Documentation for the AH-RJMCMC code (Age-Hyperparameter Reverse-jump Markov Chain Monte Carlo)

Version1: March 27th 2018

Authors: 
Phil Livermore (University of Leeds)
Alex Fournier (Institut de Physique du Globe de Paris, Paris)
Thomas Bodin (Universite Lyon, Lyon)

**********
Overview
**********
This directory contains all input files, datasets and plotting scripts necessary to reproduce all figures in the manuscript submitted to Geophysical Journal International (March 2018).

There are two versions of the code, one written in Python, and one written in Fortran.
Each version of the code uses a common input file, and writes the same format of output. They will differ however in the random sampling of the posterior distribution, so the outputs will not be identical. However, provided the distributions are properly sampled, they will both converge and graphically they will be indistinguishable.


*********************
Required programmes
*********************
The AH-RJMCMC code either requires Python (tested only in version 3) or Fortran; Python is used for plotting.

*****************************
Basic execution in Linux/Mac
*****************************

To run the code using the Fortran version, type

./Make_all_figures Fortran

or, to run using the Python version, type

./Make_all_figures Python


If running in Fortran, the Makefile is sourced to "make" the Fortran version of the code. By default this uses the Intel Fortran compiler ifort, but it is straightforward to edit line one to use a different compiler (the gfortran compiler option is given, but commented out).

****************************
Using an individual inputfile
****************************
If not already compiled, you can compile the Fortran code by typing

make all

(see the notes above about using a preferred compiler).

To run the Fortran code on a single inputfile, type

AH inputfile

or (using the Python version)

python AH_RJMCMC.py inputfile


********************
Inputfile structure
********************
Each inputfile contains all the information needed to run the AH-RJMCMC (or IL-RJMCMC) model.

Comment lines begin with a '#' symbol and are ignored.

Each line of information begins with a keyword (in either upper or lower case) followed by the relevant information.

There are several input files in the Github directory that you can use as templates.

Keywords:

Data_file - The relative path of the dataset to be read in 

File_format - Columns of the datafile to be read in, in the order:  age, delta age, intensity, delta intensity, stratification
A value of stratification (in the datafile) of means that the datum is not constrained by stratification; a value of 1 means that it is.
The number of the columns begins at 0.

Note: if stratification is not relevant, set the stratification parameter to be -1: no stratification column is then read in (and so the data is assumed unstratified) 

A typical structure of File_format is:
File_format 2 3 4 5 -1   (i.e. the age is in the 3rd column - since the ordering begins at 0).

Burn_in - The number of samples for the burn-in period

Nsamples - The number of samples overall (including burn-in)

model_discretisation - the resolution in time of the returned posterior distribution
 
Chain_parameters - show, thinning
show: the frequency of writing information to the screen
thinning: the frequency of keeping samples in the Markov chain
 
Age distribution -  U = uniform, N = normal

running mode - switch between posterior and prior sampling. 1 - normal, 0 - set all likelihoods to 1 to recover the prior distributions.

Age_bounds - Age Min and Max for model parametrisation

sigma_move, sigma_change, sigma_birth - Parameters that describe the various model perturbations

Age_frac - (parameter beta in the manuscript): fraction of ages to change in a proposal

Intensity_prior - min/max bounds on intensity prior for vertices: I_min, I_max in micro Tesla

Num_change_points - K_min, K_max 
Number of internal vertices is bounded between [K_min, K_max]

Outputs_directory - Directory for all outputs
This directory is created if it does not exist

Credible - Credible interval or any non-positive number for none

Nbins - Number of bins for posterior marginals

output_model - Name and write frequency of model file
For example, the line below causes the code to write every 10th model (after thinning) to the models file here "models.dat". 

output_model models.dat 10

Enter a frequency of -1 for no model output.

output_joint_distribution_freq  - Joint distribution output frequency
If frequency is -1 then nothing is output; otherwise this defines the frequency of the joint distribution after thinning.

Optional parameters:
True_data - relative path of "true" underlying evolution
If this is present, some of the plotting code will add this data to figures.

# 
Plotting_intensity_range - min/max of plotting range of intensity
Used only by the plotting scripts. e.g.
Plotting_intensity_range 50 80