mcmc fits

code/RDGW150914.py

+#!/usr/bin/env python
+# coding: utf-8
+
+# In[379]:
+
+
+import numpy as np
+import corner
+import matplotlib.pyplot as plt
+import emcee
+import math
+import h5py
+import pandas as pd
+import json
+import qnm
+import random
+from scipy.optimize import minimize
+
+
+# In[410]:
+
+
+rootpath="/work/francisco.jimenez/sio/git/rdstackingproject"
+npoints=1000000
+nmax=1
+
+
+# In[411]:
+
+
+gw = {}
+gw["SXS:BBH:0305"] = h5py.File(rootpath+"/SXS/BBH_SKS_d14.3_q1.22_sA_0_0_0.330_sB_0_0_-0.440/Lev6/rhOverM_Asymptotic_GeometricUnits_CoM.h5", 'r')
+gw_sxs_bbh_0305 = gw["SXS:BBH:0305"]["Extrapolated_N2.dir"]["Y_l2_m2.dat"]
+
+
+# In[412]:
+
+
+metadata = {}
+with open(rootpath+"/SXS/BBH_SKS_d14.3_q1.22_sA_0_0_0.330_sB_0_0_-0.440/Lev6/metadata.json") as file:
+    metadata["SXS:BBH:0305"] = json.load(file)
+
+af = metadata["SXS:BBH:0305"]['remnant_dimensionless_spin'][-1]
+mf = metadata["SXS:BBH:0305"]['remnant_mass']
+
+
+# In[413]:
+
+
+times = gw_sxs_bbh_0305[:,0]
+
+
+# In[414]:
+
+
+position = np.argmax(times >= (3702.75))
+gw_sxs_bbh_0305rd=gw_sxs_bbh_0305[position:-1]
+timesrd=gw_sxs_bbh_0305[position:-1][:,0]
+
+
+# In[415]:
+
+
+def EradUIB2017(eta,chi1,chi2):
+    
+    m1=0.5*(1+(1-4*eta)**0.5)
+    m2=0.5*(1-(1-4*eta)**0.5)
+    S= (m1**2*chi1 + m2**2*chi2)/(m1*m1 + m2*m2)
+    
+    erad=(((1-(2*(2)**0.5)/3)*eta+0.5609904135313374*eta**2-0.84667563764404*eta**3+3.145145224278187*eta**4)*(1+S**3 *(-0.6320191645391563+ 4.952698546796005*eta-10.023747993978121*eta**2)+S**2*(-0.17762802148331427+ 2.176667900182948*eta**2)+S*(-0.13084389181783257- 1.1387311580238488*eta+5.49074464410971*eta**2)))/(1+S*(-0.9919475346968611+ 0.367620218664352*eta+4.274567337924067*eta**2))-0.01978238971523653*S*(1-4.91667749015812*eta)*(1-4*eta)**0.5 *eta *(chi1-chi2)-0.09803730445895877*(1-4*eta)**0.5*(1-3.2283713377939134*eta)*eta**2 *(chi1-chi2)+0.01118530335431078*eta**3 *(chi1-chi2)**2
+    return  erad
+
+
+# In[416]:
+
+
+omegas = []
+for i in range (0,nmax+1):
+    grav_220 = qnm.modes_cache(s=-2,l=2,m=2,n=i)
+    omega = grav_220(a=af)[0]
+    omegas.append(omega)
+
+
+# In[417]:
+
+
+#deprecated
+def model(theta, t, y):
+    x0, y0, a0, b0 = theta
+    #x0, y0= theta
+    w = (np.real(omegas))[0]/mf
+    tau=-(np.imag(omegas))[0]*mf
+
+    # -1j to agree with SXS convention
+    return (x0+y0*1j)*(np.exp(-(t-timesrd[0])/(tau*(1+b0))))*(np.cos((1+a0)*w*t)-1j*np.sin((1+a0)*w*t))
+
+
+# In[418]:
+
+
+def model_dv(theta, t, y):
+    #x0, y0= theta
+    w = (np.real(omegas))/mf
+    tau=-(np.imag(omegas))*mf
+    dim =int(len(theta)/4)
+    
+    xvars = []
+    yvars = []
+    avars = []
+    bvars = []
+    for i in range (0,dim):
+        xvar = theta[4*i]
+        xvars.append(xvar)
+        yvar = theta[4*i+1]
+        yvars.append(yvar)
+        avar = theta[4*i+2]
+        avars.append(avar)
+        bvar = theta[4*i+3]  
+        bvars.append(bvar)
+        
+    
+    ansatz = 0
+    for i in range (0,dim):
+        ansatz = ansatz + (xvars[i]+1j*yvars[i])*np.exp(-(t-timesrd[0])/(tau[i]*(1+bvars[i])))        *(np.cos((1+avars[i])*w[i]*t)-1j*np.sin((1+avars[i])*w[i]*t))
+    # -1j to agree with SXS convention
+    return ansatz
+
+
+# In[419]:
+
+
+def log_likelihood(theta, t, y):
+    modelev = model_dv(theta,t,y)
+    return  -np.sum((y[:,1] - (modelev.real))**2+(y[:,2] - (modelev.imag))**2)
+
+
+# In[420]:
+
+
+def log_likelihood_match(theta, t, y):
+    
+    #model and data
+    modelev = model_dv(theta,t,y)
+    data = y[:,1] + 1j*y[:,2]
+    
+    #norms
+    norm1=np.sum(modelev*np.conj(modelev))
+    norm2=np.sum(data*np.conj(data))
+
+    #mismatch
+    myTable = data*np.conj(modelev);    
+    return -(1-(np.sum(myTable)).real/np.sqrt(norm1*norm2)).real
+
+
+# In[421]:
+
+
+def log_prior(theta):
+    x_s = theta[0::4]
+    y_s = theta[1::4]
+    a_s = theta[2::4]
+    b_s = theta[3::4]
+    if all(-10 <= t <= 10 for t in x_s) and all(-10 <= t <= 10 for t in y_s) and all(-0.1 <= t <= 0.1 for t in a_s) and all(-0.1 <= t <= 0.1 for t in b_s):
+        return 0.0
+    return -np.inf
+
+
+# In[422]:
+
+
+def log_probability(theta, t, y):
+    lp = log_prior(theta)
+    if not np.isfinite(lp):
+        return -np.inf
+    return lp + log_likelihood(theta, t, y)
+
+
+# In[423]:
+
+
+def log_probability_match(theta, t, y):
+    lp = log_prior(theta)
+    if not np.isfinite(lp):
+        return -np.inf
+    return lp + log_likelihood_match(theta, t, y)
+
+
+# Maximum estimator Fitting
+
+# In[424]:
+
+
+np.random.seed(42)
+nll = lambda *args: -log_likelihood(*args)
+initial = np.array([-3,4,0.03,-0.1])
+soln = minimize(nll, initial, args=(timesrd, gw_sxs_bbh_0305rd))
+x0_ml, y0_ml, a0_ml, b0_ml = soln.x
+print("Maximum likelihood estimates:")
+print("x0 = {0:.3f}".format(x0_ml))
+print("y0 = {0:.3f}".format(y0_ml))
+print("a0 = {0:.3f}".format(a0_ml))
+print("b0 = {0:.3f}".format(b0_ml))
+
+
+# In[359]:
+
+
+# with mismatch function
+np.random.seed(42)
+nll = lambda *args: -log_likelihood_match(*args)
+initial = np.array([-3,4,0.03,-0.1])
+soln = minimize(nll, initial, args=(timesrd, gw_sxs_bbh_0305rd))
+x0_ml, y0_ml, a0_ml, b0_ml = soln.x
+print("Maximum likelihood estimates:")
+print("x0 = {0:.3f}".format(x0_ml))
+print("y0 = {0:.3f}".format(y0_ml))
+print("a0 = {0:.3f}".format(a0_ml))
+print("b0 = {0:.3f}".format(b0_ml))
+
+
+# In[425]:
+
+
+plt.plot(timesrd, gw_sxs_bbh_0305rd[:,2], "r", alpha=0.3, lw=3, label="NR_im")
+modelfit = model([x0_ml,y0_ml,a0_ml,b0_ml],timesrd,gw_sxs_bbh_0305rd)
+plt.plot(timesrd, modelfit.imag,"b", alpha=0.3, lw=3, label="Fit_im")
+#plt.plot(x0, np.dot(np.vander(x0, 2), w), "--k", label="LS")
+plt.legend(fontsize=14)
+plt.xlim(timesrd[0], timesrd[0]+80)
+plt.xlabel("t")
+plt.ylabel("h");
+
+
+# mcmc Fitting
+
+# In[430]:
+
+
+nwalkers = 32
+ndim = int(4*(nmax+1))
+
+pos = [random.uniform(-10,10) ,random.uniform(-10,10) ,random.uniform(-0.1,0.1) ,random.uniform(-0.1,0.1)]
+for i in range (1,nmax+1):
+       pos_aux = [random.uniform(-10,10) ,random.uniform(-10,10) ,random.uniform(-0.1,0.1) ,random.uniform(-0.1,0.1)]
+       pos = pos + pos_aux
+
+pos = pos + 1e-1 * np.random.randn(nwalkers, ndim)
+
+
+# In[431]:
+
+
+sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability, args=(timesrd, gw_sxs_bbh_0305rd))
+sampler.run_mcmc(pos,npoints, progress=True);
+
+
+# In[432]:
+
+
+label = ['x0','y0','a0','b0']
+for i in range (1,nmax+1):
+    label2 = ['x' + str(i),'y' + str(i),'a' + str(i),'b' + str(i)]
+    label = label + label2
+
+
+# In[433]:
+
+
+flat_samples = sampler.get_chain(discard=100, thin=15, flat=True)
+median=np.median(sampler.flatchain, axis=0)
+
+fig = corner.corner(
+    flat_samples, labels=label, truths=median,quantiles=(0.1, 0.9)
+);
+
+
+# In[436]:
+
+
+fig.savefig(rootpath+'/plots/fit_chi2_'+str(nmax)+'.pdf')
+
+
+# In[716]:
+
+
+#sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability_match, args=(timesrd, gw_sxs_bbh_0305rd))
+#sampler.run_mcmc(pos,npoints, progress=True);
+
+
+# In[ ]:
+
+
+#labels = ["x0", "y0","a0","b0"]
+#flat_samples = sampler.get_chain(discard=100, thin=15, flat=True)
+#fig = corner.corner(
+#    flat_samples, labels=labels, truths=[x0_ml, y0_ml,a0_ml,b0_ml]
+#);
+
+
+# In[ ]:
+
+
+#fig.savefig(rootpath+'/plots/fit_match.pdf')
+

code/condor_submit_RdownPy.sub

+universe   = vanilla
+getenv     = true
+# run script -- make sure that condor has execute permission for this file (chmod a+x script.py)
+executable = RDGW150914.py
+# file to dump stdout (this directory should exist)
+output     = RDGW150914.out
+# file to dump stderr
+error     = RDGW150914.err
+# condor logs
+log     = RDGW150914.log
+initialdir = .
+notify_user = frjifo@aei.mpg.de
+notification = Complete
+arguments  = "-processid $(Process)" 
+request_memory = 16GB
+request_cpus = 1
+on_exit_remove = (ExitBySignal == False) || ((ExitBySignal == True) && (ExitSignal != 11))
+accounting_group = aei.dev.test_dynesty
+queue 1
+