![\"Launch](\"https://static.mybinder.org/badge_logo.svg\"){kind=icon}
\n",
"If running this notebook as an online Binder notebook and the webpage times out, click the badge at left to relaunch (refreshing will not work). Note that any changes will be lost!
\n", "\n", "Hint: `shift`-`enter` to run a single cell, or from the `Cell` menu select `Run All` to run the whole file. Any code from this notebook can also be copied and pasted into the Julia REPL or a `.jl` script.\n", "***\n", "\n", "## Load required Julia packages" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Load (and install if necessary) the Chron.jl package\n", "try\n", " using Chron\n", "catch\n", " using Pkg\n", " Pkg.add(\"Chron\")\n", " using Chron\n", "end\n", "\n", "using Statistics, StatsBase, DelimitedFiles, SpecialFunctions\n", "using Plots; gr(); default(fmt = :png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "## Enter sample information\n", "Paste your data in here!" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Input the number of samples we wish to model (must match below)\n", "nSamples = 4\n", "\n", "# Make an instance of a ChronSection object for nSamples\n", "smpl = NewChronAgeData(nSamples)\n", "smpl.Name = (\"Sample 1\", \"Sample 2\", \"Sample 3\", \"Sample 4\") # Et cetera\n", "smpl.Age_14C = [ 6991, 7088, 7230, 7540,] # Measured ages\n", "smpl.Age_14C_sigma = [ 30, 70, 50, 50,] # Measured 1-σ uncertainties\n", "smpl.Height[:] = [ -355, -380,-397.0,-411.5,] # Depths below surface should be negative\n", "smpl.Height_sigma[:] = fill(0.01, nSamples) # Usually assume little or no sample height uncertainty\n", "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\n", "smpl.Path = \"MyData/\" # Where do you want output files to be stored\n", "\n", "smpl.Age_Unit = \"Years BP\"; # Unit of measurement for ages\n", "smpl.Height_Unit = \"m\"; # Unit of measurement for Height and Height_sigma" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that smpl.Height *must* increase with increasing stratigraphic height -- i.e., stratigraphically younger samples must be more positive. For this reason, it is convenient to represent depths below surface as negative numbers.\n", "\n", "***\n", "## Calculate calendar age PDFs for each sample\n", "Using the [IntCal 2013](https://doi.org/10.2458/azu_js_rc.55.16947) radiocarbon calibration curve" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "smpl.Params = fill(NaN, length(intcal13[\"Age_Calendar\"]), nSamples)\n", "hdl = plot(layout=nSamples)\n", "for i = 1:nSamples\n", " # The likelihood that a measured 14C age could result from a sample of\n", " # a given calendar age is proportional to the intergral of the product\n", " # of the two respective distributions\n", " likelihood = normproduct.(smpl.Age_14C[i], smpl.Age_14C_sigma[i], intcal13[\"Age_14C\"], intcal13[\"Age_sigma\"])\n", " likelihood ./= sum(likelihood) * intcal13[\"dt\"] # Normalize\n", "\n", " # Estimate mean and standard deviation by drawing samples from distribution\n", " samples = draw_from_distribution(likelihood, 10^6) .* maximum(intcal13[\"Age_Calendar\"])\n", " smpl.Age[i] = mean(samples)\n", " smpl.Age_sigma[i] = std(samples)\n", " smpl.Age_025CI[i] = percentile(samples,2.5)\n", " smpl.Age_975CI[i] = percentile(samples,97.5)\n", "\n", " # Populate smpl.Params with log likelihood for each sample\n", " smpl.Params[:,i] = normproduct_ll.(smpl.Age_14C[i], smpl.Age_14C_sigma[i], intcal13[\"Age_14C\"], intcal13[\"Age_sigma\"])\n", "\n", " # Plot likelihood vector for each sample\n", " t = (intcal13[\"Age_Calendar\"] .> smpl.Age[i] - 5*smpl.Age_sigma[i]) .& (intcal13[\"Age_Calendar\"] .< smpl.Age[i] + 5*smpl.Age_sigma[i])\n", " plot!(hdl[i], intcal13[\"Age_Calendar\"][t], likelihood[t], label=\"\", xlabel=\"Calendar Age\", ylabel=\"Likelihood\", title=\"$(smpl.Name[i])\")\n", "end\n", "display(hdl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "## Configure and run stratigraphic model\n", "note: to spare Binder's servers, this demo uses\n", "```\n", "config.nsteps = 3000 \n", "config.burnin = 2000*npoints_approx\n", "```\n", "However, you probably want higher numbers for a publication-quality result, for instance \n", "```\n", "config.nsteps = 30000 # Number of steps to run in distribution MCMC\n", "config.burnin = 10000*npoints_approx # Number to discard\n", "```\n", "and examine the log likelihood plot to make sure you've converged.\n", "\n", "To run the stratigraphic MCMC model, we call the `StratMetropolis14C` function. If you're not using radiocarbon ages and want to run stratigraphic model with Gaussian mean age and standard deviation for some number of stratigraphic horizons, then you can set `smpl.Age` and `smpl.Age_sigma` directly, but then you'll need to call `StratMetropolis` instead of `StratMetropolis14C`" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Generating stratigraphic age-depth model...\n", "Burn-in: 680000 steps\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32mBurn-in...100%|█████████████████████████████████████████| Time: 0:12:54\u001b[39m\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Collecting sieved stationary distribution: 1020000 steps\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32mCollecting...100%|██████████████████████████████████████| Time: 0:18:54\u001b[39m\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Configure the stratigraphic Monte Carlo model\n", "config = NewStratAgeModelConfiguration()\n", "# If in doubt, you can probably leave these parameters as-is\n", "config.resolution = 0.2 # Same units as sample height. Smaller is slower!\n", "config.bounding = 0.5 # how far away do we place runaway bounds, as a fraction of total section height\n", "(bottom, top) = extrema(smpl.Height)\n", "npoints_approx = round(Int,length(bottom:config.resolution:top) * (1 + 2*config.bounding))\n", "config.nsteps = 3000 # Number of steps to run in distribution MCMC\n", "config.burnin = 2000*npoints_approx # Number to discard\n", "config.sieve = round(Int,npoints_approx) # Record one out of every nsieve steps\n", "\n", "# Run the stratigraphic MCMC model\n", "(mdl, agedist, lldist) = StratMetropolis14C(smpl, config)\n", "\n", "# Write the results to file\n", "run(`mkdir -p $(smpl.Path)`) # Make sure that the path exists\n", "writedlm(smpl.Path*\"agedist.csv\", agedist, ',') # Stationary distribution of the age-depth model\n", "writedlm(smpl.Path*\"height.csv\", mdl.Height, ',') # Stratigraphic heights corresponding to each row of agedist\n", "writedlm(smpl.Path*\"lldist.csv\", lldist, ',') # Log likelihood distribution (to check for stationarity)\n", "\n", "# Plot the log likelihood to make sure we're converged (n.b burnin isn't recorded)\n", "hdl = plot(lldist,xlabel=\"Step number\",ylabel=\"Log likelihood\",label=\"\",line=(0.85,:darkblue))\n", "savefig(hdl,smpl.Path*\"lldist.pdf\")\n", "display(hdl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "## Plot results" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot results (mean and 95% confidence interval for both model and data)\n", "hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label=\"model\")\n", "plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label=\"\")\n", "plot!(hdl, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label=\"data\",seriestype=:scatter,color=:black)\n", "plot!(hdl, xlabel=\"Age ($(smpl.Age_Unit))\", ylabel=\"Height ($(smpl.Height_Unit))\")\n", "savefig(hdl,smpl.Path*\"AgeDepthModel.svg\")\n", "savefig(hdl,smpl.Path*\"AgeDepthModel.pdf\")\n", "display(hdl)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Interpolated age at height=-404: 8213.878377486633 +145.27607341744988/-152.82619915759005 Years BP" ] }, { "data": { "image/svg+xml": [ "\n", "\n" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Stratigraphic height at which to interpolate\n", "height = -404\n", "\n", "age_at_height = linterp1s(mdl.Height,mdl.Age,height)\n", "age_at_height_min = linterp1s(mdl.Height,mdl.Age_025CI,height)\n", "age_at_height_max = linterp1s(mdl.Height,mdl.Age_975CI,height)\n", "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)\")\n", "\n", "# Optional: interpolate full age distribution\n", "interpolated_distribution = Array{Float64}(undef,size(agedist,2))\n", "for i=1:size(agedist,2)\n", " interpolated_distribution[i] = linterp1s(mdl.Height,agedist[:,i],height)\n", "end\n", "hdl = histogram(interpolated_distribution, nbins=50, fill=(0.75,:darkblue), label=\"\")\n", "plot!(hdl, xlabel=\"Age ($(smpl.Age_Unit)) at height=$height\", ylabel=\"Likelihood (unnormalized)\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are other things we can plot as well, such as deposition rate:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ " 0.197671 seconds (600.74 k allocations: 54.665 MiB, 11.02% gc time)\n" ] } ], "source": [ "# Set bin width and spacing\n", "binwidth = round(nanrange(mdl.Age)/10,sigdigits=1) # Can also set manually, commented out below\n", "# binwidth = 100 # Same units as smpl.Age\n", "binoverlap = 10\n", "ages = collect(minimum(mdl.Age):binwidth/binoverlap:maximum(mdl.Age))\n", "bincenters = ages[1+Int(binoverlap/2):end-Int(binoverlap/2)]\n", "spacing = binoverlap\n", "\n", "# Calculate rates for the stratigraphy of each markov chain step\n", "dhdt_dist = Array{Float64}(undef, length(ages)-binoverlap, config.nsteps)\n", "@time for i=1:config.nsteps\n", " heights = linterp1(reverse(agedist[:,i]), reverse(mdl.Height), ages)\n", " dhdt_dist[:,i] .= abs.(heights[1:end-spacing] - heights[spacing+1:end]) ./ binwidth\n", "end\n", "\n", "# Find mean and 1-sigma (68%) CI\n", "dhdt = nanmean(dhdt_dist,dim=2)\n", "dhdt_50p = nanmedian(dhdt_dist,dim=2)\n", "dhdt_16p = nanpctile(dhdt_dist,15.865,dim=2) # Minus 1-sigma (15.865th percentile)\n", "dhdt_84p = nanpctile(dhdt_dist,84.135,dim=2) # Plus 1-sigma (84.135th percentile)\n", "# Other confidence intervals (10:10:50)\n", "dhdt_20p = nanpctile(dhdt_dist,20,dim=2)\n", "dhdt_80p = nanpctile(dhdt_dist,80,dim=2)\n", "dhdt_25p = nanpctile(dhdt_dist,25,dim=2)\n", "dhdt_75p = nanpctile(dhdt_dist,75,dim=2)\n", "dhdt_30p = nanpctile(dhdt_dist,30,dim=2)\n", "dhdt_70p = nanpctile(dhdt_dist,70,dim=2)\n", "dhdt_35p = nanpctile(dhdt_dist,35,dim=2)\n", "dhdt_65p = nanpctile(dhdt_dist,65,dim=2)\n", "dhdt_40p = nanpctile(dhdt_dist,40,dim=2)\n", "dhdt_60p = nanpctile(dhdt_dist,60,dim=2)\n", "dhdt_45p = nanpctile(dhdt_dist,45,dim=2)\n", "dhdt_55p = nanpctile(dhdt_dist,55,dim=2)\n", "\n", "# Plot results\n", "hdl = plot(bincenters,dhdt, label=\"Mean\", color=:black, linewidth=2)\n", "plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_16p; reverse(dhdt_84p)], fill=(0,0.2,:darkblue), linealpha=0, label=\"68% CI\")\n", "plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_20p; reverse(dhdt_80p)], fill=(0,0.2,:darkblue), linealpha=0, label=\"\")\n", "plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_25p; reverse(dhdt_75p)], fill=(0,0.2,:darkblue), linealpha=0, label=\"\")\n", "plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_30p; reverse(dhdt_70p)], fill=(0,0.2,:darkblue), linealpha=0, label=\"\")\n", "plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_35p; reverse(dhdt_65p)], fill=(0,0.2,:darkblue), linealpha=0, label=\"\")\n", "plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_40p; reverse(dhdt_60p)], fill=(0,0.2,:darkblue), linealpha=0, label=\"\")\n", "plot!(hdl,[bincenters; reverse(bincenters)],[dhdt_45p; reverse(dhdt_55p)], fill=(0,0.2,:darkblue), linealpha=0, label=\"\")\n", "plot!(hdl,bincenters,dhdt_50p, label=\"Median\", color=:grey, linewidth=1)\n", "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)\n", "savefig(hdl,smpl.Path*\"DepositionRateModelCI.pdf\")\n", "display(hdl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "## Stratigraphic model including hiatuses\n", "We can also deal with discrete hiatuses in the stratigraphic section if we know where they are and about how long they lasted. We need some different models and methods though. In particular, in addition to the `ChronAgeData` struct, we also need a `HiatusData` struct now, and will have a new `hiatusdist` output." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Generating stratigraphic age-depth model...\n", "Burn-in: 680000 steps\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32mBurn-in...100%|█████████████████████████████████████████| Time: 0:13:39\u001b[39m\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Collecting sieved stationary distribution: 1020000 steps\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[32mCollecting...100%|██████████████████████████████████████| Time: 0:21:35\u001b[39m\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "config.nsteps = 3000 # Number of steps to run in distribution MCMC\n", "config.burnin = 2000*npoints_approx # Number to discard\n", "# Data about hiatuses\n", "nHiatuses = 2 # The number of hiatuses you have data for\n", "hiatus = NewHiatusData(nHiatuses) # Struct to hold data\n", "hiatus.Height = [-371.5, -405.0 ]\n", "hiatus.Height_sigma = [ 0.0, 0.0 ]\n", "hiatus.Duration = [ 100.0, 123.0]\n", "hiatus.Duration_sigma = [ 30.5, 20.0]\n", "\n", "# Run the model. Note the new `hiatus` argument and `hiatusdist` result\n", "(mdl, agedist, hiatusdist, lldist) = StratMetropolis14C(smpl, hiatus, config); sleep(0.5)\n", "\n", "# Plot results (mean and 95% confidence interval for both model and data)\n", "hdl = plot([mdl.Age_025CI; reverse(mdl.Age_975CI)],[mdl.Height; reverse(mdl.Height)], fill=(minimum(mdl.Height),0.5,:blue), label=\"model\")\n", "plot!(hdl, mdl.Age, mdl.Height, linecolor=:blue, label=\"\", fg_color_legend=:white)\n", "plot!(hdl, smpl.Age, smpl.Height, xerror=(smpl.Age-smpl.Age_025CI,smpl.Age_975CI-smpl.Age),label=\"data\",seriestype=:scatter,color=:black)\n", "plot!(hdl, xlabel=\"Age ($(smpl.Age_Unit))\", ylabel=\"Height ($(smpl.Height_Unit))\")\n", "# savefig(hdl,smpl.Path*\"AgeDepthModelHiatus.pdf\");\n", "display(hdl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "## Getting your data out\n", "As shown below, we have access to the system command line, so we can use the unix command `ls` to see all the files we have written. We'll try this first using `;` to access Julia's command-line shell mode:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AgeDepthModel.pdf\n", "AgeDepthModel.svg\n", "AgeDepthModelHiatus.pdf\n", "DepositionRateModelCI.pdf\n", "agedist.csv\n", "height.csv\n", "lldist.csv\n", "lldist.pdf\n" ] } ], "source": [ "; ls MyData" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively, we could do the same thing using the `run` function:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "AgeDepthModel.pdf\n", "AgeDepthModel.svg\n", "AgeDepthModelHiatus.pdf\n", "DepositionRateModelCI.pdf\n", "agedist.csv\n", "height.csv\n", "lldist.csv\n", "lldist.pdf\n" ] }, { "data": { "text/plain": [ "Process(`\u001b[4mls\u001b[24m \u001b[4mMyData\u001b[24m`, ProcessExited(0))" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "run(`ls MyData`)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you're running this as an online notebook, getting these files out of here may be harder though. For SVG files, one option is to view them in markdown cells, which you should then be able to right click and download as real vector graphics. e.g. pasting something like\n", "```\n", "![](\"AgeDepthModel.svg\")
\n",
"```\n",
"in a markdown cell such as this one"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Meanwhile, for the csv files we could try something like `; cat agedist.csv`, but agedist is probably too big to print. Let's try using ffsend instead, which should give you a download link. In fact, while we're at it, we might as well archive and zip the entire directory!"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# Make gzipped tar archive of the the whole MyData directory\n",
"run(`tar -zcf archive.tar.gz ./MyData`);"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
" % Total % Received % Xferd Average Speed Time Time Time Current\n",
" Dload Upload Total Spent Left Speed\n",
"100 622 0 622 0 0 1718 0 --:--:-- --:--:-- --:--:-- 1713\n",
"100 10.8M 100 10.8M 0 0 3750k 0 0:00:02 0:00:02 --:--:-- 6418k\n"
]
}
],
"source": [
"# Download prebuilt ffsend linux binary\n",
"download(\"https://github.com/timvisee/ffsend/releases/download/v0.2.65/ffsend-v0.2.65-linux-x64-static\", \"ffsend\")\n",
"\n",
"# Make ffsend executable\n",
"run(`chmod +x ffsend`);"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"; ./ffsend upload archive.tar.gz"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You could alternatively use the ffsend command in this way to transfer individual files, for instance `; ./ffsend upload MyData/agedist.csv`"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Keep in mind that any changes you make to this online notebook won't persist after you close the tab (or after it times out) even if you save your changes! You have to either copy-paste or `file`>`Download as` a copy.\n",
"***\n",
"[](https://doi.org/10.17605/osf.io/TQX3F) "
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"kernelspec": {
"display_name": "Julia 1.5.0-rc1",
"language": "julia",
"name": "julia-1.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}