Let atlas test the stepsize

code/Is_it_inherent.py

+"""
+What happens if I increase the number density of points in the GW signal, for the n=1 mode fitting the n=1 data?
+
+--Rayne Liu, 09/08/2020
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import rc
+plt.rcParams.update({'font.size': 19})
+import dynesty
+from dynesty import plotting as dyplot
+import qnm
+import random
+import h5py
+import json
+
+#This cell mimicks the (2, 2, 0) and (2, 2, 1) superposition, using the 0.01 stepsize
+tstep = 1
+ndim = 8
+rootpath = "/Users/RayneLiu/git/rdstackingproject"
+
+t = np.arange(0, 80, tstep)
+w0, tau0, x0, y0 = [0.55578191, 11.74090006, 0.98213669, -4.81250993] 
+#Can get the w and tau from example nb and amplitude and phase from the 1910 paper
+w1, tau1, x1, y1 = [0.54, 3.88312743, 4.29386867, -0.79472571]
+print('The fundamental tone frequency, damping time, amplitude and phase:')
+print(w0, tau0, x0, y0)
+print('The n=1 overtone frequency, damping time, amplitude and phase:')
+print(w1, tau1, x1, y1)
+mockdata = x0*np.exp(1j*y0)*np.exp(-t/(tau0)) * (np.cos(w0*t)-1j*np.sin(w0*t)) + \
+            x1*np.exp(1j*y1)*np.exp(-t/(tau1)) * (np.cos(w1*t)-1j*np.sin(w1*t))
+print('The mock data:')
+figdata = plt.figure(figsize = (12, 8))
+plt.plot(t, mockdata.real, label = r'Real')
+plt.plot(t, mockdata.imag, label = r'Imag')
+plt.legend()
+#plt.show()
+figdata.savefig(rootpath + '/plotsmc/n=1_mockdata.png', format='png', bbox_inches='tight', dpi=300)
+
+def modelmock(theta):
+    """
+    theta: comprised of alpha0, alpha1, beta0, beta1, x0, x1, and y0, y1
+    """ 
+    
+    alpha0, alpha1, beta0, beta1, xvar0, xvar1, yvar0, yvar1 = theta
+    #alpha0, beta0, xvar0, yvar0 = theta
+
+    tauvar0 = tau0*(1+beta0)
+    wvar0 = w0*(1+alpha0)
+    tauvar1 = tau1*(1+beta1)
+    wvar1 = w1*(1+alpha1)
+    ansatz = (xvar0*np.exp(1j*yvar0))*np.exp(-t/tauvar0)*(np.cos(wvar0*t)-1j*np.sin(wvar0*t)) +\
+             (xvar1*np.exp(1j*yvar1))*np.exp(-t/tauvar1)*(np.cos(wvar1*t)-1j*np.sin(wvar1*t))
+    # -1j to agree with SXS convention
+    return ansatz
+
+# LogLikelihood function. It is just a Gaussian loglikelihood based on computing the residuals^2
+def log_likelihood(theta):
+    model_mock = modelmock(theta)
+    
+    return  -np.sum((mockdata.real - model_mock.real)**2+(mockdata.imag - model_mock.imag)**2)
+
+def prior_transform(cube):
+    cube[0] = -0.4 + cube[0]*0.8
+    cube[1] = -0.4 + cube[1]*0.8
+    cube[2] = -1 + cube[2]*3
+    cube[3] = -1 + cube[3]*3
+    cube[4] = 0 + cube[4]*6
+    cube[5] = 0 + cube[5]*6
+    cube[6] = -np.pi + cube[6]*2*np.pi
+    cube[7] = -np.pi + cube[7]*2*np.pi
+    return cube
+
+sampler=dynesty.NestedSampler(log_likelihood, prior_transform, ndim, nlive=1000)
+sampler.run_nested()
+res = sampler.results
+res.samples_u.shape
+
+dim = 2
+paramlabels_a = [r'$\alpha_'+str(i)+'$' for i in range (dim)]
+paramlabels_b = [r'$\beta_'+str(i)+'$' for i in range (dim)]
+paramlabels_x = [r'$x_'+str(i)+'$' for i in range (dim)]
+paramlabels_y = [r'$y_'+str(i)+'$' for i in range (dim)] 
+
+paramlabels = paramlabels_a + paramlabels_b + paramlabels_x + paramlabels_y
+
+print('Our constraints:')
+fg, ax = dyplot.cornerplot(res, color='red', 
+                           show_titles=True, labels = paramlabels,
+                           quantiles=None)
+
+fg.savefig(rootpath + '/plotsmc/n=1_mockfit_tstep=' + str(tstep) +'.png', format='png', bbox_inches='tight', dpi=300)
\ No newline at end of file