Skip to content
Snippets Groups Projects
Commit b6f9cb22 authored by Rayne Liu's avatar Rayne Liu
Browse files

The test case

parent 33a3c90f
No related branches found
No related tags found
No related merge requests found
"""
In this script I use a mock data generated from combining n=0 and n=1 overtones together, and sample it with the ptemcee sampler under the n=1 scenario. It seems that either there are some fundamental internal problems that prevent us from separating alpha_0 from alpha_1, or that ptemcee doesn't work as well for this problem. Either way, there should be something going on with our approach for this entire project...
--- R^a_{yne} L^i_u
"""
import numpy as np
import corner
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from matplotlib import rc
plt.rcParams['font.family'] = 'DejaVu Sans'
rc('text', usetex=True)
plt.rcParams.update({'font.size': 19})
import pandas as pd
import ptemcee
import qnm
import random
ntemps = 20
nwalkers = 100
npoints = 150
ndim = 8
burnin = 75
numbins = 42
datacolor = '#105670' #'#4fa3a7'
pkcolor = '#f2c977' #'#ffb45f'
mediancolor = '#f7695c' #'#9b2814'
rootpath = "/Users/RayneLiu"
t = np.arange(0, 80, 0.1)
w0, tau0, x0, y0 = [0.55578191, 11.74090006, 0.978518, -2.11289]
#Can get the w and tau from example nb and amplitude and phase from the 1910 paper
w1, tau1, x1, y1 = [0.54351639, 3.88312743, 4.29435, 1.38519]
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))
def modelmock(theta):
"""
theta: comprised of alpha0, alpha1, beta0, beta1, x0, x1, and y0, y1
"""
alpha0, alpha1, beta0, beta1, xvar0, xvar1, yvar0, yvar1 = 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
def log_prior(theta):
alpha0, alpha1, beta0, beta1, xvar0, xvar1, yvar0, yvar1 = theta
if all([-0.4 <= alpha0 <= 0.4, -1.0 <= beta0 <= 2.0, 0 <= xvar0 <= 2.4, -2*np.pi <= yvar0 <= 0, \
-0.4 <= alpha1 <= 0.4, -1.0 <= beta1 <= 2.0, 0 <= xvar1 <= 6.0, -np.pi <= yvar1 <= np.pi]):
return 0.0
return -np.inf
# 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)
# Logposterior distribution for the residuals case.
# The evidence is just a normalization factor
def log_probability(theta):
lp = log_prior(theta)
if not np.isfinite(lp):
return -np.inf
return lp + log_likelihood(theta)
#Define labels used in plotting
paramlabels_a = [r'$\alpha_0$', r'$\alpha_1$']
paramlabels_b = [r'$\beta_0$', r'$\beta_1$']
paramlabels_x = [r'$x_0$', r'$x_1$']
paramlabels_y = [r'$y_0$', r'$y_1$']
paramlabels = paramlabels_a + paramlabels_b + paramlabels_x + paramlabels_y
np.random.seed(42)
pos = [random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), random.uniform(-0.1,0.1), \
random.uniform(0.8, 1.0), random.uniform(4, 5), random.uniform(-3, -2), random.uniform(1, 2)]
pos += 1e-5 * np.random.randn(ntemps, nwalkers, ndim)
sampler = ptemcee.Sampler(nwalkers, ndim, log_likelihood, log_prior, ntemps=ntemps)
sampler.run_mcmc(pos,npoints);
#Chain plot
fig, axes = plt.subplots(ndim, 1, sharex=True, figsize=(12, 4*(4)))
for i in range(ndim):
axes[i].plot(sampler.chain[0,:, :, i].T, color="k", alpha=0.4, rasterized=True)
axes[i].yaxis.set_major_locator(MaxNLocator(5))
axes[i].set_ylabel(paramlabels[i])
axes[-1].set_xlabel('Iterations')
fig.savefig(rootpath+'/git/rdstackingproject/plotsmc/Test_with_mockdatan=1'+ str(npoints)+'pts_chainplot.png', format = 'png', bbox_inches = 'tight')
#Output chain
samples = sampler.chain[0,:, burnin:, :].reshape((-1, ndim))
df0 = pd.DataFrame(samples, columns=paramlabels)
df0.to_csv(rootpath+'/git/rdstackingproject/plotsmc/Test_with_mockdatan=1'+ str(npoints)+'pts_chain.csv', index = False)
#Corner plot
lglk = np.array([log_likelihood(samples[i]) for i in range(len(samples))])
pk = samples[np.argmax(lglk)]
median = np.median(samples, axis=0)
figcorn = corner.corner(samples, bins = numbins, hist_bin_factor = 5, color = datacolor, truths=pk, truth_color = pkcolor, plot_contours = True, labels = paramlabels, quantiles=(0.05, 0.16, 0.5, 0.84, 0.95), levels=[1-np.exp(-0.5), 1-np.exp(-1.64 ** 2/2)], show_titles=True)
#Extract the axes in order to add more important line plots
naxes = len(pk)
axes = np.array(figcorn.axes).reshape((naxes, naxes))
# Loop over the diagonal
for i in range(naxes):
ax = axes[i, i]
ax.axvline(median[i], color=mediancolor)
# Loop over the histograms
for yi in range(naxes):
for xi in range(yi):
ax = axes[yi, xi]
ax.axvline(median[xi], color=mediancolor)
ax.axhline(median[yi], color=mediancolor)
ax.plot(median[xi], median[yi], color = mediancolor, marker = 's')
figcorn.savefig(rootpath+'/git/rdstackingproject/plotsmc/Test_with_mockdatan=1'+ str(npoints)+'pts_corner.pdf', format = 'pdf')
\ No newline at end of file
universe = vanilla
getenv = true
# run script -- make sure that condor has execute permission for this file (chmod a+x script.py)
executable = A_test_with_n=1.py
# file to dump stdout (this directory should exist)
output = A_test_with_n=1.py
# file to dump stderr
error = A_test_with_n=1.py
# condor logs
log = A_test_with_n=1.py
initialdir = .
notify_user = rl746@cornell.edu
notification = Complete
arguments = "-processid $(Process)"
request_memory = 232GB
request_cpus = 16
on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
accounting_group = aei.dev.test_dynesty
queue 1
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment