search

nanograv / enterprise_extensions

A set of extension codes, utilities, and scripts for the enterprise PTA analysis framework.
MIT License
26 stars 59 forks source link

Keep parameters as float128 values in tm_delay #228

Open mattpitkin opened 6 months ago

mattpitkin commented 6 months ago

In the tm_delay function, the parameters to be varied are set as np.double precision values. The original values in the t2pulsar object (in particular the frequency parameters) are stored as np.float128, so the recasting as double's loses some precision. This isn't generally problematic, but it can be problematic if you are doing simulations with very low noise and truncation of, e.g., the frequency, caused by converting to double's throws off other parameters.

I will show a concrete example. I generate some fake TOAs for J2145-0750 using the Tempo2 fake plugin, with the very low noise of 0.1ns:

tempo2 -gr fake -f J2145-0750.par -ndobs 10 -nobsd 1 -randha n -start 51545 -end 58861 -format tempo2 -rms 1e-7
# make sure file extension is .tim
cp J2145-0750.simulate J2145-0750.tim 
# add some fake -pta and -f flags
sed -i 's|0.00010 7|0.00010 7 -pta None -f None|g' J2145-0750.tim

where the .par file contained:

PSRJ           J2145-0750
ELONG          326.0246148221227000252309125   
ELAT           5.3130536066555000000895259   
F0             62.2958887973344346 1  0.0000000000009236
F1             -1.1561411786e-16 1  8.991138848391D-21
F2             9.3513157514079237246e-34
PEPOCH         55597                       
POSEPOCH       55597                       
DMEPOCH        55597                       
DM             9.001883     
TZRFRQ         1440                     
TZRSITE        7        
EPHVER         5                           
CLK            TT(BIPM2017)
UNITS          TDB
EPHEM          DE405
TIMEEPH        IF99
DILATEFREQ     Y
PLANET_SHAPIRO Y

Now, if I assume that F2 is actually unknown and I want to estimate its value from the TOAs, I could set up a .par, called, say, J2145-0750_no_f2.par file containing:

PSRJ           J2145-0750
ELONG          326.0246148221227000252309125   
ELAT           5.3130536066555000000895259   
F0             62.2958887973344346 1  0.0000000000009236
F1             -1.1561411786e-16 1  8.991138848391D-21
F2             0 1 1
PEPOCH         55597                       
POSEPOCH       55597                       
DMEPOCH        55597                       
DM             9.001883     
TZRFRQ         1440                     
TZRSITE        7        
EPHVER         5                           
CLK            TT(BIPM2017)
UNITS          TDB
EPHEM          DE405
TIMEEPH        IF99
DILATEFREQ     Y
PLANET_SHAPIRO Y

and run (the below assumes https://github.com/nanograv/enterprise_extensions/pull/226 has been applied to allow it to run with white_vary=False):

import numpy as np
from matplotlib import pyplot as plt

from enterprise.pulsar import Pulsar
from enterprise_extensions.models import model_singlepsr_noise

psrname = "J2145-0750"

parfile = "J2145-0750_no_f2.par"
timfile = "J2145-0750.tim"

psr = Pulsar(parfile, timfile, drop_t2pulsar=False)
psr.flags["pta"] = None

# set to estimate F0, F1 and F2 (but actually hold F0 and F1 fixed)
plist = ["F0", "F1", "F2"]

pta = model_singlepsr_noise(
    psr,
    tmparam_list=plist,
    tm_linear=False,
    tm_var=True,
    tm_marg=False,
    red_var=False,
    white_vary=False,
)

# get log-likelihood over f2
tmparamname = "J2145-0750_timing model_tmparams"
curparams = {}
f2 = np.linspace(-3e-33, 3e-33, 100)
ll = np.zeros_like(f2)
for i in range(len(f2)):
   curparams[tmparamname] = [0.0, 0.0, f2[i]]
   ll[i] = pta.get_lnlikelihood(curparams)

fig, ax = plt.subplots()

ax.plot(f2, np.exp(ll - ll.max()))
ax.set_ylabel("likelihood")
ax.set_xlabel("$\ddot{f}$ (Hz/$s^2$)")

fig.show()

it gives:

test_wrong_sign

where the F2 value is completely wrong (note that a fit with Tempo2 itself gets the correct value).

From digging into why this is happening, I find that when the F0 value gets converted from float128 to double in timing.py, it gets truncated from:

62.2958887973344346

to

62.295888797334435

With the very low noise, this difference in frequency is enough to throw off the estimation of F2. If I apply the patch in this PR and run the above code again, I get:

test_right_sign

which has the correct F2 value.

vallis commented 6 months ago

Note that for some platforms (most notably Apple Silicon a.k.a. arm64) float128 it's either unavailable or defaults to a regular double. So on some platforms the change may break the code, or not fix the problem. On those platforms, tempo2 will also operate with double precision internally.

So I'm not sure what the way forward is. Is it a problem if this code breaks on arm64? It would signal that simulation is not doing what one thinks it is, but should we ship broken code? Will it even CI?

mattpitkin commented 6 months ago

I've made a slight change to workaround np.float128 not being present on ARM64 systems.

codecov-commenter commented 6 months ago

Codecov Report

Attention: Patch coverage is 33.33333% with 4 lines in your changes are missing coverage. Please review.

Project coverage is 37.30%. Comparing base (d91baf2) to head (d4afb9d).

Additional details and impacted files [![Impacted file tree graph](https://app.codecov.io/gh/nanograv/enterprise_extensions/pull/228/graphs/tree.svg?width=650&height=150&src=pr&token=HTKG5WS0VZ&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nanograv)](https://app.codecov.io/gh/nanograv/enterprise_extensions/pull/228?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nanograv) ```diff @@ Coverage Diff @@ ## master #228 +/- ## ========================================== + Coverage 37.29% 37.30% +0.01% ========================================== Files 20 20 Lines 3974 3978 +4 ========================================== + Hits 1482 1484 +2 - Misses 2492 2494 +2 ``` | [Files](https://app.codecov.io/gh/nanograv/enterprise_extensions/pull/228?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nanograv) | Coverage Δ | | |---|---|---| | [enterprise\_extensions/timing.py](https://app.codecov.io/gh/nanograv/enterprise_extensions/pull/228?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nanograv#diff-ZW50ZXJwcmlzZV9leHRlbnNpb25zL3RpbWluZy5weQ==) | `32.00% <33.33%> (+3.42%)` | :arrow_up: | ------ [Continue to review full report in Codecov by Sentry](https://app.codecov.io/gh/nanograv/enterprise_extensions/pull/228?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nanograv). > **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nanograv) > `Δ = absolute (impact)`, `ø = not affected`, `? = missing data` > Powered by [Codecov](https://app.codecov.io/gh/nanograv/enterprise_extensions/pull/228?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nanograv). Last update [d91baf2...d4afb9d](https://app.codecov.io/gh/nanograv/enterprise_extensions/pull/228?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nanograv). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=nanograv).