-
Notifications
You must be signed in to change notification settings - Fork 101
How To
Quick solutions to common tasks (taken from #pint and elsewhere). It's also worth checking the PINT HOWTOs, particularly as some things that used to be here were moved there.
With pip:
pip install -U pint-pulsar
With conda:
conda update pint-pulsar
With pip:
pip install -U pint-pulsar==0.8.4
or similar.
With conda:
conda install pint-pulsar=0.8.4
or similar.
The data files (par and tim) associated with the tutorials and other examples can be located via pint.config.examplefile() (available via the pint.config module):
import pint.config
fullfilename = pint.config.examplefile(filename)For example, the file NGC6440E.par from the Time a Pulsar notebook can be found via:
import pint
fullfilename = pint.config.examplefile("NGC6440E.par")from pint.models import get_model
m = get_model(parfile)from pint.toa import get_TOAs
t = get_TOAs(timfile)from pint.models import get_model_and_toas
m, t = get_model_and_toas(parfile, timfile)You can index by column name into the TOAs object, so you can do toas["observatory"] or whatever the column is called; and that's an array, so you can do toas["observatory"]=="arecibo" to get a Boolean array; and you can index with boolean arrays, so you can do toas[toas["observatory"]=="arecibo"] to get a new TOAs object referencing a subset.
The TOAs have a table with mjd, mjd_float, tdb, and tdbld columns. To modify them all safely and consistently the best way is to use:
t.adjust_TOAs(dt)where dt is an astropy.time.TimeDelta object. This function does not change the pulse numbers column, if present, but does recompute mjd_float, the TDB times, and the observatory positions and velocities.
You need to have the TOAs object compute the positions of the planets and add them to the table.
ts.compute_posvels(ephem,planets=True)This should be done automatically if you load your TOAs with the get_TOAs() or get_model_and_TOAs()
model.as_ICRS(epoch=epoch)which will give it to you as a model with AstrometryEquatorial components at the requested epoch. Similarly,
model.as_ECL(epoch=epoch)does the same for AstrometryEcliptic (with an optional specification of the obliquity).
PINT can handle jumps in the model outside a par file. An example is:
import numpy as np
from astropy import units as u, constants as c
from pint.models import get_model, get_model_and_toas, parameter
from pint import fitter
from pint.models import PhaseJump
import pint.config
m, t = get_model_and_toas(pint.config.examplefile("NGC6440E.par"),
pint.config.examplefile("NGC6440E.tim"))
# fit the nominal model
f = fitter.WLSFitter(toas=t, model=m)
f.fit_toas()
# group TOAs: find clusters with gaps of <2h
clusters = t.get_clusters(add_column=True)
# put in the pulse numbers based on the previous fit
t.compute_pulse_numbers(f.model)
# just for a test, add an offset to a set of TOAs
t['delta_pulse_number'][clusters==3]+=3
# now fit without a jump
fnojump = fitter.WLSFitter(toas=t, model=m, track_mode="use_pulse_numbers")
fnojump.fit_toas()
# add the Jump Component to the model
m.add_component(PhaseJump(), validate=False)
# now add the actual jump
# it can be keyed on any parameter that maskParameter will accept
# here we will use a range of MJDs
par = parameter.maskParameter(
"JUMP",
key="mjd",
value=0.0,
key_value=[t[clusters==3].get_mjds().min().value,
t[clusters==3].get_mjds().max().value],
units=u.s,
frozen=False,
)
m.components['PhaseJump'].add_param(par, setup=True)
# you can also do it indirectly through the flags as:
# m.components["PhaseJump"].add_jump_and_flags(t.table["flags"][clusters == 3])
# and fit with a jump
fjump = fitter.WLSFitter(toas=t, model=m, track_mode="use_pulse_numbers")
fjump.fit_toas()
print(f"Original chi^2 = {f.resids.calc_chi2():.2f} for {f.resids.dof} DOF")
print(f"After adding 3 rotations to some TOAs, chi^2 = {fnojump.resids.calc_chi2():.2f} for {fnojump.resids.dof} DOF")
print(f"Then after adding a jump to those TOAs, chi^2 = {fjump.resids.calc_chi2():.2f} for {fjump.resids.dof} DOF")
print(f"Best-fit value of the jump is {fjump.model.JUMP1.quantity} +/- {fjump.model.JUMP1.uncertainty} ({(fjump.model.JUMP1.quantity*fjump.model.F0.quantity).decompose():.3f} +/- {(fjump.model.JUMP1.uncertainty*fjump.model.F0.quantity).decompose():.3f} rotations)")which returns:
Original chi^2 = 59.57 for 56 DOF
After adding 3 rotations to some TOAs, chi^2 = 19136746.30 for 56 DOF
Then after adding a jump to those TOAs, chi^2 = 56.60 for 55 DOF
Best-fit value of the jump is -0.048772786677935796 s +/- 1.114921182802775e-05 s (-2.999 +/- 0.001 rotations)
showing that the offset we applied has been absorbed by the jump (plus a little extra, so chi^2 has actually improved).
See maskParamter documentation on the ways to select the TOAs.
See Fitter documentation, but to summarize:
| Wideband | Narrowband | ||
|---|---|---|---|
| Correlated errors (ECORR or RN) | No correlated errors | ||
| Downhill | WidebandDownhillFitter | DownhillGLSFitter | DownhillWLSFitter |
| Not downhill | WidebandTOAFitter | GLSFitter | WLSFitter |
PINT now uses loguru for its logging. To get this working within a script, try:
import pint.logging
from loguru import logger as log
pint.logging.setup(sink=sys.stderr, level="WARNING", usecolors=True)That sets up the logging and ensures it will play nicely with the rest of PINT. You can customize the level, the destination (e.g., file, stderr, ...) and format. The LogFilter suppresses some INFO/DEBUG messages that can clog up your screen: you can make a custom filter as well to add/remove messages.
If you want to include a standard way to control the level using command line arguments, you can do:
parser.add_argument(
"--log-level",
type=str,
choices=("TRACE", "DEBUG", "INFO", "WARNING", "ERROR"),
default=pint.logging.script_level,
help="Logging level",
dest="loglevel",
)
...
pint.logging.setup(level=args.loglevel, ...)assuming you are using argparse. Note that loguru doesn't let you change existing loggers: you should just remove and add (which the setup() function does).