In [2]:
from IPython.display import Image
Image("./_static/oop_modelling_diagram_v1.png")
Out[2]:
This notebook is a simplified version of the more extended modelling tutorial notebook shown in the X-PSI documentation page, made for the INT 2023 workshop on Neutron Rich Matter on Heaven and Earth. There are also other tutorial notebooks that we recommend trying, once familiar with this basic notebook on the modelling techniques. You could gather an idea of which tutorials might be relevant for your use-case in this page.
To run this tutorial, you should install X-PSI with the default (blackbody) atmosphere extension module, which is compiled when installing X-PSI with default settings. The installation instructions can be found here.
Some extra data files not available on the GitHub repository are needed to run this tutorial. They can be found in the Zenodo and should be placed in the directory examples/examples_modeling_tutorial/model_data/ of your local xpsi directory. These are the required files:
nicer_v1.01_arf.txtnicer_v1.01_rmf_matrix.txtnicer_v1.01_rmf_energymap.txt%matplotlib inline
import os
import numpy as np
import math
import time
from matplotlib import pyplot as plt
from matplotlib import rcParams
from matplotlib.ticker import MultipleLocator, AutoLocator, AutoMinorLocator
from matplotlib import gridspec
from matplotlib import cm
import xpsi
from xpsi.global_imports import _c, _G, _dpr, gravradius, _csq, _km, _2pi
/=============================================\ | X-PSI: X-ray Pulse Simulation and Inference | |---------------------------------------------| | Version: 2.0.2 | |---------------------------------------------| | https://xpsi-group.github.io/xpsi | \=============================================/ Imported GetDist version: 1.4.3 Imported nestcheck version: 0.2.1
In order to implement an X-PSI likelihood function it is first necessary to define an underlying model parameter space. In general the parameter space is $\mathbb{R}^{d}$, where $d\in\mathbb{N}$ is the total number of free model parameters. Each free parameter is an instance of the Parameter class. We also support fixed (or frozen) variables and derived variables that are deterministic functions of other variables (notably, free parameters).
Our aim is to construct a callable likelihood object to feed as a callback function to a posterior sampler and a posterior integrator. We illustrate this object in the diagram below; all nodes in the diagram below represent objects (instances of some class). In some cases there are multiple instances of a particular class, in other cases there is only a single instance of a particular class.
from IPython.display import Image
Image("./_static/oop_modelling_diagram_v1.png")
The arrows in the diagram denote an object that is stored as an attribute of another encapsulating object. When a call is processed to evaluate the likelihood given a vector of model parameter values, encapsulating objects get attributes (including method calls) of the wrapped objects.
We do not need to subclass data container, but you can if you wish. X-PSI is designed this way because there is a clear common usage pattern that can be concretely implemented whilst preserving the freedom and the scope of applicability of the source code. The container instance will be available as an underscore instance method of Signal, and thus available in a derived class where we will later write code for likelihood evaluation.
Hereafter we will write our custom derived classes in the notebook itself, but in practice it is best if your derived classes are written in distinct modules within a project directory, so they can be imported by a script for use with an MPI command within a shell (because in general we want to exploit parallelism for expensive likelihood evaluations).
Let us load a synthetic data set that we generated in advance, and know the fictitious exposure time for:
settings = dict(counts = np.loadtxt('../../examples/examples_modeling_tutorial/model_data/example_synthetic_realisation.dat', dtype=np.double),
channels=np.arange(20,201),
phases=np.linspace(0.0, 1.0, 33),
first=0, last=180,
exposure_time=984307.6661)
data = xpsi.Data(**settings)
Setting channels for event data... Channels set.
Let's take a look at the data that we aim to model. First we define some settings and helper functions:
rcParams['text.usetex'] = False
rcParams['font.size'] = 14.0
def veneer(x, y, axes, lw=1.0, length=8, yticks=None):
""" Make the plots a little more aesthetically pleasing. """
if x is not None:
if x[1] is not None:
axes.xaxis.set_major_locator(MultipleLocator(x[1]))
if x[0] is not None:
axes.xaxis.set_minor_locator(MultipleLocator(x[0]))
else:
axes.xaxis.set_major_locator(AutoLocator())
axes.xaxis.set_minor_locator(AutoMinorLocator())
if y is not None:
if y[1] is not None:
axes.yaxis.set_major_locator(MultipleLocator(y[1]))
if y[0] is not None:
axes.yaxis.set_minor_locator(MultipleLocator(y[0]))
else:
axes.yaxis.set_major_locator(AutoLocator())
axes.yaxis.set_minor_locator(AutoMinorLocator())
axes.tick_params(which='major', colors='black', length=length, width=lw)
axes.tick_params(which='minor', colors='black', length=int(length/2), width=lw)
plt.setp(axes.spines.values(), linewidth=lw, color='black')
if yticks:
axes.set_yticks(yticks)
def plot_one_pulse(pulse, x, label=r'Counts',
cmap=cm.magma, vmin=None, vmax=None,
yticks=[20,50,100,200]):
""" Plot a pulse resolved over a single rotational cycle. """
fig = plt.figure(figsize = (7,7))
gs = gridspec.GridSpec(1, 2, width_ratios=[50,1])
ax = plt.subplot(gs[0])
ax_cb = plt.subplot(gs[1])
# Calculate the centre of phase bins (required by pcolormesh instead of phase edges)
_phase_bincenter = x[:-1]+0.5*(x[1]-x[0])
profile = ax.pcolormesh(_phase_bincenter,
data.channels,
pulse,
vmin = vmin,
vmax = vmax,
cmap = cmap,
linewidth = 0,
rasterized = True)
profile.set_edgecolor('face')
ax.set_xlim([0.0, 1.0])
ax.set_yscale('log')
ax.set_ylabel(r'Channel')
ax.set_xlabel(r'Phase')
cb = plt.colorbar(profile,
cax = ax_cb)
cb.set_label(label=label, labelpad=25)
cb.solids.set_edgecolor('face')
veneer((0.05, 0.2), (None, None), ax, yticks=yticks)
plt.subplots_adjust(wspace = 0.025)
if yticks is not None:
ax.set_yticklabels(yticks)
Now for the data:
plot_one_pulse(data.counts, data.phases)
We require a model instrument object to transform incident specific flux signals into a form which enters directly in the sampling distribution of the data.
class CustomInstrument(xpsi.Instrument):
""" A model of the NICER telescope response. """
def __call__(self, signal, *args):
""" Overwrite base just to show it is possible.
We loaded only a submatrix of the total instrument response
matrix into memory, so here we can simplify the method in the
base class.
"""
matrix = self.construct_matrix()
self._folded_signal = np.dot(matrix, signal)
return self._folded_signal
@classmethod
def from_response_files(cls, ARF, RMF, channel_edges, max_input, min_input=0):
""" Constructor which converts response files into :class:`numpy.ndarray`s.
:param str ARF: Path to ARF which is compatible with
:func:`numpy.loadtxt`.
:param str RMF: Path to RMF which is compatible with
:func:`numpy.loadtxt`.
:param str channel_edges: Path to edges which is compatible with
:func:`numpy.loadtxt`.
"""
if min_input != 0:
min_input = int(min_input)
max_input = int(max_input)
try:
ARF = np.loadtxt(ARF, dtype=np.double, skiprows=3)
RMF = np.loadtxt(RMF, dtype=np.double)
channel_edges = np.loadtxt(channel_edges, dtype=np.double, skiprows=3)[:,1:]
except:
print('A file could not be loaded.')
raise
matrix = np.ascontiguousarray(RMF[min_input:max_input,20:201].T, dtype=np.double)
edges = np.zeros(ARF[min_input:max_input,3].shape[0]+1, dtype=np.double)
edges[0] = ARF[min_input,1]; edges[1:] = ARF[min_input:max_input,2]
for i in range(matrix.shape[0]):
matrix[i,:] *= ARF[min_input:max_input,3]
channels = np.arange(20, 201)
return cls(matrix, edges, channels, channel_edges[20:202,-2])
Let's construct an instance.
NICER = CustomInstrument.from_response_files(ARF = '../../examples/examples_modeling_tutorial/model_data/nicer_v1.01_arf.txt',
RMF = '../../examples/examples_modeling_tutorial/model_data/nicer_v1.01_rmf_matrix.txt',
channel_edges = '../../examples/examples_modeling_tutorial/model_data/nicer_v1.01_rmf_energymap.txt',
max_input = 500,
min_input = 0)
Setting channels for loaded instrument response (sub)matrix... Channels set. No parameters supplied... empty subspace created.
The NICER v1.01 response matrix:
fig = plt.figure(figsize = (14,7))
ax = fig.add_subplot(111)
veneer((25, 100), (10, 50), ax)
_ = ax.imshow(NICER.matrix,
cmap = cm.viridis,
rasterized = True)
ax.set_ylabel('Channel $-\;20$')
_ = ax.set_xlabel('Energy interval')
Summed over channel subset $[20,200]$:
fig = plt.figure(figsize = (7,7))
ax = fig.add_subplot(111)
veneer((0.1, 0.5), (50,250), ax)
ax.plot((NICER.energy_edges[:-1] + NICER.energy_edges[1:])/2.0, np.sum(NICER.matrix, axis=0), 'k-')
ax.set_ylabel('Effective area [cm$^{-2}$]')
_ = ax.set_xlabel('Energy [keV]')
We can now combine the dataset and model instrument into a Signal object. The source code for this class has methods and attributes which simplify communication between the aforementioned model objects and another object representing our model star (created below). The surface radiation field of the model star is integrated over based on energies relayed from a Signal object based on the properties of the instrument and the dataset (which are tightly coupled).
We are forced to inherit from Signal and write a method that evaluates the logarithm of the likelihood conditional on a parametrised sampling distribution for the data. There is much freedom in constructing this sampling distribution, so the design strategy for X-PSI was to leave this part of the modelling process entirely to a user, guided by a number of examples.
from xpsi.likelihoods.default_background_marginalisation import eval_marginal_likelihood
from xpsi.likelihoods.default_background_marginalisation import precomputation
class CustomSignal(xpsi.Signal):
"""
A custom calculation of the logarithm of the likelihood.
We extend the :class:`~xpsi.Signal.Signal` class to make it callable.
We overwrite the body of the __call__ method. The docstring for the
abstract method is copied.
"""
def __init__(self, workspace_intervals = 1000, epsabs = 0, epsrel = 1.0e-8,
epsilon = 1.0e-3, sigmas = 10.0, support = None, **kwargs):
""" Perform precomputation.
:params ndarray[m,2] support:
Prior support bounds for background count rate variables in the
:math:`m` instrument channels, where the lower bounds must be zero
or positive, and the upper bounds must be positive and greater than
the lower bound. Alternatively, setting the an upper bounds as
negative means the prior support is unbounded and the flat prior
density functions per channel are improper. If ``None``, the lower-
bound of the support for each channel is zero but the prior is
unbounded.
"""
super(CustomSignal, self).__init__(**kwargs)
try:
self._precomp = precomputation(self._data.counts.astype(np.int32))
except AttributeError:
print('Warning: No data... can synthesise data but cannot evaluate a '
'likelihood function.')
else:
self._workspace_intervals = workspace_intervals
self._epsabs = epsabs
self._epsrel = epsrel
self._epsilon = epsilon
self._sigmas = sigmas
if support is not None:
self._support = support
else:
self._support = -1.0 * np.ones((self._data.counts.shape[0],2))
self._support[:,0] = 0.0
def __call__(self, *args, **kwargs):
self.loglikelihood, self.expected_counts, self.background_signal, self.background_given_support = \
eval_marginal_likelihood(self._data.exposure_time,
self._data.phases,
self._data.counts,
self._signals,
self._phases,
self._shifts,
self._precomp,
self._support,
self._workspace_intervals,
self._epsabs,
self._epsrel,
self._epsilon,
self._sigmas,
kwargs.get('llzero'),
allow_negative=(False, False))
In the first part of this notebook we define a marginal likelihood function. That is, instead of invoking the true background model that in this case is known to us, we invoke a default treatment whereby we marginalise over a set of channel-by-channel background count-rate parameters instead.
We wrote our __call__ method as a wrapper for a extension module to improve speed. In general, if you wish to change the model for likelihood evaluation given expected signals, you can archive the Cython extensions in, e.g., the xpsi/likelihoods directory, and compile these when X-PSI is compiled and installed (by editing the setup.py script). Alternatively, you can compile your extension elsewhere and call those compiled binaries from your custom class derived from Signal.
Let's construct and instantiate a Signal object. The bounds of the background parameter have already been specified above.
xpsi.Signal?
signal = CustomSignal(data = data,
instrument = NICER,
background = None,
interstellar = None,
workspace_intervals = 1000,
cache = True,
epsrel = 1.0e-8,
epsilon = 1.0e-3,
sigmas = 10.0,
support = None)
Creating parameter:
> Named "phase_shift" with fixed value 0.000e+00.
> The phase shift for the signal, a periodic parameter [cycles].
We now need to build our star. The basic units for building a star are:
For this demonstration we will assume that the surface radiation field elsewhere (other than the hot regions) can be ignored in the soft X-ray regime our model instrument is sensitive to. Let's start with the Spacetime class.
xpsi.Spacetime?
We can instansiate the spacetime by defining all parameters with a given value. In this case, all parameters will be fixed because the bounds are not specified (empty dictionary). Note that all parameters must be defined at least once in bounds or values.
values = dict(distance =0.150, # (Earth) distance
mass = 1.5, # mass
radius = 12.0, # equatorial radius
cos_inclination = 0.5, # (Earth) inclination to rotation axis
frequency = 300. ) # spin frequency
spacetime = xpsi.Spacetime(bounds={}, values=values)
Creating parameter:
> Named "frequency" with fixed value 3.000e+02.
> Spin frequency [Hz].
Creating parameter:
> Named "mass" with fixed value 1.500e+00.
> Gravitational mass [solar masses].
Creating parameter:
> Named "radius" with fixed value 1.200e+01.
> Coordinate equatorial radius [km].
Creating parameter:
> Named "distance" with fixed value 1.500e-01.
> Earth distance [kpc].
Creating parameter:
> Named "cos_inclination" with fixed value 5.000e-01.
> Cosine of Earth inclination to rotation axis.
Alternatively we can specify bounds manually for the free parameters. Using (None, None) will set the default bounds. The fixed parameters are defined in the values dictionary. If both the bounds and a value are defined for a parameters, the value will not be fixed but instead will be set as an initial value.
bounds = dict(distance = (None, None), # Default bounds for (Earth) distance
mass = (None, None), # Default bounds for mass
radius = (None, None), # Default bounds for equatorial radius
cos_inclination = (None, None)) # Default bounds for (Earth) inclination to rotation axis
values = dict(frequency=300.0, # Fixed spin frequency
mass=1.4) # mass with initial value
spacetime = xpsi.Spacetime(bounds=bounds, values=values)
Creating parameter:
> Named "frequency" with fixed value 3.000e+02.
> Spin frequency [Hz].
Creating parameter:
> Named "mass" with bounds [1.000e-03, 3.000e+00] and initial value 1.400e+00.
> Gravitational mass [solar masses].
Creating parameter:
> Named "radius" with bounds [1.000e+00, 2.000e+01].
> Coordinate equatorial radius [km].
Creating parameter:
> Named "distance" with bounds [1.000e-02, 3.000e+01].
> Earth distance [kpc].
Creating parameter:
> Named "cos_inclination" with bounds [-1.000e+00, 1.000e+00].
> Cosine of Earth inclination to rotation axis.
To display the free parameters:
for p in spacetime:
print(p)
Gravitational mass [solar masses]. Coordinate equatorial radius [km]. Earth distance [kpc]. Cosine of Earth inclination to rotation axis.
For the following example in this tutorial, we will specify the bounds of all parameters.
bounds = dict(distance = (0.1, 1.0), # (Earth) distance
mass = (1.0, 3.0), # mass with default bounds
radius = (3.0 * gravradius(1.0), 16.0), # equatorial radius
cos_inclination = (0.0, 1.0)) # (Earth) inclination to rotation axis
values = dict(frequency=300.0 ) # Fixed spin frequency
spacetime = xpsi.Spacetime(bounds=bounds, values=values)
Creating parameter:
> Named "frequency" with fixed value 3.000e+02.
> Spin frequency [Hz].
Creating parameter:
> Named "mass" with bounds [1.000e+00, 3.000e+00].
> Gravitational mass [solar masses].
Creating parameter:
> Named "radius" with bounds [4.430e+00, 1.600e+01].
> Coordinate equatorial radius [km].
Creating parameter:
> Named "distance" with bounds [1.000e-01, 1.000e+00].
> Earth distance [kpc].
Creating parameter:
> Named "cos_inclination" with bounds [0.000e+00, 1.000e+00].
> Cosine of Earth inclination to rotation axis.
It is not necessary for us to write a custom derived class for the photosphere object, so we will simply instantiate a Photosphere object. However, we first need an instance of HotRegion to instantiate the photosphere, and we need to implement a low-level parametrised model for the specific intensity emergent from the photosphere in a local comoving frame.
The neutron star atmosphere is assumed to be geometrically thin. In the applications thus far, the surface local-comoving-frame radiation field as being described by a single free parameter: the effective temperature. The radiation field is also dependent on the local effective gravitational acceleration, however this is a derived parameter in the model. The parametrised radiation field as a function of energy and angle subtended to the normal to the (plane-parallel) atmosphere in a local comoving frame is provided as numerical model data for multi-dimensional interpolation.
In X-PSI, integration over the surface radiation field is performed via calls to low-level C routines. To reduce likelihood evaluation times, the atmosphere interpolator is written in C, and calls to that interpolator are from C routine. In other words, in X-PSI, we do not use Python callback functions for evaluation of specific intensities, but C functions which are compiled when the X-PSI package is built. Unfortunately this means that to change the parametrised surface radiation field you need to get your hands a little dirty; on the bright side, the body of these functions can be implemented almost completely in the Cython language, so syntactically there is some similarity to Python because the language syntax is somewhat of a hybrid/superset. Beware, however, that the body of these functions must not contain calls to the Python API, and only to external C libraries if required: the code must evaluate to pure C, and not require the Python/C API.
The following is the contents of the hot.pyx which by default includes an isotropic blackbody model:
#cython: cdivision=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
from libc.math cimport exp, pow
from xpsi.global_imports import _keV, _k_B
cdef int SUCCESS = 0
cdef double erg = 1.0e-7
cdef double Planck_dist_const = 5.040366110812353e22
cdef double k_B = _k_B
cdef double keV = _keV
cdef double k_B_over_keV = k_B / keV
#----------------------------------------------------------------------->>>
# >>> User modifiable functions.
# >>> Note that the user is entirely free to wrap thread-safe and
# ... non-parallel external C routines from an external library.
# >>> Thus the bodies of the following need not be written explicitly in
# ... the Cython language.
#----------------------------------------------------------------------->>>
cdef void* init_hot(size_t numThreads, const _preloaded *const preloaded) nogil:
# This function must match the free management routine free_hot()
# in terms of freeing dynamically allocated memory. This is entirely
# the user's responsibility to manage.
# Return NULL if dynamic memory is not required for the model.
return NULL
cdef int free_hot(size_t numThreads, void *const data) nogil:
# This function must match the initialisation routine init_hot()
# in terms of freeing dynamically allocated memory. This is entirely
# the user's responsibility to manage.
# The void pointer must be appropriately cast before memory is freed --
# only the user can know this at compile time.
# Just use free(<void*> data) iff no memory was dynamically
# allocated in the function:
# init_local_hot()
# because data is expected to be NULL in this case
#printf("\nNo data to be freed.")
return SUCCESS
cdef double eval_hot(size_t THREAD,
double E,
double mu,
const double *const VEC,
void *const data) nogil:
# Arguments:
# E = photon energy in keV
# mu = cosine of ray zenith angle (i.e., angle to surface normal)
# VEC = variables such as temperature, effective gravity, ...
# data = numerical model data required for intensity evaluation
cdef double temp = k_B_over_keV * pow(10.0, VEC[0])
return E * E * E / ( exp(E / temp) - 1.0 )
cdef double eval_hot_norm() nogil:
# Source radiation field normalisation which is independent of the
# parameters of the parametrised model -- i.e. cell properties, energy,
# and angle.
# Writing the normalisation here reduces the number of operations required
# during integration.
# The units of the specific intensity need to be J/cm^2/s/keV/steradian.
return erg * Planck_dist_const
In most use-cases we need to modify these functions to enable handling of the numerical atmosphere data. An extension for such a case may be found as an example, which contains the extension used by Riley et al. (2019) to implement a numerical atmosphere generated by the NSX atmosphere code (Ho & Lai (2001); Ho & Heinke (2009))
We now instantiate hot region objects. We can find the required and optional parameter names in the class docstring:
xpsi.HotRegion?
The names can also be found as class attributes as follows:
xpsi.HotRegion.required_names
['super_colatitude', 'super_radius', 'phase_shift', 'super_temperature (if no custom specification)']
The condition if no custom specification means that this name is required if we do not supply custom parameters for the radiation field in the superseding member of the hot region. If we supply custom parameters, we also need to subclass xpsi.HotRegion and overwrite the __compute_cellParamVecs method to handle our parameters.
xpsi.HotRegion.optional_names
['omit_colatitude', 'omit_radius', 'omit_azimuth', 'cede_colatitude', 'cede_radius', 'cede_azimuth', 'cede_temperature']
Let's now instantiate our first spot and call it primary:
bounds = dict(super_colatitude = (None, None),
super_radius = (None, None),
phase_shift = (-0.25, 0.75),
super_temperature = (None, None))
# a simple circular, simply-connected spot
primary = xpsi.HotRegion(bounds=bounds,
values={}, # no initial values and no derived/fixed
symmetry=True,
omit=False,
cede=False,
concentric=False,
sqrt_num_cells=32,
min_sqrt_num_cells=10,
max_sqrt_num_cells=64,
num_leaves=100,
num_rays=200,
prefix='p') # unique prefix needed because >1 instance
Creating parameter:
> Named "super_colatitude" with bounds [0.000e+00, 3.142e+00].
> The colatitude of the centre of the superseding region [radians].
Creating parameter:
> Named "super_radius" with bounds [0.000e+00, 1.571e+00].
> The angular radius of the (circular) superseding region [radians].
Creating parameter:
> Named "phase_shift" with bounds [-2.500e-01, 7.500e-01].
> The phase of the hot region, a periodic parameter [cycles].
Creating parameter:
> Named "super_temperature" with bounds [3.000e+00, 7.600e+00].
> log10(superseding region effective temperature [K]).
The arguments omit, cede, and concentric can be used to more complex geometries such as rings, crescents, dual temperature regions, etc. All other arguments determine the numerical resolution, and have defaults which have been (somewhat arbitrarily) chosen to be result in a likelihood evaluation time of $\mathcal{O}(1)$ s.
Now let's instantiate the second hot region and call it secondary:
class derive(xpsi.Derive):
def __init__(self):
"""
We can pass a reference to the primary here instead
and store it as an attribute if there is risk of
the global variable changing.
This callable can for this simple case also be
achieved merely with a function instead of a magic
method associated with a class.
"""
pass
def __call__(self, boundto, caller = None):
# one way to get the required reference
global primary # unnecessary, but for clarity
return primary['super_temperature'] - 0.2
#bounds['super_temperature'] = None # declare fixed/derived variable
secondary = xpsi.HotRegion(bounds=bounds, # can otherwise use same bounds
values={},#{'super_temperature': derive()}, # create a callable value
symmetry=True,
omit=False,
cede=False,
concentric=False,
sqrt_num_cells=32,
min_sqrt_num_cells=10,
max_sqrt_num_cells=100,
num_leaves=100,
num_rays=200,
is_antiphased=True,
prefix='s') # unique prefix needed because >1 instance
Creating parameter:
> Named "super_colatitude" with bounds [0.000e+00, 3.142e+00].
> The colatitude of the centre of the superseding region [radians].
Creating parameter:
> Named "super_radius" with bounds [0.000e+00, 1.571e+00].
> The angular radius of the (circular) superseding region [radians].
Creating parameter:
> Named "phase_shift" with bounds [-2.500e-01, 7.500e-01].
> The phase of the hot region, a periodic parameter [cycles].
Creating parameter:
> Named "super_temperature" with bounds [3.000e+00, 7.600e+00].
> log10(superseding region effective temperature [K]).
We now need to encapsulate the hot region instances in a container with properties expected by the Photosphere class.
from xpsi import HotRegions
hot = HotRegions((primary, secondary))
Let's check out the hot regions:
hot.objects[0] # 'p'
Free parameters --------------- p__phase_shift: The phase of the hot region, a periodic parameter [cycles]. p__super_colatitude: The colatitude of the centre of the superseding region [radians]. p__super_radius: The angular radius of the (circular) superseding region [radians]. p__super_temperature: log10(superseding region effective temperature [K]). Derived/fixed parameters ------------------------ p__cede_colatitude: The colatitude of the centre of the ceding region [radians]. p__cede_radius: The angular radius of the (circular) ceding region [radians]. p__cede_azimuth: The azimuth of the centre of the ceding region relative to the centre of the superseding region [radians]. p__omit_colatitude: The colatitude of the centre of the omission region [radians]. p__omit_radius: The angular radius of the (circular) omission region [radians]. p__omit_azimuth: The azimuth of the centre of the omission region relative to the centre of the superseding region [radians].
hot.objects[1] # 's'
Free parameters --------------- s__phase_shift: The phase of the hot region, a periodic parameter [cycles]. s__super_colatitude: The colatitude of the centre of the superseding region [radians]. s__super_radius: The angular radius of the (circular) superseding region [radians]. s__super_temperature: log10(superseding region effective temperature [K]). Derived/fixed parameters ------------------------ s__cede_colatitude: The colatitude of the centre of the ceding region [radians]. s__cede_radius: The angular radius of the (circular) ceding region [radians]. s__cede_azimuth: The azimuth of the centre of the ceding region relative to the centre of the superseding region [radians]. s__omit_colatitude: The colatitude of the centre of the omission region [radians]. s__omit_radius: The angular radius of the (circular) omission region [radians]. s__omit_azimuth: The azimuth of the centre of the omission region relative to the centre of the superseding region [radians].
A list of names, with the prefix, can also be accessed as follows:
h = hot.objects[0]
h.names
['p__phase_shift', 'p__super_colatitude', 'p__super_radius', 'p__super_temperature', 'p__cede_colatitude', 'p__cede_radius', 'p__cede_azimuth', 'p__omit_colatitude', 'p__omit_radius', 'p__omit_azimuth']
h.prefix
'p'
h.get_param('phase_shift')
The phase of the hot region, a periodic parameter [cycles]
Let's set a value for the temperature of the primary hot region:
hot['p__super_temperature'] = 6.0 # equivalent to ``primary['super_temperature'] = 6.0``
We can now instantitate the photosphere:
photosphere = xpsi.Photosphere(hot = hot, elsewhere = None,
values=dict(mode_frequency = spacetime['frequency']))
Creating parameter:
> Named "mode_frequency" with fixed value 3.000e+02.
> Coordinate frequency of the mode of radiative asymmetry in the
photosphere that is assumed to generate the pulsed signal [Hz].
We do not model the surface radiation field elsewhere, and we thus leave the elsewhere keyword as None (the default). Elsewhere means on the surface, exterior to radiating hot regions that are typically expected to span a smaller angular extent; in the current version, the radiation from elsewhere, if explicitly computed is assumed to be time-invariant supposing the hot regions were not present. To account for radiation from elsewhere, a time-invariant signal is first computed, meaning an axisymmetric surface local-comoving-frame radiation field is assumed. The time-dependent signals from the hot regions are then computed, and modified by subtracting the specific intensity that would otherwise be generated by the surface local-comoving-frame radiation field from elsewhere (i.e., in place of the hot regions).
We can now combine the ambient spacetime, spacetime, and the embedded photosphere, photosphere, into a model star represented by an instance of Star. We do not need to extend this class, so we can simply construct and instantiate a star as follows:
star = xpsi.Star(spacetime = spacetime, photospheres = photosphere)
Let's check out the star object, which merged parameter subspaces associated with objects lower in the hierarchy:
star
Free parameters --------------- mass: Gravitational mass [solar masses]. radius: Coordinate equatorial radius [km]. distance: Earth distance [kpc]. cos_inclination: Cosine of Earth inclination to rotation axis. p__phase_shift: The phase of the hot region, a periodic parameter [cycles]. p__super_colatitude: The colatitude of the centre of the superseding region [radians]. p__super_radius: The angular radius of the (circular) superseding region [radians]. p__super_temperature: log10(superseding region effective temperature [K]). s__phase_shift: The phase of the hot region, a periodic parameter [cycles]. s__super_colatitude: The colatitude of the centre of the superseding region [radians]. s__super_radius: The angular radius of the (circular) superseding region [radians]. s__super_temperature: log10(superseding region effective temperature [K]).
Given the objects constructed above and the relevant pre-compiled low-level code, we can now construct and instantiate a callable likelihood object. We do not need extend (via inheritance) the Likelihood class found the source code: this class simply combines all of the model objects defined above, performs some automatic operations given the properties of the those objects, and facilitates communication of those objects when it recieves a call to evaluate the likelihood.
likelihood = xpsi.Likelihood(star = star, signals = signal,
num_energies=128,
threads=1,
externally_updated=False)
Let's retrieve the total number of free model parameters merged into the full parameter space:
likelihood
Free parameters --------------- mass: Gravitational mass [solar masses]. radius: Coordinate equatorial radius [km]. distance: Earth distance [kpc]. cos_inclination: Cosine of Earth inclination to rotation axis. p__phase_shift: The phase of the hot region, a periodic parameter [cycles]. p__super_colatitude: The colatitude of the centre of the superseding region [radians]. p__super_radius: The angular radius of the (circular) superseding region [radians]. p__super_temperature: log10(superseding region effective temperature [K]). s__phase_shift: The phase of the hot region, a periodic parameter [cycles]. s__super_colatitude: The colatitude of the centre of the superseding region [radians]. s__super_radius: The angular radius of the (circular) superseding region [radians]. s__super_temperature: log10(superseding region effective temperature [K]).
len(likelihood)
12
Let's call the likelihood object with the true model parameter values that we injected to generate the synthetic data rendered above, omitting background parameters.
p = [1.4,
12.5,
0.2,
math.cos(1.25),
0.0,
1.0,
0.075,
6.2,
0.025,
math.pi - 1.0,
0.2,
6.0]
likelihood.clear_cache()
t = time.time()
# source code changes since model was applied, so let's be a
# bit lenient when checking the likelihood function
likelihood.check(None, [-26713.6136777], 1.0e-6,
physical_points=[p])
print('time = %.3f s' % (time.time() - t))
Checking likelihood and prior evaluation before commencing sampling... Not using ``allclose`` function from NumPy. Using fallback implementation instead. Checking closeness of likelihood arrays: -2.6713601352e+04 | -2.6713613678e+04 ..... Closeness evaluated. Log-likelihood value checks passed on root process. Checks passed. time = 0.542 s
External sampling software will interact with a likelihood object in this way. That is, it will pass some ordered container of parameter values: a vector. The external sampler will also interact with likelihood object to acces the prior object to inverse sample from the joint prior distribution. Our prior object can also interact with our likelihood object outside of a sampling process.
Calling the likelihood object also modified the signals property of the photosphere object. Let's plot the signals by summing the count-rate over output instrument channels. We first define a helper function.
def plot_pulse(pulse1, pulse2, pulse3, pulse4, phase1, phase2):
""" Plot hot region signals before and after telescope operation. """
fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111)
ax.set_ylabel('Signal [arbitrary normalisation]')
ax.set_xlabel('Phase [cycles]')
temp = np.sum(pulse1, axis=0)
ax.plot(phase1, temp/np.max(temp), '-', color='k', lw=0.5)
temp = np.sum(pulse2, axis=0)
ax.plot(phase2, temp/np.max(temp), '-', color='r', lw=0.5)
temp = np.sum(pulse3, axis=0)
ax.plot(phase1, temp/np.max(temp), 'o-', color='k', lw=0.5, markersize=2)
temp = np.sum(photosphere.signal[1][0], axis=0)
ax.plot(phase2, temp/np.max(temp), 'o-', color='r', lw=0.5, markersize=2)
veneer((0.05,0.2), (0.05,0.2), ax)
likelihood(p, reinitialise=False)
_ = plot_pulse(signal.signals[0], signal.signals[1],
photosphere.signal[0][0], photosphere.signal[1][0],
signal.phases[0], signal.phases[1])
The pulse-profiles with markers are the signals incident on the telescope, before operating on them with the response model. The markers, linearly spaced in phase, denote the phase resolution.
The likelihood object calls the photosphere.integrate method, passing the energies stored as the property signal.energies. We can do this manually if we wish to integrate signals but not calculate likelihoods. Here we sum over incident specific photon flux pulses as an approximation to integrating over energy. Note that we do not change the signal.signals traced by the solid curves without markers.
Below we print crude representations of the cell meshes spanning each hot region.
fig = plt.figure(figsize = (14,7))
gs = gridspec.GridSpec(1, 3, width_ratios=[50,50,1], wspace=0.2)
ax = plt.subplot(gs[0])
veneer((1,5), (1, 5), ax)
# primary (lower colatitude) hot region
h = hot.objects[0]
z = h._HotRegion__cellArea[0]/np.max(h._HotRegion__cellArea[0])
patches = plt.pcolormesh(z,
vmin = np.min(z),
vmax = np.max(z),
cmap = cm.magma,
linewidth = 1.0,
rasterized = True,
edgecolor='black')
ax = plt.subplot(gs[1])
veneer((1,5), (1, 5), ax)
# secondary (higher colatitude) hot region
h = hot.objects[1]
z = h._HotRegion__cellArea[0]/np.max(h._HotRegion__cellArea[0])
_ = plt.pcolormesh(z,
vmin = np.min(z),
vmax = np.max(z),
cmap = cm.magma,
linewidth = 1.0,
rasterized = True,
edgecolor='black')
ax_cb = plt.subplot(gs[2])
cb = plt.colorbar(patches,
cax = ax_cb,
ticks = MultipleLocator(0.2))
cb.set_label(label = r'cell area (normalised by maximum)', labelpad=25)
cb.solids.set_edgecolor('face')
veneer((None, None), (0.05, None), ax_cb)
cb.outline.set_linewidth(1.0)
Note that the lowest colatitude row is at zero on the y-axis. The elements of a mesh cell-area array which are finite are not all identical: at the boundary of a hot region the proper area elements are smaller because of partial coverage by radiating material. The sum of all finite proper areas effectively equals the total proper area within a hot-region boundary.
We can also use projection tool to visualize the geometry of the neutron star hot regions (see the separate documenation page for more information about the tool):
from xpsi.utilities import ProjectionTool
ProjectionTool.plot_projection_general((likelihood),"ST-U","I","NP")
YOU ARE USING A 2 HOT SPOT MODEL
<AxesSubplot:>
#%tb
likelihood.params
[Gravitational mass [solar masses] = 1.400e+00, Coordinate equatorial radius [km] = 1.250e+01, Earth distance [kpc] = 2.000e-01, Cosine of Earth inclination to rotation axis = 3.153e-01, The phase of the hot region, a periodic parameter [cycles] = 0.000e+00, The colatitude of the centre of the superseding region [radians] = 1.000e+00, The angular radius of the (circular) superseding region [radians] = 7.500e-02, log10(superseding region effective temperature [K]) = 6.200e+00, The phase of the hot region, a periodic parameter [cycles] = 2.500e-02, The colatitude of the centre of the superseding region [radians] = 2.142e+00, The angular radius of the (circular) superseding region [radians] = 2.000e-01, log10(superseding region effective temperature [K]) = 6.000e+00]
Let's plot a pulse in two dimensions. Also note that we can interpolate the signal in phase as follows.
from xpsi.tools import phase_interpolator
def plot_2D_pulse(z, x, shift, y, ylabel,
num_rotations=5.0, res=5000, figsize=(12,6),
cm=cm.viridis,
yticks=None):
""" Helper function to plot a phase-energy pulse.
:param array-like z:
A pair of *ndarray[m,n]* objects representing the signal at
*n* phases and *m* values of an energy variable.
:param ndarray[n] x: Phases the signal is resolved at.
:param tuple shift: Hot region phase parameters.
:param ndarray[m] x: Energy values the signal is resolved at.
"""
fig = plt.figure(figsize = figsize)
gs = gridspec.GridSpec(1, 2, width_ratios=[50,1], wspace=0.025)
ax = plt.subplot(gs[0])
ax_cb = plt.subplot(gs[1])
new_phases = np.linspace(0.0, num_rotations, res)
interpolated = phase_interpolator(new_phases,
x,
z[0], shift[0])
interpolated += phase_interpolator(new_phases,
x,
z[1], shift[1])
profile = ax.pcolormesh(new_phases,
y,
interpolated/np.max(interpolated),
cmap = cm,
linewidth = 0,
rasterized = True)
profile.set_edgecolor('face')
ax.set_xlim([0.0, num_rotations])
ax.set_yscale('log')
ax.set_ylabel(ylabel)
ax.set_xlabel(r'Phase')
veneer((0.1, 0.5), (None,None), ax, yticks=yticks)
if yticks is not None:
ax.set_yticklabels(yticks)
cb = plt.colorbar(profile,
cax = ax_cb,
ticks = MultipleLocator(0.2))
cb.set_label(label=r'Signal (normalised by maximum)', labelpad=25)
cb.solids.set_edgecolor('face')
veneer((None, None), (0.05, None), ax_cb)
cb.outline.set_linewidth(1.0)
Let's compute the incident specific flux signal normalised to the maximum. Note that if want to obtain the specific flux in units of photons/cm$^{2}$/s/keV instead, the output of photosphere.signal needs to be divided by distance squared, where distance is measured in meters.
plot_2D_pulse((photosphere.signal[0][0], photosphere.signal[1][0]),
x=signal.phases[0],
shift=signal.shifts,
y=signal.energies,
ylabel=r'Energy (keV)',
yticks=[0.1, 0.5,1.0,2.0])
The count rate signal in each channel:
plot_2D_pulse((signal.signals[0], signal.signals[1]),
x=signal.phases[0],
shift=signal.shifts,
y=NICER.channels,
ylabel=r'Channels',
cm=cm.magma,
yticks=[20,50,100,200])
Now we increase the phase resolution, and plot a single rotational pulse:
for obj in hot.objects:
obj.set_phases(num_leaves = 1024)
# the current relationship between objects requires that we reinitialise
# if we wish to automatically communicate the updated settings between objects
p[3] = math.cos(1.0) # set inclination
_ = likelihood(p, reinitialise = True)
plot_2D_pulse((signal.signals[0], signal.signals[1]),
x=signal.phases[0],
shift=signal.shifts,
y=NICER.channels,
ylabel=r'Channels',
num_rotations=1.0,
cm=cm.magma,
yticks=[20,50,100,200],
figsize=(7,7))
Finally, let's play with some geometric parameters:
fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111)
ax.set_ylabel('photons/cm$^2$/s/keV (normalised by maxima)')
ax.set_xlabel('Phase [cycles]')
for obj in hot.objects:
obj.set_phases(num_leaves = 256)
# let's play with the angular radius of the primary hot region
angular_radii = np.linspace(0.01, 1.0, 10)
likelihood.externally_updated = True
likelihood['cos_inclination'] = math.cos(1.25)
for angular_radius in angular_radii:
likelihood['p__super_radius'] = angular_radius # play time
star.update()
photosphere.integrate(energies=signal.energies, threads=3)
temp = np.sum(photosphere.signal[0][0] + photosphere.signal[1][0], axis=0)
_ = ax.plot(hot.phases_in_cycles[0], temp/np.max(temp), 'k-', lw=0.5)
likelihood['cos_inclination'] = math.cos(1.0)
for angular_radius in angular_radii:
likelihood['p__super_radius'] = angular_radius
star.update()
photosphere.integrate(energies=signal.energies, threads=3)
temp = np.sum(photosphere.signal[0][0] + photosphere.signal[1][0], axis=0)
_ = ax.plot(hot.phases_in_cycles[0], temp/np.max(temp), 'r-', lw=0.5)
likelihood['cos_inclination'] = math.cos(0.5)
for angular_radius in angular_radii:
likelihood['p__super_radius'] = angular_radius
star.update()
photosphere.integrate(energies=signal.energies, threads=3)
temp = np.sum(photosphere.signal[0][0] + photosphere.signal[1][0], axis=0)
_ = ax.plot(hot.phases_in_cycles[0], temp/np.max(temp), 'b-', lw=0.5)
veneer((0.05,0.2), (0.05,0.2), ax)
Let us now construct a callable object representing a joint prior density distribution on the space $\mathbb{R}^{d}$. We need to extend the base class to implement our distribution, which with respect to some parameters is separable, but for others it is uniform on a joint space, and compactly supported according to non-trivial constraint equations.
As an example gravitational mass and equatorial radius: a joint constraint is imposed to assign zero prior density to stars which are too compact: the polar radius, in units of the gravitational radius, of the rotationally deformed stellar 2-surface is too small.
from scipy.stats import truncnorm
class CustomPrior(xpsi.Prior):
""" A custom (joint) prior distribution.
Source: Fictitious
Model variant: ST-U
Two single-temperature, simply-connected circular hot regions with
unshared parameters.
"""
__derived_names__ = ['compactness', 'phase_separation',]
__draws_from_support__ = 3
def __init__(self):
""" Nothing to be done.
A direct reference to the spacetime object could be put here
for use in __call__:
.. code-block::
self.spacetime = ref
Instead we get a reference to the spacetime object through the
a reference to a likelihood object which encapsulates a
reference to the spacetime object.
"""
super(CustomPrior, self).__init__() # not strictly required if no hyperparameters
def __call__(self, p = None):
""" Evaluate distribution at ``p``.
:param list p: Model parameter values.
:returns: Logarithm of the distribution evaluated at ``p``.
"""
temp = super(CustomPrior, self).__call__(p)
if not np.isfinite(temp):
return temp
# based on contemporary EOS theory
if not self.parameters['radius'] <= 16.0:
return -np.inf
ref = self.parameters.star.spacetime # shortcut
# limit polar radius to try to exclude deflections >= \pi radians
# due to oblateness this does not quite eliminate all configurations
# with deflections >= \pi radians
R_p = 1.0 + ref.epsilon * (-0.788 + 1.030 * ref.zeta)
if R_p < 1.76 / ref.R_r_s:
return -np.inf
# polar radius at photon sphere for ~static star (static ambient spacetime)
#if R_p < 1.5 / ref.R_r_s:
# return -np.inf
mu = math.sqrt(-1.0 / (3.0 * ref.epsilon * (-0.788 + 1.030 * ref.zeta)))
# 2-surface cross-section have a single maximum in |z|
# i.e., an elliptical surface; minor effect on support, if any,
# for high spin frequenies
if mu < 1.0:
return -np.inf
ref = self.parameters # redefine shortcut
# enforce order in hot region colatitude
if ref['p__super_colatitude'] > ref['s__super_colatitude']:
return -np.inf
phi = (ref['p__phase_shift'] - 0.5 - ref['s__phase_shift']) * _2pi
ang_sep = xpsi.HotRegion.psi(ref['s__super_colatitude'],
phi,
ref['p__super_colatitude'])
# hot regions cannot overlap
if ang_sep < ref['p__super_radius'] + ref['s__super_radius']:
return -np.inf
return 0.0
def inverse_sample(self, hypercube=None):
""" Draw sample uniformly from the distribution via inverse sampling. """
to_cache = self.parameters.vector
if hypercube is None:
hypercube = np.random.rand(len(self))
# the base method is useful, so to avoid writing that code again:
_ = super(CustomPrior, self).inverse_sample(hypercube)
ref = self.parameters # shortcut
idx = ref.index('distance')
ref['distance'] = truncnorm.ppf(hypercube[idx], -2.0, 7.0, loc=0.3, scale=0.1)
# flat priors in cosine of hot region centre colatitudes (isotropy)
# support modified by no-overlap rejection condition
idx = ref.index('p__super_colatitude')
a, b = ref.get_param('p__super_colatitude').bounds
a = math.cos(a); b = math.cos(b)
ref['p__super_colatitude'] = math.acos(b + (a - b) * hypercube[idx])
idx = ref.index('s__super_colatitude')
a, b = ref.get_param('s__super_colatitude').bounds
a = math.cos(a); b = math.cos(b)
ref['s__super_colatitude'] = math.acos(b + (a - b) * hypercube[idx])
# restore proper cache
for parameter, cache in zip(ref, to_cache):
parameter.cached = cache
# it is important that we return the desired vector because it is
# automatically written to disk by MultiNest and only by MultiNest
return self.parameters.vector
def transform(self, p, **kwargs):
""" A transformation for post-processing. """
p = list(p) # copy
# used ordered names and values
ref = dict(zip(self.parameters.names, p))
# compactness ratio M/R_eq
p += [gravradius(ref['mass']) / ref['radius']]
# phase separation between hot regions
# first some temporary variables:
if ref['p__phase_shift'] < 0.0:
temp_p = ref['p__phase_shift'] + 1.0
else:
temp_p = ref['p__phase_shift']
temp_s = 0.5 + ref['s__phase_shift']
if temp_s > 1.0:
temp_s = temp_s - 1.0
# now append:
if temp_s >= temp_p:
p += [temp_s - temp_p]
else:
p += [1.0 - temp_p + temp_s]
return p
We can now construct and instantiate a callable prior object, pass it to the likelihood object, which stores a mutual reference to itself in prior.
prior = CustomPrior()
No parameters supplied... empty subspace created.
likelihood.prior = prior
Let's draw samples from the prior and plot them:
samps, _ = prior.draw(int(1e4))
Drawing samples from the joint prior... Samples drawn.
def plot_samps(samps, x, y, xlabel, ylabel, s=1.0, color='k', **kwargs):
""" Plot samples as 2D scatter plot. """
fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111)
ax.scatter(samps[:,x], samps[:,y], s=s, color=color, **kwargs)
veneer(None, None, ax)
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
return plt.gca()
Let's first plot the $(M, R_{\rm eq})$ samples:
ax = plot_samps(samps,
likelihood.index('radius'), # we can find the index for __getitem__ magic
likelihood.index('mass'),
likelihood.get_param('radius').symbol + ' [km]', # can use symbols like so
r'$M$ [M$_{\odot}$]') # or manual entry
# R_eq = 1.76 x r_s(M)
_ = ax.plot(2.0*1.76*gravradius(np.linspace(1.0,3.0,100)), np.linspace(1.0,3.0,100), 'r--')
Note that the prior support is defined with a constraint that the polar radius $R_{\rm p}(R_{\rm eq}, M, \Omega)\geq 1.76r_{s}(M)$, which is why there is a region devoid of samples between the prior support and the dashed line $R_{\rm eq} = 1.76r_s(M)$. This condition is not necessary though because as of v0.6, X-PSI can include higher-order images, with the secondary and tertiary expected to be entirely sufficient in practice.
Let's now plot the hot region (circular spot) colatitudes:
ax = plot_samps(samps,
likelihood.index('p__super_colatitude'),
likelihood.index('s__super_colatitude'),
r'$\Theta_{p}$ [radians]', r'$\Theta_{s}$ [radians]')
# enforce colatitude order to distinguish hot regions as primary and secondary
_ = ax.plot(np.array([0.0,math.pi]), np.array([0.0,math.pi]), 'r--')
Note that the samples, marginalised over other region geometry parameters, are sparser when both hot regions approach the poles because we exclude overlapping configurations from the prior support. This is because the hot regions are by convention defined as disjoint, and cannot merge. If one wanted a more complex hot region, one would not invoke multiple hot regions that are permitted to overlap, but one would instead handle the extra complexity within the HotRegion class or a subclass.
Let's plot the angular radii of the spots:
_ = plot_samps(samps,
likelihood.index('p__super_radius'),
likelihood.index('s__super_radius'),
r'$\zeta_{p}$ [radians]',
r'$\zeta_{s}$ [radians]')
Note that the prior density is greater for hot regions that subtend smaller solid angles at the centre of the star, which also derives from the non-overlapping criterion for prior support.
Finally, let's take a look at the phases:
_ = plot_samps(samps,
likelihood.index('p__phase_shift'),
likelihood.index('s__phase_shift'),
r'$\phi_{p}$ [cycles]',
r'$\phi_{s}$ [cycles]')
Note that again because the hot regions cannot overlap, rarefaction occurs in the vicinity of lines of minimal phase separation. Note that the boundaries are all periodic, so this pattern tesselates. Because we implemented a transformation in our CustomPrior subclass, we can actually draw the samples and transform them, which is useful in post-processing contexts. We defined the intervals [-0.25, 0.75] for the inverse sampling so that the posterior mode(s) will not be near a boundary. The nested sampling algorithm can handle periodic boundaries by defining wrapped parameters; however, this can be trivially avoided altogether by rough inspection of the phases of the subpulses in the data, which we can see above are at around $-0.1$ and $0.4$ given the respective ground truth (injected) phases of $\phi_{p}=0.0$ and $\phi_{s}=0.025$.
Transformations for the purpose of likelihood evaluation must be handled in the inverse_sample method of an instance of the Prior class, but additional transformations that extend the parameter vector are written in the transform method.
If we wanted to transform automatically upon drawing the samples, thereby extending the parameter vectors passed to the __call__ method (so be careful with wrap-around indices when evaluating prior support conditions), we would do the following:
samps_plus_transformed, _ = prior.draw(int(1e4), transform=True)
Drawing samples from the joint prior... Samples drawn.
We defined a transformation from the hot region centre phases to the phase separation:
_ = plot_samps(samps_plus_transformed,
likelihood.index('p__phase_shift'),
likelihood.prior.index('phase_separation'),
r'$\phi_{p}$ [cycles]',
r'$\Delta\phi$ [cycles]')
We can see the rarefaction occurs for $\Delta\phi\sim0.0=1.0$.
The marginal one-dimensional prior distributions are overplotted, by the PostProcessing module, with the posterior distributions.
It is recommended to carefully inspect joint prior samples for pairs of parameters before commencing a sampling run, especially if there is a non-trivial constraint equation imposed on the prior support.
We have constructed and instantiated both a callable likelihood object and a callable prior object. We could proceed, for example, to apply the open-source sampler pymultinest or emcee to the joint posterior distribution proportional to the product of the (exponentiated) calls to the likelihood and prior objects.
We interface with the nested sampler MultiNest in a similar manner, by defining some runtime settings, and then passing those settings together with likelihood and prior objects to a wrapper from the Sample module. We will run the sampler for a specified number (1000) of nested replacements (iterations).
The environment variable LD_LIBRARY_PATH must be set before launching Jupyter as follows:
$ export LD_LIBRARY_PATH=<path/to/multinest>/libGet the parameters that are periodic or wrapped to inform MultiNest of boundary properties:
wrapped_params = [0]*len(likelihood)
wrapped_params[likelihood.index('p__phase_shift')] = 1
wrapped_params[likelihood.index('s__phase_shift')] = 1
runtime_params = {'resume': False,
'importance_nested_sampling': False,
'multimodal': False,
'n_clustering_params': None,
'outputfiles_basename': './run/run', # make ./run directory manually
'n_iter_before_update': 50,
'n_live_points': 50, #100
'sampling_efficiency': 0.8,
'const_efficiency_mode': False,
'wrapped_params': wrapped_params,
'evidence_tolerance': 0.5,
'max_iter': 100, # 10000 # manual termination condition for short test
'verbose': True}
for h in hot.objects:
h.set_phases(num_leaves = 100)
likelihood.threads = 3
likelihood.reinitialise()
likelihood.clear_cache()
# inform source code that parameter objects updated when inverse sampling
likelihood.externally_updated = True
p = [1.4,
12.5,
0.2,
math.cos(1.25),
0.0,
1.0,
0.075,
6.2,
0.025,
math.pi - 1.0,
0.2,
6.0]
# let's require that checks pass before starting to sample
check_kwargs = dict(hypercube_points = None,
physical_points = p, # externally_updated preserved
loglikelihood_call_vals = [-26713.6136777], # from above
rtol_loglike = 1.0e-6) # choose a tolerance
# note that mutual refs are already stored in the likelihood and prior
# objects to facilitate communication externally of the sampling process
xpsi.Sample.nested(likelihood, prior, check_kwargs, **runtime_params)
Imported PyMultiNest. Checking likelihood and prior evaluation before commencing sampling... Not using ``allclose`` function from NumPy. Using fallback implementation instead. Checking closeness of likelihood arrays: -2.6713601352e+04 | -2.6713613678e+04 ..... Closeness evaluated. Log-likelihood value checks passed on root process. Checks passed. Commencing integration... Estimating fractional hypervolume of the unit hypercube with finite prior density: Requiring 1E+03 draws from the prior support for Monte Carlo estimation... Drawing samples from the joint prior... Samples drawn. The support occupies an estimated 11.4% of the hypervolume within the unit hypercube... Fractional hypervolume estimated. Sampling efficiency set to: 7.0200. ***************************************************** MultiNest v3.12 Copyright Farhan Feroz & Mike Hobson Release Nov 2019 no. of live points = 50 dimensionality = 12 ***************************************************** Starting MultiNest generating live points live points generated, starting sampling Acceptance Rate: 0.847458 Replacements: 100 Total Samples: 118 Nested Sampling ln(Z): -120724.880426 Acceptance Rate: 0.600000 Replacements: 150 Total Samples: 250 Nested Sampling ln(Z): -115567.940017 analysing data from ./run/run.txt ln(ev)= -115567.94001695007 +/- NaN Total Likelihood Evaluations: 250 Sampling finished. Exiting MultiNest Integration completed.
The verbose output of the MultiNest program is by default directed to the host terminal session. Instead of trying to redirect that output to that of the above cell, we simply copy and paste the output from the terminal below. For the purpose of illustration, this output was for a 12-dimensional model variant (as for the emcee plot above but using a longer test run).
*****************************************************
MultiNest v3.11
Copyright Farhan Feroz & Mike Hobson
Release Apr 2018
no. of live points = 100
dimensionality = 12
*****************************************************
Starting MultiNest
generating live points
live points generated, starting sampling
Acceptance Rate: 0.724638
Replacements: 100
Total Samples: 138
Nested Sampling ln(Z): **************
Acceptance Rate: 0.649351
Replacements: 150
Total Samples: 231
Nested Sampling ln(Z): -116670.287917
Acceptance Rate: 0.569801
Replacements: 200
Total Samples: 351
Nested Sampling ln(Z): -115291.669431
Acceptance Rate: 0.449640
Replacements: 250
Total Samples: 556
Nested Sampling ln(Z): -108499.449911
Acceptance Rate: 0.408719
Replacements: 300
Total Samples: 734
Nested Sampling ln(Z): -95430.022790
Acceptance Rate: 0.367261
Replacements: 350
Total Samples: 953
Nested Sampling ln(Z): -77360.112633
Acceptance Rate: 0.319744
Replacements: 400
Total Samples: 1251
Nested Sampling ln(Z): -66119.380404
Acceptance Rate: 0.263930
Replacements: 450
Total Samples: 1705
Nested Sampling ln(Z): -57607.930990
Acceptance Rate: 0.213675
Replacements: 500
Total Samples: 2340
Nested Sampling ln(Z): -53505.956949
Acceptance Rate: 0.173119
Replacements: 550
Total Samples: 3177
Nested Sampling ln(Z): -50428.177797
Acceptance Rate: 0.147893
Replacements: 600
Total Samples: 4057
Nested Sampling ln(Z): -47108.755667
Acceptance Rate: 0.132653
Replacements: 650
Total Samples: 4900
Nested Sampling ln(Z): -43437.007007
Acceptance Rate: 0.125381
Replacements: 700
Total Samples: 5583
Nested Sampling ln(Z): -39888.092691
Acceptance Rate: 0.113533
Replacements: 750
Total Samples: 6606
Nested Sampling ln(Z): -36841.337131
Acceptance Rate: 0.100251
Replacements: 800
Total Samples: 7980
Nested Sampling ln(Z): -34450.919514
Acceptance Rate: 0.088450
Replacements: 850
Total Samples: 9610
Nested Sampling ln(Z): -32545.531967
Acceptance Rate: 0.080121
Replacements: 900
Total Samples: 11233
Nested Sampling ln(Z): -31270.147897
Acceptance Rate: 0.069674
Replacements: 950
Total Samples: 13635
Nested Sampling ln(Z): -30103.155016
Acceptance Rate: 0.064201
Replacements: 1000
Total Samples: 15576
Nested Sampling ln(Z): -29365.169148
Acceptance Rate: 0.058427
Replacements: 1050
Total Samples: 17971
Nested Sampling ln(Z): -28879.280235
ln(ev)= -28879.280235090871 +/- NaN
Total Likelihood Evaluations: 17971
Sampling finished. Exiting MultiNestTo prove that the objects constructed above can be fed to the emcee sampler, let's run a number of iterations using a single Python process. We will initialise the ensemble by drawing from a multivariate Gaussian with mean vector equal to the ground truth (injected) vector.
std = [1.0,
0.05,
1.0,
0.01,
0.05,
0.0025,
0.01,
0.05,
0.05,
0.01,
0.01,
0.05]
runtime_params = {'resume': False,
'root_dir': './',
'nwalkers': 50,
'nsteps': 100,
'walker_dist_moments': zip(p, std)} # if resume then ``None``
for h in hot.objects:
h.set_phases(num_leaves = 100) # revert to original phase resolution above
likelihood.threads = 3
likelihood.reinitialise() # because we played with the object above
# Use MPI=False for testing purposes
backend = xpsi.Sample.ensemble(likelihood, prior,
MPI=False, **runtime_params)
Imported emcee version: 3.1.4 Warning: a ``samples.h5`` file exists in ``./``. Attempting to move ``samples.h5`` to a subdirectory of ``./``. File archived in subdirectory ``./old_samples``. Commencing posterior sampling.
100%|█████████████████████████████████████████| 100/100 [15:49<00:00, 9.50s/it]
Sampling complete.
# clean up the docs/source directory
#!rm samples.h5; rm -r old_samples
Note that we could also try initialising the ensemble by inverse sampling the joint prior distribution.
Let's quickly plot the evolution of the ensemble Markov chains to prove that the sampling process commenced and is behaving in a somewhat reasonable manner:
try:
backend
except NameError:
import emcee
backend = emcee.backends.HDFBackend('samples.h5')
chain = backend.get_chain()
# These chains plotted were generated using
# v0.2 so the labels here are in a different
# order. The model also had two free
# temperature parameters instead of just one.
# Also, inclination was the parameter, not
# cos(inclination), and the prior was flat
# in inclination, not cos(inclination).
labels = [r'$D$',
r'$M$',
r'$R_{\rm eq}$',
r'$\cos(i)$',
r'$\phi_{p}$',
r'$\Theta_{p}$',
r'$\zeta_{p}$',
r'$T_{p}$',
r'$\phi_{s}$'
r'$\Theta_{s}$',
r'$\zeta_{s}$']
fig = plt.figure(figsize=(8,32))
gs = gridspec.GridSpec(12, 1, hspace=0.15)
for i in range(len(labels)):
ax = plt.subplot(gs[i,0])
ax.set_ylabel(labels[i])
for j in range(50):
plt.plot(chain[:,j,i], 'k-', lw=0.5, alpha=0.5)
if i < 11:
ax.tick_params(axis='x', labelbottom=False)
plt.setp(ax.get_yticklabels()[0], visible=False)
plt.setp(ax.get_yticklabels()[-1], visible=False)
else: ax.set_xlabel('Steps')
veneer((250, 1000), None, ax)
The chains rendered in the documentation were run on a desktop machine in about a day of wall-time. It is visually discernable that the ensemble distribution has not yet evolved to a stationary state: a rigourous application of ensemble MCMC would cover convergence criteria, auto-correlation, and examination of sensitivity to initial conditions and the transition kernel. In fact, based on the analysis with nested sampling on path xpsi/examples/default_background, we know that the posterior mode in the vicinity of the above ensemble is rather non-linear in the space being sampled, so ensemble MCMC may require many steps in order to argue for convergence.
In this notebook thus far we have not generated sythetic data. However, we did condition on synthetic data. Below we outline how that data was generated.
The background radiation field incident on the model instrument for the purpose of generating synthetic data was a time-invariant powerlaw spectrum, and was transformed into a count-rate in each output channel using the response matrix for synthetic data generation. We would reproduce this background here by writing a custom subclass as follows.
class CustomBackground(xpsi.Background):
""" The background injected to generate synthetic data. """
def __init__(self, bounds=None, value=None):
# first the parameters that are fundemental to this class
doc = """
Powerlaw spectral index.
"""
index = xpsi.Parameter('powerlaw_index',
strict_bounds = (-3.0, -1.01),
bounds = bounds,
doc = doc,
symbol = r'$\Gamma$',
value = value)
super(CustomBackground, self).__init__(index)
def __call__(self, energy_edges, phases):
""" Evaluate the incident background field. """
G = self['powerlaw_index']
temp = np.zeros((energy_edges.shape[0] - 1, phases.shape[0]))
temp[:,0] = (energy_edges[1:]**(G + 1.0) - energy_edges[:-1]**(G + 1.0)) / (G + 1.0)
for i in range(phases.shape[0]):
temp[:,i] = temp[:,0]
self._incident_background = temp
Note that the analytic background is integrated over energy intervals, as required by a Signal instance, which would then straightforwardly apply the model instrument response to the background.
We can now construct and instantiate a background object. We need to specify the number of background parameters, and define the hard bounds on those parameters; in this case we have only a single parameter, the powerlaw index.
We would then instantiate as follows:
background = CustomBackground(bounds=(None, None)) # use strict bounds, but do not fix/derive
Creating parameter:
> Named "powerlaw_index" with bounds [-3.000e+00, -1.010e+00].
> Powerlaw spectral index.
We are also in need of a simpler data object.
class SynthesiseData(xpsi.Data):
""" Custom data container to enable synthesis. """
def __init__(self, channels, phases, first, last):
self.channels = channels
self._phases = phases
try:
self._first = int(first)
self._last = int(last)
except TypeError:
raise TypeError('The first and last channels must be integers.')
if self._first >= self._last:
raise ValueError('The first channel number must be lower than the '
'the last channel number.')
Instantiate:
_data = SynthesiseData(np.arange(20,201), np.linspace(0.0, 1.0, 33), 0, 180)
Setting channels for event data... Channels set.
We are in need of a synthesise method, which in this implementation wraps an extension module. Let's check what the extension module offers:
from xpsi.tools.synthesise import synthesise_given_total_count_number as _synthesise
_synthesise?
When creating synthetic data, we can fix the random seed for noise by providing it as an input parameter called 'gsl_seed'. If no input is provided, the seed is based on the clock time.
def synthesise(self,
require_source_counts,
require_background_counts,
name='synthetic',
directory='./data',
**kwargs):
""" Synthesise data set.
"""
self._expected_counts, synthetic, _, _ = _synthesise(self._data.phases,
require_source_counts,
self._signals,
self._phases,
self._shifts,
require_background_counts,
self._background.registered_background,
gsl_seed=0) #seed=0
try:
if not os.path.isdir(directory):
os.mkdir(directory)
except OSError:
print('Cannot create write directory.')
raise
np.savetxt(os.path.join(directory, name+'_realisation.dat'),
synthetic,
fmt = '%u')
self._write(self.expected_counts,
filename = os.path.join(directory, name+'_expected_hreadable.dat'),
fmt = '%.8e')
self._write(synthetic,
filename = os.path.join(directory, name+'_realisation_hreadable.dat'),
fmt = '%u')
def _write(self, counts, filename, fmt):
""" Write to file in human readable format. """
rows = len(self._data.phases) - 1
rows *= len(self._data.channels)
phases = self._data.phases[:-1]
array = np.zeros((rows, 3))
for i in range(counts.shape[0]):
for j in range(counts.shape[1]):
array[i*len(phases) + j,:] = self._data.channels[i], phases[j], counts[i,j]
np.savetxt(filename, array, fmt=['%u', '%.6f'] + [fmt])
Add unbound methods:
CustomSignal.synthesise = synthesise
CustomSignal._write = _write
Instantiate, and reconfigure the likelihood object:
signal = CustomSignal(data = _data,
instrument = NICER,
background = background,
interstellar = None,
prefix='NICER')
for h in hot.objects:
h.set_phases(num_leaves = 100)
likelihood = xpsi.Likelihood(star = star, signals = signal, threads=1)
Creating parameter:
> Named "phase_shift" with fixed value 0.000e+00.
> The phase shift for the signal, a periodic parameter [cycles].
Warning: No data... can synthesise data but cannot evaluate a likelihood function.
Check write path:
!pwd
/home/tuomo/xpsi/xpsi_group/xpsi/docs/source
likelihood
Free parameters --------------- mass: Gravitational mass [solar masses]. radius: Coordinate equatorial radius [km]. distance: Earth distance [kpc]. cos_inclination: Cosine of Earth inclination to rotation axis. p__phase_shift: The phase of the hot region, a periodic parameter [cycles]. p__super_colatitude: The colatitude of the centre of the superseding region [radians]. p__super_radius: The angular radius of the (circular) superseding region [radians]. p__super_temperature: log10(superseding region effective temperature [K]). s__phase_shift: The phase of the hot region, a periodic parameter [cycles]. s__super_colatitude: The colatitude of the centre of the superseding region [radians]. s__super_radius: The angular radius of the (circular) superseding region [radians]. s__super_temperature: log10(superseding region effective temperature [K]). NICER__powerlaw_index: Powerlaw spectral index.
p = [1.4,
12.5,
0.2,
math.cos(1.25),
0.0,
1.0,
0.075,
6.2,
0.025,
math.pi - 1.0,
0.2,
6.0,
-2.0]
NICER_kwargs = dict(require_source_counts=2.0e6,
require_background_counts=2.0e6,
name='new_synthetic',
directory='./data')
likelihood.synthesise(p, force=True, NICER=NICER_kwargs)
Exposure time: 984282.608809 [s] Background normalisation: 1.89132786e-05
plot_one_pulse(np.loadtxt('data/new_synthetic_realisation.dat', dtype=np.double), _data.phases)
In this notebook we constructed a model including a likelihood and prior objects. We also looked at the sampling interface, and concluded by synthesising the pre-prepared data set that was loaded at the beginning of the notebook.