Getting started

After installing kima, we'll go through one of the examples in order to start playing with the package.

Let's find 51 Peg b!
In the kima/examples/51Peg folder you will find a dataset of radial-velocity observations and the code to analyse it with kima.

51 Peg b was the first exoplanet discovered around a solar-type star. The original paper, by Mayor & Queloz (1995), used data from the ELODIE spectrograph. Here we will use another dataset, obtained with the Hamilton echelle spectrograph, situated in the Lick Observatory, California.

The data is in the 51Peg.rv file. It has 256 observations separated in three columns: time (in Modified Julian Days), measured radial velocity (in m/s) and its uncertainty (also in m/s). This type of file is the typical input to kima.

# 51Peg.rv
50002.665695	-52.9	 4.1
50002.684340	-45.8	 4.8
50002.800220	-60.8	 4.6
50002.815961	-53.3	 5.0
...

Our goal is to fit a Keplerian model to these radial-velocity (RV) observations.
To do this with kima, we prepare a kima_setup.cpp file (there's one already in the folder, we'll just go through it).

This C++ file starts with some general includes

#include "DNest4.h"
#include "Data.h"
#include "RVmodel.h"
using namespace DNest4;

Afterwards, we include kima's default priors, some of which we will change later on

#include "default_priors.h"

and we define general settings for the model

const bool obs_after_HARPS_fibers = false;
const bool GP = false;
const bool hyperpriors = false;
const bool trend = false;

These settings tell kima that we don't need to fit an RV offset due to the HARPS fibers (see here), we will use a Keplerian model without considering correlated noise (instead of a Gaussian process model), and we will not consider hyperpriors for the orbital parameters nor a RV linear trend.

Next comes the RVmodel constructor, where we set a few more options for the model and the priors.

RVmodel::RVmodel()
:planets(5, 1, false, RVConditionalPrior())
,mu(Data::get_instance().N())
,C(Data::get_instance().N(), Data::get_instance().N())
{
    Cprior = new Uniform(-10, 10);
    Jprior = new ModifiedLogUniform(1.0, 2000.);
    Pprior = new LogUniform(0.2, 1E2);
    Kprior = new ModifiedLogUniform(1.0, 2E3);
    eprior = new Uniform(0., 1.);
    wprior = new Uniform(0.0, 2*M_PI);
    phiprior = new Uniform(0.0, 2*M_PI);
    save_setup();
}

The line

planets(5, 1, false, RVConditionalPrior())

first tells kima how many orbital parameters there are in a Keplerian model (5). These are the orbital period P, semi-amplitude K, eccentricity e, and two angles w and phi.
This number should not be changed.

Next we set the maximum number of planets in the model to 1, and tell kima not to fix it during the analysis (with the third argument fix=false). That is, the number of planets is a free parameter, having a uniform prior between 0 and 1, in this case.

We then define priors for all the parameters

// systemic velocity; m/s
Cprior = new Uniform(-10, 10); 
// additional white noise; m/s
Jprior = new ModifiedLogUniform(1.0, 2000.);
// orbital period; days
Pprior = new LogUniform(0.2, 1E2); // days
// semi-amplitude; m/s
Kprior = new ModifiedLogUniform(1.0, 2E3);
// eccentricity
eprior = new Uniform(0., 1.);
// orbital angles, w and phi; radians
wprior = new Uniform(0.0, 2*M_PI);
phiprior = new Uniform(0.0, 2*M_PI);

We're almost ready. All that's left is to read the data file and run the model, which is done inside a main function

int main(int argc, char** argv)
{
    /* set the RV data file */
    char* datafile = "51Peg.rv";

    // the second argument is the units of the RVs
    // the third (optional) argument, tells kima 
    // not to skip any line in the header of the file
    Data::get_instance().load(datafile, "ms", 0);
    
    // set the sampler and run it!
    Sampler<RVmodel> sampler = setup<RVmodel>(argc, argv);
    sampler.run(10); // print to terminal every 10 saves
    return 0;
}

And that's it, we can now compile and run the model

make
./run

kima will start (using 4 threads, by default) and print some information to the terminal. The full example, with the default options, will take around 5 min to finish.

In the terminal, kima will say that it read the RVs from the file, is creating levels at given likelihood values and saving particles. The default option is to stop after saving 5000 particles.

But we don't have to wait for it to finish!

After opening another terminal in the same folder (kima/examples/51Peg), we can use the kima-showresults script included in the pykima package to look at the results.

You can run

kima-showresults -h

to see all the available options.

Running

kima-showresults rv planets orbital gp extra

will display a number of plots with the output of the run.
Note: there's no problem if you run kima-showresults while kima is still running!