Welcome to enterprise’s documentation!

ENTERPRISE (Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE) is a pulsar timing analysis code, aimed at noise analysis, gravitational-wave searches, and timing model analysis.

Installation

From sources

The sources for enterprise can be downloaded from the Github repo.

You can either clone the public repository:

$ git clone git://github.com/nanograv/enterprise

Or download the tarball:

$ curl  -OL https://github.com/nanograv/enterprise/tarball/master

Once you have a copy of the source, you can install it with:

$ pip install numpy
$ pip install -r requirements.txt
$ pip install git+https://github.com/vallis/libstempo.git --install-option="--with-tempo2=$TEMPO2"
$ python setup.py install

If you want to run tests or do any other development then also run:

$ pip install -r requirements_dev.txt

Contributing

Contributions are welcome, and they are greatly appreciated! Every little bit helps, and credit will always be given.

You can contribute in many ways:

Types of Contributions

Report Bugs

Report bugs at https://github.com/nanograv/enterprise/issues.

If you are reporting a bug, please include:

• Your operating system name and version.

• Any details about your local setup that might be helpful in troubleshooting.

• Detailed steps to reproduce the bug.

Fix Bugs

Look through the GitHub issues for bugs. Anything tagged with “bug” and “help wanted” is open to whoever wants to implement it.

Implement Features

Look through the GitHub issues for features. Anything tagged with “enhancement” and “help wanted” is open to whoever wants to implement it.

Write Documentation

enterprise could always use more documentation, whether as part of the official enterprise docs, in docstrings, or even on the web in blog posts, articles, and such.

Submit Feedback

The best way to send feedback is to file an issue at https://github.com/nanograv/enterprise/issues.

If you are proposing a feature:

• Explain in detail how it would work.

• Keep the scope as narrow as possible, to make it easier to implement.

• Remember that this is a volunteer-driven project, and that contributions are welcome :)

Get Started!

Ready to contribute? Here’s how to set up enterprise for local development.

1. Fork the enterprise repo on GitHub.

2. Clone your fork locally:

   $ git clone git@github.com:your_name_here/enterprise.git
   

3. Set enterprise/master as upstream remote:

   $ git remote add upstream https://github.com/nanograv/enterprise.git
   

   You can then pull changes from the upstream master branch with:

   $ git pull upstream master
   

4. This is how you set up your fork for local development:

   Note

   You will need to have tempo and suitesparse installed before running the commands below.

$ cd enterprise/
$ make init
$ source .enterprise/bin/activate

5. Create a branch for local development:

   $ git checkout -b name-of-your-bugfix-or-feature
   

   Now you can make your changes locally.

6. When you’re done making changes, check that your changes pass flake8 and the tests, including testing other Python versions with tox (tox not implemented yet). Also check that any new docs are formatted correctly:

   $ make test
   $ make lint
   $ make docs
   

   To get flake8 and tox, just pip install them into your virtualenv.

7. Commit your changes and push your branch to GitHub:

   $ git add .
   $ git commit -m "Your detailed description of your changes."
   $ git push origin name-of-your-bugfix-or-feature
   

8. Submit a pull request through the GitHub website.

Pull Request Guidelines

Before you submit a pull request, check that it meets these guidelines:

1. The pull request should include tests.

2. If the pull request adds functionality, the docs should be updated. Put your new functionality into a function with a docstring.

3. The pull request should work for Python 2.6, 2.7, 3.3, 3.4 and 3.5, and for PyPy. Check https://travis-ci.org/nanograv/enterprise/pull_requests and make sure that the tests pass for all supported Python versions.

Tips

To run a subset of tests:

$ python -m unittest tests.test_enterprise

To track and checkout another user’s branch:

$ git remote add other-user-username https://github.com/other-user-username/enterprise.git
$ git fetch other-user-username
$ git checkout --track -b branch-name other-user-username/branch-name

Note

This tutorial was generated from a Jupyter notebook that can be downloaded here.

enterprise Data Structures

This guide will give an introduction to the unique data structures used in enterprise. These are all designed with the goal of making this code as user-friendly as possible, both for the end user and the developer.

Class Factories

The enterprise code makes heavy use of so-called class factories. Class factories are functions that return classes (not objects of class instances). A simple example is as follows:

def A(farg1, farg2):
    class A(object):

        def __init__(self, iarg):
            self.iarg = iarg

        def print_info(self):
            print('Object instance {}\nInstance argument: {}\nFunction args: {} {}\n'.format(
                self, self.iarg, farg1, farg2))
    return A

# define class A with arguments that can be seen within the class
a = A('arg1', 'arg2')

# instantiate 2 instances of class A with different arguments
a1 = a('iarg1')
a2 = a('iarg2')

# call print_info method
a1.print_info()
a2.print_info()

Object instance <__main__.A object at 0x10bb64290>
Instance argument: iarg1
Function args: arg1 arg2

Object instance <__main__.A object at 0x10bb642d0>
Instance argument: iarg2
Function args: arg1 arg2

In the example above we see that the arguments arg1 and arg2 are seen by both instances a1 and a2; however these instances were intantiated with different input arguments iarg1 and iarg2. So we see that class-factories are great when we want to give “global” parameters to a class without having to pass them on initialization. This also allows us to mix and match classes, as we will do in enterprise before we instantiate them.

The Pulsar class

The Pulsar class is a simple data structure that stores all of the important information about a pulsar that is obtained from a timing package such as the TOAs, residuals, error-bars, flags, design matrix, etc.

This class is instantiated with a par and a tim file. Full documentation on this class can be found here.

psr = Pulsar(datadir+'/B1855+09_NANOGrav_9yv1.gls.par', datadir+'/B1855+09_NANOGrav_9yv1.tim')

This Pulsar object is then passed to other enterprise data structures in a loosley coupled way in order to interact with the pulsar data.

The Parameter class

In enterprise signal parameters are set by specifying a prior distribution (i.e., Uniform, Normal, etc.). These Parameters are how enterprise builds signals. Below we will give an example of this functionality.

# lets define an efac parameter with a uniform prior from [0.5, 5]
efac = parameter.Uniform(0.5, 5)
print(efac)

<class 'enterprise.signals.parameter.Uniform'>

Uniform is a class factory that returns a class. The parameter is then intialized via a name. This way, a single parameter class can be initialized for multiple signal parameters with different names (i.e. EFAC per observing backend, etc). Once the parameter is initialized then you then have access to many useful methods.

# initialize efac parameter with name "efac_1"
efac1 = efac('efac_1')
print(efac1)

# return parameter name
print(efac1.name)

# get pdf at a point (log pdf is access)
print(efac1.get_pdf(1.3), efac1.get_logpdf(1.3))

# return 5 samples from this prior distribution
print(efac1.sample(n=5))

"efac_1":Uniform(0.5,5)
efac_1
(0.22222222222222221, -1.5040773967762742)
[ 4.15875031  4.02527174  0.86093696  2.29835222  3.14076572]

The Function structure

In enterprise we have defined a special data structure called Function. This data structure provides the user with a way to use and combine several different enterprise components in a user friendly way. More explicitly, it converts and standard function into an enterprise Function which can extract information from the Pulsar object and can also interact with enterprise Parameters.

[put reference to docstring here]

For example, consider the function:

@signal_base.function
def sine_wave(toas, log10_A=-7, log10_f=-8):
    return 10**log10_A * np.sin(2*np.pi*toas*10**log10_f)

Notice that the first positional argument of the function is toas, which happens to be a name of an attribute in the Pulsar class and the keyword arguments specify the default parameters for this function.

The decorator converts this standard function to a Function which can be used in two ways: the first way is to treat it like any other function.

# treat it just as a standard function with a vector input
sw = sine_wave(np.array([1,2,3]), log10_A=-8, log10_f=-7.5)
print(sw)

[  1.98691765e-15   3.97383531e-15   5.96075296e-15]

the second way is to use it as a Function:

# or use it as an enterprise function
sw_function = sine_wave(log10_A=parameter.Uniform(-10,-5), log10_f=parameter.Uniform(-9, -7))
print(sw_function)

<class 'enterprise.signals.signal_base.Function'>

Here we see that Function is actually a class factory, that is, when initialized with enterprise Parameters it returns a class that is initialized with a name and a Pulsar object as follows:

sw2 = sw_function('sine_wave', psr=psr)
print(sw2)

<enterprise.signals.signal_base.Function object at 0x109da10d0>

Now this Function object carries around instances of the Parameter classes given above for this particular function and Pulsar

print(sw2.params)

["sine_wave_log10_A":Uniform(-10,-5), "sine_wave_log10_f":Uniform(-9,-7)]

Most importantly it can be called in three different ways: If given without parameters it will fall back on the defaults given in the original function definition

print(sw2())

[  5.97588901e-08   5.97588901e-08   5.97588901e-08 ...,  -5.80521219e-08
  -5.80521219e-08  -5.80521219e-08]

or we can give it new fixed parameters

print(sw2(log10_A=-8, log10_f=-6.5))

[ -7.23515356e-09  -7.23515356e-09  -7.23515356e-09 ...,   5.93768399e-09
   5.93768399e-09   5.93768399e-09]

or most importantly we can give it a parameter dictionary with the Parameter names as keys. This is how Functions are use internally inside enterprise.

params = {'sine_wave_log10_A':-8, 'sine_wave_log10_f':-6.5}
print(sw2(params=params))

[ -7.23515356e-09  -7.23515356e-09  -7.23515356e-09 ...,   5.93768399e-09
   5.93768399e-09   5.93768399e-09]

Notice that the last two methods give the same answer since we gave it the same values just in different ways. So you may be thinking: “Why did we pass the Pulsar object on initialization?” or “Wait. How does it know about the toas?!”. Well the first question answers the second. By passing the pulsar object it grabs the toas attribute internally. This feature, combined with the ability to recognize Parameters and the ability to call the original function as we always would are the main strengths of Function, which is used heavily in enterprise.

Note that if we define a function without the decorator then we can still obtain a Function via:

def sine_wave(toas, log10_A=-7, log10_f=-8):
    return 10**log10_A * np.sin(2*np.pi*toas*10**log10_f)

sw3 = signal_base.Function(sine_wave, log10_A=parameter.Uniform(-10,-5),
                           log10_f=parameter.Uniform(-9, -7))

print(sw3)

<class 'enterprise.signals.signal_base.Function'>

Make your own Function

To define your own Function all you have to do is to define a function with these rules in mind.

1. If you want to use Pulsar attributes, define them as positional arguments with the same name as used in the Pulsar class (see here for more information.

2. Any arguments that you may use as Parameters must be keyword arguments (although you can have others that aren’t Parameters)

3. Add the @function decorator.

And thats it! You can now define your own Functions with minimal overhead and use them in enterprise or for tests and simulations or whatever you want!

The Selection structure

In the course of our analysis it is useful to split different signals into pieces. The most common flavor of this is to split the white noise parameters (i.e., EFAC, EQUAD, and ECORR) by observing backend system. The Selection structure is here to make this as smooth and versatile as possible.

The Selection structure is also a class-factory that returns a specific selection dictionary with keys and Boolean arrays as values.

This will become more clear with an example. Lets say that you want to split our parameters between the first and second half of the dataset, then we can define the following function:

def cut_half(toas):
    midpoint = (toas.max() + toas.min()) / 2
    return dict(zip(['t1', 't2'], [toas <= midpoint, toas > midpoint]))

This function will return a dictionary with keys (i.e. the names of the different subsections) t1 and t2 and boolean arrays corresponding to the first and second halves of the data span, respectively. So for a simple input we have:

toas = np.array([1,2,3,4])
print(cut_half(toas))

{'t2': array([False, False,  True,  True], dtype=bool), 't1': array([ True,  True, False, False], dtype=bool)}

To pass this to enterprise we turn it into a Selection via:

ch = Selection(cut_half)
print(ch)

<class 'enterprise.signals.selections.Selection'>

As we have stated, this is class factory that will be initialized inside enterprise signals with a Pulsar object in a very similar way to Functions.

ch1 = ch(psr)
print(ch1)
print(ch1.masks)

<enterprise.signals.selections.Selection object at 0x1048212d0>
{'t2': array([False, False, False, ...,  True,  True,  True], dtype=bool), 't1': array([ True,  True,  True, ..., False, False, False], dtype=bool)}

The Selection object has a method masks that uses the Pulsar object to evaluate the arguments of cut_half (these can be any number of Pulsar attributes, not just toas). The Selection object can also be called to return initialized Parameters with the split names as follows:

# make efac class factory
efac = parameter.Uniform(0.1, 5.0)

# now give it to selection
params, masks = ch1('efac', efac)

# named parameters
print(params)

# named masks
print(masks)

{u't1_efac': "B1855+09_t1_efac":Uniform(0.1,5.0), u't2_efac': "B1855+09_t2_efac":Uniform(0.1,5.0)}
{u't1_efac': array([ True,  True,  True, ..., False, False, False], dtype=bool), u't2_efac': array([False, False, False, ...,  True,  True,  True], dtype=bool)}

Make your own Selection

To define your own Selection all you have to do is to define a function with these rules in mind.

1. If you want to use Pulsar attributes, define them as positional arguments with the same name as used in the Pulsar class (see here for more information.

2. Make sure the return value is a dictionary with the names you want for the different segments and values as boolean arrays specifying which points to apply the split to.

3. A selection does not have to apply to all points. You can make it apply to only single points or single segments if you wish.

And thats it! You can now define your own Selections with minimal overhead and use them in enterprise or for tests and simulations or whatever you want!

Signals, SignalCollections, and PTAs oh my!

• The data (residuals) are modeled as the sum of Signal components which have their own Parameters.

• The sum of all Signal components is a SignalCollection.

• Each pulsar’s model is a SignalCollection that are combined to form a PTA.

• Common Signals are shared across pulsars

• Likelihoods act on PTAs.

Anatomy of an enterprise Signal

• \(\delta\tau = \sum_{i} X(\phi_{\rm basis})_{(i)}w_{(i)} + s(\phi_{\rm det}) + n(\phi_{\rm white})\)

• \(w_{(i)} | K_{(i)} = \mathrm{Normal}(0, K(\phi_{\rm gp})_{(i)})\)

class Signal(object):
    """Base class for Signal objects."""

    def get_ndiag(self, params):
        """Returns the diagonal of the white noise vector `N`.
        This method also supports block diagaonal sparse matrices.
        """
        return None

    def get_delay(self, params):
        """Returns the waveform of a deterministic signal."""
        return 0

    def get_basis(self, params=None):
        """Returns the basis array of shape N_toa x N_basis."""
        return None

    def get_phi(self, params):
        """Returns a diagonal or full rank covaraince matrix
        of the basis amplitudes."""
        return None

    def get_phiinv(self, params):
        """Returns inverse of the covaraince of basis amplitudes."""
        return None

Note

This tutorial was generated from a Jupyter notebook that can be downloaded here.

Analyzing Open MDC data

In this tutorial we will use enterprise to analyze open MDC dataset 1.

Get par and tim files

The first step in the process is getting the open MDC1 par and tim files in the tests directory.

parfiles = sorted(glob.glob(datadir + '/*.par'))
timfiles = sorted(glob.glob(datadir + '/*.tim'))

Load pulsars into Pulsar objects

enterprise uses a specific Pulsar object to store all of the relevant pulsar information (i.e. TOAs, residuals, error bars, flags, etc) from the timing package. Eventually enterprise will support both PINT and tempo2; however, for the moment it only supports tempo2 through the libstempo package. This object is then used to initalize Signals that define the generative model for the pulsar residuals. This is in keeping with the overall enterprise philosophy that the pulsar data should be as loosley coupled as possible to the pulsar model.

psrs = []
for p, t in zip(parfiles, timfiles):
    psr = Pulsar(p, t)
    psrs.append(psr)

Setup and run a simple noise model on a single pulsar

Here we will demonstrate how to do a simple noise run on a single pulsar. In this analysis we will simply model the noise via a single EFAC parameter and a power-law red noise process.

Set up model

Here we see the basic enterprise model building steps:

1. Define parameters and priors (This makes use of the Parameter class factory)

2. Set up the signals making use of the Signal class factories.

3. Define the model by summing the individual Signal classes.

4. Define a PTA by initializing the signal model with a Pulsar object.

Notice that powerlaw is uses as a Function here.

##### parameters and priors #####

# Uniform prior on EFAC
efac = parameter.Uniform(0.1, 5.0)

# red noise parameters
# Uniform in log10 Amplitude and in spectral index
log10_A = parameter.Uniform(-18,-12)
gamma = parameter.Uniform(0,7)

##### Set up signals #####

# white noise
ef = white_signals.MeasurementNoise(efac=efac)

# red noise (powerlaw with 30 frequencies)
pl = utils.powerlaw(log10_A=log10_A, gamma=gamma)
rn = gp_signals.FourierBasisGP(spectrum=pl, components=30)

# timing model
tm = gp_signals.TimingModel()

# full model is sum of components
model = ef + rn + tm

# initialize PTA
pta = signal_base.PTA([model(psrs[0])])

We can see which parameters we are going to be searching over with:

print(pta.params)

["J0030+0451_efac":Uniform(0.1,5.0), "J0030+0451_gamma":Uniform(0,7), "J0030+0451_log10_A":Uniform(-18,-12)]

Get initial parameters

We will start our MCMC chain at a random point in parameter space. We accomplish this by setting up a parameter dictionary using the name and sample methods for each Parameter.

xs = {par.name: par.sample() for par in pta.params}
print(xs)

{u'J0030+0451_efac': 4.7352650698633516, u'J0030+0451_gamma': 3.8216965873513029, u'J0030+0451_log10_A': -15.161366939011094}

Note that the rest of the analysis here is dependent on the sampling method and not on enterprise itself.

Set up sampler

Here we are making use of the PTMCMCSampler package for sampling. For this sampler, as in many others, it requires a function to compute the log-likelihood and log-prior given a vector of parameters. Here, these are supplied by PTA as pta.get_lnlikelihood and pta.get_lnprior.

# dimension of parameter space
ndim = len(xs)

# initial jump covariance matrix
cov = np.diag(np.ones(ndim) * 0.01**2)

# set up jump groups by red noise groups
ndim = len(xs)
groups  = [range(0, ndim)]
groups.extend([[1,2]])

# intialize sampler
sampler = ptmcmc(ndim, pta.get_lnlikelihood, pta.get_lnprior, cov, groups=groups,
                 outDir='chains/mdc/open1/')

Sample!

# sampler for N steps
N = 100000
x0 = np.hstack(p.sample() for p in pta.params)
sampler.sample(x0, N, SCAMweight=30, AMweight=15, DEweight=50)

Finished 10.00 percent in 7.578883 s Acceptance rate = 0.27876Adding DE jump with weight 50
Finished 99.00 percent in 77.849424 s Acceptance rate = 0.404505
Run Complete

Examine chain output

We see here that we have indeed recovered the injected values!

chain = np.loadtxt('chains/mdc/open1/chain_1.txt')
pars = sorted(xs.keys())
burn = int(0.25 * chain.shape[0])

truths = [1.0, 4.33, np.log10(5e-14)]
corner.corner(chain[burn:,:-4], 30, truths=truths, labels=pars);

Run full PTA GWB analysis

Here we will use the full 36 pulsar PTA to conduct a search for the GWB. In this analysis we fix the EFAC=1 for simplicity (and since we already know the answer!). This shows an example of how to use Constant parameters in enterprise.

Here you notice some of the simplicity of enterprise. For the most part, setting up the model for the full PTA is identical to that for one pulsar. In this case the only differences are that we are specifying the timespan to use when setting the GW and red noise frequencies and we are including a FourierBasisCommonGP signal, which models the GWB spectrum and spatial correlations.

After this setup, the rest is nearly identical to the single pulsar run above.

# find the maximum time span to set GW frequency sampling
tmin = [p.toas.min() for p in psrs]
tmax = [p.toas.max() for p in psrs]
Tspan = np.max(tmax) - np.min(tmin)

##### parameters and priors #####

# white noise parameters
# in this case we just set the value here since all efacs = 1
# for the MDC data
efac = parameter.Constant(1.0)

# red noise parameters
log10_A = parameter.Uniform(-18,-12)
gamma = parameter.Uniform(0,7)

##### Set up signals #####

# white noise
ef = white_signals.MeasurementNoise(efac=efac)

# red noise (powerlaw with 30 frequencies)
pl = utils.powerlaw(log10_A=log10_A, gamma=gamma)
rn = gp_signals.FourierBasisGP(spectrum=pl, components=30, Tspan=Tspan)

# gwb
# We pass this signal the power-law spectrum as well as the standard
# Hellings and Downs ORF
orf = utils.hd_orf()
crn = gp_signals.FourierBasisCommonGP(pl, orf, components=30, name='gw', Tspan=Tspan)

# timing model
tm = gp_signals.TimingModel()

# full model is sum of components
model = ef + rn + tm  + crn

# initialize PTA
pta = signal_base.PTA([model(psr) for psr in psrs])

Set up sampler

# initial parameters
xs = {par.name: par.sample() for par in pta.params}

# dimension of parameter space
ndim = len(xs)

# initial jump covariance matrix
cov = np.diag(np.ones(ndim) * 0.01**2)

# set up jump groups by red noise groups
ndim = len(xs)
groups  = [range(0, ndim)]
groups.extend(map(list, zip(range(0,ndim,2), range(1,ndim,2))))

sampler = ptmcmc(ndim, pta.get_lnlikelihood, pta.get_lnprior, cov, groups=groups,
                 outDir='chains/mdc/open1_gwb/')

# sampler for N steps
N = 100000
x0 = np.hstack(p.sample() for p in pta.params)
sampler.sample(x0, N, SCAMweight=30, AMweight=15, DEweight=50)

Plot output

chain = np.loadtxt('chains/mdc/open1_gwb/chain_1.txt')
pars = sorted(xs.keys())
burn = int(0.25 * chain.shape[0])

corner.corner(chain[burn:,-6:-4], 40, labels=pars[-2:], smooth=True, truths=[4.33, np.log10(5e-14)]);

Note

This tutorial was generated from a Jupyter notebook that can be downloaded here.

Analyzing 9-year NANOGrav data

In this tutorial we will use enterprise to analyze the NANOGrav 9-year data release for a stochastic GW background. We will reproduce the power-law GWB limit from this paper.

Function to convert PAL2 noise parameters to enterprise parameter dict

def get_noise_from_pal2(noisefile):
    psrname = noisefile.split('/')[-1].split('_noise.txt')[0]
    fin = open(noisefile, 'r')
    lines = fin.readlines()
    params = {}
    for line in lines:
        ln = line.split()
        if 'efac' in line:
            par = 'efac'
            flag = ln[0].split('efac-')[-1]
        elif 'equad' in line:
            par = 'log10_equad'
            flag = ln[0].split('equad-')[-1]
        elif 'jitter_q' in line:
            par = 'log10_ecorr'
            flag = ln[0].split('jitter_q-')[-1]
        elif 'RN-Amplitude' in line:
            par = 'log10_A'
            flag = ''
        elif 'RN-spectral-index' in line:
            par = 'gamma'
            flag = ''
        else:
            break
        if flag:
            name = [psrname, flag, par]
        else:
            name = [psrname, par]
        pname = '_'.join(name)
        params.update({pname: float(ln[1])})
    return params

Get par, tim, and noise files

Here we collect the tim and par files as well as noise files made from the PAL2 code. These are the same par, tim, and noise files used in the 9-year analysis papers. We use the convienience function above to convert from PAL2 noise files to enterprise parameter dictionaries.

parfiles = sorted(glob.glob(datadir + '/*.par'))
timfiles = sorted(glob.glob(datadir + '/*.tim'))
noisefiles = sorted(glob.glob(datadir + '/*noise.txt'))

# 18 pulsars used in 9 year analysis
p9 = np.loadtxt(datadir+'/9yr_pulsars.txt', dtype='S42')

# filter
parfiles = [x for x in parfiles if x.split('/')[-1].split('_')[0] in p9]
timfiles = [x for x in timfiles if x.split('/')[-1].split('_')[0] in p9]
noisefiles = [x for x in noisefiles if x.split('/')[-1].split('_')[0] in p9]

Load into Pulsar class list

psrs = []
for p, t in zip(parfiles, timfiles):
    psr = Pulsar(p, t, ephem='DE421')
    psrs.append(psr)

Get parameter dict from noisefiles

params = {}
for nfile in noisefiles:
    params.update(get_noise_from_pal2(nfile))

Set up model

When setting up the model for our upper limit run we fix all of the white noise (EFAC, EQUAD, and ECORR) parameters to the values obtained from the noise files. This is done by using Constant parameters. In this case we do not specify a default value for all instances of that parameter but instead will set them, based on their initialized pulsar and backend specific name, later via the set_default_params method of PTA.

Speaking of white noise parameters here, we also use the Selection object.

Another feature to notice is that we do not use a uniform prior on the log of the red noise or GWB amplitude. Instead we use a LinearExp prior (short for linear-exponent prior), that is a prior of the form \(p(x)\propto 10^x\). This is how we can still use the log of the parameter to sample but place a uniform prior on the parameter itself. We do this for both the red noise and GWB amplitude parameters.

Next, in order to save on computing time we do not include spatial correlations here. Instead we model the GWB as a common red process across all pulsars. In enterprise we can do this with a simple trick. We pre-initialize the parameters before passing them to the Signal model. In this way the same parameter instance is used for all pulsars. Lastly, we fixt the spectral index of the GWB to be 13/3 (4.33) using the Constant parameter.

# find the maximum time span to set GW frequency sampling
tmin = [p.toas.min() for p in psrs]
tmax = [p.toas.max() for p in psrs]
Tspan = np.max(tmax) - np.min(tmin)

# selection class to break white noise by backend
selection = selections.Selection(selections.by_backend)

##### parameters and priors #####

# white noise parameters
# since we are fixing these to values from the noise file we set
# them as constant parameters
efac = parameter.Constant()
equad = parameter.Constant()
ecorr = parameter.Constant()

# red noise parameters
log10_A = parameter.LinearExp(-20,-12)
gamma = parameter.Uniform(0,7)

# GW parameters (initialize with names here to use parameters in common across pulsars)
log10_A_gw = parameter.LinearExp(-18,-12)('log10_A_gw')
gamma_gw = parameter.Constant(4.33)('gamma_gw')

##### Set up signals #####

# white noise
ef = white_signals.MeasurementNoise(efac=efac, selection=selection)
eq = white_signals.EquadNoise(log10_equad=equad, selection=selection)
ec = white_signals.EcorrKernelNoise(log10_ecorr=ecorr, selection=selection)

# red noise (powerlaw with 30 frequencies)
pl = utils.powerlaw(log10_A=log10_A, gamma=gamma)
rn = gp_signals.FourierBasisGP(spectrum=pl, components=30, Tspan=Tspan)

# gwb (no spatial correlations)
cpl = utils.powerlaw(log10_A=log10_A_gw, gamma=gamma_gw)
gw = gp_signals.FourierBasisGP(spectrum=cpl, components=30, Tspan=Tspan)

# for spatial correltions you can do...
#orf = utils.hd_orf()
#crn = gp_signals.FourierBasisCommonGP(cpl, orf, components=30, name='gw', Tspan=Tspan)

# timing model
tm = gp_signals.TimingModel()

# to add solar system ephemeris modeling...
#eph = deterministic_signals.PhysicalEphemerisSignal(use_epoch_toas=True)

# full model is sum of components
model = ef + eq + ec + rn + tm + gw

# intialize PTA
pta = signal_base.PTA([model(psr) for psr in psrs])

Set white noise parameters

pta.set_default_params(params)

Set initial parameters drawn from prior and evaluate likelihood to fill caches

Evaluating the likelihood is not necessary, the caches will be filled the first time it is called within the sampler if not called here.

xs = {par.name: par.sample() for par in pta.params}
print pta.get_lnlikelihood(xs);
print pta.get_lnprior(xs);

1396202.32558
-32.2501201076

Set up sampler

# dimension of parameter space
ndim = len(xs)

# initial jump covariance matrix
cov = np.diag(np.ones(ndim) * 0.01**2)

# set up jump groups by red noise groups
ndim = len(xs)
groups  = [range(0, ndim)]
groups.extend(map(list, zip(range(0,ndim,2), range(1,ndim,2))))
groups.extend([[36]])

sampler = ptmcmc(ndim, pta.get_lnlikelihood, pta.get_lnprior, cov, groups=groups,
                 outDir='chains/nano_9_gwb/')

Sample!

# sampler for N steps
N = 1000000
x0 = np.hstack(p.sample() for p in pta.params)
sampler.sample(x0, N, SCAMweight=30, AMweight=15, DEweight=50, )

Plot output

chain = np.loadtxt('chains/nano_9_gwb/chain_1.txt)
pars = sorted(xs.keys())
burn = int(0.25 * chain.shape[0])

plt.hist(chain[burn:,-5], 50, normed=True, histtype='step', lw=2);
plt.xlabel(pars[-1]);

Upper limit value

We see that the upper limit agrees perfectly with the published value.

upper = 10**np.percentile(chain[burn:, -5], q=0.95)
print(upper)

1.49899289556e-15

enterprise package

Subpackages

enterprise.signals package

Submodules

enterprise.signals.anis_coefficients module

enterprise.signals.anis_coefficients.almFromClm(clm)

Given an array of clm values, return an array of complex alm valuex

Note: There is a bug in healpy for the negative m values. This function just takes the imaginary part of the abs(m) alm index.

enterprise.signals.anis_coefficients.anis_basis(psr_locs, lmax, nside=32)

Calculate the correlation basis matrices using the pixel-space transormations

@param psr_locs: Location of the pulsars [phi, theta] @param lmax: Maximum l to go up to @param nside: What nside to use in the pixelation [32]

Note: GW directions are in direction of GW propagation

enterprise.signals.anis_coefficients.clmFromAlm(alm)

Given an array of clm values, return an array of complex alm valuex

Note: There is a bug in healpy for the negative m values. This function just takes the imaginary part of the abs(m) alm index.

enterprise.signals.anis_coefficients.clmFromMap(h, lmax)

Given a pixel map, and a maximum l-value, return the corresponding C_{lm} values.

@param h: Sky power map @param lmax: Up to which order we’ll be expanding

return: clm values

Use real_sph_harm for the map

enterprise.signals.anis_coefficients.clmFromMap_fast(h, lmax)

Given a pixel map, and a maximum l-value, return the corresponding C_{lm} values.

@param h: Sky power map @param lmax: Up to which order we’ll be expanding

return: clm values

Use Healpix spherical harmonics for computational efficiency

enterprise.signals.anis_coefficients.createSignalResponse(pphi, ptheta, gwphi, gwtheta)

Create the signal response matrix. All parameters are assumed to be of the same dimensionality.

@param pphi: Phi of the pulsars @param ptheta: Theta of the pulsars @param gwphi: Phi of GW propagation direction @param gwtheta: Theta of GW propagation direction

@return: Signal response matrix of Earth-term

enterprise.signals.anis_coefficients.createSignalResponse_pol(pphi, ptheta, gwphi, gwtheta, plus=True, norm=True)

Create the signal response matrix. All parameters are assumed to be of the same dimensionality.

@param pphi: Phi of the pulsars @param ptheta: Theta of the pulsars @param gwphi: Phi of GW propagation direction @param gwtheta: Theta of GW propagation direction @param plus: Whether or not this is the plus-polarization @param norm: Normalise the correlations to equal Jenet et. al (2005)

@return: Signal response matrix of Earth-term

enterprise.signals.anis_coefficients.getCov(clm, nside, F_e)

Given a vector of clm values, construct the covariance matrix

@param clm: Array with Clm values @param nside: Healpix nside resolution @param F_e: Signal response matrix

@return: Cross-pulsar correlation for this array of clm values

enterprise.signals.anis_coefficients.mapFromClm(clm, nside)

Given an array of C_{lm} values, produce a pixel-power-map (non-Nested) for healpix pixelation with nside

@param clm: Array of C_{lm} values (inc. 0,0 element) @param nside: Nside of the healpix pixelation

return: Healpix pixels

Use real_sph_harm for the map

enterprise.signals.anis_coefficients.mapFromClm_fast(clm, nside)

Given an array of C_{lm} values, produce a pixel-power-map (non-Nested) for healpix pixelation with nside

@param clm: Array of C_{lm} values (inc. 0,0 element) @param nside: Nside of the healpix pixelation

return: Healpix pixels

Use Healpix spherical harmonics for computational efficiency

enterprise.signals.anis_coefficients.orfFromMap_fast(psr_locs, usermap, response=None)

Calculate an ORF from a user-defined sky map.

@param psr_locs: Location of the pulsars [phi, theta] @param usermap: Provide a healpix map for GW power

Note: GW directions are in direction of GW propagation

enterprise.signals.anis_coefficients.real_sph_harm(mm, ll, phi, theta)

The real-valued spherical harmonics.

enterprise.signals.anis_coefficients.signalResponse_fast(ptheta_a, pphi_a, gwtheta_a, gwphi_a)

Create the signal response matrix FAST

enterprise.signals.deterministic_signals module

Contains class factories for deterministic signals. Deterministic signals are defined as the class of signals that have a delay that is to be subtracted from the residuals.

enterprise.signals.deterministic_signals.Deterministic(waveform, selection=<class 'enterprise.signals.selections.Selection.<locals>.Selection'>, name='')

Class factory for generic deterministic signals.

enterprise.signals.deterministic_signals.PhysicalEphemerisSignal(frame_drift_rate=True, d_jupiter_mass=True, d_saturn_mass=True, d_uranus_mass=True, d_neptune_mass=True, jup_orb_elements=True, sat_orb_elements=False, model='setIII', use_epoch_toas=True, name='')

Class factory for physical ephemeris model signal.

This function implements a physically motivated ephemeris delay model. It is parameterized by an overall ecliptic-z frame drift rate, by the masses of gas giants, by the six orbital elements of Jupiter, and by the six orbital elements of Saturn. All these contributions can be disabled individually (Saturn orbit corrections are disabled by default).

Note

This signal is only compatible with a tempo2 Pulsar object.

The user can implement their own priors (e.g., by setting frame_drift_rate = parameter.Uniform(-1e-10,1e-10)(‘frame_drift_rate’) but we have set reasonable defaults (see below).

Parameters

• frame_drift_rate – ecliptic z-drift rate in units of rad/year referred to offset 1/1/2010. Default prior is Uniform(-1e-9, 1e-9).

• d_jupiter_mass – Mass deviation of Jupiter in solar masses. Default prior taken from IAU mass measurement uncertainty - Normal(0, 1.54976690e-11)

• d_saturn_mass – Mass deviation of Saturn in solar masses. Default prior taken from IAU mass measurement uncertainty - Normal(0, 8.17306184e-12)

• d_uranus_mass – Mass deviation of Uranus in solar masses. Default prior taken from IAU mass measurement uncertainty - Normal(0, 5.71923361e-11)

• d_neptune_mass – Mass deviation of Neptune in solar masses. Default prior taken from IAU mass measurement uncertainty - Normal(0, 7.96103855e-11)

• jup_orb_elements – Jupiter orbital-element perturbation coefficients. Default prior is Uniform(-0.05, 0.05) for each element, appropriate for standard SVD basis.

• sat_orb_elements – Saturn orbital-element perturbations coefficient. Default prior is Uniform(-0.5, 0.5) for each element, appropriate for standard SVD basis.

• model – Sets the vector basis used by Jupiter and Saturn orbital-element perturbations. Currently can be set to ‘orbel’, ‘orbel-v2’, ‘setIII’. Default: ‘setIII’

• use_epoch_toas – Use interpolation from epoch to full TOAs. This option reduces computational cost for large multi-channel TOA data sets. Default: True

enterprise.signals.gp_bases module

Utilities module containing various useful functions for use in other modules.

enterprise.signals.gp_bases.createfourierdesignmatrix_dm(toas, freqs, nmodes=30, Tspan=None, pshift=False, fref=1400, logf=False, fmin=None, fmax=None, modes=None)

Construct DM-variation fourier design matrix. Current normalization expresses DM signal as a deviation [seconds] at fref [MHz]

Parameters

• toas – vector of time series in seconds

• freqs – radio frequencies of observations [MHz]

• nmodes – number of fourier coefficients to use

• Tspan – option to some other Tspan

• pshift – option to add random phase shift

• fref – reference frequency [MHz]

• logf – use log frequency spacing

• fmin – lower sampling frequency

• fmax – upper sampling frequency

• modes – option to provide explicit list or array of sampling frequencies

Returns

F: DM-variation fourier design matrix

Returns

f: Sampling frequencies

enterprise.signals.gp_bases.createfourierdesignmatrix_env(toas, log10_Amp=-7, log10_Q=2.4771212547196626, t0=4579200000, nmodes=30, Tspan=None, logf=False, fmin=None, fmax=None, modes=None)

Construct fourier design matrix with gaussian envelope.

Parameters

• toas – vector of time series in seconds

• nmodes – number of fourier coefficients to use

• freqs – radio frequencies of observations [MHz]

• freq – option to output frequencies

• Tspan – option to some other Tspan

• logf – use log frequency spacing

• fmin – lower sampling frequency

• fmax – upper sampling frequency

• log10_Amp – log10 of the Amplitude [s]

• t0 – mean of gaussian envelope [s]

• log10_Q – log10 of standard deviation of gaussian envelope [days]

• modes – option to provide explicit list or array of sampling frequencies

Returns

F: fourier design matrix with gaussian envelope

Returns

f: Sampling frequencies

enterprise.signals.gp_bases.createfourierdesignmatrix_eph(t, nmodes, phi, theta, freq=False, Tspan=None, logf=False, fmin=None, fmax=None, modes=None)

enterprise.signals.gp_bases.createfourierdesignmatrix_ephem(toas, pos, nmodes=30, Tspan=None)

Construct ephemeris perturbation Fourier design matrix and frequencies. The matrix contains nmodes*6 columns, ordered as by frequency first, Cartesian coordinate second:

sin(f0) [x], sin(f0) [y], sin(f0) [z], cos(f0) [x], cos(f0) [y], cos(f0) [z], sin(f1) [x], sin(f1) [y], sin(f1) [z], …

The corresponding frequency vector repeats every entry six times. This design matrix should be used with monopole_orf and with a powerlaw that specifies components=6.

Parameters

• toas – vector of time series in seconds

• pos – pulsar position as Cartesian vector

• nmodes – number of Fourier coefficients

• Tspan – Tspan used to define Fourier bins

Returns

F: Fourier design matrix of shape (len(toas),6*nmodes)

Returns

f: Sampling frequencies (6*nmodes)

enterprise.signals.gp_bases.createfourierdesignmatrix_red(toas, nmodes=30, Tspan=None, logf=False, fmin=None, fmax=None, pshift=False, modes=None, pseed=None)

Construct fourier design matrix from eq 11 of Lentati et al, 2013 :param toas: vector of time series in seconds :param nmodes: number of fourier coefficients to use :param freq: option to output frequencies :param Tspan: option to some other Tspan :param logf: use log frequency spacing :param fmin: lower sampling frequency :param fmax: upper sampling frequency :param pshift: option to add random phase shift :param pseed: option to provide phase shift seed :param modes: option to provide explicit list or array of

sampling frequencies

Returns

F: fourier design matrix

Returns

f: Sampling frequencies

enterprise.signals.gp_priors module

Utilities module containing various useful functions for use in other modules.

enterprise.signals.gp_priors.InvGamma(alpha=1, gamma=1, size=None)

Class factory for Inverse Gamma parameters.

enterprise.signals.gp_priors.InvGammaPrior(value, alpha=1, gamma=1)

Prior function for InvGamma parameters.

enterprise.signals.gp_priors.InvGammaSampler(alpha=1, gamma=1, size=None)

Sampling function for Uniform parameters.

enterprise.signals.gp_priors.broken_powerlaw(f, log10_A, gamma, delta, log10_fb, kappa=0.1)

Generic broken powerlaw spectrum. :param f: sampling frequencies :param A: characteristic strain amplitude [set for gamma at f=1/yr] :param gamma: negative slope of PSD for f > f_break [set for comparison

at f=1/yr (default 13/3)]

Parameters

• delta – slope for frequencies < f_break

• log10_fb – log10 transition frequency at which slope switches from gamma to delta

• kappa – smoothness of transition (Default = 0.1)

enterprise.signals.gp_priors.free_spectrum(f, log10_rho=None)

Free spectral model. PSD amplitude at each frequency is a free parameter. Model is parameterized by S(f_i) =

ho_i^2 * T,

where

ho_i is the free parameter and T is the observation

length.

enterprise.signals.gp_priors.infinitepower(f)

enterprise.signals.gp_priors.powerlaw(f, log10_A=-16, gamma=5, components=2)

enterprise.signals.gp_priors.powerlaw_genmodes(f, log10_A=-16, gamma=5, components=2, wgts=None)

enterprise.signals.gp_priors.t_process(f, log10_A=-15, gamma=4.33, alphas=None)

t-process model. PSD amplitude at each frequency is a fuzzy power-law.

enterprise.signals.gp_priors.t_process_adapt(f, log10_A=-15, gamma=4.33, alphas_adapt=None, nfreq=None)

t-process model. PSD amplitude at each frequency is a fuzzy power-law.

enterprise.signals.gp_priors.turnover(f, log10_A=-15, gamma=4.33, lf0=-8.5, kappa=3.3333333333333335, beta=0.5)

enterprise.signals.gp_priors.turnover_knee(f, log10_A, gamma, lfb, lfk, kappa, delta)

Generic turnover spectrum with a high-frequency knee. :param f: sampling frequencies of GWB :param A: characteristic strain amplitude at f=1/yr :param gamma: negative slope of PSD around f=1/yr (usually 13/3) :param lfb: log10 transition frequency at which environment dominates GWs :param lfk: log10 knee frequency due to population finiteness :param kappa: smoothness of turnover (10/3 for 3-body stellar scattering) :param delta: slope at higher frequencies

enterprise.signals.gp_signals module

Contains class factories for Gaussian Process (GP) signals. GP signals are defined as the class of signals that have a basis function matrix and basis prior vector..

class enterprise.signals.gp_signals.MarginalizingNmat(Mmat, Nmat=0)

Bases: object

MNF(arg)

MNMMNF(arg)

MNr(arg)

solve(right, left_array=None, logdet=False)

property cf

enterprise.signals.gp_signals.BasisCommonGP(priorFunction, basisFunction, orfFunction, coefficients=False, combine=True, name='')

enterprise.signals.gp_signals.BasisGP(priorFunction, basisFunction, coefficients=False, combine=True, selection=<class 'enterprise.signals.selections.Selection.<locals>.Selection'>, name='')

Class factory for generic GPs with a basis matrix.

enterprise.signals.gp_signals.EcorrBasisModel(log10_ecorr=<class 'enterprise.signals.parameter.Uniform.<locals>.Uniform'>, coefficients=False, selection=<class 'enterprise.signals.selections.Selection.<locals>.Selection'>, name='basis_ecorr')

Convenience function to return a BasisGP class with a quantized ECORR basis.

enterprise.signals.gp_signals.FourierBasisCommonGP(spectrum, orf, coefficients=False, combine=True, components=20, Tspan=None, modes=None, name='common_fourier', pshift=False, pseed=None)

enterprise.signals.gp_signals.FourierBasisCommonGP_ephem(spectrum, components, Tspan, name='ephem_gp')

enterprise.signals.gp_signals.FourierBasisCommonGP_physicalephem(frame_drift_rate=1e-09, d_jupiter_mass=1.5497669e-11, d_saturn_mass=8.17306184e-12, d_uranus_mass=5.71923361e-11, d_neptune_mass=7.96103855e-11, jup_orb_elements=0.05, sat_orb_elements=0.5, model='setIII', coefficients=False, name='phys_ephem_gp')

Class factory for physical ephemeris corrections as a common GP. Individual perturbations can be excluded by setting the corresponding prior sigma to None.

Parameters

• frame_drift_rate – Gaussian sigma for frame drift rate

• d_jupiter_mass – Gaussian sigma for Jupiter mass perturbation

• d_saturn_mass – Gaussian sigma for Saturn mass perturbation

• d_uranus_mass – Gaussian sigma for Uranus mass perturbation

• d_neptune_mass – Gaussian sigma for Neptune mass perturbation

• jup_orb_elements – Gaussian sigma for Jupiter orbital elem. perturb.

• sat_orb_elements – Gaussian sigma for Saturn orbital elem. perturb.

• model – vector basis used by Jupiter and Saturn perturb.; see PhysicalEphemerisSignal, defaults to “setIII”

• coefficients – if True, treat GP coefficients as enterprise parameters; if False, marginalize over them

Returns

BasisCommonGP representing ephemeris perturbations

enterprise.signals.gp_signals.FourierBasisGP(spectrum, coefficients=False, combine=True, components=20, selection=<class 'enterprise.signals.selections.Selection.<locals>.Selection'>, Tspan=None, modes=None, name='red_noise', pshift=False, pseed=None)

Convenience function to return a BasisGP class with a fourier basis.

enterprise.signals.gp_signals.MarginalizingTimingModel(name='marginalizing_linear_timing_model', use_svd=False, normed=True)

Class factory for marginalizing (fast-likelihood) linear timing model signals.

enterprise.signals.gp_signals.TimingModel(coefficients=False, name='linear_timing_model', use_svd=False, normed=True)

Class factory for marginalized linear timing model signals.

enterprise.signals.gp_signals.WidebandTimingModel(dmefac=<class 'enterprise.signals.parameter.Uniform.<locals>.Uniform'>, log10_dmequad=<class 'enterprise.signals.parameter.Uniform.<locals>.Uniform'>, dmjump=<class 'enterprise.signals.parameter.Uniform.<locals>.Uniform'>, selection=<class 'enterprise.signals.selections.Selection.<locals>.Selection'>, dmjump_selection=<class 'enterprise.signals.selections.Selection.<locals>.Selection'>, dmjump_ref=None, name='wideband_timing_model')

Class factory for marginalized linear timing model signals that take wideband TOAs and DMs. Works in tandem with DMX parameters, by effectively setting their prior mean and variance. DM noise can be adjusted in analogy to the tempo/tempo2/PINT parameter convention, with variance = dmefac^2 (dmerr^2 + dmequad^2). DM noise can be segmented with selection (e.g., backends). DM jumps can also be added for each dmjump_selection; setting dmjump_ref to one of the selection labels keeps the corresponding jump to zero.

enterprise.signals.gp_signals.ecorr_basis_prior(weights, log10_ecorr=-8)

Returns the ecorr prior. :param weights: A vector or weights for the ecorr prior.

enterprise.signals.gp_signals.get_timing_model_basis(use_svd=False, normed=True)

enterprise.signals.parameter module

Contains parameter types for use in enterprise Signal classes.

class enterprise.signals.parameter.ConstantParameter(name)

Bases: object

Constant Parameter base class.

property value

class enterprise.signals.parameter.FunctionBase

Bases: object

class enterprise.signals.parameter.Parameter(name)

Bases: object

get_logpdf(value=None, **kwargs)

get_pdf(value=None, **kwargs)

sample(**kwargs)

property params

property size

enterprise.signals.parameter.Constant(val=None)

Class factory for Constant parameters. Leave val=None to set value later, for example with signal_base.PTA.set_default_params().

enterprise.signals.parameter.Function(func, name='', **func_kwargs)

enterprise.signals.parameter.GPCoefficients(logprior, size)

Class factory for GP coefficients, which are usually created inside gp_signals.BasisGP.

enterprise.signals.parameter.LinearExp(pmin, pmax, size=None)

Class factory for LinearExp parameters (with pdf(x) ~ 10^x, and 0 outside [pmin,``max``]). Handles vectors correctly if pmin and pmax are scalars or if size == len(pmin) == len(pmax)

Parameters

• pmin – minimum of range

• pmax – maximum of range

• size – length for vector parameter (default None)

Returns

LinearExp parameter class

enterprise.signals.parameter.LinearExpPrior(value, pmin, pmax)

Prior function for LinearExp parameters.

enterprise.signals.parameter.LinearExpSampler(pmin, pmax, size=None)

Sampling function for LinearExp parameters.

enterprise.signals.parameter.Normal(mu=0, sigma=1, size=None)

Class factory for Normal parameters (with pdf(x) ~ N(mu,``sigma``)). Handles vectors correctly if size == len(mu) == len(sigma), in which case sigma is taken as the sqrt of the diagonal of the covariance matrix; sigma can also be given passed as the size x size covariance matrix.

Parameters

• mu – center of normal distribution

• sigma – standard deviation of normal distribution

• size – length for vector parameter

Returns

Normal parameter class

enterprise.signals.parameter.NormalPrior(value, mu, sigma)

Prior function for Normal parameters. Handles scalar mu and sigma, compatible vector value/mu/sigma, vector value/mu and compatible covariance matrix sigma.

enterprise.signals.parameter.NormalSampler(mu, sigma, size=None)

Sampling function for Normal parameters. Handles scalar mu and sigma, compatible vector value/mu/sigma, vector value/mu and compatible covariance matrix sigma.

enterprise.signals.parameter.TruncNormal(mu=0, sigma=1, pmin=-2, pmax=2, size=None)

Class factory for TruncNormal parameters (with pdf(x) ~ N(mu,``sigma``) on a finite, truncated domain). Handles vectors correctly if size == len(mu) == len(sigma), in which case sigma is taken as the sqrt of the diagonal of the covariance matrix. :param mu: center of normal distribution :param sigma: standard deviation of normal distribution :param pmin: lower bound of domain :param pmax: upper bound of domain :param size: length for vector parameter :return: TruncNormal` parameter class

enterprise.signals.parameter.TruncNormalPrior(value, mu, sigma, pmin, pmax, norm=None)

Prior function for TruncNormal parameters. Handles scalar mu/sigma/pmin/pmax or compatible vector value/mu/sigma/pmin/pmax/.

enterprise.signals.parameter.TruncNormalSampler(mu, sigma, pmin, pmax, norm=None, size=None)

Sampling function for TruncNormal parameters. Handles scalar mu/sigma/pmin/pmax or compatible vector value/mu/sigma/pmin/pmax/. Implements rejection sampling on a normal distribution.

enterprise.signals.parameter.Uniform(pmin, pmax, size=None)

Class factory for Uniform parameters (with pdf(x) ~ 1/[pmax - pmin] inside [pmin,pmax], 0 outside. Handles vectors correctly, if pmin and pmax are scalars, or if len(size) == len(pmin) == len(pmax)

Parameters

• pmin – minimum of uniform range

• pmax – maximum of uniform range

• size – length for vector parameter

Returns

Uniform parameter class

enterprise.signals.parameter.UniformPrior(value, pmin, pmax)

Prior function for Uniform parameters. Accepts vector value and boundaries.

enterprise.signals.parameter.UniformSampler(pmin, pmax, size=None)

Sampling function for Uniform parameters.

enterprise.signals.parameter.UserParameter(prior=None, logprior=None, sampler=None, size=None)

Class factory for UserParameter, implementing Enterprise parameters with arbitrary priors. The prior is specified by way of an Enterprise Function of the form prior(value, [par1, par2]). Optionally, include sampler (a function with the same parameters as prior), to allow random sampling of the parameter through enterprise.signals.parameter.sample.

Parameters

• prior – parameter prior pdf, given as Enterprise Function

• sampler – function returning a randomly sampled parameter according to prior

• size – length for vector parameter

Returns

UserParameter class

enterprise.signals.parameter.function(func)

Decorator for Function.

enterprise.signals.parameter.get_funcargs(func)

Convenience function to get args and kwargs of any function.

enterprise.signals.parameter.sample(parlist)

Sample a list of Parameters consistently (i.e., keeping track of hyperparameters).

enterprise.signals.selections module

Contains various selection functions to mask parameters by backend flags, time-intervals, etc.

enterprise.signals.selections.Selection(func)

Class factory for residual selection.

enterprise.signals.selections.by_backend(backend_flags)

Selection function to split by backend flags.

enterprise.signals.selections.by_band(flags)

Selection function to split by PPTA frequency band under -B flag

enterprise.signals.selections.by_frontend(flags)

Selection function to split by frontend under -fe flag

enterprise.signals.selections.by_telescope(telescope)

Selection function to split by telescope

enterprise.signals.selections.call_me_maybe(obj)

See here for description.

enterprise.signals.selections.cut_half(toas)

Selection function to split by data segment

enterprise.signals.selections.nanograv_backends(backend_flags)

Selection function to split by NANOGRav backend flags only.

enterprise.signals.selections.no_selection(toas)

Default selection with no splitting.

enterprise.signals.selections.selection_func(func)

enterprise.signals.signal_base module

Defines the signal base classes and metaclasses. All signals will then be derived from these base classes.

class enterprise.signals.signal_base.BlockMatrix(blocks, slices, nvec=0)

Bases: object

solve(other, left_array=None, logdet=False)

class enterprise.signals.signal_base.CommonSignal(psr)

Bases: Signal

Base class for CommonSignal objects.

classmethod get_phicross(signal1, signal2, params)

get_phiinv(params)

Returns inverse of the covaraince of basis amplitudes.

class enterprise.signals.signal_base.LogLikelihood(pta, cholesky_sparse=True)

Bases: object

class enterprise.signals.signal_base.MetaCollection

Bases: type

Metaclass for Signal collections. Allows addition of SignalCollection classes.

class enterprise.signals.signal_base.MetaSignal

Bases: type

Metaclass for Signals. Allows addition of Signal classes.

class enterprise.signals.signal_base.PTA(init, lnlikelihood=<class 'enterprise.signals.signal_base.LogLikelihood'>)

Bases: object

get_TNT(params)

get_TNr(params)

get_basis(params={})

get_delay(params={})

get_lnlikelihood(params, **kwargs)

get_lnprior(params)

get_logsignalprior(params)

get_ndiag(params={})

get_phi(params, cliques=False)

get_phiinv(params, logdet=False, method='cliques')

get_phiinv_byfreq_cliques(params, logdet=False, cholesky=False)

get_phiinv_byfreq_partition(params, logdet=False)

get_phiinv_sparse(params, logdet=False)

get_rNr_logdet(params)

get_residuals()

get_signal(name)

Returns Signal instance given the signal name.

items()

keys()

map_params(xs)

set_default_params(params)

summary(include_params=True, to_stdout=False)

generate summary string for PTA model

Parameters

• include_params – [bool] list all parameters for each signal

• to_stdout – [bool] print summary to stdout instead of returning it

Returns

[string]

values()

property param_names

property params

property pulsarmodels

property pulsars

property signals

Return signal dictionary.

class enterprise.signals.signal_base.ShermanMorrison(jvec, slices, nvec=0.0)

Bases: object

Custom container class for Sherman-morrison array inversion.

solve(other, left_array=None, logdet=False)

class enterprise.signals.signal_base.Signal(psr)

Bases: object

Base class for Signal objects.

get(parname, params={})

get_basis(params=None)

Returns the basis array of shape N_toa x N_basis.

get_delay(params)

Returns the waveform of a deterministic signal.

get_logsignalprior(params)

Returns an additional prior/likelihood terms associated with a signal.

get_ndiag(params)

Returns the diagonal of the white noise vector N.

This method also supports block diagonal sparse matrices.

get_phi(params)

Returns a diagonal covariance matrix of the basis amplitudes.

get_phiinv(params)

Returns inverse of the covaraince of basis amplitudes.

set_default_params(params)

Set default parameters.

property param_names

property params

class enterprise.signals.signal_base.csc_matrix_alt(arg1, shape=None, dtype=None, copy=False)

Bases: csc_matrix

Sub-class of scipy.sparse.csc_matrix with custom add and solve methods.

solve(other, left_array=None, logdet=False)

class enterprise.signals.signal_base.ndarray_alt(inputarr)

Bases: ndarray

Sub-class of np.ndarray with custom solve method.

solve(other, left_array=None, logdet=False)

enterprise.signals.signal_base.LogLikelihoodDenseCholesky(pta)

enterprise.signals.signal_base.SignalCollection(metasignals)

Class factory for SignalCollection objects.

enterprise.signals.signal_base.cache_call(attrs, limit=2)

This decorator caches the output of a class method that takes a single parameter ‘params’. It saves the cache in the instance attributes _cache_<methodname> and _cache_list_<methodname>.

The cache keys are listed in the class attribute (or attributes) specified in the initial decorator call. For instance, if the decorator is applied as @cache_call(‘basis_params’), then the parameters listed in self.basis_params (together with their values) will be used as the key.

The parameter ‘limit’ specifies the number of entries saved in the cache.

enterprise.signals.signal_base.simplememobyid(method)

This decorator caches the last call of a class method that takes a single parameter arg. It holds a reference to the last arg as key, and uses the cached value if arg is key. If arg is a Sequence, then the decorator uses the cached value if the is relation is true element by element.

enterprise.signals.utils module

Utilities module containing various useful functions for use in other modules.

class enterprise.signals.utils.ConditionalGP(pta, phiinv_method='cliques')

Bases: object

This class allows the computation of conditional means and random draws for all GP coefficients/realizations in a model, given a vector of hyperparameters. It currently requires combine=False for all GPs (or otherwise distinct bases) and does not work with MarginalizingTimingModel (fast new-style likelihood).

get_mean_coefficients(params)

get_mean_processes(params)

sample_coefficients(params, n=1)

sample_processes(params, n=1)

class enterprise.signals.utils.KernelMatrix(init)

Bases: ndarray

add(other, idx)

inv(logdet=False)

set(other, idx)

enterprise.signals.utils.anis_orf(pos1, pos2, params, **kwargs)

Anisotropic GWB spatial correlation function.

enterprise.signals.utils.bwm_delay(toas, pos, log10_h=-14.0, cos_gwtheta=0.0, gwphi=0.0, gwpol=0.0, t0=55000, antenna_pattern_fn=None)

Function that calculates the Earth-term gravitational-wave Burst-With-Memory signal, as described in: Seto et al., van Haasteren and Levin, Pshirkov et al., Cordes and Jenet. This version uses the F+/Fx polarization modes, and matches the continuous-wave and anisotropy papers.

Parameters

• toas – time-of-arrival measurements [s]

• pos – unit vector from Earth to pulsar

• log10_h – log10 of GW strain

• cos_gwtheta – cosine of GW polar angle

• gwphi – GW azimuthal polar angle [rad]

• gwpol – GW polarization angle

• t0 – burst central time [day]

• antenna_pattern_fn – user defined function that takes pos, gwtheta, gwphi as arguments and returns (fplus, fcross)

Returns

the waveform as induced timing residuals (seconds)

enterprise.signals.utils.calculate_splus_scross(nmax, mc, dl, h0, F, e, t, l0, gamma, gammadot, inc)

Calculate splus and scross for a CGW summed over all harmonics. This waveform differs slightly from that in Taylor et al (2016) in that it includes the time dependence of the advance of periastron.

Parameters

• nmax – Total number of harmonics to use

• mc – Chirp mass of binary [Solar Mass]

• dl – Luminosity distance [Mpc]

• F – Orbital frequency of binary [Hz]

• e – Orbital Eccentricity

• t – TOAs [s]

• l0 – Initial eccentric anomoly [rad]

• gamma – Angle of periastron advance [rad]

• gammadot – Time derivative of angle of periastron advance [rad/s]

• inc – Inclination angle [rad]

Return splus, scross

plus and cross time-domain waveforms for a CGW

enterprise.signals.utils.create_gw_antenna_pattern(pos, gwtheta, gwphi)

Function to create pulsar antenna pattern functions as defined in Ellis, Siemens, and Creighton (2012). :param pos: Unit vector from Earth to pulsar :param gwtheta: GW polar angle in radians :param gwphi: GW azimuthal angle in radians

Returns

(fplus, fcross, cosMu), where fplus and fcross are the plus and cross antenna pattern functions and cosMu is the cosine of the angle between the pulsar and the GW source.

enterprise.signals.utils.create_quantization_matrix(toas, dt=1, nmin=2)

Create quantization matrix mapping TOAs to observing epochs.

enterprise.signals.utils.create_stabletimingdesignmatrix(designmat, fastDesign=True)

Stabilize the timing-model design matrix.

Parameters

• designmat – Pulsar timing model design matrix

• fastDesign – Stabilize the design matrix the fast way [True]

Returns

Mm: Stabilized timing model design matrix

enterprise.signals.utils.createfourierdesignmatrix_physicalephem(toas, planetssb, pos_t, frame_drift_rate=1e-09, d_jupiter_mass=1.5497669e-11, d_saturn_mass=8.17306184e-12, d_uranus_mass=5.71923361e-11, d_neptune_mass=7.96103855e-11, jup_orb_elements=0.05, sat_orb_elements=0.5, model='setIII')

Construct physical ephemeris perturbation design matrix and ‘frequencies’. Parameters can be excluded by setting the corresponding prior sigma to None

Parameters

• toas – vector of time series in seconds

• pos – pulsar position as Cartesian vector

• frame_drift_rate – normal sigma for frame drift rate

• d_jupiter_mass – normal sigma for Jupiter mass perturbation

• d_saturn_mass – normal sigma for Saturn mass perturbation

• d_uranus_mass – normal sigma for Uranus mass perturbation

• d_neptune_mass – normal sigma for Neptune mass perturbation

• jup_orb_elements – normal sigma for Jupiter orbital elem. perturb.

• sat_orb_elements – normal sigma for Saturn orbital elem. perturb.

• model – vector basis used by Jupiter and Saturn perturb.; see PhysicalEphemerisSignal, defaults to “setIII”

Returns

F: Fourier design matrix of shape (len(toas), nvecs)

Returns

sigmas: Phi sigmas (nvecs, to be passed to physicalephem_spectrum)

enterprise.signals.utils.dipole_orf(pos1, pos2)

Dipole spatial correlation function.

enterprise.signals.utils.dmass(planet, dm_over_Msun)

enterprise.signals.utils.ecl2eq_vec(x)

Rotate (n,3) vector time series from ecliptic to equatorial.

enterprise.signals.utils.eq2ecl_vec(x)

Rotate (n,3) vector time series from equatorial to ecliptic.

enterprise.signals.utils.euler_vec(z, y, x, n)

Return (n,3,3) tensor with each (3,3) block containing an Euler rotation with angles z, y, x. Optionally each of z, y, x can be a vector of length n.

enterprise.signals.utils.get_Fdot(F, mc, e)

Compute frequency derivative from Taylor et al. (2016)

Parameters

• F – Orbital frequency [Hz]

• mc – Chirp mass of binary [Solar Mass]

• e – Eccentricity of binary

Returns

dF/dt

enterprise.signals.utils.get_an(n, mc, dl, h0, F, e)

Compute a_n from Eq. 22 of Taylor et al. (2016).

Parameters

• n – Harmonic number

• mc – Chirp mass of binary [Solar Mass]

• dl – Luminosity distance [Mpc]

• F – Orbital frequency of binary [Hz]

• e – Orbital Eccentricity

Returns

a_n

enterprise.signals.utils.get_bn(n, mc, dl, h0, F, e)

Compute b_n from Eq. 22 of Taylor et al. (2015).

Parameters

• n – Harmonic number

• mc – Chirp mass of binary [Solar Mass]

• dl – Luminosity distance [Mpc]

• F – Orbital frequency of binary [Hz]

• e – Orbital Eccentricity

Returns

b_n

enterprise.signals.utils.get_cn(n, mc, dl, h0, F, e)

Compute c_n from Eq. 22 of Taylor et al. (2016).

Parameters

• n – Harmonic number

• mc – Chirp mass of binary [Solar Mass]

• dl – Luminosity distance [Mpc]

• F – Orbital frequency of binary [Hz]

• e – Orbital Eccentricity

Returns

c_n

enterprise.signals.utils.get_coefficients(pta, params, n=1, phiinv_method='cliques', variance=True, common_sparse=False)

enterprise.signals.utils.get_coupled_constecc_eqns(y, t, mc, e0)

Computes the coupled system of differential equations from Peters (1964) and Barack & Cutler (2004). This is a system of three variables:

F: Orbital frequency [Hz] phase0: Orbital phase [rad]

Parameters

• y – Vector of input parameters [F, e, gamma]

• t – Time [s]

• mc – Chirp mass of binary [Solar Mass]

Returns

array of derivatives [dF/dt, dphase/dt]

enterprise.signals.utils.get_coupled_ecc_eqns(y, t, mc, q)

Computes the coupled system of differential equations from Peters (1964) and Barack & Cutler (2004). This is a system of three variables:

F: Orbital frequency [Hz] e: Orbital eccentricity gamma: Angle of precession of periastron [rad] phase0: Orbital phase [rad]

Parameters

• y – Vector of input parameters [F, e, gamma]

• t – Time [s]

• mc – Chirp mass of binary [Solar Mass]

• q – Mass ratio of binary

Returns

array of derivatives [dF/dt, de/dt, dgamma/dt, dphase/dt]

enterprise.signals.utils.get_edot(F, mc, e)

Compute eccentricity derivative from Taylor et al. (2016)

Parameters

• F – Orbital frequency [Hz]

• mc – Chirp mass of binary [Solar Mass]

• e – Eccentricity of binary

Returns

de/dt

enterprise.signals.utils.get_gammadot(F, mc, q, e)

Compute gamma dot from Barack and Cutler (2004)

Parameters

• F – Orbital frequency [Hz]

• mc – Chirp mass of binary [Solar Mass]

• q – Mass ratio of binary

• e – Eccentricity of binary

Returns

dgamma/dt

enterprise.signals.utils.get_planet_orbital_elements(model='setIII')

Grab physical ephemeris model files

enterprise.signals.utils.hd_orf(pos1, pos2)

Hellings & Downs spatial correlation function.

enterprise.signals.utils.linear_interp_basis(toas, dt=2592000)

Provides a basis for linear interpolation.

Parameters

• toas – Pulsar TOAs in seconds

• dt – Linear interpolation step size in seconds.

Returns

Linear interpolation basis and nodes

enterprise.signals.utils.make_ecc_interpolant()

Make interpolation function from eccentricity file to determine number of harmonics to use for a given eccentricity.

Returns

interpolant

enterprise.signals.utils.monopole_orf(pos1, pos2)

Monopole spatial correlation function.

enterprise.signals.utils.normed_tm_basis(Mmat, norm=None)

enterprise.signals.utils.physical_ephem_delay(toas, planetssb, pos_t, frame_drift_rate=0, d_jupiter_mass=0, d_saturn_mass=0, d_uranus_mass=0, d_neptune_mass=0, jup_orb_elements=array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), sat_orb_elements=array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0]), times=None, jup_orbit=None, sat_orbit=None, equatorial=True)

enterprise.signals.utils.physicalephem_spectrum(sigmas)

enterprise.signals.utils.quant2ind(U)

Use quantization matrix to return slices of non-zero elements.

Parameters

U – quantization matrix

Returns

list of `slice`s for non-zero elements of U

Note

This function assumes that the pulsar TOAs were sorted by time.

enterprise.signals.utils.solve_coupled_constecc_solution(F0, e0, phase0, mc, t)

Compute the solution to the coupled system of equations from from Peters (1964) and Barack & Cutler (2004) at a given time.

Parameters

• F0 – Initial orbital frequency [Hz]

• mc – Chirp mass of binary [Solar Mass]

• t – Time at which to evaluate solution [s]

Returns

(F(t), phase(t))

enterprise.signals.utils.solve_coupled_ecc_solution(F0, e0, gamma0, phase0, mc, q, t)

Compute the solution to the coupled system of equations from from Peters (1964) and Barack & Cutler (2004) at a given time.

Parameters

• F0 – Initial orbital frequency [Hz]

• e0 – Initial orbital eccentricity

• gamma0 – Initial angle of precession of periastron [rad]

• mc – Chirp mass of binary [Solar Mass]

• q – Mass ratio of binary

• t – Time at which to evaluate solution [s]

Returns

(F(t), e(t), gamma(t), phase(t))

enterprise.signals.utils.ss_framerotate(mjd, planet, x, y, z, dz, offset=None, equatorial=False)

Rotate planet trajectory given as (n,3) tensor, by ecliptic Euler angles x, y, z, and by z rate dz. The rate has units of rad/year, and is referred to offset 2010/1/1. dates must be given in MJD.

enterprise.signals.utils.svd_tm_basis(Mmat)

enterprise.signals.utils.tm_prior(weights)

enterprise.signals.utils.unnormed_tm_basis(Mmat)

enterprise.signals.white_signals module

Contains class factories for white noise signals. White noise signals are defined as the class of signals that only modifies the white noise matrix N.

enterprise.signals.white_signals.EcorrKernelNoise(log10_ecorr=<class 'enterprise.signals.parameter.Uniform.<locals>.Uniform'>, selection=<class 'enterprise.signals.selections.Selection.<locals>.Selection'>, method='sherman-morrison', name='')

Class factory for ECORR type noise.

Parameters

• log10_ecorr – Parameter type for log10 or ecorr parameter.

• selection – Selection object specifying masks for backends, time segments, etc.

• method – Method for computing noise covariance matrix. Options include sherman-morrison, sparse, and block

Returns

EcorrKernelNoise class.

ECORR is a noise signal that is used for data with multi-channel TOAs that are nearly simultaneous in time. It is a white noise signal that is uncorrelated epoch to epoch but completely correlated for TOAs in a given observing epoch.

For this implementation we use this covariance matrix as part of the white noise covariance matrix \(N\). It can be seen from above that this covariance is block diagonal, thus allowing us to exploit special methods to make matrix manipulations easier.

In this signal implementation we offer three methods of performing these matrix operations:

sherman-morrison

Uses the Sherman-Morrison forumla to compute the matrix inverse and other matrix operations. Note: This method can only be used for covariances that can be constructed by the outer product of two vectors, \(uv^T\).

sparse

Uses Scipy Sparse matrices to construct the block diagonal covariance matrix and perform matrix operations.

block

Uses a custom scheme that uses the individual blocks from the block diagonal matrix to perform fast matrix inverse and other solve operations.

Note

The sherman-morrison method is the fastest, followed by the block and then sparse methods, however; the block and sparse methods are more general and should be used if sub-classing this signal for more complicated blocks.

enterprise.signals.white_signals.EquadNoise(*args, **kwargs)

enterprise.signals.white_signals.MeasurementNoise(efac=<class 'enterprise.signals.parameter.Uniform.<locals>.Uniform'>, log10_t2equad=None, selection=<class 'enterprise.signals.selections.Selection.<locals>.Selection'>, name='')

Class factory for EFAC+EQUAD measurement noise (with tempo/tempo2/pint parameter convention, variance = efac^2 (toaerr^2 + t2equad^2)). Leave out log10_t2equad to use EFAC noise only.

enterprise.signals.white_signals.TNEquadNoise(log10_tnequad=<class 'enterprise.signals.parameter.Uniform.<locals>.Uniform'>, selection=<class 'enterprise.signals.selections.Selection.<locals>.Selection'>, name='')

Class factory for TNEQUAD type measurement noise (legacy, not multiplied by EFAC).

enterprise.signals.white_signals.WhiteNoise(varianceFunction, selection=<class 'enterprise.signals.selections.Selection.<locals>.Selection'>, name='')

Class factory for generic white noise signals.

enterprise.signals.white_signals.combined_ndiag(toaerrs, efac=1.0, log10_t2equad=-8)

enterprise.signals.white_signals.efac_ndiag(toaerrs, efac=1.0)

enterprise.signals.white_signals.tnequad_ndiag(toas, log10_tnequad=-8)

Submodules

enterprise.constants module

Declares physical constants for use in enterprise. Depends on numpy for base mathematical constants, and scipy.constants for physical constants.

enterprise.pulsar module

Class containing pulsar data from timing package [tempo2/PINT].

class enterprise.pulsar.BasePulsar

Bases: object

Abstract Base Class for Pulsar objects.

filter_data(start_time=None, end_time=None)

Filter data to create a time-slice of overall dataset.

set_flags(flagname, values)

Set value of existing or new flags.

sort_data()

Sort data by time.

to_pickle(outdir=None)

Save object to pickle file.

property Mmat

Return ntoa x npar design matrix.

property backend_flags

Return array of backend flags.

Not all TOAs have the same flags for all data sets. In order to facilitate this we have a ranked ordering system that will look for flags. The order is group, g, sys, i, f, fe`+`be.

property dm

Return DM parameter from parfile.

property dmx

Return a dictionary of DMX-parameter values and stoa ranges from parfile.

property flags

Return a dictionary of tim-file flags.

property freqs

Return array of radio frequencies in MHz.

property iisort

Return inverse of sorting indices.

property isort

Return sorting indices.

property pdist

Return tuple of pulsar distance and uncertainty in kpc.

property phi

Return azimuthal angle of pulsar in radians.

property planetssb

Return planetary position vectors at all timestamps

property pos

Return unit vector from SSB to pulsar at fiducial POSEPOCH.

property pos_t

Return unit vector from SSB to pulsar as function of time.

property residuals

Return array of residuals in seconds.

property stoas

Return array of observatory TOAs in seconds.

property sunssb

Return sun position vector at all timestamps

property telescope

Return telescope vector at all timestamps

property theta

Return polar angle of pulsar in radians.

property toaerrs

Return array of TOA errors in seconds.

property toas

Return array of TOAs in seconds.

class enterprise.pulsar.PintPulsar(toas, model, sort=True, drop_pintpsr=True, planets=True)

Bases: BasePulsar

class enterprise.pulsar.Tempo2Pulsar(t2pulsar, sort=True, drop_t2pulsar=True, planets=True)

Bases: BasePulsar

deflate()

destroy()

inflate()

enterprise.pulsar.Pulsar(*args, **kwargs)

enterprise.pulsar.get_maxobs(timfile)

Utility function to return number of lines in tim file. :param timfile:

Full path to tim-file. For tim-files that use INCLUDEs this should be the base tim file.

Returns

Number of lines in tim-file

enterprise.pulsar_inflate module

Defines PulsarInflater class: instances copy a numpy array to shared memory, and (after pickling) will reinflate to a numpy array that refers to the shared data.

class enterprise.pulsar_inflate.PulsarInflater(array)

Bases: object

destroy()

inflate()

class enterprise.pulsar_inflate.memmap

Bases: ndarray

tests package

Submodules

tests.enterprise_test_data module

This module provides the location of data files for tests as datadir. Currently they are in tests/data, and they are based on the 9-yr data release.

tests.test_deterministic_signals module

Deterministic Signals are those that provide a get_delay() method, and are created by the class factories in enterprise.signals.deterministic_signals. All tests in this module are run on B1855+09_NANOGrav_9yv1.

class tests.test_deterministic_signals.TestDeterministicSignals(methodName='runTest')

Bases: TestCase

Tests deterministic signals with a tempo2 Pulsar object.

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Set up the enterprise.Pulsar() object used in tests (tempo2 version).

test_bwm()

Tests enterprise.signals.deterministic_signals.Deterministic() using the burst-with-memory function enterprise.signals.utils.bwm_delay(). The test instantiates a deterministic Signal on our test pulsar, and compares the array returned by calling get_delay() on the Signal with a fixed dictionary of parameters, with the result of calling the function directly with those parameters.

test_delay()

Same as TestDeterministicSignals.test_bwm(), but for a simple sine wave signal.

test_delay_backend()

Same as TestDeterministicSignals.test_delay(), but instantiates the Signal with enterprise.signals.selections.by_backend(), which creates separated named parameters for 430_ASP, 430_PUPPI, L-wide_ASP, L-wide_PUPPI. The parameters are automatically accounted for in get_delay(), but they need to be used explicitly when calling the function directly. The tests therefore reconstructs the delay vector by building selection masks from enterprise.Pulsar.backend_flags().

test_physical_ephem_model()

Tests physical ephemeris model (which is implemented as a deterministic signal) four ways:

• computed directly with enterprise.signals.utils.physical_ephem_delay();

• computed with enterprise.signals.deterministic_signals.PhysicalEphemerisSignal.get_delay() with use_epoch_toas=True (the default), which reduces computation by evaluating ephemeris corrections once per measurement epoch, and then interpolating to the full toas vector;

• computed with enterprise.signals.deterministic_signals.PhysicalEphemerisSignal.get_delay(), setting use_epoch_toas=False;

• loaded from a golden copy.

test_physical_ephem_model_setIII()

Test physical ephemeris model

class tests.test_deterministic_signals.TestDeterministicSignalsPint(methodName='runTest')

Bases: TestDeterministicSignals

Tests deterministic signals with a PINT Pulsar object.

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Set up the enterprise.Pulsar() object used in tests (PINT version).

tests.test_deterministic_signals.sine_wave(toas, log10_A=-7, log10_f=-8, phase=0.0)

A simple sine wave Enterprise function object. When instantiated, it will create named Parameters for log10_A, log10_f, phase, and it will automatically extract toas from the linked Pulsar object.

tests.test_gp_coefficients module

test_gp_coefficients

Tests for GP signals used with deterministic coefficients.

class tests.test_gp_coefficients.TestGPCoefficients(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_common_red_noise()

Test of a coefficient-based common GP.

test_conditional_gp()

test_ecorr_backend()

Test that ecorr-backend signal returns correct values.

test_ephemeris()

Test physical-ephemeris delay, made three ways: from marginalized GP, from coefficient-based GP, from deterministic model.

test_formalism()

test_fourier_red_noise()

Test that implicit and explicit GP delays are the same.

class tests.test_gp_coefficients.TestGPCoefficientsPint(methodName='runTest')

Bases: TestGPCoefficients

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_ephemeris()

Test physical-ephemeris delay, made three ways: from marginalized GP, from coefficient-based GP, from deterministic model.

tests.test_gp_coefficients.create_quant_matrix(toas, dt=1)

tests.test_gp_coefficients.se_kernel(etoas, log10_sigma=-7, log10_lam=6.413634997198556)

tests.test_gp_priors module

test_gp_priors

Tests for GP priors and bases.

class tests.test_gp_priors.TestGPSignals(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_adapt_t_process_prior()

Test that red noise signal returns correct values.

test_broken_powerlaw_prior()

Test that red noise signal returns correct values.

test_free_spec_prior()

Test that red noise signal returns correct values.

test_powerlaw_genmodes_prior()

Test that red noise signal returns correct values.

test_t_process_prior()

Test that red noise signal returns correct values.

test_turnover_knee_prior()

Test that red noise signal returns correct values.

test_turnover_prior()

Test that red noise signal returns correct values.

tests.test_gp_signals module

test_gp_signals

Tests for GP signal modules.

class tests.test_gp_signals.TestGPSignals(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_combine_signals()

Test for combining different signals.

test_ecorr()

Test that ecorr signal returns correct values.

test_ecorr_backend()

Test that ecorr-backend signal returns correct values.

test_fourier_red_noise()

Test that red noise signal returns correct values.

test_fourier_red_noise_backend()

Test that red noise-backend signal returns correct values.

test_fourier_red_noise_pshift()

Test that red noise signal returns correct values.

test_fourier_red_user_freq_array()

Test that red noise signal returns correct values with user defined frequency array.

test_gp_parameter()

Test GP basis model with parameterized basis.

test_gp_timing_model()

Test that the timing model signal returns correct values.

test_kernel()

test_kernel_backend()

test_pshift_fourier()

Test Fourier basis with prescribed phase shifts.

test_red_noise_add()

Test that red noise addition only returns independent columns.

test_red_noise_add_backend()

Test that red noise with backend addition only returns independent columns.

class tests.test_gp_signals.TestGPSignalsPint(methodName='runTest')

Bases: TestGPSignals

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

tests.test_gp_signals.create_quant_matrix(toas, dt=1)

tests.test_gp_signals.se_kernel(etoas, log10_sigma=-7, log10_lam=6.413634997198556)

tests.test_gp_wideband module

test_gp_wideband

Tests for WidebandTimingModel.

class tests.test_gp_wideband.TestGPSignalsPint(methodName='runTest')

Bases: TestWidebandTimingModel

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_wideband()

class tests.test_gp_wideband.TestWidebandTimingModel(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_wideband()

tests.test_hierarchical_parameter module

test_hierarchical_parameter

Tests for hierarchical parameter functionality

class tests.test_hierarchical_parameter.TestHierarchicalParameter(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_Function_of_Function()

test_enterprise_Function()

test_enterprise_Parameter()

test_powerlaw()

test_powerlaw_equad()

tests.test_likelihood module

test_likelihood

Tests of likelihood module

class tests.test_likelihood.TestLikelihood(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

compute_like(npsrs=1, inc_corr=False, inc_kernel=False, cholesky_sparse=True, marginalizing_tm=False)

classmethod setUpClass()

Setup the Pulsar object.

test_compare_ecorr_likelihood()

Compare basis and kernel ecorr methods.

test_like_corr()

Test likelihood with spatial correlations.

test_like_corr_kernel()

Test likelihood with spatial correlations and kernel.

test_like_corr_one_kernel()

Test likelihood with spatial correlations and one kernel.

test_like_nocorr()

Test likelihood with no spatial correlations.

test_like_nocorr_kernel()

Test likelihood with no spatial correlations and kernel.

test_like_nocorr_one_kernel()

Test likelihood with no spatial correlations and one kernel.

class tests.test_likelihood.TestLikelihoodPint(methodName='runTest')

Bases: TestLikelihood

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

tests.test_likelihood.create_quant_matrix(toas, dt=1)

tests.test_likelihood.get_noise_from_pal2(noisefile)

tests.test_likelihood.se_kernel(etoas, log10_sigma=-7, log10_lam=6.413634997198556)

tests.test_parameter module

test_parameter

Tests Uniform and Normal parameter priors and sampling functions

class tests.test_parameter.TestParameter(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

test_linearexp()

Test LinearExp parameter prior and sampler.

test_metaparam()

test_normal()

Test Normal parameter prior and sampler for various combinations of scalar and vector arguments.

test_normalandtruncnormal()

test_truncnormal()

Test TruncNormal parameter prior and sampler for various combinations of scalar and vector arguments.

test_uniform()

Test Uniform parameter prior and sampler for various combinations of scalar and vector arguments.

tests.test_pta module

test_pta

Tests for common signal and PTA class modules.

class tests.test_pta.TestPTASignals(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_parameterized_orf()

test_pta_phi()

test_pta_phiinv_methods()

test_summary()

Test PTA summary table as well as its str representation and dict-like interface.

class tests.test_pta.TestPTASignalsPint(methodName='runTest')

Bases: TestPTASignals

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

tests.test_pta.hd_orf_generic(pos1, pos2, a=1.5, b=0.25, c=0.25)

tests.test_pta.hd_powerlaw(f, pos1, pos2, log10_A=-15, gamma=4.3)

tests.test_pulsar module

tests.test_selections module

test_pulsar

Tests for signals/selections module.

class tests.test_selections.TestSelections(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_selections()

class tests.test_selections.TestSelectionsPint(methodName='runTest')

Bases: TestSelections

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

tests.test_set_parameter module

test_set_parameter

Tests of setting constant parameters

class tests.test_set_parameter.TestSetParameters(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_pta()

test_single_pulsar()

class tests.test_set_parameter.TestSetParametersPint(methodName='runTest')

Bases: TestSetParameters

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

tests.test_set_parameter.get_noise_from_pal2(noisefile)

tests.test_utils module

test_utils

Tests for utils module.

class tests.test_utils.TestUtils(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_createfourierdesignmatrix_dm(nf=30)

Check DM-variation Fourier design matrix shape.

test_createfourierdesignmatrix_ephem(nf=30)

Check x-axis ephemeris Fourier design matrix shape.

test_createfourierdesignmatrix_red(nf=30)

Check Fourier design matrix shape.

test_createstabletimingdesignmatrix()

Timing model design matrix shape.

test_ecc_cw_waveform()

Check eccentric wafeform generation.

test_fplus_fcross()

Check fplus, fcross generation.

test_numerical_ecc_integration()

Test numerical integration of eccentric GW.

test_orf()

Test ORF functions.

test_psd()

Test PSD functions.

test_quantization_matrix()

Test quantization matrix generation.

tests.test_vector_parameter module

test_vector_parameter

Tests for vector parameter functionality

class tests.test_vector_parameter.TestVectorParameter(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_phi()

Test vector parameter on signal level.

test_vector_parameter_like()

Test vector parameter in a likelihood

class tests.test_vector_parameter.TestVectorParameterPint(methodName='runTest')

Bases: TestVectorParameter

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

tests.test_vector_parameter.free_spectrum(f, log10_rho=None)

tests.test_white_signals module

test_white_signals

Tests for white signal modules.

class tests.test_white_signals.TestWhiteSignals(methodName='runTest')

Bases: TestCase

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

test_add_efac_tnequad()

Test that addition of efac and tnequad signal returns correct covariance.

test_add_efac_tnequad_backend()

Test that addition of efac-backend and tnequad-backend signal returns correct covariance.

test_ecorr_block()

Test of block matrix ecorr signal and solve methods.

test_ecorr_block_ipta()

Test of block matrix ecorr signal and solve methods.

test_ecorr_sherman_morrison()

Test of sherman-morrison ecorr signal and solve methods.

test_ecorr_sherman_morrison_ipta()

Test of sherman-morrison ecorr signal and solve methods.

test_ecorr_sparse()

Test of sparse ecorr signal and solve methods.

test_ecorr_sparse_ipta()

Test of sparse ecorr signal and solve methods.

test_efac()

Test that efac signal returns correct covariance.

test_efac_backend()

Test that backend-efac signal returns correct covariance.

test_efac_equad_combined_backend()

Test that the combined EFAC + EQUAD noise (tempo/tempo2/pint definition) returns the correct covariance.

test_equad()

Test that the deprecated EquadNoise is not available.

test_tnequad()

Test that tnequad signal returns correct covariance.

test_tnequad_backend()

Test that backend-equad signal returns correct covariance.

class tests.test_white_signals.TestWhiteSignalsPint(methodName='runTest')

Bases: TestWhiteSignals

Create an instance of the class that will use the named test method when executed. Raises a ValueError if the instance does not have a method with the specified name.

classmethod setUpClass()

Setup the Pulsar object.

class tests.test_white_signals.Woodbury(N, U, J)

Bases: object

logdet()

solve(other)

Development Leads

• Justin A. Ellis (@jellis18) <https://github.com/jellis18>

• Stephen R. Taylor (@stevertaylor) <https://github.com/stevertaylor>

• Paul T. Baker (@paulthebaker) <https://github.com/paulthebaker>

• Michele Vallisneri (@vallis) <https://github.com/vallis>