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

upload notebook

parent d0dc1bd6
No related branches found
No related tags found
No related merge requests found
Expand all Show whitespace changes
Inline Side-by-side
notebook/prod_plot_whitendata.ipynb 0 → 100644
+ 523
− 0
View file @ 657e874a
%% Cell type:code id:b2885394 tags:
``` python
from pycbc.waveform import get_fd_waveform
from pycbc.waveform import get_td_waveform
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pycbc
from pycbc.inference import io, models
import lal
from pycbc import conversions
from pycbc import cosmology
from pycbc import transforms
from pycbc import coordinates
from pycbc import waveform
from pycbc.detector import Detector
################################################################################
#
# PLOTTING OPTIONS
#
################################################################################
# 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:14a8fdad tags:
``` python
params['figure.figsize']
```
%% Output
[10.211706102117061, 6.311181454233049]
%% Cell type:code id:438f6f50 tags:
``` python
xphm_path = '/work/yifan.wang/gw190521_mpvinverse/second/wednesdy_xphm_gr_v2data/nongr_run/xphm_morecpu_output/config_files/extract_gw190521_nongr_xphm.hdf'
```
%% Cell type:code id:df2cf8c6 tags:
``` python
gr_path = '/work/yifan.wang/gw190521_mpvinverse/second/wednesdy_xphm_gr_v2data/gr_run/gr_run_output/config_files/extract_prod_gr.hdf'
```
%% Cell type:code id:f17009cc tags:
``` python
def fn_to_waveform(fn):
# read in the raw posterior file and extract the max loglikelihood parameter value
with io.loadfile(fn, 'r') as fp:
data = fp.read_data()
psds = fp.read_psds()
cp = fp.read_config_file()
samples = fp.read_samples(list(fp['samples'].keys()))
idx = samples['loglikelihood'].argmax()
params = {p: samples[p][idx] for p in samples.fieldnames}
# set up the model
model = models.read_from_config(cp, data=data, psds=psds)
model.update(**params)
_ = model.loglikelihood
pol = model.current_stats['maxl_polarization']
#compute the max loglikelihood waveform
wfs = model.waveform_generator.generate(polarization=pol,**model.current_params)
data = {}
whdata = {}
for det in model.data:
d = model.data[det]
data[det] = d.to_timeseries()
d = d.copy()
d *= model.weight[det] / (4*model.psds[det].delta_f)**0.5
whdata[det] = d.to_timeseries()
waveforms = {}
whitened_waveforms = {}
for detname in wfs:
hs = {}
whitend_hs = {}
det = Detector(detname)
fp, fc = det.antenna_pattern(model.current_params['ra'],
model.current_params['dec'],
pol,
model.current_params['tc'])
hp, hc = wfs[detname]
h = fp*hp + fc*hc
h.resize(len(model.weight[detname]))
flow = model.current_params['f_lower']
hs[detname] = pycbc.waveform.fd_to_td(h, left_window=(flow, flow+5))
wh = h.copy()
wh *= model.weight[detname] / (4*model.psds[detname].delta_f)**0.5
whitend_hs[detname] = wh.to_timeseries()
waveforms.update(hs)
whitened_waveforms.update(whitend_hs)
return whdata, whitened_waveforms,model
```
%% Cell type:code id:41a0a746 tags:
``` python
xphm_data_whiten, xphm_waveform_whiten,model = fn_to_waveform(xphm_path)
```
%% Output
WARNING:root:No sampling_params section read from config file
WARNING:root:WARNING: The following args are not being used for waveform generation: baseapprox, srcmchirp, redshift, q, comoving_volume, trigger_time, delta_tc, parity_mpvinverse
%% Cell type:code id:26fa693a tags:
``` python
_, gr_wf,gr_model = fn_to_waveform(gr_path)
```
%% Output
WARNING:root:No sampling_params section read from config file
WARNING:root:WARNING: The following args are not being used for waveform generation: srcmchirp, redshift, q, comoving_volume, trigger_time, delta_tc
%% Cell type:code id:d50f9fd8 tags:
``` python
det = list(gr_wf.keys())
```
%% Cell type:code id:5b8fe3fd tags:
``` python
whdata,whitened_waveforms = xphm_data_whiten,xphm_waveform_whiten
xlim=(-0.1, 0.1)
fig = plt.figure(figsize=(10.211706102117061, 6.311181454233049*1.6))
fig.subplots_adjust(hspace=.2)
ax = fig.add_subplot(311);
wh = whdata[det[0]]
ax.plot(wh.sample_times-model.static_params['trigger_time'], wh, c='k', alpha=0.2, label='whitened data')
#ax.plot(gr_data[det[0]].sample_times-gr_model.static_params['trigger_time'],gr_data[det[0]],c='red',alpha=0.2)
wh = whitened_waveforms[det[0]]
ax.plot(wh.sample_times-model.static_params['trigger_time'], wh,label='birefringence')
ax.plot(gr_wf[det[0]].sample_times-gr_model.static_params['trigger_time'], gr_wf[det[0]],label='GR')
ax.grid(True)
ax.set_xticklabels([])
#ax.set_ylabel('H1')
#ax.text(0.8,0.8,'Hanford')
ax.set_xlim(*xlim)
ax.set_ylim(-100,100)
#ax.text(10,8,'Hanford')
#ax.set_xticks([])
# for minor ticks
#ax.set_xticks([], minor=True)
#ax.axes.get_xaxis().set_ticks([])
#plt.gca().xaxis.set_major_locator(plt.NullLocator())
ax.text(0.054, 70, 'LIGO Hanford', style='italic',fontsize=20,
bbox={'facecolor':'white', 'alpha':0.5, 'pad':5,'edgecolor':'none'})
################
bx = fig.add_subplot(312);
wh = whdata[det[1]]
bx.plot(wh.sample_times-model.static_params['trigger_time'], wh, c='k', alpha=0.2, label='whitened data')
wh = whitened_waveforms[det[1]]
bx.plot(wh.sample_times-model.static_params['trigger_time'], wh,label='birefringence')
bx.plot(gr_wf[det[1]].sample_times-model.static_params['trigger_time'], gr_wf[det[1]],label='GR')
bx.set_ylabel('Whitened Strain')
bx.set_xlim(*xlim)
bx.set_ylim(-100,100)
#bx.legend(loc=[0.66,0.71])
#bx.legend(loc='lower right')
bx.text(0.0515, 70, 'LIGO Livingston', style='italic',fontsize=20,
bbox={'facecolor':'white', 'alpha':0.5, 'pad':0,'edgecolor':'none'})
bx.grid(True)
bx.set_xticklabels([])
################
cx = fig.add_subplot(313);
wh = whdata[det[2]]
cx.plot(wh.sample_times-model.static_params['trigger_time'], wh, c='k', alpha=0.2, label='whitened data')
wh = whitened_waveforms[det[2]]
cx.plot(wh.sample_times-model.static_params['trigger_time'], wh,label='birefringence')
cx.plot(gr_wf[det[2]].sample_times-model.static_params['trigger_time'], gr_wf[det[2]],label='GR')
cx.text(0.07, 70, 'Virgo', style='italic',fontsize=20,
bbox={'facecolor':'white', 'alpha':0.5, 'pad':5,'edgecolor':'none'})
#cx.set_ylabel('V1')
cx.set_xlim(*xlim)
cx.set_ylim(-100,100)
cx.set_xlabel('Time [s]')
cx.legend(loc='lower right', ncol=3,framealpha=0)
fig.savefig('0918whitendata.pdf',bbox_inches='tight')
```
%% Output
%% Cell type:markdown id:eb804396 tags:
# Bayesian evidence
%% Cell type:code id:7c729ffb tags:
``` python
import h5py
grfile = h5py.File(gr_path,'r')
nongrfile = h5py.File(xphm_path)
```
%% Cell type:code id:34978275 tags:
``` python
baygr = grfile.attrs['log_evidence'] - grfile['samples'].attrs['lognl']
```
%% Cell type:code id:49553a12 tags:
``` python
baynongr = nongrfile.attrs['log_evidence'] - nongrfile['samples'].attrs['lognl']
```
%% Cell type:code id:8ccc540f tags:
``` python
baynongr,baygr,baynongr-baygr
```
%% Output
(95.5876921692543, 87.74573551728099, 7.841956651973305)
%% Cell type:markdown id:d7e44f8b tags:
# maxL SNR
%% Cell type:code id:6733a640 tags:
``` python
ii = np.argmax(grfile['samples']['loglikelihood'])
```
%% Cell type:code id:3e26e204 tags:
``` python
conversions.snr_from_loglr(grfile['samples']['loglikelihood'][ii] - grfile['samples'].attrs['lognl'])
```
%% Output
15.289043026035419
%% Cell type:code id:05e15a1f tags:
``` python
ii = np.argmax(nongrfile['samples']['loglikelihood'])
```
%% Cell type:code id:59575cf1 tags:
``` python
conversions.snr_from_loglr(nongrfile['samples']['loglikelihood'][ii] - nongrfile['samples'].attrs['lognl'])
```
%% Output
15.473953014212633
%% Cell type:markdown id:975880a9 tags:
# GR maxL parameters
%% Cell type:code id:7869a554 tags:
``` python
gr_model.current_params
```
%% Output
{'spin2_polar': 2.51825060232281,
'chi_p': 0.7386566141564779,
'inclination': 1.0490370243369231,
'srcmchirp': 54.44781356503804,
'srcmass1': 148.33013483168034,
'delta_tc': 0.003219054716161532,
'distance': 2037.5553540826304,
'srcmass2': 29.66232709467983,
'dec': 0.5304731815876179,
'loglikelihood': -87302.86681521917,
'redshift': 0.3691579561315956,
'q': 5.000623665102948,
'ra': 3.3158950374809972,
'spin1_azimuthal': 4.4461024300156655,
'logwt': -87343.31265112931,
'comoving_volume': 13805537798.745253,
'spin2_azimuthal': 0.8901187496181328,
'spin1_polar': 2.2912715299376387,
'coa_phase': 1.7398043595152775,
'chi_eff': -0.5769320542453128,
'spin2_a': 0.2699502327574271,
'spin1_a': 0.9829197339775988,
'approximant': 'IMRPhenomXPHM',
'f_lower': 20.0,
'f_ref': 20.0,
'trigger_time': 1242442967.4473,
'spin1x': -0.19437799854703947,
'spin1y': -0.7126224718025379,
'spin1z': -0.6484732915128284,
'spin2x': 0.09917070415014632,
'spin2y': 0.1224658286771233,
'spin2z': -0.21918125013462872,
'tc': 1242442967.450519,
'mass1': 203.08738423886746,
'mass2': 40.612411139058686}
%% Cell type:code id:a28763a5 tags:
``` python
```
notebook/reweightmpv_plot.ipynb 0 → 100644
+ 1932
− 0
View file @ 657e874a
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

AEI Privacy Policy Imprint