Skip to content

API documentation for pykima

kima is written in C++, to leverage the increased performance. But a more dynamic language like Python is better suited to analyse the results and create figures. The pykima package was created to provide a bridge between the two languages. It can be used in an IPython shell, a Jupyter notebook, or the standard Python interpreter. Simply type

import pykima as pk

The pykima.results module contains the KimaResults class for storing, analysing, and plotting the results from a kima run. Users will typically not instantiate the class directly, but instead use the high-level load function

>>> res = pk.load()
>>> res
KimaResults(lnZ=..., ESS=...)

This class has a number of useful attributes and methods, described below.

pykima.results API

This module defines the KimaResults class to hold results from a run.

KimaResults(options, save_plots=False, return_figs=True, verbose=False)

A class to hold, analyse, and display the results from kima

Attributes:

Name Type Description
model str

The type of kima model

priors dict

A dictionary with the priors used in the model

ESS int

Effective sample size

evidence float

The log-evidence (\(\ln Z\)) of the model

information float

The Kulback-Leibler divergence between prior and posterior

data data_holder

The data

posteriors posterior_holder

The marginal posterior samples

apply_cuts_period(samples, pmin=None, pmax=None, return_mask=False)

apply cuts in orbital period

burst_model(sample, t=None, v=None)

For models with multiple instruments, this function "bursts" the computed RV into n_instruments individual arrays. This is mostly useful for plotting the RV model together with the original data.

Parameters:

Name Type Description Default
sample ndarray

One posterior sample, with shape (npar,)

required
t ndarray

Times at which to evaluate the model

None
v ndarray

Pre-computed RVs. If None, calls self.eval_model

None

eval_model(sample, t=None, include_planets=True, include_known_object=True, include_trend=True, single_planet=None, except_planet=None)

Evaluate the deterministic part of the model at one posterior sample. If t is None, use the observed times. Instrument offsets are only added if t is None, but the systemic velocity is always added. To evaluate at all posterior samples, consider using np.apply_along_axis(self.eval_model, 1, self.posterior_sample)

Note: this function does not evaluate the GP component of the model.

Parameters:

Name Type Description Default
sample ndarray

One posterior sample, with shape (npar,)

required
t ndarray

Times at which to evaluate the model, or None to use observed times

None
include_planets bool

Whether to include the contribution from the planets

True
include_known_object bool

Whether to include the contribution from the known object planets

True
include_trend bool

Whether to include the contribution from the trend

True
single_planet int

Index of a single planet to include in the model, starting at 1. Use positive values (1, 2, ...) for the Np planets and negative values (-1, -2, ...) for the known object planets.

None
except_planet Union[int, List]

Index (or list of indices) of a single planet to exclude from the model, starting at 1. Use positive values (1, 2, ...) for the Np planets and negative values (-1, -2, ...) for the known object planets.

None

from_prior(n=1)

Generate n samples from the priors for all parameters.

full_model(sample, t=None, **kwargs)

Evaluate the full model at one posterior sample, including the GP. If t is None, use the observed times. Instrument offsets are only added if t is None, but the systemic velocity is always added. To evaluate at all posterior samples, consider using

np.apply_along_axis(self.full_model, 1, self.posterior_sample)

Parameters:

Name Type Description Default
sample ndarray

One posterior sample, with shape (npar,)

required
t ndarray

Times at which to evaluate the model, or None to use observed times

None
**kwargs

Keyword arguments passed directly to eval_model

{}

get_marginals()

Get the marginal posteriors from the posterior_sample matrix. They go into self.T, self.A, self.E, etc

get_medians()

return the median values of all the parameters

load(filename=None, diagnostic=False, **kwargs) classmethod

Load a KimaResults object from the current directory, a pickle file, or a zip file.

Parameters:

Name Type Description Default
filename str

If given, load the model from this file. Can be a zip or pickle file. Defaults to None.

None
diagnostic bool

Whether to plot the DNest4 diagnotics. Defaults to False.

False
**kwargs

Extra keyword arguments passed to showresults

{}

Returns:

Name Type Description
res KimaResults

An object holding the results

maximum_likelihood_sample(from_posterior=False, Np=None, printit=True, mask=None)

Get the maximum likelihood sample. By default, this is the highest likelihood sample found by DNest4. If from_posterior is True, this returns instead the highest likelihood sample from those that represent the posterior. The latter may change, due to random choices, between different calls to "showresults". If Np is given, select only samples with that number of planets.

median_sample(Np=None, printit=True)

Get the median posterior sample. If Np is given, select only from posterior samples with that number of planets.

parameter_priors() property

A list of priors which can be indexed using self.indices

planet_model(sample, t=None, include_known_object=True, single_planet=None, except_planet=None)

Evaluate the planet part of the model at one posterior sample. If t is None, use the observed times. To evaluate at all posterior samples, consider using

np.apply_along_axis(self.planet_model, 1, self.posterior_sample)
Note

this function does not evaluate the GP component of the model nor the systemic velocity and instrument offsets.

Parameters:

Name Type Description Default
sample ndarray

One posterior sample, with shape (npar,)

required
t ndarray

Times at which to evaluate the model, or None to use observed times

None
include_known_object bool

Whether to include the contribution from the known object planets

True
single_planet int

Index of a single planet to include in the model, starting at 1. Use positive values (1, 2, ...) for the Np planets and negative values (-1, -2, ...) for the known object planets.

None
except_planet Union[int, List]

Index (or list of indices) of a single planet to exclude from the model, starting at 1. Use positive values (1, 2, ...) for the Np planets and negative values (-1, -2, ...) for the known object planets.

None

save_pickle(filename, verbose=True)

Pickle this KimaResults object into a file.

Parameters:

Name Type Description Default
filename str

The name of the file where to save the model

required
verbose bool

Print a message. Defaults to True.

True

save_zip(filename, verbose=True)

Save this KimaResults object and the text files into a zip.

Parameters:

Name Type Description Default
filename str

The name of the file where to save the model

required
verbose bool

Print a message. Defaults to True.

True

stochastic_model(sample, t=None, return_std=False, derivative=False, **kwargs)

Evaluate the stochastic part of the model (GP) at one posterior sample. If t is None, use the observed times. Instrument offsets are only added if t is None, but the systemic velocity is always added. To evaluate at all posterior samples, consider using

np.apply_along_axis(self.stochastic_model, 1, self.posterior_sample)

Parameters:

Name Type Description Default
sample ndarray

One posterior sample, with shape (npar,)

required
t ndarray

Times at which to evaluate the model, or None to use observed times

None
return_std bool

Whether to return the standard deviation of the predictive. Default is False.

False
derivative bool

Return the first time derivative of the GP prediction instead

False

data_holder dataclass

A simple class to hold the datasets used in kima

Attributes:

Name Type Description
t ndarray

The observation times

y ndarray

The observed radial velocities

e ndarray

The radial velocity uncertainties

obs ndarray

Identifier for the instrument of each observation

N int

Total number of observations

posterior_holder dataclass

A simple class to hold the posterior samples

Attributes:

Name Type Description
P ndarray

The orbital period(s)

K ndarray

The semi-amplitude(s)

e ndarray

The orbital eccentricities(s)

ω ndarray

The argument(s) of pericenter

φ ndarray

The mean anomaly(ies) at the epoch


A number of functions for more complicated analysis of the results are implemented in the pykima.analysis module. Most of these functions take a KimaResults instance as argument.

pykima.analysis API

detection_limits(results, star_mass=1.0, Np=None, bins=200, plot=True, semi_amplitude=False, semi_major_axis=False, logX=True, return_mask=False, remove_nan=True, show_eccentricity=False, K_lines=(5, 3, 1), smooth=False, smooth_degree=3)

Calculate detection limits using samples with more than Np planets. By default, this function uses the value of Np which passes the posterior probability threshold.

Parameters:

Name Type Description Default
star_mass Union[float, Tuple]

Stellar mass and optionally its uncertainty [in solar masses].

1.0
Np int

Consider only posterior samples with more than Np planets.

None
bins int

Number of bins at which to calculate the detection limits. The period ranges from the minimum to the maximum orbital period in the posterior.

200
plot bool

Whether to plot the detection limits

True
semi_amplitude bool

Show the detection limits for semi-amplitude, instead of planet mass

False
semi_major_axis bool

Show semi-major axis in the x axis, instead of orbital period

False
logX bool

Should X-axis bins be logarithmic?

True
return_mask bool False
remove_nan bool

remove bins with no counts (no posterior samples)

True
smooth bool

Smooth the binned maximum with a polynomial

False
smooth_degree int

Degree of the polynomial used for smoothing

3

Returns:

Name Type Description

P,K,E,M (ndarray): Orbital periods, semi-amplitudes, eccentricities, and planet masses used in the calculation of the detection limits. These correspond to all the posterior samples with more than Np planets

s namedtuple

Detection limits result, with attributes max and bins. The max array is in units of Earth masses, bins is in days.

find_outliers(results, sample, threshold=10, verbose=False)

Estimate which observations are outliers, for a model with a Student t likelihood. This function first calculates the residuals given the parameters in sample. Then it computes the relative probability of each residual point given a Student-t (Td) and a Gaussian (Gd) likelihoods. If the probability Td is larger than Gd (by a factor of threshold), the point is flagged as an outlier. The function returns an "outlier mask".

get_planet_mass(P, K, e, star_mass=1.0, full_output=False)

Calculate the planet (minimum) mass, \(M_p \sin i\), given orbital period P, semi-amplitude K, eccentricity e, and stellar mass. If star_mass is a tuple with (estimate, uncertainty), this (Gaussian) uncertainty will be taken into account in the calculation.

Parameters:

Name Type Description Default
P Union[float, ndarray]

orbital period [days]

required
K Union[float, ndarray]

semi-amplitude [m/s]

required
e Union[float, ndarray]

orbital eccentricity

required
star_mass Union[float, Tuple]

stellar mass, or (mass, uncertainty) [Msun]

1.0

This function returns different results depending on the inputs.

If P, K, and e are floats and star_mass is a float

Returns:

Name Type Description
Msini float

planet mass, in \(M_{\rm Jup}\)

Msini float

planet mass, in \(M_{\rm Earth}\)

If P, K, and e are floats and star_mass is a tuple

Returns:

Name Type Description
Msini tuple

planet mass and uncertainty, in \(M_{\rm Jup}\)

Msini tuple

planet mass and uncertainty, in \(M_{\rm Earth}\)

If P, K, and e are arrays and full_output=True

Returns:

Name Type Description
m_Msini float

posterior mean for the planet mass, in \(M_{\rm Jup}\)

s_Msini float

posterior standard deviation for the planet mass, in \(M_{\rm Jup}\)

Msini array

posterior samples for the planet mass, in \(M_{\rm Jup}\)

If P, K, and e are arrays and full_output=False

Returns:

Name Type Description
m_Msini float

posterior mean for the planet mass, in \(M_{\rm Jup}\)

s_Msini float

posterior standard deviation for the planet mass, in \(M_{\rm Jup}\)

m_Msini float

posterior mean for the planet mass, in \(M_{\rm Earth}\)

s_Msini float

posterior standard deviation for the planet mass, in \(M_{\rm Earth}\)

get_planet_mass_and_semimajor_axis(P, K, e, star_mass=1.0, full_output=False, verbose=False)

Calculate the planet (minimum) mass Msini and the semi-major axis given orbital period P, semi-amplitude K, eccentricity e, and stellar mass. If star_mass is a tuple with (estimate, uncertainty), this (Gaussian) uncertainty will be taken into account in the calculation.

Units

P [days] K [m/s] e [] star_mass [Msun]

Returns:

Type Description

(M, A) where M is the output of get_planet_mass A is the output of get_planet_semimajor_axis

get_planet_semimajor_axis(P, K, star_mass=1.0, full_output=False)

Calculate the semi-major axis of the planet's orbit given orbital period P, semi-amplitude K, and stellar mass.

Parameters:

Name Type Description Default
P Union[float, ndarray]

orbital period [days]

required
K Union[float, ndarray]

semi-amplitude [m/s]

required
star_mass Union[float, Tuple]

stellar mass, or (mass, uncertainty) [Msun]

1.0

This function returns different results depending on the inputs.

If P and K are floats and star_mass is a float

Returns:

Name Type Description
a float

planet semi-major axis, in AU

If P and K are floats and star_mass is a tuple

Returns:

Name Type Description
a tuple

semi-major axis and uncertainty, in AU

If P and K are arrays and full_output=True

Returns:

Name Type Description
m_a float

posterior mean for the semi-major axis, in AU

s_a float

posterior standard deviation for the semi-major axis, in AU

a array

posterior samples for the semi-major axis, in AU

If P and K are arrays and full_output=False

Returns:

Name Type Description
m_a float

posterior mean for the semi-major axis, in AU

s_a float

posterior standard deviation for the semi-major axis, in AU

np_bayes_factor_threshold(results, threshold=150)

Return the value of Np supported by the data considering a posterior ratio (Bayes factor) threshold.

Parameters:

Name Type Description Default
results KimaResults

A results instance

required
threshold float

Posterior ratio threshold.

150

np_most_probable(results)

Return the value of Np with the highest posterior probability.

Parameters:

Name Type Description Default
results KimaResults

A results instance

required

np_posterior_threshold(results, threshold=0.9)

Return the value of Np supported by the data considering an absolute posterior probability threshold.

Parameters:

Name Type Description Default
results KimaResults

A results instance

required
threshold float

Posterior probability threshold

0.9


The pykima.utils module contains a few less commonly used utility functions:

pykima.utils API

Small (but sometimes important) utility functions

Fixed(value)

Not a real distribution, just a parameter fixed at a given value

logpdf(x)

Logarithm of the probability density function

pdf(x)

Probability density function: 1.0 if x==value, otherwise 0.0

rvs()

Random sample (always returns value)

ZeroDist()

A dummy probability distribution which always returns 0.0

apply_argsort(arr1, arr2, axis=-1)

Apply arr1.argsort() on arr2, along axis.

chdir(dir)

A simple context manager to switch directories temporarily

clipped_mean(arr, min, max)

Mean of arr between min and max

clipped_std(arr, min, max)

Standard deviation of arr between min and max

date_to_jd(year, month, day)

Convert a date (year, month, day) to Julian Day.

Algorithm from 'Practical Astronomy with your Calculator or Spreadsheet', 4th ed., Duffet-Smith and Zwart, 2011.

Parameters:

Name Type Description Default
year int

Year as integer. Years preceding 1 A.D. should be 0 or negative. The year before 1 A.D. is 0, 10 B.C. is year -9.

required
month int

Month as integer, Jan = 1, Feb. = 2, etc.

required
day float

Day, may contain fractional part.

required

Returns:

Name Type Description
jd float

Julian Day

Examples:

Convert 6 a.m., February 17, 1985 to Julian Day

>>> date_to_jd(1985, 2, 17.25)
2446113.75

datetime_to_jd(date)

Convert a datetime.datetime object to Julian day.

Parameters:

Name Type Description Default
date datetime

the date to be converted to Julian day

required

Returns:

Name Type Description
jd float

Julian day

Examples:

>>> d = datetime.datetime(1985, 2, 17, 6)
>>> d
datetime.datetime(1985, 2, 17, 6, 0)
>>> jdutil.datetime_to_jd(d)
2446113.75

days_to_hmsm(days)

Convert fractional days to hours, minutes, seconds, and microseconds. Precision beyond microseconds is rounded to the nearest microsecond.

Parameters:

Name Type Description Default
days float

A fractional number of days. Must be less than 1.

required

Returns:

Name Type Description
hour int

Hour number

min int

Minute number

sec int

Second number

micro int

Microsecond number

Raises:

Type Description
ValueError

If days is >= 1.

Examples:

>>> days_to_hmsm(0.1)
(2, 24, 0, 0)

find_prior_limits(prior)

Find lower and upper limits of a prior from the kima_model_setup.txt file.

find_prior_parameters(prior)

Find parameters of a prior from the kima_model_setup.txt file.

get_instrument_name(data_file)

Find instrument name

get_planet_teq(Tstar=5777, Rstar=1, a=1, A=0, f=1)

Calculate the planet's equlibrium temperature

Parameters:

Name Type Description Default
Tstar float

Stellar effective temperature [K]. Defaults to 5777.

5777
Rstar float

Stellar radius [Rsun]. Defaults to 1.0.

1
a float

Planet semi-major axis [AU]. Defaults to 1.0.

1
A float

Bond albedo. Defaults to 0.

0
f float

Redistribution factor. Defaults to 1.

1

Returns:

Name Type Description
Teq float

Planet equiolibrium temperature [K]

get_prior(prior)

Return a scipt.stats-like prior from a kima_model_setup.txt prior

get_star_name(data_file)

Find star name (usually works for approx. standard filenames)

get_transit_probability(Rstar=1.0, a=1.0)

Transit probability, simple. Eq 6.1 in Exoplanet Handbook Rstar: stellar radius [Rsun] a: semi-major axis [au]

hmsm_to_days(hour=0, min=0, sec=0, micro=0)

Convert hours, minutes, seconds, and microseconds to fractional days.

Parameters:

Name Type Description Default
hour int

Hour number. Defaults to 0.

0
min int

Minute number. Defaults to 0.

0
sec int

Second number. Defaults to 0.

0
micro int

Microsecond number. Defaults to 0.

0

Returns:

Name Type Description
days float

Fractional days

Examples:

>>> hmsm_to_days(hour=6)
0.25

jd_to_date(jd)

Convert Julian Day to date.

Algorithm from 'Practical Astronomy with your Calculator or Spreadsheet', 4th ed., Duffet-Smith and Zwart, 2011.

Parameters:

Name Type Description Default
jd float

Julian Day

required

Returns:

Name Type Description
year int

Year as integer. Years preceding 1 A.D. should be 0 or negative. The year before 1 A.D. is 0, 10 B.C. is year -9.

month int

Month as integer, Jan = 1, Feb. = 2, etc.

day float

Day, may contain fractional part.

Examples:

Convert Julian Day 2446113.75 to year, month, and day.

>>> jd_to_date(2446113.75)
(1985, 2, 17.25)

jd_to_mjd(jd)

Convert Julian Day to Modified Julian Day

Parameters:

Name Type Description Default
jd float

Julian Day

required

Returns mjd (float): Modified Julian Day

lighten_color(color, amount=0.5)

Lightens the given color by multiplying (1-luminosity) by the given amount. Input can be a matplotlib color string, hex string, or RGB tuple.

Examples:

lighten_color('g', 0.3) lighten_color('#F034A3', 0.6) lighten_color((.3, .55, .1), 0.5)

mjd_to_jd(mjd)

Convert Modified Julian Day to Julian Day.

Parameters:

Name Type Description Default
mjd float

Modified Julian Day

required

Returns jd (float): Julian Day

read_datafile(datafile, skip)

Read data from datafile for multiple instruments. Can be str, in which case the 4th column is assumed to contain an integer identifier of the instrument. Or list, in which case each element will be one different filename containing three columns each.

read_datafile_rvfwhm(datafile, skip)

Read data from datafile for multiple instruments and RV-FWHM data. Can be str, in which case the 4th column is assumed to contain an integer identifier of the instrument. Or list, in which case each element will be one different filename containing three columns each.

rms(array)

Root mean square of array

show_tips()

Show a few tips on how to use kima (but only sometimes)

wrms(array, weights)

Weighted root mean square of array, given weights


Last update: 2022-08-19