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