Move utility functions out of StratMetropolis.jl

brenhinkeller · Jul 2, 2024 · a816d4f · a816d4f
1 parent ec24946
commit a816d4f
Show file tree

Hide file tree

Showing 5 changed files with 36 additions and 34 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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