[set-up], imports¶

See README for more details about installing necessary packages and setting up your environment.

In [ ]:
 
In [1]:
import pint_pal as pp
import pint_pal.lite_utils as lu
import pint_pal.plot_utils as pu
from pint_pal.timingconfiguration import TimingConfiguration
from pint_pal.ftester import run_Ftests
from astropy import log
from pint.fitter import ConvergenceFailure
import pint.fitter
import os
import copy
from astropy.visualization import quantity_support
quantity_support()

# notebook gives interactive plots but not until the kernel is done
%matplotlib notebook
# inline gives non-interactive plots right away
#%matplotlib inline

# Set logging level (PINT uses loguru)
log.setLevel("INFO") # Set desired verbosity of log statements (DEBUG/INFO/WARNING/ERROR)
pint.logging.setup(level="WARNING", usecolors=True)
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Out[1]:
1
Fatal Python error: config_get_locale_encoding: failed to get the locale encoding: nl_langinfo(CODESET) failed
Python runtime state: preinitialized

Pre or Post Noise solution?¶

In [3]:
hasnoise = False
if hasnoise == False:
    ext = '_prenoise'
else:
    ext = '_postnoise'
In [ ]:
 

Load/update timing solution¶

Load configuration (.yaml) file, get TOAs and timing model; if you're running from the root of the git distribution, simply edit the .yaml file name, otherwise include relevant paths to the .yaml file, and .par/.tim directories as kwargs (see commented example).

In [5]:
config = "sample.dr3.yaml"  # fill in actual path
par_directory = None   # default location
tim_directory = None   # default location
tc = TimingConfiguration(config, par_directory=par_directory, tim_directory=tim_directory)

# To combine TOAs, assumption is that cuts have already been applied properly
mo,to = tc.get_model_and_toas(apply_initial_cuts=False,usepickle=False)

# Uncomment this line to manually excise TOAs.
#tc.manual_cuts(to)

# Computing pulse numbers ensures param changes in the model will not break phase connection
to.compute_pulse_numbers(mo)

# Set non-binary epochs to the center of the data span
lu.center_epochs(mo,to)

# Summarize TOAs present
to.print_summary()
WARNING  (pint.models.noise_model       ): 'TNEQ1 -sys ['JBO.DFB.1520']' is provided by parameter EQUAD, using EQUAD instead. 
WARNING  (pint.models.noise_model       ): 'TNEQ2 -sys ['EFF.EBPP.1360']' is provided by parameter EQUAD, using EQUAD instead. 
WARNING  (pint.models.noise_model       ): 'TNEQ3 -sys ['EFF.EBPP.1410']' is provided by parameter EQUAD, using EQUAD instead. 
WARNING  (pint.models.noise_model       ): 'TNEQ4 -sys ['EFF.EBPP.2639']' is provided by parameter EQUAD, using EQUAD instead. 
WARNING  (pint.models.noise_model       ): 'TNEQ5 -sys ['NRT.BON.1400']' is provided by parameter EQUAD, using EQUAD instead. 
WARNING  (pint.models.noise_model       ): 'TNEQ6 -sys ['NRT.BON.1600']' is provided by parameter EQUAD, using EQUAD instead. 
WARNING  (pint.models.noise_model       ): 'TNEQ7 -sys ['NRT.BON.2000']' is provided by parameter EQUAD, using EQUAD instead. 
WARNING  (pint.models.noise_model       ): 'TNEQ8 -sys ['WSRT.P1.323.C']' is provided by parameter EQUAD, using EQUAD instead. 
WARNING  (pint.models.noise_model       ): 'TNEQ9 -sys ['WSRT.P1.367.C']' is provided by parameter EQUAD, using EQUAD instead. 
INFO: Par file created: 2023-03-17T15:20:55.686266 [pint_pal.timingconfiguration]
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 0: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 1: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 2: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 3: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 4: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 5: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 6: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 7: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 8: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 11: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 12: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 13: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 14: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 15: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 16: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 17: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 18: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 19: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 20: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 21: inserting NaNs
WARNING  (pint.toa                      ): 'pulse_number' not present in data set 22: inserting NaNs
Number of TOAs:  26105
Number of commands:  [3, 5, 7, 1, 1, 1, 1, 2, 2, 2, 1, 2, 2, 1, 2, 3, 2, 2, 2, 3, 3, 3, 1]
Number of observatories: 12 ['parkes', 'leap', 'effelsberg_asterix', 'wsrt', 'gbt', 'ncyobs', 'effelsberg', 'gmrt', 'meerkat', 'nancay', 'jbroach', 'jodrell']
MJD span:  50460.444 to 59644.026
Date span: 1997-01-12 10:39:32.749631210 to 2022-03-06 00:37:27.443521010
effelsberg TOAs (137):
  Min freq:      1348.029 MHz
  Max freq:      2634.249 MHz
  Min error:     0.093 us
  Max error:     23.2 us
  Median error:  1.02 us
effelsberg_asterix TOAs (230):
  Min freq:      1292.221 MHz
  Max freq:      4881.875 MHz
  Min error:     0.099 us
  Max error:     15.5 us
  Median error:  0.964 us
gbt TOAs (17745):
  Min freq:      724.687 MHz
  Max freq:      1880.945 MHz
  Min error:     0.04 us
  Max error:     14.3 us
  Median error:  1.01 us
gmrt TOAs (193):
  Min freq:      305.469 MHz
  Max freq:      1410.108 MHz
  Min error:     0.268 us
  Max error:     15.8 us
  Median error:  1.18 us
jbroach TOAs (174):
  Min freq:      1410.091 MHz
  Max freq:      1676.311 MHz
  Min error:     0.06 us
  Max error:     40.3 us
  Median error:  1.12 us
jodrell TOAs (25):
  Min freq:      1520.000 MHz
  Max freq:      1520.000 MHz
  Min error:     0.42 us
  Max error:     21.7 us
  Median error:  1.01 us
leap TOAs (87):
  Min freq:      1388.000 MHz
  Max freq:      1420.000 MHz
  Min error:     0.031 us
  Max error:     4.56 us
  Median error:  0.221 us
meerkat TOAs (817):
  Min freq:      910.395 MHz
  Max freq:      1648.321 MHz
  Min error:     0.057 us
  Max error:     8.1 us
  Median error:  0.567 us
nancay TOAs (183):
  Min freq:      1397.000 MHz
  Max freq:      2298.305 MHz
  Min error:     0.046 us
  Max error:     7.36 us
  Median error:  0.698 us
ncyobs TOAs (1002):
  Min freq:      1289.836 MHz
  Max freq:      2668.086 MHz
  Min error:     0.045 us
  Max error:     24.4 us
  Median error:  1.19 us
parkes TOAs (5401):
  Min freq:      662.350 MHz
  Max freq:      3852.400 MHz
  Min error:     0.029 us
  Max error:     4.63 us
  Median error:  0.886 us
wsrt TOAs (111):
  Min freq:      323.750 MHz
  Max freq:      1400.000 MHz
  Min error:     0.08 us
  Max error:     15 us
  Median error:  2.6 us

In [6]:
# Define the fitter object and plot pre-fit residuals
fo = tc.construct_fitter(to,mo)
pu.plot_residuals_time(fo, restype='prefit', legend=False)
if mo.is_binary:
    pu.plot_residuals_orb(fo, restype='prefit', legend=False)
In [ ]:
# Set free params based on list in the config file (want to update JUMP handling differently soon)
fo.model.free_params = tc.get_free_params(fo)

# Do the fit
try:
    fo.fit_toas(maxiter=tc.get_niter())
    fo.model.CHI2.value = fo.resids.chi2
except ConvergenceFailure:
    run_Ftest = False
    log.warning('Failed to converge; moving on with best result, but should address before final version.')
In [ ]:
# Plot post-fit residuals, print summary of results, write prenoise solution
pu.plot_residuals_time(fo, restype='postfit', legend=True)
if mo.is_binary:
    pu.plot_residuals_orb(fo, restype='postfit', legend=True)
    
fo.print_summary()
lu.check_convergence(fo)

lu.write_par(fo,toatype=tc.get_toa_type(),addext=ext)
In [1]:
# Look at the whitened residuals as well
if hasnoise == True:
    pu.plot_residuals_time(fo, restype='postfit', whitened=True, legend=False)
    if mo.is_binary:
        pu.plot_residuals_orb(fo, restype='postfit', whitened=True, legend=False)
  Cell In[1], line 7
    continue
    ^
SyntaxError: 'continue' not properly in loop

Compare to one other model¶

You can use this cell to quickly check if/which any parameters have changed by more than 3 sigma (number of sigma is tunable).

Best use case: compare to a pre-noise model or a previous run. I'm working on a more comprehensive comparison notebook for things like combined vs. individual PTA.

In [ ]:
compmodel = lu.compare_models(fo,
               model_to_compare=tc.get_compare_model(),
               verbosity='min',
               nodmx=True,
               threshold_sigma=3.)

Check for additional significant parameters (F-test)¶

This is not likely to be meaningful without noise parameters and it's time consuming -- so we really don't recommend running it until you have some noise parameters. Without that info, it'll be pretty meaningless.

In [3]:
run_Ftest = False
Ftest_dict = None
if run_Ftest:
    savedLevel = log.getEffectiveLevel()
    try:
        log.setLevel("WARNING")
        Ftest_dict = run_Ftests(fo, 
                                alpha=0.0027, 
                                NITS=tc.get_niter())
    finally:
        log.setLevel(savedLevel)

[changelog] entries¶

New changelog entries in the .yaml file should follow a specific format and are only added for specified reasons (excising TOAs, adding/removing params, changing binary models, etc.). For more detailed instructions, run lu.new_changelog_entry? in a new cell. This function can be used to format your entry, which should be added to the bottom of the appropriate .yaml file. Note: make sure your git user.email is properly configured, since this field is used to add your name to the entry.

In [ ]:
lu.new_changelog_entry?
In [ ]: