add rational fits

gw_eccentricity/gw_eccentricity.py

@@ -314,6 +314,11 @@ def measure_eccentricity(tref_in=None,
             omega_gw_apocenters(t).
             Defaults are set using utils.get_default_spline_kwargs
 
+        rational_fit_kwargs: dict
+            Dictionary of arguments to be passed to the rational
+            fit function. Defaults are set using
+            `utils.get_default_rational_fit_kwargs`
+
         extrema_finding_kwargs:
             Dictionary of arguments to be passed to the extrema finder,
             scipy.signal.find_peaks.
@@ -412,6 +417,29 @@ def measure_eccentricity(tref_in=None,
             mean anomaly to zero.
             USE THIS WITH CAUTION!
 
+        omega_gw_extrema_interpolation_method : str, default="rational_fit"
+            Specifies the method used to build the interpolations for 
+            `omega_gw_pericenters_interp(t)` or `omega_gw_apocenters_interp(t)`.
+            The available options are:
+
+            - `spline`: Uses `scipy.interpolate.InterpolatedUnivariateSpline`.
+                - Best suited for cleaner data, such as when waveform modes are generated
+                    using models like SEOB or TEOB.
+                - Faster to construct and evaluate.
+                - Since it fits through every data point, it may exhibit oscillatory 
+                    behavior, particularly near the merger.
+
+            - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation`.
+                - Can handle both clean and noisy data, e.g., waveform modes
+                    from numerical simulations.
+                - Better monotonic behaviour, particularly near the merger.
+                - Significantly slower compared to the `spline` method. This is because
+                    finding optimal numerator and denominator degree needs several iterations
+                - Can suppress pathologies in the waveform that might be visible with 
+                    `spline`.
+
+            Default value: `"rational_fit"`.
+
     Returns
     -------
     A dictionary containing the following keys

gw_eccentricity/utils.py

@@ -1,6 +1,7 @@
 """Useful functions for the project."""
 import numpy as np
 import argparse
+from polyrat import StabilizedSKRationalApproximation
 from scipy.interpolate import InterpolatedUnivariateSpline
 from scipy.interpolate import PchipInterpolator
 import warnings
@@ -317,6 +318,42 @@ def get_interpolant(oldX,
     return interpolant
 
 
+def get_default_rational_fit_kwargs():
+    """Get default kwargs for rational fit."""
+    default_rational_fit_kwargs = {
+        "num_degree": None,
+        "denom_degree": None,
+        "norm": 2,
+        "maxiter": 20,
+        "verbose": False,
+        "xtol": 1e-07,
+    }
+    return default_rational_fit_kwargs
+
+
+def get_rational_fit(x, y, rational_fit_kwargs=None, check_kwargs=True):
+    """Get rational fit for the data set (x, y).
+
+    We use `polyrat.StabilizedSKRationalApproximation` to obtain
+    rational fit to the data.
+    """
+    if check_kwargs:
+        rational_fit_kwargs = check_kwargs_and_set_defaults(
+            rational_fit_kwargs,
+            get_default_rational_fit_kwargs(),
+            "rational_fit_kwargs",
+            "utils.get_default_rational_fit_kwargs"
+        )
+    # We use a rational approximation based on Stabilized Sanathanan-Koerner Iteration
+    # described in arXiv:2009.10803 and implemented in `polyrat.StabilizedSKRationalApproximation`.
+    rat = StabilizedSKRationalApproximation(**rational_fit_kwargs)
+    # The input x-axis data must be 2-dimensional
+    rat.fit(x.reshape(-1, 1), y)
+    # Using a lambda function to change the input 1-d data to
+    # float64 and reshape to 2-d as required by the fit.
+    return lambda t: rat(t.astype('float64').reshape(-1, 1))
+
+
 def debug_message(message, debug_level, important=True,
                   point_to_verbose_output=False):
     """Show message based on debug_level.

pytest.ini

@@ -0,0 +1,3 @@
+[pytest]  
+filterwarnings =  
+    ignore::DeprecationWarning:polyrat.*

setup.py

@@ -31,7 +31,8 @@ def get_version():
         'h5py',
         'lalsuite',
         'sxs',
-        'scri'
+        'scri',
+        'polyrat'
     ],
     classifiers=[
         "Intended Audience :: Science/Research",

test/generate_regression_data.py

@@ -17,10 +17,15 @@
     type=str,
     required=True,
     help="EccDefinition method to save the regression data for.")
+parser.add_argument(
+    "--interp_method",
+    type=str,
+    required=True,
+    help="omega_gw_extrema_interpolation_method to save the regression data for.")
 args = parser.parse_args()
 
 
-def generate_regression_data(method):
+def generate_regression_data(method, interp_method):
     """Generate data for regression test using a method."""
     # Load test waveform
     lal_kwargs = {"approximant": "EccentricTD",
@@ -42,7 +47,7 @@ def generate_regression_data(method):
         raise Exception(f"method {method} is not available. Must be one of "
                         f"{available_methods}")
 
-    extra_kwargs = {}
+    extra_kwargs = {"omega_gw_extrema_interpolation_method": interp_method}
     user_kwargs = extra_kwargs.copy()
     regression_data.update({"extra_kwargs": extra_kwargs})
     # Try evaluating at an array of times
@@ -56,9 +61,9 @@ def generate_regression_data(method):
     meanano_ref = gwecc_dict["mean_anomaly"]
     # We save the measured data 3 reference times
     n = len(tref_out)
-    dict_tref =  {"time": [tref_out[0], tref_out[n//4], tref_out[n//2]],
-                  "eccentricity": [ecc_ref[0], ecc_ref[n//4], ecc_ref[n//2]],
-                  "mean_anomaly": [meanano_ref[0], meanano_ref[n//4], meanano_ref[n//2]]}
+    dict_tref =  {"time": [tref_out[n//8], tref_out[n//4], tref_out[n//2]],
+                  "eccentricity": [ecc_ref[n//8], ecc_ref[n//4], ecc_ref[n//2]],
+                  "mean_anomaly": [meanano_ref[n//8], meanano_ref[n//4], meanano_ref[n//2]]}
 
     # Try evaluating at an array of frequencies
     gwecc_dict = measure_eccentricity(
@@ -70,18 +75,18 @@ def generate_regression_data(method):
     ecc_ref = gwecc_dict["eccentricity"]
     meanano_ref = gwecc_dict["mean_anomaly"]
     n = len(fref_out)
-    dict_fref = {"frequency": [fref_out[0], fref_out[n//4], fref_out[n//2]],
-                 "eccentricity": [ecc_ref[0], ecc_ref[n//4], ecc_ref[n//2]],
-                 "mean_anomaly": [meanano_ref[0], meanano_ref[n//4], meanano_ref[n//2]]}
+    dict_fref = {"frequency": [fref_out[n//8], fref_out[n//4], fref_out[n//2]],
+                 "eccentricity": [ecc_ref[n//8], ecc_ref[n//4], ecc_ref[n//2]],
+                 "mean_anomaly": [meanano_ref[n//8], meanano_ref[n//4], meanano_ref[n//2]]}
     regression_data.update({"tref": dict_tref,
                             "fref": dict_fref})
 
     if not os.path.exists(data_dir):
         os.mkdir(data_dir)
     # save to a json file
-    fl = open(f"{data_dir}/{method}_regression_data.json", "w")
+    fl = open(f"{data_dir}/{method}_{interp_method}_regression_data.json", "w")
     json.dump(regression_data, fl)
     fl.close()
 
 # generate regression data
-generate_regression_data(args.method)
+generate_regression_data(args.method, args.interp_method)

test/regression_data/AmplitudeFits_spline_regression_data.json

@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12981.119163243457, -11192.719163243455, -7616.019163243456], "eccentricity": [0.14125503147618124, 0.1347862810667858, 0.11943002063131647], "mean_anomaly": [3.271739184693624, 0.9679656489638688, 5.4841245849525535]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11926882656303317, 0.11542822840106648, 0.10530937888140646], "mean_anomaly": [5.832592108883006, 1.494348360501732, 2.830869757710083]}}

test/regression_data/Amplitude_spline_regression_data.json

@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12997.219163243457, -11225.819163243457, -7683.019163243456], "eccentricity": [0.1412963090542808, 0.13489390800343737, 0.11972197672456475], "mean_anomaly": [3.118093974362708, 0.642446628947404, 4.757738663295697]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11923643778091675, 0.11538856688636778, 0.10524430745877145], "mean_anomaly": [5.8083670411936055, 1.4699121058751956, 2.7932679992189193]}}

test/regression_data/FrequencyFits_spline_regression_data.json

@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12981.319163243457, -11193.019163243456, -7616.419163243456], "eccentricity": [0.1412559990968395, 0.13478783801299066, 0.11943256094180654], "mean_anomaly": [3.271739184693624, 0.9671579875661678, 5.4821375349876575]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11926984282942188, 0.1154273089414749, 0.10530750015791723], "mean_anomaly": [5.833828675147657, 1.5003684452927075, 2.839055270348034]}}

test/regression_data/Frequency_spline_regression_data.json

@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12997.819163243457, -11226.819163243457, -7684.919163243456], "eccentricity": [0.14129511303046816, 0.13489405452014136, 0.11972728556745593], "mean_anomaly": [3.114478792943153, 0.634809819093654, 4.7407656825570825]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11923247746634258, 0.11538428023179614, 0.10524044651877096], "mean_anomaly": [5.811428452219403, 1.473195504466105, 2.798061691253885]}}

test/regression_data/ResidualAmplitude_spline_regression_data.json

@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12981.319163243457, -11193.019163243456, -7616.519163243456], "eccentricity": [0.14125596891213044, 0.13478776432643125, 0.11943299066201574], "mean_anomaly": [3.2708353893387336, 0.9662032412803043, 5.4800155541341695]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11926975624561775, 0.1154272021945365, 0.10530996061768227], "mean_anomaly": [5.832818126814786, 1.4992647336341633, 2.8323037712231454]}}

test/regression_data/ResidualFrequency_spline_regression_data.json

@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12981.319163243457, -11193.019163243456, -7616.519163243456], "eccentricity": [0.14125598209018864, 0.1347878224760185, 0.1194330394418408], "mean_anomaly": [3.271739184693624, 0.967011047732143, 5.481076544560906]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11926981413331783, 0.11542727032728273, 0.10531052353279924], "mean_anomaly": [5.833866315156342, 1.5003849616201421, 2.8327681616289127]}}