The approach taken by Moore and Gair is to use the waveform differences between the NR and PN waveforms to train a Gaussian process. Here I investigate whether it is possible to by-pass the PN approximant, and predict a waveform using a Gaussian process trained only from the NR data.

import gwgpr.nr as nr
import numpy as np
import gwgpr.predict as gwpred
import astropy
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')

For computational simplicity we will work with just the “S-series-v2” waveforms.

cat = nr.NRCatalogue('/scratch/aries/gt_bbh/')
cat_f = nr.NRCatalogue('/scratch/aries/gt_bbh/')
cat.waveforms = cat.waveforms[cat.waveforms['series']=='S-series-v2']
cat.waveforms[['q', 'a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL']] = 15*cat.waveforms[['q','a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL']]/np.array(cat.waveforms[['q', 'a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL']].max())

We can have a quick look at how the training set is distributed through parameter space, which should give us some idea of the range of waveforms the final GP should be able to replicate.

import seaborn as sns
cat.waveforms = cat.waveforms[['q','a1','a2', 'th1L', 'th2L', 'ph1', 'ph2', 'th12', 'thSL', 'thJL']]#.to_pandas()
g = sns.pairplot(cat.waveforms,  dropna=True);
 /home/daniel/.virtualenvs/jupyter/local/lib/python2.7/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))
/home/daniel/.virtualenvs/jupyter/local/lib/python2.7/site-packages/matplotlib/__init__.py:892: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))

png

This graph reveals that at least two of the parameters are redundant: they always have the same value, those are ph2 and th2L. (We could very nearly remove q as a parameter too, but it has a few values which don’t lie within its main range). So we can reduce the number of parameters we train with by two (which is good, it makes the computation easier!), and the distribution will look something like this:

import seaborn as sns
pdcat = cat.waveforms[['q','a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL']]
g = sns.pairplot(pdcat,  dropna=True);
png

We now train the GPR with the waveform data and the parameters of the waveform.

predictor = gwpred.GPWaveform(cat, parameters = ['q', 'a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL'], times=(-50, 50))
#skpredictor = gwpred.SKWaveform(cat, parameters = ['q', 'a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL'], times=(-50, 50))
 0 483
1 483
2 483
3 483
4 483
5 483
6 483
7 483
8 483
9 483
10 483
11 483
12 483
13 483
14 483
15 483
16 483
17 483
18 483
{} failed. Continuing.
{} failed. Continuing.
21 483
22 483
23 483
24 483
25 483
26 483
27 483
28 483
29 483
{} failed. Continuing.
31 483
32 483
33 483
{} failed. Continuing.
35 483
36 483
37 483
{} failed. Continuing.
{} failed. Continuing.
{} failed. Continuing.
41 483
42 483
43 483
44 483
45 483
46 483
47 483
48 483
49 483
50 483
51 483
52 483
53 483
54 483
55 483
56 483
57 483
58 483
59 483
60 483
61 483
62 483
63 483
64 483
65 483
66 483
67 483
68 483
69 483
70 483
71 483
72 483
73 483
74 483
75 483
76 483
77 483
78 483
79 483
80 483
81 483
82 483
83 483
84 483
85 483
Training complete.
log Likelihood: 3.67600768525e+16
0 483
1 483
2 483
3 483
4 483
5 483
6 483
7 483
8 483
9 483
10 483
11 483
12 483
13 483
14 483
15 483
16 483
17 483
18 483
{} failed. Continuing.
{} failed. Continuing.
21 483
22 483
23 483
24 483
25 483
26 483
27 483
28 483
29 483
{} failed. Continuing.
31 483
32 483
33 483
{} failed. Continuing.
35 483
36 483
37 483
{} failed. Continuing.
{} failed. Continuing.
{} failed. Continuing.
41 483
42 483
43 483
44 483
45 483
46 483
47 483
48 483
49 483
50 483
51 483
52 483
53 483
54 483
55 483
56 483
57 483
58 483
59 483
60 483
61 483
62 483
63 483
64 483
65 483
66 483
67 483
68 483
69 483
70 483
71 483
72 483
73 483
74 483
75 483
76 483
77 483
78 483
79 483
80 483
81 483
82 483
83 483
84 483
85 483

---------------------------------------------------------------------------

Exception                                 Traceback (most recent call last)

<ipython-input-38-8e9c1ef5a405> in <module>()
      1 predictor = gwpred.GPWaveform(cat, parameters = ['q', 'a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL'], times=(-50, 50))
----> 2 skpredictor = gwpred.SKWaveform(cat, parameters = ['q', 'a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL'], times=(-50, 50))


/scratch/aries/daniel/repositories/gwgpr/gwgpr/predict.py in __init__(self, catalogue, parameters, times)
    125         self.gp = gp = gaussian_process.GaussianProcess()
    126         #self.gp = gp = george.GP(kernel)
--> 127         self.gp.fit(self.train_x, self.train_y)
    128         print("Training complete.")
    129 


/usr/lib/python2.7/dist-packages/sklearn/gaussian_process/gaussian_process.pyc in fit(self, X, y)
    307         if (np.min(np.sum(D, axis=1)) == 0.
    308                 and self.corr != correlation.pure_nugget):
--> 309             raise Exception("Multiple input features cannot have the same"
    310                             " value.")
    311 


Exception: Multiple input features cannot have the same value.

We can now input a set of parameters and see the resulting waveform.

predictor.optimise()
 3.67600768525e+16

---------------------------------------------------------------------------

MemoryError                               Traceback (most recent call last)

<ipython-input-39-7d6cee686a71> in <module>()
----> 1 predictor.optimise()


/scratch/aries/daniel/repositories/gwgpr/gwgpr/predict.py in optimise(self)
     69         # Run the optimization routine.
     70         p0 = self.gp.kernel.vector
---> 71         results = op.minimize(self.nll, p0, jac=self.grad_nll)
     72 
     73         self.gp.kernel[:] = results.x


/usr/lib/python2.7/dist-packages/scipy/optimize/_minimize.pyc in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    370         return _minimize_cg(fun, x0, args, jac, callback, **options)
    371     elif meth == 'bfgs':
--> 372         return _minimize_bfgs(fun, x0, args, jac, callback, **options)
    373     elif meth == 'newton-cg':
    374         return _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,


/usr/lib/python2.7/dist-packages/scipy/optimize/optimize.pyc in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, **unknown_options)
    830     else:
    831         grad_calls, myfprime = wrap_function(fprime, args)
--> 832     gfk = myfprime(x0)
    833     k = 0
    834     N = len(x0)


/usr/lib/python2.7/dist-packages/scipy/optimize/optimize.pyc in function_wrapper(*wrapper_args)
    279     def function_wrapper(*wrapper_args):
    280         ncalls[0] += 1
--> 281         return function(*(wrapper_args + args))
    282 
    283     return ncalls, function_wrapper


/scratch/aries/daniel/repositories/gwgpr/gwgpr/predict.py in grad_nll(self, p)
     59         # Update the kernel parameters and compute the likelihood.
     60         self.gp.kernel[:] = p
---> 61         return -self.gp.grad_lnlikelihood(self.train_y, quiet=True)
     62 
     63     def optimise(self):


/home/daniel/.virtualenvs/jupyter/local/lib/python2.7/site-packages/george/gp.pyc in grad_lnlikelihood(self, y, quiet)
    261 
    262         # Calculate the gradient.
--> 263         A = np.outer(self._alpha, self._alpha) - K_inv
    264         g = 0.5 * np.einsum('ijk,ij', Kg, A)
    265 


MemoryError: 
a = 4
x,y,z = cat.load(a, ['q', 'a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL']).output_full()
times = np.linspace(-50, 50, 500)
pars = x[2][:-1] #/ np.array(cat.waveforms[['q', 'a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL']].max())
pars[4] += 3.9
py, ps = predictor.predict(pars, times)
#ps /= ps.sum()
py /= py.sum()
x = np.array(x)
f, ax2 = plt.subplots(1, sharex=True, figsize=(15,6))
f.suptitle('GPR Predicted waveform')
std = np.sqrt(np.diag(ps/ps.sum()))
ax2.fill(np.concatenate([times, times[::-1]]),
         np.concatenate([py - 1.9600 * std, (py + 1.9600 * std)[::-1]]),
         alpha=.3);
ax2.set_xlim([-50,50]);
ax2.plot(times, py)
ax2.set_xlabel('Time [sec since maximum]');
ax2.set_ylabel('h+ Strain');
f.subplots_adjust(hspace=0.02);
plt.setp([a.get_xticklabels() for a in f.axes[:-1]], visible=False);
plt.tight_layout();
png

We can also explore around a given point in the parameter space and have a look at the uncertainties in the prediction.

from matplotlib import cm

plots = ['q', 'a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL']
f, ax = plt.subplots(len(plots), len(plots), figsize=(13,13))


#plt.imshow(yerr, extent=[0, 1, 1, 0], aspect=1)
#plt.colorbar(label="Standard deviation")
#plt.plot(np.array(cat.waveforms[['a1', 'a2']]), '+')
#plt.plot(np.array(cat.waveforms[['a1']]),np.array(cat.waveforms[['a2']]), '.')

res = 10
yerr = np.zeros((res,res))    

#ax[1,0].contour(x1, x2,yerr);
    
f.subplots_adjust(hspace=0.01, wspace=0.01);
plt.setp([a.get_xticklabels() for a in f.axes[::]], visible=False);
plt.setp([a.get_yticklabels() for a in f.axes[::]], visible=False);

for i,val in enumerate(plots):
    ax[-1,i].set_xlabel(val);
    plt.setp(ax[-1,i].get_xticklabels(), visible=True, rotation='vertical');
    ax[i, 0].set_ylabel(val);
    plt.setp(ax[i, 0].get_yticklabels(), visible=True)

levels = np.linspace(0, 1, 100)

#pars = [0.5, 0.5, 60, 90, 90, 10, 12.5]
# x[0][:-1]
#diffs = diffs**2
#diffs = diffs.sum(axis=1)
#diffs = np.sqrt(diffs)    

wv = cat.waveforms[['q','a1','a2', 'th1L', 'ph1', 'th12', 'thSL', 'thJL']]  
diffs = np.array(wv / wv.max()) - pars/np.array(wv.max())



for i in xrange(len(plots)):
    for j in xrange(len(plots)):
        
        pars = [ 14.99565126,  11.25      ,  11.25      ,   2.5       ,
        15.        ,   5.        ,   5.23255814,   7.64367816]
        diffs = np.array(wv / wv.max()) - pars/np.array(wv.max())
        if j<i: 
            ax[j,i].axis('off')
            continue
        elif j==i:
            #ax[j,i].patch.set_visible(False)
            #ax[j,i].plot(np.array(cat.waveforms[[plots[i]]]), np.sqrt((diffs**2).sum(axis=1)), '.')
            ax[j,i].hist2d(np.squeeze(np.array(cat.waveforms[[plots[i]]])), np.sqrt((diffs**2).sum(axis=1)), bins=20, cmap=cm.get_cmap('Greys'));
            plt.setp(ax[j, i].get_yticklabels(), visible=True)
            ax[j,i].yaxis.tick_right()
            ax[j,i].grid(b=False);
            continue
            
        lowlim1 = np.float(cat.waveforms[[ plots[i] ]].min())
        lowlim2 = np.float(cat.waveforms[[ plots[j] ]].min())
        higlim1 = np.float(cat.waveforms[[ plots[i] ]].max())
        higlim2 = np.float(cat.waveforms[[ plots[j] ]].max())
        #lowlim1, higlim1 = pars[i]-2, pars[i]+2
        #lowlim2, higlim2 = pars[j]-2, pars[j]+2
        
        x1, x2 = np.meshgrid(np.linspace(lowlim1, higlim1, res),
                             np.linspace(lowlim2, higlim2, res))
                             #np.linspace(0, 0.4, res))


       
        times = np.linspace(-10,10, 30)
        


        for l in range(res):
            for m in range(res):
                pars = [ 14.99565126,  11.25      ,  11.25      ,   2.5       ,
        15.        ,   5.        ,   5.23255814,   7.64367816]
                pars[i]=x1[l,m]
                pars[j]=x2[l,m]
                py, ps = predictor.predict(pars, times)
                yerr[l,m] = np.mean(np.sqrt(np.diag(ps)))

        im=ax[j,i].contourf(x1, x2,yerr, levels=levels, cmap=cm.get_cmap('cubehelix'))#cmap=cm.get_cmap('Oranges', len(levels) - 1));
        #ax[j,i].plot(np.array(cat.waveforms[[ plots[i] ]]),np.array(cat.waveforms[[ plots[j] ]]), '.', alpha=0.1, label="Training points")
        
cbar_ax = f.add_axes([0.87, 0.3, 0.05, 0.5])
f.colorbar(im, cax=cbar_ax, label="Standard deviation");


f.text(0.5, 0.95, 'Exploring the area arround \n \
                    $q$ \t {0[0]} \n \
                    $a_1$ \t {0[1]} \n \
                    $a_2$ \t {0[2]} \n \
                    $\\theta_$ \t {0[3]} \n \
                    $\phi_1$ \t {0[4]} \n \
                    $\\theta_1$\t {0[5]} \n \
                    $\\theta_$ \t {0[6]} \n \
                    $\\theta_$ \t {0[7]} \n '.format(pars),
        verticalalignment='top', horizontalalignment='left',
        color='black', fontsize=15)

plt.legend()
plt.tight_layout();
#f.savefig('100_resultscorner.png', dpi=300)
png