From d770df3d8cac396d7afbc95b6b91fd2158b48c14 Mon Sep 17 00:00:00 2001
From: "C. Brenhin Keller" <cbkeller@dartmouth.edu>
Date: Mon, 11 Nov 2024 23:19:25 -0500
Subject: [PATCH] Move birth to separate function, perturb an interpolation of
 the current line

---
 src/inversion.jl |  3 +--
 src/utilities.jl | 29 +++++++++++++++++++++++++++++
 2 files changed, 30 insertions(+), 2 deletions(-)

diff --git a/src/inversion.jl b/src/inversion.jl
index 15a75e2..050181a 100644
--- 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
index 8a6ca07..de8c30d 100644
--- 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]