Unit tests failing about pint

[paulray]

Currently, tests/test_bt.py is failing because the delays computed by the PINT BT model are the opposite sign to those saved in tests/J1955_ltdelays.dat. This was written by @matteobachetti Which sign is correct? Maybe @vhaasteren knows?

[vhaasteren]

The class BT(TimingModel) from pint/models/bt.py is what is being used. That delay function produces delays so that the residuals are ok (see the notebook PINT_example_20150720.ipynb) in the root directory (EDIT no, it's not there. Get it here: http://www.tapir.caltech.edu/~vhaasteren/PINT_example_20150720.ipynb). So the delay function in BT(TimingModel) is ok.

I don't know how the file J1955_ltdelays.dat was made? I assume that lt stands for libstempo? In that case, the value 'torb' of tempo2 is returned in the function 'binarydelay' of libstempo. In tempo2, the torb is applied to the bbat as follows:

        psr[p].obsn[i].torb = torb; // save for plotting etc
        deltaT = (psr[p].obsn[i].bbat-psr[p].param[param_pepoch].val[0])*86400.0+torb;
        fct = psr[p].obsn[i].bbat-(int)psr[p].obsn[i].bbat - (psr[p].param[param_pepoch].val[0] - 
                (int)psr[p].param[param_pepoch].val[0]);       
        ftpd = fct+torb/86400.0;       
        psr[p].obsn[i].pet = psr[p].obsn[i].bbat - torb/SECDAY;

(line 239 of formResiduals.C)

In PINT, the delay is applied to the residuals in the phase function. For spindown, for instance:

    dt_tzrmjd = (toas['tdbld'] - self.TZRMJDld) * SECS_PER_DAY - delay

(line 104 of pint/models/spindown.py). Seems the same to me. So it kind of makes no sense to me why the sign would be different. Would be great if @matteobachetti could comment...

EDIT: Oh, I see that the file PINT_example_20150720.ipynb is not in the master branch yet. Get it here: http://www.tapir.caltech.edu/~vhaasteren/PINT_example_20150720.ipynb)

[paulray]

I was having a failure in the test_model.py program that turned out to be due to a segfault in the general2 plugin for Tempo2 that it was using. I fixed the bug (introduced recently by Mike Keith, I think) and submitted a pull request to the tempo2 repo on bitbucket.

I'll keep this issue open until ./run-tests.sh goes through with no errors, then we should make sure not to accept any pull requests that break the tests. Sadly, it is easy for changes to numpy, astropy, TEMPO, Tempo2 to break our tests without us doing anything...

[matteobachetti]

Hi,
I wrote that test with the old version of the BT model. I'm not sure what changed in the meantime, it's been several months since I last touched that test. I suspect I just used the wrong sign convention in the old version of bt, and used the same wrong convention in libstempo.

[paulray]

That's what I was guessing @matteobachetti . Do you have the code that generated J1955_ltdelays.dat? Maybe you can add that to the repo, or post it here so we can generate a new version with the sign convention matching what @vhaasteren specified above?

[matteobachetti]

Hi,
I lost the script in a disk failure, but I found this in an old backup. I can't test it right now because my libstempo installation has gone nuts. Could you try it out?

def produce_libstempo_delays():
    """Use simulated data of J1955 with TEMPO2 and test if reproducible"""

    parfile = 'tests/J1955.par'
    timfile = 'tests/J1955.tim'

    toas = toa.get_TOAs(timfile, planets=False, usepickle=False)

    SECS_PER_DAY = 86400

    toasld = toas.table['tdbld']
    t0 = toasld[0]
    toaslds = (toasld - t0) * SECS_PER_DAY

    newmodel = bt.BT()

    newmodel.read_parfile(parfile)

    phases = newmodel.phase(toas.table)

    t0 = time.time()
    delays = newmodel.delay(toas.table)

    psr = lt.tempopulsar(parfile, timfile)

    btdelay = psr.binarydelay()
    bttoasld = psr.toas()
    btt0 = bttoasld[0]
    bttoaslds = (bttoasld - btt0) * SECS_PER_DAY

    with open('libstempodelays.txt', 'w') as fobj:
        for i in range(len(delays)):
            print("{:.20} {:.20}".format(toaslds[i], btdelay[i]), file=fobj)

    plt.figure("Delays")
    plt.plot(toaslds, delays, label='PINT')
    plt.scatter(toaslds, delays)

    plt.plot(bttoaslds, btdelay, label='T2')
    plt.scatter(bttoaslds, btdelay)

    plt.legend()
    plt.xlabel('Time (d)')
    plt.ylabel('Delay (s)')

    plt.figure("Diff")

    plt.plot(toaslds, delays - btdelay)
    plt.scatter(toaslds, delays - btdelay)

    plt.legend()
    plt.xlabel('Time (d)')
    plt.ylabel('Diff (s)')
    plt.show()

[paulray]

I was able to run the above code, after I commented out the t0 = time.time() line that I didn't understand. I also added the following header:

from __future__ import print_function, division
import matplotlib.pyplot as plt
import libstempo as lt
import pint.toa as toa
import pint.models.bt as bt

When I run it, PINT produces a nice set of delays, but libstempo just produces 0.0s.

Testing it piece by piece, it seems that libstempo is disabling delays for some reason:

In [20]:  psr = lt.tempopulsar(parfile, timfile)
WARNING [TIM1]: Please place MODE flags in the parameter file 
Delay correction turned off for psr 0

Anyone know why that might be? I've never actually used libstempo.

[vhaasteren]

Paul, you need to have tempo2 first run formResiduals. Before running binarydelay, run:

_ = psr.residuals()

That should make it work.

[paulray]

That indeed fixed it. However, it is the case that that the pint line
delays = newmodel.delay(toas.table)
and the libstempo line
btdelay = psr.binarydelay()
are returning values with opposite sign. Here are the first 10 entries of each:

In [25]: btdelay[:10]
Out[25]: 
array([ 31.271713,  31.065,  30.769283,  30.385405,  29.914461,  29.357797,
        28.717004,  27.993915,  27.190596,  26.309347], dtype=float128)

In [26]: delays[:10]
Out[26]: 
array([-31.27171272, -31.06499991, -30.76928284, -30.38540461,
       -29.91446101, -29.35779727, -28.7170044 , -27.99391453,
       -27.190596  , -26.309347  ])

So, we just need to sort out why the sign conventions are different.

[paulray]

OK, at least on my computer pull request #61 fixes these issues, with the exception of pickling TOAs. That fix is coming in astropy 1.0.7. I don't want to install the dev version of astropy on my machine since astropy is installed by MacPorts and there are several other packages I use that depend on astropy.
The fix for the btdelay sign convention is just a hack. We should understand, and document somewhere, what the proper sign convention is for binary delays in libstempo and PINT.

[matteobachetti]

Oh, sorry for that time.time() call, I had some execution time measurements that I didn't clean up properly.
For these astropy version issues: I usually use anaconda and create virtual environments with different combinations of library versions, it's very convenient. Have you tried it?

[paulray]

I keep hearing about virtualenv. I should give it a shot. But, it looks like astropy 1.1 is about to release and it will fix the pickling bug, as well as adding better support for units. In particular, matplotlib will be able to plot arrays with units. Very nice!
http://astropy.readthedocs.org/en/v1.1rc1/whatsnew/1.1.html

[paulray]

Pull request #61 was merged, so I'm closing this issue.