Commit abe7dfc0 authored by Rayne Liu's avatar Rayne Liu
Browse files
parents 0238221e 072a4c73
This diff is collapsed.
This diff is collapsed.
......@@ -4144,7 +4144,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.10"
"version": "3.7.3"
}
},
"nbformat": 4,
#!/usr/bin/env python
# coding: utf-8
"""
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 = 0.01
ndim = 8
rootpath = '/work/rayne.liu/git/rdstackingproject'#"/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)
#!/usr/bin/env python
# coding: utf-8
"""
What happens if I fit the mock data simply using omega and tau, and don't look at alpha and beta?
--Rayne Liu, 09/09/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 = 0.01
ndim = 8
rootpath = '/work/rayne.liu/git/rdstackingproject'#"/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 wvar0, wvar1, tauvar0, tauvar1, xvar0, xvar1, and yvar0, yvar1
"""
wvar0, wvar1, tauvar0, tauvar1, xvar0, xvar1, yvar0, yvar1 = theta
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] = (1.-0.4)*w0 + cube[0]*0.8*w0
cube[1] = (1.-0.4)*w1 + cube[1]*0.8*w1
cube[2] = cube[2]*3*tau0
cube[3] = cube[3]*3*tau1
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=2000)
sampler.run_nested()
res = sampler.results
res.samples_u.shape
dim = 2
paramlabels_w = [r'$\omega_'+str(i)+'$' for i in range (dim)]
paramlabels_t = [r'$\tau_'+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_w + paramlabels_t + paramlabels_x + paramlabels_y
print('Our constraints:')
fg, ax = dyplot.cornerplot(res, color='blue',
show_titles=True, labels = paramlabels,
quantiles=None)
fg.savefig(rootpath + '/plotsmc/n=1_mockfit_tstep=' + str(tstep) +'wandt.png', format='png', bbox_inches='tight', dpi=300)
This diff is collapsed.
universe = vanilla
getenv = true
# run script -- make sure that condor has execute permission for this file (chmod a+x script.py)
executable = Is_it_inherent.py
# file to dump stdout (this directory should exist)
output = Is_it_inherent.out
# file to dump stderr
error = Is_it_inherent.err
# condor logs
log = Is_it_inherent.log
initialdir = .
notify_user = rl746@cornell.edu
notification = Complete
arguments = "-processid $(Process)"
request_memory = 232GB
request_cpus = 1
on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
accounting_group = aei.dev.test_dynesty
queue 1
universe = vanilla
getenv = true
# run script -- make sure that condor has execute permission for this file (chmod a+x script.py)
executable = Is_it_inherent_wandt.py
# file to dump stdout (this directory should exist)
output = Is_it_inherent_wandt.out
# file to dump stderr
error = Is_it_inherent_wandt.err
# condor logs
log = Is_it_inherent_wandt.log
initialdir = .
notify_user = rl746@cornell.edu
notification = Complete
arguments = "-processid $(Process)"
request_memory = 64GB
request_cpus = 1
on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
accounting_group = aei.dev.test_dynesty
queue 1
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment