diff --git a/Project.toml b/Project.toml index 8f77394..a72495c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Chron" uuid = "68885b1f-77b5-52a7-b2e7-6a8014c56b98" authors = ["C. Brenhin Keller "] -version = "0.6.1" +version = "0.6.2" [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" diff --git a/examples/Chron1.0StratOnlyGeneral.jl b/examples/Chron1.0StratOnlyGeneral.jl index 085aa11..2c6e2ba 100644 --- a/examples/Chron1.0StratOnlyGeneral.jl +++ b/examples/Chron1.0StratOnlyGeneral.jl @@ -18,10 +18,10 @@ nsamples = 5 # Make an instance of a GeneralAgeData object for n samples smpl = GeneralAgeData(nsamples) - smpl.Name = ( "Sample 1", "Sample 2", "Sample 3", "Sample 4", "Sample 5",) # Et cetera - smpl.Age_Distribution = [Normal(39.5,0.1), Uniform(37, 38), Normal(36.3, 0.1), Uniform(33.5,34), Normal(32.1, 0.1),] # Measured ages - smpl.Height = [ 100, 200, 300, 400, 500,] # Depths below surface should be negative - 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.Name = ( "Sample 1", "Sample 2", "Sample 3", "Sample 4", "Sample 5",) # Et cetera + smpl.Age = [Normal(39.5,0.1), Uniform(37, 38), Normal(36.3, 0.1), Uniform(33.5,34), Normal(32.1, 0.1),] # Measured ages + smpl.Height = [ 100, 200, 300, 400, 500,] # Depths below surface should be negative + 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 = "Ma" # Unit of measurement for ages smpl.Height_Unit = "m" # Unit of measurement for Height and Height_sigma @@ -60,13 +60,13 @@ plot!(hdl, [mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(round(Int,minimum(mdl.Height)),0.5,:blue), label="model") plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="") # Center line t = smpl.Age_Sidedness .== 0 # Two-sided constraints (plot in black) - any(t) && plot!(hdl, mean.(smpl.Age_Distribution[t]), smpl.Height[t], xerror=2*std.(smpl.Age_Distribution[t]),label="data",seriestype=:scatter,color=:black) + any(t) && plot!(hdl, mean.(smpl.Age[t]), smpl.Height[t], xerror=2*std.(smpl.Age[t]),label="data",seriestype=:scatter,color=:black) t = smpl.Age_Sidedness .== 1 # Minimum ages (plot in cyan) - any(t) && plot!(hdl, mean.(smpl.Age_Distribution[t]), smpl.Height[t], xerror=(2*std.(smpl.Age_Distribution[t]),zeros(count(t))),label="",seriestype=:scatter,color=:cyan,msc=:cyan) - any(t) && zip(mean.(smpl.Age_Distribution[t]), mean.(smpl.Age_Distribution[t]).+nanmean(std.(smpl.Age_Distribution[t]))*4, smpl.Height[t]) .|> x-> plot!([x[1],x[2]],[x[3],x[3]], arrow=true, label="", color=:cyan) + any(t) && plot!(hdl, mean.(smpl.Age[t]), smpl.Height[t], xerror=(2*std.(smpl.Age[t]),zeros(count(t))),label="",seriestype=:scatter,color=:cyan,msc=:cyan) + any(t) && zip(mean.(smpl.Age[t]), mean.(smpl.Age[t]).+nanmean(std.(smpl.Age[t]))*4, smpl.Height[t]) .|> x-> plot!([x[1],x[2]],[x[3],x[3]], arrow=true, label="", color=:cyan) t = smpl.Age_Sidedness .== -1 # Maximum ages (plot in orange) - any(t) && plot!(hdl, mean.(smpl.Age_Distribution[t]), smpl.Height[t], xerror=(zeros(count(t)),2*std.(smpl.Age_Distribution[t])),label="",seriestype=:scatter,color=:orange,msc=:orange) - any(t) && zip(mean.(smpl.Age_Distribution[t]), mean.(smpl.Age_Distribution[t]).-nanmean(std.(smpl.Age_Distribution[t]))*4, smpl.Height[t]) .|> x-> plot!([x[1],x[2]],[x[3],x[3]], arrow=true, label="", color=:orange) + any(t) && plot!(hdl, mean.(smpl.Age[t]), smpl.Height[t], xerror=(zeros(count(t)),2*std.(smpl.Age[t])),label="",seriestype=:scatter,color=:orange,msc=:orange) + any(t) && zip(mean.(smpl.Age[t]), mean.(smpl.Age[t]).-nanmean(std.(smpl.Age[t]))*4, smpl.Height[t]) .|> x-> plot!([x[1],x[2]],[x[3],x[3]], arrow=true, label="", color=:orange) savefig(hdl,"AgeDepthModel.pdf") display(hdl) @@ -147,7 +147,7 @@ # Plot results (mean and 95% confidence interval for both model and data) hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label="model") plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label="", fg_color_legend=:white) - plot!(hdl, mean.(smpl.Age_Distribution), smpl.Height, xerror=std.(smpl.Age_Distribution)*2,label="data",seriestype=:scatter,color=:black) + plot!(hdl, mean.(smpl.Age), smpl.Height, xerror=std.(smpl.Age)*2,label="data",seriestype=:scatter,color=:black) plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Height ($(smpl.Height_Unit))") savefig(hdl, "Interpolated age distribution.pdf") display(hdl) diff --git a/src/Objects.jl b/src/Objects.jl index e7b5228..bf65600 100644 --- a/src/Objects.jl +++ b/src/Objects.jl @@ -24,22 +24,22 @@ Sidedness_Method::Symbol end - function ChronAgeData(nSamples::Integer) - ChronAgeData{nSamples}( - ntuple(i->"Sample name", nSamples), - collect(1.0:nSamples), # Sample Height - zeros(nSamples), # Sample Height_sigma - fill(NaN,nSamples), # Sample ages - fill(NaN,nSamples), # Sample age uncertainty - fill(NaN,nSamples), # Sample age 2.5% CI - fill(NaN,nSamples), # Sample age 97.5% CI - fill(NaN,nSamples), # Sample 14C ages - fill(NaN,nSamples), # Sample 14C uncertainties - 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 - zeros(nSamples), # DistType (for Distribution-fitting only: 0=distribution to be fit, 1=single Gaussian) - Vector{Vector{Float64}}(undef,nSamples), # Stationary distribution of eruption age - ntuple(i->:Chronometer, nSamples), # Age Types (e.g., :UPb or :ArAr) - fill(NaN,5,nSamples), # Sample age distribution parameters + function ChronAgeData(nsamples::Integer) + ChronAgeData{nsamples}( + ntuple(i->"Sample name", nsamples), + collect(1.0:nsamples), # Sample Height + zeros(nsamples), # Sample Height_sigma + fill(NaN,nsamples), # Sample ages + fill(NaN,nsamples), # Sample age uncertainty + fill(NaN,nsamples), # Sample age 2.5% CI + fill(NaN,nsamples), # Sample age 97.5% CI + fill(NaN,nsamples), # Sample 14C ages + fill(NaN,nsamples), # Sample 14C uncertainties + 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 + zeros(nsamples), # DistType (for Distribution-fitting only: 0=distribution to be fit, 1=single Gaussian) + Vector{Vector{Float64}}(undef,nsamples), # Stationary distribution of eruption age + ntuple(i->:Chronometer, nsamples), # Age Types (e.g., :UPb or :ArAr) + fill(NaN,5,nsamples), # Sample age distribution parameters "./", # Relative path where we can find .csv data files 2, # i.e., are the data files 1-sigma or 2-sigma "Ma", @@ -52,7 +52,7 @@ Name::NTuple{N, String} Height::Vector{Float64} Height_sigma::Vector{Float64} - Age_Distribution::Vector{<:Union{<:Distribution{Univariate, Continuous}}} + Age::Vector{<:Union{<:Distribution{Univariate, Continuous}}} Age_Sidedness::Vector{Float64} Chronometer::NTuple{N, Symbol} Age_Unit::String @@ -60,14 +60,14 @@ Sidedness_Method::Symbol end - function GeneralAgeData(nSamples::Integer) - GeneralAgeData{nSamples}( - ntuple(i->"Sample name", nSamples), - collect(1.0:nSamples), # Sample Height - zeros(nSamples), # Sample Height_sigma - fill(Uniform(0,4567),nSamples), # Sample ages - 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 - ntuple(i->:Chronometer, nSamples), # Age Types (e.g., :UPb or :ArAr) + function GeneralAgeData(nsamples::Integer) + GeneralAgeData{nsamples}( + ntuple(i->"Sample name", nsamples), + collect(1.0:nsamples), # Sample Height + zeros(nsamples), # Sample Height_sigma + fill(Uniform(0,4567), nsamples), # Sample ages + 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 + ntuple(i->:Chronometer, nsamples), # Age Types (e.g., :UPb or :ArAr) "Ma", "m", :cdf, @@ -85,12 +85,12 @@ Height_Unit::String end - function HiatusData(nHiatuses::Integer) + function HiatusData(nhiatuses::Integer) HiatusData( - fill(NaN,nHiatuses), # Height - fill(NaN,nHiatuses), # Height_sigma - fill(NaN,nHiatuses), # Duration - fill(NaN,nHiatuses), # Duration_sigma + fill(NaN, nhiatuses), # Height + fill(NaN, nhiatuses), # Height_sigma + fill(NaN, nhiatuses), # Duration + fill(NaN, nhiatuses), # Duration_sigma "Ma", "m", ) diff --git a/src/StratMetropolis.jl b/src/StratMetropolis.jl index f5260c9..3c2bec3 100644 --- a/src/StratMetropolis.jl +++ b/src/StratMetropolis.jl @@ -164,7 +164,7 @@ bounding = config.bounding # Stratigraphic age constraints - ages = unionize(smpl.Age_Distribution)::Vector{<:Union{<:Distribution{Univariate, Continuous}}} + ages = unionize(smpl.Age)::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 @@ -227,7 +227,7 @@ bounding = config.bounding # Stratigraphic age constraints - ages = unionize(smpl.Age_Distribution)::Vector{<:Union{<:Distribution{Univariate, Continuous}}} + ages = unionize(smpl.Age)::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 diff --git a/test/testStratOnlyGeneral.jl b/test/testStratOnlyGeneral.jl index 48276ce..25bc9f3 100644 --- a/test/testStratOnlyGeneral.jl +++ b/test/testStratOnlyGeneral.jl @@ -2,10 +2,10 @@ nsamples = 5 # Make an instance of a GeneralAgeData object for n samples smpl = GeneralAgeData(nsamples) -smpl.Name = ( "Sample 1", "Sample 2", "Sample 3", "Sample 4", "Sample 5",) # Et cetera -smpl.Age_Distribution = [Normal(39.5,0.1), Uniform(37, 38), Normal(36.3, 0.1), Uniform(33.5,34), Normal(32.1, 0.1),] # Measured ages -smpl.Height = [ 100, 200, 300, 400, 500,] # Depths below surface should be negative -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.Name = ( "Sample 1", "Sample 2", "Sample 3", "Sample 4", "Sample 5",) # Et cetera +smpl.Age = [Normal(39.5,0.1), Uniform(37, 38), Normal(36.3, 0.1), Uniform(33.5,34), Normal(32.1, 0.1),] # Measured ages +smpl.Height = [ 100, 200, 300, 400, 500,] # Depths below surface should be negative +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 = "Ma" # Unit of measurement for ages smpl.Height_Unit = "m" # Unit of measurement for Height and Height_sigma @test smpl isa GeneralAgeData