gwnrtools / nr-catalog-tools Goto Github PK

View Code? Open in Web Editor NEW

1.0 3.0 5.0 8.64 MB

A unified interface to various catalogs of Numerical Relativity simulations of compact binary mergers.

Home Page: https://github.com/gwnrtools/nr-catalog-tools

License: GNU General Public License v3.0

Python 100.00%

binary-black-holes gravitational-waveforms gravitational-waves numerical-relativity binary-neutron-stars

nr-catalog-tools's Introduction

Interface to Numerical Relativity Catalogs

The nrcatalogtools python package provides a unified high-level interface to multiple catalogs of data products from Numerical Relativity simulations of compact-object binary mergers. At the moment of writing, different research groups have separate formats of data and/or tools to interface with them. This package will be a convenience layer atop those for downstream research applications. This package derives where it can from the sxs package of the Simulating eXtreme Spacetimes collaboration, as that has much of the capabilities for handling Numerical Relativity data that we need in general for all catalogs. Much of the interface here is deliberately identical to that in the sxs package to facilitate their interuse.

We currently support the following catalogs:

Usage

>>> from nrcatalogtools import RITCatalog
>>> rcatalog = RITCatalog.load()
>>> print(rcatalog.simulations_dataframe.index)
Index(['RIT:BBH:0001-n100-id3', 'RIT:BBH:0002-n100-id0',
       'RIT:BBH:0003-n100-id0', 'RIT:BBH:0004-n100-id0',
       'RIT:BBH:0005-n100-id0', 'RIT:BBH:0006-n100-id3',
       'RIT:BBH:0007-n100-id0', 'RIT:BBH:0008-n100-id0',
       'RIT:BBH:0009-n100-id0', 'RIT:BBH:0010-n100-id0',
       ...
       'RIT:BBH:1914-n144-id1', 'RIT:BBH:1915-n144-id1',
       'RIT:BBH:1916-n100-id1', 'RIT:BBH:1917-n100-id1',
       'RIT:BBH:1918-n100-id1', 'RIT:BBH:1919-n100-id1',
       'RIT:BBH:1920-n100-id1', 'RIT:BBH:1921-n100-id1',
       'RIT:BBH:1922-n100-id1', 'RIT:BBH:1923-n100-id1'],
      dtype='object', length=1879)

Now, if one needs a particular simulation, they can do:

>>> rwf = rcatalog.get('RIT:BBH:0003-n100-id0')

To check which modes are available for this simulation:

>>> print(rwf.LM)
[[ 2 -2]
 [ 2 -1]
 [ 2  0]
 [ 2  1]
 [ 2  2]
 [ 3 -3]
 [ 3 -2]
 [ 3 -1]
 [ 3  0]
 [ 3  1]
 [ 3  2]
 [ 3  3]
 [ 4 -4]
 [ 4 -3]
 [ 4 -2]
 [ 4 -1]
 [ 4  0]
 [ 4  1]
 [ 4  2]
 [ 4  3]
 [ 4  4]]

To extract a single mode from this:

>>> rwf.get_mode(2, 2)
array([[-1.18175000e+03,  8.41055081e-02,  6.60652456e-04],
       [-1.18150000e+03,  8.41034759e-02, -7.94687302e-04],
       [-1.18125000e+03,  8.40763695e-02, -2.25019642e-03],
       ...,
       [ 3.61000000e+02,  2.56889323e-12, -3.97799029e-25],
       [ 3.61250000e+02,  1.30444912e-12, -1.64922275e-25],
       [ 3.61500000e+02,  0.00000000e+00,  0.00000000e+00]])

To get polarizations for the same simulation:

>>> pols = rwf.get_td_waveform(total_mass = 40, # solar masses
                               distance = 100., # Megaparsecs
                               inclination = 0.2, # radians
                               coa_phase = 0.3) # radians
>>> hp, hc = pols.real(), -1 * pols.imag()

which can subsequently be plotted easily:

>>> import matplotlib.pyplot as plt
>>> plt.plot(hp.sample_times, hp, label='+')
>>> plt.plot(hc.sample_times, hc, label='x')
>>> plt.legend()
>>> plt.show()

which should give the following figure:

nr-catalog-tools's People

Contributors

Stargazers

Watchers

Forkers

vaishakp adivijaykumar prayush divyajyoti09 anuj137

nr-catalog-tools's Issues

Set up CI with basic tests

Need this to avoid issues like #19. For starters, could just include flake8 and black checks.

Possible inconsistency in loading of SXS data.

On calling

sxs_cat = SXSCatalog.load(download=True);
sxs_wf =  sxs1.get(sim_name, extrapolation_order=2)

The loader reports:

Found the following files to load from the SXS catalog:
    SXS:BBH:0001v6/Lev5/rhOverM_Asymptotic_GeometricUnits_CoM.h5

Further, on calling for the metadata:

sxs_cat.get_metadata(sim_name)['metadata_path']

It returns the metadata file path:

/home/user/.cache/sxs/SXS:BBH:0001v6/Lev5/metadata.json

However, it was found that the directory SXS:BBH:0001v6 contains only the metadata.json file and the directory SXS:BBH:0001v3 contains only the rhOverM_Asymptotic_GeometricUnits_CoM.h5 file.

On computing the md5 checksum of this h5 file, it was found to match with that of the v3 and v6 version of the simulation hosted on the SXS web catalog.

Is this correct behaviour?

Occasional freezes during GT catalog load

I have found that the execution freezes occasionally when attempting to load the GT catalog after displaying the message:
Downloading to {HOME}/.cache/MAYA/catalog.temp.txt:,
even if it has previously been downloaded. In this case, MAYACatalog.load shows no error but doesn't finish executing. It sometimes takes multiple restarts to get it to load.

Slow SXS catalog load

The execution of
sxs1 = SXSCatalog.load(download=False)
takes upwards of 290s, even if the previously downloaded catalog metadata is the latest.

Further, fetching modes data for a simulation that has been previously downloaded,
sxsw = sxs1.get(sim_name, extrapolation_order=2)
throws an error:
RuntimeError: Failed to download 'https://data.black-holes.org/catalog.json'; try setting download=False.

Store simulation file paths in simulation metadata.

It would be convenient to have the full path of the simulation metadata (used by the get_metadata method) and the file containing the modes data as an entry in the simulation metadata. It is recommended that this entry is set when the simulation metadata dictionary is created, and not computed again for the sake of consistency.

This would enable consistency checking between the H5 files and the simulation metadata loaded.

Speedup first time loading of RIT catalog

I was wondering if its possible to speedup the first time download/ load of RIT catalog using parallel requests. As I understand, the catalog presently downloads the metadata serially. It takes about half a day to download the catalog.

Would it be possible to cut down this time using wget's parallel download (using a constructed list of urls) or multiple processes?

Retire legacy code in `nrcatalogtools/legacy.py`

The code currently in nrcatalogtools/legacy.py was written to use waveform data from NR simulations, rotate and combine them appropriately to get observable GW polarizations. However, all of this functionality is now being derived from the sxs package which is more complete and computationally efficient than the implementation here. So we need to retire all the code in this file in order to avoid any confusion to users.

Translate inclination/ref_phase to LAL frame when calculating polarizations

At the moment the inclination and coa_phase angles used to construct the plus and cross polarizations by this member method of WaveformModes, are being interpreted in whatever frame the NR waveform modes data is available. In order to ensure uniformity of frames across NR catalogs, we should translate to the LAL frame that is documented here and currently used in lalsimulation. We choose this frame because lalsuite has the largest collection of GW waveform models from compact binary coalescences, and so NR waveform data is most likely to be used in conjunction with tools from lalsuite.

Test `get_td_waveform` method.

Test the functionality to compute polarizations from modes data.

Observations:

Apart from a coordinate rotation by pi/2 (this has not been added implicitly yet as it is immediately not clear if the rotation angle is the same across NR catalogs.), the waveform polarizations are differing with a mismatch of up to 1.5% w.r.t. the lal counterparts.

Case study: GT0001

The raw individual waveform modes loaded by nrcatalogtools, waveformtools and LAL are the same, upto numerical precision errors.
The individual waveform modes after resampling at user specified time-step, but before obtaining polarizations differ with a mismatch of about 0.7%. The plots (1) below show the differences in the (2,2) modes, after they have been suitably aligned in time, phase.

(if the image labels are not visible, please try turning off dark mode.)
3. The mismatches after obtaining the polarizations are about 1.3%.

Long SXS catalog loading times

The following code snippet takes a long time to execute every time I run it:

from nrcatalogtools import SXSCatalog
sxscatalog = SXSCatalog.load()

(Just to be clear, it's the second line that's taking long time to execute)

I already have catalog.zip stored in my ~/.cache/sxs directory, and still it takes a lot of time to load. I also tried to explicitly disable the downloading like so:

sxscatalog = SXSCatalog.load(download=False)

and it still takes long to execute.

Can this be cured or is this some "fundamental" I/O speed limitation in reading the catalog.zip file itself?

Interface with the new MAYA catalog data

Data is available here.

An interfacing package is available here.

Remove print statements, use logger instead

logger allows good control over what is displayed, and will also eliminate the need to use print_verbosity.

Uniform Formatting and PEP8 consistency

I saw a lot of flake8 failures and also noted that there is no uniform formatting across the package. For the latter, I prefer black, but I know people have major disagreements there. In any case, thoughts on these?

I can volunteer to put in a PR implementing these (possibly also adding flake8 checks to CI).

Pinging @prayush on this.

Add functionality to calculate error norm between two waveforms

PR #10 proposes to remove all legacy code, including calculation of the error norms between waveforms. However, it would be nice to have a version of the error norm compatible with this package. Will have to think a bit more about the best way to design this function, so creating this issue to make a note, and also to solicit input. The old code is pasted below:

    def get_L2_norm(self):
        sum_num = 0
        norm = 0
        for mode in self.modes:
            times_maya, hlm_maya = self.get_mode(mode=mode)
            times_sur, hlm_sur = self.get_mode_surrogate(mode=mode)
            phase_maya = np.unwrap(np.angle(hlm_maya))
            phase_maya = phase_maya - phase_maya[0]
            phase_sur = np.unwrap(np.angle(hlm_sur))
            phase_sur = phase_sur - phase_sur[0]
            amp_maya = np.abs(hlm_maya)
            amp_maya_interp = InterpolatedUnivariateSpline(
                times_maya, amp_maya)
            phase_maya_interp = InterpolatedUnivariateSpline(
                times_maya, phase_maya)
            amp_sur = np.abs(hlm_sur)
            aligned_strain_maya = (amp_maya_interp(times_sur) *
                                   np.exp(1j * phase_maya_interp(times_sur)))
            aligned_strain_sur = amp_sur * np.exp(1j * phase_sur)

            diff = np.sum(
                np.abs(aligned_strain_maya - aligned_strain_sur)**2 *
                self.delta_t)
            sum_num += diff
            norm += np.sum(np.abs(aligned_strain_maya)**2 * self.delta_t)
        return 0.5 * (sum_num / norm)

Need better documentation!

....everywhere.

Serious flake8 issues

Here's a summary of the flake8 errors introduced by #18. This needs to be fixed on priority.

./nrcatalogtools/lvc.py:454:5: F841 local variable 'LNhat' is assigned to but never used
./nrcatalogtools/lvc.py:458:5: F841 local variable 'nhat' is assigned to but never used
./nrcatalogtools/lvc.py:462:19: F821 undefined name 'LNhatx'
./nrcatalogtools/lvc.py:463:19: F821 undefined name 'LNhaty'
./nrcatalogtools/lvc.py:464:19: F821 undefined name 'LNhatz'
./nrcatalogtools/lvc.py:465:18: F821 undefined name 'nhatx'
./nrcatalogtools/lvc.py:466:18: F821 undefined name 'nhaty'
./nrcatalogtools/lvc.py:467:18: F821 undefined name 'nhatz'
./nrcatalogtools/lvc.py:496:5: F841 local variable 'IntParams' is assigned to but never used
./nrcatalogtools/lvc.py:496:17: F821 undefined name 'GetInterpRefValuesFromH5File'
./nrcatalogtools/lvc.py:499:13: F821 undefined name 'RefParams'
./nrcatalogtools/lvc.py:499:47: F821 undefined name 'RefParams'
./nrcatalogtools/lvc.py:500:13: F821 undefined name 'RefParams'
./nrcatalogtools/lvc.py:500:47: F821 undefined name 'RefParams'
./nrcatalogtools/lvc.py:501:13: F821 undefined name 'RefParams'
./nrcatalogtools/lvc.py:501:47: F821 undefined name 'RefParams'
./nrcatalogtools/lvc.py:505:5: F841 local variable 'nhat' is assigned to but never used
./nrcatalogtools/lvc.py:509:15: F821 undefined name 'RefParams'
./nrcatalogtools/lvc.py:510:15: F821 undefined name 'RefParams'
./nrcatalogtools/lvc.py:511:15: F821 undefined name 'RefParams'
./nrcatalogtools/lvc.py:519:19: F821 undefined name 'LNhatx'
./nrcatalogtools/lvc.py:520:19: F821 undefined name 'LNhaty'
./nrcatalogtools/lvc.py:521:19: F821 undefined name 'LNhatz'
./nrcatalogtools/lvc.py:522:18: F821 undefined name 'nhatx'
./nrcatalogtools/lvc.py:523:18: F821 undefined name 'nhaty'
./nrcatalogtools/lvc.py:524:18: F821 undefined name 'nhatz'
./nrcatalogtools/lvc.py:670:5: F841 local variable 'tol' is assigned to but never used
./nrcatalogtools/lvc.py:698:5: F841 local variable 'ReqDefAttrsSXS' is assigned to but never used
./nrcatalogtools/lvc.py:720:24: F821 undefined name 'GetRefTimeFromRefFreq'
./nrcatalogtools/lvc.py:722:42: F821 undefined name 'CheckInterpReq'
./nrcatalogtools/lvc.py:730:34: F821 undefined name 'CheckInterpReq'
./nrcatalogtools/lvc.py:737:16: F821 undefined name 'GetRefTimeFromMetadata'
./nrcatalogtools/lvc.py:744:40: F821 undefined name 'CheckNRAttrs'
./nrcatalogtools/lvc.py:754:48: F821 undefined name 'CheckNRAttrs'
./nrcatalogtools/lvc.py:760:33: F821 undefined name 'ComputeLALSourceFrameFromSXSMetadata'
./nrcatalogtools/lvc.py:768:29: F821 undefined name 'GetRefVals'
./nrcatalogtools/lvc.py:771:25: F821 undefined name 'GetRefVals'
./nrcatalogtools/lvc.py:792:48: F821 undefined name 'CheckNRAttrs'
./nrcatalogtools/lvc.py:800:25: F821 undefined name 'ComputeLALSourceFrameByInterp'
./nrcatalogtools/lvc.py:913:5: F841 local variable 'salpha' is assigned to but never used