Phase rounding feels inconsistent around 0.5

[theXYZT]

See this code snippet:

In [1]: x = np.arange(10) + 0.5

In [2]: x
Out[2]: array([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5])

In [3]: ph = Phase(x)

In [4]: ph.int
Out[4]: <Angle [ 0.,  2.,  2.,  4.,  4.,  6.,  6.,  8.,  8., 10.] cycle>

In [5]: ph.frac
Out[5]: <FractionalPhase [ 0.5, -0.5,  0.5, -0.5,  0.5, -0.5,  0.5, -0.5,  0.5, -0.5] cycle>

This is an edge case that feels inconsistent, if not incorrect, to me. The 0.5 values should reliably be part of a specific "integer cycle". However, since np.round rounds to the even integer, we get even integer cycles having both -0.5 and +0.5 entries, and the odd integers have none.

For simply folding pulse profiles we only care about fractional phases, so we can fix this via:

In [6]: FractionalPhase(ph)
Out[6]: <FractionalPhase [-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5] cycle>

which solves half the problem because FractionalPhase uses a range of [-0.5, 0.5) consistently while wrapping values.

However, it doesn't fix the integer part being inconsistently distributed. So, if you are trying to fold pulse stacks or individual pulses or doing anything where you require consistent management of whole cycles and well-distributed fractional cycles within them, then this can be a problem.

I've fixed this by editing the implementation of astropy.time.utils.day_frac to use np.floor(x + 0.5) instead of np.round(x) which consistently rounds towards +inf. This gives more consistent results (that are much more desirable in the use-case that the Phase class was written for - pulsar folding).

Now, I get:

In [1]: x = np.arange(10) + 0.5

In [2]: x
Out[2]: array([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5])

In [3]: ph = Phase(x)

In [4]: ph.int
Out[4]: <Angle [ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.] cycle>

In [5]: ph.frac
Out[5]: <FractionalPhase [-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5] cycle>

This seems consistent to me and I don't think I have introduced any additional errors by doing this. What do you think?

Unless baseband-tasks specifically needs the rounding behavior of np.round, I think this would be a good change.

[mhvk]

My overall feeling would be that there should be a separate method to get the fraction, which includes arguments for the lower edge, units, and rounding method. Could be somewhat like the .fold() method of astropy.timeseries: https://docs.astropy.org/en/latest/api/astropy.timeseries.TimeSeries.html#astropy.timeseries.TimeSeries.fold.