This Jupyter notebook can be downloaded from build_model_from_scratch.ipynb, or viewed as a python script at build_model_from_scratch.py.
Building a timing model from scratch
This example includes: * Constructing a timing model object from scratch * Adding and deleting components * Assigning parameter values * Adding prefix-able parameters
[1]:
import astropy.units as u # Astropy units is a very useful module.
import pint.logging
try:
from IPython.display import display
except ImportError:
# Older IPython
from IPython.core.display_functions import display
# setup logging
pint.logging.setup(level="INFO")
from pint.models import (
parameter as p,
) # We would like to add parameters to the model, so we need parameter module.
from pint.models.timing_model import (
TimingModel,
Component,
) # Interface for timing model
import pint
from astropy.time import Time # PINT uses astropy Time objects to represent times
Typically, timing models are built by reading a par file with the get_model() function, but it is possible to construct them entirely programmatically from scratch. Also, once you have a TimingModel object (no matter how you built it), you can modify it by adding or removing parameters or entire components. This example show how this is done.
We are going to build the model for “NGC6440E.par” from scratch
First let us see all the possible components we can use
All built-in component classes can be viewed from Component class, which uses the meta-class to collect the built-in component class. For how to make a component class, see example “make_component_class” (in preparation)
[2]:
# list all the existing components
# all_components is a dictionary, with the component name as the key and component class as the value.
all_components = Component.component_types
# Print the component class names.
_ = [print(x) for x in all_components] # The "_ =" just suppresses excess output
AbsPhase
AstrometryEquatorial
AstrometryEcliptic
BinaryBT
BinaryBTPiecewise
BinaryDD
BinaryDDS
BinaryDDGR
BinaryDDH
BinaryDDK
BinaryELL1
BinaryELL1H
BinaryELL1k
DispersionDM
DispersionDMX
DispersionJump
DMWaveX
FD
Glitch
PhaseOffset
PiecewiseSpindown
IFunc
PhaseJump
ScaleToaError
ScaleDmError
EcorrNoise
PLDMNoise
PLRedNoise
SolarSystemShapiro
SolarWindDispersion
SolarWindDispersionX
Spindown
FDJump
TroposphereDelay
Wave
WaveX
Choose your components
Let’s start from a relatively simple model, with AbsPhase: The absolute phase of the pulsar, typical parameters, TZRMJD, TZRFREQ… AstrometryEquatorial: The ICRS equatorial coordinate, parameters, RAJ, DECJ, PMRA, PMDEC… Spindown: The pulsar spin-down model, parameters, F0, F1…
We will add a dispersion model as a demo.
[3]:
selected_components = ["AbsPhase", "AstrometryEquatorial", "Spindown"]
component_instances = []
# Initiate the component instances
for cp_name in selected_components:
component_class = all_components[cp_name] # Get the component class
component_instance = component_class() # Instantiate a component object
component_instances.append(component_instance)
Construct timing model (i.e., TimingModel instance)
TimingModel class provides the storage and interface for the components. It also manages the components internally.
[4]:
# Construct timing model instance, given a name and a list of components to include (that we just created above)
tm = TimingModel("NGC6400E", component_instances)
View the components in the timing model instance
To view all the components in TimingModel instance, we can use the property .components, which returns a dictionary (name as the key, component instance as the value).
Internally, the components are stored in a list(ordered list, you will see why this is important below) according to their types. All the delay type of components (subclasses of DelayComponent class) are stored in the DelayComponent_list, and the phase type of components (subclasses of PhaseComponent class) in the PhaseComponent_list.
[5]:
# print the components in the timing model
for cp_name, cp_instance in tm.components.items():
print(cp_name, cp_instance)
AbsPhase AbsPhase(
MJDParameter( TZRMJD UNSET,
strParameter( TZRSITE UNSET,
floatParameter( TZRFRQ UNSET)
Spindown Spindown(
floatParameter( F0 0.0 (Hz) frozen=True),
MJDParameter( PEPOCH UNSET)
AstrometryEquatorial AstrometryEquatorial(
MJDParameter( POSEPOCH UNSET,
floatParameter( PX 0.0 (mas) frozen=True),
AngleParameter( RAJ UNSET,
AngleParameter( DECJ UNSET,
floatParameter( PMRA 0.0 (mas / yr) frozen=True),
floatParameter( PMDEC 0.0 (mas / yr) frozen=True))
Useful methods of TimingModel
TimingModel.components(): List all the components in the timing model.TimingModel.add_component(): Add component into the timing model.TimingModel.remove_component(): Remove a component from the timing model.TimingModel.params(): List all the parameters in the timing model from all components.TimingModel.setup(): Setup the components (e.g., register the derivatives).TimingModel.validate(): Validate the components see if the parameters are setup properly.TimingModel.delay(): Compute the total delay.TimingModel.phase(): Compute the total phase.TimingModel.delay_funcs(): List all the delay functions from all the delay components.TimingModel.phase_funcs(): List all the phase functions from all the phase components.TimingModel.get_component_type(): Get all the components from one categoryTimingModel.map_component(): Get the component location. It returns the component’s instance, order in the list, host list and its type.TimingModel.get_params_mapping(): Report which component each parameter comes from.TimingModel.get_prefix_mapping(): Get the index mapping for one prefix parameters.TimingModel.param_help(): Print the help line for all available parameters.
Component order
Since the times that are given to a delay component include all the delays from the previously-evaluted delay components, the order of delay components is important. For example, the solar system delays need to be applied to get to barycentric time, which is needed to evaluate the binary delays, then the binary delays must be applied to get to pulsar proper time.
PINT provides a default ordering for the components. In most cases this should be correct, but can be modified by expert users for a particular purpose.
Here is the default order:
[6]:
from pint.models.timing_model import DEFAULT_ORDER
_ = [print(order) for order in DEFAULT_ORDER]
astrometry
jump_delay
troposphere
solar_system_shapiro
solar_wind
dispersion_constant
dispersion_dmx
dispersion_jump
pulsar_system
frequency_dependent
absolute_phase
spindown
phase_jump
wave
wavex
Add parameter values
Initially, the parameters have no values or the default values, so we must add them before validating the model.
Please note, PINT’s convention for fitting flag is defined in the Parameter.frozen attribute. Parameter.frozen = True means “do not fit this parameter”. This is the opposite of TEMPO/TEMPO2 .par file flag where “1” means the parameter is fitted.
[7]:
# We build a dictionary with a key for each parameter we want to set.
# The dictionary entries can be either
# {'pulsar name': (parameter value, TEMPO_Fit_flag, uncertainty)} akin to a TEMPO par file form
# or
# {'pulsar name': (parameter value, )} for parameters that can't be fit
# NOTE: The values here are assumed to be in the default units for each parameter
# Notice that we assign values with units, and pint defines a special hourangle_second unit that can be use for
# right ascensions. Also, angles can be specified as strings that will be parsed by astropy.
params = {
"PSR": ("1748-2021E",),
"RAJ": ("17:48:52.75", 1, 0.05 * pint.hourangle_second),
"DECJ": ("-20:21:29.0", 1, 0.4 * u.arcsec),
"F0": (61.485476554 * u.Hz, 1, 5e-10 * u.Hz),
"PEPOCH": (Time(53750.000000, format="mjd", scale="tdb"),),
"POSEPOCH": (Time(53750.000000, format="mjd", scale="tdb"),),
"TZRMJD": (Time(53801.38605120074849, format="mjd", scale="tdb"),),
"TZRFRQ": (1949.609 * u.MHz,),
"TZRSITE": (1,),
}
# Assign the parameters
for name, info in params.items():
par = getattr(tm, name) # Get parameter object from name
par.quantity = info[0] # set parameter value
if len(info) > 1:
if info[1] == 1:
par.frozen = False # Frozen means not fit.
par.uncertainty = info[2]
Set up and Validating the model
Setting up the model builds the necessary model attributes, and validating model checks if there is any important parameter values missing, and if the parameters are assigned correctly. If there is anything not assigned correctly, it will raise an exception.
[8]:
tm.setup()
tm.validate()
# You should see all the assigned parameters.
# Printing a TimingModel object shows the parfile representation
print(tm)
# Created: 2024-04-02T05:45:16.035674
# PINT_version: 0.9.8+832.g080d65c
# User: docs
# Host: build-23937247-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-aws-x86_64-with-glibc2.35
# Python: 3.11.6 (main, Feb 1 2024, 16:47:41) [GCC 11.4.0]
# Format: pint
PSR 1748-2021E
DILATEFREQ N
DMDATA N
NTOA 0
RAJ 17:48:52.75000000 1 0.04999999999999999584
DECJ -20:21:29.00000000 1 0.40000000000000002220
PMRA 0.0
PMDEC 0.0
PX 0.0
POSEPOCH 53750.0000000000000000
F0 61.48547655400000167 1 5.0000000000000003114e-10
PEPOCH 53750.0000000000000000
TZRMJD 53801.3860512007449870
TZRSITE 1
TZRFRQ 1949.609
The validate function is also integrated into the add_component() function. When adding a component it will validate the timing model by default; however, it can be switched off by setting flag validate=False. We will use this flag in the next section.
Add a component to the timing model
We will add the dispersion component to the timing model. The steps are: 1. Instantiate the Dispersion class 2. Add dispersion instance into the timing model, with validation as False.
Since the dispersion model’s parameter have not set yet, validation would fail. We will validate it after the parameters filled in.
Add parameters and set values
Validate the timing model.
[9]:
dispersion_class = all_components["DispersionDM"]
dispersion = dispersion_class() # Make the dispersion instance.
# Using validate=False here allows a component being added first and validate later.
tm.add_component(dispersion, validate=False)
Let us examine the components in the timing model.
[10]:
# print the components out, DispersionDM should be there.
print("All components in timing model:")
display(tm.components)
print("\n")
print("Delay components in the DelayComponent_list (order matters!):")
# print the delay component order, dispersion should be after the astrometry
display(tm.DelayComponent_list)
All components in timing model:
{'AbsPhase': AbsPhase(
MJDParameter( TZRMJD 53801.3860512007449870 (d) frozen=True),
strParameter( TZRSITE 1 frozen=True),
floatParameter( TZRFRQ 1949.609 (MHz) frozen=True)),
'Spindown': Spindown(
floatParameter( F0 61.48547655400000167 (Hz) +/- 5e-10 Hz frozen=False),
MJDParameter( PEPOCH 53750.0000000000000000 (d) frozen=True)),
'AstrometryEquatorial': AstrometryEquatorial(
MJDParameter( POSEPOCH 53750.0000000000000000 (d) frozen=True),
floatParameter( PX 0.0 (mas) frozen=True),
AngleParameter( RAJ 17:48:52.75000000 (hourangle) +/- 0h00m00.05s frozen=False),
AngleParameter( DECJ -20:21:29.00000000 (deg) +/- 0d00m00.4s frozen=False),
floatParameter( PMRA 0.0 (mas / yr) frozen=True),
floatParameter( PMDEC 0.0 (mas / yr) frozen=True)),
'DispersionDM': DispersionDM(
floatParameter( DM 0.0 (pc / cm3) frozen=True),
floatParameter( DM1 UNSET,
MJDParameter( DMEPOCH UNSET)}
Delay components in the DelayComponent_list (order matters!):
[AstrometryEquatorial(
MJDParameter( POSEPOCH 53750.0000000000000000 (d) frozen=True),
floatParameter( PX 0.0 (mas) frozen=True),
AngleParameter( RAJ 17:48:52.75000000 (hourangle) +/- 0h00m00.05s frozen=False),
AngleParameter( DECJ -20:21:29.00000000 (deg) +/- 0d00m00.4s frozen=False),
floatParameter( PMRA 0.0 (mas / yr) frozen=True),
floatParameter( PMDEC 0.0 (mas / yr) frozen=True)),
DispersionDM(
floatParameter( DM 0.0 (pc / cm3) frozen=True),
floatParameter( DM1 UNSET,
MJDParameter( DMEPOCH UNSET)]
The DM value can be set as we set the parameters above.
[11]:
tm.DM.quantity = 223.9 * u.pc / u.cm**3
tm.DM.frozen = False # Frozen means not fit.
tm.DM.uncertainty = 0.3 * u.pc / u.cm**3
Run validate again and just make sure everything is setup good.
[12]:
tm.validate() # If this fails, that means the DM model was not setup correctly.
Now the dispersion model component is added and you are now set for your analysis.
Delete a component
Deleting a component will remove the component from component list.
[13]:
# Remove by name
tm.remove_component("DispersionDM")
[14]:
display(tm.components)
{'AbsPhase': AbsPhase(
MJDParameter( TZRMJD 53801.3860512007449870 (d) frozen=True),
strParameter( TZRSITE 1 frozen=True),
floatParameter( TZRFRQ 1949.609 (MHz) frozen=True)),
'Spindown': Spindown(
floatParameter( F0 61.48547655400000167 (Hz) +/- 5e-10 Hz frozen=False),
MJDParameter( PEPOCH 53750.0000000000000000 (d) frozen=True)),
'AstrometryEquatorial': AstrometryEquatorial(
MJDParameter( POSEPOCH 53750.0000000000000000 (d) frozen=True),
floatParameter( PX 0.0 (mas) frozen=True),
AngleParameter( RAJ 17:48:52.75000000 (hourangle) +/- 0h00m00.05s frozen=False),
AngleParameter( DECJ -20:21:29.00000000 (deg) +/- 0d00m00.4s frozen=False),
floatParameter( PMRA 0.0 (mas / yr) frozen=True),
floatParameter( PMDEC 0.0 (mas / yr) frozen=True))}
Dispersion model should disappear from the timing model.
Add prefix-style parameters
Prefix style parameters are used in certain models (e.g., DMX_nnnn or Fn).
Let us use the DispersionDMX model to demonstrate how it works.
[15]:
tm.add_component(all_components["DispersionDMX"]())
[16]:
_ = [print(cp) for cp in tm.components]
# "DispersionDMX" should be there.
AbsPhase
Spindown
AstrometryEquatorial
DispersionDMX
Display the existing DMX parameters
What do we have in DMX model.
Note, Component class also has the attribute params, which is only for the parameters in the component.
[17]:
print(tm.components["DispersionDMX"].params)
['DMX', 'DMX_0001', 'DMXR1_0001', 'DMXR2_0001']
Add DMX parameters
Since we already have DMX_0001, we will add DMX_0003, just want to show that for DMX model, DMX index(‘0001’ part) does not have to follow the consecutive order.
The prefix type of parameters have to use prefixParameter class from pint.models.parameter module.
[18]:
# Add prefix parameters
dmx_0003 = p.prefixParameter(
parameter_type="float", name="DMX_0003", value=None, units=u.pc / u.cm**3
)
tm.components["DispersionDMX"].add_param(dmx_0003, setup=True)
# tm.add_param_from_top(dmx_0003, "DispersionDMX", setup=True)
# # Component should given by its name string. use setup=True make sure new parameter get registered.
Check if the parameter and component setup correctly.
[19]:
display(tm.params)
display(tm.delay_deriv_funcs.keys()) # the derivative function should be added.
['PSR',
'TRACK',
'EPHEM',
'CLOCK',
'UNITS',
'START',
'FINISH',
'RM',
'INFO',
'TIMEEPH',
'T2CMETHOD',
'BINARY',
'DILATEFREQ',
'DMDATA',
'NTOA',
'CHI2',
'CHI2R',
'TRES',
'DMRES',
'POSEPOCH',
'PX',
'RAJ',
'DECJ',
'PMRA',
'PMDEC',
'F0',
'PEPOCH',
'TZRMJD',
'TZRSITE',
'TZRFRQ',
'DMX',
'DMX_0001',
'DMXR1_0001',
'DMXR2_0001',
'DMX_0003']
dict_keys(['PX', 'RAJ', 'DECJ', 'PMRA', 'PMDEC', 'DMX_0001', 'DMX_0003'])
However only adding DMX_0003 is not enough, since one DMX parameter also need a DMX range, DMXR1_0003, DMXR2_0003 in this case. Without them, the validation will fail. So let us add them as well.
[20]:
dmxr1_0003 = p.prefixParameter(
parameter_type="mjd", name="DMXR1_0003", value=None, units=u.day
) # DMXR1 is a type of MJD parameter internally.
dmxr2_0003 = p.prefixParameter(
parameter_type="mjd", name="DMXR2_0003", value=None, units=u.day
) # DMXR1 is a type of MJD parameter internally.
tm.components["DispersionDMX"].add_param(dmxr1_0003, setup=True)
tm.components["DispersionDMX"].add_param(dmxr2_0003, setup=True)
[21]:
tm.params
[21]:
['PSR',
'TRACK',
'EPHEM',
'CLOCK',
'UNITS',
'START',
'FINISH',
'RM',
'INFO',
'TIMEEPH',
'T2CMETHOD',
'BINARY',
'DILATEFREQ',
'DMDATA',
'NTOA',
'CHI2',
'CHI2R',
'TRES',
'DMRES',
'POSEPOCH',
'PX',
'RAJ',
'DECJ',
'PMRA',
'PMDEC',
'F0',
'PEPOCH',
'TZRMJD',
'TZRSITE',
'TZRFRQ',
'DMX',
'DMX_0001',
'DMXR1_0001',
'DMXR2_0001',
'DMX_0003',
'DMXR1_0003',
'DMXR2_0003']
Then validate it.
[22]:
tm.validate()
Remove a parameter
Remove a parameter is just use the remove_param() function.
[23]:
tm.remove_param("DMX_0003")
tm.remove_param("DMXR1_0003")
tm.remove_param("DMXR2_0003")
display(tm.params)
['PSR',
'TRACK',
'EPHEM',
'CLOCK',
'UNITS',
'START',
'FINISH',
'RM',
'INFO',
'TIMEEPH',
'T2CMETHOD',
'BINARY',
'DILATEFREQ',
'DMDATA',
'NTOA',
'CHI2',
'CHI2R',
'TRES',
'DMRES',
'POSEPOCH',
'PX',
'RAJ',
'DECJ',
'PMRA',
'PMDEC',
'F0',
'PEPOCH',
'TZRMJD',
'TZRSITE',
'TZRFRQ',
'DMX',
'DMX_0001',
'DMXR1_0001',
'DMXR2_0001']
Add higher order derivatives of spin frequency to timing model
Adding higher order derivatives of spin frequency (e.g., F2, F3, F4) is a common use case. Fn is a prefixParameter, but unlike the DMX_ parameters, all indexes up to the maximum order must exist, since it represents the coefficients of a Taylor expansion.
Let us list the current spindown model parameters:
[24]:
display(tm.components["Spindown"].params)
['F0', 'PEPOCH']
Let us add F1 and F2 to the model. Both F1 and F2 needs a very high precision, we use longdouble=True flag to specify the F2 value to be a longdouble type.
Note, if we add F3 directly without F2, the validation will fail.
[25]:
f1 = p.prefixParameter(
parameter_type="float", name="F1", value=0.0, units=u.Hz / (u.s), longdouble=True
)
f2 = p.prefixParameter(
parameter_type="float",
name="F2",
value=0.0,
units=u.Hz / (u.s) ** 2,
longdouble=True,
)
[26]:
tm.components["Spindown"].add_param(f1, setup=True)
tm.components["Spindown"].add_param(f2, setup=True)
[27]:
tm.validate()
[28]:
display(tm.params)
['PSR',
'TRACK',
'EPHEM',
'CLOCK',
'UNITS',
'START',
'FINISH',
'RM',
'INFO',
'TIMEEPH',
'T2CMETHOD',
'BINARY',
'DILATEFREQ',
'DMDATA',
'NTOA',
'CHI2',
'CHI2R',
'TRES',
'DMRES',
'POSEPOCH',
'PX',
'RAJ',
'DECJ',
'PMRA',
'PMDEC',
'F0',
'PEPOCH',
'F1',
'F2',
'TZRMJD',
'TZRSITE',
'TZRFRQ',
'DMX',
'DMX_0001',
'DMXR1_0001',
'DMXR2_0001']
Now F2 can be used in the timing model.
[29]:
tm.F1.quantity = -1.181e-15 * u.Hz / u.s
tm.F1.uncertainty = 1e-18 * u.Hz / u.s
tm.F2.quantity = 2e-10 * u.Hz / u.s**2
display(tm.F2)
floatParameter( F2 2e-10 (Hz / s2) frozen=True)
[ ]: