Disable TZRMJD corrections if TZRSITE is bat

[gcerib]

This is a pull request following the discussion with @abhisrkckl , aimed at resolving #1691.
It introduces a check on TZRSITE in absolute_phase.py, and disables clock corrections on TZRMJD if the former is one of the aliases for the SSB. This is necessary for PINT to consistently reproduce results by tempo2, which implements such a check.

I managed to run tox to test the changes, it somehow took 32 hours on my server. It fails on tests/test_precision.py::test_pulsar_mjd_close_to_mjd_on_leap_second_days . However, so does the version I originally checked out, prior to applying any change to absolute_phase.py . Please find attached the short summary I obtained from both the original (unchanged) version and the one with changes.

tox_result.txt
tox_result__orig.txt

Cheers,
GC

[abhisrkckl]

Can you please add a test to verify that your original issue has been resolved?

[abhisrkckl]

The barycenter has more aliases than "@" and bat. It may be better to change this comparison to
self.TZRSITE.value.lower() not in [get_observatory("@" ).name] + get_observatory("@" ).aliases

[abhisrkckl]

Also, please add a similar condition for the geocenter ("geo") as well.

[gcerib]

Hi. In tempo2's https://bitbucket.org/psrsoft/tempo2/src/master/refphs.C, this is what it does:

 21             if (obs->telID[0]=='@'            /* Have barycentric arrival time */
 22                     || strcasecmp(obs->telID,"bat")==0)
 23             {
 24                 obs->clockCorr=0;  /* therefore don't do clock corrections */
 25                 obs->delayCorr=0;
 26             }
 27             else if (strcmp(obs->telID,"STL")==0)
 28             {
 29                 obs->clockCorr=0;  /* don't do clock corrections */
 30                 obs->delayCorr=1;
 31             }
 32             else if (strcmp(obs->telID,"STL_FBAT")==0)
 33             {
 34                 obs->clockCorr=0;  /* don't do clock corrections */
 35                 obs->delayCorr=1;
 36             }
 37             else
 38             {
 39                 obs->clockCorr=1;
 40                 obs->delayCorr=1;
 41             }

It doesn't involve geocentric TOAs at all, and indeed, a par file with "COE" as TZRSITE will get the clock corrections in both PINT and tempo2. If we now disable them in PINT, I am afraid we will end up in the mirrored situation in which PINT doesn't apply the clock corrections to geocentric TZRMJD, whereas tempo2 does, producing again absolute phase differences between the two packages.

Regarding the aliases for "@" and "BAT", this can be easily changed. I have little experience with pytest, but I can try and build a test by looking at other ones in the tests folder.

[abhisrkckl]

OK... Then you can ignore the part about the geocenter.

About STL and STL_FBAT, I don't know what they are. They are not supported by PINT in any case.

[abhisrkckl]

I dug into this a bit more. The include_bipm and include_gps attributes of the BarycenterObs class are explicitly set to False in https://github.com/nanograv/PINT/blob/master/src/pint/observatory/special_locations.py

So this bug should not have occurred in the first place.

I think the bug is in the get_observatory() function. When I do print(get_observatory("@").include_bipm), I get True. This is because the get_observatory() function wrongly sets this from default arguments. I think the correction should be in get_observatory(). This special case treatment is not actually needed.

[gcerib]

Thanks. If I set defaults to None and reset them only if they are provided, as you suggest, it doesn't actually work, because in pint.toa get_observatory() gets called with explicit arguments:

2206         for obs, grp in self.get_obs_groups():
2207             site = get_observatory(
2208                 obs,
2209                 include_gps=include_gps,
2210                 include_bipm=include_bipm,
2211                 bipm_version=bipm_version,
2212             )

Any idea?

[abhisrkckl]

I meant changing the get_observatory function itself to something like this:

def get_observatory(
    name, include_gps=True, include_bipm=True, bipm_version=bipm_default
):
    """Convenience function to get observatory object with options.

    This function will simply call the ``Observatory.get`` method but
    will manually modify the global observatory object after the method is called.
    Name-matching is case-insensitive.

    If the observatory is not present in the PINT list, will fallback to astropy.

    Parameters
    ----------
    name : str
        The name of the observatory
    include_gps : bool, optional
        Set False to disable UTC(GPS)->UTC clock correction.
    include_bipm : bool, optional
        Set False to disable TAI TT(BIPM) clock correction.
    bipm_version : str, optional
        Set the version of TT BIPM clock correction files.

    .. note:: This function can and should be expanded if more clock
        file switches/options are added at a public API level.

    """
    site = Observatory.get(name)
    if include_gps is not None:
        site.include_gps = include_gps
    if include_bipm is not None:
        site.include_bipm = include_bipm
        site.bipm_version = bipm_version
    return site

[gcerib]

Hi,
Yes. That is precisely what I did. It still fails, even if I set the defaults of get_observatory to None. This is because toa.py:get_TOAs() sets include_bipm and include_gps to True even when it is called without them (and should therefore default to None), if all TOAs provided do not have yet the clock corrections (as is the case):

 283     if all("clkcorr" not in f for f in t.table["flags"]):
 284         if include_gps is None:
 285             include_gps = True
 286         if bipm_version is None:
 287             bipm_version = bipm_default
 288         if include_bipm is None:
 289             include_bipm = True
 290         # FIXME: should we permit existing clkcorr flags?
 291         t.apply_clock_corrections(
 292             include_gps=include_gps,
 293             include_bipm=include_bipm,
 294             bipm_version=bipm_version,
 295             limits=limits,
 296         )

The conditions in lines 284 and 288 are met by default, i.e. when calling get_TOAs without explicitly setting include_gps and include_bipm, and those parameters are passed over to apply_clock_corrections, which in turn passes them to get_observatory. Therefore, even if a TOA in a tim file is already defined in the barycenter, whenever it is read with get_TOAs it will get the clock corrections, unless the user manually sets get_TOAs(..., include_gps=False, include_bipm=False). This would, however, affect all TOAs in the tim file, not only the barycentered ones!

The behavior with TZRMJD when TZRSITE is BAT follows from the same origin. While building the TZR TOA, get_TZR_toa(self,toas) in absolute_phase.py reads the clock correction information from the provided toas, and calls get_TOAs_array() explicitly, passing the include_gps and include_bipm parameters that were applied on the toas. Since these were set to True when reading them for the first time, it ends up applying the correction also on TZRMJD, even if it is defined in BAT..

At least, this was my attempt to follow the code around. It involves an interplay of get_observatory() [observatory/init.py], get_TOAs(), get_TOAs_array(), apply_clock_corrections.py [toa.py], and get_TZR_toa.py [absolute_phase.py].

[gcerib]

Hi @abhisrkckl ,

Unfortunately I am on a technical shift to the telescopes in La Palma, I don't have much time these days to continue digging on this. Did you have a chance to have a look at the problem?

So far, the only solution that worked consistently is adding a check for the BarycenterObs class in get_observatory() (#1707 (comment)), and it passes the tests I've implemented. but it also has some drawbacks. Namely, if apply_clock_corrections() is called on some TOAs multiple times, and there is a barycentered TOA in the list, all calls after the first one will fail.

Cheers,
Giovanni

[abhisrkckl]

No problem. I'll close this PR and open a separate one.

[gcerib]

@abhisrkckl I think in here self.tz_clck_info should also reflect the fact that if TZRMJD is barycentered, no corrections have been applied. What do you think?

[gcerib]

Sorry, it would be line 124/136, it is just outside the commit and cannot comment directly on it.

[abhisrkckl]

Please see my comment above. I think if we fix the bug in get_observatory, this should get fixed automatically.

[gcerib]

Can you please add a test to verify that your original issue has been resolved?

@abhisrkckl I've added the test you requested, hope it helps!