Bug in xpphase (met0 command line switch ignored)

[lucabaldini]

This was reported by Fei

Dear Luca,

I find a strange thing in ixpeobssim pipeline.xpphase.

For line 112 in /ixpeobssim/bin/xpphase.py,  phase is calculated with time reference of 'TSTART', not the given reference time met0.

I tested them with commands attached (line 112 and line 113), they will result in different peak phase positions, while the profiles are both correct based on the fact that met_to_phase() function in /srcmodel/ephemeris.py is using the met0 correctly.

Is this on purpose or a bug?

If this is on purpose, it may result in a difficult in dealing with the phase lag.

Best wishes,

Fei

[lucabaldini]

Actually, looking further into this, I now remember why this implemented this way in the first place---the first argument to ephemeris.fold() is not intended to be the reference met from the ephemeris, but rather a convenient time origin to improve the numerical accuracy of the calculation, see the docstring

    def fold(self, met, start_met, phi0=0.):
        """Fold the MET values to return the corresponding arrays of pulse phase.

        .. warning::

            This is creating a new xEphemeris object with the reference MET set
            to start_met, and the frequency derivatives changed accordingly, and
            then calling the met_to_phase() method of this new object. This is
            supposed to roundtrip with the rvs() method to a decent accuracy,
            and is therefore the interface used in the xpphase application.
            Note, however, that there is a lot of numerical (more or less random)
            rounding going on under the hood, and this is not intended
            for "professional" use. (In a sense, the function internals are
            tweaked specifically for the numerical noise to play in such a way
            that we're making approximately the same errors in rvs() and fold()).
        """
        met0 = start_met
        nu0 = self.nu(start_met)
        nudot0 = self.nudot(start_met)
        phase = xEphemeris(met0, nu0, nudot0, self.nuddot).met_to_phase(met)
        pulse_phase = numpy.mod(phase + phi0, 1)
        return pulse_phase

Now, the current implementation might very well be wrong, but I would need more context (e.g., a complete setup, including the input files and the desired output) to reproduce the problem and fix this.

@s-silvestri and Fei, can you help me out, here?

[s-silvestri]

Are you sure that this is for accuracy and not just for convenience (all positive numbers)? This could ultimately be the reason why the results are so TSTART-dependent in the main branch.

I can provide some file for testing as soon as i get back to my computer

Il giorno sab 8 ott 2022 alle ore 15:49 Luca Baldini < @.***> ha scritto:

Actually, looking further into this, I now remember why this implemented this way in the first place---the first argument to ephemeris.fold() is not intended to be the reference met from the ephemeris, but rather a convenient time origin to improve the numerical accuracy of the calculation, see the docstring

def fold(self, met, start_met, phi0=0.):
    """Fold the MET values to return the corresponding arrays of pulse phase.

    .. warning::

        This is creating a new xEphemeris object with the reference MET set
        to start_met, and the frequency derivatives changed accordingly, and
        then calling the met_to_phase() method of this new object. This is
        supposed to roundtrip with the rvs() method to a decent accuracy,
        and is therefore the interface used in the xpphase application.
        Note, however, that there is a lot of numerical (more or less random)
        rounding going on under the hood, and this is not intended
        for "professional" use. (In a sense, the function internals are
        tweaked specifically for the numerical noise to play in such a way
        that we're making approximately the same errors in rvs() and fold()).
    """
    met0 = start_met
    nu0 = self.nu(start_met)
    nudot0 = self.nudot(start_met)
    phase = xEphemeris(met0, nu0, nudot0, self.nuddot).met_to_phase(met)
    pulse_phase = numpy.mod(phase + phi0, 1)
    return pulse_phase

Now, the current implementation might very well be wrong, but I would need more context (e.g., a complete setup, including the input files and the desired output) to reproduce the problem and fix this.

@s-silvestri https://github.com/s-silvestri and Fei, can you help me out, here?

— Reply to this email directly, view it on GitHub https://github.com/lucabaldini/ixpeobssim/issues/651#issuecomment-1272324846, or unsubscribe https://github.com/notifications/unsubscribe-auth/AQSBKVI3K3FD47SRLD3KW6TWCF3WTANCNFSM6AAAAAAQ7XUJBE . You are receiving this because you were mentioned.Message ID: @.***>

[lucabaldini]

Are you sure that this is for accuracy and not just for convenience (all positive numbers)?

Well, I am sure that my intent when I coded this was for accuracy :-)

if you look at test_roundtrip in test_ephemeris.py

    def test_roundtrip(self, duration=1.e7, padding=1.e10):
        """
        """
        met0 = self.crab_ephemeris.met0
        for start_met in (met0, met0 + padding, met0 - padding):
            met1 = numpy.linspace(start_met, start_met + duration, 100000)
            phase = self.crab_ephemeris.met_to_phase(met1)
            met2 = self.crab_ephemeris.phase_spline_inverse(start_met, duration)(phase)
            self.assertTrue(numpy.allclose(met2, met1))
            met3 = self.crab_ephemeris.phase_to_met(phase, start_met, duration)
            self.assertTrue(numpy.allclose(met3, met1))

the basic idea that I had on my mind was to avoid operating with exceedingly large numbers (e.g., if met0 is much earlier than the observation), and the reasoning was: start from the given ephemeris, create a new ephemeris referred to the start of the observation, and then use that for the folding.

Again, this might be foolish, but before we change this I'd like to understand exacty which problem we are solving.

[s-silvestri]

the basic idea that I had on my mind was to avoid operating with exceedingly large numbers (e.g., if met0 is much earlier than the observation), and the reasoning was: start from the given ephemeris, create a new ephemeris referred to the start of the observation, and then use that for the folding.

But I don't understand why met0 at TSTART would always give better accuracy with respect to a user defined value which is usually the epoch at the middle of the observation and is already optimized for that, or at least that is what i usually choose to get a symmetric interval for the taylor expansion to be valid at both extrema.

[Jaffexf]

I basically agree with Stefano. Firstly, I don't think that you could improve the accuracy with this method. We have no intention to optimize the ephemeris with our pipeline, so the given ephemeris is usually having the epoch close to the observation time or within the duration. Secondly, this is misleading considering the conventional definition of met0. For now, I have no further requirement on xpphase, I use it only for a quick adding of phase column with the correct ephemeris derived from softwares like 'stingray' or 'hendrics'.

[s-silvestri]

Ok, if we all agree on this, a simple change to xpphase.py

-135 phase = eph.fold(time_, kwargs.get('met0'), kwargs.get('phi0'))
+135 if kwargs.get('met0') is None:
+136    met0 = evt_header['TSTART']
+137    logger.warning('Using TSTART as epoch since no met0 was input')
+138 else:
+139    met0 = kwargs.get('met0')
+140 phase = eph.fold(time_, met0, kwargs.get('phi0'))

fixes the issue of unit tests not running correctly when no met0 is given as input and everything looks fine. In this case the results might not be as accurate as by giving met0, so we send a warning to the user

[lucabaldini]

I am afraid I did not make myself clear :-)

I am sorry to insist, but this is not how the current code is working. I am not trying to improve the ephemeris by realigning the calculation to the start time of the observations, nor to change the met0 of the ephemeris. I am just trying to avoid working with exceedingly large numbers, in which case the term in dt**2 might explode. In a nutshell, the edge case I was trying to handle here was the one where the user was passing a set of ephemeris with met0 being a loooong time in the past. All the current code is intended to preserve the full information (including met0) in the ephemeris.

[And, I gather from your finding, it does not, so there is obviously a bug in there that we need to fix.]

Note that met0 is used in the definition of the ephemeris within the _get_ephemeris() function, and it should not be used outside there. (One possible edge case that I see is if you pass a par file and leave the --met0 command line switch to None, in which case met0 is read from file, and correctly stores in the ephemeris object, but then you are using the irrelevant fact that the value of the command-line argument is None later in the main function body.)

The warning that you are suggesting, in particular, is misleading since the first argument to ephemeris.fold() is not intented to represent the met0 of the ephemeris (which, as the name says, is already contained in the ephemeris) but rather a conventional start_time to realign the calculation, which should not have any effect on the output phase, if not for numerical rounding errors when met0 is too far away in the past from the start of the observation.

[Which, again, is not what's happening, and we do have to understand why.]

Can I ask again a full setup reproducing the issue, so that I can see where things go wrong, exactly?

[raestew]

Hello all, I recently ran into an issue that may be related to this thread. I am working with a timing solution that phase aligns an IXPE observation to a NuSTAR observation so that I can perform phase-resolved spectroscopy on the source. Upon inputting the timing solution into ixpeobssim using xpphase, I encountered an unexpected phase shift in the pulse profiles that resulted in the two observations no longer being aligned.  After making sure the timing solution was indeed correct and updating ixpeobssim to version 30.6.4, I examined xpphase.py and saw that line 112 read:       

phase = eph.fold(time_, evt_header['TSTART'], kwargs.get('phi0'))

I believe this had the effect of overwriting the provided t0 with the start time of the observation, which may not be appropriate for all types of timing analysis.  I replaced line 112 in xpphase.py with:          

met0 = kwargs.get('met0')
phase = eph.fold(time_, met0, kwargs.get('phi0'))

and subsequently obtained the phase alignment I had expected. In case xpphase.py will be updated in future releases I wanted to highlight this issue as this may impact timing analysis results performed by ixpeobssim. Sorry if this is not the appropriate channel to bring up this issue, I will be happy to resubmit the same information via a different channel if that is more beneficial to workflow.