Merge pull request #257 from vijayvarma392/add_tests

add more tests

Author: vijayvarma392

.github/workflows/test.yml

@@ -19,7 +19,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [macos-latest, ubuntu-latest]
-        python-version: ['3.9']
+        python-version: ['3.9', '3.10']
 
     steps:
       - name: Check out repository code
@@ -98,6 +98,18 @@ jobs:
         run: |
           pip install .
 
+      - name: Download SXS NR data for precessing tests
+        if: matrix.python-version != '3.9'
+        shell: bash -l {0}
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+        run: |
+          python - <<'EOF'
+          from gw_eccentricity.load_data import download_sxs_waveform
+          download_sxs_waveform("SXS:BBH:2859", 4, "/tmp/gwecc_test_data",
+                                catalog_tag="v3.0.0")
+          EOF
+
       - name: Run test suite
         shell: bash -l {0}
         run: |

gw_eccentricity/eccDefinition.py

@@ -2379,6 +2379,21 @@ def measure_ecc(self):
         self.tmax = min(self.t_pericenters[-1], self.t_apocenters[-1])
         self.tmin = max(self.t_pericenters[0], self.t_apocenters[0])
 
+        # t_for_checks depends only on tmin/tmax (already set above).
+        # Build it early so get_omega_gw_extrema_interpolant can use it when
+        # converting fref_in → tref_in via mean_of_extrema_interpolants.
+        self.t_for_checks = self.dataDict["t"][
+            np.logical_and(self.dataDict["t"] >= self.tmin,
+                           self.dataDict["t"] <= self.tmax)]
+
+        # Build omega_gw extrema interpolants before the fref_in conversion so
+        # that compute_mean_of_extrema_interpolants is available as an
+        # averaging method in compute_tref_in_and_fref_out_from_fref_in.
+        self.omega_gw_pericenters_interp \
+            = self.get_omega_gw_extrema_interpolant("pericenters")
+        self.omega_gw_apocenters_interp \
+            = self.get_omega_gw_extrema_interpolant("apocenters")
+
         if self.domain == "frequency":
             # get the tref_in and fref_out from fref_in
             self.tref_in, self.fref_out \
@@ -2387,10 +2402,6 @@ def measure_ecc(self):
         self.tref_out = self.tref_in[
             np.logical_and(self.tref_in <= self.tmax,
                            self.tref_in >= self.tmin)]
-        # set time for checks and diagnostics
-        self.t_for_checks = self.dataDict["t"][
-            np.logical_and(self.dataDict["t"] >= self.tmin,
-                           self.dataDict["t"] <= self.tmax)]
 
         # Sanity checks
         # check that fref_out and tref_out are of the same length
@@ -2427,11 +2438,6 @@ def measure_ecc(self):
            or self.tref_out[-1] > self.t_pericenters[-1]:
             raise Exception("Reference time must be within two pericenters.")
 
-        # Build omega_gw extrema interpolants
-        self.omega_gw_pericenters_interp \
-            = self.get_omega_gw_extrema_interpolant("pericenters")
-        self.omega_gw_apocenters_interp \
-            = self.get_omega_gw_extrema_interpolant("apocenters")
         # compute eccentricity at self.tref_out
         self.eccentricity = self.compute_eccentricity(self.tref_out)
         # Compute mean anomaly at tref_out

gw_eccentricity/load_data.py

@@ -743,6 +743,87 @@ def load_lvcnr_waveform(**kwargs):
     return return_dict
 
 
+def download_sxs_waveform(sxs_id, lev, data_dir, extrap_order=2,
+                          catalog_tag=""):
+    """Download SXS catalog waveform files to a local directory.
+
+    Downloads the files required by ``load_waveform(origin="SXSCatalog")``:
+
+    - ``Strain_N{extrap_order}.h5``
+    - ``Strain_N{extrap_order}.json``
+    - ``metadata.json``
+    - ``Horizons.h5``
+
+    Parameters
+    ----------
+    sxs_id : str
+        SXS simulation identifier, e.g. ``"SXS:BBH:3726"``.
+    lev : int
+        Resolution level, e.g. ``3``.
+    data_dir : str
+        Root directory.  Files are placed inside
+        ``<data_dir>/<sxs_id_safe>/Lev<lev>/``, where ``sxs_id_safe``
+        is ``sxs_id`` with ``":"`` replaced by ``"_"``
+        (e.g. ``"SXS_BBH_3726"``).
+    extrap_order : int, optional
+        Waveform extrapolation order.  Default is 2, which selects the
+        ``Strain_N2.h5`` / ``Strain_N2.json`` file pair.
+    catalog_tag : str, optional
+        Git tag of the SXS catalog release to use, e.g. ``"v3.0.0"``.
+        When provided the catalog is fetched directly from that tagged
+        release, bypassing the GitHub API call that discovers the latest
+        release.  This avoids GitHub API rate-limit errors in CI.
+        Defaults to ``""`` which fetches the latest release.
+
+    Returns
+    -------
+    str
+        Full path to the directory that now contains the downloaded
+        files.  Pass this as ``data_dir`` to
+        ``load_waveform(origin="SXSCatalog", data_dir=...)``.
+
+    Notes
+    -----
+    Files that already exist on disk are skipped, that is not re-downloaded.
+    The SXS simulations catalogue is loaded via the ``sxs`` Python package to
+    retrieve direct download URLs for each file.
+    """
+    from sxs.utilities import download_file
+
+    # Build target directory
+    sxs_id_safe = sxs_id.replace(":", "_")
+    target_dir = os.path.join(
+        os.path.expanduser(data_dir), sxs_id_safe, f"Lev{lev}")
+    os.makedirs(target_dir, exist_ok=True)
+
+    # Load the simulations catalogue to resolve file URLs.
+    # Pass `tag` when provided to skip the GitHub API call that discovers
+    # the latest release (avoids rate-limit errors in CI environments).
+    simulations = sxs.Simulations.load(download=True, tag=catalog_tag)
+    if sxs_id not in simulations:
+        raise ValueError(
+            f"'{sxs_id}' not found in the SXS simulations catalogue.  "
+            "Check the SXS ID and that the catalogue is up to date.")
+
+    sim_files = simulations[sxs_id].get("files", {})
+
+    filenames = [
+        f"Strain_N{extrap_order}.h5",
+        f"Strain_N{extrap_order}.json",
+        "metadata.json",
+        "Horizons.h5",
+    ]
+
+    for filename in filenames:
+        dest = os.path.join(target_dir, filename)
+        if os.path.exists(dest):
+            continue  # already present, skip
+        url = sim_files[f"Lev{lev}:{filename}"]["link"]
+        download_file(url, dest, progress=False)
+
+    return target_dir
+
+
 def load_sxs_catalogformat(**kwargs):
     """Load modes from sxs waveform files in sxs catalog format.
 

gw_eccentricity/rational_fit.py

@@ -104,7 +104,8 @@ def _eval_arnoldi_basis(x, H, degree, v0_norm):
     # Build higher degree basis functions using the stored Hessenberg
     # matrix H
     for k in range(degree):
-        v = V[:, k] * x - V[:, :k+1] @ H[:k+1, k]
+        with np.errstate(divide='ignore', over='ignore', invalid='ignore'):
+            v = V[:, k] * x - V[:, :k+1] @ H[:k+1, k]
         if H[k + 1, k] < 1e-14:
             break
         V[:, k + 1] = v / H[k + 1, k]
@@ -405,8 +406,9 @@ def predict(self, x_new):
 
         # Compute the numerator and denominator at x_new using the
         # fitted coefficients
-        p = P @ self.a
-        q = Q @ self.b
+        with np.errstate(divide='ignore', over='ignore', invalid='ignore'):
+            p = P @ self.a
+            q = Q @ self.b
 
         # Return the rational function values at x_new
         r = p / q

pytest.ini

@@ -1 +1,5 @@
-[pytest]
+[pytest]
+filterwarnings =
+    ignore::DeprecationWarning:pyseobnr
+    ignore::DeprecationWarning:scipy
+    ignore::DeprecationWarning:qnm

test/conftest.py

@@ -0,0 +1,37 @@
+"""Shared test helpers and fixtures for the gw_eccentricity test suite."""
+import numpy as np
+from pyseobnr.generate_waveform import generate_modes_opt
+
+
+def get_seob_datadict(q=5.0, chi1z=0.4, chi2z=0.3,
+                      omega_start=0.02, ecc=0.1, mean_ano=1.7,
+                      include_zero_ecc=False):
+    """Return a dataDict built from a SEOBNRv5EHM waveform.
+
+    Parameters
+    ----------
+    q            : mass ratio
+    chi1z, chi2z : aligned spin components
+    omega_start  : initial orbital frequency in geometric units
+    ecc          : initial eccentricity
+    mean_ano     : initial mean anomaly
+    include_zero_ecc : if True, also add t_zeroecc/hlm_zeroecc keys built
+                       from a zero-eccentricity waveform starting slightly
+                       earlier (omega_start * 0.9), required by Residual* methods
+    """
+    t, modes = generate_modes_opt(
+        q, chi1z, chi2z, omega_start,
+        eccentricity=ecc, rel_anomaly=mean_ano,
+        approximant="SEOBNRv5EHM")
+    hlm = {tuple(int(x) for x in k.split(",")): v for k, v in modes.items()}
+    dataDict = {"t": t, "hlm": hlm}
+
+    if include_zero_ecc:
+        t_z, modes_z = generate_modes_opt(
+            q, chi1z, chi2z, omega_start * 0.9,
+            eccentricity=0, rel_anomaly=0,
+            approximant="SEOBNRv5EHM")
+        hlm_z = {tuple(int(x) for x in k.split(",")): v for k, v in modes_z.items()}
+        dataDict.update({"t_zeroecc": t_z, "hlm_zeroecc": hlm_z})
+
+    return dataDict

test/test_mks_vs_dimless_units.py

@@ -5,6 +5,13 @@
 from lal import MTSUN_SI
 
 
+# rational_fit introduces ~2x more floating-point noise than spline due to its
+# iterative QR algorithm, pushing unit-conversion differences to ~2e-7.
+# rtol=1e-6 comfortably covers both methods while still catching real unit bugs.
+# (numpy assert_allclose default is 1e-7)
+UNIT_CONSISTENCY_RTOL = 1e-6
+
+
 def test_mks_vs_dimless_units():
     """ Tests that the measure_eccentricity interface is working for both
     MKS and dimensionless units.
@@ -38,9 +45,7 @@ def test_mks_vs_dimless_units():
                        "D": 1})
     dataDictMKS = load_data.load_waveform(**lal_kwargs)
 
-    # use different omega_gw_extrema_interpolation methods
-    # TODO: Add rational_fit to the list
-    omega_gw_extrema_interpolation_methods = ["spline"]
+    omega_gw_extrema_interpolation_methods = ["spline", "rational_fit"]
 
     # List of all available methods
     available_methods = gw_eccentricity.get_available_methods()
@@ -74,19 +79,22 @@ def test_mks_vs_dimless_units():
             np.testing.assert_allclose(
                 [tref_out],
                 [tref_out_MKS * sec_to_dimless],
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("tref_out at a single dimensionless and MKS"
                          " time are inconsistent.\n"
                          "x = Dimensionless, y = MKS converted to dimless"))
             # Check if the measured ecc an mean ano are the same from the two units
             np.testing.assert_allclose(
                 [ecc_ref],
                 [ecc_ref_MKS],
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("Eccentricity at a single dimensionless and MKS"
                          " time gives different results.\n"
                          "x = Dimensionless, y = MKS"))
             np.testing.assert_allclose(
                 [meanano_ref],
                 [meanano_ref_MKS],
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("Mean anomaly at a single dimensionless and MKS"
                          " time gives different results.\n"
                          "x = Dimensionless, y = MKS"))
@@ -115,13 +123,15 @@ def test_mks_vs_dimless_units():
             np.testing.assert_allclose(
                 [tref_out],
                 [tref_out_MKS * sec_to_dimless],
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("tref_out array for dimensionless and MKS"
                          " tref_in are inconsistent.\n"
                          "x = Dimensionless, y = MKS converted to dimless"))
             # Check if the measured ecc an mean ano are the same from the two units
             np.testing.assert_allclose(
                 ecc_ref,
                 ecc_ref_MKS,
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("Eccentricity at dimensionless and MKS array of"
                          " times are different\n."
                          "x = Dimensionless, y = MKS"))
@@ -130,6 +140,7 @@ def test_mks_vs_dimless_units():
             np.testing.assert_allclose(
                 np.unwrap(meanano_ref),
                 np.unwrap(meanano_ref_MKS),
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("Mean anomaly at dimensionless and MKS array of"
                          " times are different.\n"
                          "x = Dimensionless, y = MKS"))
@@ -160,19 +171,22 @@ def test_mks_vs_dimless_units():
             np.testing.assert_allclose(
                 [fref_out],
                 [fref_out_MKS / sec_to_dimless],
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("fref_out for a single dimensionless and MKS"
                          " fref_in are inconsistent.\n"
                          "x = Dimensionless, y = MKS converted to dimless"))
             # Check if the measured ecc an mean ano are the same from the two units
             np.testing.assert_allclose(
                 [ecc_ref],
                 [ecc_ref_MKS],
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("Eccentricity at a single dimensionless and MKS"
                          " frequency gives different results.\n"
                          "x = Dimensionless, y = MKS"))
             np.testing.assert_allclose(
                 [meanano_ref],
                 [meanano_ref_MKS],
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("Mean anomaly at a single dimensionless and MKS"
                          " frequency gives different results.\n"
                          "x = Dimensionless, y = MKS"))
@@ -203,19 +217,22 @@ def test_mks_vs_dimless_units():
             np.testing.assert_allclose(
                 [fref_out],
                 [fref_out_MKS / sec_to_dimless],
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("fref_out for an array of dimensionless and MKS"
                          " fref_in are inconsistent.\n"
                          "x = Dimensionless, y = MKS converted to dimless"))
             # Check if the measured ecc an mean ano are the same from the two units
             np.testing.assert_allclose(
                 ecc_ref,
                 ecc_ref_MKS,
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("Eccentricity at dimensionless and MKS array of"
                          " frequencies are different.\n"
                          "x = Dimensionless, y = MKS"))
             np.testing.assert_allclose(
                 np.unwrap(meanano_ref),
                 np.unwrap(meanano_ref_MKS),
+                rtol=UNIT_CONSISTENCY_RTOL,
                 err_msg=("Mean anomaly at dimensionless and MKS array of"
                          " frequencies are different.\n"
                          "x = Dimensionless, y = MKS"))