From 07cd500a4580aaa9164ab50b95b87eb2fded29a2 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Wed, 25 Sep 2024 15:03:11 +0530 Subject: [PATCH 01/34] add rational fits --- gw_eccentricity/eccDefinition.py | 70 ++++++++++++++++++++++++++++++++ gw_eccentricity/utils.py | 32 +++++++++++++++ setup.py | 3 +- 3 files changed, 104 insertions(+), 1 deletion(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 8bd0f019..8fb4dfa7 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -14,6 +14,8 @@ from .utils import interpolate from .utils import get_interpolant from .utils import get_default_spline_kwargs +from .utils import get_rational_fit +from .utils import get_default_rational_fit_kwargs from .utils import debug_message from .plot_settings import use_fancy_plotsettings, colorsDict, labelsDict from .plot_settings import figWidthsTwoColDict, figHeightsDict @@ -310,9 +312,15 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, get_default_spline_kwargs(), "spline_kwargs", "utils.get_default_spline_kwargs()") + self.rational_fit_kwargs = check_kwargs_and_set_defaults( + self.extra_kwargs["rational_fit_kwargs"], + get_default_rational_fit_kwargs(), + "rational_fit_kwargs", + "utils.get_default_rational_fit_kwargs()") self.available_averaging_methods \ = self.get_available_omega_gw_averaging_methods() self.debug_level = self.extra_kwargs["debug_level"] + self.rational_fit_kwargs["verbose"] = True if self.debug_level >=1 else False self.debug_plots = self.extra_kwargs["debug_plots"] self.return_zero_if_small_ecc_failure = self.extra_kwargs["return_zero_if_small_ecc_failure"] # check if there are unrecognized keys in the dataDict @@ -823,6 +831,7 @@ def get_default_extra_kwargs(self): """Defaults for additional kwargs.""" default_extra_kwargs = { "spline_kwargs": {}, + "rational_fit_kwargs": {}, "extrema_finding_kwargs": {}, # Gets overridden in methods like # eccDefinitionUsingAmplitude "debug_level": 0, @@ -1174,6 +1183,53 @@ def interp_extrema(self, extrema_type="pericenters"): raise Exception( f"Sufficient number of {extrema_type} are not found." " Can not create an interpolant.") + + def rational_fit(self, x, y): + return get_rational_fit(x, y, + rational_fit_kwargs=self.rational_fit_kwargs, + check_kwargs=False) + + def rational_fit_extrema(self, extrema_type="pericenters"): + if extrema_type == "pericenters": + extrema = self.pericenters_location + elif extrema_type == "apocenters": + extrema = self.apocenters_location + else: + raise Exception("extrema_type must be either " + "'pericenrers' or 'apocenters'.") + if len(extrema) >= 2: + # assign degree + # TODO: Use an optimal degree based on the wavefrom length + self.rational_fit_kwargs["num_degree"],\ + self.rational_fit_kwargs["denom_degree"] = 2, 2 + rational_fit = self.rational_fit(self.t[extrema], + self.omega_gw[extrema]) + t = self.t_for_checks.copy() + omega = rational_fit(t) + while self.check_domega_dt(t, omega, 1.0): + self.rational_fit_kwargs["num_degree"] -= 1 + self.rational_fit_kwargs["denom_degree"] -= 1 + debug_message(f"degree is lowered to {self.rational_fit_kwargs['num_degree']}", + debug_level=self.debug_level, + important=True) + rational_fit = self.rational_fit(self.t[extrema], + self.omega_gw[extrema]) + omega = rational_fit(t) + return rational_fit + + else: + raise Exception( + f"Sufficient number of {extrema_type} are not found." + " Can not create an interpolant.") + + def check_domega_dt(self, t, omega, tol=1.0): + """Check first derivative of omega. + + The first derivative of interpolant/fit of omega at the extrema + should be monotonically increasing. + """ + domega_dt = np.gradient(omega, t[1] - t[0]) + return any(domega_dt[1:]/domega_dt[:-1] < tol) def check_num_extrema(self, extrema, extrema_type="extrema"): """Check number of extrema. @@ -1518,6 +1574,7 @@ def measure_ecc(self, tref_in=None, fref_in=None): self.t_apocenters = self.t[self.apocenters_location] self.tmax = min(self.t_pericenters[-1], self.t_apocenters[-1]) self.tmin = max(self.t_pericenters[0], self.t_apocenters[0]) + if self.domain == "frequency": # get the tref_in and fref_out from fref_in self.tref_in, self.fref_out \ @@ -1531,6 +1588,19 @@ def measure_ecc(self, tref_in=None, fref_in=None): np.logical_and(self.dataDict["t"] >= self.tmin, self.dataDict["t"] <= self.tmax)] + # Check monotonicity of the first derivative of omega_gw interpolant + # and build rational fit if necessary. + if (self.check_domega_dt(self.t_for_checks, + self.omega_gw_pericenters_interp(self.t_for_checks)) + or + self.check_domega_dt(self.t_for_checks, + self.omega_gw_apocenters_interp(self.t_for_checks))): + debug_message("The spline interpolant through extrema has non " + "monotonic time derivative. Trying rational fit instead.", + debug_level=self.debug_level, important=True) + self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") + self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") + # Sanity checks # check that fref_out and tref_out are of the same length if self.domain == "frequency": diff --git a/gw_eccentricity/utils.py b/gw_eccentricity/utils.py index e66b3647..29eed8cc 100644 --- a/gw_eccentricity/utils.py +++ b/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,37 @@ 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" + ) + rat = StabilizedSKRationalApproximation(**rational_fit_kwargs) + rat.fit(x.reshape(-1, 1), y) + return lambda x: rat(x.astype('float64').reshape(-1, 1)) + + def debug_message(message, debug_level, important=True, point_to_verbose_output=False): """Show message based on debug_level. diff --git a/setup.py b/setup.py index e3c66413..3f6dd633 100644 --- a/setup.py +++ b/setup.py @@ -30,7 +30,8 @@ def get_version(): 'scipy', 'h5py', 'lalsuite', - 'sxs' + 'sxs', + 'polyrat' ], classifiers=[ "Intended Audience :: Science/Research", From 107fbf2e27cf1a5068360f59d1db932b26235d8d Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Thu, 26 Sep 2024 08:51:46 +0530 Subject: [PATCH 02/34] add documentation and restructure rational_fit --- gw_eccentricity/eccDefinition.py | 69 +++++++++++++++++++++++--------- 1 file changed, 49 insertions(+), 20 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 8fb4dfa7..56a45bc8 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -1184,12 +1184,58 @@ def interp_extrema(self, extrema_type="pericenters"): f"Sufficient number of {extrema_type} are not found." " Can not create an interpolant.") - def rational_fit(self, x, y): + def get_rat_fit(self, x, y): + """Get interpolant. + + A wrapper of `utils.get_rational_fit` with check_kwargs=False. This is + to make sure that the checking of kwargs is not performed everytime the + rational fit function is called. Instead, the kwargs are checked once + in the init and passed to the rational fit function without repeating + checks. + """ return get_rational_fit(x, y, rational_fit_kwargs=self.rational_fit_kwargs, check_kwargs=False) + def rational_fit(self, x, y): + """Get rational fit with adaptive numerator and denominator degree. + + This function ensures that the rational fit we obtain is using the + optimal degree for the numerator and the denominator. We start with an + initial estimate of what these degrees should be based on the length of + the waveform. We check the first derivative of the resultant fit to + check for any nonmonotonicity. In case of nonmonocity detected, we + lower the degree by 1 and repeat until the check passes successfully. + """ + # assign degree + # TODO: Use an optimal degree based on the wavefrom length + self.rational_fit_kwargs["num_degree"],\ + self.rational_fit_kwargs["denom_degree"] = 2, 2 + rat_fit = self.get_rat_fit(x, y) + t = np.arange(x[0], x[-1], self.t[1] - self.t[0]) + omega = rat_fit(t) + while self.check_domega_dt(t, omega, 1.0): + self.rational_fit_kwargs["num_degree"] -= 1 + self.rational_fit_kwargs["denom_degree"] -= 1 + debug_message(f"degree is lowered to {self.rational_fit_kwargs['num_degree']}", + debug_level=self.debug_level, + important=True) + rat_fit = self.get_rat_fit(x, y) + omega = rational_fit(t) + return rat_fit + def rational_fit_extrema(self, extrema_type="pericenters"): + """Build rational fit through extrema. + + parameters: + ----------- + extrema_type: + Either "pericenters" or "apocenters". + + returns: + ------ + Rational fit through extrema + """ if extrema_type == "pericenters": extrema = self.pericenters_location elif extrema_type == "apocenters": @@ -1198,25 +1244,8 @@ def rational_fit_extrema(self, extrema_type="pericenters"): raise Exception("extrema_type must be either " "'pericenrers' or 'apocenters'.") if len(extrema) >= 2: - # assign degree - # TODO: Use an optimal degree based on the wavefrom length - self.rational_fit_kwargs["num_degree"],\ - self.rational_fit_kwargs["denom_degree"] = 2, 2 - rational_fit = self.rational_fit(self.t[extrema], - self.omega_gw[extrema]) - t = self.t_for_checks.copy() - omega = rational_fit(t) - while self.check_domega_dt(t, omega, 1.0): - self.rational_fit_kwargs["num_degree"] -= 1 - self.rational_fit_kwargs["denom_degree"] -= 1 - debug_message(f"degree is lowered to {self.rational_fit_kwargs['num_degree']}", - debug_level=self.debug_level, - important=True) - rational_fit = self.rational_fit(self.t[extrema], - self.omega_gw[extrema]) - omega = rational_fit(t) - return rational_fit - + return self.rational_fit(self.t[extrema], + self.omega_gw[extrema]) else: raise Exception( f"Sufficient number of {extrema_type} are not found." From 7b527fd5918d0381016b101e2b7336074d6ee58c Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Thu, 26 Sep 2024 11:28:05 +0530 Subject: [PATCH 03/34] add omega_gw_extrema_interpolation_method --- gw_eccentricity/eccDefinition.py | 88 ++++++++++++++++++++++++++------ 1 file changed, 73 insertions(+), 15 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 56a45bc8..a37fd5a9 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -176,6 +176,10 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, 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. + extrema_finding_kwargs: dict Dictionary of arguments to be passed to the extrema finder, scipy.signal.find_peaks. @@ -274,6 +278,42 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, eccentricity is the cause, and set the returned eccentricity and mean anomaly to zero. USE THIS WITH CAUTION! + + "omega_gw_extrema_interpolation_method" : str, default="spline" + Method to use for building `omega_gw_pericenters_interp(t)` or + `omega_gw_apocenters_interp(t)`. Available options are: + + - `spline`: Uses `scipy.interpolate.InterpolatedUnivariateSline`. + - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation` + + `spline` is useful when the data is cleaner. For example, when the + waveform modes are generated using an waveform model like SEOB/TEOB. + It is also faster to build and evaluate. Since spline tries to go through + every data point, it may show some oscillatory behaviour, specially, near + the merger. + + `rational_fit` is useful when the data is noisy. For example, when the + waveform modes are generated from numerical simulations. Rational fit + minimizes least square error making the overall trend smooth and less + oscillatory. However, it could be orders of magnitude slower compared to + `spline` method. Also, since it emphasizes to fit the overall trend, + it may suppress pathologies in the waveform which could be checked using + `spine` method. + + Therefore, we use `spline` as default with `rational_fit` as fall back. + This means that, if the extrema interpolant using `spline` shows + nonmonotonicity in its first derivative then we fall back to `rational_fit`. + + Default is `spline`. + + use_rational_fit_as_fallback : bool, default=True + Use rational fit for interpolant of omega at extrema when the + interpolant built using spline shows nonmonotonicity in its + first derivative. + + This is used only when `omega_gw_extrema_interpolation_method` is `spline`. + If `omega_gw_extrema_interpolation_method` is `rational_fit` and it has + no use. """ self.precessing = precessing # Get data necessary for eccentricity measurement @@ -323,6 +363,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, self.rational_fit_kwargs["verbose"] = True if self.debug_level >=1 else False self.debug_plots = self.extra_kwargs["debug_plots"] self.return_zero_if_small_ecc_failure = self.extra_kwargs["return_zero_if_small_ecc_failure"] + self.use_rational_fit_as_fallback = self.extra_kwargs["use_rational_fit_as_fallback"] # check if there are unrecognized keys in the dataDict self.recognized_dataDict_keys = self.get_recognized_dataDict_keys() for kw in dataDict.keys(): @@ -841,6 +882,8 @@ def get_default_extra_kwargs(self): "refine_extrema": False, "kwargs_for_fits_methods": {}, # Gets overriden in fits methods "return_zero_if_small_ecc_failure": False, + "use_rational_fit_as_fallback": True, + "omega_gw_extrema_interpolation_method": "spline" } return default_extra_kwargs @@ -1185,7 +1228,7 @@ def interp_extrema(self, extrema_type="pericenters"): " Can not create an interpolant.") def get_rat_fit(self, x, y): - """Get interpolant. + """Get rational fit. A wrapper of `utils.get_rational_fit` with check_kwargs=False. This is to make sure that the checking of kwargs is not performed everytime the @@ -1221,7 +1264,7 @@ def rational_fit(self, x, y): debug_level=self.debug_level, important=True) rat_fit = self.get_rat_fit(x, y) - omega = rational_fit(t) + omega = rat_fit(t) return rat_fit def rational_fit_extrema(self, extrema_type="pericenters"): @@ -1596,8 +1639,15 @@ def measure_ecc(self, tref_in=None, fref_in=None): "apocenters") # Build the interpolants of omega_gw at the extrema - self.omega_gw_pericenters_interp = self.interp_extrema("pericenters") - self.omega_gw_apocenters_interp = self.interp_extrema("apocenters") + if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": + self.omega_gw_pericenters_interp = self.interp_extrema("pericenters") + self.omega_gw_apocenters_interp = self.interp_extrema("apocenters") + elif self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": + self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") + self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") + else: + raise Exception("Unknown method for `omega_gw_extrema_interpolation_method`. " + "Must be one of `spline` or `rational_fit`.") self.t_pericenters = self.t[self.pericenters_location] self.t_apocenters = self.t[self.apocenters_location] @@ -1618,17 +1668,25 @@ def measure_ecc(self, tref_in=None, fref_in=None): self.dataDict["t"] <= self.tmax)] # Check monotonicity of the first derivative of omega_gw interpolant - # and build rational fit if necessary. - if (self.check_domega_dt(self.t_for_checks, - self.omega_gw_pericenters_interp(self.t_for_checks)) - or - self.check_domega_dt(self.t_for_checks, - self.omega_gw_apocenters_interp(self.t_for_checks))): - debug_message("The spline interpolant through extrema has non " - "monotonic time derivative. Trying rational fit instead.", - debug_level=self.debug_level, important=True) - self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") - self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") + # in case it is built with spine and fallback to build rational fit if necessary. + if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": + if (self.check_domega_dt(self.t_for_checks, + self.omega_gw_pericenters_interp(self.t_for_checks)) + or + self.check_domega_dt(self.t_for_checks, + self.omega_gw_apocenters_interp(self.t_for_checks))): + if self.use_rational_fit_as_fallback: + debug_message("The spline interpolant through extrema has non " + "monotonic time derivative. Trying rational fit instead.", + debug_level=self.debug_level, important=True) + self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") + self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") + else: + debug_message("The spline interpolant through extrema has non " + "monotonic time derivative. To avoid this, try rational " + "fit by setting `use_rational_fit_as_fallback` to `True` " + "in `extra_kwargs`.", + debug_level=self.debug_level, important=True) # Sanity checks # check that fref_out and tref_out are of the same length From 59eb2bef6245206e55c5f9042370ab3617a48a40 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Thu, 26 Sep 2024 12:45:54 +0530 Subject: [PATCH 04/34] use approximate orbits to guess degree for rational fit --- gw_eccentricity/eccDefinition.py | 30 +++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index a37fd5a9..1fc18390 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -298,7 +298,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, oscillatory. However, it could be orders of magnitude slower compared to `spline` method. Also, since it emphasizes to fit the overall trend, it may suppress pathologies in the waveform which could be checked using - `spine` method. + `spline` method. Therefore, we use `spline` as default with `rational_fit` as fall back. This means that, if the extrema interpolant using `spline` shows @@ -312,7 +312,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, first derivative. This is used only when `omega_gw_extrema_interpolation_method` is `spline`. - If `omega_gw_extrema_interpolation_method` is `rational_fit` and it has + If `omega_gw_extrema_interpolation_method` is `rational_fit` then it has no use. """ self.precessing = precessing @@ -1253,7 +1253,8 @@ def rational_fit(self, x, y): # assign degree # TODO: Use an optimal degree based on the wavefrom length self.rational_fit_kwargs["num_degree"],\ - self.rational_fit_kwargs["denom_degree"] = 2, 2 + self.rational_fit_kwargs["denom_degree"] \ + = self.get_optimal_degree_for_rational_fit() rat_fit = self.get_rat_fit(x, y) t = np.arange(x[0], x[-1], self.t[1] - self.t[0]) omega = rat_fit(t) @@ -1266,6 +1267,19 @@ def rational_fit(self, x, y): rat_fit = self.get_rat_fit(x, y) omega = rat_fit(t) return rat_fit + + def get_optimal_degree_for_rational_fit(self): + """Get optimal degree based on the approximate number of orbits.""" + approximate_num_orbits = ((self.phase_gw[-1] - self.phase_gw[0]) + / (4 * np.pi)) + if approximate_num_orbits <= 5: + return 1, 1 + elif (approximate_num_orbits > 5) and (approximate_num_orbits <= 20): + return 2, 2 + elif (approximate_num_orbits > 20) and (approximate_num_orbits <= 50): + return 3, 3 + else: + return 4, 4 def rational_fit_extrema(self, extrema_type="pericenters"): """Build rational fit through extrema. @@ -1683,10 +1697,12 @@ def measure_ecc(self, tref_in=None, fref_in=None): self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") else: debug_message("The spline interpolant through extrema has non " - "monotonic time derivative. To avoid this, try rational " - "fit by setting `use_rational_fit_as_fallback` to `True` " - "in `extra_kwargs`.", - debug_level=self.debug_level, important=True) + "monotonic time derivative. To avoid this, you may try one of " + "the following: \n" + " - Use rational fit as fall back option by setting" + " `use_rational_fit_as_fallback` to `True` in `extra_kwargs`\n" + " - Use rational fit as `omega_gw_extrema_interpolation_method`.", + debug_level=self.debug_level, important=True) # Sanity checks # check that fref_out and tref_out are of the same length From 37ec2b35ba7bae8e946eea1ec2a584c5c5290251 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Thu, 26 Sep 2024 14:47:21 +0530 Subject: [PATCH 05/34] use optimal degree function only user provided degree is None --- gw_eccentricity/eccDefinition.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 1fc18390..effdf6a2 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -1250,11 +1250,14 @@ def rational_fit(self, x, y): check for any nonmonotonicity. In case of nonmonocity detected, we lower the degree by 1 and repeat until the check passes successfully. """ - # assign degree - # TODO: Use an optimal degree based on the wavefrom length - self.rational_fit_kwargs["num_degree"],\ - self.rational_fit_kwargs["denom_degree"] \ - = self.get_optimal_degree_for_rational_fit() + # assign degree based on approximate number of orbits if user provided + # degree is None. + # TODO: Improve this code for optimal degree + if (self.rational_fit_kwargs["num_degree"] is None) or ( + self.rational_fit_kwargs["denom_degree"] is None): + self.rational_fit_kwargs["num_degree"],\ + self.rational_fit_kwargs["denom_degree"] \ + = self.get_optimal_degree_for_rational_fit() rat_fit = self.get_rat_fit(x, y) t = np.arange(x[0], x[-1], self.t[1] - self.t[0]) omega = rat_fit(t) From 546b406c318ef3b038feca6a92ba65cc99c2328e Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Thu, 26 Sep 2024 18:10:32 +0530 Subject: [PATCH 06/34] improve debug message --- gw_eccentricity/eccDefinition.py | 46 +++++++++++++++++++------------- 1 file changed, 28 insertions(+), 18 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index effdf6a2..2376c603 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -1264,7 +1264,9 @@ def rational_fit(self, x, y): while self.check_domega_dt(t, omega, 1.0): self.rational_fit_kwargs["num_degree"] -= 1 self.rational_fit_kwargs["denom_degree"] -= 1 - debug_message(f"degree is lowered to {self.rational_fit_kwargs['num_degree']}", + debug_message(f"Rational fit with degree {self.rational_fit_kwargs['num_degree'] + 1} " + " has non-monotonic time derivative. Lowering degree to " + f"{self.rational_fit_kwargs['num_degree']} and trying again.", debug_level=self.debug_level, important=True) rat_fit = self.get_rat_fit(x, y) @@ -1684,28 +1686,36 @@ def measure_ecc(self, tref_in=None, fref_in=None): np.logical_and(self.dataDict["t"] >= self.tmin, self.dataDict["t"] <= self.tmax)] - # Check monotonicity of the first derivative of omega_gw interpolant - # in case it is built with spine and fallback to build rational fit if necessary. + # Verify the monotonicity of the first derivative of the omega_gw interpolant. + # If a spline is used for interpolation (as specified by 'omega_gw_extrema_interpolation_method'), + # non-monotonicity may occur in the first derivative. + # If 'use_rational_fit_as_fallback' is set to True, the spline interpolant + # will be replaced with a rational fit to ensure monotonic behavior. if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": - if (self.check_domega_dt(self.t_for_checks, - self.omega_gw_pericenters_interp(self.t_for_checks)) - or - self.check_domega_dt(self.t_for_checks, - self.omega_gw_apocenters_interp(self.t_for_checks))): + # Check if the first derivative of omega_gw at pericenters or apocenters is non-monotonic + has_non_monotonicity = ( + self.check_domega_dt(self.t_for_checks, self.omega_gw_pericenters_interp(self.t_for_checks)) or + self.check_domega_dt(self.t_for_checks, self.omega_gw_apocenters_interp(self.t_for_checks)) + ) + + if has_non_monotonicity: if self.use_rational_fit_as_fallback: - debug_message("The spline interpolant through extrema has non " - "monotonic time derivative. Trying rational fit instead.", - debug_level=self.debug_level, important=True) + debug_message( + "Non-monotonic time derivative detected in the spline interpolant through extrema. " + "Switching to rational fit.", + debug_level=self.debug_level, important=True + ) + # Use rational fit for both pericenters and apocenters self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") else: - debug_message("The spline interpolant through extrema has non " - "monotonic time derivative. To avoid this, you may try one of " - "the following: \n" - " - Use rational fit as fall back option by setting" - " `use_rational_fit_as_fallback` to `True` in `extra_kwargs`\n" - " - Use rational fit as `omega_gw_extrema_interpolation_method`.", - debug_level=self.debug_level, important=True) + debug_message( + "Non-monotonic time derivative detected in the spline interpolant through extrema. " + "Consider the following options to avoid this: \n" + " - Set 'use_rational_fit_as_fallback' to True in 'extra_kwargs' to switch to rational fit.\n" + " - Use rational fit directly by setting 'omega_gw_extrema_interpolation_method' to 'rational_fit'.", + debug_level=self.debug_level, important=True + ) # Sanity checks # check that fref_out and tref_out are of the same length From 89fc1814fa796ce5c7fb281761e386e6109d3bcf Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Fri, 27 Sep 2024 12:35:55 +0530 Subject: [PATCH 07/34] more polishing --- gw_eccentricity/eccDefinition.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 2376c603..bbdeafac 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -1230,7 +1230,7 @@ def interp_extrema(self, extrema_type="pericenters"): def get_rat_fit(self, x, y): """Get rational fit. - A wrapper of `utils.get_rational_fit` with check_kwargs=False. This is + A wrapper of `utils.get_rational_fit` with check_kwargs=False. This is to make sure that the checking of kwargs is not performed everytime the rational fit function is called. Instead, the kwargs are checked once in the init and passed to the rational fit function without repeating @@ -1250,9 +1250,6 @@ def rational_fit(self, x, y): check for any nonmonotonicity. In case of nonmonocity detected, we lower the degree by 1 and repeat until the check passes successfully. """ - # assign degree based on approximate number of orbits if user provided - # degree is None. - # TODO: Improve this code for optimal degree if (self.rational_fit_kwargs["num_degree"] is None) or ( self.rational_fit_kwargs["denom_degree"] is None): self.rational_fit_kwargs["num_degree"],\ @@ -1260,8 +1257,7 @@ def rational_fit(self, x, y): = self.get_optimal_degree_for_rational_fit() rat_fit = self.get_rat_fit(x, y) t = np.arange(x[0], x[-1], self.t[1] - self.t[0]) - omega = rat_fit(t) - while self.check_domega_dt(t, omega, 1.0): + while self.check_domega_dt(t, rat_fit(t), 1.0): self.rational_fit_kwargs["num_degree"] -= 1 self.rational_fit_kwargs["denom_degree"] -= 1 debug_message(f"Rational fit with degree {self.rational_fit_kwargs['num_degree'] + 1} " @@ -1270,11 +1266,20 @@ def rational_fit(self, x, y): debug_level=self.debug_level, important=True) rat_fit = self.get_rat_fit(x, y) - omega = rat_fit(t) return rat_fit - + def get_optimal_degree_for_rational_fit(self): - """Get optimal degree based on the approximate number of orbits.""" + """Get optimal degree based on the approximate number of orbits. + + Assign degree based on approximate number of orbits if user provided + degree is None. The number of orbits is approximated by the change in + phase_gw divided by 4*pi. The degree of the rational fit is then chosen + based on this approximate number of orbits. The degree is increased as + the number of orbits increases. + """ + # TODO: Optimize this. + # assign degree based on approximate number of orbits if user provided + # degree is None. approximate_num_orbits = ((self.phase_gw[-1] - self.phase_gw[0]) / (4 * np.pi)) if approximate_num_orbits <= 5: From b17483cd0c62027a05acb8f01fdfbfca8e5c5be8 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Fri, 27 Sep 2024 13:01:47 +0530 Subject: [PATCH 08/34] improve omega interpolant documentation --- gw_eccentricity/eccDefinition.py | 56 +++++++++++++++++++------------- 1 file changed, 33 insertions(+), 23 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index bbdeafac..2aefb249 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -280,32 +280,42 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, USE THIS WITH CAUTION! "omega_gw_extrema_interpolation_method" : str, default="spline" - Method to use for building `omega_gw_pericenters_interp(t)` or - `omega_gw_apocenters_interp(t)`. Available options are: + 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.InterpolatedUnivariateSline`. - - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation` + - `spline`: Uses `scipy.interpolate.InterpolatedUnivariateSpline`. + - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation`. + + ### When to Use: + + - **`spline`** (default): + - 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. - `spline` is useful when the data is cleaner. For example, when the - waveform modes are generated using an waveform model like SEOB/TEOB. - It is also faster to build and evaluate. Since spline tries to go through - every data point, it may show some oscillatory behaviour, specially, near - the merger. - - `rational_fit` is useful when the data is noisy. For example, when the - waveform modes are generated from numerical simulations. Rational fit - minimizes least square error making the overall trend smooth and less - oscillatory. However, it could be orders of magnitude slower compared to - `spline` method. Also, since it emphasizes to fit the overall trend, - it may suppress pathologies in the waveform which could be checked using - `spline` method. - - Therefore, we use `spline` as default with `rational_fit` as fall back. - This means that, if the extrema interpolant using `spline` shows - nonmonotonicity in its first derivative then we fall back to `rational_fit`. + - **`rational_fit`**: + - More appropriate for noisy data, e.g., waveform modes from numerical + simulations. + - Minimizes least squares error, resulting in a smoother overall trend + with less oscillation. + - Significantly slower compared to the `spline` method. + - Can suppress pathologies in the waveform that might be visible with + `spline`. + + ### Fallback Behavior: - Default is `spline`. - + With `omega_gw_extrema_interpolation_method` set to `spline`, if + `use_rational_fit_as_fallback` is set to `True`, the method will + initially use `spline`. If the first derivative of the spline interpolant + exhibits non-monotonicity, the code will automatically fall back to the + `rational_fit` method. This ensures a more reliable fit when the spline + method shows undesirable behavior. + + Default value: `"spline"`. + use_rational_fit_as_fallback : bool, default=True Use rational fit for interpolant of omega at extrema when the interpolant built using spline shows nonmonotonicity in its From 92ddab01e525178f0967277a372b3c2ec3386659 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Fri, 27 Sep 2024 13:32:58 +0530 Subject: [PATCH 09/34] minor --- gw_eccentricity/eccDefinition.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 2aefb249..07131c60 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -279,7 +279,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, and mean anomaly to zero. USE THIS WITH CAUTION! - "omega_gw_extrema_interpolation_method" : str, default="spline" + omega_gw_extrema_interpolation_method : str, default="spline" Specifies the method used to build the interpolations for `omega_gw_pericenters_interp(t)` or `omega_gw_apocenters_interp(t)`. The available options are: From 30f737623f5245ce6fba5d700bcb3edef21b24b8 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Sat, 28 Sep 2024 11:34:31 +0530 Subject: [PATCH 10/34] copy documentation to gw_eccentricity --- gw_eccentricity/eccDefinition.py | 3 +- gw_eccentricity/gw_eccentricity.py | 50 ++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 1 deletion(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 07131c60..a7d5e647 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -178,7 +178,8 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, rational_fit_kwargs: dict Dictionary of arguments to be passed to the rational - fit function. + fit function. Defaults are set using + `utils.get_default_rational_fit_kwargs` extrema_finding_kwargs: dict Dictionary of arguments to be passed to the extrema finder, diff --git a/gw_eccentricity/gw_eccentricity.py b/gw_eccentricity/gw_eccentricity.py index 69372742..ffece41c 100644 --- a/gw_eccentricity/gw_eccentricity.py +++ b/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. @@ -411,6 +416,51 @@ def measure_eccentricity(tref_in=None, eccentricity is the cause, and set the returned eccentricity and mean anomaly to zero. USE THIS WITH CAUTION! + omega_gw_extrema_interpolation_method : str, default="spline" + 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`. + - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation`. + + ### When to Use: + + - **`spline`** (default): + - 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`**: + - More appropriate for noisy data, e.g., waveform modes from numerical + simulations. + - Minimizes least squares error, resulting in a smoother overall trend + with less oscillation. + - Significantly slower compared to the `spline` method. + - Can suppress pathologies in the waveform that might be visible with + `spline`. + + ### Fallback Behavior: + + With `omega_gw_extrema_interpolation_method` set to `spline`, if + `use_rational_fit_as_fallback` is set to `True`, the method will + initially use `spline`. If the first derivative of the spline interpolant + exhibits non-monotonicity, the code will automatically fall back to the + `rational_fit` method. This ensures a more reliable fit when the spline + method shows undesirable behavior. + + Default value: `"spline"`. + + use_rational_fit_as_fallback : bool, default=True + Use rational fit for interpolant of omega at extrema when the + interpolant built using spline shows nonmonotonicity in its + first derivative. + + This is used only when `omega_gw_extrema_interpolation_method` is `spline`. + If `omega_gw_extrema_interpolation_method` is `rational_fit` then it has + no use. Returns ------- From bf5619ccb281ebff50c39b87d232c0ddd726303f Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Sun, 29 Sep 2024 15:44:24 +0530 Subject: [PATCH 11/34] minor --- gw_eccentricity/eccDefinition.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index a7d5e647..1abf9eb9 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -371,7 +371,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, self.available_averaging_methods \ = self.get_available_omega_gw_averaging_methods() self.debug_level = self.extra_kwargs["debug_level"] - self.rational_fit_kwargs["verbose"] = True if self.debug_level >=1 else False + self.rational_fit_kwargs["verbose"] = True if self.debug_level >= 1 else False self.debug_plots = self.extra_kwargs["debug_plots"] self.return_zero_if_small_ecc_failure = self.extra_kwargs["return_zero_if_small_ecc_failure"] self.use_rational_fit_as_fallback = self.extra_kwargs["use_rational_fit_as_fallback"] From 326b10a1cb3894df69d47ae182ace3e66cfcd1ab Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Mon, 14 Oct 2024 18:50:23 +0530 Subject: [PATCH 12/34] raise exception if degree=1 is attempted to be lowered --- gw_eccentricity/eccDefinition.py | 27 ++++++++++++++++++++------- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index fbfbf080..bfc77265 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -1279,13 +1279,26 @@ def rational_fit(self, x, y): rat_fit = self.get_rat_fit(x, y) t = np.arange(x[0], x[-1], self.t[1] - self.t[0]) while self.check_domega_dt(t, rat_fit(t), 1.0): - self.rational_fit_kwargs["num_degree"] -= 1 - self.rational_fit_kwargs["denom_degree"] -= 1 - debug_message(f"Rational fit with degree {self.rational_fit_kwargs['num_degree'] + 1} " - " has non-monotonic time derivative. Lowering degree to " - f"{self.rational_fit_kwargs['num_degree']} and trying again.", - debug_level=self.debug_level, - important=True) + if self.rational_fit_kwargs["num_degree"] == 1: + raise Exception("`num_degree` for rational fit is already = 1. " + "Can not be lowered further.") + else: + self.rational_fit_kwargs["num_degree"] -= 1 + debug_message(f"Rational fit with `num_degree` {self.rational_fit_kwargs['num_degree'] + 1} " + " has non-monotonic time derivative. Lowering degree to " + f"{self.rational_fit_kwargs['num_degree']} and trying again.", + debug_level=self.debug_level, + important=True) + if self.rational_fit_kwargs["denom_degree"] == 1: + raise Exception("`denmm_degree` for rational fit is already = 1. " + "Can not be lowered further.") + else: + self.rational_fit_kwargs["denom_degree"] -= 1 + debug_message(f"Rational fit with `denom_degree` {self.rational_fit_kwargs['num_degree'] + 1} " + " has non-monotonic time derivative. Lowering degree to " + f"{self.rational_fit_kwargs['num_degree']} and trying again.", + debug_level=self.debug_level, + important=True) rat_fit = self.get_rat_fit(x, y) return rat_fit From 42a67555a1076bede1be4c9238e6018fdb208cb0 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Thu, 17 Oct 2024 13:26:35 +0900 Subject: [PATCH 13/34] fix nonmonotic egw by increasing rational fit degree --- gw_eccentricity/eccDefinition.py | 129 ++++++++++++++++++++----------- 1 file changed, 83 insertions(+), 46 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index bfc77265..787c382a 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -1290,7 +1290,7 @@ def rational_fit(self, x, y): debug_level=self.debug_level, important=True) if self.rational_fit_kwargs["denom_degree"] == 1: - raise Exception("`denmm_degree` for rational fit is already = 1. " + raise Exception("`denom_degree` for rational fit is already = 1. " "Can not be lowered further.") else: self.rational_fit_kwargs["denom_degree"] -= 1 @@ -1696,17 +1696,6 @@ def measure_ecc(self, tref_in=None, fref_in=None): = self.check_extrema_separation(self.apocenters_location, "apocenters") - # Build the interpolants of omega_gw at the extrema - if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": - self.omega_gw_pericenters_interp = self.interp_extrema("pericenters") - self.omega_gw_apocenters_interp = self.interp_extrema("apocenters") - elif self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": - self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") - self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") - else: - raise Exception("Unknown method for `omega_gw_extrema_interpolation_method`. " - "Must be one of `spline` or `rational_fit`.") - self.t_pericenters = self.t[self.pericenters_location] self.t_apocenters = self.t[self.apocenters_location] self.tmax = min(self.t_pericenters[-1], self.t_apocenters[-1]) @@ -1725,37 +1714,6 @@ def measure_ecc(self, tref_in=None, fref_in=None): np.logical_and(self.dataDict["t"] >= self.tmin, self.dataDict["t"] <= self.tmax)] - # Verify the monotonicity of the first derivative of the omega_gw interpolant. - # If a spline is used for interpolation (as specified by 'omega_gw_extrema_interpolation_method'), - # non-monotonicity may occur in the first derivative. - # If 'use_rational_fit_as_fallback' is set to True, the spline interpolant - # will be replaced with a rational fit to ensure monotonic behavior. - if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": - # Check if the first derivative of omega_gw at pericenters or apocenters is non-monotonic - has_non_monotonicity = ( - self.check_domega_dt(self.t_for_checks, self.omega_gw_pericenters_interp(self.t_for_checks)) or - self.check_domega_dt(self.t_for_checks, self.omega_gw_apocenters_interp(self.t_for_checks)) - ) - - if has_non_monotonicity: - if self.use_rational_fit_as_fallback: - debug_message( - "Non-monotonic time derivative detected in the spline interpolant through extrema. " - "Switching to rational fit.", - debug_level=self.debug_level, important=True - ) - # Use rational fit for both pericenters and apocenters - self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") - self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") - else: - debug_message( - "Non-monotonic time derivative detected in the spline interpolant through extrema. " - "Consider the following options to avoid this: \n" - " - Set 'use_rational_fit_as_fallback' to True in 'extra_kwargs' to switch to rational fit.\n" - " - Use rational fit directly by setting 'omega_gw_extrema_interpolation_method' to 'rational_fit'.", - debug_level=self.debug_level, important=True - ) - # Sanity checks # check that fref_out and tref_out are of the same length if self.domain == "frequency": @@ -1791,19 +1749,45 @@ def measure_ecc(self, tref_in=None, fref_in=None): or self.tref_out[-1] > self.t_pericenters[-1]: raise Exception("Reference time must be within two pericenters.") + # Build omega_gw extrema interpolants + self.build_omega_gw_extrema_interpolants() + # compute eccentricity at self.tref_out self.eccentricity = self.compute_eccentricity(self.tref_out) # Compute mean anomaly at tref_out self.mean_anomaly = self.compute_mean_anomaly(self.tref_out) + # check if eccentricity is nonmonotonic and try to fix. + # If eccentricity is nonmonotonic and `omega_gw_extrema_interpolation_method` + # is `rational_fit`, sometimes increasing the degree might help. + # Note that increasing the degree may also break the monotonicity of the + # omega_gw extrema interpolants (high degree my lead to divergences) + # which will automatically decrease the degree. + # Therefore, we check the degree before and after trying with increased degree. + # If the degree is the same as before, it implies that increasing degree did not + # help and we stop attempting to increase the degree any further. + if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": + while self.check_monotonicity_and_convexity()["monotonic"] == False: + # store old degrees to compare later + num_degree_old = np.copy(self.rational_fit_kwargs["num_degree"]) + debug_message("Trying to fix nonmonotic egw evolution by increasing rational fit " + "degree by 1.", self.debug_level) + self.rational_fit_kwargs["num_degree"] += 1 + self.rational_fit_kwargs["denom_degree"] += 1 + self.build_omega_gw_extrema_interpolants() + if self.rational_fit_kwargs["num_degree"] == num_degree_old: + break + self.eccentricity = self.compute_eccentricity(self.tref_out) + # update ecc for checks + self.ecc_for_checks = self.compute_eccentricity(self.t_for_checks) + else: + self.check_monotonicity_and_convexity() + # check if eccentricity is positive if any(self.eccentricity < 0): debug_message("Encountered negative eccentricity.", self.debug_level, point_to_verbose_output=True) - # check if eccentricity is monotonic and convex - self.check_monotonicity_and_convexity() - if self.debug_plots: # make a plot for diagnostics fig, axes = self.make_diagnostic_plots() @@ -1812,6 +1796,49 @@ def measure_ecc(self, tref_in=None, fref_in=None): # return measured eccentricity, mean anomaly and reference time or # frequency where these are measured. return self.make_return_dict_for_eccentricity_and_mean_anomaly() + + def build_omega_gw_extrema_interpolants(self): + # Build the interpolants of omega_gw at the extrema + if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": + self.omega_gw_pericenters_interp = self.interp_extrema("pericenters") + self.omega_gw_apocenters_interp = self.interp_extrema("apocenters") + elif self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": + self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") + self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") + else: + raise Exception("Unknown method for `omega_gw_extrema_interpolation_method`. " + "Must be one of `spline` or `rational_fit`.") + + # Verify the monotonicity of the first derivative of the omega_gw interpolant. + # If a spline is used for interpolation (as specified by 'omega_gw_extrema_interpolation_method'), + # non-monotonicity may occur in the first derivative. + # If 'use_rational_fit_as_fallback' is set to True, the spline interpolant + # will be replaced with a rational fit to ensure monotonic behavior. + if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": + # Check if the first derivative of omega_gw at pericenters or apocenters is non-monotonic + has_non_monotonicity = ( + self.check_domega_dt(self.t_for_checks, self.omega_gw_pericenters_interp(self.t_for_checks)) or + self.check_domega_dt(self.t_for_checks, self.omega_gw_apocenters_interp(self.t_for_checks)) + ) + + if has_non_monotonicity: + if self.use_rational_fit_as_fallback: + debug_message( + "Non-monotonic time derivative detected in the spline interpolant through extrema. " + "Switching to rational fit.", + debug_level=self.debug_level, important=True + ) + # Use rational fit for both pericenters and apocenters + self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") + self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") + else: + debug_message( + "Non-monotonic time derivative detected in the spline interpolant through extrema. " + "Consider the following options to avoid this: \n" + " - Set 'use_rational_fit_as_fallback' to True in 'extra_kwargs' to switch to rational fit.\n" + " - Use rational fit directly by setting 'omega_gw_extrema_interpolation_method' to 'rational_fit'.", + debug_level=self.debug_level, important=True + ) def set_eccentricity_and_mean_anomaly_to_zero(self): """Set eccentricity and mean_anomaly to zero.""" @@ -2119,6 +2146,9 @@ def check_monotonicity_and_convexity(self, self.decc_dt_for_checks = self.derivative_of_eccentricity( self.t_for_checks, n=1) + # Create an empty dictionary to store the check results + check_dict = {} + # Is ecc(t) a monotonically decreasing function? if any(self.decc_dt_for_checks > 0): idx = np.where(self.decc_dt_for_checks > 0)[0] @@ -2127,6 +2157,9 @@ def check_monotonicity_and_convexity(self, f"{'at' if len(idx) == 1 else 'in the range'} {range}") debug_message(message, self.debug_level, point_to_verbose_output=True) + check_dict.update({"monotonic": False}) + else: + check_dict.update({"monotonic": True}) # Is ecc(t) a convex function? That is, is the second # derivative always negative? @@ -2142,6 +2175,10 @@ def check_monotonicity_and_convexity(self, debug_message(f"{message} expected to be always negative", self.debug_level, point_to_verbose_output=True) + check_dict.update({"convex": False}) + else: + check_dict.update({"convex": True}) + return check_dict def get_range_from_indices(self, indices, times): """Get the range of time from indices for gives times array.""" From 4a5e4c961af4936a21e4d22f87a1b9124ad66ae6 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Thu, 17 Oct 2024 10:29:19 +0530 Subject: [PATCH 14/34] add informative message --- gw_eccentricity/eccDefinition.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 787c382a..334bd216 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -1770,8 +1770,9 @@ def measure_ecc(self, tref_in=None, fref_in=None): while self.check_monotonicity_and_convexity()["monotonic"] == False: # store old degrees to compare later num_degree_old = np.copy(self.rational_fit_kwargs["num_degree"]) - debug_message("Trying to fix nonmonotic egw evolution by increasing rational fit " - "degree by 1.", self.debug_level) + debug_message("Trying to fix nonmonotic egw evolution (current degree = " + f"{num_degree_old}) by increasing rational fit " + "degree by 1.", self.debug_level, important=False) self.rational_fit_kwargs["num_degree"] += 1 self.rational_fit_kwargs["denom_degree"] += 1 self.build_omega_gw_extrema_interpolants() From 95e60e17a115a17cf211a91d6a6b0fc5f9f8bcca Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Fri, 18 Oct 2024 07:44:14 +0530 Subject: [PATCH 15/34] turn on debug --- gw_eccentricity/eccDefinition.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 334bd216..5755baad 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -1772,7 +1772,7 @@ def measure_ecc(self, tref_in=None, fref_in=None): num_degree_old = np.copy(self.rational_fit_kwargs["num_degree"]) debug_message("Trying to fix nonmonotic egw evolution (current degree = " f"{num_degree_old}) by increasing rational fit " - "degree by 1.", self.debug_level, important=False) + "degree by 1.", self.debug_level, important=True) self.rational_fit_kwargs["num_degree"] += 1 self.rational_fit_kwargs["denom_degree"] += 1 self.build_omega_gw_extrema_interpolants() From b4b1ab364fdf81298a6653800dcd1ff8b7789239 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Fri, 18 Oct 2024 21:54:04 +0900 Subject: [PATCH 16/34] refactor egw checking and fixing with rational fits --- gw_eccentricity/eccDefinition.py | 95 ++++++++++++++++++++++---------- 1 file changed, 66 insertions(+), 29 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 5755baad..73bd5e26 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -372,6 +372,8 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, = self.get_available_omega_gw_averaging_methods() self.debug_level = self.extra_kwargs["debug_level"] self.rational_fit_kwargs["verbose"] = True if self.debug_level >= 1 else False + # keep history of rational fit degree and nonmonotonicity of the corresponding fits + self.rational_fit_nonmonotonicity_history = {} self.debug_plots = self.extra_kwargs["debug_plots"] self.return_zero_if_small_ecc_failure = self.extra_kwargs["return_zero_if_small_ecc_failure"] self.use_rational_fit_as_fallback = self.extra_kwargs["use_rational_fit_as_fallback"] @@ -1303,27 +1305,31 @@ def rational_fit(self, x, y): return rat_fit def get_optimal_degree_for_rational_fit(self): - """Get optimal degree based on the approximate number of orbits. + """Get optimal degree based on the number of extrema. - Assign degree based on approximate number of orbits if user provided - degree is None. The number of orbits is approximated by the change in - phase_gw divided by 4*pi. The degree of the rational fit is then chosen - based on this approximate number of orbits. The degree is increased as - the number of orbits increases. + Assign degree based on number of extrema found. The degree is increased as + the number of extrema increases. """ # TODO: Optimize this. - # assign degree based on approximate number of orbits if user provided + # assign degree based on the number of extrema if user provided # degree is None. - approximate_num_orbits = ((self.phase_gw[-1] - self.phase_gw[0]) - / (4 * np.pi)) + approximate_num_orbits = max(len(self.pericenters_location), + len(self.apocenters_location)) if approximate_num_orbits <= 5: - return 1, 1 + num_degree, denom_degree = 1, 1 elif (approximate_num_orbits > 5) and (approximate_num_orbits <= 20): - return 2, 2 + num_degree, denom_degree = 2, 2 elif (approximate_num_orbits > 20) and (approximate_num_orbits <= 50): - return 3, 3 + num_degree, denom_degree = 3, 3 + elif (approximate_num_orbits > 50) and (approximate_num_orbits <= 100): + num_degree, denom_degree = 4, 4 + elif (approximate_num_orbits > 100) and (approximate_num_orbits <= 150): + num_degree, denom_degree = 5, 5 + elif (approximate_num_orbits > 150) and (approximate_num_orbits <= 200): + num_degree, denom_degree = 6, 6 else: - return 4, 4 + num_degree, denom_degree = 5 + int(np.log10(approximate_num_orbits)), 5 + int(np.log10(approximate_num_orbits)) + return num_degree, denom_degree def rational_fit_extrema(self, extrema_type="pericenters"): """Build rational fit through extrema. @@ -1359,7 +1365,12 @@ def check_domega_dt(self, t, omega, tol=1.0): should be monotonically increasing. """ domega_dt = np.gradient(omega, t[1] - t[0]) - return any(domega_dt[1:]/domega_dt[:-1] < tol) + is_nonmonotonic = any(domega_dt[1:]/domega_dt[:-1] < tol) + # update history + if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": + self.rational_fit_nonmonotonicity_history.update({ + f"{self.rational_fit_kwargs['num_degree']}": is_nonmonotonic}) + return is_nonmonotonic def check_num_extrema(self, extrema, extrema_type="extrema"): """Check number of extrema. @@ -1751,7 +1762,6 @@ def measure_ecc(self, tref_in=None, fref_in=None): # Build omega_gw extrema interpolants self.build_omega_gw_extrema_interpolants() - # compute eccentricity at self.tref_out self.eccentricity = self.compute_eccentricity(self.tref_out) # Compute mean anomaly at tref_out @@ -1767,20 +1777,7 @@ def measure_ecc(self, tref_in=None, fref_in=None): # If the degree is the same as before, it implies that increasing degree did not # help and we stop attempting to increase the degree any further. if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": - while self.check_monotonicity_and_convexity()["monotonic"] == False: - # store old degrees to compare later - num_degree_old = np.copy(self.rational_fit_kwargs["num_degree"]) - debug_message("Trying to fix nonmonotic egw evolution (current degree = " - f"{num_degree_old}) by increasing rational fit " - "degree by 1.", self.debug_level, important=True) - self.rational_fit_kwargs["num_degree"] += 1 - self.rational_fit_kwargs["denom_degree"] += 1 - self.build_omega_gw_extrema_interpolants() - if self.rational_fit_kwargs["num_degree"] == num_degree_old: - break - self.eccentricity = self.compute_eccentricity(self.tref_out) - # update ecc for checks - self.ecc_for_checks = self.compute_eccentricity(self.t_for_checks) + self.check_egw_and_try_to_fix_if_nonmonotonic() else: self.check_monotonicity_and_convexity() @@ -1798,6 +1795,46 @@ def measure_ecc(self, tref_in=None, fref_in=None): # frequency where these are measured. return self.make_return_dict_for_eccentricity_and_mean_anomaly() + def check_egw_and_try_to_fix_if_nonmonotonic(self): + while self.check_monotonicity_and_convexity()["monotonic"] == False: + # store old degrees to compare later + num_degree_old = np.copy(self.rational_fit_kwargs["num_degree"]) + denom_degree_old = np.copy(self.rational_fit_kwargs["denom_degree"]) + # Increase degree by one + debug_message("Attempting to resolve nonmonotic egw evolution " + f"(current degree: {num_degree_old}) by increasing rational fit " + "degree by 1.", self.debug_level, important=True) + self.rational_fit_kwargs["num_degree"] += 1 + self.rational_fit_kwargs["denom_degree"] += 1 + # check if this degree is already tried and corresponding monotonicity + # If it was already tried and resulted in nonmonotic fits, then we abandon + # this attempt and set the rational fit degrees to its previous values. + if not self.check_whether_to_try_new_degree(self.rational_fit_kwargs["num_degree"]): + # reset the correct value + self.rational_fit_kwargs["num_degree"] = int(num_degree_old) + self.rational_fit_kwargs["denom_degree"] = int(denom_degree_old) + debug_message("Final rational fits were built with " + f"`num_degree`={self.rational_fit_kwargs['num_degree']} and " + f"`denom_degree`={self.rational_fit_kwargs['denom_degree']}.", + debug_level=self.debug_level, important=True) + break + # build interpolants with updated degree + self.build_omega_gw_extrema_interpolants() + # compute eccentricity using updated interpolants + self.eccentricity = self.compute_eccentricity(self.tref_out) + # update ecc for checks + self.ecc_for_checks = self.compute_eccentricity(self.t_for_checks) + + def check_whether_to_try_new_degree(self, new_degree): + if f"{new_degree}" in self.rational_fit_nonmonotonicity_history: + if self.rational_fit_nonmonotonicity_history[f"{new_degree}"]: + debug_message(f"Rational fit was already built with degree: {new_degree} and " + f"was found to be nonmonotonic. Abandoning attempt with degree: {new_degree}", + debug_level=self.debug_level, important=True) + return not self.rational_fit_nonmonotonicity_history[f"{new_degree}"] + else: + return True + def build_omega_gw_extrema_interpolants(self): # Build the interpolants of omega_gw at the extrema if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": From 5a0569000de1f436d0e2bc1ca9780764faeb243e Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Fri, 18 Oct 2024 22:04:05 +0900 Subject: [PATCH 17/34] fix typos --- gw_eccentricity/eccDefinition.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 73bd5e26..6f1989a7 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -1282,22 +1282,22 @@ def rational_fit(self, x, y): t = np.arange(x[0], x[-1], self.t[1] - self.t[0]) while self.check_domega_dt(t, rat_fit(t), 1.0): if self.rational_fit_kwargs["num_degree"] == 1: - raise Exception("`num_degree` for rational fit is already = 1. " + raise Exception("Current `num_degree` for rational fit is already 1. " "Can not be lowered further.") else: self.rational_fit_kwargs["num_degree"] -= 1 - debug_message(f"Rational fit with `num_degree` {self.rational_fit_kwargs['num_degree'] + 1} " - " has non-monotonic time derivative. Lowering degree to " + debug_message(f"Rational fit with `num_degree`={self.rational_fit_kwargs['num_degree'] + 1} " + " has nonmonotonic time derivative. Lowering degree to " f"{self.rational_fit_kwargs['num_degree']} and trying again.", debug_level=self.debug_level, important=True) if self.rational_fit_kwargs["denom_degree"] == 1: - raise Exception("`denom_degree` for rational fit is already = 1. " + raise Exception("Current `denom_degree` for rational fit is already 1. " "Can not be lowered further.") else: self.rational_fit_kwargs["denom_degree"] -= 1 - debug_message(f"Rational fit with `denom_degree` {self.rational_fit_kwargs['num_degree'] + 1} " - " has non-monotonic time derivative. Lowering degree to " + debug_message(f"Rational fit with `denom_degree`={self.rational_fit_kwargs['num_degree'] + 1} " + " has nonmonotonic time derivative. Lowering degree to " f"{self.rational_fit_kwargs['num_degree']} and trying again.", debug_level=self.debug_level, important=True) @@ -1802,13 +1802,13 @@ def check_egw_and_try_to_fix_if_nonmonotonic(self): denom_degree_old = np.copy(self.rational_fit_kwargs["denom_degree"]) # Increase degree by one debug_message("Attempting to resolve nonmonotic egw evolution " - f"(current degree: {num_degree_old}) by increasing rational fit " + f"(current degree={num_degree_old}) by increasing rational fit " "degree by 1.", self.debug_level, important=True) self.rational_fit_kwargs["num_degree"] += 1 self.rational_fit_kwargs["denom_degree"] += 1 # check if this degree is already tried and corresponding monotonicity - # If it was already tried and resulted in nonmonotic fits, then we abandon - # this attempt and set the rational fit degrees to its previous values. + # If it was already tried and resulted in nonmonotonic fits, then we abandon + # this attempt and set the rational fit degrees to their previous values. if not self.check_whether_to_try_new_degree(self.rational_fit_kwargs["num_degree"]): # reset the correct value self.rational_fit_kwargs["num_degree"] = int(num_degree_old) @@ -1828,7 +1828,7 @@ def check_egw_and_try_to_fix_if_nonmonotonic(self): def check_whether_to_try_new_degree(self, new_degree): if f"{new_degree}" in self.rational_fit_nonmonotonicity_history: if self.rational_fit_nonmonotonicity_history[f"{new_degree}"]: - debug_message(f"Rational fit was already built with degree: {new_degree} and " + debug_message(f"Rational fit was already built with degree={new_degree} and " f"was found to be nonmonotonic. Abandoning attempt with degree: {new_degree}", debug_level=self.debug_level, important=True) return not self.rational_fit_nonmonotonicity_history[f"{new_degree}"] From 83b7d08ac9198721aebb1e62908388d2d5ee37c4 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Fri, 25 Oct 2024 16:33:49 +0530 Subject: [PATCH 18/34] add variable to store final fit degrees --- gw_eccentricity/eccDefinition.py | 182 +++++++++++++++---------------- 1 file changed, 85 insertions(+), 97 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 6f1989a7..82f52dae 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -20,7 +20,6 @@ from .plot_settings import use_fancy_plotsettings, colorsDict, labelsDict from .plot_settings import figWidthsTwoColDict, figHeightsDict - class eccDefinition: """Base class to define eccentricity for given waveform data dictionary.""" @@ -373,7 +372,9 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, self.debug_level = self.extra_kwargs["debug_level"] self.rational_fit_kwargs["verbose"] = True if self.debug_level >= 1 else False # keep history of rational fit degree and nonmonotonicity of the corresponding fits - self.rational_fit_nonmonotonicity_history = {} + self.rational_fit_nonmonotonicity_history = {"pericenters": {}, "apocenters": {}} + # Update the degrees used to construct the final fits + self.rational_fit_degrees = {"pericenters": None, "apocenters": None} self.debug_plots = self.extra_kwargs["debug_plots"] self.return_zero_if_small_ecc_failure = self.extra_kwargs["return_zero_if_small_ecc_failure"] self.use_rational_fit_as_fallback = self.extra_kwargs["use_rational_fit_as_fallback"] @@ -1263,49 +1264,83 @@ def get_rat_fit(self, x, y): rational_fit_kwargs=self.rational_fit_kwargs, check_kwargs=False) - def rational_fit(self, x, y): + def rational_fit(self, x, y, description=None): """Get rational fit with adaptive numerator and denominator degree. - This function ensures that the rational fit we obtain is using the - optimal degree for the numerator and the denominator. We start with an - initial estimate of what these degrees should be based on the length of - the waveform. We check the first derivative of the resultant fit to - check for any nonmonotonicity. In case of nonmonocity detected, we - lower the degree by 1 and repeat until the check passes successfully. + This function ensures that the rational fit uses the optimal degree for + the numerator and denominator. It first checks for nonmonotonicity in the + fit's derivative and lowers the degrees if necessary. Only if the initial + degree was monotonic (no reduction needed), it attempts to increase the + degrees to find a higher-degree fit. If increasing the degrees causes + nonmonotonicity, it reverts to the previous valid monotonic fit. """ - if (self.rational_fit_kwargs["num_degree"] is None) or ( - self.rational_fit_kwargs["denom_degree"] is None): - self.rational_fit_kwargs["num_degree"],\ - self.rational_fit_kwargs["denom_degree"] \ - = self.get_optimal_degree_for_rational_fit() + # make sure that description is not None + if description is None: + raise Exception("Please provide a description. `description` can not be None.") + # Set initial degrees if not already specified + if not (self.rational_fit_kwargs["num_degree"] and self.rational_fit_kwargs["denom_degree"]): + self.rational_fit_kwargs["num_degree"], self.rational_fit_kwargs["denom_degree"] \ + = self.get_approximate_degree_for_rational_fit() + rat_fit = self.get_rat_fit(x, y) t = np.arange(x[0], x[-1], self.t[1] - self.t[0]) - while self.check_domega_dt(t, rat_fit(t), 1.0): - if self.rational_fit_kwargs["num_degree"] == 1: - raise Exception("Current `num_degree` for rational fit is already 1. " - "Can not be lowered further.") - else: + old_num_degree = self.rational_fit_kwargs["num_degree"] + old_denom_degree = self.rational_fit_kwargs["denom_degree"] + + degrees_were_reduced = False + + # Check for nonmonotonicity and lower degrees if needed + while self.check_if_domega_dt_is_nonmonotonic(t, rat_fit(t), 1.0, description): + degrees_were_reduced = True # Mark that reduction occurred + if self.rational_fit_kwargs["num_degree"] > 1: self.rational_fit_kwargs["num_degree"] -= 1 - debug_message(f"Rational fit with `num_degree`={self.rational_fit_kwargs['num_degree'] + 1} " - " has nonmonotonic time derivative. Lowering degree to " - f"{self.rational_fit_kwargs['num_degree']} and trying again.", - debug_level=self.debug_level, - important=True) - if self.rational_fit_kwargs["denom_degree"] == 1: - raise Exception("Current `denom_degree` for rational fit is already 1. " - "Can not be lowered further.") - else: + if self.rational_fit_kwargs["denom_degree"] > 1: self.rational_fit_kwargs["denom_degree"] -= 1 - debug_message(f"Rational fit with `denom_degree`={self.rational_fit_kwargs['num_degree'] + 1} " - " has nonmonotonic time derivative. Lowering degree to " - f"{self.rational_fit_kwargs['num_degree']} and trying again.", - debug_level=self.debug_level, - important=True) + if self.rational_fit_kwargs["num_degree"] == 1 and self.rational_fit_kwargs["denom_degree"] == 1: + raise Exception("Both numerator and denominator degrees cannot be lowered further.") + + debug_message(f"Lowering degrees to num_degree={self.rational_fit_kwargs['num_degree']}, " + f"denom_degree={self.rational_fit_kwargs['denom_degree']} and retrying.", + debug_level=self.debug_level, important=True) + rat_fit = self.get_rat_fit(x, y) + # If no degrees were reduced, try increasing the degree for better fit + if not degrees_were_reduced: + last_monotonic_rat_fit = rat_fit # Track last monotonic fit + # last_monotonic_rat_fit_vals = rat_fit(t) + last_monotonic_num_degree = old_num_degree + last_monotonic_denom_degree = old_denom_degree + + while not self.check_if_domega_dt_is_nonmonotonic(t, rat_fit(t), 1.0, description): + # Increase the degrees for both numerator and denominator + new_num_degree = self.rational_fit_kwargs["num_degree"] + 1 + new_denom_degree = self.rational_fit_kwargs["denom_degree"] + 1 + self.rational_fit_kwargs["num_degree"], self.rational_fit_kwargs["denom_degree"] \ + = new_num_degree, new_denom_degree + new_rat_fit = self.get_rat_fit(x, y) + if self.check_if_domega_dt_is_nonmonotonic(t, new_rat_fit(t), 1.0, description): + # Revert to previous fit and degrees if nonmonotonicity is detected + debug_message(f"Increasing degrees caused nonmonotonicity. Reverting to " + f"last monotonic fit with num_degree={last_monotonic_num_degree} " + f"and denom_degree={last_monotonic_denom_degree}.", + debug_level=self.debug_level, important=True) + self.rational_fit_kwargs["num_degree"] = last_monotonic_num_degree + self.rational_fit_kwargs["denom_degree"] = last_monotonic_denom_degree + rat_fit = last_monotonic_rat_fit + break + + last_monotonic_rat_fit = new_rat_fit + last_monotonic_num_degree = new_num_degree + last_monotonic_denom_degree = new_denom_degree + rat_fit = new_rat_fit + + # update final degrees used to build the fit + self.rational_fit_degrees[description] = (self.rational_fit_kwargs["num_degree"], + self.rational_fit_kwargs["denom_degree"]) return rat_fit - def get_optimal_degree_for_rational_fit(self): - """Get optimal degree based on the number of extrema. + def get_approximate_degree_for_rational_fit(self): + """Get approximate degree based on the number of extrema. Assign degree based on number of extrema found. The degree is increased as the number of extrema increases. @@ -1352,14 +1387,15 @@ def rational_fit_extrema(self, extrema_type="pericenters"): "'pericenrers' or 'apocenters'.") if len(extrema) >= 2: return self.rational_fit(self.t[extrema], - self.omega_gw[extrema]) + self.omega_gw[extrema], + description=extrema_type) else: raise Exception( f"Sufficient number of {extrema_type} are not found." " Can not create an interpolant.") - def check_domega_dt(self, t, omega, tol=1.0): - """Check first derivative of omega. + def check_if_domega_dt_is_nonmonotonic(self, t, omega, tol=1.0, description=None): + """Check if the first derivative of omega is nonmonotonic The first derivative of interpolant/fit of omega at the extrema should be monotonically increasing. @@ -1368,8 +1404,9 @@ def check_domega_dt(self, t, omega, tol=1.0): is_nonmonotonic = any(domega_dt[1:]/domega_dt[:-1] < tol) # update history if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": - self.rational_fit_nonmonotonicity_history.update({ - f"{self.rational_fit_kwargs['num_degree']}": is_nonmonotonic}) + if description in self.rational_fit_nonmonotonicity_history: + self.rational_fit_nonmonotonicity_history[description].update({ + (self.rational_fit_kwargs['num_degree'], self.rational_fit_kwargs['num_degree']): is_nonmonotonic}) return is_nonmonotonic def check_num_extrema(self, extrema, extrema_type="extrema"): @@ -1767,19 +1804,8 @@ def measure_ecc(self, tref_in=None, fref_in=None): # Compute mean anomaly at tref_out self.mean_anomaly = self.compute_mean_anomaly(self.tref_out) - # check if eccentricity is nonmonotonic and try to fix. - # If eccentricity is nonmonotonic and `omega_gw_extrema_interpolation_method` - # is `rational_fit`, sometimes increasing the degree might help. - # Note that increasing the degree may also break the monotonicity of the - # omega_gw extrema interpolants (high degree my lead to divergences) - # which will automatically decrease the degree. - # Therefore, we check the degree before and after trying with increased degree. - # If the degree is the same as before, it implies that increasing degree did not - # help and we stop attempting to increase the degree any further. - if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": - self.check_egw_and_try_to_fix_if_nonmonotonic() - else: - self.check_monotonicity_and_convexity() + # check if eccentricity is nonmonotonic + self.check_monotonicity_and_convexity() # check if eccentricity is positive if any(self.eccentricity < 0): @@ -1795,54 +1821,14 @@ def measure_ecc(self, tref_in=None, fref_in=None): # frequency where these are measured. return self.make_return_dict_for_eccentricity_and_mean_anomaly() - def check_egw_and_try_to_fix_if_nonmonotonic(self): - while self.check_monotonicity_and_convexity()["monotonic"] == False: - # store old degrees to compare later - num_degree_old = np.copy(self.rational_fit_kwargs["num_degree"]) - denom_degree_old = np.copy(self.rational_fit_kwargs["denom_degree"]) - # Increase degree by one - debug_message("Attempting to resolve nonmonotic egw evolution " - f"(current degree={num_degree_old}) by increasing rational fit " - "degree by 1.", self.debug_level, important=True) - self.rational_fit_kwargs["num_degree"] += 1 - self.rational_fit_kwargs["denom_degree"] += 1 - # check if this degree is already tried and corresponding monotonicity - # If it was already tried and resulted in nonmonotonic fits, then we abandon - # this attempt and set the rational fit degrees to their previous values. - if not self.check_whether_to_try_new_degree(self.rational_fit_kwargs["num_degree"]): - # reset the correct value - self.rational_fit_kwargs["num_degree"] = int(num_degree_old) - self.rational_fit_kwargs["denom_degree"] = int(denom_degree_old) - debug_message("Final rational fits were built with " - f"`num_degree`={self.rational_fit_kwargs['num_degree']} and " - f"`denom_degree`={self.rational_fit_kwargs['denom_degree']}.", - debug_level=self.debug_level, important=True) - break - # build interpolants with updated degree - self.build_omega_gw_extrema_interpolants() - # compute eccentricity using updated interpolants - self.eccentricity = self.compute_eccentricity(self.tref_out) - # update ecc for checks - self.ecc_for_checks = self.compute_eccentricity(self.t_for_checks) - - def check_whether_to_try_new_degree(self, new_degree): - if f"{new_degree}" in self.rational_fit_nonmonotonicity_history: - if self.rational_fit_nonmonotonicity_history[f"{new_degree}"]: - debug_message(f"Rational fit was already built with degree={new_degree} and " - f"was found to be nonmonotonic. Abandoning attempt with degree: {new_degree}", - debug_level=self.debug_level, important=True) - return not self.rational_fit_nonmonotonicity_history[f"{new_degree}"] - else: - return True - def build_omega_gw_extrema_interpolants(self): # Build the interpolants of omega_gw at the extrema if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": - self.omega_gw_pericenters_interp = self.interp_extrema("pericenters") self.omega_gw_apocenters_interp = self.interp_extrema("apocenters") + self.omega_gw_pericenters_interp = self.interp_extrema("pericenters") elif self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": - self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") + self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") else: raise Exception("Unknown method for `omega_gw_extrema_interpolation_method`. " "Must be one of `spline` or `rational_fit`.") @@ -1855,8 +1841,10 @@ def build_omega_gw_extrema_interpolants(self): if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": # Check if the first derivative of omega_gw at pericenters or apocenters is non-monotonic has_non_monotonicity = ( - self.check_domega_dt(self.t_for_checks, self.omega_gw_pericenters_interp(self.t_for_checks)) or - self.check_domega_dt(self.t_for_checks, self.omega_gw_apocenters_interp(self.t_for_checks)) + self.check_if_domega_dt_is_nonmonotonic( + self.t_for_checks, self.omega_gw_pericenters_interp(self.t_for_checks)) or + self.check_if_domega_dt_is_nonmonotonic( + self.t_for_checks, self.omega_gw_apocenters_interp(self.t_for_checks)) ) if has_non_monotonicity: From 0fca25fb3ee50fdca9e6aba5058ae668d5078071 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Fri, 25 Oct 2024 17:03:49 +0530 Subject: [PATCH 19/34] remove fallback option --- gw_eccentricity/eccDefinition.py | 56 +++++------------------------- gw_eccentricity/gw_eccentricity.py | 18 ---------- 2 files changed, 8 insertions(+), 66 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 82f52dae..61ca31cf 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -305,25 +305,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, - Can suppress pathologies in the waveform that might be visible with `spline`. - ### Fallback Behavior: - - With `omega_gw_extrema_interpolation_method` set to `spline`, if - `use_rational_fit_as_fallback` is set to `True`, the method will - initially use `spline`. If the first derivative of the spline interpolant - exhibits non-monotonicity, the code will automatically fall back to the - `rational_fit` method. This ensures a more reliable fit when the spline - method shows undesirable behavior. - Default value: `"spline"`. - - use_rational_fit_as_fallback : bool, default=True - Use rational fit for interpolant of omega at extrema when the - interpolant built using spline shows nonmonotonicity in its - first derivative. - - This is used only when `omega_gw_extrema_interpolation_method` is `spline`. - If `omega_gw_extrema_interpolation_method` is `rational_fit` then it has - no use. """ self.precessing = precessing # Get data necessary for eccentricity measurement @@ -377,7 +359,6 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, self.rational_fit_degrees = {"pericenters": None, "apocenters": None} self.debug_plots = self.extra_kwargs["debug_plots"] self.return_zero_if_small_ecc_failure = self.extra_kwargs["return_zero_if_small_ecc_failure"] - self.use_rational_fit_as_fallback = self.extra_kwargs["use_rational_fit_as_fallback"] # check if there are unrecognized keys in the dataDict self.recognized_dataDict_keys = self.get_recognized_dataDict_keys() for kw in dataDict.keys(): @@ -906,7 +887,6 @@ def get_default_extra_kwargs(self): "refine_extrema": False, "kwargs_for_fits_methods": {}, # Gets overriden in fits methods "return_zero_if_small_ecc_failure": False, - "use_rational_fit_as_fallback": True, "omega_gw_extrema_interpolation_method": "spline" } return default_extra_kwargs @@ -1833,38 +1813,18 @@ def build_omega_gw_extrema_interpolants(self): raise Exception("Unknown method for `omega_gw_extrema_interpolation_method`. " "Must be one of `spline` or `rational_fit`.") - # Verify the monotonicity of the first derivative of the omega_gw interpolant. - # If a spline is used for interpolation (as specified by 'omega_gw_extrema_interpolation_method'), - # non-monotonicity may occur in the first derivative. - # If 'use_rational_fit_as_fallback' is set to True, the spline interpolant - # will be replaced with a rational fit to ensure monotonic behavior. + # Verify the monotonicity of the first derivative of the omega_gw interpolant with spline. if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": # Check if the first derivative of omega_gw at pericenters or apocenters is non-monotonic - has_non_monotonicity = ( - self.check_if_domega_dt_is_nonmonotonic( + if (self.check_if_domega_dt_is_nonmonotonic( self.t_for_checks, self.omega_gw_pericenters_interp(self.t_for_checks)) or self.check_if_domega_dt_is_nonmonotonic( - self.t_for_checks, self.omega_gw_apocenters_interp(self.t_for_checks)) - ) - - if has_non_monotonicity: - if self.use_rational_fit_as_fallback: - debug_message( - "Non-monotonic time derivative detected in the spline interpolant through extrema. " - "Switching to rational fit.", - debug_level=self.debug_level, important=True - ) - # Use rational fit for both pericenters and apocenters - self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") - self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") - else: - debug_message( - "Non-monotonic time derivative detected in the spline interpolant through extrema. " - "Consider the following options to avoid this: \n" - " - Set 'use_rational_fit_as_fallback' to True in 'extra_kwargs' to switch to rational fit.\n" - " - Use rational fit directly by setting 'omega_gw_extrema_interpolation_method' to 'rational_fit'.", - debug_level=self.debug_level, important=True - ) + self.t_for_checks, self.omega_gw_apocenters_interp(self.t_for_checks))): + debug_message( + "Nonmonotonic time derivative detected in the spline interpolant through extrema. " + "Using rational fit by setting 'omega_gw_extrema_interpolation_method' to 'rational_fit' " + "may provide better result.", + debug_level=self.debug_level, important=True) def set_eccentricity_and_mean_anomaly_to_zero(self): """Set eccentricity and mean_anomaly to zero.""" diff --git a/gw_eccentricity/gw_eccentricity.py b/gw_eccentricity/gw_eccentricity.py index ffece41c..af442bbd 100644 --- a/gw_eccentricity/gw_eccentricity.py +++ b/gw_eccentricity/gw_eccentricity.py @@ -442,25 +442,7 @@ def measure_eccentricity(tref_in=None, - Can suppress pathologies in the waveform that might be visible with `spline`. - ### Fallback Behavior: - - With `omega_gw_extrema_interpolation_method` set to `spline`, if - `use_rational_fit_as_fallback` is set to `True`, the method will - initially use `spline`. If the first derivative of the spline interpolant - exhibits non-monotonicity, the code will automatically fall back to the - `rational_fit` method. This ensures a more reliable fit when the spline - method shows undesirable behavior. - Default value: `"spline"`. - - use_rational_fit_as_fallback : bool, default=True - Use rational fit for interpolant of omega at extrema when the - interpolant built using spline shows nonmonotonicity in its - first derivative. - - This is used only when `omega_gw_extrema_interpolation_method` is `spline`. - If `omega_gw_extrema_interpolation_method` is `rational_fit` then it has - no use. Returns ------- From 065776e6e349f75e1357917eb06713320d783b92 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Sat, 26 Oct 2024 11:26:59 +0530 Subject: [PATCH 20/34] improve documentation, change to rational_fit --- gw_eccentricity/eccDefinition.py | 73 ++++++++++++++++++------------ gw_eccentricity/gw_eccentricity.py | 38 +++++++--------- 2 files changed, 60 insertions(+), 51 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 61ca31cf..f05a64ac 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -279,33 +279,28 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, and mean anomaly to zero. USE THIS WITH CAUTION! - omega_gw_extrema_interpolation_method : str, default="spline" + 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`. - - ### When to Use: - - - **`spline`** (default): - - 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`**: - - More appropriate for noisy data, e.g., waveform modes from numerical - simulations. - - Minimizes least squares error, resulting in a smoother overall trend - with less oscillation. - - Significantly slower compared to the `spline` method. - - Can suppress pathologies in the waveform that might be visible with - `spline`. - - Default value: `"spline"`. + - Can handle both 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"`. """ self.precessing = precessing # Get data necessary for eccentricity measurement @@ -352,7 +347,9 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, self.available_averaging_methods \ = self.get_available_omega_gw_averaging_methods() self.debug_level = self.extra_kwargs["debug_level"] - self.rational_fit_kwargs["verbose"] = True if self.debug_level >= 1 else False + # set verbose to debug_level. If verbose is True, then it prints information + # of each iteration for rational fits for a given degree. + self.rational_fit_kwargs["verbose"] = self.debug_level # keep history of rational fit degree and nonmonotonicity of the corresponding fits self.rational_fit_nonmonotonicity_history = {"pericenters": {}, "apocenters": {}} # Update the degrees used to construct the final fits @@ -887,7 +884,7 @@ def get_default_extra_kwargs(self): "refine_extrema": False, "kwargs_for_fits_methods": {}, # Gets overriden in fits methods "return_zero_if_small_ecc_failure": False, - "omega_gw_extrema_interpolation_method": "spline" + "omega_gw_extrema_interpolation_method": "rational_fit" } return default_extra_kwargs @@ -1254,19 +1251,33 @@ def rational_fit(self, x, y, description=None): degrees to find a higher-degree fit. If increasing the degrees causes nonmonotonicity, it reverts to the previous valid monotonic fit. """ - # make sure that description is not None + # make sure that description is not None. A description is needed to + # update the optimal values of the numerator and denominator degrees + # used to build the final rational fit if description is None: - raise Exception("Please provide a description. `description` can not be None.") - # Set initial degrees if not already specified + raise Exception("Please provide a description for which to build a " + "rational fit. `description` can not be None. For example, " + "it can be 'apocenters', 'pericenters' or 'omega_gw_average'.") + # Set initial degrees if not already specified. + # The optimal degrees for the numerator and denominator changes depending on the + # eccentricity and duration of the waveform. We set these values by counting the + # number of approximate orbits and using some prior mapping between number of orbits + # and the optimal degrees. The mapping is done via `get_approximate_degree_for_rational_fit`. + # Usually this is a good starting point but not the optimal one + # as it varies as a function of eccentricity which we do not know apriori. The optimal + # value is found by doing some iterations. if not (self.rational_fit_kwargs["num_degree"] and self.rational_fit_kwargs["denom_degree"]): self.rational_fit_kwargs["num_degree"], self.rational_fit_kwargs["denom_degree"] \ = self.get_approximate_degree_for_rational_fit() rat_fit = self.get_rat_fit(x, y) t = np.arange(x[0], x[-1], self.t[1] - self.t[0]) + # save the degrees for checks at each step of iterations for finding the optimal degrees old_num_degree = self.rational_fit_kwargs["num_degree"] old_denom_degree = self.rational_fit_kwargs["denom_degree"] + # use a flag to know whether the initial degrees produced monotonic + # time derivative of the fits or it was required to be reduced. degrees_were_reduced = False # Check for nonmonotonicity and lower degrees if needed @@ -1277,17 +1288,19 @@ def rational_fit(self, x, y, description=None): if self.rational_fit_kwargs["denom_degree"] > 1: self.rational_fit_kwargs["denom_degree"] -= 1 if self.rational_fit_kwargs["num_degree"] == 1 and self.rational_fit_kwargs["denom_degree"] == 1: - raise Exception("Both numerator and denominator degrees cannot be lowered further.") + debug_message("Both numerator and denominator degrees are equal to 1 " + "and cannot be lowered further.", + debug_level=self.debug_level, important=False) + break debug_message(f"Lowering degrees to num_degree={self.rational_fit_kwargs['num_degree']}, " f"denom_degree={self.rational_fit_kwargs['denom_degree']} and retrying.", - debug_level=self.debug_level, important=True) + debug_level=self.debug_level, important=False) rat_fit = self.get_rat_fit(x, y) # If no degrees were reduced, try increasing the degree for better fit if not degrees_were_reduced: last_monotonic_rat_fit = rat_fit # Track last monotonic fit - # last_monotonic_rat_fit_vals = rat_fit(t) last_monotonic_num_degree = old_num_degree last_monotonic_denom_degree = old_denom_degree @@ -1303,7 +1316,7 @@ def rational_fit(self, x, y, description=None): debug_message(f"Increasing degrees caused nonmonotonicity. Reverting to " f"last monotonic fit with num_degree={last_monotonic_num_degree} " f"and denom_degree={last_monotonic_denom_degree}.", - debug_level=self.debug_level, important=True) + debug_level=self.debug_level, important=False) self.rational_fit_kwargs["num_degree"] = last_monotonic_num_degree self.rational_fit_kwargs["denom_degree"] = last_monotonic_denom_degree rat_fit = last_monotonic_rat_fit diff --git a/gw_eccentricity/gw_eccentricity.py b/gw_eccentricity/gw_eccentricity.py index af442bbd..58228601 100644 --- a/gw_eccentricity/gw_eccentricity.py +++ b/gw_eccentricity/gw_eccentricity.py @@ -416,33 +416,29 @@ def measure_eccentricity(tref_in=None, eccentricity is the cause, and set the returned eccentricity and mean anomaly to zero. USE THIS WITH CAUTION! - omega_gw_extrema_interpolation_method : str, default="spline" + + 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`. - - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation`. + - 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. - ### When to Use: - - - **`spline`** (default): - - 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`**: - - More appropriate for noisy data, e.g., waveform modes from numerical - simulations. - - Minimizes least squares error, resulting in a smoother overall trend - with less oscillation. - - Significantly slower compared to the `spline` method. - - Can suppress pathologies in the waveform that might be visible with - `spline`. - - Default value: `"spline"`. + - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation`. + - Can handle both 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 ------- From 671c54a9b4c71989ba2f055697acad008cd01045 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Sat, 26 Oct 2024 11:51:28 +0530 Subject: [PATCH 21/34] make tests pass --- gw_eccentricity/eccDefinition.py | 4 ++-- gw_eccentricity/gw_eccentricity.py | 4 ++-- test/test_mks_vs_dimless_units.py | 28 ++++++++++++++++++++-------- test/test_regression.py | 8 +++++++- 4 files changed, 31 insertions(+), 13 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index f05a64ac..9d6fd35c 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -292,8 +292,8 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, behavior, particularly near the merger. - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation`. - - Can handle both noisy data, e.g., waveform modes from numerical - simulations. + - 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 diff --git a/gw_eccentricity/gw_eccentricity.py b/gw_eccentricity/gw_eccentricity.py index 58228601..27d8969a 100644 --- a/gw_eccentricity/gw_eccentricity.py +++ b/gw_eccentricity/gw_eccentricity.py @@ -430,8 +430,8 @@ def measure_eccentricity(tref_in=None, behavior, particularly near the merger. - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation`. - - Can handle both noisy data, e.g., waveform modes from numerical - simulations. + - 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 diff --git a/test/test_mks_vs_dimless_units.py b/test/test_mks_vs_dimless_units.py index c48a4fd6..b915b3c4 100644 --- a/test/test_mks_vs_dimless_units.py +++ b/test/test_mks_vs_dimless_units.py @@ -38,6 +38,10 @@ def test_mks_vs_dimless_units(): "D": 1}) dataDictMKS = load_data.load_waveform(**lal_kwargs) + # set omega_gw_extrema_interpolation method + # TODO: Need to do the same for `rational_fit` + extra_kwargs = {"omega_gw_extrema_interpolation_method": "spline"} + # List of all available methods available_methods = gw_eccentricity.get_available_methods() for method in available_methods: @@ -46,7 +50,8 @@ def test_mks_vs_dimless_units(): gwecc_dict = measure_eccentricity( tref_in=dataDict["t"][idx], method=method, - dataDict=dataDict) + dataDict=dataDict, + extra_kwargs=extra_kwargs) tref_out = gwecc_dict["tref_out"] ecc_ref = gwecc_dict["eccentricity"] meanano_ref = gwecc_dict["mean_anomaly"] @@ -55,7 +60,8 @@ def test_mks_vs_dimless_units(): gwecc_dict_MKS = measure_eccentricity( tref_in=dataDictMKS["t"][idx], method=method, - dataDict=dataDictMKS) + dataDict=dataDictMKS, + extra_kwargs=extra_kwargs) tref_out_MKS = gwecc_dict_MKS["tref_out"] ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] @@ -88,7 +94,8 @@ def test_mks_vs_dimless_units(): gwecc_dict = measure_eccentricity( tref_in=dataDict["t"][idx_start: idx_end], method=method, - dataDict=dataDict) + dataDict=dataDict, + extra_kwargs=extra_kwargs) tref_out = gwecc_dict["tref_out"] ecc_ref = gwecc_dict["eccentricity"] meanano_ref = gwecc_dict["mean_anomaly"] @@ -96,7 +103,8 @@ def test_mks_vs_dimless_units(): gwecc_dict_MKS = measure_eccentricity( tref_in=dataDictMKS["t"][idx_start: idx_end], method=method, - dataDict=dataDictMKS) + dataDict=dataDictMKS, + extra_kwargs=extra_kwargs) tref_out_MKS = gwecc_dict_MKS["tref_out"] ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] @@ -129,7 +137,8 @@ def test_mks_vs_dimless_units(): gwecc_dict = measure_eccentricity( fref_in=fref_in, method=method, - dataDict=dataDict) + dataDict=dataDict, + extra_kwargs=extra_kwargs) fref_out = gwecc_dict["fref_out"] ecc_ref = gwecc_dict["eccentricity"] meanano_ref = gwecc_dict["mean_anomaly"] @@ -139,7 +148,8 @@ def test_mks_vs_dimless_units(): gwecc_dict_MKS = measure_eccentricity( fref_in=fref_in, method=method, - dataDict=dataDictMKS) + dataDict=dataDictMKS, + extra_kwargs=extra_kwargs) fref_out_MKS = gwecc_dict_MKS["fref_out"] ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] @@ -170,7 +180,8 @@ def test_mks_vs_dimless_units(): gwecc_dict = measure_eccentricity( fref_in=fref_in, method=method, - dataDict=dataDict) + dataDict=dataDict, + extra_kwargs=extra_kwargs) fref_out = gwecc_dict["fref_out"] ecc_ref = gwecc_dict["eccentricity"] meanano_ref = gwecc_dict["mean_anomaly"] @@ -180,7 +191,8 @@ def test_mks_vs_dimless_units(): gwecc_dict_MKS = measure_eccentricity( fref_in=fref_in, method=method, - dataDict=dataDictMKS) + dataDict=dataDictMKS, + extra_kwargs=extra_kwargs) fref_out_MKS = gwecc_dict_MKS["fref_out"] ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] diff --git a/test/test_regression.py b/test/test_regression.py index ca4e3345..74220991 100644 --- a/test/test_regression.py +++ b/test/test_regression.py @@ -6,7 +6,13 @@ from gw_eccentricity import measure_eccentricity # locally it passs without problem but on github action, we need to set this tolerance -atol = 1e-9 +atol = 1e-5 +#TODO: Changed atol from 1e-9 to 1e-5, so that tests can pass. This was required because +# the regression data was created using `spline`, whereas now we are using `rational_fit` +# as default. These two methods are expected to produdce slightly different values of eccenrtcity. +# We should created separate data for `spline` and `rational_fit` and then compare data with the +# corresponding method. + def test_regression(): """Regression test using all methods.""" From ae873cf8aea4256405408c7453e2f5a9a2b39740 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Sat, 26 Oct 2024 17:07:24 +0530 Subject: [PATCH 22/34] change egw monotonicity check to previous version --- gw_eccentricity/eccDefinition.py | 10 ---------- test/test_regression.py | 4 ++-- 2 files changed, 2 insertions(+), 12 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 9d6fd35c..533ff31a 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -2145,9 +2145,6 @@ def check_monotonicity_and_convexity(self, self.decc_dt_for_checks = self.derivative_of_eccentricity( self.t_for_checks, n=1) - # Create an empty dictionary to store the check results - check_dict = {} - # Is ecc(t) a monotonically decreasing function? if any(self.decc_dt_for_checks > 0): idx = np.where(self.decc_dt_for_checks > 0)[0] @@ -2156,9 +2153,6 @@ def check_monotonicity_and_convexity(self, f"{'at' if len(idx) == 1 else 'in the range'} {range}") debug_message(message, self.debug_level, point_to_verbose_output=True) - check_dict.update({"monotonic": False}) - else: - check_dict.update({"monotonic": True}) # Is ecc(t) a convex function? That is, is the second # derivative always negative? @@ -2174,10 +2168,6 @@ def check_monotonicity_and_convexity(self, debug_message(f"{message} expected to be always negative", self.debug_level, point_to_verbose_output=True) - check_dict.update({"convex": False}) - else: - check_dict.update({"convex": True}) - return check_dict def get_range_from_indices(self, indices, times): """Get the range of time from indices for gives times array.""" diff --git a/test/test_regression.py b/test/test_regression.py index 74220991..eb2d6c04 100644 --- a/test/test_regression.py +++ b/test/test_regression.py @@ -7,10 +7,10 @@ # locally it passs without problem but on github action, we need to set this tolerance atol = 1e-5 -#TODO: Changed atol from 1e-9 to 1e-5, so that tests can pass. This was required because +# TODO: Changed atol from 1e-9 to 1e-5, so that tests can pass. This was required because # the regression data was created using `spline`, whereas now we are using `rational_fit` # as default. These two methods are expected to produdce slightly different values of eccenrtcity. -# We should created separate data for `spline` and `rational_fit` and then compare data with the +# We should create separate data for `spline` and `rational_fit` and then compare data with the # corresponding method. From 336227dfc64a1a3f439f56c4005adc93397ddf3e Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Sat, 26 Oct 2024 19:13:59 +0530 Subject: [PATCH 23/34] add docs under get_rational_fit --- gw_eccentricity/eccDefinition.py | 1 + gw_eccentricity/utils.py | 7 ++++++- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 533ff31a..446c6bca 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -20,6 +20,7 @@ from .plot_settings import use_fancy_plotsettings, colorsDict, labelsDict from .plot_settings import figWidthsTwoColDict, figHeightsDict + class eccDefinition: """Base class to define eccentricity for given waveform data dictionary.""" diff --git a/gw_eccentricity/utils.py b/gw_eccentricity/utils.py index 29eed8cc..3411b19a 100644 --- a/gw_eccentricity/utils.py +++ b/gw_eccentricity/utils.py @@ -344,9 +344,14 @@ def get_rational_fit(x, y, rational_fit_kwargs=None, check_kwargs=True): "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) - return lambda x: rat(x.astype('float64').reshape(-1, 1)) + # 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, From c3b4eb87f9ccada99ffd801e1a17eebcf6a036ef Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Sun, 27 Oct 2024 09:38:53 +0530 Subject: [PATCH 24/34] update regression data and test based on interp method --- test/generate_regression_data.py | 13 ++- .../AmplitudeFits_regression_data.json | 1 - .../AmplitudeFits_spline_regression_data.json | 1 + .../Amplitude_regression_data.json | 1 - .../Amplitude_spline_regression_data.json | 1 + .../FrequencyFits_regression_data.json | 1 - .../FrequencyFits_spline_regression_data.json | 1 + .../Frequency_regression_data.json | 1 - .../Frequency_spline_regression_data.json | 1 + .../ResidualAmplitude_regression_data.json | 1 - ...idualAmplitude_spline_regression_data.json | 1 + .../ResidualFrequency_regression_data.json | 1 - ...idualFrequency_spline_regression_data.json | 1 + test/test_regression.py | 109 +++++++++--------- 14 files changed, 68 insertions(+), 66 deletions(-) delete mode 100644 test/regression_data/AmplitudeFits_regression_data.json create mode 100644 test/regression_data/AmplitudeFits_spline_regression_data.json delete mode 100644 test/regression_data/Amplitude_regression_data.json create mode 100644 test/regression_data/Amplitude_spline_regression_data.json delete mode 100644 test/regression_data/FrequencyFits_regression_data.json create mode 100644 test/regression_data/FrequencyFits_spline_regression_data.json delete mode 100644 test/regression_data/Frequency_regression_data.json create mode 100644 test/regression_data/Frequency_spline_regression_data.json delete mode 100644 test/regression_data/ResidualAmplitude_regression_data.json create mode 100644 test/regression_data/ResidualAmplitude_spline_regression_data.json delete mode 100644 test/regression_data/ResidualFrequency_regression_data.json create mode 100644 test/regression_data/ResidualFrequency_spline_regression_data.json diff --git a/test/generate_regression_data.py b/test/generate_regression_data.py index e17cbd0d..61bdd1bf 100644 --- a/test/generate_regression_data.py +++ b/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 @@ -79,9 +84,9 @@ def generate_regression_data(method): 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) diff --git a/test/regression_data/AmplitudeFits_regression_data.json b/test/regression_data/AmplitudeFits_regression_data.json deleted file mode 100644 index faeab019..00000000 --- a/test/regression_data/AmplitudeFits_regression_data.json +++ /dev/null @@ -1 +0,0 @@ -{"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": {}, "tref": {"time": [-14769.419163230703, -11192.719163230704, -7616.019163230703], "eccentricity": [0.14715861688347398, 0.1347862810647561, 0.11943002062852026], "mean_anomaly": [0.0, 0.9679656489638546, 5.4841245849525535]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338650648955707, 0.11542822840300593, 0.10530937888307923], "mean_anomaly": [2.9766378492619907, 1.494348360501732, 2.830869757710353]}} \ No newline at end of file diff --git a/test/regression_data/AmplitudeFits_spline_regression_data.json b/test/regression_data/AmplitudeFits_spline_regression_data.json new file mode 100644 index 00000000..3ab3cd16 --- /dev/null +++ b/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": [-14769.419163243456, -11192.719163243455, -7616.019163243456], "eccentricity": [0.1471586168829805, 0.1347862810667858, 0.11943002063131647], "mean_anomaly": [0.0, 0.9679656489638688, 5.4841245849525535]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338650648871319, 0.11542822840106648, 0.10530937888140646], "mean_anomaly": [2.9766378492612873, 1.494348360501732, 2.830869757710083]}} \ No newline at end of file diff --git a/test/regression_data/Amplitude_regression_data.json b/test/regression_data/Amplitude_regression_data.json deleted file mode 100644 index ec603f0b..00000000 --- a/test/regression_data/Amplitude_regression_data.json +++ /dev/null @@ -1 +0,0 @@ -{"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": {}, "tref": {"time": [-14768.619163230704, -11225.819163230703, -7683.019163230703], "eccentricity": [0.14714463101268582, 0.1348939080003978, 0.11972197672186036], "mean_anomaly": [0.0, 0.6424466289474253, 4.757738663295697]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12336087182171862, 0.11538856688498489, 0.10524430745756652], "mean_anomaly": [2.955433309853852, 1.4699121058753946, 2.7932679992191325]}} \ No newline at end of file diff --git a/test/regression_data/Amplitude_spline_regression_data.json b/test/regression_data/Amplitude_spline_regression_data.json new file mode 100644 index 00000000..f067d341 --- /dev/null +++ b/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": [-14768.619163243457, -11225.819163243457, -7683.019163243456], "eccentricity": [0.14714463101287545, 0.13489390800343737, 0.11972197672456475], "mean_anomaly": [0.0, 0.642446628947404, 4.757738663295697]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.123360871821657, 0.11538856688636778, 0.10524430745877145], "mean_anomaly": [2.955433309853831, 1.4699121058751956, 2.7932679992189193]}} \ No newline at end of file diff --git a/test/regression_data/FrequencyFits_regression_data.json b/test/regression_data/FrequencyFits_regression_data.json deleted file mode 100644 index 3a0e061e..00000000 --- a/test/regression_data/FrequencyFits_regression_data.json +++ /dev/null @@ -1 +0,0 @@ -{"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": {}, "tref": {"time": [-14769.519163230703, -11193.019163230703, -7616.419163230703], "eccentricity": [0.14715912645792362, 0.13478783801420957, 0.11943256094553112], "mean_anomaly": [0.0, 0.9671579875661607, 5.4821375349876575]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338317931464393, 0.11542730894258924, 0.1053075001590753], "mean_anomaly": [2.988070979817415, 1.5003684452927502, 2.8390552703477496]}} \ No newline at end of file diff --git a/test/regression_data/FrequencyFits_spline_regression_data.json b/test/regression_data/FrequencyFits_spline_regression_data.json new file mode 100644 index 00000000..11371e40 --- /dev/null +++ b/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": [-14769.519163243456, -11193.019163243456, -7616.419163243456], "eccentricity": [0.1471591264575487, 0.13478783801299066, 0.11943256094180654], "mean_anomaly": [0.0, 0.9671579875661678, 5.4821375349876575]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338317931488396, 0.1154273089414749, 0.10530750015791723], "mean_anomaly": [2.988070979816996, 1.5003684452927075, 2.839055270348034]}} \ No newline at end of file diff --git a/test/regression_data/Frequency_regression_data.json b/test/regression_data/Frequency_regression_data.json deleted file mode 100644 index 8eab7988..00000000 --- a/test/regression_data/Frequency_regression_data.json +++ /dev/null @@ -1 +0,0 @@ -{"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": {}, "tref": {"time": [-14768.719163230704, -11226.819163230703, -7684.919163230703], "eccentricity": [0.14714131326030366, 0.13489405451806347, 0.11972728556646606], "mean_anomaly": [0.0, 0.634809819093654, 4.7407656825570825]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12335278321610421, 0.11538428023363356, 0.10524044651268605], "mean_anomaly": [2.9664207187406504, 1.4731955044661476, 2.7980616912547234]}} \ No newline at end of file diff --git a/test/regression_data/Frequency_spline_regression_data.json b/test/regression_data/Frequency_spline_regression_data.json new file mode 100644 index 00000000..8bd97f73 --- /dev/null +++ b/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": [-14768.719163243457, -11226.819163243457, -7684.919163243456], "eccentricity": [0.14714131326021462, 0.13489405452014136, 0.11972728556745593], "mean_anomaly": [0.0, 0.634809819093654, 4.7407656825570825]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12335278321687171, 0.11538428023179614, 0.10524044651877096], "mean_anomaly": [2.9664207187394354, 1.473195504466105, 2.798061691253885]}} \ No newline at end of file diff --git a/test/regression_data/ResidualAmplitude_regression_data.json b/test/regression_data/ResidualAmplitude_regression_data.json deleted file mode 100644 index 75ef4a34..00000000 --- a/test/regression_data/ResidualAmplitude_regression_data.json +++ /dev/null @@ -1 +0,0 @@ -{"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": {}, "tref": {"time": [-14769.519163230703, -11193.019163230703, -7616.519163230703], "eccentricity": [0.14715912350724414, 0.134787764326868, 0.11943299066186941], "mean_anomaly": [0.0, 0.9662032412803185, 5.4800155541341695]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.123383193388651, 0.11542720219146674, 0.10530996061674747], "mean_anomaly": [2.9867814403803123, 1.499264733634135, 2.8323037712234154]}} \ No newline at end of file diff --git a/test/regression_data/ResidualAmplitude_spline_regression_data.json b/test/regression_data/ResidualAmplitude_spline_regression_data.json new file mode 100644 index 00000000..e0953d17 --- /dev/null +++ b/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": [-14769.519163243456, -11193.019163243456, -7616.519163243456], "eccentricity": [0.14715912350736482, 0.13478776432643125, 0.11943299066201574], "mean_anomaly": [0.0, 0.9662032412803043, 5.4800155541341695]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338319338698023, 0.1154272021945365, 0.10530996061768227], "mean_anomaly": [2.9867814403809874, 1.4992647336341633, 2.8323037712231454]}} \ No newline at end of file diff --git a/test/regression_data/ResidualFrequency_regression_data.json b/test/regression_data/ResidualFrequency_regression_data.json deleted file mode 100644 index b6df16de..00000000 --- a/test/regression_data/ResidualFrequency_regression_data.json +++ /dev/null @@ -1 +0,0 @@ -{"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": {}, "tref": {"time": [-14769.519163230703, -11193.019163230703, -7616.519163230703], "eccentricity": [0.14715912219591287, 0.13478782247732968, 0.11943303944129302], "mean_anomaly": [0.0, 0.9670110477321288, 5.481076544560906]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338313657709132, 0.11542727032631284, 0.10531052353292425], "mean_anomaly": [2.9880704763720587, 1.500384961620071, 2.8327681616288416]}} \ No newline at end of file diff --git a/test/regression_data/ResidualFrequency_spline_regression_data.json b/test/regression_data/ResidualFrequency_spline_regression_data.json new file mode 100644 index 00000000..30e3134c --- /dev/null +++ b/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": [-14769.519163243456, -11193.019163243456, -7616.519163243456], "eccentricity": [0.14715912219591198, 0.1347878224760185, 0.1194330394418408], "mean_anomaly": [0.0, 0.967011047732143, 5.481076544560906]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338313657401168, 0.11542727032728273, 0.10531052353279924], "mean_anomaly": [2.9880704763715116, 1.5003849616201421, 2.8327681616289127]}} \ No newline at end of file diff --git a/test/test_regression.py b/test/test_regression.py index eb2d6c04..59f92210 100644 --- a/test/test_regression.py +++ b/test/test_regression.py @@ -4,14 +4,10 @@ import gw_eccentricity from gw_eccentricity import load_data from gw_eccentricity import measure_eccentricity +import glob # locally it passs without problem but on github action, we need to set this tolerance -atol = 1e-5 -# TODO: Changed atol from 1e-9 to 1e-5, so that tests can pass. This was required because -# the regression data was created using `spline`, whereas now we are using `rational_fit` -# as default. These two methods are expected to produdce slightly different values of eccenrtcity. -# We should create separate data for `spline` and `rational_fit` and then compare data with the -# corresponding method. +atol = 1e-9 def test_regression(): @@ -20,55 +16,56 @@ def test_regression(): available_methods = gw_eccentricity.get_available_methods() for method in available_methods: # Load regression data - regression_data_file = f"test/regression_data/{method}_regression_data.json" - # Load the regression data - fl = open(regression_data_file, "r") - regression_data = json.load(fl) - fl.close() - # waveform kwargs - lal_kwargs = regression_data["waveform_kwargs"] - # load waveform data - dataDict = load_data.load_waveform(**lal_kwargs) - extra_kwargs = regression_data["extra_kwargs"] + regression_data_files = glob.glob(f"test/regression_data/{method}_*_regression_data.json") + for regression_data_file in regression_data_files: + # Load the regression data + fl = open(regression_data_file, "r") + regression_data = json.load(fl) + fl.close() + # waveform kwargs + lal_kwargs = regression_data["waveform_kwargs"] + # load waveform data + dataDict = load_data.load_waveform(**lal_kwargs) + extra_kwargs = regression_data["extra_kwargs"] - # Try evaluating at times where regression data are saved - regression_data_at_tref = regression_data["tref"] - tref_in = regression_data_at_tref["time"] - gwecc_dict = measure_eccentricity( - tref_in=tref_in, - method=method, - dataDict=dataDict, - extra_kwargs=extra_kwargs) - tref_out = gwecc_dict["tref_out"] - ecc_ref = gwecc_dict["eccentricity"] - meanano_ref = gwecc_dict["mean_anomaly"] - # Compare the measured data with the saved data - np.testing.assert_allclose( - ecc_ref, regression_data_at_tref["eccentricity"], - atol=atol, - err_msg="measured and saved eccentricity at saved times do not match.") - np.testing.assert_allclose( - meanano_ref, regression_data_at_tref["mean_anomaly"], - atol=atol, - err_msg="measured and saved mean anomaly at saved times do not match.") + # Try evaluating at times where regression data are saved + regression_data_at_tref = regression_data["tref"] + tref_in = regression_data_at_tref["time"] + gwecc_dict = measure_eccentricity( + tref_in=tref_in, + method=method, + dataDict=dataDict, + extra_kwargs=extra_kwargs) + tref_out = gwecc_dict["tref_out"] + ecc_ref = gwecc_dict["eccentricity"] + meanano_ref = gwecc_dict["mean_anomaly"] + # Compare the measured data with the saved data + np.testing.assert_allclose( + ecc_ref, regression_data_at_tref["eccentricity"], + atol=atol, + err_msg="measured and saved eccentricity at saved times do not match.") + np.testing.assert_allclose( + meanano_ref, regression_data_at_tref["mean_anomaly"], + atol=atol, + err_msg="measured and saved mean anomaly at saved times do not match.") - # Try evaluating at frequencies where regression data are saved - regression_data_at_fref = regression_data["fref"] - fref_in = regression_data_at_fref["frequency"] - gwecc_dict = measure_eccentricity( - fref_in=fref_in, - method=method, - dataDict=dataDict, - extra_kwargs=extra_kwargs) - fref_out = gwecc_dict["fref_out"] - ecc_ref = gwecc_dict["eccentricity"] - meanano_ref = gwecc_dict["mean_anomaly"] - # Compare the measured data with the saved data - np.testing.assert_allclose( - ecc_ref, regression_data_at_fref["eccentricity"], - atol=atol, - err_msg="measured and saved eccentricity at saved frequencies do not match.") - np.testing.assert_allclose( - meanano_ref, regression_data_at_fref["mean_anomaly"], - atol=atol, - err_msg="measured and saved mean anomaly at saved frequencies do not match.") + # Try evaluating at frequencies where regression data are saved + regression_data_at_fref = regression_data["fref"] + fref_in = regression_data_at_fref["frequency"] + gwecc_dict = measure_eccentricity( + fref_in=fref_in, + method=method, + dataDict=dataDict, + extra_kwargs=extra_kwargs) + fref_out = gwecc_dict["fref_out"] + ecc_ref = gwecc_dict["eccentricity"] + meanano_ref = gwecc_dict["mean_anomaly"] + # Compare the measured data with the saved data + np.testing.assert_allclose( + ecc_ref, regression_data_at_fref["eccentricity"], + atol=atol, + err_msg="measured and saved eccentricity at saved frequencies do not match.") + np.testing.assert_allclose( + meanano_ref, regression_data_at_fref["mean_anomaly"], + atol=atol, + err_msg="measured and saved mean anomaly at saved frequencies do not match.") From 9c376dbd3ef745d9f64aa6e3f93d908dfe7aa926 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Sun, 27 Oct 2024 09:57:30 +0530 Subject: [PATCH 25/34] fix regression test --- test/generate_regression_data.py | 12 ++++++------ .../AmplitudeFits_spline_regression_data.json | 2 +- .../Amplitude_spline_regression_data.json | 2 +- .../FrequencyFits_spline_regression_data.json | 2 +- .../Frequency_spline_regression_data.json | 2 +- .../ResidualAmplitude_spline_regression_data.json | 2 +- .../ResidualFrequency_spline_regression_data.json | 2 +- 7 files changed, 12 insertions(+), 12 deletions(-) diff --git a/test/generate_regression_data.py b/test/generate_regression_data.py index 61bdd1bf..6b309835 100644 --- a/test/generate_regression_data.py +++ b/test/generate_regression_data.py @@ -61,9 +61,9 @@ def generate_regression_data(method, interp_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( @@ -75,9 +75,9 @@ def generate_regression_data(method, interp_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}) diff --git a/test/regression_data/AmplitudeFits_spline_regression_data.json b/test/regression_data/AmplitudeFits_spline_regression_data.json index 3ab3cd16..002b5202 100644 --- a/test/regression_data/AmplitudeFits_spline_regression_data.json +++ b/test/regression_data/AmplitudeFits_spline_regression_data.json @@ -1 +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": [-14769.419163243456, -11192.719163243455, -7616.019163243456], "eccentricity": [0.1471586168829805, 0.1347862810667858, 0.11943002063131647], "mean_anomaly": [0.0, 0.9679656489638688, 5.4841245849525535]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338650648871319, 0.11542822840106648, 0.10530937888140646], "mean_anomaly": [2.9766378492612873, 1.494348360501732, 2.830869757710083]}} \ No newline at end of file +{"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]}} \ No newline at end of file diff --git a/test/regression_data/Amplitude_spline_regression_data.json b/test/regression_data/Amplitude_spline_regression_data.json index f067d341..38be1630 100644 --- a/test/regression_data/Amplitude_spline_regression_data.json +++ b/test/regression_data/Amplitude_spline_regression_data.json @@ -1 +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": [-14768.619163243457, -11225.819163243457, -7683.019163243456], "eccentricity": [0.14714463101287545, 0.13489390800343737, 0.11972197672456475], "mean_anomaly": [0.0, 0.642446628947404, 4.757738663295697]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.123360871821657, 0.11538856688636778, 0.10524430745877145], "mean_anomaly": [2.955433309853831, 1.4699121058751956, 2.7932679992189193]}} \ No newline at end of file +{"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]}} \ No newline at end of file diff --git a/test/regression_data/FrequencyFits_spline_regression_data.json b/test/regression_data/FrequencyFits_spline_regression_data.json index 11371e40..d532f681 100644 --- a/test/regression_data/FrequencyFits_spline_regression_data.json +++ b/test/regression_data/FrequencyFits_spline_regression_data.json @@ -1 +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": [-14769.519163243456, -11193.019163243456, -7616.419163243456], "eccentricity": [0.1471591264575487, 0.13478783801299066, 0.11943256094180654], "mean_anomaly": [0.0, 0.9671579875661678, 5.4821375349876575]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338317931488396, 0.1154273089414749, 0.10530750015791723], "mean_anomaly": [2.988070979816996, 1.5003684452927075, 2.839055270348034]}} \ No newline at end of file +{"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]}} \ No newline at end of file diff --git a/test/regression_data/Frequency_spline_regression_data.json b/test/regression_data/Frequency_spline_regression_data.json index 8bd97f73..9025b76f 100644 --- a/test/regression_data/Frequency_spline_regression_data.json +++ b/test/regression_data/Frequency_spline_regression_data.json @@ -1 +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": [-14768.719163243457, -11226.819163243457, -7684.919163243456], "eccentricity": [0.14714131326021462, 0.13489405452014136, 0.11972728556745593], "mean_anomaly": [0.0, 0.634809819093654, 4.7407656825570825]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12335278321687171, 0.11538428023179614, 0.10524044651877096], "mean_anomaly": [2.9664207187394354, 1.473195504466105, 2.798061691253885]}} \ No newline at end of file +{"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]}} \ No newline at end of file diff --git a/test/regression_data/ResidualAmplitude_spline_regression_data.json b/test/regression_data/ResidualAmplitude_spline_regression_data.json index e0953d17..9caccb28 100644 --- a/test/regression_data/ResidualAmplitude_spline_regression_data.json +++ b/test/regression_data/ResidualAmplitude_spline_regression_data.json @@ -1 +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": [-14769.519163243456, -11193.019163243456, -7616.519163243456], "eccentricity": [0.14715912350736482, 0.13478776432643125, 0.11943299066201574], "mean_anomaly": [0.0, 0.9662032412803043, 5.4800155541341695]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338319338698023, 0.1154272021945365, 0.10530996061768227], "mean_anomaly": [2.9867814403809874, 1.4992647336341633, 2.8323037712231454]}} \ No newline at end of file +{"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]}} \ No newline at end of file diff --git a/test/regression_data/ResidualFrequency_spline_regression_data.json b/test/regression_data/ResidualFrequency_spline_regression_data.json index 30e3134c..1dc9abd0 100644 --- a/test/regression_data/ResidualFrequency_spline_regression_data.json +++ b/test/regression_data/ResidualFrequency_spline_regression_data.json @@ -1 +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": [-14769.519163243456, -11193.019163243456, -7616.519163243456], "eccentricity": [0.14715912219591198, 0.1347878224760185, 0.1194330394418408], "mean_anomaly": [0.0, 0.967011047732143, 5.481076544560906]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338313657401168, 0.11542727032728273, 0.10531052353279924], "mean_anomaly": [2.9880704763715116, 1.5003849616201421, 2.8327681616289127]}} \ No newline at end of file +{"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]}} \ No newline at end of file From 03c75ae1a0caca397a9d0c6a7132b5e613981525 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Sun, 27 Oct 2024 11:04:11 +0530 Subject: [PATCH 26/34] ignore deprecationwarning from polyrat --- pytest.ini | 3 +++ test/test_regression.py | 1 - 2 files changed, 3 insertions(+), 1 deletion(-) create mode 100644 pytest.ini diff --git a/pytest.ini b/pytest.ini new file mode 100644 index 00000000..8a199c10 --- /dev/null +++ b/pytest.ini @@ -0,0 +1,3 @@ +[pytest] +filterwarnings = + ignore::DeprecationWarning:polyrat.* \ No newline at end of file diff --git a/test/test_regression.py b/test/test_regression.py index 59f92210..2de08489 100644 --- a/test/test_regression.py +++ b/test/test_regression.py @@ -9,7 +9,6 @@ # locally it passs without problem but on github action, we need to set this tolerance atol = 1e-9 - def test_regression(): """Regression test using all methods.""" # List of all available methods From c2c361f390cbb1a061b0d2e3483139812db113bd Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Sun, 27 Oct 2024 11:47:37 +0530 Subject: [PATCH 27/34] test units for spline for now --- gw_eccentricity/eccDefinition.py | 6 +- test/test_mks_vs_dimless_units.py | 345 +++++++++++++++--------------- 2 files changed, 178 insertions(+), 173 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 446c6bca..d3be637e 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -1824,8 +1824,10 @@ def build_omega_gw_extrema_interpolants(self): self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") else: - raise Exception("Unknown method for `omega_gw_extrema_interpolation_method`. " - "Must be one of `spline` or `rational_fit`.") + raise Exception( + f"Unknown method {self.extra_kwargs['omega_gw_extrema_interpolation_method']} " + "for `omega_gw_extrema_interpolation_method`. " + "Must be one of `spline` or `rational_fit`.") # Verify the monotonicity of the first derivative of the omega_gw interpolant with spline. if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": diff --git a/test/test_mks_vs_dimless_units.py b/test/test_mks_vs_dimless_units.py index b915b3c4..769abf2e 100644 --- a/test/test_mks_vs_dimless_units.py +++ b/test/test_mks_vs_dimless_units.py @@ -38,181 +38,184 @@ def test_mks_vs_dimless_units(): "D": 1}) dataDictMKS = load_data.load_waveform(**lal_kwargs) - # set omega_gw_extrema_interpolation method - # TODO: Need to do the same for `rational_fit` - extra_kwargs = {"omega_gw_extrema_interpolation_method": "spline"} + # use different omega_gw_extrema_interpolation methods + # TODO: Add rational_fit to the list + omega_gw_extrema_interpolation_methods = ["spline"] # List of all available methods available_methods = gw_eccentricity.get_available_methods() for method in available_methods: - # Try evaluating at a single dimless time t = -12000 - idx = np.argmin(np.abs(dataDict["t"] - (-12000))) - gwecc_dict = measure_eccentricity( - tref_in=dataDict["t"][idx], - method=method, - dataDict=dataDict, - extra_kwargs=extra_kwargs) - tref_out = gwecc_dict["tref_out"] - ecc_ref = gwecc_dict["eccentricity"] - meanano_ref = gwecc_dict["mean_anomaly"] - gwecc_object = gwecc_dict["gwecc_object"] - # Try evaluating at a single MKS time - gwecc_dict_MKS = measure_eccentricity( - tref_in=dataDictMKS["t"][idx], - method=method, - dataDict=dataDictMKS, - extra_kwargs=extra_kwargs) - tref_out_MKS = gwecc_dict_MKS["tref_out"] - ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] - meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] - gwecc_object_MKS = gwecc_dict_MKS["gwecc_object"] - # Check that the tref_out is the same as rescaled (to dimless) tref_out_MKS - sec_to_dimless = 1/lal_kwargs["M"]/MTSUN_SI - np.testing.assert_allclose( - [tref_out], - [tref_out_MKS * sec_to_dimless], - 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], - 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], - err_msg=("Mean anomaly at a single dimensionless and MKS" - " time gives different results.\n" - "x = Dimensionless, y = MKS")) + for omega_gw_extrema_interpolation_method in omega_gw_extrema_interpolation_methods: + extra_kwargs = { + "omega_gw_extrema_interpolation_method": omega_gw_extrema_interpolation_method} + # Try evaluating at a single dimless time t = -12000 + idx = np.argmin(np.abs(dataDict["t"] - (-12000))) + gwecc_dict = measure_eccentricity( + tref_in=dataDict["t"][idx], + method=method, + dataDict=dataDict, + extra_kwargs=extra_kwargs) + tref_out = gwecc_dict["tref_out"] + ecc_ref = gwecc_dict["eccentricity"] + meanano_ref = gwecc_dict["mean_anomaly"] + gwecc_object = gwecc_dict["gwecc_object"] + # Try evaluating at a single MKS time + gwecc_dict_MKS = measure_eccentricity( + tref_in=dataDictMKS["t"][idx], + method=method, + dataDict=dataDictMKS, + extra_kwargs=extra_kwargs) + tref_out_MKS = gwecc_dict_MKS["tref_out"] + ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] + meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] + gwecc_object_MKS = gwecc_dict_MKS["gwecc_object"] + # Check that the tref_out is the same as rescaled (to dimless) tref_out_MKS + sec_to_dimless = 1/lal_kwargs["M"]/MTSUN_SI + np.testing.assert_allclose( + [tref_out], + [tref_out_MKS * sec_to_dimless], + 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], + 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], + err_msg=("Mean anomaly at a single dimensionless and MKS" + " time gives different results.\n" + "x = Dimensionless, y = MKS")) - # Try evaluating at an array of dimless times - idx_start = np.argmin(np.abs(dataDict["t"] - (-12000))) - idx_end = np.argmin(np.abs(dataDict["t"] - (-5000))) - gwecc_dict = measure_eccentricity( - tref_in=dataDict["t"][idx_start: idx_end], - method=method, - dataDict=dataDict, - extra_kwargs=extra_kwargs) - tref_out = gwecc_dict["tref_out"] - ecc_ref = gwecc_dict["eccentricity"] - meanano_ref = gwecc_dict["mean_anomaly"] - # Try evaluating at an array of MKS times - gwecc_dict_MKS = measure_eccentricity( - tref_in=dataDictMKS["t"][idx_start: idx_end], - method=method, - dataDict=dataDictMKS, - extra_kwargs=extra_kwargs) - tref_out_MKS = gwecc_dict_MKS["tref_out"] - ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] - meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] - # check that the tref_out times are consistent - np.testing.assert_allclose( - [tref_out], - [tref_out_MKS * sec_to_dimless], - 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, - err_msg=("Eccentricity at dimensionless and MKS array of" - " times are different\n." - "x = Dimensionless, y = MKS")) - # using unwrapped mean anomaly since 0 and 2pi should be treated as - # the same - np.testing.assert_allclose( - np.unwrap(meanano_ref), - np.unwrap(meanano_ref_MKS), - err_msg=("Mean anomaly at dimensionless and MKS array of" - " times are different.\n" - "x = Dimensionless, y = MKS")) + # Try evaluating at an array of dimless times + idx_start = np.argmin(np.abs(dataDict["t"] - (-12000))) + idx_end = np.argmin(np.abs(dataDict["t"] - (-5000))) + gwecc_dict = measure_eccentricity( + tref_in=dataDict["t"][idx_start: idx_end], + method=method, + dataDict=dataDict, + extra_kwargs=extra_kwargs) + tref_out = gwecc_dict["tref_out"] + ecc_ref = gwecc_dict["eccentricity"] + meanano_ref = gwecc_dict["mean_anomaly"] + # Try evaluating at an array of MKS times + gwecc_dict_MKS = measure_eccentricity( + tref_in=dataDictMKS["t"][idx_start: idx_end], + method=method, + dataDict=dataDictMKS, + extra_kwargs=extra_kwargs) + tref_out_MKS = gwecc_dict_MKS["tref_out"] + ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] + meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] + # check that the tref_out times are consistent + np.testing.assert_allclose( + [tref_out], + [tref_out_MKS * sec_to_dimless], + 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, + err_msg=("Eccentricity at dimensionless and MKS array of" + " times are different\n." + "x = Dimensionless, y = MKS")) + # using unwrapped mean anomaly since 0 and 2pi should be treated as + # the same + np.testing.assert_allclose( + np.unwrap(meanano_ref), + np.unwrap(meanano_ref_MKS), + err_msg=("Mean anomaly at dimensionless and MKS array of" + " times are different.\n" + "x = Dimensionless, y = MKS")) - # Try evaluating at single dimensionless frequency - fref_in = gwecc_object.compute_orbit_averaged_omega_gw_between_extrema( - dataDict["t"][idx]) / (2 * np.pi) - gwecc_dict = measure_eccentricity( - fref_in=fref_in, - method=method, - dataDict=dataDict, - extra_kwargs=extra_kwargs) - fref_out = gwecc_dict["fref_out"] - ecc_ref = gwecc_dict["eccentricity"] - meanano_ref = gwecc_dict["mean_anomaly"] - # Try evaluating at single MKS frequency - fref_in = gwecc_object_MKS.compute_orbit_averaged_omega_gw_between_extrema( - dataDictMKS["t"][idx]) / (2 * np.pi) - gwecc_dict_MKS = measure_eccentricity( - fref_in=fref_in, - method=method, - dataDict=dataDictMKS, - extra_kwargs=extra_kwargs) - fref_out_MKS = gwecc_dict_MKS["fref_out"] - ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] - meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] - # Check the fref_out frequencies are consistent - np.testing.assert_allclose( - [fref_out], - [fref_out_MKS / sec_to_dimless], - 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], - 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], - err_msg=("Mean anomaly at a single dimensionless and MKS" - " frequency gives different results.\n" - "x = Dimensionless, y = MKS")) + # Try evaluating at single dimensionless frequency + fref_in = gwecc_object.compute_orbit_averaged_omega_gw_between_extrema( + dataDict["t"][idx]) / (2 * np.pi) + gwecc_dict = measure_eccentricity( + fref_in=fref_in, + method=method, + dataDict=dataDict, + extra_kwargs=extra_kwargs) + fref_out = gwecc_dict["fref_out"] + ecc_ref = gwecc_dict["eccentricity"] + meanano_ref = gwecc_dict["mean_anomaly"] + # Try evaluating at single MKS frequency + fref_in = gwecc_object_MKS.compute_orbit_averaged_omega_gw_between_extrema( + dataDictMKS["t"][idx]) / (2 * np.pi) + gwecc_dict_MKS = measure_eccentricity( + fref_in=fref_in, + method=method, + dataDict=dataDictMKS, + extra_kwargs=extra_kwargs) + fref_out_MKS = gwecc_dict_MKS["fref_out"] + ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] + meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] + # Check the fref_out frequencies are consistent + np.testing.assert_allclose( + [fref_out], + [fref_out_MKS / sec_to_dimless], + 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], + 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], + err_msg=("Mean anomaly at a single dimensionless and MKS" + " frequency gives different results.\n" + "x = Dimensionless, y = MKS")) - # Try evaluating at an array of dimensionless frequencies - fref_in = gwecc_object.compute_orbit_averaged_omega_gw_between_extrema( - dataDict["t"][idx: idx+500]) / (2 * np.pi) - gwecc_dict = measure_eccentricity( - fref_in=fref_in, - method=method, - dataDict=dataDict, - extra_kwargs=extra_kwargs) - fref_out = gwecc_dict["fref_out"] - ecc_ref = gwecc_dict["eccentricity"] - meanano_ref = gwecc_dict["mean_anomaly"] - # Try evaluating at an array of MKS frequencies - fref_in = gwecc_object_MKS.compute_orbit_averaged_omega_gw_between_extrema( - dataDictMKS["t"][idx: idx+500]) / (2 * np.pi) - gwecc_dict_MKS = measure_eccentricity( - fref_in=fref_in, - method=method, - dataDict=dataDictMKS, - extra_kwargs=extra_kwargs) - fref_out_MKS = gwecc_dict_MKS["fref_out"] - ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] - meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] - # Check fref_out - np.testing.assert_allclose( - [fref_out], - [fref_out_MKS / sec_to_dimless], - 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, - 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), - err_msg=("Mean anomaly at dimensionless and MKS array of" - " frequencies are different.\n" - "x = Dimensionless, y = MKS")) + # Try evaluating at an array of dimensionless frequencies + fref_in = gwecc_object.compute_orbit_averaged_omega_gw_between_extrema( + dataDict["t"][idx: idx+500]) / (2 * np.pi) + gwecc_dict = measure_eccentricity( + fref_in=fref_in, + method=method, + dataDict=dataDict, + extra_kwargs=extra_kwargs) + fref_out = gwecc_dict["fref_out"] + ecc_ref = gwecc_dict["eccentricity"] + meanano_ref = gwecc_dict["mean_anomaly"] + # Try evaluating at an array of MKS frequencies + fref_in = gwecc_object_MKS.compute_orbit_averaged_omega_gw_between_extrema( + dataDictMKS["t"][idx: idx+500]) / (2 * np.pi) + gwecc_dict_MKS = measure_eccentricity( + fref_in=fref_in, + method=method, + dataDict=dataDictMKS, + extra_kwargs=extra_kwargs) + fref_out_MKS = gwecc_dict_MKS["fref_out"] + ecc_ref_MKS = gwecc_dict_MKS["eccentricity"] + meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"] + # Check fref_out + np.testing.assert_allclose( + [fref_out], + [fref_out_MKS / sec_to_dimless], + 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, + 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), + err_msg=("Mean anomaly at dimensionless and MKS array of" + " frequencies are different.\n" + "x = Dimensionless, y = MKS")) From 07336b0d38836232ff1274313e95df1dfc6c4064 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Tue, 5 Nov 2024 19:51:17 +0530 Subject: [PATCH 28/34] address comments --- gw_eccentricity/eccDefinition.py | 159 ++++++++++++++++++------------- 1 file changed, 95 insertions(+), 64 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index d3be637e..dd9db313 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -290,7 +290,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, 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. + behavior, particularly near the merger, especially for noisy NR data. - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation`. - Can handle both clean and noisy data, e.g., waveform @@ -1229,7 +1229,7 @@ def interp_extrema(self, extrema_type="pericenters"): f"Sufficient number of {extrema_type} are not found." " Can not create an interpolant.") - def get_rat_fit(self, x, y): + def get_rational_fit_wrapper(self, x, y): """Get rational fit. A wrapper of `utils.get_rational_fit` with check_kwargs=False. This is @@ -1242,48 +1242,63 @@ def get_rat_fit(self, x, y): rational_fit_kwargs=self.rational_fit_kwargs, check_kwargs=False) - def rational_fit(self, x, y, description=None): + def rational_fit(self, x, y, data_name=None): """Get rational fit with adaptive numerator and denominator degree. This function ensures that the rational fit uses the optimal degree for - the numerator and denominator. It first checks for nonmonotonicity in the - fit's derivative and lowers the degrees if necessary. Only if the initial - degree was monotonic (no reduction needed), it attempts to increase the - degrees to find a higher-degree fit. If increasing the degrees causes - nonmonotonicity, it reverts to the previous valid monotonic fit. + the numerator and denominator. It first checks for nonmonotonicity in + the fit's derivative and lowers the degrees if necessary. Only if the + fit for the initial degree was monotonic (no reduction of degree + needed), it attempts to increase the degrees to find a higher-degree + fit. If increasing the degrees causes nonmonotonicity, it reverts to + the previous valid monotonic fit. + + The initial degrees for the rational fit can be specified through + `rational_fit_kwargs` in `extra_kwargs`. Default values are provided by + `get_default_rational_fit_kwargs`, where both degrees are set to + `None`. If both degrees remain `None`, appropriate starting values are + determined using `self.get_approximate_degree_for_rational_fit`. """ - # make sure that description is not None. A description is needed to + # make sure that data_name is not None. A data_name is needed to # update the optimal values of the numerator and denominator degrees # used to build the final rational fit - if description is None: - raise Exception("Please provide a description for which to build a " - "rational fit. `description` can not be None. For example, " + if data_name is None: + raise Exception("Please provide a data_name for which to build a " + "rational fit. `data_name` can not be None. For example, " "it can be 'apocenters', 'pericenters' or 'omega_gw_average'.") - # Set initial degrees if not already specified. - # The optimal degrees for the numerator and denominator changes depending on the - # eccentricity and duration of the waveform. We set these values by counting the - # number of approximate orbits and using some prior mapping between number of orbits - # and the optimal degrees. The mapping is done via `get_approximate_degree_for_rational_fit`. - # Usually this is a good starting point but not the optimal one - # as it varies as a function of eccentricity which we do not know apriori. The optimal - # value is found by doing some iterations. - if not (self.rational_fit_kwargs["num_degree"] and self.rational_fit_kwargs["denom_degree"]): - self.rational_fit_kwargs["num_degree"], self.rational_fit_kwargs["denom_degree"] \ - = self.get_approximate_degree_for_rational_fit() - - rat_fit = self.get_rat_fit(x, y) - t = np.arange(x[0], x[-1], self.t[1] - self.t[0]) + # Set initial degrees for the rational fit if not already specified. + # If either `num_degree` or `denom_degree` is specified, use that value + # to set the other. If both are unspecified, we use the following + # procedure to set appropriate values for `num_degree` and `denom_degree`. + + # The optimal degrees for the numerator and denominator can vary with + # the eccentricity and duration of the waveform. To estimate these values, + # we first count the approximate number of orbits, then use a predefined + # mapping from orbit count to optimal degrees, provided by + # `get_approximate_degree_for_rational_fit`. + + # While this provides a good initial estimate, it may not be optimal, + # as the ideal degree can change based on the eccentricity which we do + # not know apriori. The true optimal values are determined through + # further iterations. + if self.rational_fit_kwargs["num_degree"] is None and self.rational_fit_kwargs["denom_degree"] is None: + self.rational_fit_kwargs["num_degree"], self.rational_fit_kwargs["denom_degree"] = \ + self.get_approximate_degree_for_rational_fit() + elif self.rational_fit_kwargs["num_degree"] is None: + self.rational_fit_kwargs["num_degree"] = self.rational_fit_kwargs["denom_degree"] + elif self.rational_fit_kwargs["denom_degree"] is None: + self.rational_fit_kwargs["denom_degree"] = self.rational_fit_kwargs["num_degree"] + + rat_fit = self.get_rational_fit_wrapper(x, y) + x_test = np.arange(x[0], x[-1], self.t[1] - self.t[0]) # save the degrees for checks at each step of iterations for finding the optimal degrees old_num_degree = self.rational_fit_kwargs["num_degree"] old_denom_degree = self.rational_fit_kwargs["denom_degree"] - # use a flag to know whether the initial degrees produced monotonic - # time derivative of the fits or it was required to be reduced. - degrees_were_reduced = False - # Check for nonmonotonicity and lower degrees if needed - while self.check_if_domega_dt_is_nonmonotonic(t, rat_fit(t), 1.0, description): - degrees_were_reduced = True # Mark that reduction occurred + fit_is_nonmonotonic = self.check_if_first_derivative_is_not_strictly_monotonic( + x_test, rat_fit(x_test), 1.0, data_name) + while fit_is_nonmonotonic: if self.rational_fit_kwargs["num_degree"] > 1: self.rational_fit_kwargs["num_degree"] -= 1 if self.rational_fit_kwargs["denom_degree"] > 1: @@ -1298,21 +1313,27 @@ def rational_fit(self, x, y, description=None): f"denom_degree={self.rational_fit_kwargs['denom_degree']} and retrying.", debug_level=self.debug_level, important=False) - rat_fit = self.get_rat_fit(x, y) - # If no degrees were reduced, try increasing the degree for better fit - if not degrees_were_reduced: + # build new fit and check monotonicity + rat_fit = self.get_rational_fit_wrapper(x, y) + fit_is_nonmonotonic = self.check_if_first_derivative_is_not_strictly_monotonic( + x_test, rat_fit(x_test), 1.0, data_name) + # If fit with initial degree is monotonic, try increasing the degree for a better fit + if not fit_is_nonmonotonic: last_monotonic_rat_fit = rat_fit # Track last monotonic fit last_monotonic_num_degree = old_num_degree last_monotonic_denom_degree = old_denom_degree - while not self.check_if_domega_dt_is_nonmonotonic(t, rat_fit(t), 1.0, description): + while not fit_is_nonmonotonic: # Increase the degrees for both numerator and denominator new_num_degree = self.rational_fit_kwargs["num_degree"] + 1 new_denom_degree = self.rational_fit_kwargs["denom_degree"] + 1 self.rational_fit_kwargs["num_degree"], self.rational_fit_kwargs["denom_degree"] \ = new_num_degree, new_denom_degree - new_rat_fit = self.get_rat_fit(x, y) - if self.check_if_domega_dt_is_nonmonotonic(t, new_rat_fit(t), 1.0, description): + # build new fit and check monotonicity + new_rat_fit = self.get_rational_fit_wrapper(x, y) + fit_is_nonmonotonic = self.check_if_first_derivative_is_not_strictly_monotonic( + x_test, new_rat_fit(x_test), 1.0, data_name) + if fit_is_nonmonotonic: # Revert to previous fit and degrees if nonmonotonicity is detected debug_message(f"Increasing degrees caused nonmonotonicity. Reverting to " f"last monotonic fit with num_degree={last_monotonic_num_degree} " @@ -1329,7 +1350,7 @@ def rational_fit(self, x, y, description=None): rat_fit = new_rat_fit # update final degrees used to build the fit - self.rational_fit_degrees[description] = (self.rational_fit_kwargs["num_degree"], + self.rational_fit_degrees[data_name] = (self.rational_fit_kwargs["num_degree"], self.rational_fit_kwargs["denom_degree"]) return rat_fit @@ -1346,15 +1367,15 @@ def get_approximate_degree_for_rational_fit(self): len(self.apocenters_location)) if approximate_num_orbits <= 5: num_degree, denom_degree = 1, 1 - elif (approximate_num_orbits > 5) and (approximate_num_orbits <= 20): + elif approximate_num_orbits <= 20: num_degree, denom_degree = 2, 2 - elif (approximate_num_orbits > 20) and (approximate_num_orbits <= 50): + elif approximate_num_orbits <= 50: num_degree, denom_degree = 3, 3 - elif (approximate_num_orbits > 50) and (approximate_num_orbits <= 100): + elif approximate_num_orbits <= 100: num_degree, denom_degree = 4, 4 - elif (approximate_num_orbits > 100) and (approximate_num_orbits <= 150): + elif approximate_num_orbits <= 150: num_degree, denom_degree = 5, 5 - elif (approximate_num_orbits > 150) and (approximate_num_orbits <= 200): + elif approximate_num_orbits <= 200: num_degree, denom_degree = 6, 6 else: num_degree, denom_degree = 5 + int(np.log10(approximate_num_orbits)), 5 + int(np.log10(approximate_num_orbits)) @@ -1382,26 +1403,23 @@ def rational_fit_extrema(self, extrema_type="pericenters"): if len(extrema) >= 2: return self.rational_fit(self.t[extrema], self.omega_gw[extrema], - description=extrema_type) + data_name=extrema_type) else: raise Exception( f"Sufficient number of {extrema_type} are not found." " Can not create an interpolant.") - def check_if_domega_dt_is_nonmonotonic(self, t, omega, tol=1.0, description=None): - """Check if the first derivative of omega is nonmonotonic - - The first derivative of interpolant/fit of omega at the extrema - should be monotonically increasing. + def check_if_first_derivative_is_not_strictly_monotonic(self, x, y, tol=1.0, data_name=None): + """Check if the first derivative of data is not strictly monotonic. """ - domega_dt = np.gradient(omega, t[1] - t[0]) - is_nonmonotonic = any(domega_dt[1:]/domega_dt[:-1] < tol) + dy_dx = np.gradient(y, x[1] - x[0]) + is_not_strictly_monotonic = any(dy_dx[1:]/dy_dx[:-1] < tol) # update history if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": - if description in self.rational_fit_nonmonotonicity_history: - self.rational_fit_nonmonotonicity_history[description].update({ - (self.rational_fit_kwargs['num_degree'], self.rational_fit_kwargs['num_degree']): is_nonmonotonic}) - return is_nonmonotonic + if data_name in self.rational_fit_nonmonotonicity_history: + self.rational_fit_nonmonotonicity_history[data_name].update({ + (self.rational_fit_kwargs['num_degree'], self.rational_fit_kwargs['num_degree']): is_not_strictly_monotonic}) + return is_not_strictly_monotonic def check_num_extrema(self, extrema, extrema_type="extrema"): """Check number of extrema. @@ -1816,6 +1834,20 @@ def measure_ecc(self, tref_in=None, fref_in=None): return self.make_return_dict_for_eccentricity_and_mean_anomaly() def build_omega_gw_extrema_interpolants(self): + """Build interpolant through omega_gw extrema. + + Calculating the eccentricity using the gravitational wave frequency + (`omega_gw`), requires building an interpolant through its extrema + (pericenters and apocenters). + + The method for constructing this interpolant can be specified by the + user via `omega_gw_extrema_interpolation_method` in + `extra_kwargs`. Currently, two options are supported: `"spline"` and + `"rational_fit"`. + + Based on the chosen method, an interpolation will be constructed for + `omega_gw` at the pericenter and apocenter points. + """ # Build the interpolants of omega_gw at the extrema if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": self.omega_gw_apocenters_interp = self.interp_extrema("apocenters") @@ -1832,9 +1864,9 @@ def build_omega_gw_extrema_interpolants(self): # Verify the monotonicity of the first derivative of the omega_gw interpolant with spline. if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": # Check if the first derivative of omega_gw at pericenters or apocenters is non-monotonic - if (self.check_if_domega_dt_is_nonmonotonic( + if (self.check_if_first_derivative_is_not_strictly_monotonic( self.t_for_checks, self.omega_gw_pericenters_interp(self.t_for_checks)) or - self.check_if_domega_dt_is_nonmonotonic( + self.check_if_first_derivative_is_not_strictly_monotonic( self.t_for_checks, self.omega_gw_apocenters_interp(self.t_for_checks))): debug_message( "Nonmonotonic time derivative detected in the spline interpolant through extrema. " @@ -2438,16 +2470,15 @@ def compute_orbit_averaged_omega_gw_between_extrema(self, t): def check_monotonicity_of_omega_gw_average(self, omega_gw_average, - description="omega_gw average"): + data_name="omega_gw average"): """Check that omega average is monotonically increasing. Parameters ---------- omega_gw_average : array-like 1d array of omega_gw averages to check for monotonicity. - description : str - String to describe what the the which omega_gw average we are - looking at. + data_name : str + String to describe which omega_gw average we are looking at. """ idx_non_monotonic = np.where( np.diff(omega_gw_average) <= 0)[0] @@ -2503,7 +2534,7 @@ def check_monotonicity_of_omega_gw_average(self, fig.tight_layout() figName = ( "./gwecc_" - f"{self.method}_{description.replace(' ', '_')}.pdf") + f"{self.method}_{data_name.replace(' ', '_')}.pdf") # fig.savefig(figName) self.save_debug_fig(fig, figName) plt.close(fig) @@ -2513,7 +2544,7 @@ def check_monotonicity_of_omega_gw_average(self, "for diagnostic plot use `debug_plots=True` in " "extra_kwargs") raise Exception( - f"{description} are non-monotonic.\n" + f"{data_name} are non-monotonic.\n" f"First non-monotonicity occurs at peak number {first_idx}," f" where omega_gw drops from {omega_gw_average[first_idx]} to" f" {omega_gw_average[first_idx+1]}, a decrease by" From 6f2d3aeccf68ab9d26f644004b753598a92ef573 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Tue, 5 Nov 2024 19:55:43 +0530 Subject: [PATCH 29/34] update -> store --- gw_eccentricity/eccDefinition.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index dd9db313..9a0fdb17 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -353,7 +353,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, self.rational_fit_kwargs["verbose"] = self.debug_level # keep history of rational fit degree and nonmonotonicity of the corresponding fits self.rational_fit_nonmonotonicity_history = {"pericenters": {}, "apocenters": {}} - # Update the degrees used to construct the final fits + # Store degrees used to construct the final fits self.rational_fit_degrees = {"pericenters": None, "apocenters": None} self.debug_plots = self.extra_kwargs["debug_plots"] self.return_zero_if_small_ecc_failure = self.extra_kwargs["return_zero_if_small_ecc_failure"] From 0e89ed4a289bffc0c706a4ac142b80ffd0555eb1 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Wed, 13 Nov 2024 21:41:53 +0530 Subject: [PATCH 30/34] separate interpolation for extrema and non extrema data --- gw_eccentricity/eccDefinition.py | 480 +++++++++++++++++++------------ 1 file changed, 292 insertions(+), 188 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 9a0fdb17..ea117d1b 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -168,18 +168,39 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, extra_kwargs: dict A dictionary of any extra kwargs to be passed. Allowed kwargs are: - spline_kwargs: dict + extrema_interp_kwargs: dict + A dictionary of dictionaries where each key corresponds to a + specific interpolation method from the available + `available_omega_gw_extrema_interpolation_methods`. Each + top-level key maps to a nested dictionary containing keyword + arguments (`kwargs`) specific to that interpolation method. The + inner dictionary is passed as arguments to the respective + `omega_gw_extrema_interpolation_method`. + + Example structure: + { + "method_1": {"param1": value1, "param2": value2}, + "method_2": {"paramA": valueA, "paramB": valueB} + }, + "method_1", "method_2" should be available methods in + `available_omega_gw_extrema_interpolation_methods`. + + currently, the available methods are: + + - "spline": default kwargs are set using + `utils.get_default_spline_kwargs` + - "rational_fit": default kwargs are set using + `utils.get_default_rational_fit_kwargs` + + See under `omega_gw_extrema_interpolation_method` for more + details on these interpolation methods. + + non_extrema_interp_kwargs: dict Dictionary of arguments to be passed to the spline interpolation routine (scipy.interpolate.InterpolatedUnivariateSpline) used to - compute quantities like omega_gw_pericenters(t) and - 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` + interpolate data other than the omega_gw extrema. + Defaults are set using `utils.get_default_spline_kwargs` extrema_finding_kwargs: dict Dictionary of arguments to be passed to the extrema finder, @@ -281,25 +302,27 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, 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: + 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. + - 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, especially for noisy NR data. + - Since it fits through every data point, it may exhibit + oscillatory behavior, particularly near the merger, + especially for noisy NR data. - `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`. + - 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"`. """ @@ -330,33 +353,47 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, extra_kwargs, self.get_default_extra_kwargs(), "extra_kwargs", "eccDefinition.get_default_extra_kwargs()") + self.debug_level = self.extra_kwargs["debug_level"] self.extrema_finding_kwargs = check_kwargs_and_set_defaults( self.extra_kwargs['extrema_finding_kwargs'], self.get_default_extrema_finding_kwargs(min_width_for_extrema), "extrema_finding_kwargs", "eccDefinition.get_default_extrema_finding_kwargs()") - self.spline_kwargs = check_kwargs_and_set_defaults( - self.extra_kwargs["spline_kwargs"], + # set omega extrema interpolation method + self.omega_gw_extrema_interpolation_method = \ + self.extra_kwargs["omega_gw_extrema_interpolation_method"] + # check extrema interpolation method + self.available_omega_gw_extrema_interpolation_methods = \ + self.get_available_omega_gw_extrema_interpolation_methods() + if self.omega_gw_extrema_interpolation_method \ + not in self.available_omega_gw_extrema_interpolation_methods: + raise Exception( + "Unknown omega_gw_extrema_interpolation_method " + f"`{self.omega_gw_extrema_interpolation_method}`. Should be " + f"of {self.available_omega_gw_extrema_interpolation_methods.keys()}") + # set extrema interpolation kwargs + self.extrema_interp_kwargs = check_kwargs_and_set_defaults( + self.extra_kwargs[ + "extrema_interp_kwargs"][ + self.omega_gw_extrema_interpolation_method], + self.get_default_extrema_interp_kwargs(), + "extrema_interp_kwargs", + "eccDefinition.get_default_extrema_interp_kwargs()") + # set other variables required for the omega_gw extrema interpolation + # method + self.set_other_variables_for_extrema_interpolation() + # set spline interpolation kwargs to be used for interpolating data + # other than the extrema. + self.non_extrema_interp_kwargs = check_kwargs_and_set_defaults( + self.extra_kwargs["non_extrema_interp_kwargs"], get_default_spline_kwargs(), - "spline_kwargs", + "non_extrema_interp_kwargs", "utils.get_default_spline_kwargs()") - self.rational_fit_kwargs = check_kwargs_and_set_defaults( - self.extra_kwargs["rational_fit_kwargs"], - get_default_rational_fit_kwargs(), - "rational_fit_kwargs", - "utils.get_default_rational_fit_kwargs()") self.available_averaging_methods \ = self.get_available_omega_gw_averaging_methods() - self.debug_level = self.extra_kwargs["debug_level"] - # set verbose to debug_level. If verbose is True, then it prints information - # of each iteration for rational fits for a given degree. - self.rational_fit_kwargs["verbose"] = self.debug_level - # keep history of rational fit degree and nonmonotonicity of the corresponding fits - self.rational_fit_nonmonotonicity_history = {"pericenters": {}, "apocenters": {}} - # Store degrees used to construct the final fits - self.rational_fit_degrees = {"pericenters": None, "apocenters": None} self.debug_plots = self.extra_kwargs["debug_plots"] - self.return_zero_if_small_ecc_failure = self.extra_kwargs["return_zero_if_small_ecc_failure"] + self.return_zero_if_small_ecc_failure \ + = self.extra_kwargs["return_zero_if_small_ecc_failure"] # check if there are unrecognized keys in the dataDict self.recognized_dataDict_keys = self.get_recognized_dataDict_keys() for kw in dataDict.keys(): @@ -870,12 +907,53 @@ def get_default_extrema_finding_kwargs(self, width): "rel_height": 0.5, "plateau_size": None} return default_extrema_finding_kwargs + + def get_default_extrema_interp_kwargs(self): + """Get default kwargs to be passed to extrema interpolating method. + + Accurate interpolation of omega_gw extrema is crucial for obtaining + monotonic eccentricity evoultion with time. Depending on the + interpolation method, this function returns the default kwargs to be + passed to the interpolating function. + """ + allowed_methods = self.extra_kwargs["extrema_interp_kwargs"].keys() + if self.omega_gw_extrema_interpolation_method == "spline": + kwargs = get_default_spline_kwargs() + elif self.omega_gw_extrema_interpolation_method == "rational_fit": + kwargs = get_default_rational_fit_kwargs() + else: + raise Exception("Unknown omega_gw extrema interpolation method " + f"`{self.omega_gw_extrema_interpolation_method}`. " + f"Allowed methods are {allowed_methods}") + return kwargs + + def set_other_variables_for_extrema_interpolation(self): + """Set other variables required for extrema interpolation. + + Depending on the omega_gw_extrema_interpolation_method, additional + varibles may be required for obtaining good interpolant and/or for + debugging purpose. These variables may vary from method to method. + """ + if self.omega_gw_extrema_interpolation_method == "rational_fit": + # set verbose to debug_level. If verbose is True, then it prints + # information of each iteration for rational fits for a given + # degree. + self.extrema_interp_kwargs["verbose"] = self.debug_level + # keep history of rational fit degree and nonmonotonicity of the + # corresponding fits + self.rational_fit_nonmonotonicity_history = { + "pericenters": {}, "apocenters": {}} + # Store degrees used to construct the final fits + self.rational_fit_degrees = { + "pericenters": None, "apocenters": None} def get_default_extra_kwargs(self): """Defaults for additional kwargs.""" default_extra_kwargs = { - "spline_kwargs": {}, - "rational_fit_kwargs": {}, + "non_extrema_interp_kwargs": {}, + "extrema_interp_kwargs": { + "spline": {}, + "rational_fit": {}}, "extrema_finding_kwargs": {}, # Gets overridden in methods like # eccDefinitionUsingAmplitude "debug_level": 0, @@ -1177,9 +1255,17 @@ def get_good_extrema(self, pericenters, apocenters, pericenters, apocenters) return pericenters, apocenters - def get_interp(self, oldX, oldY, allowExtrapolation=False, - interpolator="spline"): - """Get interpolant. + def get_spline_interpolant_for_non_extrema( + self, oldX, oldY, allowExtrapolation=False, interpolator="spline"): + """Get spline interpolant for data other than the omega_gw extrema. + + Interpolating the omega_gw extrema could be challenging, particularly + near the merger due to the rapid change in the omega_gw values and + limited number of data points. Thus, it requires special treatment. For + other data, this problem is less likely to occur. Therefore, we use + this function to interpolate data other than omega_gw extrema. The + spline kwargs for this function is provided by `extra_kwargs` using the + key "non_extrema_interp_kwargs". A wrapper of utils.get_interpolant with check_kwargs=False. This is to make sure that the checking of kwargs is not performed @@ -1188,46 +1274,33 @@ def get_interp(self, oldX, oldY, allowExtrapolation=False, function without repeating checks. """ return get_interpolant(oldX, oldY, allowExtrapolation, interpolator, - spline_kwargs=self.spline_kwargs, + spline_kwargs=self.non_extrema_interp_kwargs, check_kwargs=False) - def interp(self, newX, oldX, oldY, allowExtrapolation=False, - interpolator="spline"): - """Get interpolated values. + def get_spline_interpolant_for_extrema( + self, oldX, oldY, allowExtrapolation=False, interpolator="spline"): + """Get spline interpolant for omega_gw extrema. - A wrapper of utils.interpolate with check_kwargs=False for - reasons explained in the documentation of get_interp function. + Same as `get_spline_interpolant_for_non_extrema` but uses the + `spline_kwargs` provided via `self.extrema_interp_kwargs`. For more, + see `get_spline_interpolant_for_non_extrema`. """ - return interpolate(newX, oldX, oldY, allowExtrapolation, interpolator, - spline_kwargs=self.spline_kwargs, - check_kwargs=False) - - def interp_extrema(self, extrema_type="pericenters"): - """Build interpolant through extrema. + return get_interpolant(oldX, oldY, allowExtrapolation, interpolator, + spline_kwargs=self.extrema_interp_kwargs, + check_kwargs=False) - parameters: - ----------- - extrema_type: - Either "pericenters" or "apocenters". + def spline_interpolate_non_extrema(self, newX, oldX, oldY, allowExtrapolation=False, + interpolator="spline"): + """Get interpolated values using spline for data other than omega_gw + extrema. - returns: - ------ - Interpolant through extrema + A wrapper of utils.interpolate with check_kwargs=False for reasons + explained in the documentation of + `get_spline_interpolant_for_non_extrema` function. """ - if extrema_type == "pericenters": - extrema = self.pericenters_location - elif extrema_type == "apocenters": - extrema = self.apocenters_location - else: - raise Exception("extrema_type must be either " - "'pericenrers' or 'apocenters'.") - if len(extrema) >= 2: - return self.get_interp(self.t[extrema], - self.omega_gw[extrema]) - else: - raise Exception( - f"Sufficient number of {extrema_type} are not found." - " Can not create an interpolant.") + return interpolate(newX, oldX, oldY, allowExtrapolation, interpolator, + spline_kwargs=self.non_extrema_interp_kwargs, + check_kwargs=False) def get_rational_fit_wrapper(self, x, y): """Get rational fit. @@ -1239,10 +1312,10 @@ def get_rational_fit_wrapper(self, x, y): checks. """ return get_rational_fit(x, y, - rational_fit_kwargs=self.rational_fit_kwargs, + rational_fit_kwargs=self.extrema_interp_kwargs, check_kwargs=False) - def rational_fit(self, x, y, data_name=None): + def get_rational_fit_for_extrema(self, x, y, data_name=None): """Get rational fit with adaptive numerator and denominator degree. This function ensures that the rational fit uses the optimal degree for @@ -1254,70 +1327,81 @@ def rational_fit(self, x, y, data_name=None): the previous valid monotonic fit. The initial degrees for the rational fit can be specified through - `rational_fit_kwargs` in `extra_kwargs`. Default values are provided by - `get_default_rational_fit_kwargs`, where both degrees are set to - `None`. If both degrees remain `None`, appropriate starting values are - determined using `self.get_approximate_degree_for_rational_fit`. + `extrema_interp_kwargs` in `extra_kwargs` with the key "rational_fit". + (See `extrema_interp_kwargs` under `extra_kwargs` for more details) + Default values are provided by `get_default_rational_fit_kwargs`, where + both degrees are set to `None`. If both degrees remain `None`, + appropriate starting values are determined using + `self.get_approximate_degree_for_rational_fit`. """ # make sure that data_name is not None. A data_name is needed to # update the optimal values of the numerator and denominator degrees # used to build the final rational fit if data_name is None: - raise Exception("Please provide a data_name for which to build a " - "rational fit. `data_name` can not be None. For example, " - "it can be 'apocenters', 'pericenters' or 'omega_gw_average'.") - # Set initial degrees for the rational fit if not already specified. - # If either `num_degree` or `denom_degree` is specified, use that value - # to set the other. If both are unspecified, we use the following - # procedure to set appropriate values for `num_degree` and `denom_degree`. + raise Exception( + "Please provide a data_name for which to build a " + "rational fit. `data_name` can not be None. For example, " + "it can be 'apocenters', 'pericenters' or 'omega_gw_average'.") + # Set initial degrees for the rational fit if not already specified. If + # either `num_degree` or `denom_degree` is specified, use that value to + # set the other. If both are unspecified, we use the following + # procedure to set appropriate values for `num_degree` and + # `denom_degree`. - # The optimal degrees for the numerator and denominator can vary with - # the eccentricity and duration of the waveform. To estimate these values, - # we first count the approximate number of orbits, then use a predefined - # mapping from orbit count to optimal degrees, provided by + # The optimal degrees for the numerator and denominator can vary with + # the eccentricity and duration of the waveform. To estimate these + # values, we first count the approximate number of orbits, then use a + # predefined mapping from orbit count to optimal degrees, provided by # `get_approximate_degree_for_rational_fit`. # While this provides a good initial estimate, it may not be optimal, # as the ideal degree can change based on the eccentricity which we do # not know apriori. The true optimal values are determined through # further iterations. - if self.rational_fit_kwargs["num_degree"] is None and self.rational_fit_kwargs["denom_degree"] is None: - self.rational_fit_kwargs["num_degree"], self.rational_fit_kwargs["denom_degree"] = \ + if self.extrema_interp_kwargs["num_degree"] is None \ + and self.extrema_interp_kwargs["denom_degree"] is None: + self.extrema_interp_kwargs["num_degree"], self.extrema_interp_kwargs["denom_degree"] = \ self.get_approximate_degree_for_rational_fit() - elif self.rational_fit_kwargs["num_degree"] is None: - self.rational_fit_kwargs["num_degree"] = self.rational_fit_kwargs["denom_degree"] - elif self.rational_fit_kwargs["denom_degree"] is None: - self.rational_fit_kwargs["denom_degree"] = self.rational_fit_kwargs["num_degree"] + elif self.extrema_interp_kwargs["num_degree"] is None: + self.extrema_interp_kwargs["num_degree"] = self.extrema_interp_kwargs["denom_degree"] + elif self.extrema_interp_kwargs["denom_degree"] is None: + self.extrema_interp_kwargs["denom_degree"] = self.extrema_interp_kwargs["num_degree"] rat_fit = self.get_rational_fit_wrapper(x, y) x_test = np.arange(x[0], x[-1], self.t[1] - self.t[0]) - # save the degrees for checks at each step of iterations for finding the optimal degrees - old_num_degree = self.rational_fit_kwargs["num_degree"] - old_denom_degree = self.rational_fit_kwargs["denom_degree"] + # save the degrees for checks at each step of iterations for finding + # the optimal degrees + old_num_degree = self.extrema_interp_kwargs["num_degree"] + old_denom_degree = self.extrema_interp_kwargs["denom_degree"] # Check for nonmonotonicity and lower degrees if needed fit_is_nonmonotonic = self.check_if_first_derivative_is_not_strictly_monotonic( x_test, rat_fit(x_test), 1.0, data_name) while fit_is_nonmonotonic: - if self.rational_fit_kwargs["num_degree"] > 1: - self.rational_fit_kwargs["num_degree"] -= 1 - if self.rational_fit_kwargs["denom_degree"] > 1: - self.rational_fit_kwargs["denom_degree"] -= 1 - if self.rational_fit_kwargs["num_degree"] == 1 and self.rational_fit_kwargs["denom_degree"] == 1: - debug_message("Both numerator and denominator degrees are equal to 1 " - "and cannot be lowered further.", - debug_level=self.debug_level, important=False) + if self.extrema_interp_kwargs["num_degree"] > 1: + self.extrema_interp_kwargs["num_degree"] -= 1 + if self.extrema_interp_kwargs["denom_degree"] > 1: + self.extrema_interp_kwargs["denom_degree"] -= 1 + if self.extrema_interp_kwargs["num_degree"] == 1 \ + and self.extrema_interp_kwargs["denom_degree"] == 1: + debug_message( + "Both numerator and denominator degrees are equal to 1 " + "and cannot be lowered further.", + debug_level=self.debug_level, important=False) break - debug_message(f"Lowering degrees to num_degree={self.rational_fit_kwargs['num_degree']}, " - f"denom_degree={self.rational_fit_kwargs['denom_degree']} and retrying.", - debug_level=self.debug_level, important=False) + debug_message( + "Lowering degrees to " + f"num_degree={self.extrema_interp_kwargs['num_degree']}, " + f"denom_degree={self.extrema_interp_kwargs['denom_degree']} " + "and retrying.", debug_level=self.debug_level, important=False) # build new fit and check monotonicity rat_fit = self.get_rational_fit_wrapper(x, y) fit_is_nonmonotonic = self.check_if_first_derivative_is_not_strictly_monotonic( x_test, rat_fit(x_test), 1.0, data_name) - # If fit with initial degree is monotonic, try increasing the degree for a better fit + # If fit with initial degree is monotonic, try increasing the degree + # for a better fit if not fit_is_nonmonotonic: last_monotonic_rat_fit = rat_fit # Track last monotonic fit last_monotonic_num_degree = old_num_degree @@ -1325,22 +1409,27 @@ def rational_fit(self, x, y, data_name=None): while not fit_is_nonmonotonic: # Increase the degrees for both numerator and denominator - new_num_degree = self.rational_fit_kwargs["num_degree"] + 1 - new_denom_degree = self.rational_fit_kwargs["denom_degree"] + 1 - self.rational_fit_kwargs["num_degree"], self.rational_fit_kwargs["denom_degree"] \ + new_num_degree = self.extrema_interp_kwargs["num_degree"] + 1 + new_denom_degree = self.extrema_interp_kwargs["denom_degree"] + 1 + self.extrema_interp_kwargs["num_degree"], self.extrema_interp_kwargs["denom_degree"] \ = new_num_degree, new_denom_degree # build new fit and check monotonicity new_rat_fit = self.get_rational_fit_wrapper(x, y) - fit_is_nonmonotonic = self.check_if_first_derivative_is_not_strictly_monotonic( - x_test, new_rat_fit(x_test), 1.0, data_name) + fit_is_nonmonotonic \ + = self.check_if_first_derivative_is_not_strictly_monotonic( + x_test, new_rat_fit(x_test), 1.0, data_name) if fit_is_nonmonotonic: - # Revert to previous fit and degrees if nonmonotonicity is detected - debug_message(f"Increasing degrees caused nonmonotonicity. Reverting to " - f"last monotonic fit with num_degree={last_monotonic_num_degree} " - f"and denom_degree={last_monotonic_denom_degree}.", - debug_level=self.debug_level, important=False) - self.rational_fit_kwargs["num_degree"] = last_monotonic_num_degree - self.rational_fit_kwargs["denom_degree"] = last_monotonic_denom_degree + # Revert to previous fit and degrees if nonmonotonicity is + # detected + debug_message( + "Increasing degrees caused nonmonotonicity. Reverting to " + f"last monotonic fit with num_degree={last_monotonic_num_degree} " + f"and denom_degree={last_monotonic_denom_degree}.", + debug_level=self.debug_level, important=False) + self.extrema_interp_kwargs["num_degree"] \ + = last_monotonic_num_degree + self.extrema_interp_kwargs["denom_degree"] \ + = last_monotonic_denom_degree rat_fit = last_monotonic_rat_fit break @@ -1350,15 +1439,16 @@ def rational_fit(self, x, y, data_name=None): rat_fit = new_rat_fit # update final degrees used to build the fit - self.rational_fit_degrees[data_name] = (self.rational_fit_kwargs["num_degree"], - self.rational_fit_kwargs["denom_degree"]) + self.rational_fit_degrees[data_name] = ( + self.extrema_interp_kwargs["num_degree"], + self.extrema_interp_kwargs["denom_degree"]) return rat_fit def get_approximate_degree_for_rational_fit(self): """Get approximate degree based on the number of extrema. - Assign degree based on number of extrema found. The degree is increased as - the number of extrema increases. + Assign degree based on number of extrema found. The degree is increased + as the number of extrema increases. """ # TODO: Optimize this. # assign degree based on the number of extrema if user provided @@ -1378,11 +1468,49 @@ def get_approximate_degree_for_rational_fit(self): elif approximate_num_orbits <= 200: num_degree, denom_degree = 6, 6 else: - num_degree, denom_degree = 5 + int(np.log10(approximate_num_orbits)), 5 + int(np.log10(approximate_num_orbits)) + num_degree, denom_degree \ + = (5 + int(np.log10(approximate_num_orbits)), + 5 + int(np.log10(approximate_num_orbits))) return num_degree, denom_degree - def rational_fit_extrema(self, extrema_type="pericenters"): - """Build rational fit through extrema. + def check_if_first_derivative_is_not_strictly_monotonic( + self, x, y, tol=1.0, data_name=None): + """Check if the first derivative of data is not strictly monotonic. + """ + dy_dx = np.gradient(y, x[1] - x[0]) + is_not_strictly_monotonic = any(dy_dx[1:]/dy_dx[:-1] < tol) + # update history + if self.extra_kwargs["omega_gw_extrema_interpolation_method"] \ + == "rational_fit": + if data_name in self.rational_fit_nonmonotonicity_history: + self.rational_fit_nonmonotonicity_history[data_name].update( + {(self.extrema_interp_kwargs['num_degree'], + self.extrema_interp_kwargs['num_degree']): + is_not_strictly_monotonic}) + return is_not_strictly_monotonic + + def get_available_omega_gw_extrema_interpolation_methods(self): + """Return available omega_gw extrema interpolation methods.""" + available_methods = { + "spline": self.get_spline_interpolant_for_extrema, + "rational_fit": self.get_rational_fit_for_extrema, + } + return available_methods + + def get_omega_gw_extrema_interpolant(self, extrema_type="pericenters"): + """Build interpolant through omega_gw extrema. + + Calculating the eccentricity using the gravitational wave frequency + (`omega_gw`), requires building an interpolant through its extrema + (pericenters and apocenters). + + The method for constructing this interpolant can be specified by the + user via `omega_gw_extrema_interpolation_method` in + `extra_kwargs`. Currently, two options are supported: `"spline"` and + `"rational_fit"`. + + Based on the chosen method, an interpolation will be constructed for + `omega_gw` at the pericenter/apocenter points. parameters: ----------- @@ -1391,7 +1519,7 @@ def rational_fit_extrema(self, extrema_type="pericenters"): returns: ------ - Rational fit through extrema + Interpolant through extrema """ if extrema_type == "pericenters": extrema = self.pericenters_location @@ -1401,26 +1529,20 @@ def rational_fit_extrema(self, extrema_type="pericenters"): raise Exception("extrema_type must be either " "'pericenrers' or 'apocenters'.") if len(extrema) >= 2: - return self.rational_fit(self.t[extrema], - self.omega_gw[extrema], - data_name=extrema_type) + method = self.available_omega_gw_extrema_interpolation_methods[ + self.omega_gw_extrema_interpolation_method] + if self.omega_gw_extrema_interpolation_method == "rational_fit": + # rational fit method `get_rational_fit_for_extrema` takes an + # additional argument + return method(self.t[extrema], self.omega_gw[extrema], + extrema_type) + else: + return method(self.t[extrema], self.omega_gw[extrema]) else: raise Exception( f"Sufficient number of {extrema_type} are not found." " Can not create an interpolant.") - def check_if_first_derivative_is_not_strictly_monotonic(self, x, y, tol=1.0, data_name=None): - """Check if the first derivative of data is not strictly monotonic. - """ - dy_dx = np.gradient(y, x[1] - x[0]) - is_not_strictly_monotonic = any(dy_dx[1:]/dy_dx[:-1] < tol) - # update history - if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": - if data_name in self.rational_fit_nonmonotonicity_history: - self.rational_fit_nonmonotonicity_history[data_name].update({ - (self.rational_fit_kwargs['num_degree'], self.rational_fit_kwargs['num_degree']): is_not_strictly_monotonic}) - return is_not_strictly_monotonic - def check_num_extrema(self, extrema, extrema_type="extrema"): """Check number of extrema. @@ -1810,7 +1932,10 @@ def measure_ecc(self, tref_in=None, fref_in=None): raise Exception("Reference time must be within two pericenters.") # Build omega_gw extrema interpolants - self.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") + # check monotonicity of the interpolants + self.check_omega_gw_extrema_interpolants() # compute eccentricity at self.tref_out self.eccentricity = self.compute_eccentricity(self.tref_out) # Compute mean anomaly at tref_out @@ -1833,45 +1958,24 @@ def measure_ecc(self, tref_in=None, fref_in=None): # frequency where these are measured. return self.make_return_dict_for_eccentricity_and_mean_anomaly() - def build_omega_gw_extrema_interpolants(self): - """Build interpolant through omega_gw extrema. - - Calculating the eccentricity using the gravitational wave frequency - (`omega_gw`), requires building an interpolant through its extrema - (pericenters and apocenters). - - The method for constructing this interpolant can be specified by the - user via `omega_gw_extrema_interpolation_method` in - `extra_kwargs`. Currently, two options are supported: `"spline"` and - `"rational_fit"`. - - Based on the chosen method, an interpolation will be constructed for - `omega_gw` at the pericenter and apocenter points. + def check_omega_gw_extrema_interpolants(self): """ - # Build the interpolants of omega_gw at the extrema - if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": - self.omega_gw_apocenters_interp = self.interp_extrema("apocenters") - self.omega_gw_pericenters_interp = self.interp_extrema("pericenters") - elif self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "rational_fit": - self.omega_gw_apocenters_interp = self.rational_fit_extrema("apocenters") - self.omega_gw_pericenters_interp = self.rational_fit_extrema("pericenters") - else: - raise Exception( - f"Unknown method {self.extra_kwargs['omega_gw_extrema_interpolation_method']} " - "for `omega_gw_extrema_interpolation_method`. " - "Must be one of `spline` or `rational_fit`.") - - # Verify the monotonicity of the first derivative of the omega_gw interpolant with spline. + check monotonicity of the omega_gw extrema interpolants. + """ + # Verify the monotonicity of the first derivative of the omega_gw + # interpolant with spline. if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": - # Check if the first derivative of omega_gw at pericenters or apocenters is non-monotonic + # Check if the first derivative of omega_gw at pericenters or + # apocenters is non-monotonic if (self.check_if_first_derivative_is_not_strictly_monotonic( self.t_for_checks, self.omega_gw_pericenters_interp(self.t_for_checks)) or self.check_if_first_derivative_is_not_strictly_monotonic( self.t_for_checks, self.omega_gw_apocenters_interp(self.t_for_checks))): debug_message( - "Nonmonotonic time derivative detected in the spline interpolant through extrema. " - "Using rational fit by setting 'omega_gw_extrema_interpolation_method' to 'rational_fit' " - "may provide better result.", + "Nonmonotonic time derivative detected in the spline " + "interpolant through extrema. Using rational fit by " + "setting 'omega_gw_extrema_interpolation_method' to " + "'rational_fit' may provide better result.", debug_level=self.debug_level, important=True) def set_eccentricity_and_mean_anomaly_to_zero(self): @@ -2005,8 +2109,8 @@ def derivative_of_eccentricity(self, t, n=1): self.t_for_checks) if self.ecc_interp is None: - self.ecc_interp = self.get_interp(self.t_for_checks, - self.ecc_for_checks) + self.ecc_interp = self.get_spline_interpolant_for_non_extrema( + self.t_for_checks, self.ecc_for_checks) # Get derivative of ecc(t) using spline return self.ecc_interp.derivative(n=n)(t) @@ -2328,11 +2432,11 @@ def compute_res_amp_gw_and_res_omega_gw(self): # residual quantities can be computed. Above, we check that this # extrapolation does not happen before t_merger, which is where # eccentricity is normally measured. - self.amp_gw_zeroecc_interp = self.interp( + self.amp_gw_zeroecc_interp = self.spline_interpolate_non_extrema( self.t, self.t_zeroecc_shifted, self.amp_gw_zeroecc, allowExtrapolation=True) self.res_amp_gw = self.amp_gw - self.amp_gw_zeroecc_interp - self.omega_gw_zeroecc_interp = self.interp( + self.omega_gw_zeroecc_interp = self.spline_interpolate_non_extrema( self.t, self.t_zeroecc_shifted, self.omega_gw_zeroecc, allowExtrapolation=True) self.res_omega_gw = (self.omega_gw - self.omega_gw_zeroecc_interp) @@ -2465,7 +2569,7 @@ def compute_orbit_averaged_omega_gw_between_extrema(self, t): orbit_averaged_omega_gw, "omega_gw averaged [apocenter to apocenter] and " "[pericenter to pericenter]") - return self.interp( + return self.spline_interpolate_non_extrema( t, self.t_for_orbit_averaged_omega_gw, orbit_averaged_omega_gw) def check_monotonicity_of_omega_gw_average(self, @@ -2569,7 +2673,7 @@ def compute_mean_of_extrema_interpolants(self, t): def compute_omega_gw_zeroecc(self, t): """Find omega_gw from zeroecc data.""" - return self.interp( + return self.spline_interpolate_non_extrema( t, self.t_zeroecc_shifted, self.omega_gw_zeroecc) def get_available_omega_gw_averaging_methods(self): @@ -2721,7 +2825,7 @@ def compute_tref_in_and_fref_out_from_fref_in(self, fref_in): self.omega_gw_average, "Interpolated omega_gw_average") # Get tref_in using interpolation - tref_in = self.interp(fref_out, + tref_in = self.spline_interpolate_non_extrema(fref_out, self.omega_gw_average/(2 * np.pi), self.t_for_omega_gw_average) # check if tref_in is monotonically increasing From d41645e790a07273d3fc06e858ffcef1b3a76e68 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Wed, 13 Nov 2024 21:54:35 +0530 Subject: [PATCH 31/34] copy over to gw_eccentricity --- gw_eccentricity/eccDefinition.py | 4 +-- gw_eccentricity/gw_eccentricity.py | 40 +++++++++++++++++++++++------- 2 files changed, 33 insertions(+), 11 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index ea117d1b..1c581507 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -170,7 +170,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, are: extrema_interp_kwargs: dict A dictionary of dictionaries where each key corresponds to a - specific interpolation method from the available + specific interpolation method from the `available_omega_gw_extrema_interpolation_methods`. Each top-level key maps to a nested dictionary containing keyword arguments (`kwargs`) specific to that interpolation method. The @@ -182,7 +182,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, "method_1": {"param1": value1, "param2": value2}, "method_2": {"paramA": valueA, "paramB": valueB} }, - "method_1", "method_2" should be available methods in + "method_1", "method_2" should be methods from `available_omega_gw_extrema_interpolation_methods`. currently, the available methods are: diff --git a/gw_eccentricity/gw_eccentricity.py b/gw_eccentricity/gw_eccentricity.py index 27d8969a..4a5d0182 100644 --- a/gw_eccentricity/gw_eccentricity.py +++ b/gw_eccentricity/gw_eccentricity.py @@ -307,17 +307,39 @@ def measure_eccentricity(tref_in=None, Default is False which implies the system to be nonprecessing. extra_kwargs: A dict of any extra kwargs to be passed. Allowed kwargs are: - spline_kwargs: + extrema_interp_kwargs: dict + A dictionary of dictionaries where each key corresponds to a + specific interpolation method from the + `available_omega_gw_extrema_interpolation_methods`. Each top-level + key maps to a nested dictionary containing keyword arguments + (`kwargs`) specific to that interpolation method. The inner + dictionary is passed as arguments to the respective + `omega_gw_extrema_interpolation_method`. + + Example structure: + { + "method_1": {"param1": value1, "param2": value2}, + "method_2": {"paramA": valueA, "paramB": valueB} + }, + "method_1", "method_2" should be methods from + `available_omega_gw_extrema_interpolation_methods`. + + currently, the available methods are: + + - "spline": default kwargs are set using + `utils.get_default_spline_kwargs` + - "rational_fit": default kwargs are set using + `utils.get_default_rational_fit_kwargs` + + See under `omega_gw_extrema_interpolation_method` for more details + on these interpolation methods. + + non_extrema_interp_kwargs: dict Dictionary of arguments to be passed to the spline interpolation routine (scipy.interpolate.InterpolatedUnivariateSpline) used to - compute quantities like omega_gw_pericenters(t) and - 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` + interpolate data other than the omega_gw extrema. + + Defaults are set using `utils.get_default_spline_kwargs` extrema_finding_kwargs: Dictionary of arguments to be passed to the extrema finder, From ef8aaa3f1e3b11a89808865bb5962b30159fc671 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Thu, 14 Nov 2024 09:48:34 +0530 Subject: [PATCH 32/34] add checks for only one key in extrema_interp_kwargs --- gw_eccentricity/eccDefinition.py | 69 ++++++++++++++++++++++-------- gw_eccentricity/gw_eccentricity.py | 41 +++++++----------- 2 files changed, 68 insertions(+), 42 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 1c581507..00c70bac 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -8,6 +8,7 @@ import numpy as np import matplotlib.pyplot as plt +from copy import deepcopy from .utils import peak_time_via_quadratic_fit, check_kwargs_and_set_defaults from .utils import amplitude_using_all_modes from .utils import time_deriv_4thOrder @@ -169,32 +170,23 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, A dictionary of any extra kwargs to be passed. Allowed kwargs are: extrema_interp_kwargs: dict - A dictionary of dictionaries where each key corresponds to a - specific interpolation method from the - `available_omega_gw_extrema_interpolation_methods`. Each - top-level key maps to a nested dictionary containing keyword - arguments (`kwargs`) specific to that interpolation method. The - inner dictionary is passed as arguments to the respective - `omega_gw_extrema_interpolation_method`. + A dictionary with a single key matching the current + `omega_gw_extrema_interpolation_method`. See under + `omega_gw_extrema_interpolation_method` for more details on + possible values of interpolation methods. Example structure: { - "method_1": {"param1": value1, "param2": value2}, - "method_2": {"paramA": valueA, "paramB": valueB} - }, - "method_1", "method_2" should be methods from - `available_omega_gw_extrema_interpolation_methods`. + "method": {"param1": value1, "param2": value2}, + } - currently, the available methods are: + currently, the available options for "method" are: - "spline": default kwargs are set using `utils.get_default_spline_kwargs` - "rational_fit": default kwargs are set using `utils.get_default_rational_fit_kwargs` - See under `omega_gw_extrema_interpolation_method` for more - details on these interpolation methods. - non_extrema_interp_kwargs: dict Dictionary of arguments to be passed to the spline interpolation routine @@ -350,7 +342,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, = self.get_amp_phase_omega_gw(self.dataDict) # Sanity check various kwargs and set default values self.extra_kwargs = check_kwargs_and_set_defaults( - extra_kwargs, self.get_default_extra_kwargs(), + deepcopy(extra_kwargs), self.get_default_extra_kwargs(), "extra_kwargs", "eccDefinition.get_default_extra_kwargs()") self.debug_level = self.extra_kwargs["debug_level"] @@ -371,6 +363,8 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, "Unknown omega_gw_extrema_interpolation_method " f"`{self.omega_gw_extrema_interpolation_method}`. Should be " f"of {self.available_omega_gw_extrema_interpolation_methods.keys()}") + # check extrema interp kwargs + self.check_extrema_interp_kwargs(extra_kwargs) # set extrema interpolation kwargs self.extrema_interp_kwargs = check_kwargs_and_set_defaults( self.extra_kwargs[ @@ -926,6 +920,47 @@ def get_default_extrema_interp_kwargs(self): f"`{self.omega_gw_extrema_interpolation_method}`. " f"Allowed methods are {allowed_methods}") return kwargs + + def check_extrema_interp_kwargs(self, extra_kwargs): + """Check extrema interp kwargs provided in extra_kwargs. + + The `extrema_interp_kwargs` should have exactly one + key matching the current `omega_gw_extrema_interpolation_method`. + """ + + if "extrema_interp_kwargs" in extra_kwargs: + # check that kwargs for only one extrema interpolation method is + # provided. + common_message = ( + "The input `extrema_interp_kwargs` must contain a single " + "key, matching the current " + "`omega_gw_extrema_interpolation_method` " + f"'{self.omega_gw_extrema_interpolation_method}', " + "with a dictionary of kwargs for " + f"'{self.omega_gw_extrema_interpolation_method}' " + "as its value.") + extrema_interp_keys \ + = list(extra_kwargs["extrema_interp_kwargs"].keys()) + if len(extrema_interp_keys) == 0: + raise Exception( + "Dictionay provided by 'extrema_interp_kwargs' " + " in 'extra_kwargs' can not be empty.\n" + f"{common_message}" + ) + # check that the key in extrema_interp_kwargs matches + # `omega_gw_extrema_interpolation_method` + if len(extrema_interp_keys) == 1 and ( + extrema_interp_keys[0] + != self.omega_gw_extrema_interpolation_method): + raise Exception( + "`omega_gw_extrema_interpolation_method` is " + f"{self.omega_gw_extrema_interpolation_method} but the " + "kwargs in 'extrema_interp_kwargs' via 'extra_kwargs' is " + f"for {extrema_interp_keys[0]}.\n" + f"{common_message}" + ) + if len(extrema_interp_keys) > 1: + raise Exception(f"{common_message}") def set_other_variables_for_extrema_interpolation(self): """Set other variables required for extrema interpolation. diff --git a/gw_eccentricity/gw_eccentricity.py b/gw_eccentricity/gw_eccentricity.py index 4a5d0182..13292cc0 100644 --- a/gw_eccentricity/gw_eccentricity.py +++ b/gw_eccentricity/gw_eccentricity.py @@ -308,31 +308,22 @@ def measure_eccentricity(tref_in=None, extra_kwargs: A dict of any extra kwargs to be passed. Allowed kwargs are: extrema_interp_kwargs: dict - A dictionary of dictionaries where each key corresponds to a - specific interpolation method from the - `available_omega_gw_extrema_interpolation_methods`. Each top-level - key maps to a nested dictionary containing keyword arguments - (`kwargs`) specific to that interpolation method. The inner - dictionary is passed as arguments to the respective - `omega_gw_extrema_interpolation_method`. - - Example structure: - { - "method_1": {"param1": value1, "param2": value2}, - "method_2": {"paramA": valueA, "paramB": valueB} - }, - "method_1", "method_2" should be methods from - `available_omega_gw_extrema_interpolation_methods`. - - currently, the available methods are: - - - "spline": default kwargs are set using - `utils.get_default_spline_kwargs` - - "rational_fit": default kwargs are set using - `utils.get_default_rational_fit_kwargs` - - See under `omega_gw_extrema_interpolation_method` for more details - on these interpolation methods. + A dictionary with a single key matching the current + `omega_gw_extrema_interpolation_method`. See under + `omega_gw_extrema_interpolation_method` for more details on + possible values of interpolation methods. + + Example structure: + { + "method": {"param1": value1, "param2": value2}, + } + + currently, the available options for "method" are: + + - "spline": default kwargs are set using + `utils.get_default_spline_kwargs` + - "rational_fit": default kwargs are set using + `utils.get_default_rational_fit_kwargs` non_extrema_interp_kwargs: dict Dictionary of arguments to be passed to the spline interpolation From d158f311b99cc5cbdf84953a8b66fa1d33827c53 Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Thu, 14 Nov 2024 10:03:11 +0530 Subject: [PATCH 33/34] do only extra_kwargs not None --- gw_eccentricity/eccDefinition.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 00c70bac..2f01761b 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -927,8 +927,7 @@ def check_extrema_interp_kwargs(self, extra_kwargs): The `extrema_interp_kwargs` should have exactly one key matching the current `omega_gw_extrema_interpolation_method`. """ - - if "extrema_interp_kwargs" in extra_kwargs: + if extra_kwargs is not None and "extrema_interp_kwargs" in extra_kwargs: # check that kwargs for only one extrema interpolation method is # provided. common_message = ( From c916d7032b97c13b1219c401301dc4d141e5354f Mon Sep 17 00:00:00 2001 From: md-arif-shaikh Date: Wed, 27 Nov 2024 12:01:12 +0530 Subject: [PATCH 34/34] rename interp kwargs --- gw_eccentricity/eccDefinition.py | 284 +++++++++++++---------------- gw_eccentricity/gw_eccentricity.py | 4 +- 2 files changed, 131 insertions(+), 157 deletions(-) diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py index 7452d9d7..f1bf023f 100644 --- a/gw_eccentricity/eccDefinition.py +++ b/gw_eccentricity/eccDefinition.py @@ -202,7 +202,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, extra_kwargs: dict A dictionary of any extra kwargs to be passed. Allowed kwargs are: - extrema_interp_kwargs: dict + special_interp_kwargs_for_extrema: dict A dictionary with a single key matching the current `omega_gw_extrema_interpolation_method`. See under `omega_gw_extrema_interpolation_method` for more details on @@ -220,7 +220,7 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, - "rational_fit": default kwargs are set using `utils.get_default_rational_fit_kwargs` - non_extrema_interp_kwargs: dict + general_interp_kwargs: dict Dictionary of arguments to be passed to the spline interpolation routine (scipy.interpolate.InterpolatedUnivariateSpline) used to @@ -404,24 +404,24 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2, f"`{self.omega_gw_extrema_interpolation_method}`. Should be " f"of {self.available_omega_gw_extrema_interpolation_methods.keys()}") # check extrema interp kwargs - self.check_extrema_interp_kwargs(extra_kwargs) + self.check_special_interp_kwargs_for_extrema(extra_kwargs) # set extrema interpolation kwargs - self.extrema_interp_kwargs = check_kwargs_and_set_defaults( + self.special_interp_kwargs_for_extrema = check_kwargs_and_set_defaults( self.extra_kwargs[ - "extrema_interp_kwargs"][ + "special_interp_kwargs_for_extrema"][ self.omega_gw_extrema_interpolation_method], - self.get_default_extrema_interp_kwargs(), - "extrema_interp_kwargs", - "eccDefinition.get_default_extrema_interp_kwargs()") + self.get_default_special_interp_kwargs_for_extrema(), + "special_interp_kwargs_for_extrema", + "eccDefinition.get_default_special_interp_kwargs_for_extrema()") # set other variables required for the omega_gw extrema interpolation # method self.set_other_variables_for_extrema_interpolation() # set spline interpolation kwargs to be used for interpolating data # other than the extrema. - self.non_extrema_interp_kwargs = check_kwargs_and_set_defaults( - self.extra_kwargs["non_extrema_interp_kwargs"], + self.general_interp_kwargs = check_kwargs_and_set_defaults( + self.extra_kwargs["general_interp_kwargs"], get_default_spline_kwargs(), - "non_extrema_interp_kwargs", + "general_interp_kwargs", "utils.get_default_spline_kwargs()") self.available_averaging_methods \ = self.get_available_omega_gw_averaging_methods() @@ -1054,7 +1054,7 @@ def get_default_extrema_finding_kwargs(self, width): "plateau_size": None} return default_extrema_finding_kwargs - def get_default_extrema_interp_kwargs(self): + def get_default_special_interp_kwargs_for_extrema(self): """Get default kwargs to be passed to extrema interpolating method. Accurate interpolation of omega_gw extrema is crucial for obtaining @@ -1062,7 +1062,7 @@ def get_default_extrema_interp_kwargs(self): interpolation method, this function returns the default kwargs to be passed to the interpolating function. """ - allowed_methods = self.extra_kwargs["extrema_interp_kwargs"].keys() + allowed_methods = self.extra_kwargs["special_interp_kwargs_for_extrema"].keys() if self.omega_gw_extrema_interpolation_method == "spline": kwargs = get_default_spline_kwargs() elif self.omega_gw_extrema_interpolation_method == "rational_fit": @@ -1073,32 +1073,33 @@ def get_default_extrema_interp_kwargs(self): f"Allowed methods are {allowed_methods}") return kwargs - def check_extrema_interp_kwargs(self, extra_kwargs): + def check_special_interp_kwargs_for_extrema(self, extra_kwargs): """Check extrema interp kwargs provided in extra_kwargs. - The `extrema_interp_kwargs` should have exactly one + The `special_interp_kwargs_for_extrema` should have exactly one key matching the current `omega_gw_extrema_interpolation_method`. """ - if extra_kwargs is not None and "extrema_interp_kwargs" in extra_kwargs: + if (extra_kwargs is not None + and "special_interp_kwargs_for_extrema" in extra_kwargs): # check that kwargs for only one extrema interpolation method is # provided. common_message = ( - "The input `extrema_interp_kwargs` must contain a single " - "key, matching the current " + "The input `special_interp_kwargs_for_extrema` must contain " + "a single key, matching the current " "`omega_gw_extrema_interpolation_method` " f"'{self.omega_gw_extrema_interpolation_method}', " "with a dictionary of kwargs for " f"'{self.omega_gw_extrema_interpolation_method}' " "as its value.") extrema_interp_keys \ - = list(extra_kwargs["extrema_interp_kwargs"].keys()) + = list(extra_kwargs["special_interp_kwargs_for_extrema"].keys()) if len(extrema_interp_keys) == 0: raise Exception( - "Dictionay provided by 'extrema_interp_kwargs' " + "Dictionay provided via 'special_interp_kwargs_for_extrema'" " in 'extra_kwargs' can not be empty.\n" f"{common_message}" ) - # check that the key in extrema_interp_kwargs matches + # check that the key in special_interp_kwargs_for_extrema matches # `omega_gw_extrema_interpolation_method` if len(extrema_interp_keys) == 1 and ( extrema_interp_keys[0] @@ -1106,8 +1107,8 @@ def check_extrema_interp_kwargs(self, extra_kwargs): raise Exception( "`omega_gw_extrema_interpolation_method` is " f"{self.omega_gw_extrema_interpolation_method} but the " - "kwargs in 'extrema_interp_kwargs' via 'extra_kwargs' is " - f"for {extrema_interp_keys[0]}.\n" + "kwargs in 'special_interp_kwargs_for_extrema' via " + f"'extra_kwargs' is for {extrema_interp_keys[0]}.\n" f"{common_message}" ) if len(extrema_interp_keys) > 1: @@ -1124,7 +1125,7 @@ def set_other_variables_for_extrema_interpolation(self): # set verbose to debug_level. If verbose is True, then it prints # information of each iteration for rational fits for a given # degree. - self.extrema_interp_kwargs["verbose"] = self.debug_level + self.special_interp_kwargs_for_extrema["verbose"] = self.debug_level # keep history of rational fit degree and nonmonotonicity of the # corresponding fits self.rational_fit_nonmonotonicity_history = { @@ -1136,8 +1137,8 @@ def set_other_variables_for_extrema_interpolation(self): def get_default_extra_kwargs(self): """Defaults for additional kwargs.""" default_extra_kwargs = { - "non_extrema_interp_kwargs": {}, - "extrema_interp_kwargs": { + "general_interp_kwargs": {}, + "special_interp_kwargs_for_extrema": { "spline": {}, "rational_fit": {}}, "extrema_finding_kwargs": {}, # Gets overridden in methods like @@ -1441,66 +1442,6 @@ def get_good_extrema(self, pericenters, apocenters, pericenters, apocenters) return pericenters, apocenters - def get_spline_interpolant_for_non_extrema( - self, oldX, oldY, allowExtrapolation=False, interpolator="spline"): - """Get spline interpolant for data other than the omega_gw extrema. - - Interpolating the omega_gw extrema could be challenging, particularly - near the merger due to the rapid change in the omega_gw values and - limited number of data points. Thus, it requires special treatment. For - other data, this problem is less likely to occur. Therefore, we use - this function to interpolate data other than omega_gw extrema. The - spline kwargs for this function is provided by `extra_kwargs` using the - key "non_extrema_interp_kwargs". - - A wrapper of utils.get_interpolant with check_kwargs=False. - This is to make sure that the checking of kwargs is not performed - everytime the interpolation function is called. Instead, the kwargs - are checked once in the init and passed to the interpolation - function without repeating checks. - """ - return get_interpolant(oldX, oldY, allowExtrapolation, interpolator, - spline_kwargs=self.non_extrema_interp_kwargs, - check_kwargs=False) - - def get_spline_interpolant_for_extrema( - self, oldX, oldY, allowExtrapolation=False, interpolator="spline"): - """Get spline interpolant for omega_gw extrema. - - Same as `get_spline_interpolant_for_non_extrema` but uses the - `spline_kwargs` provided via `self.extrema_interp_kwargs`. For more, - see `get_spline_interpolant_for_non_extrema`. - """ - return get_interpolant(oldX, oldY, allowExtrapolation, interpolator, - spline_kwargs=self.extrema_interp_kwargs, - check_kwargs=False) - - def spline_interpolate_non_extrema(self, newX, oldX, oldY, allowExtrapolation=False, - interpolator="spline"): - """Get interpolated values using spline for data other than omega_gw - extrema. - - A wrapper of utils.interpolate with check_kwargs=False for reasons - explained in the documentation of - `get_spline_interpolant_for_non_extrema` function. - """ - return interpolate(newX, oldX, oldY, allowExtrapolation, interpolator, - spline_kwargs=self.non_extrema_interp_kwargs, - check_kwargs=False) - - def get_rational_fit_wrapper(self, x, y): - """Get rational fit. - - A wrapper of `utils.get_rational_fit` with check_kwargs=False. This is - to make sure that the checking of kwargs is not performed everytime the - rational fit function is called. Instead, the kwargs are checked once - in the init and passed to the rational fit function without repeating - checks. - """ - return get_rational_fit(x, y, - rational_fit_kwargs=self.extrema_interp_kwargs, - check_kwargs=False) - def get_rational_fit_for_extrema(self, x, y, data_name=None): """Get rational fit with adaptive numerator and denominator degree. @@ -1513,12 +1454,12 @@ def get_rational_fit_for_extrema(self, x, y, data_name=None): the previous valid monotonic fit. The initial degrees for the rational fit can be specified through - `extrema_interp_kwargs` in `extra_kwargs` with the key "rational_fit". - (See `extrema_interp_kwargs` under `extra_kwargs` for more details) - Default values are provided by `get_default_rational_fit_kwargs`, where - both degrees are set to `None`. If both degrees remain `None`, - appropriate starting values are determined using - `self.get_approximate_degree_for_rational_fit`. + `special_interp_kwargs_for_extrema` in `extra_kwargs` with the key + "rational_fit". (See `special_interp_kwargs_for_extrema` under + `extra_kwargs` for more details) Default values are provided by + `get_default_rational_fit_kwargs`, where both degrees are set to + `None`. If both degrees remain `None`, appropriate starting values are + determined using `self.get_approximate_degree_for_rational_fit`. """ # make sure that data_name is not None. A data_name is needed to # update the optimal values of the numerator and denominator degrees @@ -1544,32 +1485,37 @@ def get_rational_fit_for_extrema(self, x, y, data_name=None): # as the ideal degree can change based on the eccentricity which we do # not know apriori. The true optimal values are determined through # further iterations. - if self.extrema_interp_kwargs["num_degree"] is None \ - and self.extrema_interp_kwargs["denom_degree"] is None: - self.extrema_interp_kwargs["num_degree"], self.extrema_interp_kwargs["denom_degree"] = \ + if self.special_interp_kwargs_for_extrema["num_degree"] is None \ + and self.special_interp_kwargs_for_extrema["denom_degree"] is None: + self.special_interp_kwargs_for_extrema["num_degree"], \ + self.special_interp_kwargs_for_extrema["denom_degree"] = \ self.get_approximate_degree_for_rational_fit() - elif self.extrema_interp_kwargs["num_degree"] is None: - self.extrema_interp_kwargs["num_degree"] = self.extrema_interp_kwargs["denom_degree"] - elif self.extrema_interp_kwargs["denom_degree"] is None: - self.extrema_interp_kwargs["denom_degree"] = self.extrema_interp_kwargs["num_degree"] - - rat_fit = self.get_rational_fit_wrapper(x, y) + elif self.special_interp_kwargs_for_extrema["num_degree"] is None: + self.special_interp_kwargs_for_extrema["num_degree"] \ + = self.special_interp_kwargs_for_extrema["denom_degree"] + elif self.special_interp_kwargs_for_extrema["denom_degree"] is None: + self.special_interp_kwargs_for_extrema["denom_degree"] \ + = self.special_interp_kwargs_for_extrema["num_degree"] + + rat_fit = get_rational_fit( + x, y, rational_fit_kwargs=self.special_interp_kwargs_for_extrema, + check_kwargs=False) x_test = np.arange(x[0], x[-1], self.t[1] - self.t[0]) # save the degrees for checks at each step of iterations for finding # the optimal degrees - old_num_degree = self.extrema_interp_kwargs["num_degree"] - old_denom_degree = self.extrema_interp_kwargs["denom_degree"] + old_num_degree = self.special_interp_kwargs_for_extrema["num_degree"] + old_denom_degree = self.special_interp_kwargs_for_extrema["denom_degree"] # Check for nonmonotonicity and lower degrees if needed fit_is_nonmonotonic = self.check_if_first_derivative_is_not_strictly_monotonic( x_test, rat_fit(x_test), 1.0, data_name) while fit_is_nonmonotonic: - if self.extrema_interp_kwargs["num_degree"] > 1: - self.extrema_interp_kwargs["num_degree"] -= 1 - if self.extrema_interp_kwargs["denom_degree"] > 1: - self.extrema_interp_kwargs["denom_degree"] -= 1 - if self.extrema_interp_kwargs["num_degree"] == 1 \ - and self.extrema_interp_kwargs["denom_degree"] == 1: + if self.special_interp_kwargs_for_extrema["num_degree"] > 1: + self.special_interp_kwargs_for_extrema["num_degree"] -= 1 + if self.special_interp_kwargs_for_extrema["denom_degree"] > 1: + self.special_interp_kwargs_for_extrema["denom_degree"] -= 1 + if self.special_interp_kwargs_for_extrema["num_degree"] == 1 \ + and self.special_interp_kwargs_for_extrema["denom_degree"] == 1: debug_message( "Both numerator and denominator degrees are equal to 1 " "and cannot be lowered further.", @@ -1578,13 +1524,19 @@ def get_rational_fit_for_extrema(self, x, y, data_name=None): debug_message( "Lowering degrees to " - f"num_degree={self.extrema_interp_kwargs['num_degree']}, " - f"denom_degree={self.extrema_interp_kwargs['denom_degree']} " + "num_degree=" + f"{self.special_interp_kwargs_for_extrema['num_degree']}, " + "denom_degree=" + f"{self.special_interp_kwargs_for_extrema['denom_degree']} " "and retrying.", debug_level=self.debug_level, important=False) # build new fit and check monotonicity - rat_fit = self.get_rational_fit_wrapper(x, y) - fit_is_nonmonotonic = self.check_if_first_derivative_is_not_strictly_monotonic( + rat_fit = get_rational_fit( + x, y, + rational_fit_kwargs=self.special_interp_kwargs_for_extrema, + check_kwargs=False) + fit_is_nonmonotonic \ + = self.check_if_first_derivative_is_not_strictly_monotonic( x_test, rat_fit(x_test), 1.0, data_name) # If fit with initial degree is monotonic, try increasing the degree # for a better fit @@ -1595,12 +1547,18 @@ def get_rational_fit_for_extrema(self, x, y, data_name=None): while not fit_is_nonmonotonic: # Increase the degrees for both numerator and denominator - new_num_degree = self.extrema_interp_kwargs["num_degree"] + 1 - new_denom_degree = self.extrema_interp_kwargs["denom_degree"] + 1 - self.extrema_interp_kwargs["num_degree"], self.extrema_interp_kwargs["denom_degree"] \ + new_num_degree \ + = self.special_interp_kwargs_for_extrema["num_degree"] + 1 + new_denom_degree \ + = self.special_interp_kwargs_for_extrema["denom_degree"] + 1 + self.special_interp_kwargs_for_extrema["num_degree"],\ + self.special_interp_kwargs_for_extrema["denom_degree"] \ = new_num_degree, new_denom_degree # build new fit and check monotonicity - new_rat_fit = self.get_rational_fit_wrapper(x, y) + new_rat_fit = get_rational_fit( + x, y, + rational_fit_kwargs=self.special_interp_kwargs_for_extrema, + check_kwargs=False) fit_is_nonmonotonic \ = self.check_if_first_derivative_is_not_strictly_monotonic( x_test, new_rat_fit(x_test), 1.0, data_name) @@ -1608,13 +1566,14 @@ def get_rational_fit_for_extrema(self, x, y, data_name=None): # Revert to previous fit and degrees if nonmonotonicity is # detected debug_message( - "Increasing degrees caused nonmonotonicity. Reverting to " - f"last monotonic fit with num_degree={last_monotonic_num_degree} " - f"and denom_degree={last_monotonic_denom_degree}.", + "Increasing degrees caused nonmonotonicity. " + "Reverting to last monotonic fit with " + f"num_degree={last_monotonic_num_degree} and " + f"denom_degree={last_monotonic_denom_degree}.", debug_level=self.debug_level, important=False) - self.extrema_interp_kwargs["num_degree"] \ + self.special_interp_kwargs_for_extrema["num_degree"] \ = last_monotonic_num_degree - self.extrema_interp_kwargs["denom_degree"] \ + self.special_interp_kwargs_for_extrema["denom_degree"] \ = last_monotonic_denom_degree rat_fit = last_monotonic_rat_fit break @@ -1626,8 +1585,8 @@ def get_rational_fit_for_extrema(self, x, y, data_name=None): # update final degrees used to build the fit self.rational_fit_degrees[data_name] = ( - self.extrema_interp_kwargs["num_degree"], - self.extrema_interp_kwargs["denom_degree"]) + self.special_interp_kwargs_for_extrema["num_degree"], + self.special_interp_kwargs_for_extrema["denom_degree"]) return rat_fit def get_approximate_degree_for_rational_fit(self): @@ -1670,15 +1629,15 @@ def check_if_first_derivative_is_not_strictly_monotonic( == "rational_fit": if data_name in self.rational_fit_nonmonotonicity_history: self.rational_fit_nonmonotonicity_history[data_name].update( - {(self.extrema_interp_kwargs['num_degree'], - self.extrema_interp_kwargs['num_degree']): + {(self.special_interp_kwargs_for_extrema['num_degree'], + self.special_interp_kwargs_for_extrema['num_degree']): is_not_strictly_monotonic}) return is_not_strictly_monotonic def get_available_omega_gw_extrema_interpolation_methods(self): """Return available omega_gw extrema interpolation methods.""" available_methods = { - "spline": self.get_spline_interpolant_for_extrema, + "spline": get_interpolant, "rational_fit": self.get_rational_fit_for_extrema, } return available_methods @@ -1718,12 +1677,13 @@ def get_omega_gw_extrema_interpolant(self, extrema_type="pericenters"): method = self.available_omega_gw_extrema_interpolation_methods[ self.omega_gw_extrema_interpolation_method] if self.omega_gw_extrema_interpolation_method == "rational_fit": - # rational fit method `get_rational_fit_for_extrema` takes an - # additional argument return method(self.t[extrema], self.omega_gw[extrema], extrema_type) - else: - return method(self.t[extrema], self.omega_gw[extrema]) + if self.omega_gw_extrema_interpolation_method == "spline": + return method( + self.t[extrema], self.omega_gw[extrema], + spline_kwargs=self.special_interp_kwargs_for_extrema, + check_kwargs=False) else: raise Exception( f"Sufficient number of {extrema_type} are not found." @@ -2026,11 +1986,11 @@ def measure_ecc(self, tref_in=None, fref_in=None): = self.check_num_extrema(apocenters, "apocenters") # If the eccentricity is too small for a method to find the extrema, - # and `return_zero_if_small_ecc_failure` is true, then we set the eccentricity and - # mean anomaly to zero and return them. In this case, the rest of the - # code in this function is not executed, and therefore, many variables - # that are needed for making diagnostic plots are not computed. Thus, - # in such cases, the diagnostic plots may not work. + # and `return_zero_if_small_ecc_failure` is true, then we set the + # eccentricity and mean anomaly to zero and return them. In this case, + # the rest of the code in this function is not executed, and therefore, + # many variables that are needed for making diagnostic plots are not + # computed. Thus, in such cases, the diagnostic plots may not work. if any([insufficient_pericenters_but_long_waveform, insufficient_apocenters_but_long_waveform]) \ and self.return_zero_if_small_ecc_failure: @@ -2118,8 +2078,10 @@ def measure_ecc(self, tref_in=None, fref_in=None): 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") + self.omega_gw_pericenters_interp \ + = self.get_omega_gw_extrema_interpolant("pericenters") + self.omega_gw_apocenters_interp \ + = self.get_omega_gw_extrema_interpolant("apocenters") # check monotonicity of the interpolants self.check_omega_gw_extrema_interpolants() # compute eccentricity at self.tref_out @@ -2150,13 +2112,16 @@ def check_omega_gw_extrema_interpolants(self): """ # Verify the monotonicity of the first derivative of the omega_gw # interpolant with spline. - if self.extra_kwargs["omega_gw_extrema_interpolation_method"] == "spline": + if self.extra_kwargs["omega_gw_extrema_interpolation_method"] \ + == "spline": # Check if the first derivative of omega_gw at pericenters or # apocenters is non-monotonic if (self.check_if_first_derivative_is_not_strictly_monotonic( - self.t_for_checks, self.omega_gw_pericenters_interp(self.t_for_checks)) or + self.t_for_checks, + self.omega_gw_pericenters_interp(self.t_for_checks)) or self.check_if_first_derivative_is_not_strictly_monotonic( - self.t_for_checks, self.omega_gw_apocenters_interp(self.t_for_checks))): + self.t_for_checks, + self.omega_gw_apocenters_interp(self.t_for_checks))): debug_message( "Nonmonotonic time derivative detected in the spline " "interpolant through extrema. Using rational fit by " @@ -2295,8 +2260,9 @@ def derivative_of_eccentricity(self, t, n=1): self.t_for_checks) if self.ecc_interp is None: - self.ecc_interp = self.get_spline_interpolant_for_non_extrema( - self.t_for_checks, self.ecc_for_checks) + self.ecc_interp = get_interpolant( + self.t_for_checks, self.ecc_for_checks, + spline_kwargs=self.general_interp_kwargs, check_kwargs=False) # Get derivative of ecc(t) using spline return self.ecc_interp.derivative(n=n)(t) @@ -2619,13 +2585,15 @@ def compute_res_amp_gw_and_res_omega_gw(self): # residual quantities can be computed. Above, we check that this # extrapolation does not happen before t_merger, which is where # eccentricity is normally measured. - self.amp_gw_zeroecc_interp = self.spline_interpolate_non_extrema( + self.amp_gw_zeroecc_interp = interpolate( self.t, self.t_zeroecc_shifted, self.amp_gw_zeroecc, - allowExtrapolation=True) + allowExtrapolation=True, spline_kwargs=self.general_interp_kwargs, + check_kwargs=False) self.res_amp_gw = self.amp_gw - self.amp_gw_zeroecc_interp - self.omega_gw_zeroecc_interp = self.spline_interpolate_non_extrema( + self.omega_gw_zeroecc_interp = interpolate( self.t, self.t_zeroecc_shifted, self.omega_gw_zeroecc, - allowExtrapolation=True) + allowExtrapolation=True, spline_kwargs=self.general_interp_kwargs, + check_kwargs=False) self.res_omega_gw = (self.omega_gw - self.omega_gw_zeroecc_interp) def get_t_average_for_orbit_averaged_omega_gw(self): @@ -2756,8 +2724,9 @@ def compute_orbit_averaged_omega_gw_between_extrema(self, t): orbit_averaged_omega_gw, "omega_gw averaged [apocenter to apocenter] and " "[pericenter to pericenter]") - return self.spline_interpolate_non_extrema( - t, self.t_for_orbit_averaged_omega_gw, orbit_averaged_omega_gw) + return interpolate( + t, self.t_for_orbit_averaged_omega_gw, orbit_averaged_omega_gw, + spline_kwargs=self.general_interp_kwargs, check_kwargs=False) def check_monotonicity_of_omega_gw_average(self, omega_gw_average, @@ -2860,8 +2829,9 @@ def compute_mean_of_extrema_interpolants(self, t): def compute_omega_gw_zeroecc(self, t): """Find omega_gw from zeroecc data.""" - return self.spline_interpolate_non_extrema( - t, self.t_zeroecc_shifted, self.omega_gw_zeroecc) + return interpolate( + t, self.t_zeroecc_shifted, self.omega_gw_zeroecc, + spline_kwargs=self.general_interp_kwargs, check_kwargs=False) def get_available_omega_gw_averaging_methods(self): """Return available omega_gw averaging methods.""" @@ -3005,16 +2975,19 @@ def compute_tref_in_and_fref_out_from_fref_in(self, fref_in): # of omega_gw_average. # We get omega_gw_average by evaluating the omega_gw_average(t) # on t, from tmin_for_fref to tmax_for_fref - self.t_for_omega_gw_average, self.omega_gw_average = self.get_omega_gw_average(method) + self.t_for_omega_gw_average, self.omega_gw_average \ + = self.get_omega_gw_average(method) # check that omega_gw_average is monotonically increasing self.check_monotonicity_of_omega_gw_average( self.omega_gw_average, "Interpolated omega_gw_average") # Get tref_in using interpolation - tref_in = self.spline_interpolate_non_extrema(fref_out, + tref_in = interpolate(fref_out, self.omega_gw_average/(2 * np.pi), - self.t_for_omega_gw_average) + self.t_for_omega_gw_average, + spline_kwargs=self.general_interp_kwargs, + check_kwargs=False) # check if tref_in is monotonically increasing if any(np.diff(tref_in) <= 0): debug_message(f"tref_in from fref_in using method {method} is" @@ -3053,7 +3026,8 @@ def get_fref_bounds(self, method=None): frequency. """ if self.omega_gw_average is None: - self.t_for_omega_gw_average, self.omega_gw_average = self.get_omega_gw_average(method) + self.t_for_omega_gw_average, self.omega_gw_average \ + = self.get_omega_gw_average(method) return [min(self.omega_gw_average)/2/np.pi, max(self.omega_gw_average)/2/np.pi] diff --git a/gw_eccentricity/gw_eccentricity.py b/gw_eccentricity/gw_eccentricity.py index 8f83ea4e..79ca22c7 100644 --- a/gw_eccentricity/gw_eccentricity.py +++ b/gw_eccentricity/gw_eccentricity.py @@ -338,7 +338,7 @@ def measure_eccentricity(tref_in=None, Default value is "inertial". extra_kwargs: A dict of any extra kwargs to be passed. Allowed kwargs are: - extrema_interp_kwargs: dict + special_interp_kwargs_for_extrema: dict A dictionary with a single key matching the current `omega_gw_extrema_interpolation_method`. See under `omega_gw_extrema_interpolation_method` for more details on @@ -356,7 +356,7 @@ def measure_eccentricity(tref_in=None, - "rational_fit": default kwargs are set using `utils.get_default_rational_fit_kwargs` - non_extrema_interp_kwargs: dict + general_interp_kwargs: dict Dictionary of arguments to be passed to the spline interpolation routine (scipy.interpolate.InterpolatedUnivariateSpline) used to interpolate data other than the omega_gw extrema.