Skip to content
Snippets Groups Projects
Commit 2b512ed6 authored by Yifan Wang's avatar Yifan Wang
Browse files

add fisher matrix

parents
Branches master
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
curve_data.txt 0 → 100644
+ 0
− 0
View file @ 2b512ed6
Source diff could not be displayed: it is too large. Options to address this: view the blob.
fisher_matrix_numerical.ipynb 0 → 100644
+ 534
− 0
View file @ 2b512ed6
%% Cell type:markdown id: tags:
## Testing Parity with Fisher Matrix
%% Cell type:code id: tags:
``` python
%matplotlib inline
import numpy as np
import sys
import lal
import lalsimulation as lalsim
import pycbc.conversions as conversions
import pycbc.psd
#from lalinference.rapid_pe import lalsimutils
import matplotlib.pyplot as plt
from scipy import interpolate
from numpy import pi
from scipy import constants
import scipy.interpolate
import astropy
from astropy import cosmology
```
%% Cell type:code id: tags:
``` python
import matplotlib as mpl
# PLOTTING OPTIONS
fig_width_pt = 3*246.0 # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27 # Convert pt to inch
golden_mean = (np.sqrt(5)-1.0)/2.0 # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt # width in inches
fig_height = fig_width*golden_mean # height in inches
fig_size = [fig_width,fig_height]
params = { 'axes.labelsize': 24,
'font.family': 'serif',
'font.serif': 'Computer Modern Raman',
'font.size': 24,
'legend.fontsize': 20,
'xtick.labelsize': 24,
'ytick.labelsize': 24,
'axes.grid' : True,
'text.usetex': True,
'savefig.dpi' : 100,
'lines.markersize' : 14,
'figure.figsize': fig_size}
mpl.rcParams.update(params)
```
%% Cell type:code id: tags:
``` python
def generate_waveform(m_chirp, eta, dist, Aeff,delta_f, flow, fhigh,
approx=lalsim.IMRPhenomPv2, both=False):
"""
Generate a waveform given a chirp mass and symmetric mass ratio
for binaries with zero spin.
m_chirp -- chirp mass in solar masses
eta -- float symmetric mass ratio
dist -- distance in Mpc
Aeff -- beyond GR parameters
delta_f -- frequency step size in Hz
flow -- minimum frequency in Hz
fmax -- maximum frequency in Hz
approx -- waveform approximant (e.g. lalsimulation.IMRPhenomPv2)
both -- Boolean dictating the summation of h_plus and h_cross
Returns:
h1 -- magnitude of h_plus and h_cross in the frequency domain
"""
# Generate a waveform following the same procedure in waveform_overlap.ipynb
m1 = conversions.mass1_from_mchirp_eta(m_chirp, eta)
m2 = conversions.mass2_from_mchirp_eta(m_chirp, eta)
m1 *= lal.MSUN_SI
m2 *= lal.MSUN_SI
s1x, s1y, s1z = 0.0, 0.0, 0.0
s2x, s2y, s2z = 0.0, 0.0, 0.0
coa_phase, polarization_angle, eccentricity = 0.0, 0.0, 0.0
dist *= lal.PC_SI * 1e6
incl = 0.0
params = lal.CreateDict()
lal.DictInsertINT4Value(params,"parity_beta_mu",1)
lal.DictInsertREAL8Value(params,"parity_Aeff_phi",Aeff)
hpf, hxf = lalsim.SimInspiralChooseFDWaveform(
m1, m2,
s1x, s1y, s1z,
s2x, s2y, s2z,
dist,incl,coa_phase,
polarization_angle,eccentricity,0.0,
delta_f,
flow, fhigh, 20.0,
params,
approx)
#print hpf.f0, hpf.deltaF, len(hpf.data.data)
# For convenience, we'll reuse hpf and redefine to be h1
if both:
hpf.data.data += hxf.data.data
#print 'Combining hpf and hxf data'
h1 = hpf
#TODO: cf Eric Thrane's script to include beam pattern function
return h1.data.data
```
%% Cell type:code id: tags:
``` python
psd_file = np.loadtxt('curve_data.txt')
psd_interpolation = scipy.interpolate.interp1d(psd_file[:,0],psd_file[:,2]**2)
```
%% Cell type:markdown id: tags:
# Plot the waveform with parity violation
%% Cell type:code id: tags:
``` python
flow = 20 # Hz
fhigh = 2048 # Hz
df = 1./32
h1 = generate_waveform(30,0.25,400,0.,df,flow,fhigh)
h2 = generate_waveform(30,0.25,400,1e11,df,flow,fhigh)
h1 = h1[int(flow/df):]
h2 = h2[int(flow/df):]
freqs = np.arange(flow, fhigh + df, df)
plt.figure()
plt.plot(freqs,h1,label='$A_\mathrm{eff}=0$')
plt.plot(freqs,h2,label='$A_\mathrm{eff}=10^{11}$')
plt.xlim(20,400)
plt.legend()
plt.show()
```
%% Output
/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order)
findfont: Font family ['serif'] not found. Falling back to DejaVu Sans.
%% Cell type:markdown id: tags:
# Define the derivative functions
%% Cell type:code id: tags:
``` python
def deriv_chirp_mass(m_chirp, eta, distance, Aeff, df, flow, fhigh):
"""
Find the partial derivative of h(f) with respect to chirp mass.
m_chirp -- chirp mass in solar masses
eta -- symmetric mass ratio
redshift -- redshift value
df -- frequency step size in Hz
flow -- minimum frequency in Hz
fhigh -- maximum frequency in Hz
Returns:
Array of partial derivative of h(f) with respect to chirp mass for various frequencies
"""
derivs = []
# Calculate the waveform h(f) at this given parameter set
waveform_orig = generate_waveform(m_chirp, eta, distance, Aeff, df, flow, fhigh)
# Find the h(f) for a chirp mass very close to the original chirp mass
delta_mc = 0.001 # Solar Masses
waveform_changed = generate_waveform(m_chirp + delta_mc, eta, distance, Aeff, df, flow, fhigh)
# Clip off h(f) values corresponding to frequencies before flow -- !!! CHECK THIS
flow_index = int(flow/df)
waveform_orig = waveform_orig[flow_index:]
waveform_changed = waveform_changed[flow_index:]
# I'm concerned that I'm actually cutting off high values, not low values
# Calculate the derivative for each frequency value
for index, value in enumerate(waveform_orig):
derivs.append((waveform_changed[index] - value) / delta_mc) # !!! - CHECK THIS: is the division ok?
return np.asarray(derivs)
def deriv_eta(m_chirp, eta, distance, Aeff, df, flow, fhigh):
derivs = []
# Calculate the waveform h(f) at this given parameter set
waveform_orig = generate_waveform(m_chirp, eta, distance, Aeff, df, flow, fhigh)
# Find the h(f) for a eta very close to the original eta
delta_eta = 0.001 # Solar Masses
waveform_changed = generate_waveform(m_chirp, eta + delta_eta, distance, Aeff, df, flow, fhigh)
# Clip off h(f) values corresponding to frequencies before flow -- !!! CHECK THIS
flow_index = int(flow/df)
waveform_orig = waveform_orig[flow_index:]
waveform_changed = waveform_changed[flow_index:]
# I'm concerned that I'm actually cutting off high values, not low values
# Calculate the derivative for each frequency value
for index, value in enumerate(waveform_orig):
derivs.append((waveform_changed[index] - value) / delta_eta) # !!! - CHECK THIS: is the division ok?
return np.asarray(derivs)
def deriv_dist(m_chirp, eta, distance, Aeff, df, flow, fhigh):
derivs = []
# Calculate the waveform h(f) at this given parameter set
waveform_orig = generate_waveform(m_chirp, eta, distance, Aeff,df, flow, fhigh)
# Find the h(f) for a eta very close to the original eta
delta_dist = 1.
waveform_changed = generate_waveform(m_chirp, eta , distance+delta_dist, Aeff,df, flow, fhigh)
# Clip off h(f) values corresponding to frequencies before flow -- !!! CHECK THIS
flow_index = int(flow/df)
waveform_orig = waveform_orig[flow_index:]
waveform_changed = waveform_changed[flow_index:]
# I'm concerned that I'm actually cutting off high values, not low values
# Calculate the derivative for each frequency value
for index, value in enumerate(waveform_orig):
derivs.append((waveform_changed[index] - value) / delta_dist) # !!! - CHECK THIS: is the division ok?
return np.asarray(derivs)
def deriv_Aeff(m_chirp, eta, distance, Aeff, df, flow, fhigh):
derivs = []
# Calculate the waveform h(f) at this given parameter set
waveform_orig = generate_waveform(m_chirp, eta, distance, Aeff,df, flow, fhigh)
# Find the h(f) for an Aeff very close to the original Aeff
delta_Aeff = 0.1
waveform_changed = generate_waveform(m_chirp, eta, distance, Aeff+delta_Aeff, df, flow, fhigh)
# Clip off h(f) values corresponding to frequencies before flow -- !!! CHECK THIS
flow_index = int(flow/df)
waveform_orig = waveform_orig[flow_index:]
waveform_changed = waveform_changed[flow_index:]
# I'm concerned that I'm actually cutting off high values, not low values
# Calculate the derivative for each frequency value
for index, value in enumerate(waveform_orig):
derivs.append((waveform_changed[index] - value) / delta_Aeff) # !!! - CHECK THIS: is the division ok?
return np.asarray(derivs)
```
%% Cell type:markdown id: tags:
# Fisher Information Matrix
%% Cell type:code id: tags:
``` python
def find_Aeff_cov(m_chirp, eta, distance, Aeff, df, flow, fhigh):
"""
Calculate a Fisher Matrix given a set of masses, redshifts, and frequency values.
We assume a minimum frequency value of 30 Hz and a maximum of 2048 Hz.
m_chirp -- chirp mass in solar mass
eta -- symmetric mass ratio
redshift -- redshift value
df -- spacing between discretely sampled frequency values
Returns:
FM -- 4-by-4 matrix for fisher information matrix using chirp mass, eta, distance and Aeff derivatives
Aeff's variance
"""
# Set frequency range
freqs = np.arange(flow, fhigh + df, df)
# Find PSD at given frequency values
# map is for python2
#psd = map(lalsim.SimNoisePSDaLIGOZeroDetHighPower, freqs)
#psd = pycbc.psd.EinsteinTelescopeP1600143(len(freqs),df, flow)
psd = psd_interpolation(freqs)
# Calculate the partial of h(f) with respect to chirp mass, eta, dchi2 and dist
dMc = deriv_chirp_mass(m_chirp, eta, distance, Aeff,df, flow, fhigh)
dEta = deriv_eta(m_chirp, eta, distance, Aeff,df, flow, fhigh)
dDist = deriv_dist(m_chirp, eta, distance, Aeff,df, flow, fhigh)
dAeff = deriv_Aeff(m_chirp, eta, distance, Aeff,df, flow, fhigh)
FisherList = [dMc, dEta, dDist, dAeff]
FisherMatrix = np.zeros(shape=(4,4))
for i in range(4):
for j in range(4):
FisherMatrix[i,j] = 4*np.real(np.sum(FisherList[i]*np.conjugate(FisherList[j]) / psd * df))
return np.sqrt(np.linalg.inv(FisherMatrix)[3][3])
def find_SNR(m_chirp, eta, distance, Aeff, df, flow, fhigh):
# Set frequency range
freqs = np.arange(flow, fhigh + df, df)
# Find PSD at given frequency values
psd = psd_interpolation(freqs)
hp = generate_waveform(m_chirp, eta, distance, Aeff,df, flow, fhigh)
flow_index = int(flow/df)
hp = hp[flow_index:]
return np.sqrt(4*np.real(np.sum(hp*np.conjugate(hp) / psd * df)))
#TODO: investigate whether the selected signals' SNR passes the threshold
```
%% Cell type:code id: tags:
``` python
flow = 10 # Hz
fhigh = 2048 # Hz
df = 1./32
chirpmass = 30
eta = 0.249
distance = 2000
Aeff = 0
Aeff_constraint = find_Aeff_cov(chirpmass,eta,distance, Aeff, df,flow,fhigh)
eventSNR = find_SNR(chirpmass,eta,distance, Aeff, df,flow,fhigh)
Aeff_constraint,eventSNR
```
%% Output
(540827913.424093, 265.3088993879092)
%% Cell type:markdown id: tags:
## Convert Aeff to Mpv
%% Cell type:code id: tags:
``` python
from pycbc import cosmology
from scipy import integrate
QE = 1.60217656535e-19
def integrand_parityaeff(redshift, parity_beta):
"""
The integrand:
(1.0 + z)^parity_beta / sqrt(Omega_m (1+z)^3 + Omega_Lambda)
"""
omega_m = 0.3075 #pycbc.cosmology.get_cosmology().Om0 # matter density
omega_l = 0.6910098821161554 #pycbc.cosmology.get_cosmology().Ode0 # dark energy density
return (1.0+redshift)**(parity_beta)/ np.sqrt(omega_m*(1.0+redshift)**3.0 + omega_l)
def Aeff_to_Mpv(parity_beta, parity_aeff, redshiftval):
"""
convert Aeff to Mpv
"""
zs = np.zeros(parity_aeff.shape, dtype=float) # the output array
for (ii, val) in enumerate(parity_aeff):
zs[ii] = integrate.quad(integrand_parityaeff, 0, redshiftval[ii] ,args=(parity_beta))[0]
return (parity_aeff / zs ) ** (1.0 / parity_beta) * QE
```
%% Cell type:code id: tags:
``` python
redshift = cosmology.redshift(distance)
mpv = Aeff_to_Mpv(1,np.array([Aeff_constraint]),np.array([redshift]))
1./mpv/1e9
```
%% Output
array([4.50256769])
%% Cell type:markdown id: tags:
### Therefore, for a 30 chirpmass BH located at 2000 Mpc, Mpv can be constrained to be 4.5 GeV
%% Cell type:code id: tags:
``` python
```
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint