This Jupyter notebook can be downloaded from How_to_build_a_timing_model_component.ipynb, or viewed as a python script at How_to_build_a_timing_model_component.py.
How to compose a timing model component
PINT’s design makes it easy to add a new, custom timing model component to meet specific needs. This notebook demonstrates how to write your own timing model component with the minimal requirements so that PINT can recognize and use it in fits. Here, we implement a new spindown class, PeriodSpindown
, that is parameterized by P0
, P1
, instead of the built-in Spindown
model component, which uses F0
, F1
.
Building the timing model component from scratch
This example notebook includes the following contents * Defining a timing model component class * Necessary parts * Conventions * Use it with the TimingModel
class * Add the new component to the TimingModel
class * Use the functions in the TimingModel
class to interact with the new component.
We will build a simple model component, pulsar spindown model with spin period as parameters, instead of spin frequency.
Import the necessary modules
[1]:
# PINT uses astropy units in the internal calculation and is highly recommended for a new component
import astropy.units as u
# Import the component classes.
from pint.models.spindown import SpindownBase
import pint.models.parameter as p
import pint.config
import pint.logging
# setup the logging
pint.logging.setup(level="INFO")
[1]:
1
Define the timing model class
A timing model component should be an inheritance/subclass of pint.models.timing_model.Component
. PINT also pre-defines three component subclasses for the most used type of components and they have different attribute and functions (see: https://nanograv-pint.readthedocs.io/en/latest/api/pint.models.timing_model.html): * DelayComponent for delay type of models. * PhaseComponent for phase type of models. * NoiseComponent for noise type of models.
Here since we are making a spin-down model, we will use the PhaseComponent
.
Required parts
Model parameters, generally defined as
PINT.models.parameter.Parameter
class or its subclasses. (see https://nanograv-pint.readthedocs.io/en/latest/api/pint.models.parameter.html)Model functions, defined as methods in the component, including:
.setup(), for setting up the component(e.g., registering the derivatives).
.validate(), for checking if the parameters have the correct inputs.
Modeled quantity functions.
The derivative of modeled quantities.
Other support functions.
Conventions
To make a component work as a part of a timing model, it has to follow the following rules to interface the TimingModel
class. Using the analog of a circuit board, the TimingModel
object is the mother board, and the Component
objects are the electronic components(e.g., resistors and transistors); and the following rules are the pins of a component.
Set the class attribute
.register
to be True so that the component is in the searching space of model builderAdd the method of final result in the designated list, so the
TimingModel
’s collecting function(e.g., total delay or total phase) can collect the result. Here are the designated list for the most common component type:DelayComponent: .delay_funcs_component
PhaseComponent: .phase_funcs_component
NoiseComponent: .
.basis_funcs
.covariance_matrix_funcs
.scaled_toa_sigma_funcs
.scaled_dm_sigma_funcs
.dm_covariance_matrix_funcs_component
Register the analytical derivative functions using the
.register_deriv_funcs(derivative function, parameter name)
if any.If one wants to access the attribute in the parent
TimingModel
class or from other components, please use._parent
attribute which is a linker to theTimingModel
class and other components.
[2]:
class PeriodSpindown(SpindownBase):
"""This is an example model component of pulsar spindown but parametrized as period."""
register = True # Flags for the model builder to find this component.
# define the init function.
# Most components do not have a parameter for input.
def __init__(self):
# Get the attributes that initialized in the parent class
super().__init__()
# Add parameters using the add_params in the TimingModel
# Add spin period as parameter
self.add_param(
p.floatParameter(
name="P0",
value=None,
units=u.s,
description="Spin period",
longdouble=True,
tcb2tdb_scale_factor=u.Quantity(1),
)
)
# Add spin period derivative P1. Since it is not all rquired, we are setting the
# default value to 0.0
self.add_param(
p.floatParameter(
name="P1",
value=0.0,
units=u.s / u.s,
description="Spin period derivative",
longdouble=True,
tcb2tdb_scale_factor=u.Quantity(1),
)
)
# Add reference epoch time.
self.add_param(
p.MJDParameter(
name="PEPOCH_P0",
description="Reference epoch for spin-down",
time_scale="tdb",
tcb2tdb_scale_factor=u.Quantity(1),
)
)
# Add spindown phase model function to phase functions
self.phase_funcs_component += [self.spindown_phase_period]
# Add the d_phase_d_delay derivative to the list
self.phase_derivs_wrt_delay += [self.d_spindown_phase_period_d_delay]
def setup(self):
"""Setup the model. Register the derivative functions"""
super().setup() # This will run the setup in the Component class.
# The following lines are resgistering the derivative functions to the timingmodel.
self.register_deriv_funcs(self.d_phase_d_P0, "P0")
self.register_deriv_funcs(self.d_phase_d_P1, "P1")
def validate(self):
"""Check the parameter value."""
super().validate() # This will run the .validate() in the component class
# Check required parameters, since P1 is not required, we are not checking it here
for param in ["P0"]:
if getattr(self, param) is None:
raise ValueError(f"Spindown period model needs {param}")
# One can always setup properties for updating attributes automatically.
@property
def F0(self):
# We return F0 as parameter here since the other place of PINT code use F0
# in the format of PINT parameter.
return p.floatParameter(
name="F0",
value=1.0 / self.P0.quantity,
units="Hz",
description="Spin-frequency",
long_double=True,
tcb2tdb_scale_factor=u.Quantity(1),
)
# Defining the derivatives. In the PINT, a common format of derivative naming is
# d_xxx_d_xxxx
@property
def d_F0_d_P0(self):
return -1.0 / self.P0.quantity**2
@property
def F1(self):
return p.floatParameter(
name="F1",
value=self.d_F0_d_P0 * self.P1.quantity,
units=u.Hz / u.s,
description="Spin down frequency",
long_double=True,
tcb2tdb_scale_factor=u.Quantity(1),
)
@property
def d_F1_d_P0(self):
return self.P1.quantity * 2.0 / self.P0.quantity**3
@property
def d_F1_d_P1(self):
return self.d_F0_d_P0
def get_dt(self, toas, delay):
"""dt from the toas to the reference time."""
# toas.table['tdbld'] stores the tdb time in longdouble.
return (toas.table["tdbld"] - self.PEPOCH_P0.value) * u.day - delay
# Defining the phase function, which is added to the self.phase_funcs_component
def spindown_phase_period(self, toas, delay):
"""Spindown phase using P0 and P1"""
dt = self.get_dt(toas, delay)
return self.F0.quantity * dt + 0.5 * self.F1.quantity * dt**2
def d_spindown_phase_period_d_delay(self, toas, delay):
"""This is part of the derivative chain for the parameters in the delay term."""
dt = self.get_dt(toas, delay)
return -(self.F0.quantity + dt * self.F1.quantity)
def d_phase_d_P0(self, toas, param, delay):
dt = self.get_dt(toas, delay)
return self.d_F0_d_P0 * dt + 0.5 * self.d_F1_d_P0 * dt**2
def d_phase_d_P1(self, toas, param, delay):
dt = self.get_dt(toas, delay)
return 0.5 * self.d_F1_d_P1 * dt**2
Apply the new component to the TimingModel
Let us use this new model component in our example pulsar “NGC6440E”, which has F0
and F1
. Instead, we will use the model component above. The following .par
file string if converted from the NGC6440E.par
with P0
and P1
instead of F0
, F1
.
[3]:
par_string = """
PSR 1748-2021E
RAJ 17:48:52.75 1 0.05
DECJ -20:21:29.0 1 0.4
P0 0.016264003404474613 1 0
P1 3.123955D-19 1 0
PEPOCH_P0 53750.000000
POSEPOCH 53750.000000
DM 223.9 1 0.3
SOLARN0 0.00
EPHEM DE421
UNITS TDB
TIMEEPH FB90
CORRECT_TROPOSPHERE N
PLANET_SHAPIRO N
DILATEFREQ N
TZRMJD 53801.38605120074849
TZRFRQ 1949.609
TZRSITE 1
"""
[4]:
import io
from pint.models import get_model
Load the timing model with new parameterization.
[5]:
model = get_model(
io.StringIO(par_string)
) # PINT can take a string IO for inputing the par file
Check if the component is loaded into the timing model and make sure there is no built-in spindown model.
[6]:
print(model.components["PeriodSpindown"])
print(
"Is the built-in spin-down model in the timing model: ",
"Spindown" in model.components.keys(),
)
print("Is 'P0' in the timing model: ", "P0" in model.params)
print("Is 'P1' in the timing model: ", "P1" in model.params)
print("Is 'F0' in the timing model: ", "F0" in model.params)
print("Is 'F1' in the timing model: ", "F1" in model.params)
PeriodSpindown(
floatParameter( P0 0.016264003404474613 (s) +/- 0.0 s frozen=False),
floatParameter( P1 3.123955e-19 () +/- 0.0 frozen=False),
MJDParameter( PEPOCH_P0 53750.0000000000000000 (d) frozen=True))
Is the built-in spin-down model in the timing model: False
Is 'P0' in the timing model: True
Is 'P1' in the timing model: True
Is 'F0' in the timing model: False
Is 'F1' in the timing model: False
Load TOAs and prepare for fitting
[7]:
from pint.fitter import WLSFitter
from pint.toa import get_TOAs
[8]:
toas = get_TOAs(pint.config.examplefile("NGC6440E.tim"), ephem="DE421")
f = WLSFitter(toas, model)
Plot the residuals
[9]:
import matplotlib.pyplot as plt
Plot the prefit residuals.
[10]:
plt.errorbar(
toas.get_mjds().value,
f.resids_init.time_resids.to_value(u.us),
yerr=toas.get_errors().to_value(u.us),
fmt=".",
)
plt.title(f"{model.PSR.value} Pre-Fit Timing Residuals")
plt.xlabel("MJD")
plt.ylabel("Residual (us)")
plt.grid()
Fit the TOAs using P0
and P1
[11]:
f.fit_toas()
[11]:
np.longdouble('59.574712573297878792')
Plot the post-fit residuals
[12]:
plt.errorbar(
toas.get_mjds().value,
f.resids.time_resids.to_value(u.us),
yerr=toas.get_errors().to_value(u.us),
fmt=".",
)
plt.title(f"{model.PSR.value} Pre-Fit Timing Residuals")
plt.xlabel("MJD")
plt.ylabel("Residual (us)")
plt.grid()
Print out the summary
[13]:
f.print_summary()
Fitted model using weighted_least_square method with 5 free parameters to 62 TOAs
Prefit residuals Wrms = 1090.5803879598564 us, Postfit residuals Wrms = 21.18210858178613 us
Chisq = 59.575 for 56 d.o.f. for reduced Chisq of 1.064
PAR Prefit Postfit Units
=================== ==================== ============================ =====
PSR 1748-2021E 1748-2021E None
EPHEM DE421 DE421 None
CLOCK TT(BIPM2023) None
UNITS TDB TDB None
START 53478.3 d
FINISH 54187.6 d
TIMEEPH FB90 FB90 None
DILATEFREQ N None
DMDATA N None
NTOA 0 None
CHI2 59.5747
CHI2R 1.06383
TRES 21.1821 us
POSEPOCH 53750 d
PX 0 mas
RAJ 17h48m52.75s 17h48m52.80034691s +/- 0.00014 hourangle_second
DECJ -20d21m29s -20d21m29.38331083s +/- 0.033 arcsec
PMRA 0 mas / yr
PMDEC 0 mas / yr
CORRECT_TROPOSPHERE N None
PLANET_SHAPIRO N None
NE_SW 0 1 / cm3
SWP 2
SWM 0 None
DM 223.9 224.114(35) pc / cm3
TZRMJD 53801.4 d
TZRSITE 1 1 None
TZRFRQ 1949.61 MHz
P0 0.016264 0.016264003404376(5) s
P1 3.12396e-19 3.125(4)×10⁻¹⁹
PEPOCH_P0 53750 d
Derived Parameters:
Period = 0.01626400340437608±0 s
Pdot = (3.1248326935082824±0)×10⁻¹⁹
Characteristic age = 8.246e+08 yr (braking index = 3)
Surface magnetic field = 2.28e+09 G
Magnetic field at light cylinder = 4806 G
Spindown Edot = 2.868e+33 erg / s (I=1e+45 cm2 g)
Write out a par file for the result
[14]:
f.model.write_parfile("/tmp/output.par")
print(f.model.as_parfile())
# Created: 2024-11-05T15:45:43.086798
# PINT_version: 1.1
# User: docs
# Host: build-26181900-project-85767-nanograv-pint
# OS: Linux-5.19.0-1028-aws-x86_64-with-glibc2.35
# Python: 3.11.9 (main, Jun 18 2024, 09:40:25) [GCC 11.4.0]
# Format: pint
PSR 1748-2021E
EPHEM DE421
CLOCK TT(BIPM2023)
UNITS TDB
START 53478.2858714195382639
FINISH 54187.5873241702319097
TIMEEPH FB90
DILATEFREQ N
DMDATA N
NTOA 62
CHI2 59.57471257329788
CHI2R 1.063834153094605
TRES 21.182108581786131332
RAJ 17:48:52.80034691 1 0.00013524660915124006
DECJ -20:21:29.38331083 1 0.03285153312806998044
PMRA 0.0
PMDEC 0.0
PX 0.0
POSEPOCH 53750.0000000000000000
CORRECT_TROPOSPHERE N
PLANET_SHAPIRO N
SOLARN0 0.0
SWM 0
DM 224.11379619426721291 1 0.034938980494125089493
TZRMJD 53801.3860512007484954
TZRSITE 1
TZRFRQ 1949.609
P0 0.01626400340437608 1 4.784091376106004e-15
P1 3.1248326935082817e-19 1 3.813960677879849e-22
PEPOCH_P0 53750.0000000000000000