design matrix uses a fixed F0

[scottransom]

I'm working on updating random models to use the design matrix and while testing that, I noticed small differences with the current "exact" calculation which subtracts residuals.

In timing_model.py, the designmatrix() method, which we use in our fitters, currently converts from phase to time simply by dividing by m.F0.value. For MSPs, this is likely fine, but for young pulsars, as @paulray has mentioned in the past, this can cause fairly large differences.

Anyone know if there is a reason why the design matrix is this way? Is this compatible with Enterprise's use of the design matrix?

I'm happy to fix this (and potentially give the user an option of return the matrix in term of phase or time) if we agree that it is not the behavior that we should have.

[dlakaplan]

I seem to recall a number of issues like this. I feel like I fixed one or two, but I definitely didn't look exhaustively. There's the residuals.get_PSR_freq() routine, which has modelF0 as the default. Could we change that to something else? And I know the same math is implemented in several places (also in spindown.py, and probably elsewhere?).

[aarchiba]

There can be headaches in fitting when uncertainties are converted between phase and time using fitted parameters F0, F1, ... - for one thing our derivatives of chi-squared with respect to these parameters become wrong. Convergence can also get weird - if twiddling a parameter rescales the errors, well, maybe there's a way to make the errors super small. But of course Paul is right that when the spin frequency changes substantially over the observation span we need to be more careful which value we use.

There was a suggestion that a "reference F0" could be in the par file and used for some calculations but not modified during fitting. Possibly this could be extended to include reference values of F1, ... .

There might also be performance implications - an accurate calculation of the pulsar spin frequency at a certain time could be expensive if it's supposed to be the observed frequency - do we include the Doppler shifts from orbital motion?

I do think it makes sense to have a single function for converting phase to frequency, that accepts options for several possibilities (use a reference F0, use full observed spin frequency, ...). Then at least we could easily audit the places it was used and discuss carefully what the right option was for each of those cases.

[paulray]

There was some earlier discussion here: #866

[dlakaplan]

PhaseJump also converts to time using a fixed F0. That might be another one to change.

[aarchiba]

Also some discussion in #713 where it was also suggested that all such conversions go through a single function (get_PSR_freq?) and in #415

[dlakaplan]

But I think that function should be moved to timing_model rather than residuals

[aarchiba]

Indeed, what you need for a get_psr_freq-like function all lives in a TimingModel, in any of its versions. It's not quite clear how we're going to pick the right version, but if all the phase/time conversions go through it at least we can work towards sensible answers.

[scottransom]

Just to let people know, I'm going to move and improve get_PSR_freq() (and change its name) along with fixing this design matrix issue with my random model improvements PR. They are all linked. Hopefully you will see a PR over the next couple days.