From 5cc1b28c2301f81c6a1ebfa27a3b6a47d064232f Mon Sep 17 00:00:00 2001 From: "C. Brenhin Keller" Date: Thu, 27 Jun 2024 14:06:55 -0400 Subject: [PATCH] Update examples, add tabular input options --- .gitignore | 2 + examples/Chron1.0Coupled.jl | 16 ++- examples/Chron1.0CoupledConcordia.jl | 13 +- examples/Chron1.0CoupledSystematic.jl | 12 +- examples/Chron1.0CoupledTabularInput.jl | 49 ++++--- examples/Chron1.0Radiocarbon.jl | 15 +- examples/Chron1.0StratOnly.jl | 16 ++- examples/Chron1.0StratOnlyTabularInput.jl | 46 +++--- examples/coupled_example_height_query.csv | 147 ++++++++++++++++++++ examples/stratonly_example_height_query.csv | 58 ++++++++ 10 files changed, 312 insertions(+), 62 deletions(-) create mode 100644 examples/coupled_example_height_query.csv create mode 100644 examples/stratonly_example_height_query.csv diff --git a/.gitignore b/.gitignore index f7f3d9d..7f7e849 100644 --- a/.gitignore +++ b/.gitignore @@ -18,6 +18,8 @@ Manifest.toml examples/DenverUPbExampleData/results.csv results.csv distresults.csv +AgeDepthModel.csv +stratonly_example_height_results.csv .ipynb_check* /etc diff --git a/examples/Chron1.0Coupled.jl b/examples/Chron1.0Coupled.jl index 7aad51c..c2b3a6b 100644 --- a/examples/Chron1.0Coupled.jl +++ b/examples/Chron1.0Coupled.jl @@ -111,7 +111,7 @@ # Run the stratigraphic MCMC model @time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config) - + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") ## --- Plot stratigraphic model - - - - - - - - - - - - - - - - - - - - - - - - @@ -178,7 +178,13 @@ dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile) # Plot results - hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2) + hdl = plot( + xlabel="Age ($(smpl.Age_Unit))", + ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", + fg_color_legend=:white, + framestyle=:box, + ) + plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2) plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI") for lci in 20:5:45 dhdt_lp = nanpctile(dhdt_dist,lci,dim=2) @@ -186,8 +192,7 @@ plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="") end plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1) - plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white) - # savefig(hdl,"DepositionRateModelCI.pdf") + savefig(hdl,"DepositionRateModelCI.pdf") display(hdl) ## --- Make heatmap - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -250,12 +255,15 @@ # Run the model. Note the additional `hiatus` arguments @time (mdl, agedist, hiatusdist, lldist) = StratMetropolisDist(smpl, hiatus, config) + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") # 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="") plot!(hdl, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label="data",seriestype=:scatter,color=:black) plot!(hdl, xlabel="Age (Ma)", ylabel="Height (cm)", framestyle=:box) + savefig(hdl,"AgeDepthModel.pdf") + display(hdl) ## --- (Optional) Add systematic uncertainties for U-Pb data diff --git a/examples/Chron1.0CoupledConcordia.jl b/examples/Chron1.0CoupledConcordia.jl index 9c7619a..6355725 100644 --- a/examples/Chron1.0CoupledConcordia.jl +++ b/examples/Chron1.0CoupledConcordia.jl @@ -110,7 +110,7 @@ # Run the stratigraphic MCMC model @time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config) - + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") ## --- Plot stratigraphic model - - - - - - - - - - - - - - - - - - - - - - - - @@ -178,7 +178,13 @@ dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile) # Plot results - hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2) + hdl = plot( + xlabel="Age ($(smpl.Age_Unit))", + ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", + fg_color_legend=:white, + framestyle=:box, + ) + plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2) plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI") for lci in 20:5:45 dhdt_lp = nanpctile(dhdt_dist,lci,dim=2) @@ -186,8 +192,7 @@ plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="") end plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1) - plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white) - # savefig(hdl,"DepositionRateModelCI.pdf") + savefig(hdl,"DepositionRateModelCI.pdf") display(hdl) ## --- End of File diff --git a/examples/Chron1.0CoupledSystematic.jl b/examples/Chron1.0CoupledSystematic.jl index 4e875e3..86e85a3 100644 --- a/examples/Chron1.0CoupledSystematic.jl +++ b/examples/Chron1.0CoupledSystematic.jl @@ -124,6 +124,7 @@ # Run the stratigraphic MCMC model @time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config, systematic) + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") # # Other options: # # Youngest Zircon @@ -213,7 +214,13 @@ dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile) # Plot results - hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2) + hdl = plot( + xlabel="Age ($(smpl.Age_Unit))", + ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", + fg_color_legend=:white, + framestyle=:box, + ) + plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2) plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI") for lci in 20:5:45 dhdt_lp = nanpctile(dhdt_dist,lci,dim=2) @@ -221,8 +228,7 @@ plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="") end plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1) - plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white) - # savefig(hdl,"DepositionRateModelCI.pdf") + savefig(hdl,"DepositionRateModelCI.pdf") display(hdl) ## --- Probability that a given interval of stratigraphy was deposited entirely before/after a given time diff --git a/examples/Chron1.0CoupledTabularInput.jl b/examples/Chron1.0CoupledTabularInput.jl index 2418798..e3b09da 100644 --- a/examples/Chron1.0CoupledTabularInput.jl +++ b/examples/Chron1.0CoupledTabularInput.jl @@ -12,18 +12,19 @@ using Chron using Plots + cd(@__DIR__) # Move to script directory ## --- Define sample properties - - - - - - - - - - - - - - - - - - - - - - - - # As in Chron1.0Coupled, but here we read in the information from a csv table: - ds = importdataset(joinpath(@__DIR__, "coupled_example_input.csv"), importas=:Tuple) + ds = importdataset("coupled_example_input.csv", importas=:Tuple) smpl = ChronAgeData(length(ds.Name)) smpl.Name = Tuple(ds.Name) smpl.Height .= ds.Height smpl.Height_sigma .= ds.Height_sigma smpl.Age_Sidedness .= ds.Age_Sidedness # 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 = joinpath(@__DIR__, first(ds.Path)) # Where are the data files? + smpl.Path = first(ds.Path) # Where are the data files? smpl.inputSigmaLevel = first(ds.inputSigmaLevel) # i.e., are the data files 1-sigma or 2-sigma. Integer. smpl.Age_Unit = first(ds.Age_Unit) # Unit of measurement for ages and errors in the data files smpl.Height_Unit = first(ds.Height_Unit) # Unit of measurement for Height and Height_sigma @@ -103,7 +104,7 @@ # Run the stratigraphic MCMC model @time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config) - + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") ## --- Plot stratigraphic model - - - - - - - - - - - - - - - - - - - - - - - - @@ -126,25 +127,19 @@ display(hdl) -## --- Interpolate model age at a specific stratigraphic height - - - - - - - - - - # Stratigraphic height at which to interpolate - interp_height = 0 +## --- Interpolate results at a specific heights read from file - - - - - - - - - age_at_height = linterp1s(mdl.Height,mdl.Age,interp_height) - age_at_height_min = linterp1s(mdl.Height,mdl.Age_025CI,interp_height) - age_at_height_max = linterp1s(mdl.Height,mdl.Age_975CI,interp_height) - print("Interpolated age at height=$interp_height: $age_at_height +$(age_at_height_max-age_at_height)/-$(age_at_height-age_at_height_min) $(smpl.Age_Unit)") + # Stratigraphic heights at which to interpolate + hds = importdataset("coupled_example_height_query.csv", importas=:Tuple) - # Optional: interpolate full age distribution - interpolated_distribution = Array{Float64}(undef,size(agedist,2)) - for i=1:size(agedist,2) - interpolated_distribution[i] = linterp1s(mdl.Height,agedist[:,i],interp_height) - end - hdl = histogram(interpolated_distribution, nbins=50, label="", framestyle=:box) - plot!(hdl, xlabel="Age ($(smpl.Age_Unit)) at height=$interp_height", ylabel="Likelihood (unnormalized)") - savefig(hdl, "Interpolated age distribution.pdf") - display(hdl) + hds = (; + Height = hds.Height, + Age = linterp1s(mdl.Height, mdl.Age, hds.Height), + Age_sigma = linterp1s(mdl.Height, mdl.Age_sigma, hds.Height), + Age_025CI = linterp1s(mdl.Height, mdl.Age_025CI, hds.Height), + Age_975CI = linterp1s(mdl.Height, mdl.Age_975CI, hds.Height), + ) + exportdataset(hds, "coupled_example_height_results.csv") ## --- Calculate deposition rate binned by age - - - - - - - - - - - - - - - - @@ -170,7 +165,13 @@ dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile) # Plot results - hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2) + hdl = plot( + xlabel="Age ($(smpl.Age_Unit))", + ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", + fg_color_legend=:white, + framestyle=:box, + ) + plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2) plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI") for lci in 20:5:45 dhdt_lp = nanpctile(dhdt_dist,lci,dim=2) @@ -178,8 +179,7 @@ plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="") end plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1) - plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white) - # savefig(hdl,"DepositionRateModelCI.pdf") + savefig(hdl,"DepositionRateModelCI.pdf") display(hdl) ## --- Make heatmap - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -216,12 +216,15 @@ # Run the model. Note the additional `hiatus` arguments @time (mdl, agedist, hiatusdist, lldist) = StratMetropolisDist(smpl, hiatus, config) + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") # 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="") plot!(hdl, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label="data",seriestype=:scatter,color=:black) plot!(hdl, xlabel="Age (Ma)", ylabel="Height (cm)", framestyle=:box) + savefig(hdl,"AgeDepthModel.pdf") + display(hdl) ## --- (Optional) Add systematic uncertainties for U-Pb data diff --git a/examples/Chron1.0Radiocarbon.jl b/examples/Chron1.0Radiocarbon.jl index a0fb10c..9bde75b 100644 --- a/examples/Chron1.0Radiocarbon.jl +++ b/examples/Chron1.0Radiocarbon.jl @@ -84,6 +84,7 @@ # Run the stratigraphic MCMC model @time (mdl, agedist, lldist) = StratMetropolis14C(smpl, config) + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") # Plot results (mean and 95% confidence interval for both model and data) hdl = plot(framestyle=:box, @@ -140,7 +141,13 @@ dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile) # Plot results - hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2) + hdl = plot( + xlabel="Age ($(smpl.Age_Unit))", + ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", + fg_color_legend=:white, + framestyle=:box, + ) + plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2) plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI") for lci in 20:5:45 dhdt_lp = nanpctile(dhdt_dist,lci,dim=2) @@ -148,8 +155,7 @@ plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="") end plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1) - plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white) - # savefig(hdl,"DepositionRateModelCI.pdf") + savefig(hdl,"DepositionRateModelCI.pdf") display(hdl) ## --- Optional: Stratigraphic model including hiatuses - - - - - - - - - - - - @@ -164,6 +170,7 @@ # Run the model. Note the additional `hiatus` arguments @time (mdl, agedist, hiatusdist, lldist) = StratMetropolis14C(smpl, hiatus, config) + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") # Plot results (mean and 95% confidence interval for both model and data) hdl = plot(framestyle=:box, @@ -173,5 +180,7 @@ plot!(hdl, [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, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label="data",seriestype=:scatter,color=:black) + savefig(hdl, "Interpolated age distribution.pdf") + display(hdl) ## --- End of File - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/Chron1.0StratOnly.jl b/examples/Chron1.0StratOnly.jl index 56c24da..bf71a77 100644 --- a/examples/Chron1.0StratOnly.jl +++ b/examples/Chron1.0StratOnly.jl @@ -50,6 +50,7 @@ # Run the stratigraphic MCMC model @time (mdl, agedist, lldist) = StratMetropolis(smpl, config) + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") # Plot results (mean and 95% CI for model / 2-sigma for data) hdl = plot(framestyle=:box, @@ -114,7 +115,13 @@ dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile) # Plot results - hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2) + hdl = plot( + xlabel="Age ($(smpl.Age_Unit))", + ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", + fg_color_legend=:white, + framestyle=:box, + ) + plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2) plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI") for lci in 20:5:45 dhdt_lp = nanpctile(dhdt_dist,lci,dim=2) @@ -122,10 +129,8 @@ plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="") end plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1) - plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white) - # savefig(hdl,"DepositionRateModelCI.pdf") + savefig(hdl,"DepositionRateModelCI.pdf") display(hdl) - ## --- Optional: Stratigraphic model including hiatuses - - - - - - - - - - - - # Data about hiatuses @@ -138,11 +143,14 @@ # Run the model. Note the additional `hiatus` arguments @time (mdl, agedist, hiatusdist, lldist) = StratMetropolis(smpl, hiatus, config); sleep(0.5) + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") # 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, smpl.Age, smpl.Height, xerror=smpl.Age_sigma*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) ## --- End of File - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/Chron1.0StratOnlyTabularInput.jl b/examples/Chron1.0StratOnlyTabularInput.jl index ea09e61..b8d61c3 100644 --- a/examples/Chron1.0StratOnlyTabularInput.jl +++ b/examples/Chron1.0StratOnlyTabularInput.jl @@ -9,13 +9,14 @@ using Chron using Plots + cd(@__DIR__) # Move to script directory ## --- Define sample properties - - - - - - - - - - - - - - - - - - - - - - - - # As in Chron1.0Coupled, but here we read in the information from a csv table: # # # # # # # # # # # Enter sample information here # # # # # # # # # # # # # - ds = importdataset(joinpath(@__DIR__, "stratonly_example_input.csv"), importas=:Tuple) + ds = importdataset("stratonly_example_input.csv", importas=:Tuple) smpl = ChronAgeData(length(ds.Name)) smpl.Name = Tuple(ds.Name) smpl.Age .= ds.Age @@ -50,6 +51,7 @@ # Run the stratigraphic MCMC model @time (mdl, agedist, lldist) = StratMetropolis(smpl, config) + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") # Plot results (mean and 95% CI for model / 2-sigma for data) hdl = plot(framestyle=:box, @@ -70,25 +72,19 @@ savefig(hdl,"AgeDepthModel.pdf") display(hdl) -## --- Interpolate results at a specific height - - - - - - - - - - - - - - - - +## --- Interpolate results at a specific heights read from file - - - - - - - - - # Stratigraphic height at which to interpolate - height = -404 + # Stratigraphic heights at which to interpolate + hds = importdataset("stratonly_example_height_query.csv", importas=:Tuple) - age_at_height = linterp1s(mdl.Height,mdl.Age,height) - age_at_height_min = linterp1s(mdl.Height,mdl.Age_025CI,height) - age_at_height_max = linterp1s(mdl.Height,mdl.Age_975CI,height) - print("Interpolated age at height=$height: $age_at_height +$(age_at_height_max-age_at_height)/-$(age_at_height-age_at_height_min) $(smpl.Age_Unit)") - - # Optional: interpolate full age distribution - interpolated_distribution = Array{Float64}(undef,size(agedist,2)) - for i=1:size(agedist,2) - interpolated_distribution[i] = linterp1s(mdl.Height,agedist[:,i],height) - end - hdl = histogram(interpolated_distribution, nbins=50, label="") - plot!(hdl, xlabel="Age ($(smpl.Age_Unit)) at height=$height", ylabel="Likelihood (unnormalized)") - savefig(hdl, "Interpolated age distribution.pdf") - display(hdl) + hds = (; + Height = hds.Height, + Age = linterp1s(mdl.Height, mdl.Age, hds.Height), + Age_sigma = linterp1s(mdl.Height, mdl.Age_sigma, hds.Height), + Age_025CI = linterp1s(mdl.Height, mdl.Age_025CI, hds.Height), + Age_975CI = linterp1s(mdl.Height, mdl.Age_975CI, hds.Height), + ) + exportdataset(hds, "stratonly_example_height_results.csv") ## --- Calculate deposition rate binned by age - - - - - - - - - - - - - - - - - @@ -114,7 +110,13 @@ dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile) # Plot results - hdl = plot(agebincenters,dhdt, label="Mean", color=:black, linewidth=2) + hdl = plot( + xlabel="Age ($(smpl.Age_Unit))", + ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", + fg_color_legend=:white, + framestyle=:box, + ) + plot!(hdl, agebincenters,dhdt, label="Mean", color=:black, linewidth=2) plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label="68% CI") for lci in 20:5:45 dhdt_lp = nanpctile(dhdt_dist,lci,dim=2) @@ -122,8 +124,7 @@ plot!(hdl,[agebincenters; reverse(agebincenters)],[dhdt_lp; reverse(dhdt_up)], fill=(0,0.2,:darkblue), linealpha=0, label="") end plot!(hdl, agebincenters,dhdt_50p, label="Median", color=:grey, linewidth=1) - plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Depositional Rate ($(smpl.Height_Unit) / $(smpl.Age_Unit) over $binwidth $(smpl.Age_Unit))", fg_color_legend=:white) - # savefig(hdl,"DepositionRateModelCI.pdf") + savefig(hdl,"DepositionRateModelCI.pdf") display(hdl) ## --- Optional: Stratigraphic model including hiatuses - - - - - - - - - - - - @@ -138,11 +139,14 @@ # Run the model. Note the additional `hiatus` arguments @time (mdl, agedist, hiatusdist, lldist) = StratMetropolis(smpl, hiatus, config); sleep(0.5) + exportdataset(NamedTuple(mdl), "AgeDepthModel.csv") # 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, smpl.Age, smpl.Height, xerror=smpl.Age_sigma*2,label="data",seriestype=:scatter,color=:black) plot!(hdl, xlabel="Age ($(smpl.Age_Unit))", ylabel="Height ($(smpl.Height_Unit))") + savefig(hdl,"AgeDepthModel.pdf") + display(hdl) ## --- End of File - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/examples/coupled_example_height_query.csv b/examples/coupled_example_height_query.csv new file mode 100644 index 0000000..75405d0 --- /dev/null +++ b/examples/coupled_example_height_query.csv @@ -0,0 +1,147 @@ +Height +-52 +-51 +-50 +-49 +-48 +-47 +-46 +-45 +-44 +-43 +-42 +-41 +-40 +-39 +-38 +-37 +-36 +-35 +-34 +-33 +-32 +-31 +-30 +-29 +-28 +-27 +-26 +-25 +-24 +-23 +-22 +-21 +-20 +-19 +-18 +-17 +-16 +-15 +-14 +-13 +-12 +-11 +-10 +-9 +-8 +-7 +-6 +-5 +-4 +-3 +-2 +-1 +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +32 +33 +34 +35 +36 +37 +38 +39 +40 +41 +42 +43 +44 +45 +46 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +58 +59 +60 +61 +62 +63 +64 +65 +66 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +81 +82 +83 +84 +85 +86 +87 +88 +89 +90 +91 +92 +93 diff --git a/examples/stratonly_example_height_query.csv b/examples/stratonly_example_height_query.csv new file mode 100644 index 0000000..7ad7c5f --- /dev/null +++ b/examples/stratonly_example_height_query.csv @@ -0,0 +1,58 @@ +Height +-355 +-356 +-357 +-358 +-359 +-360 +-361 +-362 +-363 +-364 +-365 +-366 +-367 +-368 +-369 +-370 +-371 +-372 +-373 +-374 +-375 +-376 +-377 +-378 +-379 +-380 +-381 +-382 +-383 +-384 +-385 +-386 +-387 +-388 +-389 +-390 +-391 +-392 +-393 +-394 +-395 +-396 +-397 +-398 +-399 +-400 +-401 +-402 +-403 +-404 +-405 +-406 +-407 +-408 +-409 +-410 +-411