Add ability to simulate noise in Pint TOAs

[Hazboun6]

A few of us developed code that is basically the Pint rewrite of the noise simulation pieces of toasim.py in libstempo. This was developed by myself, @andrealommen and some of her students for adding noise to NICER X-ray TOAs, but should be generally applicable to any TOAs.

New methods include:

make_ideal
add_efac
add_equad
add_ecorr
add_dm_rednoise
add_rednoise

I think the names are fairly understandable to those who are familiar with these common types of pulsar timing data noise. Both red noise methods add in time-correlated noise constructed using a Fourier basis. The add_dm_rednoise method adds in that noise with a dependence on the observation frequency of 1/RF^2, but that could easily be generalized to any 1/RF^p dependence.

This does not include any gravitational wave simulation code, but that could definitely be added in the future.

I'll post the code in a separate comment below. There is a chance the code is a bit dated at this point, but I'm happy to update and add this into Pint.

Two questions:

1. Given that I'll write some tests and update the code, is this something folks would like to see a part of the package?
2. Where is the best place for this? I could follow the libstempo naming and call it toa_sim, or start a noise.py. There is the toa.make_fake_toas method, but it seems like there is plenty of code in there already.

[Hazboun6]

The current state of the code is here.

[Hazboun6]

@aarchiba is this still something that would be useful to you?

[aarchiba]

This does seem useful! I had envisioned PINT's noise model objects being able to generate noise of the kind they describe - so if you had a noise model describing red noise of a particular amplitude and spectral index, or an ECORR affecting a certain category of TOAs, you could just pass those models TOAs and they would generate adjustments to reflect an instance of the noise they describe.

It seems to me that the machinery here is more or less capable of doing that - that is, that given an appropriate description of the noise model they can generate the noise. So it might be just a matter of translating from the noise model parameters in the PINT objects to the appropriate values for these functions?

[Hazboun6]

This makes sense to me. One could add a method to the Base class, something like noise_model. NoiseComponent.add_fake_noise(), and populate that with the various methods from above.

Should the code keep track of whether one has added in fake signals of some sort? Once we start adding methods to classes that normally hold real data with real noise, I'd worry about accidental misuse. Although maybe I'm trying to do too much hand holding for the user.

Or a slightly different option would be to use an underscored method, noise_model. NoiseComponent._add_fake_noise() that rides along with the noise model and gets called when the user calls a separate method, noise_model.add_fake_noise(toas, model). That method would then look for all of the noise models in model and add them to the appropriate TOAs using the _add_fake_noise() method.

[aarchiba]

Adjusting TOAs is sort of expensive (stuff ends up getting recalculated), so I'd rather just return an array with a noise adjustment per TOA. And TimingModel objects don't actually contain any TOAs, or even have any rigid association with any. So the pattern so far of TimingModel functions that act on TOAs is that they get an argument toas, and they either modify that or return a relevant derived thing (say a list of the indices that match a selection criterion). So I'd go with model.add_fake_noise(toas) that returns an array in units of (micro)seconds. The TimingModel version of the function would call noise_component.add_fake_noise(toas) for each noise_component and return the sum.

The main drawback of this approach is that it's (currently) a bit of a pain in the neck to build a TimingModel or a collection of NoiseComponents if all you need them for is to feed in a set of noise model parameters. But I think that's probably fixable. Right now it looks something like:

n = EcorrNoise()
n.ECORR1.flag = "-fe"
n.ECORR1.flag_value = "430"
n.ECORR1.quantity = 1*u.us
n.generate_noise(toas)

(and worse if you want more ECORRs) but we need to make adding parameters and making components less painful.

[abhisrkckl]

Already implemented.