From a816d4f0f59a00deccfe69522be9e847d81b07fa Mon Sep 17 00:00:00 2001
From: "C. Brenhin Keller" <cbkeller@dartmouth.edu>
Date: Tue, 2 Jul 2024 15:00:56 -0400
Subject: [PATCH] Move utility functions out of StratMetropolis.jl

---
 src/Chron.jl           |  5 +++--
 src/Objects.jl         |  2 +-
 src/StratMetropolis.jl | 27 ---------------------------
 src/Systematic.jl      | 19 ++++++++++++++++++-
 src/Utilities.jl       | 17 ++++++++++++++---
 5 files changed, 36 insertions(+), 34 deletions(-)

diff --git a/src/Chron.jl b/src/Chron.jl
index 3b8cf7f..2c69a51 100644
--- a/src/Chron.jl
+++ b/src/Chron.jl
@@ -34,11 +34,12 @@ module Chron
     using Distributions
     using QuadGK
 
+    # Custom objects for holding Chron age data
+    include("Objects.jl")
+    # General utility functions
     include("Utilities.jl")
     # Functions for propagating systematic uncertainties
     include("Systematic.jl")
-    # Custom objects for holding Chron age data
-    include("Objects.jl")
     # Intcal2013 calibration curve for radiocarbion
     include("Intcal.jl")
     # Functions for estimating extrema of a finite-range distribution
diff --git a/src/Objects.jl b/src/Objects.jl
index f9c0efc..eee14fa 100644
--- a/src/Objects.jl
+++ b/src/Objects.jl
@@ -50,7 +50,7 @@
         Name::NTuple{N, String}
         Height::Vector{Float64}
         Height_sigma::Vector{Float64}
-        Age::Vector{Distribution{Univariate, Continuous}}
+        Age_Distribution::Vector{<:Union{<:Distribution{Univariate, Continuous}}}
         Age_Sidedness::Vector{Float64}
         Chronometer::NTuple{N, Symbol}
         Age_Unit::String
diff --git a/src/StratMetropolis.jl b/src/StratMetropolis.jl
index aef1d57..38ee59c 100644
--- a/src/StratMetropolis.jl
+++ b/src/StratMetropolis.jl
@@ -433,33 +433,6 @@
 
 ## --- # Internals of the Markov chain
 
-    # Use dispatch to let us reduce duplication
-    strat_ll(x, ages::AbstractVector{<:Radiocarbon}) = interpolate_ll(x, ages)
-    strat_ll(x, ages::AbstractVector{<:Normal}) = normpdf_ll(x, ages)
-    strat_ll(x::Real, age::Distribution) = logpdf(age, x)
-    function strat_ll(x, ages)
-        ll = zero(float(eltype(x)))
-        @inbounds for i in eachindex(x, ages)
-            ll += logpdf(ages[i], x[i])
-        end
-        return ll
-    end
-
-    adjust!(ages, chronometer, systematic::Nothing) = ages
-    function adjust!(ages, chronometer, systematic::SystematicUncertainty)
-        systUPb = randn()*systematic.UPb
-        systArAr = randn()*systematic.ArAr
-        @assert eachindex(ages)==eachindex(chronometer)
-        @inbounds for i ∈ eachindex(ages)
-            if  chronometer[i] === :UPb 
-                ages[i] += systUPb
-            elseif chronometer[i] === :ArAr
-                ages[i] += systArAr
-            end
-        end
-        return ages
-    end
-
     function stratmetropolis(Height, Height_sigma, model_heights::AbstractRange, Age_Sidedness, ages, model_ages, proposal_sigma, burnin::Integer, nsteps::Integer, sieve::Integer, Chronometer=nothing, systematic=nothing)
         resolution = step(model_heights)
         npoints = length(model_heights)
diff --git a/src/Systematic.jl b/src/Systematic.jl
index 0d8fab8..5e7017c 100644
--- a/src/Systematic.jl
+++ b/src/Systematic.jl
@@ -1,4 +1,21 @@
-## --- Dealing with systematic uncertainty
+## --- Adjust distributions for systematic uncertainty
+
+    adjust!(ages, chronometer, systematic::Nothing) = ages
+    function adjust!(ages, chronometer, systematic::SystematicUncertainty)
+        systUPb = randn()*systematic.UPb
+        systArAr = randn()*systematic.ArAr
+        @assert eachindex(ages)==eachindex(chronometer)
+        @inbounds for i ∈ eachindex(ages)
+            if  chronometer[i] === :UPb 
+                ages[i] += systUPb
+            elseif chronometer[i] === :ArAr
+                ages[i] += systArAr
+            end
+        end
+        return ages
+    end
+
+## --- Add systematic uncertainty when given an already-generated age or age-depth distribution
 
     function add_systematic_uncert_UPb(agedistmyr::Vector{<:AbstractFloat})
         λ238 = val(λ238U) # Jaffey decay constant, 1/Myr
diff --git a/src/Utilities.jl b/src/Utilities.jl
index da9de0e..ae43596 100644
--- a/src/Utilities.jl
+++ b/src/Utilities.jl
@@ -1,6 +1,3 @@
-## --- Dealing with arbitrary distributions
-
-
 ## --- "bilinear exponential" distribution type
 
     """
@@ -280,4 +277,18 @@
         return ll
     end
 
+## --- log likelihood functions allowing for arbitrary Distributions
+
+    # Use dispatch to let us reduce duplication
+    strat_ll(x, ages::AbstractVector{<:Radiocarbon}) = interpolate_ll(x, ages)
+    strat_ll(x, ages::AbstractVector{<:Normal}) = normpdf_ll(x, ages)
+    strat_ll(x::Real, age::Distribution) = logpdf(age, x)
+    function strat_ll(x, ages)
+        ll = zero(float(eltype(x)))
+        @inbounds for i in eachindex(x, ages)
+            ll += logpdf(ages[i], x[i])
+        end
+        return ll
+    end
+
 ## --- End of File