Add option

OpenThermochronology · Dec 13, 2024 · 05657d4 · 05657d4
1 parent 6468bd6
commit 05657d4
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 6 deletions.
diff --git a/src/chronometers.jl b/src/chronometers.jl
@@ -735,4 +735,18 @@ function chronometers(T::Type{<:AbstractFloat}, data, model)
 
     isempty(result) && @error "No chronometers found"
     return unionize(result)
+end
+
+## --- model uncertainties of chronometers
+
+function movesigma!(σcalc::AbstractVector{T}, chrons::AbstractVector{<:Chronometer}) where {T<:AbstractFloat}
+    for C in (ZirconHe, ApatiteHe, ApatiteFT)
+        r = abs(randn(T))
+        for i in eachindex(σcalc, chrons)
+            if chrons[i] isa C
+                σcalc[i] *= r
+            end
+        end
+    end
+    return σcalc
 end
diff --git a/src/inversion.jl b/src/inversion.jl
@@ -65,6 +65,7 @@
         npoints = (haskey(model, :npoints) ? model.npoints : minpoints)::Int
         totalpoints = maxpoints + boundary.npoints + constraint.npoints::Int
         simplified = (haskey(model, :simplified) ? model.simplified : false)::Bool
+        dynamicsigma = (haskey(model, :dynamicsigma) ? model.dynamicsigma : true)::Bool
         dynamicjumping = (haskey(model, :dynamicjumping) ? model.dynamicjumping : false)::Bool
         dTmax = T(haskey(model, :dTmax) ? model.dTmax : 10)::T
         dTmax_sigma = T(haskey(model, :dTmax_sigma) ? model.dTmax_sigma : 5)::T
@@ -76,9 +77,9 @@
         Tnow, Tinit = extrema(boundary.T₀)
         Tr = T(haskey(model, :Tr) ? model.Tr : (Tinit+Tnow)/2)::T
         dt = T(model.dt)::T
-        σmodel = T(model.σmodel)::T
-        σannealing = T(model.σannealing)::T
-        λannealing = T(model.λannealing)::T
+        σmodel = T(haskey(model, :σmodel) ? model.σmodel : 0)::T
+        σannealing = T(haskey(model, :σannealing) ? model.σannealing : 125)::T
+        λannealing = T(haskey(model, :λannealing) ? model.λannealing : 2/burnin)::T
 
         # Arrays to hold all t and T points (up to npoints=maxpoints)
         agepoints = zeros(T, maxpoints) 
@@ -155,6 +156,7 @@
             copyto!(constraint.agepointsₚ, constraint.agepoints)
             copyto!(constraint.Tpointsₚ, constraint.Tpoints)
             copyto!(boundary.Tpointsₚ, boundary.Tpoints)
+            copyto!(σcalcₚ, σcalc)
             copyto!(σⱼtₚ, σⱼt)
             copyto!(σⱼTₚ, σⱼT)
 
@@ -186,6 +188,9 @@
                 # If there's an imposed unconformity or other t-T constraint, adjust within bounds
                 movebounds!(constraint, boundary)
 
+                # Move uncertainties
+                dynamicsigma && movesigma!(σcalcₚ, data)
+
             end
 
             # Recalculate interpolated proposed t-T path
@@ -269,6 +274,7 @@
             copyto!(constraint.agepointsₚ, constraint.agepoints)
             copyto!(constraint.Tpointsₚ, constraint.Tpoints)
             copyto!(boundary.Tpointsₚ, boundary.Tpoints)
+            copyto!(σcalcₚ, σcalc)
             copyto!(σⱼtₚ, σⱼt)
             copyto!(σⱼTₚ, σⱼT)
 
@@ -300,6 +306,9 @@
                 # If there's an imposed unconformity or other t-T constraint, adjust within bounds
                 movebounds!(constraint, boundary)
 
+                # Move uncertainties
+                dynamicsigma && movesigma!(σcalcₚ, data)
+
             end
 
             # Recalculate interpolated proposed t-T path
@@ -411,6 +420,7 @@
         npoints = (haskey(model, :npoints) ? model.npoints : minpoints)::Int
         totalpoints = maxpoints + boundary.npoints + constraint.npoints::Int
         simplified = (haskey(model, :simplified) ? model.simplified : false)::Bool
+        dynamicsigma = (haskey(model, :dynamicsigma) ? model.dynamicsigma : true)::Bool
         dynamicjumping = (haskey(model, :dynamicjumping) ? model.dynamicjumping : false)::Bool
         dTmax = T(haskey(model, :dTmax) ? model.dTmax : 10)::T
         dTmax_sigma = T(haskey(model, :dTmax_sigma) ? model.dTmax_sigma : 5)::T
@@ -422,9 +432,9 @@
         Tnow, Tinit = extrema(boundary.T₀)
         Tr = T(haskey(model, :Tr) ? model.Tr : (Tinit+Tnow)/2)::T
         dt = T(model.dt)::T
-        σmodel = T(model.σmodel)::T
-        σannealing = T(model.σannealing)::T
-        λannealing = T(model.λannealing)::T
+        σmodel = T(haskey(model, :σmodel) ? model.σmodel : 0)::T
+        σannealing = T(haskey(model, :σannealing) ? model.σannealing : 125)::T
+        λannealing = T(haskey(model, :λannealing) ? model.λannealing : 2/burnin)::T
 
         # Arrays to hold all t and T points (up to npoints=maxpoints)
         agepoints = zeros(T, maxpoints) 
@@ -504,6 +514,7 @@
             copyto!(constraint.agepointsₚ, constraint.agepoints)
             copyto!(constraint.Tpointsₚ, constraint.Tpoints)
             copyto!(boundary.Tpointsₚ, boundary.Tpoints)
+            copyto!(σcalcₚ, σcalc)
             copyto!(σⱼtₚ, σⱼt)
             copyto!(σⱼTₚ, σⱼT)
 
@@ -535,6 +546,9 @@
                 # If there's an imposed unconformity or other t-T constraint, adjust within bounds
                 movebounds!(constraint, boundary)
 
+                # Move uncertainties
+                dynamicsigma && movesigma!(σcalcₚ, data)
+
             elseif (r < p_move+p_birth+p_death+p_bounds+p_kinetics)
                 # Adjust kinetic parameters, one at a time
                 haszhe && (zdmₚ = movekinetics(zdm))
@@ -630,6 +644,7 @@
             copyto!(constraint.agepointsₚ, constraint.agepoints)
             copyto!(constraint.Tpointsₚ, constraint.Tpoints)
             copyto!(boundary.Tpointsₚ, boundary.Tpoints)
+            copyto!(σcalcₚ, σcalc)
             copyto!(σⱼtₚ, σⱼt)
             copyto!(σⱼTₚ, σⱼT)
 
@@ -661,6 +676,9 @@
                 # If there's an imposed unconformity or other t-T constraint, adjust within bounds
                 movebounds!(constraint, boundary)
 
+                # Move uncertainties
+                dynamicsigma && movesigma!(σcalcₚ, data)
+
             elseif (r < p_move+p_birth+p_death+p_bounds+p_kinetics)
                 # Adjust kinetic parameters, one at a time
                 haszhe && (zdmₚ = movekinetics(zdm))