In [1]:
%pylab

Using matplotlib backend: <object object at 0x7f219c099de0>
Populating the interactive namespace from numpy and matplotlib
In [2]:
%matplotlib inline
In [3]:
from pycbc.workflow import WorkflowConfigParser
from pycbc.inference import models,io
from pycbc.inference.io import DynestyFile
import h5py,pycbc
from pycbc.waveform.utils import apply_fd_time_shift
from pycbc.types import TimeSeries
In [4]:
pycbc.init_logging(True)
In [5]:
cp = WorkflowConfigParser(['./inference-2048.ini'])
In [6]:
#fp = io.loadfile('GW150914-ringdown220.hdf','r')
#data = fp.read_data()
#psd = fp.read_psds()

#model = models.read_from_config(cp,data=data,psds=psd)
In [7]:
model = models.read_from_config(cp)

Setting up priors for each parameter
No sampling_params section read from config file
Loading waveform transforms
Determining analysis times to use
Padding H1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
Padding L1 analysis start and end times by 4 (= psd-inverse-length/2) seconds to account for PSD wrap around effects.
Checking that H1 has valid data in the requested segment
Checking that L1 has valid data in the requested segment
Reading Frames
Highpass Filtering
Resampling data
Converting to float64
Highpass Filtering
Remove Padding
Reading Frames
Highpass Filtering
Resampling data
Converting to float64
Highpass Filtering
Remove Padding
Will generate a different time series for PSD estimation
Checking that H1 has valid data in the requested segment
Checking that L1 has valid data in the requested segment
Reading Frames
Highpass Filtering
Resampling data
Converting to float64
Highpass Filtering
Remove Padding
Reading Frames
Highpass Filtering
Resampling data
Converting to float64
Highpass Filtering
Remove Padding
Applying gates to PSD data
Using Hann window to truncate inverse PSD
Setting low frequency cutoff of PSD to 0
Will highpass waveforms at 15.000000 Hz
In [8]:
with DynestyFile('GW150914-ringdown220.hdf','w') as f:
    f.write_psd(model.psds)
    f.write_stilde(model.data)
    f.write_config_file(cp)
In [9]:
#Choose an arbitrary set of parameters, for demonstration purpose
params = {'final_mass':60,
        'final_spin':0.8,
         'amp220':5e-20,
         'phi220':0}
model.update(**params)
In [10]:
_ = model.loglikelihood
In [11]:
model.loglikelihood
Out[11]:
-336914.711623026
In [12]:
model.lognl
Out[12]:
-325057.6323774979
In [13]:
model.loglr
Out[13]:
-11857.07924552809
In [14]:
print(sqrt(2*model.loglr))

nan

/work/yifan.wang/ringdown/GW150914/pycbc_run/subgating/env-subgating/lib/python3.7/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in sqrt
  """Entry point for launching an IPython kernel.
In [15]:
model.static_params
Out[15]:
{'approximant': 'TdQNMfromFinalMassSpin',
 'harmonics': 'spheroidal',
 'tref': 1126259462.4083147,
 'ra': 1.95,
 'dec': -1.27,
 'toffset': 0.0,
 'lmns': 221.0,
 'polarization': 0.82,
 'inclination': 3.14159265}

Plot data and waveforms (in Hanford)¶

In [16]:
data = model.td_data['H1']
In [17]:
h = model.get_waveforms()
ht = h['H1'].to_timeseries()
In [18]:
gate_times = model.get_gate_times()
In [19]:
start,delay = gate_times['H1']
In [26]:
# the ringdown starting time tc is t_gate_end
tc = start + delay 
# counts the index of the nearest data sample on the floor
idx = np.floor(float(tc-ht.start_time)*ht.sample_rate).astype(int)
# the new starting time, hence t_gate_end lands on a specific data sample 
disctc = ht.sample_times[idx]
# shift the waveform
shifth = apply_fd_time_shift(h['H1'], disctc-(tc-h['H1'].epoch), copy=True)
shiftht = shifth.to_timeseries()
# shift the data
fdata = data.to_frequencyseries()
newdata = apply_fd_time_shift(fdata,disctc-tc+fdata.epoch,copy=True)
newdatat = newdata.to_timeseries()
In [29]:
plt.plot(data.sample_times,data,'o-',label='data')
plt.plot(ht.sample_times,ht,'o-',label='ht')
plt.plot(shiftht.sample_times,shiftht,'o-',label='shifted ht')
plt.plot(newdatat.sample_times,newdatat,'o-',label='shift data')
tc = start+delay
plt.axvline(tc,ls='--',color='grey',label='tc')
plt.axvline(disctc,ls='--',color='grey',label='shifted tc')
plt.xlim(tc-0.001,tc+0.005)
plt.ylim(-3e-20,3e-20)
plt.legend()
plt.savefig('shiftgating.png',,bbox_inches='tight')

  File "/local/user/yifan.wang/ipykernel_404373/2935839773.py", line 11
    plt.savefig('shiftgating.png',,bbox_inches='tight')
                                  ^
SyntaxError: invalid syntax
In [ ]: