Simplify σ, revert to 'retry' method for moving changepoints

examples/ZrnHeInversionVartCryst.jl

@@ -34,19 +34,19 @@
 
 ## --- Prepare problem
 
-    burnin = 350000
+    burnin = 350_000
     model = (
-        nsteps = 500000, # How many steps of the Markov chain should we run?
+        nsteps = 500_000, # How many steps of the Markov chain should we run?
         burnin = burnin, # How long should we wait for MC to converge (become stationary)
         dr = 1.0,    # Radius step, in microns
         dt = 10.0,   # time step size in Myr
-        dTmax = 10.0, # Maximum reheating/burial per model timestep
+        dTmax = 15.0, # Maximum reheating/burial per model timestep
         TInit = 400.0, # Initial model temperature (in C) (i.e., crystallization temperature)
         ΔTInit = -50.0, # TInit can vary from TInit to TInit+ΔTinit
         TNow = 0.0, # Current surface temperature (in C)
         ΔTNow = 10.0, # TNow may vary from TNow to TNow+ΔTNow
         tInitMax = 4000.0, # Ma -- forbid anything older than this
-        minPoints = 5,  # Minimum allowed number of t-T points
+        minPoints = 10,  # Minimum allowed number of t-T points
         maxPoints = 50, # Maximum allowed number of t-T points
         simplified = false, # Prefer simpler tT paths?
         # Diffusion parameters

src/inversion.jl

@@ -116,7 +116,7 @@
         # Simulated annealing of uncertainty
         simannealmodel = (σModel=T(model.σModel), σAnnealing=T(model.σAnnealing), λAnnealing=T(model.λAnnealing))
         σₐ = simannealsigma.(1, HeAge_sigma; simannealmodel)::DenseVector{T}
-        σₙ = simannealsigma.(nsteps, HeAge_sigma; simannealmodel)::DenseVector{T}
+        σ = sqrt.(T(model.σModel)^2 .+ HeAge_sigma.^2)
 
         # Log-likelihood for initial proposal
         ll = normpdf_ll(HeAge, σₐ, calcHeAges)
@@ -195,6 +195,10 @@
 
                 # Retry unless we have satisfied the maximum reheating rate
                 (maxdiff(TStepsₚ) < dTmax) && break
+
+                # Copy last accepted solution to re-modify if we don't break
+                copyto!(agePointsₚ, agePoints)
+                copyto!(TPointsₚ, TPoints)
                 end
             elseif (r < move+birth) && (nPointsₚ < maxPoints)
                 # Birth: add a new model point
@@ -246,6 +250,11 @@
 
                 # Retry unless we have satisfied the maximum reheating rate
                 (maxdiff(TStepsₚ) < dTmax) && break
+
+                # Copy last accepted solution to re-modify if we don't break
+                copyto!(unconf.agePointsₚ, unconf.agePoints)
+                copyto!(unconf.TPointsₚ, unconf.TPoints)
+                copyto!(boundary.TPointsₚ, boundary.TPoints)
                 end
             end
 
@@ -283,13 +292,13 @@
                 copyto!(boundary.TPoints, boundary.TPointsₚ)
                 copyto!(calcHeAges, calcHeAgesₚ)
 
-                # This is saved for ouput, but not critical to the function of the MCMC loop
+                # Not critical to the function of the MCMC loop, but critical for recording stationary distribution!
                 copyto!(TSteps, TStepsₚ)
                 acceptancedist[n] = true
             end
 
             # Record results for analysis and troubleshooting
-            lldist[n] = normpdf_ll(HeAge, σₙ, calcHeAges) # Recalculated to constant baseline
+            lldist[n] = normpdf_ll(HeAge, σ, calcHeAges) # Recalculated to constant baseline
             ndist[n] = nPoints # distribution of # of points
             HeAgedist[:,n] .= calcHeAges # distribution of He ages