Move birth to separate function, perturb an interpolation of the curr…

…ent line
OpenThermochronology · Nov 12, 2024 · d770df3 · d770df3
1 parent a897f4f
commit d770df3
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 2 deletions.
diff --git a/src/inversion.jl b/src/inversion.jl
@@ -161,10 +161,9 @@
             elseif (r < p_move+p_birth) && (npoints < maxpoints)
                 # Birth: add a new model point
                 k = npointsₚ = npoints + 1
-                agepointsₚ[k] = tnow + rand()*(tinit-tnow)
-                Tpointsₚ[k] = Tnow + rand()*(Tinit-Tnow)
                 σⱼt[k] = (tinit-Tnow)/60
                 σⱼT[k] = (Tinit-Tnow)/60
+                addpoint!(agepointsₚ, Tpointsₚ, k, σⱼt[k], σⱼT[k], boundary)
 
             elseif (r < p_move+p_birth+p_death) && (r >= p_move+p_birth) && (npoints > max(minpoints, detail.minpoints))
                 # Death: remove a model point

diff --git a/src/utilities.jl b/src/utilities.jl
@@ -190,6 +190,35 @@
         return agepointsₚ[k], Tpointsₚ[k]
     end
 
+    function addpoint!(agepointsₚ::Vector{T}, Tpointsₚ::Vector{T}, k::Int, σⱼt::T, σⱼT::T, boundary::Boundary{T}) where {T}
+        @assert eachindex(agepointsₚ) == eachindex(Tpointsₚ)
+
+        tmin, tmax = extrema(boundary.agepoints)
+        Tmin, Tmax = extrema(boundary.T₀)
+
+        # Pick an age uniformly within the boundaries
+        agepointsₚ[k] = tmin + rand()*(tmax-tmin)
+
+        # Find the closest existing points (if any)
+        ages = view(agepointsₚ, 1:k-1)
+        i₋ = findclosestbelow(agepointsₚ[k], ages)
+        i₊ = findclosestabove(agepointsₚ[k], ages)
+
+        # Pick a temperature by interpolating
+        t₋ = (firstindex(ages) <= i₋ <= lastindex(ages)) ? agepointsₚ[i₋] : tmin
+        T₋ = (firstindex(ages) <= i₋ <= lastindex(ages)) ? Tpointsₚ[i₋] : Tmin
+        t₊ = (firstindex(ages) <= i₊ <= lastindex(ages)) ? agepointsₚ[i₊] : tmax
+        T₊ = (firstindex(ages) <= i₊ <= lastindex(ages)) ? Tpointsₚ[i₊] : Tmax
+
+        f = (agepointsₚ[k] - t₋) / (t₊ - t₋)
+        Tpointsₚ[k] = f*T₊ + (1-f)*T₋
+
+        # Move the point from the interpolated value
+        movepoint!(agepointsₚ, Tpointsₚ, k, σⱼt, σⱼT, boundary)
+
+        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]