GNS 1.0
Geometric nested sampling algorithm, as described in https://arxiv.org/abs/1905.09110.
License: MIT
Requires Python: >=2.7.0

# geometric nested sampling

## Version

Version: 1.0

## Overview

Nested sampling algorithms are a family of algorithms designed to sample from posterior probability distributions associated with e.g. Bayesian inference. The overall idea in nested sampling is to convert an n-dimensional integral over the parameter space into a single dimensional integral over 'prior volume' and equivalently sample from this domain. The purpose of the geometric nested sampling algorithm presented here is to efficiently sample the distributions of parameters which exhibit certain geometrical properties e.g. circular, toroidal or spherical. Such parameters are rife in the astrophysics community, where many time- or angle-periodic parameters are present.

## The algorithm

Can run Markov Chain Monte Carlo (MCMC) Metropolis Hastings (MH) nested sampling algorithm in 'vanilla' (see https://projecteuclid.org/euclid.ba/1340370944 and https://arxiv.org/abs/0704.3704) and/or 'geometric' (https://arxiv.org/abs/1905.09110) sampling modes. See the appropriate papers for detailed explanations of how these algorithms operate.

In vanilla mode, the trial points are sampled in the physical (native) space of the input parameters.
In geometric mode (main point of the algorithm), parameter values represent points on a geometric object (e.g. circle, torus or sphere), and trial points are picked by sampling from spaces parameterised by Cartesian coordinates. If the sampled points do not lay on the original object (in general they will not), then they are projected back onto the object before being transformed back to the native space for evaluation of the MH acceptance ratio.

The sampling mode for each parameter is specified in `paramGeomList` passed through `setupDict` to `nestedRun()` ("vanilla", "wrapped" or "sphere"). "wrapped" parameters use a wrapped trial distribution. "sphere" parameters are transformed into their Euclidean representation and sampled in this space. Geometric sampling works best for parameters with 'periodic' sampling domains e.g. angles, as sampling in geometric space results in wrapping parameter values at end of domains. The sampler also works well for parameters which exhibit "spherical" properties. Note for sampling spherical sampling, two physical parameters are required for each geometric transformation. See `tests.py` for examples of how to specify what parameters are sampled in what way and the below section showing `paramGeomList` examples. 

Main output of geometric nested sampling run is file of sampled points and their weights. These are stored in `text_output/file_name.txt`, whose format follows that explained in https://getdist.readthedocs.io/en/latest/intro.html. Similarly, metadata files are also generated for each run e.g. `file_name.paramnames`, see https://getdist.readthedocs.io/en/latest/intro.html for more details on these files. The format of the text outputs is such that they can be plotted very easily using `getdist`, a kernel density esimate plotter.

## Toy models

To run with toy models (defined and explained in `toy_models.py`), go to `tests/` and run `tests.py`, which calls the function `runTests()`, which can be found in the `tests.py` file. Plots for comparison can be found in `image_output/toy_models_empirical_getdist/`. Text output can be found in `text_output/`. The default toy model used when `main.py` is ran as is when cloned from this repo, should be compared with files `image_output/toy_models_empirical_getdist/circle.png` (in particular, the MHWG curve on this plot). Files created from running main.py for this toy example are saved as th_Sc_mhwg* files in `text_output/` and `image_output/toy_models_empirical_getdist`. Further toy models can be run by adding them to `tests.py` and adjusting the necessary parameters (see the commented out example for how to configure the parameters). The following plot was produced from running the default circle toy model, and shows the kernel density esimate of samples obtained by the geometric nested sampler from the von Mises distribution.


![alt text](gns_circle.png)

Examples of `setupDict` objects used to specify parameters of a nested sampling run can be found in `setup_dicts.py`. 'verbose' corresponds to whether to output text during the sampling run, 'trapezoidalFlag' specifies whether to use the trapezoidal rectangular quadrature techniques for approximating the integral, 'ZLiveType' specifies how to determine the final contribution to the nested sample run, and shouldn't be significant if the algorithm is ran for a decent amount of time, 'terminationType' is a similar parameter, as is 'terminationFactor', 'sampler' is the type of sampler used i.e. Blind nested sampling, which involves randomly sampling values from the unit hypercube (B) Metropolis Hastings nested sampling (MH) geometric Metropolis Hastings nested sampling (MHWG), or use MultiNest if it is installed (MN). 'outputFile': is the location to store the text output from the runs, 'space' corresponds to whether the likelihood or log likelihood is used in calculations (and thus the corresponding function should be passed to the nested sampling algorithm through the appropriate arguments, see `nestedRun()`). 'paramGeomList': is a list of size of the dimensionality of the problem, and specifies how each parameter should be treated by the geometric nested sampler: 'vanilla', 'wrapped' or 'sphere' (must also be specified for the MH algorithm, which only accepts 'vanilla' values) Note for the toy models, none of these really need to be changed from their default in `tests.py` apart from `paramGeomList`, only the models chosen and the prefixes/suffixes used to index them need to be altered.

The toy models (see https://arxiv.org/abs/1905.09110) are designed to make use of the key features of geometric nested sampling, i.e. to take advantage of the way it treats parameters defined on curved manifolds. The three toy models which feature in `tests.py` by default are a circular distribution (von Mises distribution), a toroidal distribution (two-dimensional von Mises distribution), and 6 spheres with Kent distributions defined on them (for more on this see https://arxiv.org/abs/1905.09110 and https://www.jstor.org/stable/2984712?seq=1#page_scan_tab_contents). Many more are available in `toy_models.py`, all that needs to be configured is there `paramGeomList` values when added to `tests.py`, and the relevant suffixes/prefixes to index them.

## Running the algorithm with custom likelihood/prior functions

If custom priors/ likelihoods are to be used, one needs to call the `nestedRun()` function in the `gns.nested_run` module, and pass the required likelihood and prior functions. All of these functions must take as an argument a `numpy` array with shape `(nLive, nDims)` where `nLive` is the number of live points in the nested run, and is defined in the functions `NestedRunLinear` and `NestedRunLog` which are in the `nested_run.py` file. `nDims` is the dimensionality of the parameter space of the inference. Furthermore, these functions should return an array of size `(nLive)`, with the respective likelihood/prior values for each live point. Note the algorithm can be run in linear or log mode (by running `NestedRunLinear` or `NestedRunLog`, which is determined by the corresponding value in the `setupDict`), the latter of which requires a function which calculates the log likelihood. See `prob_funcs.py` for examples of such functions. Sampler will also work with probability distributions (and their inverse CDFs) defined in `scipy.stats`. 

Note that because the nested sampler samples in the unit hypercube, functions for the inverse of the CDF of the prior distributions are also required. This function also takes an array of size `(nLive, nDims)`, but returns an array of the same shape. Note also that the unit hypercube and physical parameter vectors take the shape `(nLive, nDims)`. 

Note that the geometric nested sampling algorithm has many tunable parameters, namely `nLive`, `maxTrials`, the max number of points to trial in each Markov chain, and the standard distributions of the trial distributions used in the Markov chains, `trialSigma` for `wrapped` parameters which is found in `samplers.py`, and `sigmaX`, `sigmaY`, `sigmaZ` for `sphere` parameters found in `geom_sampler.py`. The 'optimum' values of these parameters is entirely domain specific (though the sampler will generally be more accurate the higher `nLive` is), but refer to https://arxiv.org/abs/1905.09110 and https://www.repository.cam.ac.uk/handle/1810/293473 for more information on these parameters.

## `paramGeomList` examples

### Example 1

Consider a problem with `nLive = 5`, with the first two parameters exhibiting circular properties (i.e. a torus), the third being a vanilla paramater, and the last pair forming a pair of sperical parameters. Then `paramGeomList = ['wrapped', 'wrapped', 'vanilla', 'sphere', 'sphere']`

### Example 2

Different types of parameters don't have to be in sequential order. The example `paramGeomList = ['vanilla', 'wrapped', 'wrapped', 'wrapped', 'sphere', 'sphere', 'sphere', 'sphere', 'wrapped']` is perfectly valid and corresponds to the first parameter being sampled in vanilla model, the next three parameters being circular (3-torus), the following four being spherical, and thus defining the points on two spheres (with the first and third being the azimuthal angle of each sphere and the second and fourth being the polar angles), and the final parameter being circular.

In general the parameters can be in an arbitrary order, the only requirements are that the order of the list corresponds to the variable's position in the `(nLive, nDims)` arrays passed to the algorithm, and that spherical parameters come in (sequential) pairs.

## Additional features

Functions to calculate theoretical Bayesian evidence & H (K-L divergence) values are also provided, but not recommended in higher dimensional problems.

There is also a function to calculate the analytical posterior across a grid of the sampling space (or random samples), which produces files which can be input into getdist to be used for comparison with empirical results (assigns weights equal to Lhood or posterior for use in `getdist`). 

Plotting functions are also provided, in the form of simple `matplotlib` functions, as well as `getdist` (https://pypi.python.org/pypi/GetDist/) plotting functions.

Can also call MultiNest (https://arxiv.org/abs/0809.3437) by calling the pyMultiNest wrapper (https://github.com/JohannesBuchner/PyMultiNest) for comparison. 

These are automatically employed in `tests.py`, but can be turned on and off by setting the appropriate flags in the file. 

## Installation

### Pip install

To install using pip, do 

`pip install GNS==1.0,` 

to install the gns algorithm and the requirements to run the algorithm. To also install the packages required for plotting, do 

`pip install GNS[plotting]==1.0.`

Finally, to install the gns package and include `MultiNest` algorithm runs (for comparison tests), do

`pip install GNS[multinest]==1.0.`

Note that the above do not install the tests found in the `tests/` directory

### git

No installation of the geometric nested sampler itself is strictly necessary, one can simply clone the gns repo:

`git clone https://github.com/SuperKam91/gns.git.`

Next, install all the modules included in the `requirements.txt` file to be able to run the geometric nested sampling algorithm. To use the plotting functions featured in the package, one must also install all the files included in `requirements_plotting.txt`. I highly recommend installing the plot-dependent packages, to compare results the results of the toy models with those in `image_output`. Finally, to run MultiNest, install all packages in `requirements_mutlinest.txt`.

Alternatively, one can run `python setup.py install` or `pip install .` from the root of the repo to install the base requirements (i.e. those in `requirements.txt`. As with the non-local `pip install` case, one can install the packages required for plotting and Multinest by doing or `pip install .[plotting]` (and similar for Multinest).

Note that by default, the tests in the `tests/` directory run the gns algorithm from the code in the repo in which `tests/` sits. To change this i.e. to use the installed package, in `tests.py`, set `run_repo = False`.