Integration in presence of Linear Control Parameters

[brownadder]

When the control is a linear combination of other functions (as discussed with Guannan)

$f(t) = \sum_j \theta_j \phi_j(t)$

Then the $\theta_j$ parameters completely specify the pulse $f(t)$. Everytime the value of $\theta_j$ changes, we do not want to compute the integrals of $f(t)$ again. Instead, note:

$\mu_0^{f}(t) = \int_0^h f(t+s) ds = \sum_j \theta_j \int_0^h \phi_j(t+s) ds = \sum_j \theta_j \mu_0^{\phi_j}(t)$

So all we have to do is compute and store $\phi_j$ integrals and linearly combine when needed. Similarly,

$\mu_1^{f}(t) = \int_0^h (s-h/2) f(t+s) ds = \sum_j \theta_j \int_0^h (s-h/2) \phi_j(t+s) ds = \sum_j \theta_j \mu_1^{\phi_j}(t)$

$\nu$ becomes more complicated, but the same idea works:

$\nu^{f,g}(t) = \int_0^h \int_0^\zeta f(t+\zeta) g(t + \xi) = \sum_j \sum_k \theta_j \beta_k \int_0^h \int_0^\zeta \phi_j(t+\zeta) \phi_j(t+\xi) = \sum_j \sum_k \theta_j \beta_k \nu^{\phi_j,\phi_k}(t) $

Lastly, even if $g$ is expressed in different basis, it doesn't change much:

$\nu^{f,g}(t) = \int_0^h \int_0^\zeta f(t+\zeta) g(t + \xi) = \sum_j \sum_k \theta_j \beta_k \int_0^h \int_0^\zeta \phi_j(t+\zeta) \psi_j(t+\xi) = \sum_j \sum_k \theta_j \beta_k \nu^{\phi_j,\psi_k}(t) $

We just need to precompute different things. If the second entity is not encoded with such a (linear combination) of basis function, this breaks down and we will have to see if anything can be done or not.

[brownadder]

How do we handle this with an appropriate data structures and design?

1. An abstract class should be created:

Class LCF % linearly controlled function - is this a good name? probably not

init should take the list of parameters, its length is N which defines the dimension (j=1,...,N)

def phi(j, t) should return value of jth basis function at time t. It should be aware of N due to init.

def val(t) should return the value i.e. $\sum_j \theta_j \phi(j,t)$

2. Concrete classes can then be derived from it, e.g. FourierLCF or WaveletLCF etc. This will only require specifying phi(j,t), which is quite simple.

3. In an optimal control function, a precompute function should be created which should do the required precomputations. It should figure out the various basis used and should pre-compute and store their integals mu0,mu1,nu in a dictionary. These should give the required $\mu_0^{\phi_j}(t_n)$ i.e. for each $j$ and $n$, and similarly $\mu_1$. $\nu$ needs more data $\nu^{\phi_j, \psi_k}(t_n)$, i.e. three dimensional data for each pair of combinations of function types. Fourier-Fourier, Wavelet-Wavelet, Fourier-Wavelet etc. Note that $N$ also matters, so Fourier10-Fourier100 would be separate.

4. The integration methods mu should be upgraded in this way:

If instead of a function handle, you get an instance of an LCF, e.g. f an instances of Fourier with N=10. Then find precomputed values of Fourier for N=10, and find the required $\mu_0^{\phi_j}(t_n)$. Use this to compute the integral (linear combination).

This only holds if the integral range agrees with the precomputed time grid, but this should be a simple matter of checking that timegrid in the mu0 function is the same as that used in precomputed dataset before doing a lookup.

4. Integration method nu

In this case, it is also possible that f has a Fourier decomposition and g has Wavelet. In that case we need to find pre-computed data for Fourier-Wavelet and read appropriate nu from there.

5. Instead of doing pre-compute, we can do lazy evaluation, cacheing etc. to only compute when needed, but store/cache when computed.

[brownadder]

• First implementation should be Fourier only
• The parameters of the LCF need to be stored as Torch tensors with requires gradient on by default so we can autodiff
• If you pass f.val as a function to mu or nu, it should behave just like a normal function and we should get the correct integrals as we do for other functions. This will also be a unit test for your fancier approach of passing f as an object to mu or nu! Also, while you fix your latter approach, you can test out the rest of the apparatus of your numerical scheme using f.val instead of f.