```python
import numpy as np

solar_mass_in_seconds = 4.92686088e-6

def phase(f, Mc):

        phase = -np.pi/4. + ( 3./( 128. * pow(Mc*np.pi*f, 5./3.) ) )
        return phase

def htilde(f, Mc):
        Mc *= solar_mass_in_seconds
        htilde = pow(f, -7./6.) * pow(Mc,5./6.) * np.exp(1j*phase(f,Mc))
        return htilde
```

In this note I will show how we use the empirical interpolant of a function to make a "reduced order quadrature" (ROQ). The exercise is to compute the inner product between a data set d (which plays the role of $A(x,y)$ in the coupling constant work) and the function h. The standard inner product is:
\begin{equation}
\langle\, d | h(M_c)\, \rangle = \Re \Delta f \sum_{i=1}^{N}\,d^{*}(f_{min} + i \Delta f)\,h(M_c;f_{min} +i \Delta f)\,,
\end{equation}
In this quadrature rule there are N terms. In general N is the number of sample points of the discretely sampled data and function h and could be very large for some problems. If we substitute the empirical interpolant of h into the above equation for the inner product then we get the ROQ. To recap, the empirical interpolant of h is
\begin{equation} 
\mathcal{I}[h](M_c;f) = \sum_{j=1}^m B_j (f) h(M_c;F_j)
\end{equation}
and in the previous note we computed $\{B_j\}_{i=1}^{n}$ and $\{F_j\}_{i=1}^{n}$. Substituting the empirical interpolant for h in the inner product we get 
\begin{equation}
\langle\,d | \mathcal{I}[h]\,\rangle = \Re \Delta f\,\sum_{j=1}^{m}\sum_{i=1}^{N} d^{*}(f_{min} + i \Delta f) B_j(f_{min} + i \Delta f) h(M_c;F_j)\,.
\end{equation}
Notice that the sum over frequency now decouples from the dependence on $M_c$! Thus we can do all the sums over frequency once-and-for-all, and evaluating the inner product should be cheaper. Thus we write the reduced order quadrature for the inner product as
\begin{equation}
\langle\, d | h(M_c)\, \rangle_{\mathrm{ROQ}} = \Re \sum_{j=1}^{m}\,w_j\,h(M_c;F_j)\,,
\end{equation}




where

\begin{equation}
w_j = \Delta f\sum_{i=1}^{N} d^{*}(f_{min} + i \Delta f) B_j(f_{min} + i \Delta f)\,.
\end{equation}

Notice that the ROQ has only m terms in the sum, whereas the original had N. Thus we get a compression of N/m. For this specific problem, the original waveforms which we built the bases for had N=5000 samples. The number of bases is around m=300, so in this case we should expect a compression of N/m = 5000/300 = 17. Let's try it out.

First I'll create a data set. It doesn't matter what it is, so I'll just let it be a random vector of length 5000. I'll start by loading in the frequency series and frequency steps $\Delta f$, and then generate the data:


```python
fseries = np.loadtxt("fseries.dat")
df = np.loadtxt("df.dat")

data = np.random.rand(len(fseries)) + 1j*np.random.rand(len(fseries))
```

Next I'll load in the $B_j$s and the empirical interpolation nodes $F_j$.


```python
B_re = np.loadtxt("B_re.dat")
B_im = np.loadtxt("B_im.dat")

B = B_re + 1j*B_im

nodes = np.loadtxt("nodes.dat")
```

Now I have everything I need to make the weights $\omega_j$ which is done below


```python
# multipy data by \Delta f
data *= df

weights = np.inner(B, data.conjugate().T)
```

Now I'll compute the regular inner product and the ROQ


```python
d_dot_h = np.vdot(data, htilde(fseries, 1.5)).real # regular inner product

ROQ = np.dot(weights, htilde(nodes, 1.5)).real # ROQ inner product

print "regular inner product = %.15e"%d_dot_h
print "ROQ = %.15e"%ROQ
```

    regular inner product = 5.805838449180775e-07
    ROQ = 5.805838449179256e-07


so the ROQ is essentially error free at round off. Next I'll time the regular inner product and the ROQ.


```python
import time

Mc = 2

t1 = time.time()
for i in range(50000):
    np.vdot(data, htilde(fseries, Mc)).real # regular inner product
e1 = time.time()

t2 = time.time()
for i in range(50000):
    np.dot(weights, htilde(nodes, Mc)).real # ROQ inner product
e2 = time.time()

print "regular inner product took %f s"%((e1-t1)/50000.)
print "ROQ took %f s"%((e2-t2)/50000.)
print "speedup = %f"%((e1-t1) / (e2-t2))
```

    regular inner product took 0.000582 s
    ROQ took 0.000068 s
    speedup = 8.501153


This seems to give a speed up between 8 and 9, which isn't quite the predicted speed up of 16: presumably there is some overhead in passing the vectors to np.vdot and np.dot which is on the order of the time it takes to compute the inner products. Anyway, hopefully this illustrates the main point which is that if you have an empirical interpolant then computing inner products of the kind we considered can be significantly sped up. This is what we should do for the coupling coefficient work!


```python

```

b509830e58e79e8fcb1bcdbeedeab27302f11c58
ipynb
Jupyter Notebook
ROQ/Building a reduced order quadrature (ROQ).ipynb
mattpitkin/random_scripts
8fcfc1d25d8ca7ef66778b7b30be564962e3add3
MIT
ROQ/Building a reduced order quadrature (ROQ).ipynb
mattpitkin/random_scripts
8fcfc1d25d8ca7ef66778b7b30be564962e3add3
MIT
ROQ/Building a reduced order quadrature (ROQ).ipynb
mattpitkin/random_scripts
8fcfc1d25d8ca7ef66778b7b30be564962e3add3
MIT
Qwen/Qwen-72B
1. YES
2. YES