Added mul and rmul to Phase.py and created test file for Phase.py

[champagne-cmd]

No description provided.

[paulray]

This is great! Thanks Chloe!

It has also made me think, and I've been looking at how Phase() is used in our codes and I have a couple of thoughts on the Phase class. Interested in @demorest, @scottransom, and @luojing1211's thoughts:

• Interestingly, the Residuals.phase_resids does not return Phase objects! Instead, it does this:

    # This converts full from a Phase object to a np.float128
    full = full.int + full.frac
    return full

Of course, residuals are not phases, but phase differences. This reminds me of Time and TimeDelta in astropy. I wonder if we should consider something similar, or alternatively just make a Phase method that returns a float128 version of a Phase, instead of having other code do full.int + full.frac?

• The other thing commonly done is to get phases in the range [0,1) instead of [-0.5,0.5). This is a simple thing, but again maybe should be done with a Phase method instead of by the user of a Phase. Maybe phase.frac_nonnegative (or some better name) could return np.where(self.frac < 0.0, self.frac + 1.0, self.frac)? This is commonly used for histogramming phases especially when the phases are photons.
• There are a few places where the code accesses ph.frac as ph[1]. This works because Phase is a tuple and the fractional part is the second element in the tuple, but it is not clear and readable. I see that @mhvk fixed some of those in his PR New phase class #695. We should do the same here. See mcmc_fitter.py lines 640 and 643; test_fermiphase.py; test_absphase.py and there may be others.
• It is interesting that Phase is not used by the model classes that compute the model phase. For example see spindown_phase(). It computes dt as a longdouble then evaluates the Taylor series and returns a longdouble Quantity. It is only converted to a Phase object in TimingModel before being summed with any other phase terms. In other places, jump_phase() just uses np.zeros(), so it is probably just a float64, but glitch_phase does use np.longdouble. I'm not sure if this is any big deal especially since at least glitch phases and jump phases are relative and thus always quite small in magnitude, but I wonder if we are really gaining any precision by using Phase? Something to think about.
• I asked @champagne-cmd to add the multiply capabilities to Phase, but now I'm wondering if that is really a good idea. I think for sure multiplying a Phase by a scalar should work, but should we support multiplying two Phase objects? Oops, my mistake. I see this isn't there.

[mhvk]

@champagne-cmd, @paulray - seeing this, I think you should definitely consider just using my phase class in #695 - it has all this already included (I now wish I had pushed it up earlier, as it obviously would have saved work).

[luojing1211]

One of the major problems is how do we handle the tow-double type. the reason we use np.float128 is that we don't have a good two double computation support (i.e. APIs and speed). If we can have a two double type in NumPy, like the complex numbers, it could be a lot easier. @mhvk Is there any two-double data type in numpy?

[mhvk]

No, there isn't anything in numpy - I think judicious using of int/frac for Time and Phase should cover it, though.

[luojing1211]

It would be useful to have a standard two-double type in Numpy. But it is out of our scope. I think if we have the mul, rmul, etc functions setup, we can compute from a two-double time object directly to a two-double phase. One thing I am not sure of is speed. Evaluating a Taylor series using two-double could be time-consuming.

[mhvk]

Agreed. Note that #695 does do multiplication and division, ensuring that precision is preserved (like TimeDelta). Plus, as mentioned there, trig functions that just use the fraction, etc.

[paulray]

I agree @luojing1211 . We can experiment with that. One thing I think @champagne-cmd will work on is a benchmarking code that will help us evaluate the performance effects of such changes.

[scottransom]

@paulray and @champagne-cmd, that's great to hear about the benchmarking. Are there plans to do some full-fledged profiling of the code as well (for pulsars with small, med, and large numbers of TOAs)? I think that would be incredibly valuable to help clean up some low-hanging fruit optimization-wise.

The other thing to remember for the residuals class, for instance, is that while we are computing phase differences, those residuals used by an external optimization program (almost always demanding double precision) for the fits. And we aren't always dealing with just fractional parts of phase for this as we need to handle PHASE terms and known pulse rotation numbers that mean that the phase differences might be huge (which is exactly why there is that
full = full.int + full.frac
thing.

[paulray]

@scottransom we don't have full fledged plans yet, so suggestions on what to do for performance testing are most welcome! I recall a previous discussion of how Quantity and Units can be done in ways that are more or less efficient. Is there a reference for that somewhere that has best practices for performance?

And I agree that residuals are typically required to be normal double precision numbers. So, I think that having some code in Phase to return the difference of two Phase objects as either float128 or float64 in as efficient a way as possible might be helpful.

[luojing1211]

Astropy unit tips https://docs.astropy.org/en/stable/units/ in the Performance Tips section. The biggest thing is to avoid making copies. By the way, there is a simple script to do the profiling in the PINT profiling directory. It may be not the best but help me quite a bit.

[scottransom]

ooh... that's a useful page. I had completely forgotten about the << operator.....