Add and test StratMetropolis functions that take GeneralAgeData inputs

Author: brenhinkeller

Project.toml

@@ -1,7 +1,7 @@
 name = "Chron"
 uuid = "68885b1f-77b5-52a7-b2e7-6a8014c56b98"
 authors = ["C. Brenhin Keller <cbkeller@dartmouth.edu>"]
-version = "0.6.0"
+version = "0.6.1"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/StratMetropolis.jl

@@ -4,6 +4,7 @@
     """
     ```julia
     StratMetropolis(smpl::ChronAgeData, [hiatus::HiatusData,] config::StratAgeModelConfiguration)
+    StratMetropolis(smpl::GeneralAgeData, [hiatus::HiatusData,] config::StratAgeModelConfiguration)
     ```
     Runs the main Chron.jl age-depth model routine for a stratigraphic set of
     samples defined by sample heights and simple Gaussian age constraints in the
@@ -39,7 +40,7 @@
         Height = copy(smpl.Height)::Vector{Float64}
         Height_sigma = smpl.Height_sigma::Vector{Float64} .+ 1E-9 # Avoid divide-by-zero issues
         Age_Sidedness = copy(smpl.Age_Sidedness)::Vector{Float64} # Bottom is a maximum age and top is a minimum age
-        Chronometer = smpl.Chronometer
+        chronometer = smpl.Chronometer
         (bottom, top) = extrema(Height)
         model_heights = bottom:resolution:top
 
@@ -60,7 +61,7 @@
             Height_sigma = [0; Height_sigma; 0] .+ 1E-9 # Avoid divide-by-zero issues
             Age_Sidedness = [-1.0; Age_Sidedness; 1.0;] # Bottom is a maximum age and top is a minimum age
             model_heights = (bottom-offset):resolution:(top+offset)
-            Chronometer = (:None, Chronometer..., :None)
+            chronometer = (:None, chronometer..., :None)
         end
         active_height_t = bottom .<= model_heights .<= top
 
@@ -77,7 +78,7 @@
 
         # Run the Markov chain
         ages = Normal.(Age, Age_sigma)
-        agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, Chronometer, systematic)
+        agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, chronometer, systematic)
 
         # Crop the result
         agedist = agedist[active_height_t,:]
@@ -86,7 +87,7 @@
 
         return mdl, agedist, lldist
     end
-    function StratMetropolis(smpl::ChronAgeData, hiatus::HiatusData, config::StratAgeModelConfiguration)
+    function StratMetropolis(smpl::ChronAgeData, hiatus::HiatusData, config::StratAgeModelConfiguration, systematic=nothing)
         # Run stratigraphic MCMC model, with hiata
         @info "Generating stratigraphic age-depth model..."
 
@@ -103,6 +104,7 @@
         Height = copy(smpl.Height)::Vector{Float64}
         Height_sigma = smpl.Height_sigma::Vector{Float64} .+ 1E-9 # Avoid divide-by-zero issues
         Age_Sidedness = copy(smpl.Age_Sidedness)::Vector{Float64} # Bottom is a maximum age and top is a minimum age
+        chronometer = smpl.Chronometer
         (bottom, top) = extrema(Height)
         model_heights = bottom:resolution:top
 
@@ -123,6 +125,7 @@
             Height_sigma = [0; Height_sigma; 0] .+ 1E-9 # Avoid divide-by-zero issues
             Age_Sidedness = [-1.0; Age_Sidedness; 1.0;] # Bottom is a maximum age and top is a minimum age
             model_heights = (bottom-offset):resolution:(top+offset)
+            chronometer = (:None, chronometer..., :None)
         end
         active_height_t = bottom .<= model_heights .<= top
         npoints = length(model_heights)
@@ -140,7 +143,133 @@
 
         # Run the Markov chain
         ages = Normal.(Age, Age_sigma)
-        agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve)
+        agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, chronometer, systematic)
+
+        # Crop the result
+        agedist = agedist[active_height_t,:]
+        model_heights = model_heights[active_height_t]
+        mdl = StratAgeModel(model_heights, agedist)
+
+        return mdl, agedist, hiatusdist, lldist
+    end
+    function StratMetropolis(smpl::GeneralAgeData, config::StratAgeModelConfiguration, systematic=nothing)
+        # Run stratigraphic MCMC model
+        @info "Generating stratigraphic age-depth model..."
+
+        # Model configuration -- read from struct
+        resolution = config.resolution
+        burnin = config.burnin
+        nsteps = config.nsteps
+        sieve = config.sieve
+        bounding = config.bounding
+
+        # Stratigraphic age constraints
+        ages = unionize(smpl.Age_Distribution)::Vector{<:Union{<:Distribution{Univariate, Continuous}}}
+        Height = copy(smpl.Height)::Vector{Float64}
+        Height_sigma = smpl.Height_sigma::Vector{Float64} .+ 1E-9 # Avoid divide-by-zero issues
+        Age_Sidedness = copy(smpl.Age_Sidedness)::Vector{Float64} # Bottom is a maximum age and top is a minimum age
+        chronometer = smpl.Chronometer
+        (bottom, top) = extrema(Height)
+        model_heights = bottom:resolution:top
+
+        aveuncert = nanmean(std.(ages))
+        absdiff = diff(sort!(mean.(ages[Age_Sidedness.==0])))
+        maxdiff = isempty(absdiff) ? 0.0 : nanmaximum(absdiff)
+        proposal_sigma = sqrt(aveuncert^2 + (maxdiff/10)^2)
+
+        if bounding>0
+            # If bounding is requested, add extrapolated top and bottom bounds to avoid
+            # issues with the stratigraphic markov chain wandering off to +/- infinity
+            (youngest, oldest) = extrema(mean.(ages))
+            dt_dH = (oldest-youngest)/(top-bottom)
+            offset = round((top-bottom)*bounding/resolution)*resolution
+            ages = unionize([Normal(oldest+offset*dt_dH, aveuncert/10); 
+                             ages; 
+                             Normal(youngest-offset*dt_dH, aveuncert/10)])
+            Height = [bottom-offset; Height; top+offset]
+            Height_sigma = [0; Height_sigma; 0] .+ 1E-9 # Avoid divide-by-zero issues
+            Age_Sidedness = [-1.0; Age_Sidedness; 1.0;] # Bottom is a maximum age and top is a minimum age
+            model_heights = (bottom-offset):resolution:(top+offset)
+            chronometer = (:None, chronometer..., :None)
+        end
+        active_height_t = bottom .<= model_heights .<= top
+
+        # Start with a linear fit as an initial proposal
+        (a,b) = hcat(fill!(similar(Height), 1), Height) \ mean.(ages)
+        model_ages = a .+ b .* collect(model_heights)
+
+        # Select sidedness method
+        sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
+            FastSidedness(Age_Sidedness)
+        else
+            CDFSidedness(Age_Sidedness)
+        end
+
+        # Run the Markov chain
+        agedist, lldist = stratmetropolis(Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, chronometer, systematic)
+
+        # Crop the result
+        agedist = agedist[active_height_t,:]
+        model_heights = model_heights[active_height_t]
+        mdl = StratAgeModel(model_heights, agedist)
+
+        return mdl, agedist, lldist
+    end
+    function StratMetropolis(smpl::GeneralAgeData, hiatus::HiatusData, config::StratAgeModelConfiguration, systematic=nothing)
+        # Run stratigraphic MCMC model
+        @info "Generating stratigraphic age-depth model..."
+
+        # Model configuration -- read from struct
+        resolution = config.resolution
+        burnin = config.burnin
+        nsteps = config.nsteps
+        sieve = config.sieve
+        bounding = config.bounding
+
+        # Stratigraphic age constraints
+        ages = unionize(smpl.Age_Distribution)::Vector{<:Union{<:Distribution{Univariate, Continuous}}}
+        Height = copy(smpl.Height)::Vector{Float64}
+        Height_sigma = smpl.Height_sigma::Vector{Float64} .+ 1E-9 # Avoid divide-by-zero issues
+        Age_Sidedness = copy(smpl.Age_Sidedness)::Vector{Float64} # Bottom is a maximum age and top is a minimum age
+        chronometer = smpl.Chronometer
+        (bottom, top) = extrema(Height)
+        model_heights = bottom:resolution:top
+
+        aveuncert = nanmean(std.(ages))
+        absdiff = diff(sort!(mean.(ages[Age_Sidedness.==0])))
+        maxdiff = isempty(absdiff) ? 0.0 : nanmaximum(absdiff)
+        proposal_sigma = sqrt(aveuncert^2 + (maxdiff/10)^2)
+
+        if bounding>0
+            # If bounding is requested, add extrapolated top and bottom bounds to avoid
+            # issues with the stratigraphic markov chain wandering off to +/- infinity
+            (youngest, oldest) = extrema(mean.(ages))
+            dt_dH = (oldest-youngest)/(top-bottom)
+            offset = round((top-bottom)*bounding/resolution)*resolution
+            ages = unionize([Normal(oldest+offset*dt_dH, aveuncert/10); 
+                             ages; 
+                             Normal(youngest-offset*dt_dH, aveuncert/10)])
+            Height = [bottom-offset; Height; top+offset]
+            Height_sigma = [0; Height_sigma; 0] .+ 1E-9 # Avoid divide-by-zero issues
+            Age_Sidedness = [-1.0; Age_Sidedness; 1.0;] # Bottom is a maximum age and top is a minimum age
+            model_heights = (bottom-offset):resolution:(top+offset)
+            chronometer = (:None, chronometer..., :None)
+        end
+        active_height_t = bottom .<= model_heights .<= top
+
+        # Start with a linear fit as an initial proposal
+        (a,b) = hcat(fill!(similar(Height), 1), Height) \ mean.(ages)
+        model_ages = a .+ b .* collect(model_heights)
+
+        # Select sidedness method
+        sidedness = if smpl.Sidedness_Method === :fast || all(iszero, smpl.Age_Sidedness)
+            FastSidedness(Age_Sidedness)
+        else
+            CDFSidedness(Age_Sidedness)
+        end
+
+        # Run the Markov chain
+        agedist, lldist, hiatusdist = stratmetropolis(hiatus, Height, Height_sigma, model_heights, sidedness, ages, model_ages, proposal_sigma, burnin, nsteps, sieve, chronometer, systematic)
 
         # Crop the result
         agedist = agedist[active_height_t,:]
@@ -401,9 +530,6 @@
 
         return mdl, agedist, lldist
     end
-
-## --- Stratigraphic MCMC model with hiatus, for radiocarbon ages # # # # # #
-
     function StratMetropolis14C(smpl::ChronAgeData, hiatus::HiatusData, config::StratAgeModelConfiguration)
         # Run stratigraphic MCMC model, with hiata
         @info "Generating stratigraphic age-depth model..."

test/runtests.jl

@@ -4,6 +4,9 @@ using Test, Statistics, Distributions
 const make_plots = get(ENV, "MAKE_PLOTS", false) == "true"
 @testset "Utilities" begin include("testUtilities.jl") end
 @testset "Eruption / deposition age distributions" begin include("testDist.jl") end
-@testset "Strat only" begin include("testStratOnly.jl") end
+@testset "Strat only" begin 
+    include("testStratOnly.jl") 
+    include("testStratOnlyGeneral.jl") 
+end
 @testset "Radiocarbon" begin include("testRadiocarbon.jl") end
 @testset "Coupled model" begin include("testCoupled.jl") end

test/testCoupled.jl

@@ -1,11 +1,11 @@
 ## --- Standard coupled eruption-deposition modelling
 
-nSamples = 5 # The number of samples you have data for
-smpl = ChronAgeData(nSamples)
+nsamples = 5 # The number of samples you have data for
+smpl = ChronAgeData(nsamples)
 smpl.Name      =  ("KJ08-157", "KJ04-75", "KJ09-66",  "KJ04-72", "KJ04-70",)
 smpl.Height   .=  [     -52.0,      44.0,      54.0,      82.0,      93.0,]
 smpl.Height_sigma .= [    3.0,       1.0,       3.0,       3.0,       3.0,]
-smpl.Age_Sidedness .= zeros(nSamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
+smpl.Age_Sidedness .= zeros(nsamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
 smpl.Path = abspath("../examples/DenverUPbExampleData/") # Where are the data files?
 smpl.inputSigmaLevel = 2 # i.e., are the data files 1-sigma or 2-sigma. Integer.
 smpl.Age_Unit = "Ma" # Unit of measurement for ages and errors in the data files
@@ -124,11 +124,11 @@ println("StratMetropolisDist with fitted Gaussians:")
 
 ## -- Test coupled with Isoplot.jl Pb-loss-aware eruption ages
 
-nSamples = 3 # The number of samples you have data for
-smpl = ChronAgeData(nSamples)
+nsamples = 3 # The number of samples you have data for
+smpl = ChronAgeData(nsamples)
 smpl.Name      =  ("KR18-04", "KR18-01", "KR18-05")
 smpl.Height   .=  [      0.0,     100.0,     200.0] # Arbitrary heights
-smpl.Age_Sidedness .= zeros(nSamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
+smpl.Age_Sidedness .= zeros(nsamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
 smpl.Path = abspath("../examples/ConcordiaExampleData/") # Where are the data files?
 smpl.inputSigmaLevel = 1 # i.e., are the data files 1-sigma or 2-sigma. Integer.
 smpl.Age_Unit = "Ma" # Unit of measurement for ages and errors in the data files

test/testRadiocarbon.jl

@@ -7,15 +7,15 @@
 ## --- Test age-depth model
 
 # Input the number of samples we wish to model (must match below)
-nSamples = 4
-# Make an instance of a ChronSection object for nSamples
-smpl = ChronAgeData(nSamples)
+nsamples = 4
+# Make an instance of a ChronSection object for nsamples
+smpl = ChronAgeData(nsamples)
 smpl.Name           = ("Sample 1", "Sample 2", "Sample 3", "Sample 4") # Et cetera
 smpl.Age_14C       .= [ 6991,  7088,  7230,  7540,] # Measured ages
 smpl.Age_14C_sigma .= [   30,    70,    50,    50,] # Measured 1-σ uncertainties
 smpl.Height        .= [ -355,  -380,-397.0,-411.5,] # Depths below surface should be negative
-smpl.Height_sigma  .= fill(0.01, nSamples) # Usually assume little or no sample height uncertainty
-smpl.Age_Sidedness .= zeros(nSamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
+smpl.Height_sigma  .= fill(0.01, nsamples) # Usually assume little or no sample height uncertainty
+smpl.Age_Sidedness .= zeros(nsamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
 smpl.Age_Unit = "Years BP" # Unit of measurement for ages
 smpl.Height_Unit = "m" # Unit of measurement for Height and Height_sigma
 
@@ -24,8 +24,8 @@ smpl.Height_Unit = "m" # Unit of measurement for Height and Height_sigma
 calibration = intcal13
 
 # Calculate calendar age PDFs for each sample
-smpl.Params = fill(NaN, length(calibration.Age_Calendar), nSamples)
-for i = 1:nSamples
+smpl.Params = fill(NaN, length(calibration.Age_Calendar), nsamples)
+for i = 1:nsamples
     # The likelihood that a measured 14C age could result from a sample of
     # a given calendar age is proportional to the intergral of the product
     # of the two respective distributions
@@ -98,8 +98,8 @@ hiatus.Duration_sigma = [  30.5,   20.0 ]
 calibration = intcal20
 
 # Calculate calendar age PDFs for each sample
-smpl.Params = fill(NaN, length(calibration.Age_Calendar), nSamples)
-for i = 1:nSamples
+smpl.Params = fill(NaN, length(calibration.Age_Calendar), nsamples)
+for i = 1:nsamples
     # The likelihood that a measured 14C age could result from a sample of
     # a given calendar age is proportional to the intergral of the product
     # of the two respective distributions

test/testStratOnly.jl

@@ -1,13 +1,13 @@
-# Make an instance of a ChronAgeData object for nSamples
-nSamples = 6
-smpl = ChronAgeData(nSamples)
+# Make an instance of a ChronAgeData object for nsamples
+nsamples = 6
+smpl = ChronAgeData(nsamples)
 @test smpl isa ChronAgeData
 smpl.Name          = ("minimum age", "Sample 1", "Sample 2", "Sample 3", "Sample 4", "maximum age") # Et cetera
 smpl.Age          .= [ 690.0,  699.1,  708.8,  723.0,  754.0,  812.0] # Measured ages
 smpl.Age_sigma    .= [   7.0,    3.0,    7.0,    5.0,    5.0,    6.0] # Measured 1-σ uncertainties
 smpl.Height       .= [-350.0, -355.0, -380.0, -397.0, -411.5, -420.0] # Depths below surface should be negative
-smpl.Height_sigma .= fill(0.01, nSamples) # Usually assume little or no sample height uncertainty
-smpl.Age_Sidedness .= zeros(nSamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
+smpl.Height_sigma .= fill(0.01, nsamples) # Usually assume little or no sample height uncertainty
+smpl.Age_Sidedness .= zeros(nsamples) # Sidedness (zeros by default: geochron constraints are two-sided). Use -1 for a maximum age and +1 for a minimum age, 0 for two-sided
 smpl.Age_Sidedness[1] = 1. # Minimum age
 smpl.Age_Sidedness[end] = -1. # Maximum age
 smpl.Age_Unit = "Years BP" # Unit of measurement for ages
@@ -60,10 +60,3 @@ hiatus.Duration_sigma = [   3.1,    2.0 ]
 @test size(hiatusdist) == (nHiatuses, config.nsteps)
 @test mean(hiatusdist, dims=2) ≈ [10.580012942504894; 18.96167245288326;;] atol=2
 @test -Inf < mean(lldist) < 0
-
-## --- Strat only, general Distributions-based case
-
-# Make an instance of a ChronAgeData object for nSamples
-nSamples = 4
-smpl = GeneralAgeData(nSamples)
-@test smpl isa GeneralAgeData