Move boundary types inside Boundary struct

OpenThermochronology · Nov 7, 2024 · 6630f94 · 6630f94
1 parent 01e5f9b
commit 6630f94
Show file tree

Hide file tree

Showing 6 changed files with 77 additions and 68 deletions.
diff --git a/examples/ZrnHeInversionVartCryst.jl b/examples/ZrnHeInversionVartCryst.jl
@@ -115,6 +115,7 @@
         ΔTinit = -50.0,             # Tinit can vary from Tinit to Tinit+ΔTinit
         Tnow = 0.0,                 # Current surface temperature (in C)
         ΔTnow = 20.0,               # Tnow may vary from Tnow to Tnow+ΔTnow
+        tnow = 0.0,                 # Today
         tinitMax = 3500.0,          # Ma -- forbid anything older than this
         minpoints = 20,             # Minimum allowed number of t-T points
         maxpoints = 50,             # Maximum allowed number of t-T points
@@ -154,10 +155,11 @@
 
     # Boundary conditions (e.g. 10C at present and 650 C at the time of zircon formation).
     boundary = Boundary(
-        agepoints = Float64[model.tnow, model.tinit],  # Ma
-        Tpoints = Float64[model.Tnow, model.Tinit],    # Degrees C
-        T₀ = Float64[model.Tnow, model.Tinit],
-        ΔT = Float64[model.ΔTnow, model.ΔTinit],
+        agepoints = [model.tnow, model.tinit],   # [Ma] Final and initial time
+        T₀ = [model.Tnow, model.Tinit],          # [C] Final and initial temperature
+        ΔT = [model.ΔTnow, model.ΔTinit],        # [C] Final and initial temperature range (positive or negative)
+        tboundary = :reflecting, # Reflecting time boundary conditions
+        Tboundary = :hard, # Hard temperature boundary conditions
     )
 
     # Default: No constraints are imposed
@@ -166,10 +168,10 @@
     # # Uncomment this section if you wish to impose an unconformity at any point in the record
     # # Uniform distributions from Age₀ to Age₀+ΔAge, T₀ to T₀+ΔT,
     # constraint = Constraint(
-    #     Age₀ = Float64[500,],         # [Ma] Minimum age
-    #     ΔAge = Float64[80,],          # [Ma] Duration of age range
-    #     T₀ = Float64[0,],             # [C] Minimum temperature
-    #     ΔT = Float64[40,],            # [C] Width of temperature range
+    #     Age₀ = [500,],         # [Ma] Age
+    #     ΔAge = [80,],          # [Ma] Age range
+    #     T₀ = [0,],             # [C] Temperature
+    #     ΔT = [40,],            # [C] Temperature range
     # )
     # name *= "_unconf"
 

diff --git a/src/inversion.jl b/src/inversion.jl
@@ -27,17 +27,13 @@
         npoints = (haskey(model, :npoints) ? model.npoints : minpoints)::Int
         totalpoints = maxpoints + boundary.npoints + constraint.npoints::Int
         simplified = (haskey(model, :simplified) ? model.simplified : false)::Bool
-        tboundary = (haskey(model, :tboundary) ? model.tboundary : :reflecting)::Symbol
-        Tboundary = (haskey(model, :Tboundary) ? model.Tboundary : :hard)::Symbol
         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
         agesteps = T.(model.agesteps)::DenseVector{T}
         tsteps = T.(model.tsteps)::DenseVector{T}
-        tinit = T(model.tinit)::T
-        tnow = zero(T) # It's always time zero today
-        Tinit = T(model.Tinit)::T
-        Tnow = T(model.Tnow)::T
+        tnow, tinit = extrema(boundary.agepoints)
+        Tnow, Tinit = extrema(boundary.T₀)
         Tr = T(haskey(model, :Tr) ? model.Tr : (Tinit+Tnow)/2)::T
         dt = T(model.dt)::T
         dr = T(model.dr)::T
@@ -169,7 +165,7 @@
             # Adjust the proposal
             if r < p_move
                 # Move one t-T point
-                movepoint!(agepointsₚ, Tpointsₚ, k, tnow+dt, tinit-dt, Tnow, Tinit, σⱼt[k], σⱼT[k], tboundary, Tboundary)
+                movepoint!(agepointsₚ, Tpointsₚ, k, σⱼt[k], σⱼT[k], boundary)
 
             elseif (r < p_move+p_birth) && (npoints < maxpoints)
                 # Birth: add a new model point
@@ -187,13 +183,10 @@
 
             elseif (r < p_move+p_birth+p_death+p_bounds)
                 # Move the temperatures of the starting and ending boundaries
-                @. boundary.Tpointsₚ = boundary.T₀ + rand()*boundary.ΔT
+                movebounds!(boundary)
+                # If there's an imposed unconformity or other t-T constraint, adjust within parameters
+                movebounds!(constraint)
 
-                # If there's an imposed unconformity, adjust within parameters
-                if constraint.npoints > 0
-                    @. constraint.agepointsₚ = constraint.Age₀ + rand()*constraint.ΔAge
-                    @. constraint.Tpointsₚ = constraint.T₀ + rand()*constraint.ΔT
-                end
             end
 
             # Recalculate interpolated proposed t-T path
@@ -296,17 +289,13 @@
         npoints = (haskey(model, :npoints) ? model.npoints : minpoints)::Int
         totalpoints = maxpoints + boundary.npoints + constraint.npoints::Int
         simplified = (haskey(model, :simplified) ? model.simplified : false)::Bool
-        tboundary = (haskey(model, :tboundary) ? model.tboundary : :reflecting)::Symbol
-        Tboundary = (haskey(model, :Tboundary) ? model.Tboundary : :hard)::Symbol
         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
         agesteps = T.(model.agesteps)::DenseVector{T}
         tsteps = T.(model.tsteps)::DenseVector{T}
-        tinit = T(model.tinit)::T
-        tnow = zero(T) # It's always time zero today
-        Tinit = T(model.Tinit)::T
-        Tnow = T(model.Tnow)::T
+        tnow, tinit = extrema(boundary.agepoints)
+        Tnow, Tinit = extrema(boundary.T₀)
         Tr = T(haskey(model, :Tr) ? model.Tr : (Tinit+Tnow)/2)::T
         dt = T(model.dt)::T
         dr = T(model.dr)::T
@@ -441,7 +430,7 @@
             # Adjust the proposal
             if r < p_move
                 # Move one t-T point
-                movepoint!(agepointsₚ, Tpointsₚ, k, tnow+dt, tinit-dt, Tnow, Tinit, σⱼt[k], σⱼT[k], tboundary, Tboundary)
+                movepoint!(agepointsₚ, Tpointsₚ, k, σⱼt[k], σⱼT[k], boundary)
 
             elseif (r < p_move+p_birth) && (npoints < maxpoints)
                 # Birth: add a new model point
@@ -459,13 +448,9 @@
 
             elseif (r < p_move+p_birth+p_death+p_bounds)
                 # Move the temperatures of the starting and ending boundaries
-                @. boundary.Tpointsₚ = boundary.T₀ + rand()*boundary.ΔT
-
-                # If there's an imposed unconformity, adjust within parameters
-                if constraint.npoints > 0
-                    @. constraint.agepointsₚ = constraint.Age₀ + rand()*constraint.ΔAge
-                    @. constraint.Tpointsₚ = constraint.T₀ + rand()*constraint.ΔT
-                end
+                movebounds!(boundary)
+                # If there's an imposed unconformity or other t-T constraint, adjust within parameters
+                movebounds!(constraint)
 
             elseif (r < p_move+p_birth+p_death+p_bounds+p_kinetics)
                 # Adjust kinetic parameters, one at a time

diff --git a/src/types.jl b/src/types.jl
@@ -50,23 +50,27 @@ end
 export RDAAM
 
 ## --- Define Boundary type to specify the working area
-struct Boundary{T<:AbstractFloat, npoints}
-    agepoints::NTuple{npoints,T} # Ma
-    Tpoints::Vector{T}   # Degrees C
-    Tpointsₚ::Vector{T}   # Degrees C
-    T₀::NTuple{npoints,T}
-    ΔT::NTuple{npoints,T}
+struct Boundary{T<:AbstractFloat}
+    agepoints::NTuple{2,T}  # Ma
+    Tpoints::Vector{T}      # Degrees C
+    Tpointsₚ::Vector{T}     # Degrees C
+    T₀::NTuple{2,T}         # Degrees C
+    ΔT::NTuple{2,T}         # Degrees C
+    tboundary::Symbol
+    Tboundary::Symbol
     npoints::Int
 end
-function Boundary(T::Type=Float64; agepoints, T₀, ΔT, Tpoints=collect(T₀))
-    @assert length(agepoints) == length(T₀) == length(ΔT) == length(Tpoints)
-    npoints = length(agepoints)
-    Boundary{T,npoints}(Tuple(T.(agepoints)), 
-        collect(T.(Tpoints)), 
-        collect(T.(Tpoints)), 
+function Boundary(T::Type=Float64; agepoints, T₀, ΔT, Tpoints=collect(T₀), tboundary=:reflecting, Tboundary=:hard)
+    @assert length(agepoints) == length(T₀) == length(ΔT) == length(Tpoints) == 2
+    Boundary{T}(
+        Tuple(T.(agepoints)), 
+        collect(T, Tpoints), 
+        collect(T, Tpoints), 
         Tuple(T.(T₀)), 
         Tuple(T.(ΔT)), 
-        npoints,
+        Symbol(tboundary),
+        Symbol(Tboundary),
+        2,
     )
 end
 export Boundary
@@ -83,9 +87,7 @@ struct Constraint{T<:AbstractFloat}
     ΔT::Vector{T}
     npoints::Int
 end
-function Constraint(T::Type=Float64; agepoints=T[], Tpoints=T[], Age₀=T[], ΔAge=T[], T₀=T[], ΔT=T[])
-    agepoints = isempty(agepoints) ? Age₀ + ΔAge/2 : agepoints
-    Tpoints = isempty(Tpoints) ? T₀ + ΔT/2 : Tpoints
+function Constraint(T::Type=Float64; Age₀=T[], ΔAge=T[], T₀=T[], ΔT=T[], agepoints=(Age₀ + ΔAge/2), Tpoints=(T₀ + ΔT/2))
     @assert length(agepoints) == length(Tpoints) == length(Age₀) == length(ΔAge) == length(T₀) == length(ΔT)
     Constraint{T}(T.(agepoints),
         T.(Tpoints),

diff --git a/src/utilities.jl b/src/utilities.jl
@@ -155,29 +155,33 @@
     end
 
     # Move a t-T point and apply boundary conditions
-    function movepoint!(agepointsₚ, Tpointsₚ, k, tmin, tmax, Tmin, Tmax, σⱼt, σⱼT, tboundary::Symbol, Tboundary::Symbol)
+    function movepoint!(agepointsₚ::Vector{T}, Tpointsₚ::Vector{T}, k::Int, σⱼt::T, σⱼT::T, boundary::Boundary{T}) where {T}
 
         # Move the age of one model point
         agepointsₚ[k] += randn() * σⱼt
 
         # Move the Temperature of one model point
         Tpointsₚ[k] += randn() * σⱼT
 
-        if tboundary === :reflecting
+        # Apply time boundary conditions
+        tmin, tmax = extrema(boundary.agepoints)
+        if boundary.tboundary === :reflecting
             agepointsₚ[k] = reflecting(agepointsₚ[k], tmin, tmax)
-        elseif tboundary === :hard
+        elseif boundary.tboundary === :hard
             agepointsₚ[k] = hard(agepointsₚ[k], tmin, tmax)
-        elseif tboundary === :periodic
+        elseif boundary.tboundary === :periodic
             agepointsₚ[k] = periodic(agepointsₚ[k], tmin, tmax)
         else
             @error "`tboundary` must be either `:reflecting`, `:hard`, or `:periodic`"
         end
 
-        if Tboundary === :reflecting
+        # Apply Temperature boundary conditions
+        Tmin, Tmax = extrema(boundary.T₀)
+        if boundary.Tboundary === :reflecting
             Tpointsₚ[k] = reflecting(Tpointsₚ[k], Tmin, Tmax)
-        elseif Tboundary === :hard
+        elseif boundary.Tboundary === :hard
             Tpointsₚ[k] = hard(Tpointsₚ[k], Tmin, Tmax)
-        elseif Tboundary === :periodic
+        elseif boundary.Tboundary === :periodic
             Tpointsₚ[k] = periodic(Tpointsₚ[k], Tmin, Tmax)
         else
             @error "`Tboundary` must be either `:reflecting`, `:hard`, or `:periodic`"
@@ -186,6 +190,20 @@
         return agepointsₚ[k], Tpointsₚ[k]
     end
 
+    function movebounds!(boundary::Boundary)
+        @inbounds for i in eachindex(boundary.Tpointsₚ)
+            boundary.Tpointsₚ[i] = boundary.T₀[i] + rand()*boundary.ΔT[i]
+        end
+        boundary
+    end
+    function movebounds!(constraint::Constraint)
+        @inbounds for i in eachindex(constraint.Tpointsₚ)
+            constraint.agepointsₚ[i] = constraint.Age₀[i] + rand()*constraint.ΔAge[i]
+            constraint.Tpointsₚ[i] = constraint.T₀[i] + rand()*constraint.ΔT[i]
+        end
+        constraint
+    end
+
     function movekinetics(zdm::ZRDAAM)
         rn = rand(1:4)
         zdmₚ = ZRDAAM(

diff --git a/test/testcomplete.jl b/test/testcomplete.jl
@@ -18,14 +18,13 @@ model = (
     Δ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
+    tnow = 0.0 ,  # Today
     tinitMax = 4000.0, # Ma -- forbid anything older than this
     minpoints = 1,  # Minimum allowed number of t-T points
     maxpoints = 40, # Maximum allowed number of t-T points
     npoints = 5, # Initial number of t-T points
     Tr = 250., # Residence temperature of initial proposal
     simplified = false, # Prefer simpler tT paths?
-    tboundary = :reflecting, # Reflecting time boundary conditions
-    Tboundary = :reflecting, # Reflecting temperature boundary conditions
     # Model uncertainty is not well known (depends on annealing parameters,
     # decay constants, diffusion parameters, etc.), but is certainly non-zero.
     # Here we add (in quadrature) a blanket model uncertainty of 25 Ma.
@@ -57,10 +56,11 @@ model = (model...,
 
 # Boundary conditions (e.g. 10C at present and 650 C at the time of zircon formation).
 boundary = Boundary(
-    agepoints = Float64[model.Tnow, model.tinit],  # Ma
-    Tpoints = Float64[model.Tnow, model.Tinit],    # Degrees C
-    T₀ = Float64[model.Tnow, model.Tinit],
-    ΔT = Float64[model.ΔTnow, model.ΔTinit],
+    agepoints = [model.tnow, model.tinit],   # [Ma] Final and initial time
+    T₀ = [model.Tnow, model.Tinit],          # [C] Final and initial temperature
+    ΔT = [model.ΔTnow, model.ΔTinit],        # [C] Final and initial temperature range (positive or negative)
+    tboundary = :reflecting, # Reflecting time boundary conditions
+    Tboundary = :reflecting, # Reflecting temperature boundary conditions
 )
 
 # Default: No unconformity is imposed

diff --git a/test/testcompletezrn.jl b/test/testcompletezrn.jl
@@ -18,6 +18,7 @@ model = (
     Δ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
+    tnow = 0.0 ,  # Today
     tinitMax = 4000.0, # Ma -- forbid anything older than this
     minpoints = 1,  # Minimum allowed number of t-T points
     maxpoints = 40, # Maximum allowed number of t-T points
@@ -56,10 +57,11 @@ model = (model...,
 
 # Boundary conditions (e.g. 10C at present and 650 C at the time of zircon formation).
 boundary = Boundary(
-    agepoints = Float64[model.Tnow, model.tinit],  # Ma
-    Tpoints = Float64[model.Tnow, model.Tinit],    # Degrees C
-    T₀ = Float64[model.Tnow, model.Tinit],
-    ΔT = Float64[model.ΔTnow, model.ΔTinit],
+    agepoints = [model.tnow, model.tinit],   # [Ma] Final and initial time
+    T₀ = [model.Tnow, model.Tinit],          # [C] Final and initial temperature
+    ΔT = [model.ΔTnow, model.ΔTinit],        # [C] Final and initial temperature range (positive or negative)
+    tboundary = :reflecting, # Reflecting time boundary conditions
+    Tboundary = :reflecting, # Reflecting temperature boundary conditions
 )
 
 # Default: No constraint is imposed