README.md

BayeZirChron.c

C version of the Bayesian zircon eruption age estimation code from Keller, Schoene, and Samperton (2018)

Installation from command line

Installation from the command line requires a working C compiler. The default makefile assumes gcc (or an alias) is available. On linux/unix/bsd this is likely already true; on Mac OS the necessary tools for compiling C source can be installed by typing xcode-select --install at the command line.

# Download
git clone https://github.com/brenhinkeller/BayeZirChron.c.git

# Move to folder containing source code
cd BayeZirChron.c/src/

# Compile
make serial

To compile the parallel code used for synthetic distribution tests, additionally run:

make parallel

or

mpicc -std=c11 -O3 -o tzircrystestaccuracy tzircrystestaccuracy.c

Compiling and running this parallel version additionally requires a working installation of MPI (either Open MPI or MPICH) A sample batchfile is provided in the example folder: runTest.pbs

Usage

A range of examples, including the application Bayesian zircon eruption age estimation code to literature datasets, is provided. A Matlab script to run all literature examples is provided in examples/literature dataset tests/RunLiteratureExamples.m

Basic command-line usage follows the pattern:

tzircrystmetropolis <nsteps> distribution.tsv  sample.tsv > output.tsv

where is the length of Markov chain to run, distribution.tsv is an ascii file describing a zircon saturation distribution, sample.tsv is a tab or comma-separated file containing a list of individual zircon ages and 2-sigma uncertainties; the resulting stationary distributions will be written to output.tsv. For example

tzircrystmetropolis 10000 MeltsTZircDistribution.tsv  sample.tsv > metropolisdata.tsv

The results of synthetic dataset tests, which compare traditional weighted-mean and youngest-zircon interpretations to Bayesian eruption age estimates (using several different crystallization distributions) are provided in examples/synthetic dataset tests/. Figures can be re-plotted using a Matlab script tzircrystestaccuracyPlots.m.

To reproduce the datafiles provided in this folder, compile the parallel code tzircrystestaccuracy.c as described above (using mpicc), and run using

mpiexec -np <number-of-tasks> ./tzircrystestaccuracy  <sims-per-task> <nsteps> <dt/sigma>  Distribution.tsv  > results.tsv

where is the number of MPI tasks to run (typically you want this to be equal to the number of CPU cores or hardware threads you are running on), is the number of simulations (at each N) to run per MPI task, is the length of Markov chain to run, <dt/sigma> is the crystallization timescale in units of sigma (analytical uncertainty), pulling synthetic data from a distribution specified in an ascii file `Distribution.tsv for example:

mpiexec -np 16 ./tzircrystestaccuracy 4 10000 1 MeltsTZircDistribution.tsv > eruptionestimates1.tsv

to run on 16 cores with 4 simulations per taks, each 10000 MCMC steps long, with a dt/sigma of 1 and using the crystallization distribution found in MeltsTZircDistribution.tsv. To run on a cluster, you will need a batch file suited to your cluster's workload manager. For instance, the example batchfile runTest.pbs:

#!/bin/bash
#PBS -l nodes=20:ppn=16,walltime=00:40:00

module load openmpi
mpiexec ./tzircrystestaccuracy 4 10000 1 MeltsTZircDistribution.tsv > eruptionestimates1.tsv

runs 4 simulations per task on 20 nodes with 16 cores each, for a total of 1280 simulations at each N

Julia version

See julia folder. Try it out in this interactive Jupyter notebook via Binder!

For future work, see Chron.jl, which can be run in standalone eruption/deposition mode using the approach developed here.