Should we organize a short meeting to discuss this and close this issue?

**Jean-Baptiste Bayle**
(017d59b1)
*
at
30 Nov 12:06
*

Fix report job

**Jean-Baptiste Bayle**
(6d3742ab)
*
at
29 Nov 14:25
*

Fix report job

**Jean-Baptiste Bayle**
(cfb039e9)
*
at
28 Nov 17:39
*

Add switch to toggle concurrency

So, if I try to summarize. We need to compute

`\delta \hat\tau_i^{t}(\tau) = \delta \tau_i^{t}(\tau) + \delta \hat\tau_i^{\tau_i}(\tau + \delta \tau_i^{t}(\tau)).`

from the known time series `\delta \tau_i^{t}(\tau)`

(from the orbit file) and `\delta \hat\tau_i^{\tau_i}(\tau)`

(from the instrumental simulation). The latter is defined as the sum of a stochastic clock noise `N^q_i(\tau)`

and deterministic drifts `\{y_{i,0}, y_{i,1}, y_{i,2}\}`

,

`\delta \hat{\tau_i}^{\tau_i}(\tau) = y_{i,0} \cdot \tau + \frac{y_{i,1}}{2} \cdot \tau^2 + \frac{y_{i,2}}{3} \cdot \tau^3 + N^q_i(\tau),`

so that we need to compute

`\delta \hat\tau_i^{t}(\tau) = \delta \tau_i^{t}(\tau) + y_{i,0}(\tau + \delta \tau_i^{t}(\tau)) + \frac{y_{i,1}}{2}(\tau + \delta \tau_i^{t}(\tau))^2 + \frac{y_{i,2}}{3}(\tau + \delta \tau_i^{t}(\tau))^3 + N^q_i(\tau + \delta \tau_i^{t}(\tau)).`

We know how to compute each of these terms, since the first one is given in the orbit files, the deterministic drifts can be analytically evaluated, and the last can be numerically interpolated (with splines, as explained). Interpolation is the only part that is (moderately) computationally expensive and (badly) memory intensive.

So the game is to **assess if we can get rid of this interpolation term (and maybe some of the deterministic parts) and still match the requirements in terms of precision**. So the questions are

- What are the requirements in terms of timing precision here?
- What are the precision of various approximations?

If I understand correctly what you did @Nik, you have computed the full expression (see above), some approximations following expansions of the full expression at various orders in `\delta \tau_i^{t}(\tau)`

. However, I'm not sure I understand the expansion.

After 10 years of mission, `\delta \tau_i^{t}(\tau)`

is of the order of several seconds (see https://gitlab.in2p3.fr/lisa-simulation/orbits/-/jobs/576897/artifacts/file/esa-leading-orbits-sc1.pdf), so it's not particularly small. However, we can probably take a look at the various expansion orders. If I didn't make any mistakes,

- 0-th order expansion is
`[\delta \hat\tau_i^{t}]^{(0)}(\tau) = y_{i,0} \tau + \frac{y_{i,1}}{2} \tau^2 + \frac{y_{i,2}}{3} \tau^3 + N^q_i(\tau)`

, - 1-st order coefficient is
`[\delta \hat\tau_i^{t}]^{(1)}(\tau) = \delta \tau_i^{t}(\tau) [ 1 + y_{i,0} + \frac{y_{i,1}}{2} \tau + \frac{y_{i,2}}{3} \tau^2 + \dot N^q_i(\tau) ]`

, - 2-nd order coefficient is
`[\delta \hat\tau_i^{t}]^{(2)}(\tau) = \frac{1}{2} \delta (\tau_i^{t}(\tau))^2 [ y_{i,1} + 2 \tau y_{i,2} + \ddot N^q_i(\tau) ]`

.

We don't want to evaluate those clock noise derivatives, so we can use the following "approximated" expansion terms,

- 0-th order approximated expansion is
`[\delta \hat\tau_i^{t}]^{(0)}(\tau) = y_{i,0} \tau + \frac{y_{i,1}}{2} \tau^2 + \frac{y_{i,2}}{3} \tau^3 + N^q_i(\tau)`

, - 1-st order approximated coefficient is
`[\delta \hat\tau_i^{t}]^{(1)}(\tau) \approx \delta \tau_i^{t}(\tau) [ 1 + y_{i,0} + \frac{y_{i,1}}{2} \tau + \frac{y_{i,2}}{3} \tau^2 ]`

, - 2-nd order approximated coefficient is
`[\delta \hat\tau_i^{t}]^{(2)}(\tau) \approx \frac{1}{2} \delta (\tau_i^{t}(\tau))^2 [ y_{i,1} + 2 \tau y_{i,2} ]`

.

I'm not sure these expansions are actually super relevant since they don't capture increasing computing complexity (even if we exclude the stochastic term derivatives, or treat it separately), so I would propose another set of simplifications,

- Frequency offset
`y_{i,0}(\tau + \delta \tau_i^{t}(\tau))`

is replaced by`y_{i,0} \tau`

, - Frequency linear drift
`\frac{y_{i,1}}{2}(\tau + \delta \tau_i^{t}(\tau))^2`

is replaced by`\frac{y_{i,1}}{2} \tau^2`

, - Frequency quadratic drift
`\frac{y_{i,2}}{3}(\tau + \delta \tau_i^{t}(\tau))^3`

is replaced by`\frac{y_{i,2}}{3}\tau^3`

, - Stochastic term
`N^q_i(\tau + \delta \tau_i^{t}(\tau))`

is replaced by`N^q_i(\tau)`

.

We can evaluate the timing error that each of these simplifications introduces by computing, e.g. for the frequency offset term, `|y_{i,0} \tau - y_{i,0}(\tau + \delta \tau_i^{t}(\tau))|`

. Similarly to @Nik, I've used the ESA trailing orbit file to evaluate `\delta\tau_1^{t}(\tau)`

.

I've used the following approximated default values for the frequency offset and drift coefficients,

`y_{i,0} \approx 10^{-7}\,\text{s}^{-1}`

`y_{i,1} \approx 10^{-14}\,\text{s}^{-2}`

`y_{i,2} \approx 10^{-22}\,\text{s}^{-3}`

Here's what I get. This is close (but not exactly what @Nik finds).

Just my opinion: Given that Martin actually only got marginal wall-time improvement, but payed for it with 244% cpu utilization.. I would maybe keep it off by default, but announce it in the usual channels (slack, telecon). We could also ask people to test it on their machines after it's on master, and maybe change it to the default later down the line if the majority sees a decent speed-up.

Alright, so giving the decent speedup we get, I would merge this. Users can still recover the previous exact behavior using `Instrument(concurrent=False)`

, so I would expect that this does not introduce any bugs.

Maybe the only question is: should this be enabled by default?

I re-ran the same test script, using either latest master or this branch. Here is the result:

Master: python3 test-performance.py 25.06s user 6.10s system 108% cpu 28.768 total New: python3 test-performance.py 28.04s user 7.58s system 217% cpu 16.399 total

So quite similar to JBs results, a decent speed-up! Which probably makes sense, as we have almost the same machine.

Hi @j2b.bayle ,

I compared the two methods in general:

THE-dev-wrt-TCB via interpolation vs THE-dev-wrt-TCB via approximation (expansion)

The plots show the residuals for the different expansion orders, i.e. they show the accumulated error of the expansion (if I computed the interpolation correctly)

Hi @Nik, thanks for this work!

I computed now the_wrt_tcb using an interpolation of the_wrt_tps via InterpolatedUnivariateSpline order=5 (is this the way you would interpolate this?)

Yes, this is how we would do it (or with some lower order).

I'm not sure what you plotted here: what did you compare the interpolated THE-to-TPS against?

Hi all,

I computed now the_wrt_tcb using an interpolation of the_wrt_tps via InterpolatedUnivariateSpline order=5 (is this the way you would interpolate this?) and compared the result to the different expansion orders shown above:

If I did not do a mistake, the deviation between the two approaches is higher than the contribution of the higher order terms. I guess then the discussion would not about which expansion terms we include but rather about should we do an expansion at all or better go with the interpolation?

Hi all,

for the case we decide to implement the expansion I tried to evaluate the contributions of the different orders:

Short description:

Upper left: TPS deviation wrt TCB simulated with LISA Orbits using an ESA orbit file

Upper right: THE deviation wrt TPS neglecting stochastic clock noise (did not want to simulate that for 9 years), i.e.

`\delta \hat{\tau}^{\tau}(\tau) = y_0 \cdot \tau + \frac{y_1}{2} \cdot \tau^2 + \frac{y_2}{3} \cdot \tau^3,`

using y_0 = 1e-7, y_1 = 1e-14 s-1, y_2 = 1e-23 s-2.

Center left: THE deviation wrt TCB 0'th order expansion

`\delta \hat{\tau}^{t}(\tau) \approx \delta \tau^{t}(\tau) + \delta \hat{\tau}^{\tau}(\tau)`

Center right: 1'st order correction

`y_{0} \cdot \delta \tau^{t}(\tau)`

Lower left: 2'nd order correction

`\frac{y_{1}}{2} \cdot (\delta \tau^{t}(\tau))^2`

Lower right: 3'rd order correction

`\frac{y_{2}}{3} \cdot (\delta \tau^{t}(\tau))^3`

I guess what still needs to be evaluated is the error we are making when using the expansion instead of the interpolation.

Hi all!

Did anyone put some thoughts into those last 3 questions?

Since @aurelien.hees is adding the first elements for orbit reconstruction and time correlation in the L0-L2 pipeline, we should solve this as soon as possible.

We can easily have two names for the same dataset (underlying data is not copied, it's just a view on the same array). So we can go for `moc_time_correlation`

and `the_wrt_tps`

(and co).

About the model, if there's an easy approximation that fulfills our requirements, and we can easily show that it does, I don't see why we should go for the computationally (more) expansive solution. So what we should check is

- What kind of requirements do we have?
- What are the relevant components that must be taken into account when computing
`\delta \hat{\tau}^{\tau_i}_i(\tau + \delta \tau^t(\tau))`

over 10 years of data? - If one 1 or 2 deterministic terms are sufficient, I'd go for the easy approximation. If not, we can use the full interpolated model.

About the model: what do you think of trying the full model (using a low order interpolation) first, and check what the performance impact is? It should be pretty simple, and if it's negligable, I don't see a reason to use a simplified model.

The expansion is an interesting idea, but would break if either the stochastic part of `\delta \hat \tau_i^\tau`

or the total value of `\delta \tau_i^t`

ever become too large. for the stochastic part, we could probably avoid that by using a band-limited noise. For the `\delta \tau_i^t`

, the main way I can think of this might happen is if we at some point re-define the time TPS-TCB offset at the time origin of the orbit file, which is arbitrarily set to 0 right now.

About the names: `moc_time_correlation`

has the advantage of being very descriptive of what the measurement represents. Having this as output in the h5 file would not be bad, to make it clear this is representing a real measurement.

Though internally, I agree with @j2b.bayle that having a consistent naming for the different time deviations would be nice. The `timer_deviations`

should then probably also be renamed to `the_wrt_tps`

or `the_asfunc_tps`

?

Indeed, that's a conversation we had a long time ago.

@Nik not sure why I used this $y_1$ term, when the dominant one is $y_0$. Your equation is the right one. Note that we can add higher order terms ($y_1$ or $y_2$) if required – but we then comes back to @aurelien.hees's question about the required accuracy here.

I like `moc_time_correlations`

, or more explicit `the_wrt_tcb`

(or variations such as `the_asfunc_tcb`

). Following this last pattern, `tps_proper_time_deviations`

would be renamed `tps_wrt_tcb`

.