gwfast TUTORIAL¶
[1]:
import os
import sys
import copy
import numpy as onp
from astropy.cosmology import Planck18
PACKAGE_PARENT = '..'
SCRIPT_DIR = os.path.dirname(os.path.realpath(os.path.join(os.getcwd())))
sys.path.append(SCRIPT_DIR)
[2]:
import gwfast.gwfastGlobals as glob
COMPLETE EXAMPLE: GW170817¶
We use the positions of existing LVC detectors¶
[3]:
alldetectors = copy.deepcopy(glob.detectors)
print('All available detectors are: '+str(list(alldetectors.keys())))
# select only LIGO and Virgo
LVdetectors = {det:alldetectors[det] for det in ['L1', 'H1', 'Virgo']}
print('Using detectors '+str(list(LVdetectors.keys())))
All available detectors are: ['L1', 'H1', 'Virgo', 'KAGRA', 'LIGOI', 'ETS', 'ETMR', 'CE1Id', 'CE2NM', 'CE2NSW']
Using detectors ['L1', 'H1', 'Virgo']
[4]:
# We use the O2 psds
LVdetectors['L1']['psd_path'] = os.path.join(glob.detPath, 'LVC_O1O2O3', '2017-08-06_DCH_C02_L1_O2_Sensitivity_strain_asd.txt')
LVdetectors['H1']['psd_path'] = os.path.join(glob.detPath, 'LVC_O1O2O3', '2017-06-10_DCH_C02_H1_O2_Sensitivity_strain_asd.txt')
LVdetectors['Virgo']['psd_path'] = os.path.join(glob.detPath, 'LVC_O1O2O3', 'Hrec_hoft_V1O2Repro2A_16384Hz.txt')
[5]:
from gwfast.waveforms import IMRPhenomD_NRTidalv2
from gwfast.signal import GWSignal
from gwfast.network import DetNet
from fisherTools import CovMatr, compute_localization_region, check_covariance, fixParams
WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)
Initialise the signals and then the network¶
[6]:
myLVSignals = {}
for d in LVdetectors.keys():
myLVSignals[d] = GWSignal(IMRPhenomD_NRTidalv2(),
psd_path=LVdetectors[d]['psd_path'],
detector_shape = LVdetectors[d]['shape'],
det_lat= LVdetectors[d]['lat'],
det_long=LVdetectors[d]['long'],
det_xax=LVdetectors[d]['xax'],
verbose=True,
useEarthMotion = False,
fmin=10.,
IntTablePath=None)
myLVNet = DetNet(myLVSignals)
Using ASD from file /Users/francesco.iacovelli/Desktop/PhD/Research/2023-04_GWFAST/gwfast/psds/LVC_O1O2O3/2017-08-06_DCH_C02_L1_O2_Sensitivity_strain_asd.txt
WARNING: the motion of Earth gives a relevant contribution for BNS signals, consider switching it on
Initializing jax...
Jax local device count: 8
Jax device count: 8
Using ASD from file /Users/francesco.iacovelli/Desktop/PhD/Research/2023-04_GWFAST/gwfast/psds/LVC_O1O2O3/2017-06-10_DCH_C02_H1_O2_Sensitivity_strain_asd.txt
WARNING: the motion of Earth gives a relevant contribution for BNS signals, consider switching it on
Initializing jax...
Jax local device count: 8
Jax device count: 8
Using ASD from file /Users/francesco.iacovelli/Desktop/PhD/Research/2023-04_GWFAST/gwfast/psds/LVC_O1O2O3/Hrec_hoft_V1O2Repro2A_16384Hz.txt
WARNING: the motion of Earth gives a relevant contribution for BNS signals, consider switching it on
Initializing jax...
Jax local device count: 8
Jax device count: 8
[7]:
from gwfastUtils import GPSt_to_LMST
# Median values of the posterior samples for all the parameters,
# except psi and the coalescence phase that are set to 0
z = onp.array([0.00980])
tGPS = onp.array([1187008882.4])
GW170817 = {'Mc':onp.array([1.1859])*(1.+z),
'dL':Planck18.luminosity_distance(z).value/1000.,
'theta':onp.array([onp.pi/2. + 0.4080839999999999]),
'phi':onp.array([3.4461599999999994]),
'iota':onp.array([2.545065595974997]),
'psi':onp.array([0.]),
'tcoal':GPSt_to_LMST(tGPS, lat=0., long=0.), # GMST is LMST computed at long = 0°
'eta':onp.array([0.24786618323504223]),
'Phicoal':onp.array([0.]),
'chi1z':onp.array([0.005136138323169717]),
'chi2z':onp.array([0.003235146993487445]),
'Lambda1':onp.array([368.17802383555687]),
'Lambda2':onp.array([586.5487031450857])
}
print('Parameters for GW170817 are:')
GW170817
Parameters for GW170817 are:
[7]:
{'Mc': array([1.19752182]),
'dL': array([0.04374755]),
'theta': array([1.97888033]),
'phi': array([3.44616]),
'iota': array([2.5450656]),
'psi': array([0.]),
'tcoal': DeviceArray([0.43432288], dtype=float64),
'eta': array([0.24786618]),
'Phicoal': array([0.]),
'chi1z': array([0.00513614]),
'chi2z': array([0.00323515]),
'Lambda1': array([368.17802384]),
'Lambda2': array([586.54870315])}
Compute the expected matched-filter SNR (true is 33)¶
[8]:
SNR = myLVNet.SNR(GW170817)
print('SNR for GW170817 is %.2f to compare with 33'%SNR[0])
SNR for GW170817 is 33.16 to compare with 33
Compute the total Fisher matrix¶
[9]:
totF = myLVNet.FisherMatr(GW170817)
print('The computed Fisher matrix has shape %s'%str(totF.shape))
Computing Fisher for L1...
Computing derivatives...
Computing Fisher for H1...
Computing derivatives...
Computing Fisher for Virgo...
Computing derivatives...
Done.
The computed Fisher matrix has shape (13, 13, 1)
[10]:
# Check e.g. that the (dL,dL) element corresponds to (SNR/dL)^2
ParNums = IMRPhenomD_NRTidalv2().ParNums
dL_Num = ParNums['dL']
print('The relative difference is %.2e !'%((1 - totF[ParNums['dL'],ParNums['dL'],:]/(SNR/GW170817['dL'])**2)[0]))
The relative difference is 2.22e-16 !
Compute the covariance and perform some checks¶
[11]:
totCov, inversion_err = CovMatr(totF)
[12]:
_ = check_covariance(totF, totCov)
Inversion errors: [5657.625]
diagonal-1 = [array([-7.49594165e-09, 5.77928176e-08, -4.95498282e-14, 9.54375468e-14,
6.81787959e-13, 2.98029830e-13, 8.52979909e-11, 1.76168555e-08,
1.16116894e-09, -5.50993718e-06, -4.76458808e-06, -4.77896833e-10,
-5.41731424e-08], dtype=float128)]
Max off diagonal: [5657.625]
mask: where F*S(off-diagonal)>1e-10 (--> problematic if True off diagonal)
[array([[ True, False, False, False, False, False, False, False, False,
False, False, False, False],
[False, True, False, True, True, False, True, True, False,
True, True, False, False],
[False, True, True, False, False, False, False, True, False,
True, False, False, False],
[False, True, False, True, False, False, False, False, False,
False, False, False, False],
[False, True, False, False, True, False, False, False, False,
False, False, False, False],
[ True, False, False, False, False, True, False, False, False,
False, False, False, False],
[ True, False, False, False, False, False, True, False, False,
False, False, False, False],
[False, True, False, False, False, False, False, True, False,
True, False, False, False],
[ True, False, False, False, False, False, False, True, True,
True, False, False, False],
[False, True, False, True, True, False, True, False, False,
True, True, False, False],
[ True, False, True, True, True, True, False, True, True,
True, True, False, False],
[ True, False, False, False, False, True, True, False, True,
False, False, True, False],
[ True, False, True, False, False, True, False, True, True,
True, False, False, True]])]
Now try to eliminate the row corresponding to \(\delta\tilde{\Lambda}\), and see that the inversion error lowers¶
[13]:
ParNums = IMRPhenomD_NRTidalv2().ParNums
newFish, newPars = fixParams(totF, ParNums, ['deltaLambda'])
print('Now the Fisher matrix has shape %s'%str(newFish.shape))
newCov, new_inversion_err = CovMatr(newFish)
_ = check_covariance(newFish, newCov)
Now the Fisher matrix has shape (12, 12, 1)
Inversion errors: [0.00099936]
diagonal-1 = [array([-2.28188846e-09, -2.06082540e-09, 1.96557180e-13, 1.75944762e-13,
-9.95434637e-13, 2.58760678e-13, -9.65491662e-12, 4.94747979e-10,
6.65625888e-12, -1.89773989e-07, 7.42870654e-09, 3.38172546e-12],
dtype=float128)]
Max off diagonal: [1.3847575246472842991e-05]
mask: where F*S(off-diagonal)>1e-10 (--> problematic if True off diagonal)
[array([[ True, False, False, False, False, False, False, False, False,
False, False, False],
[ True, True, False, False, False, False, False, True, False,
True, False, False],
[False, False, True, False, False, False, False, False, False,
False, False, False],
[False, False, False, True, False, False, False, False, False,
False, False, False],
[ True, False, False, False, True, False, False, False, False,
False, False, False],
[ True, False, False, False, False, True, False, False, False,
False, False, False],
[False, True, False, False, False, False, True, True, False,
True, True, False],
[ True, False, False, False, False, False, False, True, False,
True, False, False],
[ True, False, False, False, False, False, False, False, True,
False, False, False],
[False, True, False, True, True, False, True, False, False,
True, True, False],
[False, True, False, True, True, False, False, True, False,
True, True, False],
[False, True, False, True, True, False, True, False, False,
False, True, True]])]
Finally compute the localisation region¶
[14]:
skyArea = compute_localization_region(newCov, newPars, GW170817['theta'])
print('The estimated sky location is %.1f deg^2, to compare with 16 deg^2'%skyArea)
The estimated sky location is 19.9 deg^2, to compare with 16 deg^2
COMPLETE EXAMPLE: ET and multiple events together¶
Configure the interforometer’s properties¶
[15]:
# Configure ET and the PSD
ETdet = {'ET': copy.deepcopy(glob.detectors).pop('ETS') }
print(ETdet)
ETdet['ET']['psd_path'] = os.path.join(glob.detPath, 'ET-0000A-18.txt')
{'ET': {'lat': 40.516666666666666, 'long': 9.416666666666666, 'xax': 0.0, 'shape': 'T'}}
Build the DetNet object¶
[16]:
from gwfast.waveforms import TaylorF2_RestrictedPN
from gwfast.signal import GWSignal
from gwfast.network import DetNet
[17]:
mySignalsET = {}
for d in ETdet.keys():
mySignalsET[d] = GWSignal(TaylorF2_RestrictedPN(use_3p5PN_SpinHO=True, is_tidal=True),
psd_path= ETdet[d]['psd_path'],
detector_shape = ETdet[d]['shape'],
det_lat= ETdet[d]['lat'],
det_long=ETdet[d]['long'],
det_xax=ETdet[d]['xax'],
verbose=True,
useEarthMotion = True,
fmin=2.,
IntTablePath=None)
myNet = DetNet(mySignalsET)
Using ASD from file /Users/francesco.iacovelli/Desktop/PhD/Research/2023-04_GWFAST/gwfast/psds/ET-0000A-18.txt
Initializing jax...
Jax local device count: 8
Jax device count: 8
Sample some BNS-like events¶
[18]:
nevents=100
zs = onp.random.uniform(1e-2, 3., nevents)
dLs = Planck18.luminosity_distance(zs).value/1000
Mcs = onp.random.normal(loc=1.156, scale=0.056, size=nevents)
events_rand = {'Mc': Mcs*(1.+zs),
'eta': onp.random.uniform(0.24, 0.25, nevents),
'dL': dLs,
'theta':onp.arccos(onp.random.uniform(-1., 1., nevents)),
'phi':onp.random.uniform(0., 2.*onp.pi, nevents),
'iota':onp.arccos(onp.random.uniform(-1., 1., nevents)),
'psi':onp.random.uniform(0., 2.*onp.pi, nevents),
'tcoal':onp.random.uniform(0., 1., nevents),
'Phicoal': onp.random.uniform(0., 2.*onp.pi, nevents),
'chi1z':onp.random.uniform(-.05, .05, nevents),
'chi2z':onp.random.uniform(-.05, .05, nevents),
'Lambda1':onp.random.uniform(0., 2000., nevents),
'Lambda2':onp.random.uniform(0., 2000., nevents),
}
Then computing SNRs and Fisher matrices is as easy and fast as¶
[19]:
%%time
snrs = myNet.SNR(events_rand)
CPU times: user 1.12 s, sys: 76 ms, total: 1.2 s
Wall time: 1.14 s
[20]:
%%time
totF = myNet.FisherMatr(events_rand)
Computing Fisher for ET...
Computing derivatives...
Filling matrix for arm 1...
Computing derivatives...
Filling matrix for arm 2...
Filling matrix for arm 3...
Done.
CPU times: user 5min 23s, sys: 4min 34s, total: 9min 58s
Wall time: 3min 27s
Both the signal and the network objects also contain a function to compute the optimal location in the sky to observe a binary with the considered configuration at a given time¶
[21]:
myNet.optimal_location(0.)
[21]:
array([0.8636426, 0.1643505])
(This has some caveats, refer to the paper for discussion)
This can be used e.g. to compute the detector reach for a given kind of source, as¶
[22]:
# Consider an equal mass non-spinning BNS system of 1.4 Msun, optimally oriented
# Notice that we here include in the dictionary the source-frame chirp mass
event = {'Mc_src':onp.array([1.2187707886145736]), 'eta':onp.array([.25]), 'iota':onp.array([0.]),
'psi':onp.array([0.]), 'tcoal':onp.array([0.]), 'Phicoal':onp.array([0.]),
'chi1z':onp.array([0.]), 'chi2z':onp.array([0.]),
'Lambda1':onp.array([0.]), 'Lambda2':onp.array([0.])}
[23]:
# We use Planck 18 cosmology to convert redshifts into luminosity distances
def get_zreach_event(detnet, event, SNRth=12, mtd='Nelder-Mead'):
from scipy.optimize import minimize
# Compute the best location and use it
best_theta, best_phi = detnet.optimal_location(event['tcoal'], is_tGPS=False)
event['theta'] = best_theta
event['phi'] = best_phi
def SNRofz(z):
event['Mc'] = event['Mc_src']*(1+z)
event['dL'] = Planck18.luminosity_distance(z)/1000.
return abs(detnet.SNR(event)-SNRth)
return minimize(SNRofz, 1, method=mtd).x
[24]:
get_zreach_event(myNet, event, SNRth=10)
[24]:
array([2.67226563])