From 9cdd8c0b500aea6a5e0894f35bc216e80734da26 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Mon, 2 Jun 2025 21:52:09 +0100 Subject: [PATCH 01/48] Enable bilby-pipe style injections --- dingo/core/result.py | 18 +-- dingo/core/samplers.py | 4 +- dingo/gw/noise/asd_dataset.py | 13 ++ dingo/pipe/dag_creator.py | 30 ++-- dingo/pipe/data_generation.py | 101 +++++++----- dingo/pipe/importance_sampling.py | 47 +++++- dingo/pipe/main.py | 136 +++++++++++++--- dingo/pipe/parser.py | 257 +++++++++++++++++------------- 8 files changed, 399 insertions(+), 207 deletions(-) diff --git a/dingo/core/result.py b/dingo/core/result.py index fa7bedda5..613c25444 100644 --- a/dingo/core/result.py +++ b/dingo/core/result.py @@ -87,8 +87,8 @@ def base_metadata(self): @property def injection_parameters(self): - if self.context: - return self.context.get("parameters") + if self.event_metadata: + return self.event_metadata.get("injection_parameters") else: return None @@ -596,7 +596,6 @@ def plot_corner( self, parameters: list = None, filename: str = "corner.pdf", - truths: dict = None, **kwargs, ): """ @@ -609,8 +608,6 @@ def plot_corner( (Default: None) filename : str Where to save samples. - truths : dict - Dictionary of truth values to include. Other Parameters ---------------- @@ -631,8 +628,8 @@ def plot_corner( # User option to plot specific parameters. if parameters: theta = theta[parameters] - if truths is not None: - kwargs["truths"] = [truths.get(k) for k in theta.columns] + if self.injection_parameters is not None: + kwargs["truths"] = [self.injection_parameters.get(k) for k in theta.columns] if weights is not None: plot_corner_multi( @@ -643,12 +640,7 @@ def plot_corner( **kwargs, ) else: - plot_corner_multi( - theta, - labels=["Dingo"], - filename=filename, - **kwargs - ) + plot_corner_multi(theta, labels=["Dingo"], filename=filename, **kwargs) def plot_log_probs(self, filename="log_probs.png"): """ diff --git a/dingo/core/samplers.py b/dingo/core/samplers.py index 377882750..f5dab4649 100644 --- a/dingo/core/samplers.py +++ b/dingo/core/samplers.py @@ -106,7 +106,9 @@ def context(self): @context.setter def context(self, value): if value is not None and "parameters" in value: - self.metadata["injection_parameters"] = value.pop("parameters") + if self.event_metadata is None: + self.event_metadata = {} + self.event_metadata["injection_parameters"] = value.pop("parameters") self._context = value @property diff --git a/dingo/gw/noise/asd_dataset.py b/dingo/gw/noise/asd_dataset.py index 4c260b3d8..ae7b6c40f 100644 --- a/dingo/gw/noise/asd_dataset.py +++ b/dingo/gw/noise/asd_dataset.py @@ -1,4 +1,5 @@ import copy +from pathlib import Path from typing import Iterable, Optional from dingo.gw.domains import build_domain, UniformFrequencyDomain @@ -193,6 +194,18 @@ def sample_random_asds(self, n: Optional[int] = None) -> dict[str, np.ndarray]: else: return {k: v[np.random.choice(len(v), n)] for k, v in self.asds.items()} + def save_psd(self, directory, ifo_name, idx: Optional[int] = None): + if idx is None: + idx = np.random.choice(len(self.asds[ifo_name])) + directory = Path(directory) + directory.mkdir(exist_ok=True) + psd_path = directory / f"{ifo_name}_{idx}_psd.txt" + np.savetxt( + psd_path, + np.vstack([self.domain(), self.asds[ifo_name][idx] ** 2]).T, + ) + return psd_path + def check_domain_compatibility(data: dict, domain: BaseFrequencyDomain) -> bool: for v in data.values(): diff --git a/dingo/pipe/dag_creator.py b/dingo/pipe/dag_creator.py index 10006d642..210231897 100644 --- a/dingo/pipe/dag_creator.py +++ b/dingo/pipe/dag_creator.py @@ -20,13 +20,13 @@ def get_trigger_time_list(inputs): """Returns a list of GPS trigger times for each data segment""" - # if (inputs.gaussian_noise or inputs.zero_noise) and inputs.trigger_time is None: - # trigger_times = [0] * inputs.n_simulation - # elif (inputs.gaussian_noise or inputs.zero_noise) and isinstance( - # inputs.trigger_time, float - # ): - # trigger_times = [inputs.trigger_time] * inputs.n_simulation - if inputs.trigger_time is not None: + if (inputs.gaussian_noise or inputs.zero_noise) and inputs.trigger_time is None: + trigger_times = [0] * inputs.n_simulation + elif (inputs.gaussian_noise or inputs.zero_noise) and isinstance( + inputs.trigger_time, float + ): + trigger_times = [inputs.trigger_time] * inputs.n_simulation + elif inputs.trigger_time is not None: trigger_times = [inputs.trigger_time] elif getattr(inputs, "gpstimes", None) is not None: start_times = inputs.gpstimes @@ -44,7 +44,7 @@ def get_parallel_list(inputs): return [f"part{idx}" for idx in range(inputs.n_parallel)] -def generate_dag(inputs, model_args): +def generate_dag(inputs): inputs = copy.deepcopy(inputs) dag = Dag(inputs) trigger_times = get_trigger_time_list(inputs) @@ -86,7 +86,9 @@ def generate_dag(inputs, model_args): # # 3. Generate new data for importance sampling **if different settings requested**. # - # If injecting into simulated noise, be sure to use consistent noise realization. + # If injecting into simulated noise, be sure to use consistent noise + # realization. Currently simulated noise + importance_sampling_updates is + # prohibited in MainInput. if len(inputs.importance_sampling_updates) > 0: # Iterate over all generation nodes and store them in a list @@ -98,7 +100,9 @@ def generate_dag(inputs, model_args): # Ensures any cached files (e.g. the distance-marginalization # lookup table) are only built once. kwargs["parent"] = generation_node_list[0] - generation_node = GenerationNode(inputs, importance_sampling=True, **kwargs) + generation_node = GenerationNode( + inputs, importance_sampling=True, **kwargs + ) importance_sampling_generation_node_list.append(generation_node) else: importance_sampling_generation_node_list = generation_node_list @@ -158,9 +162,9 @@ def generate_dag(inputs, model_args): # if inputs.create_summary: - # Add the waveform approximant to inputs, so that it can be fed to PESummary. - inputs.waveform_approximant = model_args["waveform_approximant"] - PESummaryNode(inputs, merged_importance_sampling_node_list, generation_node_list, dag=dag) + PESummaryNode( + inputs, merged_importance_sampling_node_list, generation_node_list, dag=dag + ) dag.build() # create_overview( diff --git a/dingo/pipe/data_generation.py b/dingo/pipe/data_generation.py index b0af17925..32dee08f0 100644 --- a/dingo/pipe/data_generation.py +++ b/dingo/pipe/data_generation.py @@ -9,6 +9,7 @@ from dingo.gw.data.event_dataset import EventDataset from dingo.gw.domains import UniformFrequencyDomain +from dingo.gw.gwutils import get_window, get_window_factor from dingo.pipe.parser import create_parser logger.name = "dingo_pipe" @@ -33,7 +34,7 @@ def __init__(self, args, unknown_args, create_data=True): # Run index arguments self.idx = args.idx - # self.generation_seed = args.generation_seed + self.generation_seed = args.generation_seed self.trigger_time = args.trigger_time # Naming arguments @@ -46,7 +47,7 @@ def __init__(self, args, unknown_args, create_data=True): # self.phase_marginalization = args.phase_marginalization # self.prior_file = args.prior_file # self.prior_dict = args.prior_dict - # self.deltaT = args.deltaT + self.deltaT = args.deltaT # self.default_prior = args.default_prior # Whether to generate data for importance sampling. This must be done when @@ -66,7 +67,8 @@ def __init__(self, args, unknown_args, create_data=True): self.data_format = args.data_format self.allow_tape = args.allow_tape self.tukey_roll_off = args.tukey_roll_off - self.zero_noise = False # dingo mod + self.gaussian_noise = args.gaussian_noise + self.zero_noise = args.zero_noise self.resampling_method = args.resampling_method if args.timeslide_dict is not None: @@ -85,35 +87,37 @@ def __init__(self, args, unknown_args, create_data=True): self.sampling_frequency = args.sampling_frequency self.minimum_frequency = args.minimum_frequency self.maximum_frequency = args.maximum_frequency - # self.reference_frequency = args.reference_frequency + self.reference_frequency = args.reference_frequency # Waveform, source model and likelihood - # self.waveform_generator_class = args.waveform_generator - # self.waveform_approximant = args.waveform_approximant - # self.catch_waveform_errors = args.catch_waveform_errors - # self.pn_spin_order = args.pn_spin_order - # self.pn_tidal_order = args.pn_tidal_order - # self.pn_phase_order = args.pn_phase_order - # self.pn_amplitude_order = args.pn_amplitude_order - # self.mode_array = args.mode_array - # self.waveform_arguments_dict = args.waveform_arguments_dict - # self.numerical_relativity_file = args.numerical_relativity_file - # self.injection_waveform_approximant = args.injection_waveform_approximant - # self.frequency_domain_source_model = args.frequency_domain_source_model + self.waveform_generator_class = ( + "bilby.gw.waveform_generator.LALCBCWaveformGenerator" + ) + self.waveform_approximant = args.waveform_approximant + self.catch_waveform_errors = args.catch_waveform_errors + # TODO: These are set to parser defaults. Fix to set from model. + self.pn_spin_order = -1 + self.pn_tidal_order = -1 + self.pn_phase_order = -1 + self.pn_amplitude_order = 0 + self.mode_array = None + self.waveform_arguments_dict = None + self.numerical_relativity_file = args.numerical_relativity_file + self.injection_waveform_approximant = args.injection_waveform_approximant + self.frequency_domain_source_model = "lal_binary_black_hole" # self.conversion_function = args.conversion_function # self.generation_function = args.generation_function # self.likelihood_type = args.likelihood_type # self.extra_likelihood_kwargs = args.extra_likelihood_kwargs - # self.enforce_signal_duration = args.enforce_signal_duration + self.enforce_signal_duration = args.enforce_signal_duration # PSD self.psd_maximum_duration = args.psd_maximum_duration self.psd_dict = args.psd_dict - if self.psd_dict is None: - self.psd_length = args.psd_length - self.psd_fractional_overlap = args.psd_fractional_overlap - self.psd_start_time = args.psd_start_time - self.psd_method = args.psd_method + self.psd_length = args.psd_length + self.psd_fractional_overlap = args.psd_fractional_overlap + self.psd_start_time = args.psd_start_time + self.psd_method = args.psd_method # # ROQ # self.roq_folder = args.roq_folder @@ -147,19 +151,9 @@ def __init__(self, args, unknown_args, create_data=True): # # Plotting self.plot_data = args.plot_data self.plot_spectrogram = args.plot_spectrogram - # self.plot_injection = args.plot_injection + self.plot_injection = args.plot_injection if create_data: - # Added for dingo so that create_data runs. TODO: enable injections - args.injection = False - args.injection_numbers = None - args.injection_file = None - args.injection_dict = None - args.injection_waveform_arguments = None - args.injection_frequency_domain_source_model = None - self.frequency_domain_source_model = None - self.gaussian_noise = False - self.create_data(args) def save_hdf5(self): @@ -173,7 +167,6 @@ def save_hdf5(self): # PSD and strain data. data = {"waveform": {}, "asds": {}} # TODO: Rename these keys. for ifo in self.interferometers: - strain = ifo.strain_data.frequency_domain_strain frequency_array = ifo.strain_data.frequency_array asd = ifo.power_spectral_density.get_amplitude_spectral_density_array( @@ -224,24 +217,42 @@ def save_hdf5(self): "psd_method", "channel_dict", "data_dict", + "injection", + "generation_seed", + "gaussian_noise", + "zero_noise", + "numerical_relativity_file", + "injection_waveform_approximant", + "injection_frequency_domain_source_model", ]: try: v = getattr(self, k) except AttributeError: continue - if v is not None: + if v is not None and v is not False: settings[k] = v + if self.injection: + settings["injection_parameters"] = self.injection_parameters.copy() + # Dingo and Bilby have different geocent_time conventions. + settings["injection_parameters"]["geocent_time"] -= self.trigger_time + settings["optimal_SNR"] = { + k: v["optimal_SNR"] for k, v in self.interferometers.meta_data.items() + } + settings["matched_filter_SNR"] = { + k: v["matched_filter_SNR"] + for k, v in self.interferometers.meta_data.items() + } + dataset = EventDataset( dictionary={ "data": data, - # "event_metadata": event_metadata, "settings": settings, } ) dataset.to_file(self.event_data_file) - # also saving the psd as a .dat file which can be read in + # also saving the psd as a .txt file which can be read in # easily by pesummary or bilby for ifo in self.interferometers: np.savetxt( @@ -268,6 +279,22 @@ def importance_sampling_updates(self, setting): else: self._importance_sampling_updates = None + def _set_interferometers_from_gaussian_noise(self): + super()._set_interferometers_from_gaussian_noise() + # Scale the FD strain appropriately by the window factor (not done by default + # in Bilby / bilby_pipe). This is to ensure that the injection is consistent + # with TD data with a given PSD, making it consistent also with DINGO network + # training. + for ifo in self.interferometers: + # This is a hack to set the window factor. It ensures also that the SNRs + # are calculated correctly. + td_strain = ifo.time_domain_strain + ifo.strain_data.time_domain_window(roll_off=self.tukey_roll_off) + ifo.strain_data.frequency_domain_strain = ( + ifo.strain_data.frequency_domain_strain + * np.sqrt(ifo.strain_data.window_factor) + ) + def create_generation_parser(): """Data generation parser creation""" diff --git a/dingo/pipe/importance_sampling.py b/dingo/pipe/importance_sampling.py index 850bf2a4e..ceb6fe06e 100644 --- a/dingo/pipe/importance_sampling.py +++ b/dingo/pipe/importance_sampling.py @@ -6,7 +6,13 @@ import yaml from bilby_pipe.input import Input -from bilby_pipe.utils import parse_args, logger, convert_string_to_dict +from bilby_pipe.utils import ( + parse_args, + logger, + convert_string_to_dict, + convert_prior_string_input, + BilbyPipeError, +) from dingo.gw.data.event_dataset import EventDataset from dingo.gw.domains import MultibandedFrequencyDomain @@ -44,6 +50,7 @@ def __init__(self, args, unknown_args): self.prior_dict = args.prior_dict self.default_prior = "PriorDict" self.time_reference = "geocent" + self.prior_dict_updates = args.prior_dict_updates # Choices for running # self.detectors = args.detectors @@ -158,6 +165,8 @@ def importance_sampling_settings(self, settings): self._importance_sampling_settings.update( convert_string_to_dict(settings) ) + if "phase_marginalization" in self._importance_sampling_settings: + self._importance_sampling_settings.pop("synthetic_phase", None) else: self._importance_sampling_settings = dict() @@ -166,16 +175,16 @@ def run_sampler(self): "use_base_domain", False ) - if self.prior_dict: + if self.prior_dict_updates: logger.info("Updating prior from network prior. Changes:") logger.info( yaml.dump( - self.prior_dict, + self.prior_dict_updates, default_flow_style=False, sort_keys=False, ) ) - self.result.update_prior(self.prior_dict) + self.result.update_prior(self.prior_dict_updates) if "synthetic_phase" in self.importance_sampling_settings: logger.info("Sampling synthetic phase.") @@ -196,7 +205,6 @@ def run_sampler(self): calibration_marginalization_kwargs=self.calibration_marginalization_kwargs, ) - self.result.print_summary() self.result.to_file(os.path.join(self.result_directory, self.label + ".hdf5")) @@ -207,6 +215,35 @@ def priors(self): self._priors = self._get_priors(add_time=False) return self._priors + @property + def prior_dict_updates(self): + """The input prior_dict from the ini (if given) + + Note, this is not the bilby prior (see self.priors for that), this is + a key-val dictionary where the val's are strings which are converting + into bilby priors in `_get_prior + """ + return self._prior_dict_updates + + @prior_dict_updates.setter + def prior_dict_updates(self, prior_dict_updates): + if isinstance(prior_dict_updates, dict): + prior_dict_updates = prior_dict_updates + elif isinstance(prior_dict_updates, str): + prior_dict_updates = convert_prior_string_input(prior_dict_updates) + elif prior_dict_updates is None: + self._prior_dict_updates = None + return + else: + raise BilbyPipeError( + f"prior_dict_updates={prior_dict_updates} not " f"understood" + ) + + self._prior_dict_updates = { + self._convert_prior_dict_key(key): val + for key, val in prior_dict_updates.items() + } + def create_sampling_parser(): """Data analysis parser creation""" diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py index 955d51b3b..3921a7a73 100644 --- a/dingo/pipe/main.py +++ b/dingo/pipe/main.py @@ -1,8 +1,10 @@ # # Adapted from bilby_pipe. In particular, uses the bilby_pipe data generation code. # +import copy import os +from bilby.core.prior import PriorDict from bilby_pipe.input import Input from bilby_pipe.main import MainInput as BilbyMainInput from bilby_pipe.utils import ( @@ -10,6 +12,8 @@ get_command_line_arguments, logger, parse_args, + convert_prior_string_input, + BilbyPipeError, ) from dingo.core.posterior_models.build_model import build_model_from_kwargs @@ -20,11 +24,20 @@ from ..gw.domains.build_domain import build_domain_from_model_metadata from dingo.core.posterior_models.build_model import build_model_from_kwargs +from ..gw.injection import Injection +from ..gw.noise.asd_dataset import ASDDataset logger.name = "dingo_pipe" def fill_in_arguments_from_model(args): + if args.prior_dict is not None: + raise ValueError( + "Do not specify prior-dict in INI file. This is obtained from " + "the DINGO model. To update the prior, specify " + "prior-dict-updates." + ) + logger.info(f"Loading dingo model from {args.model} in order to access settings.") try: @@ -41,6 +54,15 @@ def fill_in_arguments_from_model(args): domain = build_domain_from_model_metadata(model_metadata, base=True) + prior = Injection.from_posterior_model_metadata(model_metadata).prior + deltaT = prior["geocent_time"].maximum - prior["geocent_time"].minimum + + # Dingo and Bilby have different conventions for the geocent_time prior. Dingo + # always centers it around 0.0, whereas Bilby centers it around the trigger time. + # Drop the geocent_time parameter here so that bilby_pipe re-builds it as needed + # for Bilby. + del prior["geocent_time"] + data_settings = model_metadata["train_settings"]["data"] model_args = { @@ -53,6 +75,11 @@ def fill_in_arguments_from_model(args): "waveform_approximant": model_metadata["dataset_settings"][ "waveform_generator" ]["approximant"], + "reference_frequency": model_metadata["dataset_settings"]["waveform_generator"][ + "f_ref" + ], + "deltaT": deltaT, + "prior_dict": prior, } changed_args = {} @@ -85,6 +112,21 @@ def fill_in_arguments_from_model(args): # TODO: Also check consistency between model and init_model settings. + # If an ASDDataset is specified, pull out PSDs and save them as .txt files. Use + # these for downstream tasks by specifying the psd_dict. + if args.asd_dataset: + asd_dataset = ASDDataset(file_name=args.asd_dataset) + domain_dict = copy.deepcopy(domain.domain_dict) + if "window_factor" in domain_dict: + del domain_dict["window_factor"] + asd_dataset.update_domain(domain_dict) + psd_dict = {} + for ifo_name in args.detectors: + psd_path = asd_dataset.save_psd(args.outdir, ifo_name) + psd_dict[ifo_name] = str(psd_path) + args.asd_dataset = None + args.psd_dict = str(psd_dict) + # Updates that are explicitly provided take priority. if args.importance_sampling_updates is None: importance_sampling_updates = {} @@ -99,13 +141,16 @@ def fill_in_arguments_from_model(args): class MainInput(BilbyMainInput): - def __init__(self, args, unknown_args, importance_sampling_updates): + def __init__( + self, args, unknown_args, importance_sampling_updates, perform_checks=True + ): # Settings added for dingo. self.model = args.model self.model_init = args.model_init self.num_gnpe_iterations = args.num_gnpe_iterations self.importance_sampling_updates = importance_sampling_updates + self.prior_dict_updates = args.prior_dict_updates Input.__init__(self, args, unknown_args, print_msg=False) @@ -135,10 +180,12 @@ def __init__(self, args, unknown_args, importance_sampling_updates): self.data_find_urltype = args.data_find_urltype self.n_parallel = args.n_parallel # useful when condor nodes don't have access to submit filesystem - self.transfer_files = args.transfer_files + self.transfer_files = args.transfer_files self.additional_transfer_paths = args.additional_transfer_paths self.osg = args.osg - self.desired_sites = args.cpu_desired_sites # Dummy variable so bilby_pipe doesn't complain. + self.desired_sites = ( + args.cpu_desired_sites + ) # Dummy variable so bilby_pipe doesn't complain. self.cpu_desired_sites = args.cpu_desired_sites self.gpu_desired_sites = args.gpu_desired_sites # self.analysis_executable = args.analysis_executable @@ -162,7 +209,7 @@ def __init__(self, args, unknown_args, importance_sampling_updates): self.environment_variables = args.environment_variables self.getenv = args.getenv - # self.waveform_approximant = args.waveform_approximant + self.waveform_approximant = args.waveform_approximant # # self.time_reference = args.time_reference self.time_reference = "geocent" @@ -170,11 +217,11 @@ def __init__(self, args, unknown_args, importance_sampling_updates): # self.likelihood_type = args.likelihood_type self.duration = args.duration # self.phase_marginalization = args.phase_marginalization - self.prior_file = None # Dingo update. To change prior use the priod_dict. + self.prior_file = None # Dingo update. To change prior use the prior_dict. self.prior_dict = args.prior_dict - self.default_prior = "PriorDict" + self.default_prior = "BBHPriorDict" self.minimum_frequency = args.minimum_frequency - # self.enforce_signal_duration = args.enforce_signal_duration + self.enforce_signal_duration = args.enforce_signal_duration self.run_local = args.local self.generation_pool = args.generation_pool @@ -185,23 +232,30 @@ def __init__(self, args, unknown_args, importance_sampling_updates): # self.ignore_gwpy_data_quality_check = args.ignore_gwpy_data_quality_check self.trigger_time = args.trigger_time - # self.deltaT = args.deltaT + self.deltaT = args.deltaT self.gps_tuple = args.gps_tuple self.gps_file = args.gps_file self.timeslide_file = args.timeslide_file - self.gaussian_noise = False # DINGO MOD: Cannot use different noise types. - self.zero_noise = False # DINGO MOD: does not support zero noise yet - # self.n_simulation = args.n_simulation - # - # self.injection = args.injection - # self.injection_numbers = args.injection_numbers + self.gaussian_noise = args.gaussian_noise + self.zero_noise = args.zero_noise + self.n_simulation = args.n_simulation + + self.injection = args.injection + self.injection_numbers = args.injection_numbers self.injection_file = args.injection_file - # self.injection_dict = args.injection_dict + self.injection_dict = args.injection_dict # self.injection_waveform_arguments = args.injection_waveform_arguments # self.injection_waveform_approximant = args.injection_waveform_approximant - # self.generation_seed = args.generation_seed - # if self.injection: - # self.check_injection() + # self.injection_frequency_domain_source_model = ( + # args.injection_frequency_domain_source_model + # ) + self.generation_seed = args.generation_seed + if self.importance_sampling_updates and self.gaussian_noise: + raise ValueError( + "Cannot update data for importance sampling if using " + "simulated Gaussian noise. This risks inconsistent noise " + "realizations." + ) self.importance_sample = args.importance_sample @@ -251,7 +305,12 @@ def __init__(self, args, unknown_args, importance_sampling_updates): self.psd_start_time = args.psd_start_time self.spline_calibration_envelope_dict = args.spline_calibration_envelope_dict - # self.check_source_model(args) + if perform_checks: + # self.check_source_model(args) + # self.check_calibration_prior_boundary(args) + # self.check_cpu_parallelisation() + if self.injection: + self.check_injection() self.requirements = [] self.device = args.device @@ -280,11 +339,38 @@ def request_cpus_importance_sampling(self, request_cpus_importance_sampling): self._request_cpus_importance_sampling = request_cpus_importance_sampling @property - def priors(self): - """Read in and compose the prior at run-time""" - if getattr(self, "_priors", None) is None: - self._priors = self._get_priors(add_time=False) - return self._priors + def prior_dict_updates(self): + """The input prior_dict from the ini (if given) + + Note, this is not the bilby prior (see self.priors for that), this is + a key-val dictionary where the val's are strings which are converting + into bilby priors in `_get_prior + """ + return self._prior_dict_updates + + @prior_dict_updates.setter + def prior_dict_updates(self, prior_dict_updates): + if isinstance(prior_dict_updates, dict): + prior_dict_updates = prior_dict_updates + elif isinstance(prior_dict_updates, str): + prior_dict_updates = convert_prior_string_input(prior_dict_updates) + elif prior_dict_updates is None: + self._prior_dict_updates = None + return + else: + raise BilbyPipeError( + f"prior_dict_updates={prior_dict_updates} not " f"understood" + ) + + self._prior_dict_updates = { + self._convert_prior_dict_key(key): val + for key, val in prior_dict_updates.items() + } + + def _get_priors(self, add_time=True): + priors = super()._get_priors(add_time=add_time) + priors.update(PriorDict(self.prior_dict_updates)) + return priors def write_complete_config_file(parser, args, inputs, input_cls=MainInput): @@ -357,7 +443,7 @@ def main(): # TODO: Use two sets of inputs! The first must match the network; the second is # used in importance sampling. - generate_dag(inputs, model_args) + generate_dag(inputs) if len(unknown_args) > 0: print(f"Unrecognized arguments {unknown_args}") diff --git a/dingo/pipe/parser.py b/dingo/pipe/parser.py index 554358f8b..cd17353b5 100644 --- a/dingo/pipe/parser.py +++ b/dingo/pipe/parser.py @@ -82,9 +82,11 @@ def create_parser(top_level=True): "--calibration-correction-type", type=nonestr, default="data", - help=("Type of calibration correction: can be either `data` or `template`." - " See https://bilby-dev.github.io/bilby/api/bilby.gw.detector.calibration.html " - "for more information.") + help=( + "Type of calibration correction: can be either `data` or `template`." + " See https://bilby-dev.github.io/bilby/api/bilby.gw.detector.calibration.html " + "for more information." + ), ) calibration_parser.add( @@ -233,15 +235,15 @@ def create_parser(top_level=True): "trigger time" ), ) - # data_gen_pars.add( - # "--n-simulation", - # type=int, - # default=0, - # help=( - # "Number of simulated segments to use with gaussian-noise " - # "Note, this must match the number of injections specified" - # ), - # ) + data_gen_pars.add( + "--n-simulation", + type=int, + default=0, + help=( + "Number of simulated segments to use with gaussian-noise " + "Note, this must match the number of injections specified" + ), + ) data_gen_pars.add( "--data-dict", default=None, @@ -298,17 +300,17 @@ def create_parser(top_level=True): default="osdf", help="URL type to use for datafind, default is osdf", ) - # data_type_pars = data_gen_pars.add_mutually_exclusive_group() - # data_type_pars.add( - # "--gaussian-noise", - # action="store_true", - # help="If true, use simulated Gaussian noise", - # ) - # data_type_pars.add( - # "--zero-noise", - # action="store_true", - # help="Use a zero noise realisation", - # ) + data_type_pars = data_gen_pars.add_mutually_exclusive_group() + data_type_pars.add( + "--gaussian-noise", + action="store_true", + help="If true, use simulated Gaussian noise", + ) + data_type_pars.add( + "--zero-noise", + action="store_true", + help="Use a zero noise realisation", + ) det_parser = parser.add_argument_group( title="Detector arguments", @@ -338,20 +340,28 @@ def create_parser(top_level=True): default=None, help="The duration of data around the event to use", ) - # det_parser.add( - # "--generation-seed", - # default=None, - # type=noneint, - # help=( - # "Random seed used during data generation. If no generation seed " - # "provided, a random seed between 1 and 1e6 is selected. If a seed " - # "is provided, it is used as the base seed and all generation jobs " - # "will have their seeds set as {generation_seed = base_seed + job_idx}." - # ), - # ) det_parser.add( + "--generation-seed", + default=None, + type=noneint, + help=( + "Random seed used during data generation. If no generation seed " + "provided, a random seed between 1 and 1e6 is selected. If a seed " + "is provided, it is used as the base seed and all generation jobs " + "will have their seeds set as {generation_seed = base_seed + job_idx}." + ), + ) + psd_dict_parser = det_parser.add_mutually_exclusive_group() + psd_dict_parser.add( "--psd-dict", type=nonestr, default=None, help="Dictionary of PSD files to use" ) + psd_dict_parser.add( + "--asd-dataset", + type=nonestr, + default=None, + help="DINGO ASDDataset file to be used for injections. If specified, dingo_pipe " + "will generate PSD files based on random ASDs in the dataset.", + ) det_parser.add( "--psd-fractional-overlap", # default=0.5, @@ -445,12 +455,12 @@ def create_parser(top_level=True): title="Injection arguments", description="Whether to include software injections and how to generate them.", ) - # injection_parser.add( - # "--injection", - # action="store_true", - # default=False, - # help="Create data from an injection file", - # ) + injection_parser.add( + "--injection", + action="store_true", + default=False, + help="Create data from an injection file", + ) injection_parser_input = injection_parser.add_mutually_exclusive_group() injection_parser_input.add( "--injection-dict", @@ -467,36 +477,44 @@ def create_parser(top_level=True): " for supported formats" ), ) - # injection_parser.add( - # "--injection-numbers", - # action="append", - # type=nonestr, - # default=None, - # help=( - # "Specific injections rows to use from the injection_file, e.g. " - # "`injection_numbers=[0,3] selects the zeroth and third row. Can be " - # "a list of slice-syntax values, e.g, [0, 2:4] will produce [0, 2, 3]. " - # "Repeated entries will be ignored." - # ), - # ) - # injection_parser.add( - # "--injection-waveform-approximant", - # type=nonestr, - # default=None, - # help="The name of the waveform approximant to use to create injections. " - # "If none is specified, then the `waveform-approximant` will be used" - # "as the `injection-waveform-approximant`.", - # ) - # injection_parser.add( - # "--injection-waveform-arguments", - # type=nonestr, - # default=None, - # help=( - # "A dictionary of arbitrary additional waveform-arguments to pass " - # "to the bilby waveform generator's waveform arguments for the " - # "injection only" - # ), - # ) + injection_parser.add( + "--injection-numbers", + action="append", + type=nonestr, + default=None, + help=( + "Specific injections rows to use from the injection_file, e.g. " + "`injection_numbers=[0,3] selects the zeroth and third row. Can be " + "a list of slice-syntax values, e.g, [0, 2:4] will produce [0, 2, 3]. " + "Repeated entries will be ignored." + ), + ) + injection_parser.add( + "--injection-waveform-approximant", + type=nonestr, + default=None, + help="The name of the waveform approximant to use to create injections. " + "If none is specified, then the `waveform-approximant` will be used" + "as the `injection-waveform-approximant`.", + ) + injection_parser.add( + "--injection-frequency-domain-source-model", + type=nonestr, + default=None, + help="Frequency domain source model to use for generating injections. " + "If this is None, it will default to the frequency domain source model" + "used for analysis.", + ) + injection_parser.add( + "--injection-waveform-arguments", + type=nonestr, + default=None, + help=( + "A dictionary of arbitrary additional waveform-arguments to pass " + "to the bilby waveform generator's waveform arguments for the " + "injection only" + ), + ) submission_parser = parser.add_argument_group( title="Job submission arguments", @@ -940,11 +958,11 @@ def create_parser(top_level=True): action="store_true", help="Create plot of the frequency domain data", ) - # output_parser.add_argument( - # "--plot-injection", - # action="store_true", - # help="Create time-domain plot of the injection", - # ) + output_parser.add_argument( + "--plot-injection", + action="store_true", + help="Create time-domain plot of the injection", + ) output_parser.add_argument( "--plot-spectrogram", action="store_true", @@ -1067,7 +1085,7 @@ def create_parser(top_level=True): ) # prior_parser.add( # "--default-prior", - # default="PriorDict", + # default="BBHPriorDict", # type=str, # help=( # "The name of the prior set to base the prior on. Can be one of" @@ -1075,15 +1093,30 @@ def create_parser(top_level=True): # "or a python path to a bilby prior class available in the user's installation." # ), # ) - # prior_parser.add( - # "--deltaT", - # type=float, - # default=0.2, - # help=( - # "The symmetric width (in s) around the trigger time to" - # " search over the coalescence time" - # ), - # ) + prior_parser.add( + "--deltaT", + type=float, + default=0.2, + help=( + "The symmetric width (in s) around the trigger time to" + " search over the coalescence time" + ), + ) + prior_parser.add( + "--prior-dict-updates", + type=nonestr, + default=None, + help=( + "A dictionary of priors. Multiline " + "dictionaries are supported, but each line must contain a single" + "parameter specification and finish with a comma. Dingo priors are set at " + "network training time, so the prior-dict-update is used at importance " + "sampling " + "time to re-weight to the new prior. The prior-dict does not have to be " + "a prior over the entire set of parameters, only the parameters for which " + "the prior is changed." + ), + ) prior_parser_main = prior_parser.add_mutually_exclusive_group() # prior_parser_main.add( # "--prior-file", type=nonestr, default=None, help="The prior file" @@ -1096,21 +1129,19 @@ def create_parser(top_level=True): "A dictionary of priors (alternative to prior-file). Multiline " "dictionaries are supported, but each line must contain a single" "parameter specification and finish with a comma. Dingo priors are set at " - "network training time, so the prior-dict is used at importance sampling " - "time to re-weight to the new prior. The prior-dict does not have to be " - "a prior over the entire set of parameters, only the parameters for which " - "the prior is changed." + "network training time, so the prior-dict should not be provided by the " + "user! Use 'prior-dict-update' to update the prior for importance sampling." + ), + ) + prior_parser.add( + "--enforce-signal-duration", + action=StoreBoolean, + default=True, + help=( + "Whether to require that all signals fit within the segment duration. " + "The signal duration is calculated using a post-Newtonian approximation." ), ) - # prior_parser.add( - # "--enforce-signal-duration", - # action=StoreBoolean, - # default=True, - # help=( - # "Whether to require that all signals fit within the segment duration. " - # "The signal duration is calculated using a post-Newtonian approximation." - # ), - # ) # postprocessing_parser = parser.add_argument_group( # title="Post processing arguments", @@ -1191,12 +1222,12 @@ def create_parser(top_level=True): type=nonestr, help="The name of the waveform approximant to use for PE.", ) - # waveform_parser.add( - # "--catch-waveform-errors", - # default=True, - # action=StoreBoolean, - # help="Turns on waveform error catching", - # ) + waveform_parser.add( + "--catch-waveform-errors", + default=True, + action=StoreBoolean, + help="Turns on waveform error catching", + ) # waveform_parser.add( # "--pn-spin-order", # default=-1, @@ -1222,15 +1253,15 @@ def create_parser(top_level=True): # help="Post-Newtonian order to use for the amplitude. Also " # "used to determine the waveform starting frequency.", # ) - # waveform_parser.add( - # "--numerical-relativity-file", - # default=None, - # type=nonestr, - # help=( - # "Path to a h5 numerical relativity file to inject, see" - # "https://git.ligo.org/waveforms/lvcnr-lfs for examples" - # ), - # ) + waveform_parser.add( + "--numerical-relativity-file", + default=None, + type=nonestr, + help=( + "Path to a h5 numerical relativity file to inject, see" + "https://git.ligo.org/waveforms/lvcnr-lfs for examples" + ), + ) # waveform_parser.add( # "--waveform-arguments-dict", # default=None, From 0136e143a56e2900aedf99d40a8845bab4ef0aa5 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Tue, 17 Jun 2025 16:18:05 +0000 Subject: [PATCH 02/48] allow complex data in hdf5 --- dingo/core/dataset.py | 2 +- dingo/pipe/parser.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/dingo/core/dataset.py b/dingo/core/dataset.py index 9361f0360..4b64148ec 100644 --- a/dingo/core/dataset.py +++ b/dingo/core/dataset.py @@ -18,7 +18,7 @@ def recursive_hdf5_save(group, d): group.create_dataset(k, data=v) elif isinstance(v, pd.DataFrame): group.create_dataset(k, data=v.to_records(index=False)) - elif isinstance(v, (int, float, str, list)): + elif isinstance(v, (int, float, complex, str, list)): # TODO: Set scalars as attributes? group.create_dataset(k, data=v) else: diff --git a/dingo/pipe/parser.py b/dingo/pipe/parser.py index cd17353b5..80e1f7cc4 100644 --- a/dingo/pipe/parser.py +++ b/dingo/pipe/parser.py @@ -1095,8 +1095,8 @@ def create_parser(top_level=True): # ) prior_parser.add( "--deltaT", - type=float, - default=0.2, + type=nonefloat, + default=None, help=( "The symmetric width (in s) around the trigger time to" " search over the coalescence time" From 5829f02e9c0906de56859de1a46f25ce6618794f Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Wed, 18 Jun 2025 15:56:15 +0100 Subject: [PATCH 03/48] Enable DeltaFunction time priors. Monkey-hack of bilby_pipe.create_injections.py. --- dingo/pipe/create_injections.py | 67 ++++++++++++++++++++++++ dingo/pipe/main.py | 92 +++++++++++++++++++++++++++++++-- dingo/pipe/parser.py | 9 ++++ 3 files changed, 164 insertions(+), 4 deletions(-) create mode 100644 dingo/pipe/create_injections.py diff --git a/dingo/pipe/create_injections.py b/dingo/pipe/create_injections.py new file mode 100644 index 000000000..5165b4b7c --- /dev/null +++ b/dingo/pipe/create_injections.py @@ -0,0 +1,67 @@ +import functools + +import numpy as np + +import bilby +import bilby_pipe.input as input_mod +import bilby_pipe.create_injections as ci + +_bilby_pipe_create_injection_file = ci.create_injection_file + + +@functools.wraps(_bilby_pipe_create_injection_file) +def create_injection_file(*args, **kwargs): + # Monkey-patch create_injection_file() to allow for an offset in the time. This is + # needed if the DINGO network time prior is not centered at 0.0, e.g., if we have a + # Dirac delta prior. + + Toffset = kwargs.pop("Toffset", 0.0) + + if kwargs.get("trigger_time") is not None: + kwargs["trigger_time"] += Toffset + if kwargs.get("gpstimes") is not None: + gpstimes = kwargs["gpstimes"] + kwargs["gpstimes"] = ( + (np.asarray(gpstimes) + Toffset).tolist() + if isinstance(gpstimes, list) + else gpstimes + Toffset + ) + + return _bilby_pipe_create_injection_file(*args, **kwargs) + + +create_injection_file.__doc__ = ( + "Add a `Toffset` keyword to offset the trigger and GPS times before creating injections.\n\n" + + _bilby_pipe_create_injection_file.__doc__ +) + +ci.create_injection_file = create_injection_file + +_bilby_pipe_get_time_prior = input_mod.get_time_prior + + +@functools.wraps(_bilby_pipe_get_time_prior) +def get_time_prior(time, uncertainty, name="geocent_time", latex_label="$t_c$"): + # Monkey-patch to allow for DeltaFunction priors. + if uncertainty == 0.0: + return bilby.core.prior.DeltaFunction( + peak=time, + name=name, + latex_label=latex_label, + unit="$s$", + ) + elif uncertainty > 0.0: + return _bilby_pipe_get_time_prior( + time, uncertainty, name=name, latex_label=latex_label + ) + else: + raise ValueError(f"Time uncertainty {uncertainty} < 0.0.") + + +get_time_prior.__doc__ = ( + "A bilby.core.prior.DeltaFunction for the time parameter, if uncertainty == 0.0.\n\n" + + _bilby_pipe_get_time_prior.__doc__ +) + + +input_mod.get_time_prior = get_time_prior diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py index 3921a7a73..118ffafed 100644 --- a/dingo/pipe/main.py +++ b/dingo/pipe/main.py @@ -2,8 +2,11 @@ # Adapted from bilby_pipe. In particular, uses the bilby_pipe data generation code. # import copy +import json import os +import dingo.pipe.create_injections + from bilby.core.prior import PriorDict from bilby_pipe.input import Input from bilby_pipe.main import MainInput as BilbyMainInput @@ -55,12 +58,16 @@ def fill_in_arguments_from_model(args): domain = build_domain_from_model_metadata(model_metadata, base=True) prior = Injection.from_posterior_model_metadata(model_metadata).prior - deltaT = prior["geocent_time"].maximum - prior["geocent_time"].minimum + geocent_time_prior = prior["geocent_time"] + deltaT = geocent_time_prior.maximum - geocent_time_prior.minimum + + # Offset needed in case of DeltaFunction priors not peaked at 0.0. + Toffset = 0.5 * (geocent_time_prior.minimum + geocent_time_prior.maximum) # Dingo and Bilby have different conventions for the geocent_time prior. Dingo - # always centers it around 0.0, whereas Bilby centers it around the trigger time. - # Drop the geocent_time parameter here so that bilby_pipe re-builds it as needed - # for Bilby. + # always centers it around something close to 0.0, whereas Bilby centers it around + # the trigger time. Drop the geocent_time parameter here so that bilby_pipe + # re-builds it as needed for Bilby. del prior["geocent_time"] data_settings = model_metadata["train_settings"]["data"] @@ -79,6 +86,7 @@ def fill_in_arguments_from_model(args): "f_ref" ], "deltaT": deltaT, + "Toffset": Toffset, "prior_dict": prior, } @@ -233,6 +241,7 @@ def __init__( # self.ignore_gwpy_data_quality_check = args.ignore_gwpy_data_quality_check self.trigger_time = args.trigger_time self.deltaT = args.deltaT + self.Toffset = args.Toffset self.gps_tuple = args.gps_tuple self.gps_file = args.gps_file self.timeslide_file = args.timeslide_file @@ -372,6 +381,81 @@ def _get_priors(self, add_time=True): priors.update(PriorDict(self.prior_dict_updates)) return priors + def check_injection(self): + """Check injection behaviour + + If injections are requested, either use the injection-dict, + injection-file, or create an injection-file + + """ + default_injection_file_name = "{}/{}_injection_file.dat".format( + self.data_directory, self.label + ) + if self.injection_dict is not None: + logger.debug( + "Using injection dict from ini file {}".format( + json.dumps(self.injection_dict, indent=2) + ) + ) + elif self.injection_file is not None: + logger.debug(f"Using injection file {self.injection_file}") + elif os.path.isfile(default_injection_file_name): + # This is done to avoid overwriting the injection file + logger.debug(f"Using injection file {default_injection_file_name}") + self.injection_file = default_injection_file_name + else: + logger.debug("No injection file found, generating one now") + + if self.gps_file is not None or self.gps_tuple is not None: + if self.n_simulation > 0 and self.n_simulation != len(self.gpstimes): + raise BilbyPipeError( + "gps_file/gps_tuple option and n_simulation are not matched" + ) + gpstimes = self.gpstimes + n_injection = len(gpstimes) + else: + gpstimes = None + n_injection = self.n_simulation + + if self.trigger_time is None: + trigger_time_injections = 0 + else: + trigger_time_injections = self.trigger_time + + dingo.pipe.create_injections.create_injection_file( + filename=default_injection_file_name, + prior_file=self.prior_file, + prior_dict=self.prior_dict, + n_injection=n_injection, + trigger_time=trigger_time_injections, + deltaT=self.deltaT, + Toffset=self.Toffset, # dingo-pipe mod + gpstimes=gpstimes, + duration=self.duration, + post_trigger_duration=self.post_trigger_duration, + generation_seed=self.generation_seed, + extension="dat", + default_prior=self.default_prior, + ) + self.injection_file = default_injection_file_name + + # Check the gps_file has the sample length as number of simulation + if self.gps_file is not None: + if len(self.gpstimes) != len(self.injection_df): + raise BilbyPipeError("Injection file length does not match gps_file") + + if self.n_simulation > 0: + if self.n_simulation != len(self.injection_df): + raise BilbyPipeError( + "n-simulation does not match the number of injections: " + "please check your ini file" + ) + elif self.n_simulation == 0 and self.gps_file is None: + self.n_simulation = len(self.injection_df) + logger.debug( + f"Setting n_simulation={self.n_simulation} to match injections" + ) + def write_complete_config_file(parser, args, inputs, input_cls=MainInput): args_dict = vars(args).copy() diff --git a/dingo/pipe/parser.py b/dingo/pipe/parser.py index 80e1f7cc4..221f93329 100644 --- a/dingo/pipe/parser.py +++ b/dingo/pipe/parser.py @@ -1102,6 +1102,15 @@ def create_parser(top_level=True): " search over the coalescence time" ), ) + prior_parser.add( + "--Toffset", + type=nonefloat, + default=None, + help=( + "Offset of the center of the time prior from the trigger time (in s). " + "Useful for DeltaFunction priors. Generally set by DINGO model." + ), + ) prior_parser.add( "--prior-dict-updates", type=nonestr, From e3b019d24e9b96101a108cef4802f52d80db467d Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Thu, 19 Jun 2025 17:09:24 +0100 Subject: [PATCH 04/48] Update docs and examples. --- docs/source/dingo_pipe.md | 4 ++-- examples/gnpe_model/GW150914.ini | 2 +- examples/npe_model/GW150914.ini | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/source/dingo_pipe.md b/docs/source/dingo_pipe.md index abb628f47..d9ebf3ad4 100644 --- a/docs/source/dingo_pipe.md +++ b/docs/source/dingo_pipe.md @@ -34,7 +34,7 @@ num-samples = 50000 batch-size = 50000 recover-log-prob = true importance-sample = true -prior-dict = { +prior-dict-updates = { luminosity_distance = bilby.gw.prior.UniformComovingVolume(minimum=100, maximum=2000, name='luminosity_distance'), } @@ -97,7 +97,7 @@ Since sampling uses GPU hardware, there is an additional key `sampling-requireme For importance sampling, the Result saved in the previous step is loaded. Since this contains the strain data and ASDs, as well as all settings used for training the network, the likelihood and prior can be evaluated for each sample point. If it is necessary to change data conditioning or PSD for importance sampling (i.e., if the `importance-sampling-updates` dictionary is non-empty), then a second [data generation](#data-generation) step is first carried using the new settings, and used as importance sampling context. The importance sampled result is finally saved as HDF5, including the estimated Bayesian evidence. -If a `prior-dict` is specified in the `.ini` file, then this will be used for the importance sampling prior. One example where this is useful is for the luminosity distance prior. Indeed, Dingo tends to train better using a uniform prior over luminosity distance, but physically one would prefer a uniform in volume prior. By specifying a `prior-dict` this change can be made in importance sampling. +If `prior-dict-updates` is specified in the `.ini` file, then this will be used for the importance sampling prior. One example where this is useful is for the luminosity distance prior. Indeed, Dingo tends to train better using a uniform prior over luminosity distance, but physically one would prefer a uniform in volume prior. By specifying `prior-dict-updates` this change can be made in importance sampling. ```{caution} If extending the prior support during importance sampling, be sure that the posterior does not rail up against the prior boundary being extended. diff --git a/examples/gnpe_model/GW150914.ini b/examples/gnpe_model/GW150914.ini index 0bb52091d..5f3f0acae 100644 --- a/examples/gnpe_model/GW150914.ini +++ b/examples/gnpe_model/GW150914.ini @@ -19,7 +19,7 @@ num-gnpe-iterations = 30 num-samples = 50000 batch-size = 50000 recover-log-prob = true -prior-dict = { +prior-dict-updates = { luminosity_distance = bilby.gw.prior.UniformComovingVolume(minimum=100, maximum=2000, name='luminosity_distance'), } diff --git a/examples/npe_model/GW150914.ini b/examples/npe_model/GW150914.ini index 66b8a4756..4ddfacf13 100644 --- a/examples/npe_model/GW150914.ini +++ b/examples/npe_model/GW150914.ini @@ -18,7 +18,7 @@ device = 'cuda' num-samples = 50000 batch-size = 50000 recover-log-prob = true -prior-dict = { +prior-dict-updates = { luminosity_distance = bilby.gw.prior.UniformComovingVolume(minimum=100, maximum=2000, name='luminosity_distance'), } From 643f315b2da9662f5b635ac66f752dd4ff529a6a Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Thu, 19 Jun 2025 22:11:55 +0100 Subject: [PATCH 05/48] Ensure prior-dict-updates included when sampling injection parameters. --- dingo/pipe/main.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py index 118ffafed..6ebbce03a 100644 --- a/dingo/pipe/main.py +++ b/dingo/pipe/main.py @@ -425,7 +425,9 @@ def check_injection(self): dingo.pipe.create_injections.create_injection_file( filename=default_injection_file_name, prior_file=self.prior_file, - prior_dict=self.prior_dict, + prior_dict=self._get_priors( + add_time=False + ), # ensure prior-dict-updates incorporated n_injection=n_injection, trigger_time=trigger_time_injections, deltaT=self.deltaT, From 4e33207a80434c97df910a480868c0d6ca8fc132 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Fri, 20 Jun 2025 13:23:56 +0100 Subject: [PATCH 06/48] Fix issue with Sampler.log_prob. Ensure everything but the inference parameters are stripped from theta. For samples coming from Result, existing values for log_prob were being added to new ones. --- dingo/core/samplers.py | 13 +++++--- dingo/gw/inference/gw_samplers.py | 51 +++++++++++++++---------------- 2 files changed, 34 insertions(+), 30 deletions(-) diff --git a/dingo/core/samplers.py b/dingo/core/samplers.py index f5dab4649..638904e44 100644 --- a/dingo/core/samplers.py +++ b/dingo/core/samplers.py @@ -229,23 +229,28 @@ def run_sampler( print(f"Done. This took {time.time() - t0:.1f} s.") sys.stdout.flush() - def log_prob(self, samples: pd.DataFrame) -> np.ndarray: + def log_prob(self, samples: pd.DataFrame | dict) -> np.ndarray: """ Calculate the model log probability at specific sample points. Parameters ---------- - samples : pd.DataFrame + samples : pd.DataFrame | dict Sample points at which to calculate the log probability. Returns ------- np.array of log probabilities. """ - # TODO: Check / fix this method. It is likely broken, but is not critical. if self.context is None and not self.unconditional_model: raise ValueError("Context must be set in order to calculate log_prob.") + if isinstance(samples, dict): + if np.isscalar(list(samples.values())[0]): + samples = pd.DataFrame([samples]) + else: + samples = pd.DataFrame(samples) + # This undoes any post-correction that would have been done to the samples, # before evaluating the log_prob. E.g., the t_ref / sky position correction. samples = samples.copy() @@ -267,7 +272,7 @@ def log_prob(self, samples: pd.DataFrame) -> np.ndarray: # Context is the same for each sample. Expand across batch dimension after # pre-processing. x = self.transform_pre(self.context) - x = x.expand(len(samples), *x.shape) + x = x.expand(len(samples), *x.shape) # TODO: Make this more efficient. x = [x] else: x = [] diff --git a/dingo/gw/inference/gw_samplers.py b/dingo/gw/inference/gw_samplers.py index 4e0b39966..263666d6a 100644 --- a/dingo/gw/inference/gw_samplers.py +++ b/dingo/gw/inference/gw_samplers.py @@ -4,7 +4,7 @@ import numpy as np import pandas as pd from astropy.time import Time -from bilby.core.prior import PriorDict, DeltaFunction +from bilby.core.prior import PriorDict, DeltaFunction, Constraint from bilby.gw.detector import InterferometerList from torchvision.transforms import Compose @@ -87,11 +87,9 @@ def _correct_reference_time( Whether to apply instead the inverse transformation. This is used prior to calculating the log_prob. """ - # FIXME: Update this method to make sure that for models that do not include sky - # position, it does not do anything. if self.event_metadata is not None: t_event = self.event_metadata.get("time_event") - if t_event is not None and t_event != self.t_ref: + if t_event is not None and t_event != self.t_ref and "ra" in samples: ra = samples["ra"] time_reference = Time(self.t_ref, format="gps", scale="utc") time_event = Time(t_event, format="gps", scale="utc") @@ -108,10 +106,10 @@ def _correct_reference_time( def _post_process(self, samples: Union[dict, pd.DataFrame], inverse: bool = False): """ - Post processing of parameter samples. + Post-processing of parameter samples. + * Add any fixed parameters from the prior. * Correct the sky position for a potentially fixed reference time. (see self._correct_reference_time) - * Potentially sample a synthetic phase. (see self._sample_synthetic_phase) This method modifies the samples in place. @@ -122,34 +120,35 @@ def _post_process(self, samples: Union[dict, pd.DataFrame], inverse: bool = Fals Whether to apply instead the inverse transformation. This is used prior to calculating the log_prob. """ - # Add fixed parameters from prior intrinsic_prior = self.metadata["dataset_settings"]["intrinsic_prior"] extrinsic_prior = get_extrinsic_prior_dict( self.metadata["train_settings"]["data"]["extrinsic_prior"] ) prior = build_prior_with_defaults({**intrinsic_prior, **extrinsic_prior}) - num_samples = len(samples[list(samples.keys())[0]]) - for k, p in prior.items(): - if isinstance(p, DeltaFunction) and k not in samples: - v = p.peak - print(f"Adding fixed parameter {k} = {v} from prior.") - samples[k] = p.peak * np.ones(num_samples) + + if not inverse: + # Add fixed parameters from prior. + num_samples = len(samples[list(samples.keys())[0]]) + for k, p in prior.items(): + if isinstance(p, DeltaFunction) and k not in samples: + v = p.peak + print(f"Adding fixed parameter {k} = {v} from prior.") + samples[k] = p.peak * np.ones(num_samples) + else: + # Drop non-inference parameters from samples. + # NOTE: Important to drop "log_prob" in particular before running + # Sampler.log_prob(), otherwise log probabilities are added. + drop_parameters = [ + k for k in samples.keys() if k not in self.inference_parameters + ] + if isinstance(samples, pd.DataFrame): + samples.drop(columns=drop_parameters, inplace=True, errors="ignore") + elif isinstance(samples, dict): + for k in drop_parameters: + samples.pop(k, None) if not self.unconditional_model: self._correct_reference_time(samples, inverse) - # if not inverse: - # self._correct_reference_time(samples, inverse) - # # if self.synthetic_phase_kwargs is not None: - # # print(f"Sampling synthetic phase.") - # # t0 = time.time() - # # self._sample_synthetic_phase(samples, inverse) - # # print(f"Done. This took {time.time() - t0:.2f} seconds.") - # - # # If inverting, we go in reverse order. - # else: - # # if self.synthetic_phase_kwargs is not None: - # # self._sample_synthetic_phase(samples, inverse) - # self._correct_reference_time(samples, inverse) class GWSampler(GWSamplerMixin, Sampler): From f676c5defedb1d1a3c92cca48ad7542e1f8961e8 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Mon, 23 Jun 2025 11:33:50 +0100 Subject: [PATCH 07/48] Set default trigger time to the model reference time for injections. This is needed for models with DeltaFunction ra priors. --- dingo/pipe/dag_creator.py | 2 +- dingo/pipe/main.py | 12 +++++++++++- dingo/pipe/parser.py | 7 +++++++ 3 files changed, 19 insertions(+), 2 deletions(-) diff --git a/dingo/pipe/dag_creator.py b/dingo/pipe/dag_creator.py index 210231897..8ed8c8eac 100644 --- a/dingo/pipe/dag_creator.py +++ b/dingo/pipe/dag_creator.py @@ -21,7 +21,7 @@ def get_trigger_time_list(inputs): """Returns a list of GPS trigger times for each data segment""" if (inputs.gaussian_noise or inputs.zero_noise) and inputs.trigger_time is None: - trigger_times = [0] * inputs.n_simulation + trigger_times = [inputs.model_reference_time] * inputs.n_simulation elif (inputs.gaussian_noise or inputs.zero_noise) and isinstance( inputs.trigger_time, float ): diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py index 6ebbce03a..7e8d5a135 100644 --- a/dingo/pipe/main.py +++ b/dingo/pipe/main.py @@ -40,6 +40,11 @@ def fill_in_arguments_from_model(args): "the DINGO model. To update the prior, specify " "prior-dict-updates." ) + if args.model_reference_time is not None: + raise ValueError( + "Do not specify model-reference-time in INI file. This is obtained from the " + "DINGO model." + ) logger.info(f"Loading dingo model from {args.model} in order to access settings.") @@ -88,6 +93,7 @@ def fill_in_arguments_from_model(args): "deltaT": deltaT, "Toffset": Toffset, "prior_dict": prior, + "model_reference_time": model_metadata["train_settings"]["data"]["ref_time"], } changed_args = {} @@ -159,6 +165,7 @@ def __init__( self.num_gnpe_iterations = args.num_gnpe_iterations self.importance_sampling_updates = importance_sampling_updates self.prior_dict_updates = args.prior_dict_updates + self.model_reference_time = args.model_reference_time Input.__init__(self, args, unknown_args, print_msg=False) @@ -418,7 +425,10 @@ def check_injection(self): n_injection = self.n_simulation if self.trigger_time is None: - trigger_time_injections = 0 + # Use the model reference time as the trigger time by default. This + # avoids needing to adjust sky position in post-postprocessing, + # and it avoids complications due to DeltaFunction ra priors. + trigger_time_injections = self.model_reference_time else: trigger_time_injections = self.trigger_time diff --git a/dingo/pipe/parser.py b/dingo/pipe/parser.py index 221f93329..aa54d47f9 100644 --- a/dingo/pipe/parser.py +++ b/dingo/pipe/parser.py @@ -1349,6 +1349,13 @@ def create_parser(top_level=True): help="Neural network model for generating samples to initialize Gibbs sampling." "Must be provided if the main model is a GNPE model. ", ) + sampler_parser.add( + "--model-reference-time", + type=nonefloat, + default=None, + help="Reference time for neural network. Do not add this manually as it is " + "specified by the network.", + ) sampler_parser.add( "--recover-log-prob", action=StoreBoolean, From 794d441ffcdcd957dbc54c0d6795f3bbd85bfa9e Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Mon, 23 Jun 2025 17:32:29 +0100 Subject: [PATCH 08/48] Add pp-plot function to core.Result. Adapted from Bilby. Also added parameter names to default priors, so that they pick up latex labels for plotting. --- dingo/core/result.py | 223 +++++++++++++++++++++++++++++++++++++++++++ dingo/gw/prior.py | 41 ++++---- dingo/pipe/main.py | 3 - 3 files changed, 242 insertions(+), 25 deletions(-) diff --git a/dingo/core/result.py b/dingo/core/result.py index 613c25444..4821bcb15 100644 --- a/dingo/core/result.py +++ b/dingo/core/result.py @@ -2,11 +2,14 @@ import math import tempfile import time +from collections import namedtuple +from itertools import product import numpy as np from typing import Optional import pandas as pd +import scipy from matplotlib import pyplot as plt from scipy.constants import golden from scipy.special import logsumexp @@ -705,6 +708,226 @@ def plot_weights(self, filename="weights.png"): else: print("Results not importance sampled. Cannot plot weights.") + def get_all_injection_credible_levels( + self, keys: list[str] = None, weighted: bool = False + ): + """ + Get credible levels for all parameters. + + Adapted from Bilby. + + Parameters + ========== + keys: list, optional + A list of keys for which return the credible levels, if None, + defaults to search_parameter_keys + weighted: bool, optional + Whether to use sample weights in calculating credible level. + + Returns + ======= + credible_levels: dict + The credible levels at which the injected parameters are found. + """ + if keys is None: + keys = self.search_parameter_keys + if self.injection_parameters is None: + raise ( + TypeError, + "Result object has no 'injection_parameters'. " + "Cannot compute credible levels.", + ) + credible_levels = { + key: self.get_injection_credible_level(key, weighted=weighted) + for key in keys + if isinstance(self.injection_parameters.get(key, None), float) + } + return credible_levels + + def get_injection_credible_level(self, parameter: str, weighted: bool = False): + """ + Get the credible level of the injected parameter. + + Calculated as CDF(injection value). + + Adapted from Bilby. + + Parameters + ========== + parameter: str + Parameter to get credible level for + weighted: bool, optional + Whether to use sample weights in calculating credible level. + + Returns + ======= + float: credible level + """ + if self.injection_parameters is None: + raise ( + TypeError, + "Result object has no 'injection_parameters'. " + "Cannot compute credible levels.", + ) + theta = self._cleaned_samples() + + if weighted: + weights = theta["weights"] + else: + weights = np.ones(len(theta)) + + if parameter in theta and parameter in self.injection_parameters: + credible_level = sum( + np.array(theta[parameter].values < self.injection_parameters[parameter]) + * weights + ) / (sum(weights)) + return credible_level + else: + return np.nan + + +def make_pp_plot( + results: list[Result], + filename=None, + save=True, + confidence_interval=[0.68, 0.95, 0.997], + lines=None, + legend_fontsize="x-small", + keys=None, + title=True, + confidence_interval_alpha=0.1, + weighted: bool = False, + **kwargs, +): + """ + Make a P-P plot for a set of runs with injected signals. + + Adapted from Bilby. + + Parameters + ========== + results: list[Result] + A list of Result objects, each of these should have injected_parameters + filename: str, optional + The name of the file to save, the default is "outdir/pp.png" + save: bool, optional + Whether to save the file, default=True + confidence_interval: (float, list), optional + The confidence interval to be plotted, defaulting to 1-2-3 sigma + lines: list + If given, a list of matplotlib line formats to use, must be greater + than the number of parameters. + legend_fontsize: float + The font size for the legend + keys: list + A list of keys to use, if None defaults to search_parameter_keys + title: bool + Whether to add the number of results and total p-value as a plot title + confidence_interval_alpha: float, list, optional + The transparency for the background condifence interval + weighted: bool, optional + Whether to use weighted vs unweighted samples. It is useful to make PP plots + using unweighted samples to test networks without importance sampling. + kwargs: + Additional kwargs to pass to matplotlib.pyplot.plot + + Returns + ======= + fig, pvals: + matplotlib figure and a NamedTuple with attributes `combined_pvalue`, + `pvalues`, and `names`. + """ + if keys is None: + keys = results[0].search_parameter_keys + + credible_levels = list() + for i, result in enumerate(results): + credible_levels.append( + result.get_all_injection_credible_levels(keys, weighted=weighted) + ) + credible_levels = pd.DataFrame(credible_levels) + + if lines is None: + colors = ["C{}".format(i) for i in range(8)] + linestyles = ["-", "--", ":"] + lines = ["{}{}".format(a, b) for a, b in product(linestyles, colors)] + if len(lines) < len(credible_levels.keys()): + raise ValueError("Larger number of parameters than unique linestyles") + + x_values = np.linspace(0, 1, 1001) + + N = len(credible_levels) + fig, ax = plt.subplots() + + if isinstance(confidence_interval, float): + confidence_interval = [confidence_interval] + if isinstance(confidence_interval_alpha, float): + confidence_interval_alpha = [confidence_interval_alpha] * len( + confidence_interval + ) + elif len(confidence_interval_alpha) != len(confidence_interval): + raise ValueError( + "confidence_interval_alpha must have the same length as confidence_interval" + ) + + for ci, alpha in zip(confidence_interval, confidence_interval_alpha): + edge_of_bound = (1.0 - ci) / 2.0 + lower = scipy.stats.binom.ppf(1 - edge_of_bound, N, x_values) / N + upper = scipy.stats.binom.ppf(edge_of_bound, N, x_values) / N + # The binomial point percent function doesn't always return 0 @ 0, + # so set those bounds explicitly to be sure + lower[0] = 0 + upper[0] = 0 + ax.fill_between(x_values, lower, upper, alpha=alpha, color="k") + + pvalues = [] + print("Key: KS-test p-value") + for ii, key in enumerate(credible_levels): + pp = np.array( + [ + sum(credible_levels[key].values < xx) / len(credible_levels) + for xx in x_values + ] + ) + pvalue = scipy.stats.kstest(credible_levels[key], "uniform").pvalue + pvalues.append(pvalue) + print("{}: {}".format(key, pvalue)) + + try: + name = results[0].prior[key].latex_label + if name is None: + name = key + except (AttributeError, KeyError): + name = key + label = "{} ({:2.3f})".format(name, pvalue) + plt.plot(x_values, pp, lines[ii], label=label, **kwargs) + + Pvals = namedtuple("pvals", ["combined_pvalue", "pvalues", "names"]) + pvals = Pvals( + combined_pvalue=scipy.stats.combine_pvalues(pvalues)[1], + pvalues=pvalues, + names=list(credible_levels.keys()), + ) + print("Combined p-value: {}".format(pvals.combined_pvalue)) + + if title: + ax.set_title( + "N={}, p-value={:2.4f}".format(len(results), pvals.combined_pvalue) + + (", " "IS" if weighted else "") + ) + ax.set_xlabel("p") + ax.set_ylabel("CDF(p)") + ax.legend(handlelength=2, labelspacing=0.25, fontsize=legend_fontsize) + ax.set_xlim(0, 1) + ax.set_ylim(0, 1) + fig.tight_layout() + if save: + if filename is None: + filename = "outdir/pp.pdf" + fig.savefig(fname=filename) + + return fig, pvals + def check_equal_dict_of_arrays(a, b): if type(a) != type(b): diff --git a/dingo/gw/prior.py b/dingo/gw/prior.py index 9c6763604..cd300f046 100644 --- a/dingo/gw/prior.py +++ b/dingo/gw/prior.py @@ -74,14 +74,14 @@ def mean_std(self, keys=([]), sample_size=50000, force_numerical=False): std[key] = np.sqrt((p.maximum - p.minimum) ** 2.0 / 12.0).item() elif isinstance(p, Sine) and p.minimum == 0.0 and p.maximum == np.pi: mean[key] = np.pi / 2.0 - std[key] = np.sqrt(0.25 * (np.pi ** 2) - 2).item() + std[key] = np.sqrt(0.25 * (np.pi**2) - 2).item() elif ( isinstance(p, Cosine) and p.minimum == -np.pi / 2 and p.maximum == np.pi / 2 ): mean[key] = 0.0 - std[key] = np.sqrt(0.25 * (np.pi ** 2) - 2).item() + std[key] = np.sqrt(0.25 * (np.pi**2) - 2).item() else: estimation_keys.append(key) else: @@ -99,31 +99,28 @@ def mean_std(self, keys=([]), sample_size=50000, force_numerical=False): return mean, std -# TODO: Add latex labels, names - - default_extrinsic_dict = { - "dec": "bilby.core.prior.Cosine(minimum=-np.pi/2, maximum=np.pi/2)", - "ra": 'bilby.core.prior.Uniform(minimum=0., maximum=2*np.pi, boundary="periodic")', - "geocent_time": "bilby.core.prior.Uniform(minimum=-0.1, maximum=0.1)", - "psi": 'bilby.core.prior.Uniform(minimum=0.0, maximum=np.pi, boundary="periodic")', - "luminosity_distance": "bilby.core.prior.Uniform(minimum=100.0, maximum=6000.0)", + "dec": "bilby.core.prior.Cosine(minimum=-np.pi/2, maximum=np.pi/2, name='dec')", + "ra": 'bilby.core.prior.Uniform(minimum=0., maximum=2*np.pi, boundary="periodic", name="ra")', + "geocent_time": "bilby.core.prior.Uniform(minimum=-0.1, maximum=0.1, name='geocent_time')", + "psi": 'bilby.core.prior.Uniform(minimum=0.0, maximum=np.pi, boundary="periodic", name="psi")', + "luminosity_distance": "bilby.core.prior.Uniform(minimum=100.0, maximum=6000.0, name='luminosity_distance')", } default_intrinsic_dict = { - "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", - "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", - "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)", - "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)", + "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0, name='mass_1')", + "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0, name='mass_2')", + "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0, name='mass_ratio')", + "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0, name='chirp_mass')", "luminosity_distance": 1000.0, - "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", - "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', - "a_1": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)", - "a_2": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)", - "tilt_1": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", - "tilt_2": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", - "phi_12": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', - "phi_jl": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', + "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi, name='theta_jn')", + "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic", name="phase")', + "a_1": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99, name='a_1')", + "a_2": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99, name='a_2')", + "tilt_1": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi, name='tilt_1')", + "tilt_2": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi, name='tilt_2')", + "phi_12": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic", name="phi_12")', + "phi_jl": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic", name="phi_jl")', "geocent_time": 0.0, } diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py index 7e8d5a135..f0eccd4a8 100644 --- a/dingo/pipe/main.py +++ b/dingo/pipe/main.py @@ -19,9 +19,6 @@ BilbyPipeError, ) -from dingo.core.posterior_models.build_model import build_model_from_kwargs -from dingo.gw.domains import build_domain_from_model_metadata - from .dag_creator import generate_dag from .parser import create_parser From 603038a700157283d752782939b9245b7d540689 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Tue, 24 Jun 2025 12:11:25 +0100 Subject: [PATCH 09/48] Add PP plotting to dingo_pipe. To use, set plot-pp = true in INI file. --- dingo/pipe/dag_creator.py | 5 ++- dingo/pipe/main.py | 2 + dingo/pipe/nodes/plot_node.py | 26 +++++++++++ dingo/pipe/parser.py | 5 +++ dingo/pipe/plot.py | 14 +++--- dingo/pipe/pp_test.py | 81 +++++++++++++++++++++++++++++++++++ pyproject.toml | 3 +- 7 files changed, 128 insertions(+), 8 deletions(-) create mode 100644 dingo/pipe/pp_test.py diff --git a/dingo/pipe/dag_creator.py b/dingo/pipe/dag_creator.py index 8ed8c8eac..077864265 100644 --- a/dingo/pipe/dag_creator.py +++ b/dingo/pipe/dag_creator.py @@ -12,7 +12,7 @@ from .nodes.importance_sampling_node import ImportanceSamplingNode from .nodes.merge_node import MergeNode from .nodes.pe_summary_node import PESummaryNode -from .nodes.plot_node import PlotNode +from .nodes.plot_node import PlotNode, PlotPPNode from .nodes.sampling_node import SamplingNode logger.name = "dingo_pipe" @@ -157,6 +157,9 @@ def generate_dag(inputs): if inputs.plot_node_needed: plot_nodes_list.append(PlotNode(inputs, merged_node, dag=dag)) + if inputs.plot_pp: + PlotPPNode(inputs, merged_importance_sampling_node_list, dag=dag) + # # 6. PESummary # diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py index f0eccd4a8..f648326dc 100644 --- a/dingo/pipe/main.py +++ b/dingo/pipe/main.py @@ -292,6 +292,8 @@ def __init__( setattr(self, attr, getattr(args, attr)) if getattr(self, attr): self.plot_node_needed = True + + self.plot_pp = args.plot_pp # # # Set all other plotting options # for plot_attr in [ diff --git a/dingo/pipe/nodes/plot_node.py b/dingo/pipe/nodes/plot_node.py index 094f6b29e..99c68747e 100644 --- a/dingo/pipe/nodes/plot_node.py +++ b/dingo/pipe/nodes/plot_node.py @@ -33,3 +33,29 @@ def __init__(self, inputs, merged_node, dag): @property def executable(self): return self._get_executable_path("dingo_pipe_plot") + + +class PlotPPNode(BilbyPlotNode): + def __init__(self, inputs, merged_node_list, dag): + super(BilbyPlotNode, self).__init__(inputs) + self.dag = dag + self.request_cpus = 1 + self.job_name = f"{self.inputs.label}_plot_pp" + self.setup_arguments( + add_ini=False, add_unknown_args=False, add_command_line_args=False + ) + self.arguments.add_positional_argument(self.inputs.result_directory) + + if getattr(self, "disable_hdf5_locking", None): + self.extra_lines.append('environment = "HDF5_USE_FILE_LOCKING=FALSE"') + + self.process_node() + for node in merged_node_list: + self.job.add_parent(node.job) + + if self.inputs.simple_submission: + _strip_unwanted_submission_keys(self.job) + + @property + def executable(self): + return self._get_executable_path("dingo_pipe_pp_test") diff --git a/dingo/pipe/parser.py b/dingo/pipe/parser.py index aa54d47f9..7dbb8bad9 100644 --- a/dingo/pipe/parser.py +++ b/dingo/pipe/parser.py @@ -988,6 +988,11 @@ def create_parser(top_level=True): action="store_true", help="Create scatter plot of target versus proposal log probabilities", ) + output_parser.add_argument( + "--plot-pp", + action="store_true", + help="Create PP plot based on several injections.", + ) # output_parser.add_argument( # "--plot-marginal", # action="store_true", diff --git a/dingo/pipe/plot.py b/dingo/pipe/plot.py index ff6236b25..7dc8af756 100644 --- a/dingo/pipe/plot.py +++ b/dingo/pipe/plot.py @@ -1,5 +1,3 @@ -from pathlib import Path - from bilby_pipe.bilbyargparser import BilbyArgParser from bilby_pipe.utils import get_command_line_arguments, logger, parse_args @@ -22,10 +20,14 @@ def create_parser(): parser.add("--label", type=str, default="") # parser.add("--calibration", action="store_true", help="Generate calibration plot") parser.add("--corner", action="store_true", help="Generate corner plot") - parser.add("--weights", action="store_true", help="Generate plot of importance " - "weights") - parser.add("--log_probs", action="store_true", help="Generate plot of target" - "versus proposal log probability") + parser.add( + "--weights", action="store_true", help="Generate plot of importance " "weights" + ) + parser.add( + "--log_probs", + action="store_true", + help="Generate plot of target" "versus proposal log probability", + ) # parser.add("--marginal", action="store_true", help="Generate marginal plots") # parser.add("--skymap", action="store_true", help="Generate skymap") # parser.add("--waveform", action="store_true", help="Generate waveform") diff --git a/dingo/pipe/pp_test.py b/dingo/pipe/pp_test.py new file mode 100644 index 000000000..0716b24e6 --- /dev/null +++ b/dingo/pipe/pp_test.py @@ -0,0 +1,81 @@ +import argparse +from pathlib import Path + +from bilby_pipe.utils import logger + +from dingo.core.result import make_pp_plot +from dingo.gw.result import Result + +logger.name = "dingo_pipe" + + +def create_parser(): + parser = argparse.ArgumentParser( + prog="dingo_pipe PP test", + usage="Generates a pp plot from a directory containing a set of results", + ) + parser.add_argument("directory", help="Path to the result files") + parser.add_argument( + "--outdir", help="Path to output directory, defaults to input directory " + ) + parser.add_argument("--label", help="Additional label to use for output") + parser.add_argument( + "--print", action="store_true", help="Print the list of filenames used" + ) + parser.add_argument( + "-n", type=int, help="Number of samples to truncate to", default=None + ) + parser.add_argument( + "--filter", + type=str, + help="A string to match and filtering results. If not specified, will look " + "first for '*importance_sampling.hdf5', then '*sampling.hdf5'.", + default=None, + ) + return parser + + +def get_results_filenames(args): + p = Path(args.directory) + if args.filter: + results_files = list(p.glob(f"*{args.filter}*.hdf5")) + else: + results_files = list(p.glob("*importance_sampling.hdf5")) + if len(results_files) == 0: + results_files = list(p.glob("*sampling.hdf5")) + if len(results_files) == 0: + raise FileNotFoundError(f"No results found in path {args.directory}") + + if args.n is not None: + logger.info(f"Truncating to first {args.n} results") + results_files = results_files[: args.n] + return results_files + + +def get_basename(args): + if args.outdir is None: + args.outdir = args.directory + basename = f"{args.outdir}/" + if args.label is not None: + basename += f"{args.label}_" + return basename + + +def main(args=None): + if args is None: + args = create_parser().parse_args() + results_filenames = get_results_filenames(args) + results = [Result(file_name=f) for f in results_filenames] + basename = get_basename(args) + + logger.info("Generating PP plot.") + # keys = [ + # name + # for name, p in results[0].prior.items() + # if isinstance(p, str) or p.is_fixed is False + # ] + logger.info(f"Parameters = {results[0].search_parameter_keys}") + make_pp_plot(results, filename=f"{basename}pp.pdf") + if "weights" in results[0].samples: + logger.info("Generating IS PP plot.") + make_pp_plot(results, filename=f"{basename}pp_IS.pdf", weighted=True) diff --git a/pyproject.toml b/pyproject.toml index 12f699bff..2d4b4c9b6 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -100,8 +100,9 @@ dingo_merge_asd_datasets = "dingo.gw.noise.utils:merge_datasets_cli" dingo_pipe = "dingo.pipe.main:main" dingo_pipe_importance_sampling = "dingo.pipe.importance_sampling:main" dingo_pipe_generation = "dingo.pipe.data_generation:main" -dingo_pipe_sampling = "dingo.pipe.sampling:main" dingo_pipe_plot = "dingo.pipe.plot:main" +dingo_pipe_pp_test = "dingo.pipe.pp_test:main" +dingo_pipe_sampling = "dingo.pipe.sampling:main" dingo_pt_to_hdf5 = "dingo.core.utils.pt_to_hdf5:main" dingo_result = "dingo.pipe.dingo_result:main" dingo_train = "dingo.gw.training:train_local" From 7d2ff4b7a406c13f1cc02551c97862ff778cdbb9 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Tue, 24 Jun 2025 13:43:32 +0100 Subject: [PATCH 10/48] Enable latex labels for corner plots. --- dingo/core/result.py | 38 +++++++++++++++++++++++++++--------- dingo/core/utils/plotting.py | 9 ++++++++- 2 files changed, 37 insertions(+), 10 deletions(-) diff --git a/dingo/core/result.py b/dingo/core/result.py index 4821bcb15..e557f75da 100644 --- a/dingo/core/result.py +++ b/dingo/core/result.py @@ -13,7 +13,7 @@ from matplotlib import pyplot as plt from scipy.constants import golden from scipy.special import logsumexp -from bilby.core.prior import Constraint, DeltaFunction +from bilby.core.prior import Constraint, DeltaFunction, PriorDict from dingo.core.dataset import DingoDataset from dingo.core.density import train_unconditional_density_estimator @@ -640,6 +640,7 @@ def plot_corner( weights=[None, weights.to_numpy()], labels=["Dingo", "Dingo-IS"], filename=filename, + latex_labels_dict=get_latex_labels(self.prior), **kwargs, ) else: @@ -881,6 +882,7 @@ def make_pp_plot( ax.fill_between(x_values, lower, upper, alpha=alpha, color="k") pvalues = [] + latex_labels = get_latex_labels(results[0].prior) print("Key: KS-test p-value") for ii, key in enumerate(credible_levels): pp = np.array( @@ -892,14 +894,7 @@ def make_pp_plot( pvalue = scipy.stats.kstest(credible_levels[key], "uniform").pvalue pvalues.append(pvalue) print("{}: {}".format(key, pvalue)) - - try: - name = results[0].prior[key].latex_label - if name is None: - name = key - except (AttributeError, KeyError): - name = key - label = "{} ({:2.3f})".format(name, pvalue) + label = "{} ({:2.3f})".format(latex_labels[key], pvalue) plt.plot(x_values, pp, lines[ii], label=label, **kwargs) Pvals = namedtuple("pvals", ["combined_pvalue", "pvalues", "names"]) @@ -958,3 +953,28 @@ def freeze(d): elif isinstance(d, list): return tuple(freeze(value) for value in d) return d + + +def get_latex_labels(prior: PriorDict) -> dict: + """ + Get the latex labels for prior parameters. If no latex label exists, return the + parameter key. + + Parameters + ---------- + prior : PriorDict + + Returns + ------- + dict of latex labels + """ + labels = {} + for k, v in prior.items(): + try: + l = v.latex_label + if l is None: + l = k + except (AttributeError, KeyError): + l = k + labels[k] = l + return labels diff --git a/dingo/core/utils/plotting.py b/dingo/core/utils/plotting.py index 15073d90d..c0160472f 100644 --- a/dingo/core/utils/plotting.py +++ b/dingo/core/utils/plotting.py @@ -12,6 +12,7 @@ def plot_corner_multi( weights=None, labels=None, filename: str = "corner.pdf", + latex_labels_dict: dict = None, **kwargs, ): """ @@ -28,6 +29,8 @@ def plot_corner_multi( Labels for the posteriors. filename : str Where to save samples. + latex_labels_dict : dict + Dictionary of latex labels. Other Parameters ---------------- @@ -68,6 +71,10 @@ def plot_corner_multi( for p in samples[0].columns if p in set.intersection(*(set(s.columns) for s in samples)) ] + if latex_labels_dict: + parameter_labels = [latex_labels_dict.get(p, p) for p in common_parameters] + else: + parameter_labels = common_parameters fig = None handles = [] @@ -75,7 +82,7 @@ def plot_corner_multi( color = mpl.colors.rgb2hex(plt.get_cmap(cmap)(i)) fig = corner.corner( s[common_parameters].to_numpy(), - labels=common_parameters, + labels=parameter_labels, weights=w, color=color, no_fill_contours=True, From 421ca0818d964365d85ef970204e5b3d4224dda3 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Tue, 24 Jun 2025 13:59:27 +0100 Subject: [PATCH 11/48] Small fix for latex labels. --- dingo/core/result.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/dingo/core/result.py b/dingo/core/result.py index e557f75da..73d80f718 100644 --- a/dingo/core/result.py +++ b/dingo/core/result.py @@ -644,7 +644,13 @@ def plot_corner( **kwargs, ) else: - plot_corner_multi(theta, labels=["Dingo"], filename=filename, **kwargs) + plot_corner_multi( + theta, + labels=["Dingo"], + filename=filename, + latex_labels_dict=get_latex_labels(self.prior), + **kwargs, + ) def plot_log_probs(self, filename="log_probs.png"): """ From 01251c1489a0ecafb3376d2b90b80b4e89c03673 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Tue, 24 Jun 2025 17:16:08 +0100 Subject: [PATCH 12/48] Add 'name' attribute to priors in examples. --- dingo/core/result.py | 15 ++++++--------- dingo/pipe/data_generation.py | 1 - docs/source/example_toy_npe_model.md | 12 ++++++------ docs/source/training.md | 4 ++-- docs/source/waveform_dataset.ipynb | 12 ++++++------ examples/fmpe_model/train_settings.yaml | 4 ++-- .../fmpe_model/waveform_dataset_settings.yaml | 12 ++++++------ examples/gnpe_model/train_settings.yaml | 4 ++-- examples/gnpe_model/train_settings_init.yaml | 4 ++-- .../gnpe_model/waveform_dataset_settings.yaml | 12 ++++++------ examples/npe_model/train_settings.yaml | 4 ++-- examples/npe_model/waveform_dataset_settings.yaml | 8 ++++---- examples/toy_npe_model/train_settings.yaml | 4 ++-- .../toy_npe_model/waveform_dataset_settings.yaml | 8 ++++---- 14 files changed, 50 insertions(+), 54 deletions(-) diff --git a/dingo/core/result.py b/dingo/core/result.py index 73d80f718..466e0ecd8 100644 --- a/dingo/core/result.py +++ b/dingo/core/result.py @@ -13,7 +13,7 @@ from matplotlib import pyplot as plt from scipy.constants import golden from scipy.special import logsumexp -from bilby.core.prior import Constraint, DeltaFunction, PriorDict +from bilby.core.prior import Constraint, DeltaFunction, PriorDict, Prior from dingo.core.dataset import DingoDataset from dingo.core.density import train_unconditional_density_estimator @@ -963,8 +963,8 @@ def freeze(d): def get_latex_labels(prior: PriorDict) -> dict: """ - Get the latex labels for prior parameters. If no latex label exists, return the - parameter key. + Get the latex labels for prior parameters. If no latex label exists within the + prior object, try to choose based on parameter key. Finally, return the parameter key. Parameters ---------- @@ -976,11 +976,8 @@ def get_latex_labels(prior: PriorDict) -> dict: """ labels = {} for k, v in prior.items(): - try: - l = v.latex_label - if l is None: - l = k - except (AttributeError, KeyError): - l = k + l = v.latex_label + if l is None: + l = Prior._default_latex_labels.get(k, k) labels[k] = l return labels diff --git a/dingo/pipe/data_generation.py b/dingo/pipe/data_generation.py index 32dee08f0..20db8072f 100644 --- a/dingo/pipe/data_generation.py +++ b/dingo/pipe/data_generation.py @@ -9,7 +9,6 @@ from dingo.gw.data.event_dataset import EventDataset from dingo.gw.domains import UniformFrequencyDomain -from dingo.gw.gwutils import get_window, get_window_factor from dingo.pipe.parser import create_parser logger.name = "dingo_pipe" diff --git a/docs/source/example_toy_npe_model.md b/docs/source/example_toy_npe_model.md index 02fb4ca27..530a5cf9c 100644 --- a/docs/source/example_toy_npe_model.md +++ b/docs/source/example_toy_npe_model.md @@ -82,10 +82,10 @@ waveform_generator: # Dataset only samples over intrinsic parameters. Extrinsic parameters are chosen at train time. intrinsic_prior: - mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0) - mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0) - chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=15.0, maximum=100.0) - mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0) + mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0, name='mass_1') + mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0, name='mass_2') + chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=15.0, maximum=100.0, name='chirp_mass') + mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0, name='mass_ratio') phase: default chi_1: bilby.gw.prior.AlignedSpin(name='chi_1', a_prior=Uniform(minimum=0, maximum=0.9)) chi_2: bilby.gw.prior.AlignedSpin(name='chi_2', a_prior=Uniform(minimum=0, maximum=0.9)) @@ -188,9 +188,9 @@ data: extrinsic_prior: # Sampled at train time dec: default ra: default - geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10) + geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10, name='geocent_time') psi: default - luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0) + luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0, name='luminosity_distance') ref_time: 1126259462.391 inference_parameters: - chirp_mass diff --git a/docs/source/training.md b/docs/source/training.md index 31542681c..e48413fee 100644 --- a/docs/source/training.md +++ b/docs/source/training.md @@ -28,9 +28,9 @@ data: extrinsic_prior: dec: default ra: default - geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10) + geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10, name='geocent_time') psi: default - luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0) + luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0, name='luminosity_distance') ref_time: 1126259462.391 gnpe_time_shifts: kernel: bilby.core.prior.Uniform(minimum=-0.001, maximum=0.001) diff --git a/docs/source/waveform_dataset.ipynb b/docs/source/waveform_dataset.ipynb index 7a2b0db33..09a82914e 100644 --- a/docs/source/waveform_dataset.ipynb +++ b/docs/source/waveform_dataset.ipynb @@ -548,13 +548,13 @@ "\n", "# Dataset only samples over intrinsic parameters. Extrinsic parameters are chosen at train time.\n", "intrinsic_prior:\n", - " mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)\n", - " mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)\n", - " chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)\n", - " mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)\n", + " mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0, name='mass_1')\n", + " mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0, name='mass_2')\n", + " chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0, name='chirp_mass')\n", + " mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0, name='mass_ratio')\n", " phase: default\n", - " a_1: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)\n", - " a_2: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)\n", + " a_1: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99, name='a_1')\n", + " a_2: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99, name='a_2')\n", " tilt_1: default\n", " tilt_2: default\n", " phi_12: default\n", diff --git a/examples/fmpe_model/train_settings.yaml b/examples/fmpe_model/train_settings.yaml index 74ca0c828..560c27127 100644 --- a/examples/fmpe_model/train_settings.yaml +++ b/examples/fmpe_model/train_settings.yaml @@ -16,8 +16,8 @@ data: - L1 extrinsic_prior: dec: default - geocent_time: bilby.core.prior.Uniform(minimum=-0.03, maximum=0.03) - luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0) + geocent_time: bilby.core.prior.Uniform(minimum=-0.03, maximum=0.03, name='geocent_time') + luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0, name='luminosity_distance') psi: default ra: default inference_parameters: diff --git a/examples/fmpe_model/waveform_dataset_settings.yaml b/examples/fmpe_model/waveform_dataset_settings.yaml index 77e62dff0..e93c748c8 100644 --- a/examples/fmpe_model/waveform_dataset_settings.yaml +++ b/examples/fmpe_model/waveform_dataset_settings.yaml @@ -10,13 +10,13 @@ waveform_generator: spin_conversion_phase: 0.0 intrinsic_prior: - mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0) - mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0) - chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=20.0, maximum=120.0) - mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0) + mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0, name='mass_1') + mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0, name='mass_2') + chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=20.0, maximum=120.0, name='chirp_mass') + mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0, name='mass_ratio') phase: default - a_1: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99) - a_2: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99) + a_1: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99, name='a_1') + a_2: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99, name='a_2') tilt_1: default tilt_2: default phi_12: default diff --git a/examples/gnpe_model/train_settings.yaml b/examples/gnpe_model/train_settings.yaml index 55a9076d6..9fd1d96eb 100644 --- a/examples/gnpe_model/train_settings.yaml +++ b/examples/gnpe_model/train_settings.yaml @@ -16,9 +16,9 @@ data: extrinsic_prior: # Sampled at train time dec: default ra: default - geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10) + geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10, name='geocent_time') psi: default - luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0) + luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0, name='luminosity_distance') ref_time: 1126259462.391 gnpe_time_shifts: kernel: bilby.core.prior.Uniform(minimum=-0.001, maximum=0.001) diff --git a/examples/gnpe_model/train_settings_init.yaml b/examples/gnpe_model/train_settings_init.yaml index 39e72f214..676350392 100644 --- a/examples/gnpe_model/train_settings_init.yaml +++ b/examples/gnpe_model/train_settings_init.yaml @@ -16,9 +16,9 @@ data: extrinsic_prior: # Sampled at train time dec: default ra: default - geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10) + geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10, name='geocent_time') psi: default - luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0) + luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0, name='luminosity_distance') ref_time: 1126259462.391 inference_parameters: - H1_time diff --git a/examples/gnpe_model/waveform_dataset_settings.yaml b/examples/gnpe_model/waveform_dataset_settings.yaml index 4c9bbbf8d..dc56abfb8 100644 --- a/examples/gnpe_model/waveform_dataset_settings.yaml +++ b/examples/gnpe_model/waveform_dataset_settings.yaml @@ -12,13 +12,13 @@ waveform_generator: # Dataset only samples over intrinsic parameters. Extrinsic parameters are chosen at train time. intrinsic_prior: - mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0) - mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0) - chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=15.0, maximum=150.0) - mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0) + mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0, name='mass_1') + mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0, name='mass_2') + chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=15.0, maximum=150.0, name='chirp_mass') + mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0, name='mass_ratio') phase: default - a_1: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99) - a_2: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99) + a_1: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99, name='a_1') + a_2: bilby.core.prior.Uniform(minimum=0.0, maximum=0.99, name='a_2') tilt_1: default tilt_2: default phi_12: default diff --git a/examples/npe_model/train_settings.yaml b/examples/npe_model/train_settings.yaml index e81db84fa..e73308349 100644 --- a/examples/npe_model/train_settings.yaml +++ b/examples/npe_model/train_settings.yaml @@ -16,9 +16,9 @@ data: extrinsic_prior: # Sampled at train time dec: default ra: default - geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10) + geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10, name='geocent_time') psi: default - luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0) + luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0, name='luminosity_distance') ref_time: 1126259462.391 inference_parameters: - chirp_mass diff --git a/examples/npe_model/waveform_dataset_settings.yaml b/examples/npe_model/waveform_dataset_settings.yaml index d529f63bb..8aabec304 100644 --- a/examples/npe_model/waveform_dataset_settings.yaml +++ b/examples/npe_model/waveform_dataset_settings.yaml @@ -12,10 +12,10 @@ waveform_generator: # Dataset only samples over intrinsic parameters. Extrinsic parameters are chosen at train time. intrinsic_prior: - mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0) - mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0) - chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=15.0, maximum=150.0) - mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0) + mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0, name='mass_1') + mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=120.0, name='mass_2') + chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=15.0, maximum=150.0, name='chirp_mass') + mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0, name='mass_ratio') phase: default chi_1: bilby.gw.prior.AlignedSpin(name='chi_1', a_prior=Uniform(minimum=0, maximum=0.9)) chi_2: bilby.gw.prior.AlignedSpin(name='chi_2', a_prior=Uniform(minimum=0, maximum=0.9)) diff --git a/examples/toy_npe_model/train_settings.yaml b/examples/toy_npe_model/train_settings.yaml index 60f5f282b..e331cdc78 100644 --- a/examples/toy_npe_model/train_settings.yaml +++ b/examples/toy_npe_model/train_settings.yaml @@ -12,9 +12,9 @@ data: extrinsic_prior: # Sampled at train time dec: default ra: default - geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10) + geocent_time: bilby.core.prior.Uniform(minimum=-0.10, maximum=0.10, name='geocent_time') psi: default - luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0) + luminosity_distance: bilby.core.prior.Uniform(minimum=100.0, maximum=1000.0, name='luminosity_distance') ref_time: 1126259462.391 inference_parameters: - chirp_mass diff --git a/examples/toy_npe_model/waveform_dataset_settings.yaml b/examples/toy_npe_model/waveform_dataset_settings.yaml index 9969af980..aa27b5bfe 100644 --- a/examples/toy_npe_model/waveform_dataset_settings.yaml +++ b/examples/toy_npe_model/waveform_dataset_settings.yaml @@ -11,10 +11,10 @@ waveform_generator: # Dataset only samples over intrinsic parameters. Extrinsic parameters are chosen at train time. intrinsic_prior: - mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0) - mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0) - chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=15.0, maximum=100.0) - mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0) + mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0, name='mass_1') + mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0, name='mass_2') + chirp_mass: bilby.gw.prior.UniformInComponentsChirpMass(minimum=15.0, maximum=100.0, name='chirp_mass') + mass_ratio: bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0, name='mass_ratio') phase: default chi_1: bilby.gw.prior.AlignedSpin(name='chi_1', a_prior=Uniform(minimum=0, maximum=0.9)) chi_2: bilby.gw.prior.AlignedSpin(name='chi_2', a_prior=Uniform(minimum=0, maximum=0.9)) From 1954866e55a5407ad16211ff7aa77b3d9386938c Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Wed, 25 Jun 2025 15:44:44 +0100 Subject: [PATCH 13/48] Enable dingo signal generator for dingo_pipe injections. Set dingo-injections = true in INI file. Noise is still generated using bilby_pipe. --- dingo/gw/injection.py | 4 -- dingo/pipe/data_generation.py | 101 +++++++++++++++++++++++++++++- dingo/pipe/parser.py | 9 +++ tests/gw/test_waveform_dataset.py | 12 ++-- 4 files changed, 115 insertions(+), 11 deletions(-) diff --git a/dingo/gw/injection.py b/dingo/gw/injection.py index 6663a84c4..d7960c39b 100644 --- a/dingo/gw/injection.py +++ b/dingo/gw/injection.py @@ -111,10 +111,6 @@ def use_base_domain(self, value: bool): self.data_domain = self.data_domain.base_domain self._use_base_domain = True self._initialize_transform() - else: - print( - f"{type(self.data_domain)} has no base domain. Nothing to do." - ) else: raise NotImplementedError( "Cannot recover original domain from base domain alone." diff --git a/dingo/pipe/data_generation.py b/dingo/pipe/data_generation.py index 20db8072f..ccd34d549 100644 --- a/dingo/pipe/data_generation.py +++ b/dingo/pipe/data_generation.py @@ -5,10 +5,13 @@ from bilby_pipe.main import parse_args from bilby_pipe.utils import logger, convert_string_to_dict from bilby_pipe.data_generation import DataGenerationInput as BilbyDataGenerationInput +import lalsimulation as LS import numpy as np +from dingo.core.posterior_models.build_model import build_model_from_kwargs from dingo.gw.data.event_dataset import EventDataset from dingo.gw.domains import UniformFrequencyDomain +from dingo.gw.injection import Injection from dingo.pipe.parser import create_parser logger.name = "dingo_pipe" @@ -102,6 +105,7 @@ def __init__(self, args, unknown_args, create_data=True): self.mode_array = None self.waveform_arguments_dict = None self.numerical_relativity_file = args.numerical_relativity_file + self.dingo_injection = args.dingo_injection self.injection_waveform_approximant = args.injection_waveform_approximant self.frequency_domain_source_model = "lal_binary_black_hole" # self.conversion_function = args.conversion_function @@ -153,7 +157,101 @@ def __init__(self, args, unknown_args, create_data=True): self.plot_injection = args.plot_injection if create_data: - self.create_data(args) + if self.dingo_injection: + self.create_data_dingo_injection(args) + else: + self.create_data(args) + + def create_data_dingo_injection(self, args): + """Executes create_data but turns off any requested injections. This is + intended to create noise-only datasets, to be used with the DINGO signal + generator.""" + # Save values of relevant args. + injection = args.injection + injection_file = args.injection_file + injection_dict = args.injection_dict + + args.injection = False + args.injection_file = None + args.injection_dict = None + + # Create noise. + self.create_data(args) + + # Reset args. + args.injection = injection + args.injection_file = injection_file + args.injection_dict = injection_dict + self.injection = args.injection + self.injection_file = args.injection_file + self.injection_dict = args.injection_dict + + if self.injection: + self._inject_dingo_signal(args) + + def _inject_dingo_signal(self, args): + """Generate a GW signal using the dingo.gw.injection class and add it to the + interferometer strain data. Also compute SNRs and store them.""" + try: + model = build_model_from_kwargs( + filename=args.model, device="meta", load_training_info=False + ) + except RuntimeError: + # 'meta' is not supported by older version of python / torch + model = build_model_from_kwargs( + filename=args.model, device="cpu", load_training_info=False + ) + + injection = Injection.from_posterior_model_metadata(model.metadata) + injection.use_base_domain = True # Do not generate MFD signals. + injection.t_ref = self.trigger_time + injection._initialize_transform() + + # Possibly update waveform generator based on supplied settings. Note that + # default values will have been set by dingo_pipe from the DINGO model. + waveform_arguments = self.get_injection_waveform_arguments() + injection.f_ref = waveform_arguments["reference_frequency"] + injection.waveform_generator.approximant = LS.GetApproximantFromString( + waveform_arguments["waveform_approximant"] + ) + injection.waveform_generator.approximant_str = waveform_arguments[ + "waveform_approximant" + ] + + logger.info("Injecting waveform from DINGO with ") + logger.info(f"data_domain = {injection.data_domain.domain_dict}") + for prop in [ + "t_ref", + ]: + logger.info(f"{prop} = {getattr(injection, prop)}") + for prop in [ + "approximant_str", + "f_ref", + "f_start", + ]: + logger.info(f"{prop} = {getattr(injection.waveform_generator, prop)}") + + # Generate signal + self.injection_parameters = self.injection_df.iloc[self.idx].to_dict() + theta = self.injection_parameters.copy() + theta["geocent_time"] -= self.trigger_time + signal = injection.signal(theta) + + # Add signal to interferometer data + domain = injection.data_domain + for ifo in self.interferometers: + s = signal["waveform"][ifo.name] + s = domain.time_translate_data(s, -self.post_trigger_duration) + # Interferometer data extends to Nyquist = f_max / 2, so pad signal. Will + # be truncated when saved as HDF5. + s = np.pad( + s, (0, len(ifo.strain_data.frequency_domain_strain) - len(domain)) + ) + ifo.strain_data.frequency_domain_strain += s + ifo.meta_data["optimal_SNR"] = np.sqrt( + ifo.optimal_snr_squared(signal=s).real + ) + ifo.meta_data["matched_filter_SNR"] = ifo.matched_filter_snr(signal=s) def save_hdf5(self): """ @@ -233,6 +331,7 @@ def save_hdf5(self): if self.injection: settings["injection_parameters"] = self.injection_parameters.copy() + settings["dingo_injection"] = self.dingo_injection # Dingo and Bilby have different geocent_time conventions. settings["injection_parameters"]["geocent_time"] -= self.trigger_time settings["optimal_SNR"] = { diff --git a/dingo/pipe/parser.py b/dingo/pipe/parser.py index 7dbb8bad9..be7f5d2ae 100644 --- a/dingo/pipe/parser.py +++ b/dingo/pipe/parser.py @@ -489,6 +489,15 @@ def create_parser(top_level=True): "Repeated entries will be ignored." ), ) + injection_parser.add( + "--dingo-injection", + action="store_true", + help=( + "If true, use the DINGO Injection class to generate the injection signal. " + "Otherwise, use bilby_pipe. When using DINGO for injections, the noise is " + "still generated using bilby_pipe. Defaults to false." + ), + ) injection_parser.add( "--injection-waveform-approximant", type=nonestr, diff --git a/tests/gw/test_waveform_dataset.py b/tests/gw/test_waveform_dataset.py index b51df355e..d5150f5c6 100644 --- a/tests/gw/test_waveform_dataset.py +++ b/tests/gw/test_waveform_dataset.py @@ -30,13 +30,13 @@ # settings for intrinsic prior over parameters intrinsic_prior: # prior for non-fixed parameters - mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0) - mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0) - mass_ratio: bilby.core.prior.Uniform(minimum=0.125, maximum=1.0) - chirp_mass: bilby.core.prior.Uniform(minimum=25.0, maximum=100.0) + mass_1: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0, name='mass_1') + mass_2: bilby.core.prior.Constraint(minimum=10.0, maximum=80.0, name='mass_2') + mass_ratio: bilby.core.prior.Uniform(minimum=0.125, maximum=1.0, name='mass_ratio') + chirp_mass: bilby.core.prior.Uniform(minimum=25.0, maximum=100.0, name='chirp_mass') phase: default - a_1: bilby.core.prior.Uniform(minimum=0.0, maximum=0.88) - a_2: bilby.core.prior.Uniform(minimum=0.0, maximum=0.88) + a_1: bilby.core.prior.Uniform(minimum=0.0, maximum=0.88, name='a_1') + a_2: bilby.core.prior.Uniform(minimum=0.0, maximum=0.88, name='a_2') tilt_1: default tilt_2: default phi_12: default From 91fb671a955ffe0a7aeefdd8d6da47938d309809 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Wed, 25 Jun 2025 16:17:30 +0100 Subject: [PATCH 14/48] Ensure generation-seed is used in picking ASD. --- dingo/gw/noise/asd_dataset.py | 8 ++++++-- dingo/pipe/main.py | 5 ++++- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/dingo/gw/noise/asd_dataset.py b/dingo/gw/noise/asd_dataset.py index ae7b6c40f..7006c5f05 100644 --- a/dingo/gw/noise/asd_dataset.py +++ b/dingo/gw/noise/asd_dataset.py @@ -2,6 +2,8 @@ from pathlib import Path from typing import Iterable, Optional +import numpy as np + from dingo.gw.domains import build_domain, UniformFrequencyDomain from dingo.gw.domains.base_frequency_domain import BaseFrequencyDomain from dingo.gw.gwutils import * @@ -194,9 +196,11 @@ def sample_random_asds(self, n: Optional[int] = None) -> dict[str, np.ndarray]: else: return {k: v[np.random.choice(len(v), n)] for k, v in self.asds.items()} - def save_psd(self, directory, ifo_name, idx: Optional[int] = None): + def save_psd(self, directory, ifo_name, idx: Optional[int] = None, rng=None): + if rng is None: + rng = np.random.default_rng() if idx is None: - idx = np.random.choice(len(self.asds[ifo_name])) + idx = rng.choice(len(self.asds[ifo_name])) directory = Path(directory) directory.mkdir(exist_ok=True) psd_path = directory / f"{ifo_name}_{idx}_psd.txt" diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py index f648326dc..c60011077 100644 --- a/dingo/pipe/main.py +++ b/dingo/pipe/main.py @@ -5,6 +5,8 @@ import json import os +import numpy as np + import dingo.pipe.create_injections from bilby.core.prior import PriorDict @@ -132,8 +134,9 @@ def fill_in_arguments_from_model(args): del domain_dict["window_factor"] asd_dataset.update_domain(domain_dict) psd_dict = {} + rng = np.random.default_rng(args.generation_seed) for ifo_name in args.detectors: - psd_path = asd_dataset.save_psd(args.outdir, ifo_name) + psd_path = asd_dataset.save_psd(args.outdir, ifo_name, rng=rng) psd_dict[ifo_name] = str(psd_path) args.asd_dataset = None args.psd_dict = str(psd_dict) From 60b0fe221bec7d115d1e4b5648e9f6b6b8c53ed6 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Wed, 25 Jun 2025 16:54:38 +0100 Subject: [PATCH 15/48] Update docs. --- docs/source/dingo_pipe.md | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/docs/source/dingo_pipe.md b/docs/source/dingo_pipe.md index d9ebf3ad4..b84057e9b 100644 --- a/docs/source/dingo_pipe.md +++ b/docs/source/dingo_pipe.md @@ -85,6 +85,13 @@ The prepared event data and ASD are stored in a {py:class}`dingo.gw.data.event_d Dingo models are typically trained using Welch PSDs. For this reason we do not recommend using a BayesWave PSD for initial sampling. Rather, a BayesWave PSD should be specified within the `importance_sampling_updates` dictionary, so that it will be used during importance sampling. ``` +### Injections + +Injections mirror the [approach of `bilby_pipe`](https://lscsoft.docs.ligo.org/bilby_pipe/master/injections.html) with a few modifications: +* By default, injections are generated using `bilby_pipe`. To use signal waveforms generated by the Dingo code (using the [Injections class](inference.md#injections)) set `dingo-injection = True`. Note that Dingo and Bilby injections might differ depending on the `spin_conversion_phase` setting, which is only available in Dingo. +* Setting `zero-noise = True` will create a zero-noise injection, which is out-of-distribution with respect to Dingo training data. It should therefore be used with caution. +* When running many injections `n-simulation > 1`, one may want to produce a PP plot. This can be done automatically by setting `plot-pp = true`. + ## Sampling The next step is sampling from the Dingo model. The model is loaded into a [GWSampler](dingo.gw.inference.gw_samplers.GWSampler) or [GWSamplerGNPE](dingo.gw.inference.gw_samplers.GWSamplerGNPE) object. (If using [GNPE](gnpe) it is necessary to specify a `model-init`.) The Sampler `context` is then set from the EventDataset prepared in the previous step. `num-samples` samples are then generated in batches of size `batch-size`. The samples (and context) are stored in a [Result](dingo.gw.result.Result) object and saved in HDF5 format. @@ -127,6 +134,8 @@ spline-calibration-curves The standard Result [plots](result.md#plotting) are turned on using the `plot-corner`, `plot-weights`, and `plot-log-probs` flags. +PP plots are produced using the `plot-pp` flag. Generally this should be done when analyzing many injections (`n-simulation > 1`). Dingo will produce two PP plots, one using unweighted and another using weighted samples (from importance sampling, labeled "IS"). The PP plot for unweighted samples is useful for gauging Dingo network performance, whereas the weighted PP plot tests consistency of the likelihood. + ## Additional options extra-lines From 8a9cf98a6e3db7c8acddffb7eee8182d18a09f96 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Thu, 26 Jun 2025 10:54:44 +0100 Subject: [PATCH 16/48] Improve a docstring. --- dingo/pipe/data_generation.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/dingo/pipe/data_generation.py b/dingo/pipe/data_generation.py index ccd34d549..c17cfde93 100644 --- a/dingo/pipe/data_generation.py +++ b/dingo/pipe/data_generation.py @@ -163,9 +163,14 @@ def __init__(self, args, unknown_args, create_data=True): self.create_data(args) def create_data_dingo_injection(self, args): - """Executes create_data but turns off any requested injections. This is - intended to create noise-only datasets, to be used with the DINGO signal - generator.""" + """Adaptation of create_data to use Dingo signal models rather than Bilby. + + First, executes create_data but without any requested injections. This creates + a noise-only dataset. + + Second, calls _inject_dingo_signal to generate the Dingo signal waveform and + add it to the noisy data within the interferometers. + """ # Save values of relevant args. injection = args.injection injection_file = args.injection_file From 9619135cdded2ae5dff9aaed3064f6890d826296 Mon Sep 17 00:00:00 2001 From: Stephen Green Date: Tue, 1 Jul 2025 14:29:51 +0100 Subject: [PATCH 17/48] Update dingo/pipe/data_generation.py Co-authored-by: Nihar Gupte <40392612+nihargupte-ph@users.noreply.github.com> --- dingo/pipe/data_generation.py | 1 + 1 file changed, 1 insertion(+) diff --git a/dingo/pipe/data_generation.py b/dingo/pipe/data_generation.py index c17cfde93..f8ebafb16 100644 --- a/dingo/pipe/data_generation.py +++ b/dingo/pipe/data_generation.py @@ -392,6 +392,7 @@ def _set_interferometers_from_gaussian_noise(self): # This is a hack to set the window factor. It ensures also that the SNRs # are calculated correctly. td_strain = ifo.time_domain_strain + # TODO: correct for window factor changes https://git.ligo.org/pe/pe-group-coordination/-/issues/1 ifo.strain_data.time_domain_window(roll_off=self.tukey_roll_off) ifo.strain_data.frequency_domain_strain = ( ifo.strain_data.frequency_domain_strain From d5e91748a0f0ac28fc175e3b6f4727ff0e77ba38 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 8 Jul 2025 11:51:47 +0200 Subject: [PATCH 18/48] updates to make injections compatible with gwsignal and save the data dump --- dingo/pipe/data_generation.py | 91 +++++++++++++++++++++++++++++++---- dingo/pipe/main.py | 45 ++++++++++------- dingo/pipe/parser.py | 74 ++++++++++++++++------------ dingo/pipe/utils.py | 4 ++ 4 files changed, 154 insertions(+), 60 deletions(-) diff --git a/dingo/pipe/data_generation.py b/dingo/pipe/data_generation.py index f396b538c..26b7fcc11 100644 --- a/dingo/pipe/data_generation.py +++ b/dingo/pipe/data_generation.py @@ -2,17 +2,26 @@ import sys from bilby_pipe.input import Input -from bilby_pipe.main import parse_args -from bilby_pipe.utils import logger, convert_string_to_dict from bilby_pipe.data_generation import DataGenerationInput as BilbyDataGenerationInput +from bilby.core.prior import PriorDict +from bilby_pipe.utils import ( + parse_args, + logger, + convert_string_to_dict, + convert_prior_string_input, + BilbyPipeError, +) import lalsimulation as LS import numpy as np from dingo.core.posterior_models.build_model import build_model_from_kwargs +from dingo.gw.prior import build_prior_with_defaults +from dingo.gw.gwutils import get_extrinsic_prior_dict from dingo.gw.data.event_dataset import EventDataset from dingo.gw.domains import UniformFrequencyDomain from dingo.gw.injection import Injection from dingo.pipe.parser import create_parser +from dingo.pipe.main import fill_in_arguments_from_model logger.name = "dingo_pipe" @@ -20,7 +29,6 @@ class DataGenerationInput(BilbyDataGenerationInput): def __init__(self, args, unknown_args, create_data=True): Input.__init__(self, args, unknown_args) - # Generic initialisation self.meta_data = dict( command_line_args=args.__dict__, @@ -46,12 +54,12 @@ def __init__(self, args, unknown_args, create_data=True): # Prior arguments # self.reference_frame = args.reference_frame - # self.time_reference = args.time_reference + # self.time_reference = "geocent" # DINGO mod used for saving data dump # self.phase_marginalization = args.phase_marginalization - # self.prior_file = args.prior_file - # self.prior_dict = args.prior_dict - self.deltaT = args.deltaT - # self.default_prior = args.default_prior + self.prior_dict = args.prior_dict + self.default_prior = "BBHPriorDict" + self.time_reference = "geocent" + self.prior_dict_updates = args.prior_dict_updates # Whether to generate data for importance sampling. This must be done when # desired data settings differ from those used for network training. If this is @@ -104,11 +112,23 @@ def __init__(self, args, unknown_args, create_data=True): self.pn_phase_order = -1 self.pn_amplitude_order = 0 self.mode_array = None - self.waveform_arguments_dict = None + self.waveform_arguments_dict = args.waveform_arguments_dict self.numerical_relativity_file = args.numerical_relativity_file self.dingo_injection = args.dingo_injection self.injection_waveform_approximant = args.injection_waveform_approximant - self.frequency_domain_source_model = "lal_binary_black_hole" + if args.injection_waveform_approximant in ["SEOBNRv5PHM", "SEOBNRv5EHM", "SEOBNRv5HM"]: + self.injection_frequency_domain_source_model = "gwsignal_binary_black_hole" + self.frequency_domain_source_model = "gwsignal_binary_black_hole" + else: + self.injection_frequency_domain_source_model = "lal_binary_black_hole" + self.frequency_domain_source_model = "lal_binary_black_hole" + + # DINGO mod + self.save_bilby_data_dump = args.save_bilby_data_dump + if self.save_bilby_data_dump: + self.time_reference = args.time_reference + self.deltaT = args.deltaT + # self.conversion_function = args.conversion_function # self.generation_function = args.generation_function # self.likelihood_type = args.likelihood_type @@ -267,6 +287,24 @@ def save_hdf5(self): for easy reading by pesummary and Bilby. """ + if self.save_bilby_data_dump: + # this is needed because we want bilby to use the updated DINGO + # prior + self.prior_dict = self.priors + self.likelihood_type = "GravitationalWaveTransient" + self.calibration_marginalization = False + self.phase_marginalization = False + self.time_marginalization = False + self.distance_marginalization = False + self.number_of_response_curves = 0 + self._distance_marginalization_lookup_table = None + self.reference_frame = "sky" + self.fiducial_parameters = None + self.update_fiducial_parameters = None + self.epsilon = None + self.jitter_time = True + self.save_data_dump() + # PSD and strain data. data = {"waveform": {}, "asds": {}} # TODO: Rename these keys. for ifo in self.interferometers: @@ -400,6 +438,39 @@ def _set_interferometers_from_gaussian_noise(self): * np.sqrt(ifo.strain_data.window_factor) ) + @property + def prior_dict_updates(self): + """The input prior_dict from the ini (if given) + + Note, this is not the bilby prior (see self.priors for that), this is + a key-val dictionary where the val's are strings which are converting + into bilby priors in `_get_prior + """ + return self._prior_dict_updates + + @prior_dict_updates.setter + def prior_dict_updates(self, prior_dict_updates): + if isinstance(prior_dict_updates, dict): + prior_dict_updates = prior_dict_updates + elif isinstance(prior_dict_updates, str): + prior_dict_updates = convert_prior_string_input(prior_dict_updates) + elif prior_dict_updates is None: + self._prior_dict_updates = None + return + else: + raise BilbyPipeError( + f"prior_dict_updates={prior_dict_updates} not " f"understood" + ) + + self._prior_dict_updates = { + self._convert_prior_dict_key(key): val + for key, val in prior_dict_updates.items() + } + + def _get_priors(self, add_time=True): + priors = super()._get_priors(add_time=add_time) + priors.update(PriorDict(self.prior_dict_updates)) + return priors def create_generation_parser(): """Data generation parser creation""" diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py index c60011077..a761969be 100644 --- a/dingo/pipe/main.py +++ b/dingo/pipe/main.py @@ -23,6 +23,7 @@ from .dag_creator import generate_dag from .parser import create_parser +from .utils import dict_to_string from ..gw.domains.build_domain import build_domain_from_model_metadata from dingo.core.posterior_models.build_model import build_model_from_kwargs @@ -32,18 +33,19 @@ logger.name = "dingo_pipe" -def fill_in_arguments_from_model(args): - if args.prior_dict is not None: - raise ValueError( - "Do not specify prior-dict in INI file. This is obtained from " - "the DINGO model. To update the prior, specify " - "prior-dict-updates." - ) - if args.model_reference_time is not None: - raise ValueError( - "Do not specify model-reference-time in INI file. This is obtained from the " - "DINGO model." - ) +def fill_in_arguments_from_model(args, perform_arg_checks=True): + if perform_arg_checks: + if args.prior_dict is not None: + raise ValueError( + "Do not specify prior-dict in INI file. This is obtained from " + "the DINGO model. To update the prior, specify " + "prior-dict-updates." + ) + if args.model_reference_time is not None: + raise ValueError( + "Do not specify model-reference-time in INI file. This is obtained from the " + "DINGO model." + ) logger.info(f"Loading dingo model from {args.model} in order to access settings.") @@ -91,7 +93,7 @@ def fill_in_arguments_from_model(args): ], "deltaT": deltaT, "Toffset": Toffset, - "prior_dict": prior, + "prior_dict": dict_to_string(prior), "model_reference_time": model_metadata["train_settings"]["data"]["ref_time"], } @@ -267,11 +269,18 @@ def __init__( # ) self.generation_seed = args.generation_seed if self.importance_sampling_updates and self.gaussian_noise: - raise ValueError( - "Cannot update data for importance sampling if using " - "simulated Gaussian noise. This risks inconsistent noise " - "realizations." - ) + if not ( + list(self.importance_sampling_updates.keys()) == ["minimum_frequency"] + ): + raise ValueError( + "Cannot update data for importance sampling if using " + "simulated Gaussian noise. This risks inconsistent noise " + "realizations." + ) + else: + # only allow for changing the minimum frequency (fstart) for the + # injection waveform. But this should not be an importance sampling update + del self.importance_sampling_updates["minimum_frequency"] self.importance_sample = args.importance_sample diff --git a/dingo/pipe/parser.py b/dingo/pipe/parser.py index be7f5d2ae..43a557e59 100644 --- a/dingo/pipe/parser.py +++ b/dingo/pipe/parser.py @@ -524,6 +524,15 @@ def create_parser(top_level=True): "injection only" ), ) + injection_parser.add( + "--save-bilby-data-dump", + type=bool, + default=False, + help=( + "If given, will also save a data dump consistent with the DINGO injection." + "This is useful when comparing with Bilby" + ), + ) submission_parser = parser.add_argument_group( title="Job submission arguments", @@ -841,10 +850,10 @@ def create_parser(top_level=True): "for more details." ), ) - # likelihood_parser = parser.add_argument_group( - # title="Likelihood arguments", - # description="Options for setting up the likelihood.", - # ) + likelihood_parser = parser.add_argument_group( + title="Likelihood arguments", + description="Options for setting up the likelihood.", + ) # likelihood_parser.add( # "--distance-marginalization", # action="store_true", @@ -882,12 +891,12 @@ def create_parser(top_level=True): # type=str, # help="Reference frame for the sky parameterisation, either 'sky' (default) or, e.g., 'H1L1'", # ) - # likelihood_parser.add( - # "--time-reference", - # default="geocent", - # type=str, - # help="Time parameter to sample in, either 'geocent' (default) or, e.g., 'H1'", - # ) + likelihood_parser.add( + "--time-reference", + default="geocent", + type=str, + help="Time parameter to sample in, either 'geocent' (default) or, e.g., 'H1'", + ) # likelihood_parser.add( # "--likelihood-type", # default="GravitationalWaveTransient", @@ -1097,16 +1106,16 @@ def create_parser(top_level=True): prior_parser = parser.add_argument_group( title="Prior arguments", description="Specify the prior settings." ) - # prior_parser.add( - # "--default-prior", - # default="BBHPriorDict", - # type=str, - # help=( - # "The name of the prior set to base the prior on. Can be one of" - # "[PriorDict, BBHPriorDict, BNSPriorDict, CalibrationPriorDict]" - # "or a python path to a bilby prior class available in the user's installation." - # ), - # ) + prior_parser.add( + "--default-prior", + default="BBHPriorDict", + type=str, + help=( + "The name of the prior set to base the prior on. Can be one of" + "[PriorDict, BBHPriorDict, BNSPriorDict, CalibrationPriorDict]" + "or a python path to a bilby prior class available in the user's installation." + ), + ) prior_parser.add( "--deltaT", type=nonefloat, @@ -1141,9 +1150,9 @@ def create_parser(top_level=True): ), ) prior_parser_main = prior_parser.add_mutually_exclusive_group() - # prior_parser_main.add( - # "--prior-file", type=nonestr, default=None, help="The prior file" - # ) + prior_parser_main.add( + "--prior-file", type=nonestr, default=None, help="The prior file" + ) prior_parser_main.add( "--prior-dict", type=nonestr, @@ -1285,15 +1294,16 @@ def create_parser(top_level=True): "https://git.ligo.org/waveforms/lvcnr-lfs for examples" ), ) - # waveform_parser.add( - # "--waveform-arguments-dict", - # default=None, - # type=nonestr, - # help=( - # "A dictionary of arbitrary additional waveform-arguments to pass" - # " to the bilby waveform generator's `waveform_arguments`" - # ), - # ) + waveform_parser.add( + "--waveform-arguments-dict", + default=None, + type=nonestr, + help=( + "A dictionary of arbitrary additional waveform-arguments to pass" + " to the bilby waveform generator's `waveform_arguments`. Only used " + "for injections" + ), + ) # waveform_parser.add( # "--mode-array", # default=None, diff --git a/dingo/pipe/utils.py b/dingo/pipe/utils.py index 273c9cb76..18495b8f9 100644 --- a/dingo/pipe/utils.py +++ b/dingo/pipe/utils.py @@ -6,3 +6,7 @@ def _strip_unwanted_submission_keys(job): if not line.startswith("priority") and not line.startswith("accounting_group") and not line.startswith("ENV GET HTGETTOKENOPTS"): extra_lines.append(line) job.extra_lines = extra_lines + +def dict_to_string(d): + """Convert a dictionary to a string representation.""" + return "{" + ", ".join(f"{key}:{value}" for key, value in d.items()) + "}" \ No newline at end of file From ed31be79d9dcde18a065ca65685678d440d48d7d Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 8 Jul 2025 13:35:13 +0200 Subject: [PATCH 19/48] disallowed for minimum-frequency changes to waveform and only allowing to change it via args.injection-waveform-arguments --- dingo/pipe/data_generation.py | 5 ++++- dingo/pipe/main.py | 17 +++++------------ dingo/pipe/parser.py | 20 ++++++++++---------- 3 files changed, 19 insertions(+), 23 deletions(-) diff --git a/dingo/pipe/data_generation.py b/dingo/pipe/data_generation.py index 26b7fcc11..0f35eb58f 100644 --- a/dingo/pipe/data_generation.py +++ b/dingo/pipe/data_generation.py @@ -112,7 +112,9 @@ def __init__(self, args, unknown_args, create_data=True): self.pn_phase_order = -1 self.pn_amplitude_order = 0 self.mode_array = None - self.waveform_arguments_dict = args.waveform_arguments_dict + # don't set self.waveform_arguments_dict, it will be updated later by injection_waveform_arguments + self.waveform_arguments_dict = None + self.injection_waveform_arguments = args.injection_waveform_arguments self.numerical_relativity_file = args.numerical_relativity_file self.dingo_injection = args.dingo_injection self.injection_waveform_approximant = args.injection_waveform_approximant @@ -290,6 +292,7 @@ def save_hdf5(self): if self.save_bilby_data_dump: # this is needed because we want bilby to use the updated DINGO # prior + self.waveform_arguments_dict = self.injection_waveform_arguments self.prior_dict = self.priors self.likelihood_type = "GravitationalWaveTransient" self.calibration_marginalization = False diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py index a761969be..e904d4630 100644 --- a/dingo/pipe/main.py +++ b/dingo/pipe/main.py @@ -269,18 +269,11 @@ def __init__( # ) self.generation_seed = args.generation_seed if self.importance_sampling_updates and self.gaussian_noise: - if not ( - list(self.importance_sampling_updates.keys()) == ["minimum_frequency"] - ): - raise ValueError( - "Cannot update data for importance sampling if using " - "simulated Gaussian noise. This risks inconsistent noise " - "realizations." - ) - else: - # only allow for changing the minimum frequency (fstart) for the - # injection waveform. But this should not be an importance sampling update - del self.importance_sampling_updates["minimum_frequency"] + raise ValueError( + "Cannot update data for importance sampling if using " + "simulated Gaussian noise. This risks inconsistent noise " + "realizations." + ) self.importance_sample = args.importance_sample diff --git a/dingo/pipe/parser.py b/dingo/pipe/parser.py index 43a557e59..8463eecf2 100644 --- a/dingo/pipe/parser.py +++ b/dingo/pipe/parser.py @@ -1294,16 +1294,16 @@ def create_parser(top_level=True): "https://git.ligo.org/waveforms/lvcnr-lfs for examples" ), ) - waveform_parser.add( - "--waveform-arguments-dict", - default=None, - type=nonestr, - help=( - "A dictionary of arbitrary additional waveform-arguments to pass" - " to the bilby waveform generator's `waveform_arguments`. Only used " - "for injections" - ), - ) + # waveform_parser.add( + # "--waveform-arguments-dict", + # default=None, + # type=nonestr, + # help=( + # "A dictionary of arbitrary additional waveform-arguments to pass" + # " to the bilby waveform generator's `waveform_arguments`. Only used " + # "for injections" + # ), + # ) # waveform_parser.add( # "--mode-array", # default=None, From 5cf2018c7cc9b673676680e1314e706974a1869b Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 21 Jul 2025 15:42:02 +0200 Subject: [PATCH 20/48] added data dump to main.py --- dingo/pipe/main.py | 1 + 1 file changed, 1 insertion(+) diff --git a/dingo/pipe/main.py b/dingo/pipe/main.py index e904d4630..d4e884e6a 100644 --- a/dingo/pipe/main.py +++ b/dingo/pipe/main.py @@ -262,6 +262,7 @@ def __init__( self.injection_numbers = args.injection_numbers self.injection_file = args.injection_file self.injection_dict = args.injection_dict + self.save_bilby_data_dump = args.save_bilby_data_dump # self.injection_waveform_arguments = args.injection_waveform_arguments # self.injection_waveform_approximant = args.injection_waveform_approximant # self.injection_frequency_domain_source_model = ( From 6d1c7cf6d9222042d54b9ad4a33349c6162312b5 Mon Sep 17 00:00:00 2001 From: Nihar Gupte <40392612+nihargupte-ph@users.noreply.github.com> Date: Tue, 29 Jul 2025 11:33:15 +0200 Subject: [PATCH 21/48] Update pytest.yml Adding BLAS and MKL threading options so that the scipy/numpy svd function doesn't segfault --- .github/workflows/pytest.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index bd64e0682..762903751 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -13,6 +13,10 @@ on: jobs: pytest: runs-on: ubuntu-latest + env: + OPENBLAS_NUM_THREADS: 1 + MKL_NUM_THREADS: 1 + NUMEXPR_NUM_THREADS: 1 steps: - name: Checkout code uses: actions/checkout@v3 From 78e03c349389205aca222e5d6da20ee23b4ee210 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 29 Jul 2025 11:33:54 +0200 Subject: [PATCH 22/48] replaced default SVD from randomized_svd to scipy svd --- dingo/gw/SVD.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/dingo/gw/SVD.py b/dingo/gw/SVD.py index 7379a8c7f..9852c66e4 100644 --- a/dingo/gw/SVD.py +++ b/dingo/gw/SVD.py @@ -24,7 +24,7 @@ def __init__( data_keys=["V", "s", "mismatches"], ) - def generate_basis(self, training_data: np.ndarray, n: int, method: str = "random"): + def generate_basis(self, training_data: np.ndarray, n: int, method: str = "scipy"): """Generate the SVD basis from training data and store it. The SVD decomposition takes @@ -56,8 +56,17 @@ def generate_basis(self, training_data: np.ndarray, n: int, method: str = "rando # The randomized SVD has complexity O(m n k + k^2 (m + n)), # for a m x n matrix and k is the target rank, here called n # For small k this is much faster than the standard SVD. - U, s, Vh = randomized_svd(training_data, n, random_state=0, - power_iteration_normalizer='QR') + try: + U, s, Vh = randomized_svd(training_data, n, random_state=0, + power_iteration_normalizer='QR') + except ValueError as e: + raise ValueError( + "randomized_svd failed — possibly due to complex-valued input.\n" + "randomized_svd does not support complex arrays in scikit-learn >=1.2.\n" + "To proceed, downgrade scikit-learn to version 1.1.3:\n\n" + " pip install scikit-learn==1.1.3\n\n" + f"Original error: {e}" + ) self.Vh = Vh.astype(np.complex128) # TODO: fix types self.V = self.Vh.T.conj() From 14c279b4e4b849eb53fe63f91bde52f7a1e54fff Mon Sep 17 00:00:00 2001 From: Nihar Gupte <40392612+nihargupte-ph@users.noreply.github.com> Date: Tue, 29 Jul 2025 11:33:15 +0200 Subject: [PATCH 23/48] Update pytest.yml Adding BLAS and MKL threading options so that the scipy/numpy svd function doesn't segfault --- .github/workflows/pytest.yml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index bd64e0682..762903751 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -13,6 +13,10 @@ on: jobs: pytest: runs-on: ubuntu-latest + env: + OPENBLAS_NUM_THREADS: 1 + MKL_NUM_THREADS: 1 + NUMEXPR_NUM_THREADS: 1 steps: - name: Checkout code uses: actions/checkout@v3 From 12a0a1c74dc4769c925c097d42096c62afca6548 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 29 Jul 2025 11:33:54 +0200 Subject: [PATCH 24/48] replaced default SVD from randomized_svd to scipy svd --- dingo/gw/SVD.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/dingo/gw/SVD.py b/dingo/gw/SVD.py index 7379a8c7f..9852c66e4 100644 --- a/dingo/gw/SVD.py +++ b/dingo/gw/SVD.py @@ -24,7 +24,7 @@ def __init__( data_keys=["V", "s", "mismatches"], ) - def generate_basis(self, training_data: np.ndarray, n: int, method: str = "random"): + def generate_basis(self, training_data: np.ndarray, n: int, method: str = "scipy"): """Generate the SVD basis from training data and store it. The SVD decomposition takes @@ -56,8 +56,17 @@ def generate_basis(self, training_data: np.ndarray, n: int, method: str = "rando # The randomized SVD has complexity O(m n k + k^2 (m + n)), # for a m x n matrix and k is the target rank, here called n # For small k this is much faster than the standard SVD. - U, s, Vh = randomized_svd(training_data, n, random_state=0, - power_iteration_normalizer='QR') + try: + U, s, Vh = randomized_svd(training_data, n, random_state=0, + power_iteration_normalizer='QR') + except ValueError as e: + raise ValueError( + "randomized_svd failed — possibly due to complex-valued input.\n" + "randomized_svd does not support complex arrays in scikit-learn >=1.2.\n" + "To proceed, downgrade scikit-learn to version 1.1.3:\n\n" + " pip install scikit-learn==1.1.3\n\n" + f"Original error: {e}" + ) self.Vh = Vh.astype(np.complex128) # TODO: fix types self.V = self.Vh.T.conj() From 6f6451cb68b043f738dddc471ea1f661839a298b Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 29 Jul 2025 17:09:49 +0200 Subject: [PATCH 25/48] SEOBNRv5EHM changes --- .../waveform_generator/waveform_generator.py | 32 ++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index 53b7471be..208b6b53a 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -1017,6 +1017,10 @@ def _convert_parameters( f_min = self.f_start else: f_min = self.domain.f_min + # for SEOBNRv5EHM, the starting frequency must be the same as the reference frequency + if self.approximant_str == "SEOBNRv5EHM": + f_min = self.f_ref + f_min = self.domain.f_min # parameters needed for TD waveforms delta_t = 0.5 / self.domain.f_max @@ -1041,6 +1045,31 @@ def _convert_parameters( "condition": 1, } + + # SEOBNRv5EHM doesn't support setting a reference frequency, it is the + # same as the starting frequency + if self.approximant_str == "SEOBNRv5EHM": + # eccentric parameters + if "log10_eccentricity" in p and "eccentricity" in p: + if (10 ** p["log10_eccentricity"] - p["eccentricity"]) > 1e-4: + raise ValueError( + f"""log10_eccentricity of {p['log10_eccentricity']} and eccentricity + of {p['eccentricity']} are inconsistent, check your input values""" + ) + + if "log10_eccentricity" in p: + eccentricity = np.power(10, p["log10_eccentricity"]) + else: + eccentricity = p.get("eccentricity", 0.0) + longitude_ascending_nodes = p.get("long_asc_nodes", 0.0) + mean_per_ano = p.get("mean_anomaly", 0.0) + + params_gwsignal.update({ + 'eccentricity' : eccentricity * u.dimensionless_unscaled, + 'longAscNodes' : longitude_ascending_nodes * u.rad, + 'meanPerAno' : mean_per_ano * u.rad, + }) + # SEOBNRv5 specific parameters if "postadiabatic" in p: params_gwsignal["postadiabatic"] = p["postadiabatic"] @@ -1051,7 +1080,8 @@ def _convert_parameters( if "lmax_nyquist" in p: params_gwsignal["lmax_nyquist"] = p["lmax_nyquist"] else: - params_gwsignal["lmax_nyquist"] = 2 + if not "ROM" in self.approximant_str: + params_gwsignal["lmax_nyquist"] = 2 if return_target_function: # This is a hack to make compatible with LAL version. Target functions for From 3f65567d0268585ec88540da217ee1fc1ba9c71f Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 3 Nov 2025 14:11:47 +0100 Subject: [PATCH 26/48] fixed merge issues --- .../waveform_generator/waveform_generator.py | 34 ++++++++++++++----- dingo/pipe/data_generation.py | 5 --- dingo/pipe/nodes/plot_node.py | 2 +- dingo/pipe/parser.py | 11 ------ 4 files changed, 26 insertions(+), 26 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index 9678484f7..29adfd873 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -960,6 +960,17 @@ def __init__(self, **kwargs): extra_wf_kwargs["lmax_nyquist"] = kwargs.get("lmax_nyquist", 2) self.extra_wf_kwargs = extra_wf_kwargs + allowed_extra_kwargs = { + "postadiabatic", + "postadiabatic_type", + "lmax_nyquist", + "enable_antisymmetric_modes", + "antisymmetric_modes_hm", + } + extra_wf_kwargs = {k: v for k, v in kwargs.items() if k in allowed_extra_kwargs} + extra_wf_kwargs["lmax_nyquist"] = kwargs.get("lmax_nyquist", 2) + self.extra_wf_kwargs = extra_wf_kwargs + @property def domain(self): if self._use_base_domain: @@ -1061,19 +1072,23 @@ def _convert_parameters( # same as the starting frequency if self.approximant_str == "SEOBNRv5EHM": # eccentric parameters - if "log10_eccentricity" in p and "eccentricity" in p: - if (10 ** p["log10_eccentricity"] - p["eccentricity"]) > 1e-4: + log_ecc = p.get("log10_eccentricity") + ecc = p.get("eccentricity") + + if log_ecc is not None and ecc is not None: + if abs(10**log_ecc - ecc) > 1e-4: raise ValueError( - f"""log10_eccentricity of {p['log10_eccentricity']} and eccentricity - of {p['eccentricity']} are inconsistent, check your input values""" + f"log10_eccentricity={log_ecc} and eccentricity={ecc} are inconsistent." ) - if "log10_eccentricity" in p: - eccentricity = np.power(10, p["log10_eccentricity"]) - else: - eccentricity = p.get("eccentricity", 0.0) + eccentricity = 10**log_ecc if log_ecc is not None else p.get("eccentricity", 0.0) longitude_ascending_nodes = p.get("long_asc_nodes", 0.0) - mean_per_ano = p.get("mean_anomaly", 0.0) + if "relativistic_anomaly" in p and "mean_anomaly" in p: + raise ValueError( + "Cannot specify both relativistic_anomaly and mean_anomaly." + ) + else: + mean_per_ano = p.get("mean_anomaly", p.get("relativistic_anomaly", 0.0)) params_gwsignal.update({ 'eccentricity' : eccentricity * u.dimensionless_unscaled, @@ -1388,6 +1403,7 @@ def generate_TD_modes_L0(self, parameters): """ # TD approximants that are implemented in L0 frame. Currently tested for: # 52: SEOBNRv4PHM + # SEOBNRV5EHM parameters_gwsignal = self._convert_parameters( {**parameters, "f_ref": self.f_ref} diff --git a/dingo/pipe/data_generation.py b/dingo/pipe/data_generation.py index 786026ac3..cf2a1f05d 100644 --- a/dingo/pipe/data_generation.py +++ b/dingo/pipe/data_generation.py @@ -23,7 +23,6 @@ from dingo.gw.domains import UniformFrequencyDomain, build_domain_from_model_metadata from dingo.gw.injection import Injection from dingo.pipe.parser import create_parser -from dingo.pipe.main import fill_in_arguments_from_model logger.name = "dingo_pipe" @@ -522,10 +521,6 @@ def _get_priors(self, add_time=True): priors = super()._get_priors(add_time=add_time) priors.update(PriorDict(self.prior_dict_updates)) return priors -<<<<<<< HEAD -======= - ->>>>>>> main def create_generation_parser(): """Data generation parser creation""" diff --git a/dingo/pipe/nodes/plot_node.py b/dingo/pipe/nodes/plot_node.py index 9ed4f1b69..255b6e844 100644 --- a/dingo/pipe/nodes/plot_node.py +++ b/dingo/pipe/nodes/plot_node.py @@ -94,4 +94,4 @@ def executable(self): @property def request_memory(self): - return "16 GB" + return "16 GB" \ No newline at end of file diff --git a/dingo/pipe/parser.py b/dingo/pipe/parser.py index 82cabfc82..ce22535b8 100644 --- a/dingo/pipe/parser.py +++ b/dingo/pipe/parser.py @@ -1337,16 +1337,6 @@ def create_parser(top_level=True, usage=None): ), ) # waveform_parser.add( -<<<<<<< HEAD - # "--waveform-arguments-dict", - # default=None, - # type=nonestr, - # help=( - # "A dictionary of arbitrary additional waveform-arguments to pass" - # " to the bilby waveform generator's `waveform_arguments`. Only used " - # "for injections" - # ), -======= # "--waveform-arguments-dict", # default=None, # type=nonestr, @@ -1355,7 +1345,6 @@ def create_parser(top_level=True, usage=None): # " to the bilby waveform generator's `waveform_arguments`. Only used " # "for injections" # ), ->>>>>>> main # ) # waveform_parser.add( # "--mode-array", From 64b38db5425a5d458770e58d7a7811cfee727cd5 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 3 Nov 2025 14:23:07 +0100 Subject: [PATCH 27/48] removing data_generationextra lines --- dingo/pipe/data_generation.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/dingo/pipe/data_generation.py b/dingo/pipe/data_generation.py index cf2a1f05d..338c74682 100644 --- a/dingo/pipe/data_generation.py +++ b/dingo/pipe/data_generation.py @@ -517,11 +517,13 @@ def prior_dict_updates(self, prior_dict_updates): self._convert_prior_dict_key(key): val for key, val in prior_dict_updates.items() } + def _get_priors(self, add_time=True): priors = super()._get_priors(add_time=add_time) priors.update(PriorDict(self.prior_dict_updates)) return priors + def create_generation_parser(): """Data generation parser creation""" return create_parser(top_level=False) From bd2a32301c4a9e1809f133e1d1878989085e5814 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 3 Nov 2025 14:27:07 +0100 Subject: [PATCH 28/48] clarifying docstrings --- dingo/gw/waveform_generator/waveform_generator.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index 29adfd873..e45dddb2c 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -1040,6 +1040,7 @@ def _convert_parameters( else: f_min = self.domain.f_min # for SEOBNRv5EHM, the starting frequency must be the same as the reference frequency + # note this is the orbit averaged reference frequency if self.approximant_str == "SEOBNRv5EHM": f_min = self.f_ref f_min = self.domain.f_min @@ -1068,15 +1069,13 @@ def _convert_parameters( } - # SEOBNRv5EHM doesn't support setting a reference frequency, it is the - # same as the starting frequency if self.approximant_str == "SEOBNRv5EHM": # eccentric parameters log_ecc = p.get("log10_eccentricity") ecc = p.get("eccentricity") if log_ecc is not None and ecc is not None: - if abs(10**log_ecc - ecc) > 1e-4: + if abs(10**log_ecc - ecc) > 1e-6: raise ValueError( f"log10_eccentricity={log_ecc} and eccentricity={ecc} are inconsistent." ) @@ -1289,6 +1288,7 @@ def generate_hplus_hcross_m( # Step 2: Transform modes to target domain. hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain) else: + # SEOBNRv5EHM uses this conditioning routine # assert LS.SimInspiralImplementedTDApproximants(self.approximant) # Step 1: generate waveform modes in L0 frame in native domain of # approximant (here: TD) From 96031923668e9c2c8a15012ded30b74245921ecd Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 3 Nov 2025 15:13:24 +0100 Subject: [PATCH 29/48] added tests for decimation and mode resummation for v5EHM --- tests/gw/waveform_generator/test_wfg_m.py | 33 +++++++-- tests/gw/waveform_generator/test_wfg_mfd.py | 74 +++++++++++++++------ 2 files changed, 80 insertions(+), 27 deletions(-) diff --git a/tests/gw/waveform_generator/test_wfg_m.py b/tests/gw/waveform_generator/test_wfg_m.py index cc8769939..82730c55b 100644 --- a/tests/gw/waveform_generator/test_wfg_m.py +++ b/tests/gw/waveform_generator/test_wfg_m.py @@ -37,14 +37,15 @@ def uniform_fd_domain(): return domain -@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM"]) +@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]) def approximant(request): return request.param @pytest.fixture def intrinsic_prior(approximant): - if "PHM" in approximant: + if approximant in ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM"]: + # quasi-circular precessing-spin intrinsic_dict = { "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", @@ -61,8 +62,21 @@ def intrinsic_prior(approximant): "phi_jl": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', "geocent_time": 0.0, } - else: - # Aligned spins + elif approximant in ["SEOBNRv5HM", "SEOBNRv5EHM"]: + # quasi-circular aligned-spin + intrinsic_dict = { + "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", + "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", + "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)", + "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)", + "luminosity_distance": 1000.0, + "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", + "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', + "chi_1": 'bilby.gw.prior.AlignedSpin(name="chi_1", a_prior=Uniform(minimum=0, maximum=0.99))', + "chi_2": 'bilby.gw.prior.AlignedSpin(name="chi_2", a_prior=Uniform(minimum=0, maximum=0.99))', + "geocent_time": 0.0, + } + elif approximant in ["SEOBNRv5EHM"]: intrinsic_dict = { "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", @@ -73,15 +87,20 @@ def intrinsic_prior(approximant): "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', "chi_1": 'bilby.gw.prior.AlignedSpin(name="chi_1", a_prior=Uniform(minimum=0, maximum=0.99))', "chi_2": 'bilby.gw.prior.AlignedSpin(name="chi_2", a_prior=Uniform(minimum=0, maximum=0.99))', + "eccentricity": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.3)", + "relativistic_anomaly": "bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi)", "geocent_time": 0.0, } + else: + raise ValueError(f"Unimplemented approximant {approximant}") + prior = build_prior_with_defaults(intrinsic_dict) return prior @pytest.fixture def wfg(uniform_fd_domain, approximant): - if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM"]: + if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]: wfg_class = NewInterfaceWaveformGenerator else: wfg_class = WaveformGenerator @@ -115,7 +134,7 @@ def tolerances(approximant): # get should not have a big effect in practice. return 2e-2, 1e-5 - elif approximant == "SEOBNRv4PHM": + elif approximant in ["SEOBNRv4PHM", "SEOBNRv5EHM"]: # The mismatches are typically be of order 1e-5. This is exclusively due to # different tapering. The reference polarizations are tapered and FFTed on the # level of polarizations, while for generate_hplus_hcross_m, the tapering and FFT @@ -136,7 +155,7 @@ def tolerances(approximant): try: import pyseobnr - approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM"] + approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"] except ImportError: approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM"] diff --git a/tests/gw/waveform_generator/test_wfg_mfd.py b/tests/gw/waveform_generator/test_wfg_mfd.py index 297df73bd..63da7533a 100644 --- a/tests/gw/waveform_generator/test_wfg_mfd.py +++ b/tests/gw/waveform_generator/test_wfg_mfd.py @@ -24,36 +24,70 @@ def mfd(): return domain -@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM"]) +@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]) def approximant(request): return request.param @pytest.fixture def intrinsic_prior(approximant): - intrinsic_dict = { - "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", - "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", - "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)", - "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)", - "luminosity_distance": 1000.0, - "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", - "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', - "a_1": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)", - "a_2": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)", - "tilt_1": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", - "tilt_2": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", - "phi_12": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', - "phi_jl": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', - "geocent_time": "bilby.core.prior.Uniform(minimum=-0.1, maximum=0.1)", - } + if approximant in ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM"]: + # quasi-circular precessing-spin + intrinsic_dict = { + "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", + "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", + "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)", + "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)", + "luminosity_distance": 1000.0, + "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", + "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', + "a_1": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)", + "a_2": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.99)", + "tilt_1": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", + "tilt_2": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", + "phi_12": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', + "phi_jl": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', + "geocent_time": 0.0, + } + elif approximant in ["SEOBNRv5HM", "SEOBNRv5EHM"]: + # quasi-circular aligned-spin + intrinsic_dict = { + "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", + "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", + "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)", + "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)", + "luminosity_distance": 1000.0, + "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", + "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', + "chi_1": 'bilby.gw.prior.AlignedSpin(name="chi_1", a_prior=Uniform(minimum=0, maximum=0.99))', + "chi_2": 'bilby.gw.prior.AlignedSpin(name="chi_2", a_prior=Uniform(minimum=0, maximum=0.99))', + "geocent_time": 0.0, + } + elif approximant in ["SEOBNRv5EHM"]: + intrinsic_dict = { + "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", + "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", + "mass_ratio": "bilby.gw.prior.UniformInComponentsMassRatio(minimum=0.125, maximum=1.0)", + "chirp_mass": "bilby.gw.prior.UniformInComponentsChirpMass(minimum=25.0, maximum=100.0)", + "luminosity_distance": 1000.0, + "theta_jn": "bilby.core.prior.Sine(minimum=0.0, maximum=np.pi)", + "phase": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', + "chi_1": 'bilby.gw.prior.AlignedSpin(name="chi_1", a_prior=Uniform(minimum=0, maximum=0.99))', + "chi_2": 'bilby.gw.prior.AlignedSpin(name="chi_2", a_prior=Uniform(minimum=0, maximum=0.99))', + "eccentricity": "bilby.core.prior.Uniform(minimum=0.0, maximum=0.3)", + "relativistic_anomaly": "bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi)", + "geocent_time": 0.0, + } + else: + raise ValueError(f"Unimplemented approximant {approximant}") + prior = build_prior_with_defaults(intrinsic_dict) return prior @pytest.fixture def wfg_mfd(mfd, approximant): - if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM"]: + if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]: wfg_class = NewInterfaceWaveformGenerator else: wfg_class = WaveformGenerator @@ -68,7 +102,7 @@ def wfg_mfd(mfd, approximant): @pytest.fixture def wfg_ufd(mfd, approximant): - if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM"]: + if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]: wfg_class = NewInterfaceWaveformGenerator else: wfg_class = WaveformGenerator @@ -102,7 +136,7 @@ def tolerances(approximant): try: import pyseobnr - approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM"] + approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"] except ImportError: approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM"] From 91e615e02fad87ce2d4568d38dcc9c799d2ab616 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 3 Nov 2025 15:17:47 +0100 Subject: [PATCH 30/48] removed extra allowed kwargs --- dingo/gw/waveform_generator/waveform_generator.py | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index e45dddb2c..0680632f8 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -960,17 +960,6 @@ def __init__(self, **kwargs): extra_wf_kwargs["lmax_nyquist"] = kwargs.get("lmax_nyquist", 2) self.extra_wf_kwargs = extra_wf_kwargs - allowed_extra_kwargs = { - "postadiabatic", - "postadiabatic_type", - "lmax_nyquist", - "enable_antisymmetric_modes", - "antisymmetric_modes_hm", - } - extra_wf_kwargs = {k: v for k, v in kwargs.items() if k in allowed_extra_kwargs} - extra_wf_kwargs["lmax_nyquist"] = kwargs.get("lmax_nyquist", 2) - self.extra_wf_kwargs = extra_wf_kwargs - @property def domain(self): if self._use_base_domain: From 1ae0d76a35dd453250e7a220cdd3fd5a6bd6f8a0 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Thu, 13 Nov 2025 11:56:01 +0100 Subject: [PATCH 31/48] teobresums-dali implementation --- .../waveform_generator/waveform_generator.py | 90 ++++++++++++------- dingo/gw/waveform_generator/wfg_utils.py | 83 +++++++++++++++++ pyproject.toml | 4 + tests/gw/waveform_generator/test_wfg_m.py | 20 +++-- tests/gw/waveform_generator/test_wfg_mfd.py | 16 ++-- 5 files changed, 164 insertions(+), 49 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index 0680632f8..b12b92878 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -4,6 +4,7 @@ import numpy as np import astropy.units as u +import astropy.constants as ac from typing import Dict, List, Tuple, Union, Callable from numbers import Number import warnings @@ -727,7 +728,7 @@ def generate_hplus_hcross_m( if LS.SimInspiralImplementedFDApproximants(self.approximant): # Step 1: generate waveform modes in L0 frame in native domain of # approximant (here: FD) - hlm_fd, iota = self.generate_FD_modes_LO(parameters) + hlm_fd, iota = self.generate_FD_modes_L0(parameters) # Step 2: Transform modes to target domain. # Not required here, as approximant domain and target domain are both FD. @@ -769,7 +770,7 @@ def generate_hplus_hcross_m( self._use_base_domain = True self._domain_transform = DecimateAll(self._domain) - hlm_fd, iota = self.generate_FD_modes_LO(parameters) + hlm_fd, iota = self.generate_FD_modes_L0(parameters) pol_m = wfg_utils.get_polarizations_from_fd_modes_m( hlm_fd, iota, parameters["phase"] ) @@ -800,7 +801,7 @@ def generate_hplus_hcross_m( else: return pol_m - def generate_FD_modes_LO(self, parameters): + def generate_FD_modes_L0(self, parameters): """ Generate FD modes in the L0 frame. @@ -957,7 +958,8 @@ def __init__(self, **kwargs): "antisymmetric_modes_hm", } extra_wf_kwargs = {k: v for k, v in kwargs.items() if k in allowed_extra_kwargs} - extra_wf_kwargs["lmax_nyquist"] = kwargs.get("lmax_nyquist", 2) + if "SEOB" in self.approximant_str and not "ROM" in self.approximant_str: + extra_wf_kwargs["lmax_nyquist"] = kwargs.get("lmax_nyquist", 2) self.extra_wf_kwargs = extra_wf_kwargs @property @@ -1058,7 +1060,7 @@ def _convert_parameters( } - if self.approximant_str == "SEOBNRv5EHM": + if self.approximant_str in ["SEOBNRv5EHM", "TEOBResumSDALI"]: # eccentric parameters log_ecc = p.get("log10_eccentricity") ecc = p.get("eccentricity") @@ -1084,18 +1086,7 @@ def _convert_parameters( 'meanPerAno' : mean_per_ano * u.rad, }) - # SEOBNRv5 specific parameters - if "postadiabatic" in p: - params_gwsignal["postadiabatic"] = p["postadiabatic"] - if "postadiabatic_type" in p: - params_gwsignal["postadiabatic_type"] = p["postadiabatic_type"] - - if "lmax_nyquist" in p: - params_gwsignal["lmax_nyquist"] = p["lmax_nyquist"] - else: - if not "ROM" in self.approximant_str: - params_gwsignal["lmax_nyquist"] = 2 params_gwsignal.update(self.extra_wf_kwargs) if return_target_function: @@ -1256,18 +1247,28 @@ def generate_hplus_hcross_m( generator = new_interface_get_waveform_generator(self.approximant_str) if isinstance(self.domain, UniformFrequencyDomain): # Generate FD modes in for frequencies [-f_max, ..., 0, ..., f_max]. - if generator.domain == "freq": + if self.approximant_str in ["SEOBNRv5EHM", "TEOBResumSDALI"]: + # SEOBNRv5EHM and TEOBResumSDALI use these conditioning routines + # assert LS.SimInspiralImplementedTDApproximants(self.approximant) # Step 1: generate waveform modes in L0 frame in native domain of - # approximant (here: FD) - hlm_fd, iota = self.generate_FD_modes_LO(parameters) + # approximant (here: TD) + hlm_td, iota = self.generate_TD_modes_L0(parameters) # Step 2: Transform modes to target domain. - # Not required here, as approximant domain and target domain are both FD. + # This requires tapering of TD modes, and FFT to transform to FD. + # TEOB and SEOB have slightly different conditioning routine implementations + if self.approximant_str in ["TEOBResumSDALI"]: + wfg_utils.taper_td_modes_in_place_gwsignal(hlm_td) + hlm_fd = wfg_utils.td_modes_to_fd_modes_gwsignal(hlm_td, self.domain) + else: + wfg_utils.taper_td_modes_in_place(hlm_td) + hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain) + + # TODO, why are these different? Shouldn't the gwsignal and non-gwsignal be the + # exact same implementation? - elif ( - self.approximant_str == "SEOBNRv5PHM" - or self.approximant_str == "SEOBNRv5HM" - ): + + elif self.approximant_str in ["SEOBNRv5PHM", "SEOBNRv5HM"]: # Step 1: generate waveform modes in L0 frame in native domain of # approximant (here: TD), applying standard conditioning hlm_td, iota = self.generate_TD_modes_L0_conditioned_extra_time( @@ -1276,17 +1277,18 @@ def generate_hplus_hcross_m( # Step 2: Transform modes to target domain. hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain) - else: - # SEOBNRv5EHM uses this conditioning routine - # assert LS.SimInspiralImplementedTDApproximants(self.approximant) + + elif generator.domain == "freq": # Step 1: generate waveform modes in L0 frame in native domain of - # approximant (here: TD) - hlm_td, iota = self.generate_TD_modes_L0(parameters) + # approximant (here: FD) + hlm_fd, iota = self.generate_FD_modes_L0(parameters) # Step 2: Transform modes to target domain. - # This requires tapering of TD modes, and FFT to transform to FD. - wfg_utils.taper_td_modes_in_place(hlm_td) - hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain) + # Not required here, as approximant domain and target domain are both FD. + + else: + raise NotImplementedError() + # Step 3: Separate negative and positive frequency parts of the modes, # and add contributions according to their transformation behavior under @@ -1305,7 +1307,7 @@ def generate_hplus_hcross_m( else: return pol_m - def generate_FD_modes_LO(self, parameters): # Pending to adapt + def generate_FD_modes_L0(self, parameters): # Pending to adapt """ Generate FD modes in the L0 frame. @@ -1402,6 +1404,25 @@ def generate_TD_modes_L0(self, parameters): hlm_td = gws_wfm.GenerateTDModes(parameters_gwsignal, generator) hlms_lal = {} + # TEOBREsumS returns the modes unscaled in the distance version when calling GenerateTDModes + # therefore, we rescale them here + if "TEOB" in self.approximant_str: + m1, m2 = parameters_gwsignal["mass1"], parameters_gwsignal["mass2"] + nu = m1 * m2 / (m1 + m2) ** 2 + distance_rescaling = -( + ( + nu + * (parameters_gwsignal["mass1"] + parameters_gwsignal["mass2"]) + / parameters_gwsignal["distance"] + * ac.G + / ac.c ** 2 + ) + .to(u.dimensionless_unscaled) + .value + ) + else: + distance_rescaling = 1.0 + for key, value in hlm_td.items(): if type(key) != str: hlm_lal = lal.CreateCOMPLEX16TimeSeries( @@ -1412,9 +1433,10 @@ def generate_TD_modes_L0(self, parameters): lal.DimensionlessUnit, len(value), ) - hlm_lal.data.data = value.value + hlm_lal.data.data = value.value * distance_rescaling hlms_lal[key] = hlm_lal + return hlms_lal, parameters_gwsignal["inclination"].value def generate_TD_modes_L0_conditioned_extra_time(self, parameters): diff --git a/dingo/gw/waveform_generator/wfg_utils.py b/dingo/gw/waveform_generator/wfg_utils.py index 5a40ae381..638e380a2 100644 --- a/dingo/gw/waveform_generator/wfg_utils.py +++ b/dingo/gw/waveform_generator/wfg_utils.py @@ -1,6 +1,8 @@ import numpy as np import lal import lalsimulation as LS +from lalsimulation.gwsignal.core import conditioning_subroutines as cond +import astropy.units as u def linked_list_modes_to_dict_modes(hlm_ll): @@ -69,6 +71,63 @@ def taper_td_modes_in_place(hlm_td, tapering_flag: int = 1): window = get_tapering_window_for_complex_time_series(h, tapering_flag) h.data.data *= window +def taper_td_modes_in_place_gwsignal(hlm_td, tapering_flag: str = 'start'): + """ + Taper the time domain modes in place using gwsignal conditioning routines. + There are apparently slight differences in the tapering functions + between these two and since some wfs (EG TEOB) uses the gwsignal + routine, we use the gwsignal routine sometimes. + + Parameters + ---------- + hlm_td: dict + Dictionary with (l,m) keys and the complex lal time series objects for the + corresponding modes. + + """ + lalseries_to_gwpy_timeseries_in_place(hlm_td) + for mode in hlm_td.keys(): + hlm_td[mode] = cond.taper_gwpy_timeseries(hlm_td[mode], tapering_flag) + +def lalseries_to_gwpy_timeseries_in_place(hlm_td): + """ + Convert lal time series in place to gwpy time series. + + Parameters + ---------- + hlm_td: dict + Dictionary with (l,m) keys and the complex lal time series objects for the + corresponding modes. + + """ + from gwpy.timeseries import TimeSeries + for mode in hlm_td.keys(): + times = ( + hlm_td[mode].epoch.gpsSeconds + + hlm_td[mode].epoch.gpsNanoSeconds * 1e-9 + + np.arange(hlm_td[mode].data.length) * hlm_td[mode].deltaT + ) + + hlm_td[mode] = TimeSeries( + data=hlm_td[mode].data.data, + times=times, + name=f"h_{mode[0]}_{mode[1]}", + unit=u.dimensionless_unscaled, + ) + +def gwpyseries_to_lalseries_in_place(hlm_td): + """ + Convert gwpy time series in place to lal time series. + + Parameters + ---------- + hlm_td: dict + Dictionary with (l,m) keys and the complex gwpy time series objects for the + corresponding modes. + + """ + for mode in hlm_td.keys(): + hlm_td[mode] = hlm_td[mode].to_lal() def td_modes_to_fd_modes(hlm_td, domain): """ @@ -137,6 +196,30 @@ def td_modes_to_fd_modes(hlm_td, domain): return hlm_fd +def td_modes_to_fd_modes_gwsignal(hlm_td, domain): + """ + TODO + """ + hlm_fd = {} + + delta_f = domain.delta_f + delta_t = 0.5 / domain.f_max + f_nyquist = domain.f_max # use f_max as f_nyquist + chirplen = int(2 * f_nyquist / delta_f) + + for mode in hlm_td.keys(): + hlm_td_cond = cond.resize_gwpy_timeseries(hlm_td[mode], len(hlm_td[mode]) - chirplen, chirplen) + hlm_fd_cond = hlm_td_cond.fft() + hlm_fd_cond.epoch = hlm_td_cond.t0 + hlm_fd_cond = (hlm_fd_cond / (2 * hlm_fd_cond.df)).value + + # setting up reflection around boundary to create two-sided FFT + hlm_fd_cond = np.concatenate((hlm_fd_cond[::-1], hlm_fd_cond[1:]), axis=0) + + hlm_fd_cond[-1] = hlm_fd_cond[0] + hlm_fd[mode] = hlm_fd_cond + + return hlm_fd def get_polarizations_from_fd_modes_m(hlm_fd, iota, phase): pol_m = {} diff --git a/pyproject.toml b/pyproject.toml index 7d166faf9..ef68864b4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -85,6 +85,10 @@ pyseobnr = [ "pyseobnr", ] +teobresums = [ + "teobresums", +] + [project.scripts] dingo_append_training_stage = "dingo.gw.training:append_stage" dingo_build_svd = "dingo.gw.dataset.utils:build_svd_cli" diff --git a/tests/gw/waveform_generator/test_wfg_m.py b/tests/gw/waveform_generator/test_wfg_m.py index 82730c55b..a909b2fec 100644 --- a/tests/gw/waveform_generator/test_wfg_m.py +++ b/tests/gw/waveform_generator/test_wfg_m.py @@ -37,7 +37,7 @@ def uniform_fd_domain(): return domain -@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]) +@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"]) def approximant(request): return request.param @@ -62,7 +62,7 @@ def intrinsic_prior(approximant): "phi_jl": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', "geocent_time": 0.0, } - elif approximant in ["SEOBNRv5HM", "SEOBNRv5EHM"]: + elif approximant in ["SEOBNRv5HM"]: # quasi-circular aligned-spin intrinsic_dict = { "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", @@ -76,7 +76,7 @@ def intrinsic_prior(approximant): "chi_2": 'bilby.gw.prior.AlignedSpin(name="chi_2", a_prior=Uniform(minimum=0, maximum=0.99))', "geocent_time": 0.0, } - elif approximant in ["SEOBNRv5EHM"]: + elif approximant in ["SEOBNRv5EHM", "TEOBResumSDALI"]: intrinsic_dict = { "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", @@ -100,7 +100,7 @@ def intrinsic_prior(approximant): @pytest.fixture def wfg(uniform_fd_domain, approximant): - if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]: + if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"]: wfg_class = NewInterfaceWaveformGenerator else: wfg_class = WaveformGenerator @@ -134,7 +134,7 @@ def tolerances(approximant): # get should not have a big effect in practice. return 2e-2, 1e-5 - elif approximant in ["SEOBNRv4PHM", "SEOBNRv5EHM"]: + elif approximant in ["SEOBNRv4PHM", "SEOBNRv5EHM", "TEOBResumSDALI"]: # The mismatches are typically be of order 1e-5. This is exclusively due to # different tapering. The reference polarizations are tapered and FFTed on the # level of polarizations, while for generate_hplus_hcross_m, the tapering and FFT @@ -153,10 +153,13 @@ def tolerances(approximant): # Uncomment to test only one approximant. try: - import pyseobnr + # import pyseobnr + import EOBRun_module - approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"] + # approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"] + approximant_list = ["TEOBResumSDALI"] except ImportError: + raise ValueError("TMP") approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM"] @@ -183,7 +186,7 @@ def test_generate_hplus_hcross_m(intrinsic_prior, wfg, num_evaluations, toleranc ] ) - debug = False + debug = True if debug: maxval = max(mismatches[-1]) idx = mismatches[-1].index(maxval) @@ -198,6 +201,7 @@ def test_generate_hplus_hcross_m(intrinsic_prior, wfg, num_evaluations, toleranc plt.xscale("log") plt.xlim((5, 128)) plt.title(f"{p}, mismatch={maxval}") + plt.savefig("/work/nihargupte/src/dingo/tests/gw/mm.png") plt.show() mismatches = np.array(mismatches) diff --git a/tests/gw/waveform_generator/test_wfg_mfd.py b/tests/gw/waveform_generator/test_wfg_mfd.py index 63da7533a..badd9ac59 100644 --- a/tests/gw/waveform_generator/test_wfg_mfd.py +++ b/tests/gw/waveform_generator/test_wfg_mfd.py @@ -24,7 +24,7 @@ def mfd(): return domain -@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]) +@pytest.fixture(params=["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"]) def approximant(request): return request.param @@ -49,7 +49,7 @@ def intrinsic_prior(approximant): "phi_jl": 'bilby.core.prior.Uniform(minimum=0.0, maximum=2*np.pi, boundary="periodic")', "geocent_time": 0.0, } - elif approximant in ["SEOBNRv5HM", "SEOBNRv5EHM"]: + elif approximant in ["SEOBNRv5HM"]: # quasi-circular aligned-spin intrinsic_dict = { "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", @@ -63,7 +63,7 @@ def intrinsic_prior(approximant): "chi_2": 'bilby.gw.prior.AlignedSpin(name="chi_2", a_prior=Uniform(minimum=0, maximum=0.99))', "geocent_time": 0.0, } - elif approximant in ["SEOBNRv5EHM"]: + elif approximant in ["SEOBNRv5EHM", "TEOBResumSDALI"]: intrinsic_dict = { "mass_1": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", "mass_2": "bilby.core.prior.Constraint(minimum=10.0, maximum=80.0)", @@ -87,7 +87,7 @@ def intrinsic_prior(approximant): @pytest.fixture def wfg_mfd(mfd, approximant): - if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]: + if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"]: wfg_class = NewInterfaceWaveformGenerator else: wfg_class = WaveformGenerator @@ -102,7 +102,7 @@ def wfg_mfd(mfd, approximant): @pytest.fixture def wfg_ufd(mfd, approximant): - if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"]: + if approximant in ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"]: wfg_class = NewInterfaceWaveformGenerator else: wfg_class = WaveformGenerator @@ -134,9 +134,11 @@ def tolerances(approximant): # Uncomment to test only one approximant. try: - import pyseobnr + # import pyseobnr + import EOBRun_module - approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"] + # approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"] + approximant_list = ["TEOBResumSDALI"] except ImportError: approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM"] From 01040c76ea09b87c3356e7d4a524c7de593e9d67 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Thu, 13 Nov 2025 13:28:23 +0100 Subject: [PATCH 32/48] updated f_min typo! --- dingo/gw/waveform_generator/waveform_generator.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index b12b92878..aafea2001 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -1034,7 +1034,6 @@ def _convert_parameters( # note this is the orbit averaged reference frequency if self.approximant_str == "SEOBNRv5EHM": f_min = self.f_ref - f_min = self.domain.f_min # parameters needed for TD waveforms delta_t = 0.5 / self.domain.f_max @@ -1086,7 +1085,6 @@ def _convert_parameters( 'meanPerAno' : mean_per_ano * u.rad, }) - params_gwsignal.update(self.extra_wf_kwargs) if return_target_function: From ba0988822fbc64313ceec9ca1e3e2a183e0f29d9 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 18 Nov 2025 10:33:26 +0100 Subject: [PATCH 33/48] 2 sided fft --- dingo/gw/waveform_generator/wfg_utils.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/dingo/gw/waveform_generator/wfg_utils.py b/dingo/gw/waveform_generator/wfg_utils.py index 638e380a2..9534db7e4 100644 --- a/dingo/gw/waveform_generator/wfg_utils.py +++ b/dingo/gw/waveform_generator/wfg_utils.py @@ -214,11 +214,13 @@ def td_modes_to_fd_modes_gwsignal(hlm_td, domain): hlm_fd_cond = (hlm_fd_cond / (2 * hlm_fd_cond.df)).value # setting up reflection around boundary to create two-sided FFT - hlm_fd_cond = np.concatenate((hlm_fd_cond[::-1], hlm_fd_cond[1:]), axis=0) + hlm_fd_cond = np.concatenate((hlm_fd_cond[::-1], np.conj(hlm_fd_cond[1:])), axis=0) hlm_fd_cond[-1] = hlm_fd_cond[0] hlm_fd[mode] = hlm_fd_cond + hlm_fd[mode] = hlm_fd[mode] / (2 * hlm_fd_cond.df) + return hlm_fd def get_polarizations_from_fd_modes_m(hlm_fd, iota, phase): From 10db6cd0247736c7564553294f6f729283b8fc18 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Thu, 27 Nov 2025 18:14:54 +0100 Subject: [PATCH 34/48] renamed mean_anomaly to mean_per_ano --- dingo/asimov/dingo.ini | 6 +++--- dingo/gw/waveform_generator/waveform_generator.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/dingo/asimov/dingo.ini b/dingo/asimov/dingo.ini index 2d1e6f543..cd5caea33 100644 --- a/dingo/asimov/dingo.ini +++ b/dingo/asimov/dingo.ini @@ -202,9 +202,9 @@ prior-dict-updates = { {%- assign p = priors['eccentricity'] %} eccentricity = {{ p['type'] }}(name='eccentricity', minimum={{ p['minimum'] }}, maximum={{ p['maximum'] }}), {%- endif %} -{%- if priors.keys() contains "mean_anomaly" %} -{%- assign p = priors['mean_anomaly'] %} - mean_anomaly = {{ p['type'] }}(name='mean_anomaly', minimum={{ p['minimum'] }}, maximum={{ p['maximum'] }}, boundary={{ p['boundary'] | default: 'periodic' }} ), +{%- if priors.keys() contains "mean_per_ano" %} +{%- assign p = priors['mean_per_ano'] %} + mean_per_ano = {{ p['type'] }}(name='mean_per_ano', minimum={{ p['minimum'] }}, maximum={{ p['maximum'] }}, boundary={{ p['boundary'] | default: 'periodic' }} ), {%- endif %} {%- if priors.keys() contains "luminosity distance" %} {%- assign p = priors['luminosity distance'] %} diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index aafea2001..a7e130161 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -1072,12 +1072,12 @@ def _convert_parameters( eccentricity = 10**log_ecc if log_ecc is not None else p.get("eccentricity", 0.0) longitude_ascending_nodes = p.get("long_asc_nodes", 0.0) - if "relativistic_anomaly" in p and "mean_anomaly" in p: + if "relativistic_anomaly" in p and "mean_per_ano" in p: raise ValueError( - "Cannot specify both relativistic_anomaly and mean_anomaly." + "Cannot specify both relativistic_anomaly and mean_per_ano." ) else: - mean_per_ano = p.get("mean_anomaly", p.get("relativistic_anomaly", 0.0)) + mean_per_ano = p.get("mean_per_ano", p.get("relativistic_anomaly", 0.0)) params_gwsignal.update({ 'eccentricity' : eccentricity * u.dimensionless_unscaled, From e6543456954c6138cbb01d303320868803e0c65c Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 20 Jan 2026 09:48:41 +0100 Subject: [PATCH 35/48] teob support with synthetic phase --- dingo/gw/waveform_generator/waveform_generator.py | 5 +---- dingo/gw/waveform_generator/wfg_utils.py | 10 +++++----- 2 files changed, 6 insertions(+), 9 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index a7e130161..1dfb8bcaa 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -1255,6 +1255,7 @@ def generate_hplus_hcross_m( # Step 2: Transform modes to target domain. # This requires tapering of TD modes, and FFT to transform to FD. # TEOB and SEOB have slightly different conditioning routine implementations + # TEOB uses gwsignal routines SEOB uses lal routines if self.approximant_str in ["TEOBResumSDALI"]: wfg_utils.taper_td_modes_in_place_gwsignal(hlm_td) hlm_fd = wfg_utils.td_modes_to_fd_modes_gwsignal(hlm_td, self.domain) @@ -1262,10 +1263,6 @@ def generate_hplus_hcross_m( wfg_utils.taper_td_modes_in_place(hlm_td) hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain) - # TODO, why are these different? Shouldn't the gwsignal and non-gwsignal be the - # exact same implementation? - - elif self.approximant_str in ["SEOBNRv5PHM", "SEOBNRv5HM"]: # Step 1: generate waveform modes in L0 frame in native domain of # approximant (here: TD), applying standard conditioning diff --git a/dingo/gw/waveform_generator/wfg_utils.py b/dingo/gw/waveform_generator/wfg_utils.py index 9534db7e4..fab45e64b 100644 --- a/dingo/gw/waveform_generator/wfg_utils.py +++ b/dingo/gw/waveform_generator/wfg_utils.py @@ -211,16 +211,16 @@ def td_modes_to_fd_modes_gwsignal(hlm_td, domain): hlm_td_cond = cond.resize_gwpy_timeseries(hlm_td[mode], len(hlm_td[mode]) - chirplen, chirplen) hlm_fd_cond = hlm_td_cond.fft() hlm_fd_cond.epoch = hlm_td_cond.t0 - hlm_fd_cond = (hlm_fd_cond / (2 * hlm_fd_cond.df)).value + delta_f = hlm_fd_cond.df + hlm_fd_cond = (hlm_fd_cond / (2 * delta_f)).value # setting up reflection around boundary to create two-sided FFT - hlm_fd_cond = np.concatenate((hlm_fd_cond[::-1], np.conj(hlm_fd_cond[1:])), axis=0) + # hlm_fd_cond = np.concatenate((hlm_fd_cond[::-1], np.conj(hlm_fd_cond[1:])), axis=0) - hlm_fd_cond[-1] = hlm_fd_cond[0] + # hlm_fd_cond[-1] = hlm_fd_cond[0] + # hlm_fd[mode] = hlm_fd[mode] / (2 * delta_f) hlm_fd[mode] = hlm_fd_cond - hlm_fd[mode] = hlm_fd[mode] / (2 * hlm_fd_cond.df) - return hlm_fd def get_polarizations_from_fd_modes_m(hlm_fd, iota, phase): From 966861cc10a50651b9f780a2f86368f635be2026 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Fri, 23 Jan 2026 12:37:14 +0100 Subject: [PATCH 36/48] removed special TEOB tapering --- .../waveform_generator/waveform_generator.py | 19 +++++++------------ dingo/gw/waveform_generator/wfg_utils.py | 3 ++- 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index 1dfb8bcaa..afb185572 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -1254,14 +1254,8 @@ def generate_hplus_hcross_m( # Step 2: Transform modes to target domain. # This requires tapering of TD modes, and FFT to transform to FD. - # TEOB and SEOB have slightly different conditioning routine implementations - # TEOB uses gwsignal routines SEOB uses lal routines - if self.approximant_str in ["TEOBResumSDALI"]: - wfg_utils.taper_td_modes_in_place_gwsignal(hlm_td) - hlm_fd = wfg_utils.td_modes_to_fd_modes_gwsignal(hlm_td, self.domain) - else: - wfg_utils.taper_td_modes_in_place(hlm_td) - hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain) + wfg_utils.taper_td_modes_in_place(hlm_td) + hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain) elif self.approximant_str in ["SEOBNRv5PHM", "SEOBNRv5HM"]: # Step 1: generate waveform modes in L0 frame in native domain of @@ -1395,10 +1389,6 @@ def generate_TD_modes_L0(self, parameters): {**parameters, "f_ref": self.f_ref} ) - generator = new_interface_get_waveform_generator(self.approximant_str) - hlm_td = gws_wfm.GenerateTDModes(parameters_gwsignal, generator) - hlms_lal = {} - # TEOBREsumS returns the modes unscaled in the distance version when calling GenerateTDModes # therefore, we rescale them here if "TEOB" in self.approximant_str: @@ -1418,6 +1408,11 @@ def generate_TD_modes_L0(self, parameters): else: distance_rescaling = 1.0 + generator = new_interface_get_waveform_generator(self.approximant_str) + hlm_td = gws_wfm.GenerateTDModes(parameters_gwsignal, generator) + hlms_lal = {} + + for key, value in hlm_td.items(): if type(key) != str: hlm_lal = lal.CreateCOMPLEX16TimeSeries( diff --git a/dingo/gw/waveform_generator/wfg_utils.py b/dingo/gw/waveform_generator/wfg_utils.py index fab45e64b..74144f209 100644 --- a/dingo/gw/waveform_generator/wfg_utils.py +++ b/dingo/gw/waveform_generator/wfg_utils.py @@ -46,7 +46,8 @@ def get_tapering_window_for_complex_time_series(h, tapering_flag: int = 1): ) h_tapered.data.data = h.data.data.copy().real LS.SimInspiralREAL8WaveTaper(h_tapered.data, tapering_flag) - eps = 1e-20 * np.max(np.abs(h.data.data)) + # in case h.data.data is an array of 0's this makes the window return ones + eps = 1e-20 * np.max(np.abs(h.data.data)) if np.max(np.abs(h.data.data)) > 0 else 1e-20 window = (np.abs(h_tapered.data.data) + eps) / (np.abs(h.data.data.real) + eps) # FIXME: using eps for numerical stability is not really robust here return window From 8db341d4e05b60a57d9c557d8411d35845672972 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 2 Feb 2026 22:33:48 +0100 Subject: [PATCH 37/48] updated teob functions --- .../waveform_generator/waveform_generator.py | 101 ++++- dingo/gw/waveform_generator/wfg_utils.py | 410 +++++++++++++++++- dingo/pipe/parser.py | 17 + 3 files changed, 496 insertions(+), 32 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index afb185572..c91a02039 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -3,6 +3,7 @@ from math import isclose import numpy as np +import astropy.constants as ac import astropy.units as u import astropy.constants as ac from typing import Dict, List, Tuple, Union, Callable @@ -1058,7 +1059,6 @@ def _convert_parameters( "condition": 1, } - if self.approximant_str in ["SEOBNRv5EHM", "TEOBResumSDALI"]: # eccentric parameters log_ecc = p.get("log10_eccentricity") @@ -1085,7 +1085,13 @@ def _convert_parameters( 'meanPerAno' : mean_per_ano * u.rad, }) - params_gwsignal.update(self.extra_wf_kwargs) + # TEOBResumSDALI doesn't support lmax_nyquist, so filter it out + if "TEOBResum" in self.approximant_str: + extra_wf_kwargs_filtered = {k: v for k, v in self.extra_wf_kwargs.items() + if k != "lmax_nyquist"} + params_gwsignal.update(extra_wf_kwargs_filtered) + else: + params_gwsignal.update(self.extra_wf_kwargs) if return_target_function: # This is a hack to make compatible with LAL version. Target functions for @@ -1244,8 +1250,8 @@ def generate_hplus_hcross_m( generator = new_interface_get_waveform_generator(self.approximant_str) if isinstance(self.domain, UniformFrequencyDomain): - # Generate FD modes in for frequencies [-f_max, ..., 0, ..., f_max]. - if self.approximant_str in ["SEOBNRv5EHM", "TEOBResumSDALI"]: + # Generate FD modes for frequencies [-f_max, ..., 0, ..., f_max]. + if self.approximant_str in ["SEOBNRv5EHM"]: # SEOBNRv5EHM and TEOBResumSDALI use these conditioning routines # assert LS.SimInspiralImplementedTDApproximants(self.approximant) # Step 1: generate waveform modes in L0 frame in native domain of @@ -1257,6 +1263,11 @@ def generate_hplus_hcross_m( wfg_utils.taper_td_modes_in_place(hlm_td) hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain) + # Step 3: Two-sided modes -> polarizations organized by m. + pol_m = wfg_utils.get_polarizations_from_fd_modes_m( + hlm_fd, iota, parameters["phase"] + ) + elif self.approximant_str in ["SEOBNRv5PHM", "SEOBNRv5HM"]: # Step 1: generate waveform modes in L0 frame in native domain of # approximant (here: TD), applying standard conditioning @@ -1267,6 +1278,75 @@ def generate_hplus_hcross_m( # Step 2: Transform modes to target domain. hlm_fd = wfg_utils.td_modes_to_fd_modes(hlm_td, self.domain) + # Step 3: Two-sided modes -> polarizations organized by m. + pol_m = wfg_utils.get_polarizations_from_fd_modes_m( + hlm_fd, iota, parameters["phase"] + ) + + elif "TEOB" in self.approximant_str: + # Step 1: generate waveform modes in L0 frame in native domain of + # approximant (here: TD). + # TEOB modes from GenerateTDModes are unconditioned (gwsignal's + # GenerateTDModes bypasses the conditioning pipeline). Generate + # at lower f_start for tapering buffer, then apply conditioning. + parameters_gwsignal = self._convert_parameters( + {**parameters, "f_ref": self.f_ref} + ) + f_start, t_extra, f_min, original_f_min, f_isco = ( + wfg_utils.get_conditioning_params_for_TEOB(parameters_gwsignal) + ) + parameters_gwsignal["f22_start"] = f_start * u.Hz + + generator = new_interface_get_waveform_generator( + self.approximant_str + ) + hlm_td_gwpy = gws_wfm.GenerateTDModes( + parameters_gwsignal, generator + ) + hlm_td = { + key: value + for key, value in hlm_td_gwpy.items() + if type(key) != str + } + iota = parameters_gwsignal["inclination"].value + + # GenerateTDModes returns dimensionless strain modes; + # rescale by -nu * M_total / distance * (G/c^2) to get + # physical units consistent with GenerateFDWaveform. + m1 = parameters_gwsignal["mass1"] + m2 = parameters_gwsignal["mass2"] + nu = m1 * m2 / (m1 + m2) ** 2 + distance_rescaling = -( + ( + nu + * (m1 + m2) + / parameters_gwsignal["distance"] + * ac.G + / ac.c ** 2 + ) + .to(u.dimensionless_unscaled) + .value + ) + for lm in hlm_td: + hlm_td[lm] = hlm_td[lm] * distance_rescaling + + wfg_utils.condition_td_modes_for_TEOB_in_place( + hlm_td, t_extra, f_min, original_f_min, f_isco + ) + + # Step 2: Transform modes to target domain. + hlm_fd = wfg_utils.td_modes_to_fd_modes_gwpy( + hlm_td, self.domain + ) + + # Step 3: One-sided modes -> polarizations organized by m. + # TEOB modes from gwpy FFT are one-sided (f >= 0 only), so we + # use the non-precessing symmetry h_{l,-m}(t) = (-1)^l h*_{lm}(t) + # to directly compute h+ and hx from the one-sided modes. + pol_m = wfg_utils.get_polarizations_from_onesided_fd_modes_m( + hlm_fd, iota, parameters["phase"] + ) + elif generator.domain == "freq": # Step 1: generate waveform modes in L0 frame in native domain of # approximant (here: FD) @@ -1275,17 +1355,14 @@ def generate_hplus_hcross_m( # Step 2: Transform modes to target domain. # Not required here, as approximant domain and target domain are both FD. + # Step 3: Two-sided modes -> polarizations organized by m. + pol_m = wfg_utils.get_polarizations_from_fd_modes_m( + hlm_fd, iota, parameters["phase"] + ) + else: raise NotImplementedError() - - # Step 3: Separate negative and positive frequency parts of the modes, - # and add contributions according to their transformation behavior under - # phase shifts. - pol_m = wfg_utils.get_polarizations_from_fd_modes_m( - hlm_fd, iota, parameters["phase"] - ) - else: raise NotImplementedError( f"Target domain of type {type(self.domain)} not yet implemented." diff --git a/dingo/gw/waveform_generator/wfg_utils.py b/dingo/gw/waveform_generator/wfg_utils.py index 74144f209..69f191f85 100644 --- a/dingo/gw/waveform_generator/wfg_utils.py +++ b/dingo/gw/waveform_generator/wfg_utils.py @@ -1,8 +1,10 @@ import numpy as np import lal import lalsimulation as LS -from lalsimulation.gwsignal.core import conditioning_subroutines as cond +from scipy.signal import butter, sosfiltfilt +from lalsimulation.gwsignal.core import conditioning_subroutines as cond import astropy.units as u +from gwpy.timeseries import TimeSeries def linked_list_modes_to_dict_modes(hlm_ll): @@ -197,30 +199,64 @@ def td_modes_to_fd_modes(hlm_td, domain): return hlm_fd -def td_modes_to_fd_modes_gwsignal(hlm_td, domain): - """ - TODO +def td_modes_to_fd_modes_gwpy(hlm_td, domain): + """ + Transform dict of gwpy TD modes to dict of one-sided FD modes via FFT. + The td modes are expected to be conditioned (tapered/filtered). + + Uses numpy's complex FFT (not gwpy's rfft-based .fft()) to correctly + handle complex-valued modes, then takes the positive-frequency half. + Applies a time shift correction to account for the epoch (physical start + time of the resized time series), matching the convention used by the + LAL-based td_modes_to_fd_modes and the reference FD waveform path. + + Parameters + ---------- + hlm_td: dict + Dictionary with (l,m) keys and gwpy TimeSeries objects for the + corresponding conditioned modes. + domain: dingo.gw.domains.UniformFrequencyDomain + Target domain after FFT. + + Returns + ------- + hlm_fd: dict + Dictionary with (l,m) keys and numpy arrays with the corresponding + one-sided FD modes on [0, df, ..., f_max]. """ hlm_fd = {} - + delta_f = domain.delta_f delta_t = 0.5 / domain.f_max - f_nyquist = domain.f_max # use f_max as f_nyquist + f_nyquist = domain.f_max chirplen = int(2 * f_nyquist / delta_f) - - for mode in hlm_td.keys(): - hlm_td_cond = cond.resize_gwpy_timeseries(hlm_td[mode], len(hlm_td[mode]) - chirplen, chirplen) - hlm_fd_cond = hlm_td_cond.fft() - hlm_fd_cond.epoch = hlm_td_cond.t0 - delta_f = hlm_fd_cond.df - hlm_fd_cond = (hlm_fd_cond / (2 * delta_f)).value - - # setting up reflection around boundary to create two-sided FFT - # hlm_fd_cond = np.concatenate((hlm_fd_cond[::-1], np.conj(hlm_fd_cond[1:])), axis=0) - - # hlm_fd_cond[-1] = hlm_fd_cond[0] - # hlm_fd[mode] = hlm_fd[mode] / (2 * delta_f) - hlm_fd[mode] = hlm_fd_cond + frequency_array = domain() # [0, df, ..., f_max] + + # Resize all modes and find consistent epoch from (2,2) mode + resized_modes = {} + epoch = None + for (l, m), h in hlm_td.items(): + h_resized = cond.resize_gwpy_timeseries(h, len(h) - chirplen, chirplen) + resized_modes[(l, m)] = h_resized + # Use (2,2) mode (or first mode) to determine epoch + if (l, m) == (2, 2) or epoch is None: + peak_idx = np.argmax(np.abs(np.asarray(h_resized))) + epoch = -peak_idx * delta_t + + # Precompute time shift (same for all modes) + dt_shift = 1.0 / delta_f + epoch + time_shift = np.exp(-1j * 2 * np.pi * dt_shift * frequency_array) + + for (l, m), h in resized_modes.items(): + data = np.asarray(h, dtype=complex) + # Complex FFT -> take positive-frequency half + fft_full = np.fft.fft(data, n=chirplen) + fft_pos = fft_full[: chirplen // 2 + 1] + # Standard continuous FT normalization: H(f) = dt * DFT + hf = fft_pos * delta_t + # Time shift correction to place merger at t=0 + hf *= time_shift + hlm_fd[(l, m)] = hf return hlm_fd @@ -260,6 +296,58 @@ def get_polarizations_from_fd_modes_m(hlm_fd, iota, phase): return pol_m +def get_polarizations_from_onesided_fd_modes_m(hlm_fd, iota, phase): + """Combine one-sided (f >= 0) FD modes into polarizations organized by m. + + For one-sided FD modes from non-precessing waveform models, the polarizations + are recovered using the symmetry h_{l,-m}(t) = (-1)^l h*_{lm}(t), which + relates the negative-frequency content to the positive-frequency modes: + + h+(f) = 0.5 sum_{l,m} [Y_{lm} + (-1)^l conj(Y_{l,-m})] H_{lm}(f) + hx(f) = 0.5i sum_{l,m} [Y_{lm} - (-1)^l conj(Y_{l,-m})] H_{lm}(f) + + where Y_{lm} = _{-2}Y_{lm}(iota, pi/2 - phase). + + Each pol_m[m] transforms as exp(-1j * m * phase) under phase shifts. + + Parameters + ---------- + hlm_fd : dict + Dictionary with (l, m) keys and one-sided FD mode arrays (numpy arrays + or gwpy FrequencySeries) as values, defined for f >= 0 only. + iota : float + Inclination angle. + phase : float + Reference phase. + + Returns + ------- + pol_m : dict + Dictionary with integer m as keys and + {"h_plus": array, "h_cross": array} as values. + """ + pol_m = {} + polarizations = ["h_plus", "h_cross"] + + for (l, m), h in hlm_fd.items(): + if m not in pol_m: + pol_m[m] = {k: 0.0 for k in polarizations} + + # Extract numpy array from gwpy FrequencySeries if needed + h_data = h.value if hasattr(h, "value") else h + + ylm = lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, m) + ylm_neg_conj = np.conj( + lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, -m) + ) + + sign = (-1) ** l + pol_m[m]["h_plus"] += 0.5 * (ylm + (sign * ylm_neg_conj)) * h_data + pol_m[m]["h_cross"] += 0.5j * (ylm - (sign * ylm_neg_conj)) * h_data + + return pol_m + + def get_starting_frequency_for_SEOBRNRv5_conditioning(parameters): """ Compute starting frequency needed for having 3 extra cycles for tapering the TD modes. @@ -385,3 +473,285 @@ def taper_td_modes_for_SEOBRNRv5_extra_time( # return timeseries return h_return + + +def get_conditioning_params_for_TEOB(parameters_gwsignal): + """ + Compute conditioning parameters for TEOBResumSDALI waveforms. + + Replicates the parameter computation from gwsignal's + generate_conditioned_td_waveform_from_td + (lalsimulation/gwsignal/core/waveform_conditioning.py:34-111). + + This is needed because gwsignal's GenerateTDModes + (lalsimulation/gwsignal/core/waveform.py:668) returns unconditioned modes, + while GenerateFDWaveform routes through the conditioning pipeline. To make + mode-decomposed waveforms match the reference, the same conditioning must be + applied to each mode individually. + + Parameters + ---------- + parameters_gwsignal: dict + Dictionary of gwsignal parameters (from _convert_parameters). + + Returns + ------- + f_start: float + Lower starting frequency for generating modes with buffer for tapering. + t_extra: float + Duration of the extra time to taper at the beginning (for stage 1). + f_min: float + Effective minimum frequency, possibly lowered to fisco_9 + (for stage 2 beginning taper). + original_f_min: float + The originally requested starting frequency (for stage 1 high-pass). + f_isco: float + ISCO frequency at 6M (for stage 2 end taper). + """ + # waveform_conditioning.py:51-52 + extra_time_fraction = 0.1 + extra_cycles = 3.0 + + # waveform_conditioning.py:55-61 + # Masses are stored as astropy Quantities in solar masses; .value gives + # solar masses, .si.value gives kg. LAL bound functions expect SI (kg). + f_min = parameters_gwsignal["f22_start"].value + m1_SI = parameters_gwsignal["mass1"].si.value + m2_SI = parameters_gwsignal["mass2"].si.value + s1z = parameters_gwsignal["spin1z"].value + s2z = parameters_gwsignal["spin2z"].value + original_f_min = f_min + + # waveform_conditioning.py:70-72 — clamp f_min to ISCO at 9^1.5 + fisco_9 = 1.0 / ( + np.power(9.0, 1.5) * np.pi * (m1_SI + m2_SI) * lal.MTSUN_SI / lal.MSUN_SI + ) + if f_min > fisco_9: + f_min = fisco_9 + + # waveform_conditioning.py:76 — upper bound on chirp time from f_min + tchirp = LS.SimInspiralChirpTimeBound(f_min, m1_SI, m2_SI, s1z, s2z) + + # waveform_conditioning.py:79 + s = LS.SimInspiralFinalBlackHoleSpinBound(s1z, s2z) + + # waveform_conditioning.py:82 — merger + ringdown time + tmerge = ( + LS.SimInspiralMergeTimeBound(m1_SI, m2_SI) + + LS.SimInspiralRingdownTimeBound(m1_SI + m2_SI, s) + ) + + # waveform_conditioning.py:87 — extra cycles near merger + textra = extra_cycles / f_min + + # waveform_conditioning.py:90 — lower starting frequency for tapering buffer + f_start = LS.SimInspiralChirpStartFrequencyBound( + (1.0 + extra_time_fraction) * tchirp + tmerge + textra, m1_SI, m2_SI + ) + + # waveform_conditioning.py:106 — t_extra for stage 1 cosine taper + t_extra = extra_time_fraction * tchirp + textra + + # waveform_conditioning.py:108 — ISCO at 6M for stage 2 end taper + f_isco = 1.0 / ( + np.power(6.0, 1.5) * np.pi * (m1_SI + m2_SI) * lal.MTSUN_SI / lal.MSUN_SI + ) + + return f_start, t_extra, f_min, original_f_min, f_isco + + +def _high_pass_complex(data, dt, f_min, attenuation=0.99, order=8): + """ + High-pass filter a complex time series using a Butterworth IIR filter. + + Replicates gwsignal's high_pass_time_series + (lalsimulation/gwsignal/core/conditioning_subroutines.py:10-46). + + Parameters + ---------- + data: np.ndarray (complex) + Time series data array. + dt: float + Sampling interval in seconds. + f_min: float + Minimum frequency for high-pass. + attenuation: float + Attenuation at the low-frequency cutoff (default 0.99). + order: int + Order of Butterworth filter (default 8). + + Returns + ------- + filtered: np.ndarray (complex) + High-pass filtered data. + """ + # conditioning_subroutines.py:33 — sampling frequency + fs = 1.0 / dt + + # conditioning_subroutines.py:37-39 — bilinear transform to compute cutoff + w1 = np.tan(np.pi * f_min * dt) + wc = w1 * (1.0 / attenuation**0.5 - 1) ** (1.0 / (2.0 * order)) + fc = fs * np.arctan(wc) / np.pi + + # conditioning_subroutines.py:42-43 — forward-backward Butterworth filter + sos = butter(order, fc, btype="highpass", output="sos", fs=fs) + filtered_re = sosfiltfilt(sos, data.real) + filtered_im = sosfiltfilt(sos, data.imag) + return filtered_re + 1j * filtered_im + + +def _condition_stage1_complex(data, dt, t_extra, f_min): + """ + Stage 1 conditioning: cosine taper at the beginning + high-pass filter. + + Replicates gwsignal's time_array_condition_stage1 + (lalsimulation/gwsignal/core/conditioning_subroutines.py:50-87). + + Parameters + ---------- + data: np.ndarray (complex) + Time series data array. + dt: float + Sampling interval in seconds. + t_extra: float + Duration of extra time at the beginning to taper. + f_min: float + Minimum frequency for high-pass filter. + + Returns + ------- + data: np.ndarray (complex) + Conditioned data. + """ + # conditioning_subroutines.py:71-77 — cosine (Hann) taper at beginning + Ntaper = int(np.round(t_extra / dt)) + if Ntaper > 0 and Ntaper < len(data): + taper_array = np.arange(Ntaper) + w = 0.5 - 0.5 * np.cos(taper_array * np.pi / Ntaper) + data[:Ntaper] *= w + + # conditioning_subroutines.py:80-81 — high-pass filter + data = _high_pass_complex(data, dt, f_min) + + # conditioning_subroutines.py:84-85 — trim trailing zeros + data = np.trim_zeros(data, trim='b') + + return data + + +def _condition_stage2_complex(data, dt, f_min, f_isco): + """ + Stage 2 conditioning: end taper (1 cycle at f_isco) + beginning taper + (1 cycle at f_min). + + Replicates gwsignal's time_array_condition_stage2 + (lalsimulation/gwsignal/core/conditioning_subroutines.py:90-144). + + Parameters + ---------- + data: np.ndarray (complex) + Time series data array. + dt: float + Sampling interval in seconds. + f_min: float + Minimum frequency for beginning taper. + f_isco: float + ISCO frequency for end taper. + + Returns + ------- + data: np.ndarray (complex) + Conditioned data. + """ + # conditioning_subroutines.py:111-114 + min_taper_samples = 4 + Nsize = len(data) + if Nsize < 2 * min_taper_samples: + return data + + # conditioning_subroutines.py:119-129 — end taper: 1 cycle at f_isco + ntaper_end = max(int(np.round(1.0 / (f_isco * dt))), min_taper_samples) + taper_array = np.arange(1, ntaper_end) + w_end = 0.5 - 0.5 * np.cos(taper_array * np.pi / ntaper_end) + data[Nsize - ntaper_end + 1 :] *= w_end[::-1] + + # conditioning_subroutines.py:133-142 — beginning taper: 1 cycle at f_min + ntaper_begin = max(int(np.round(1.0 / (f_min * dt))), min_taper_samples) + taper_array = np.arange(ntaper_begin) + w_begin = 0.5 - 0.5 * np.cos(taper_array * np.pi / ntaper_begin) + data[:ntaper_begin] *= w_begin + + return data + + +def condition_td_modes_for_TEOB_in_place(hlm_td, t_extra, f_min, original_f_min, f_isco): + """ + Apply gwsignal-style conditioning to TD modes in place. + + This replicates the conditioning that gwsignal's + generate_conditioned_td_waveform_from_td + (lalsimulation/gwsignal/core/waveform_conditioning.py:34-111) applies to the + combined TD polarizations, but applies it to individual complex modes. This + is necessary because GenerateTDModes (gwsignal/core/waveform.py:668) returns + unconditioned modes — it bypasses the conditioning pipeline entirely. + + The conditioning consists of two stages (following the C implementation at + XLALSimInspiralTDConditionStage1 / Stage2): + + Stage 1 (conditioning_subroutines.py:50-87): + - Cosine (Hann) taper over t_extra seconds at the beginning + - Order-8 Butterworth high-pass filter at original_f_min (attenuation 0.99) + + Stage 2 (conditioning_subroutines.py:90-144): + - End taper: 1 cycle at f_isco (ISCO at 6M) + - Beginning taper: 1 cycle at f_min + + Usage + ----- + The modes must be generated at the lower starting frequency f_start returned + by get_conditioning_params_for_TEOB(), so that there is sufficient buffer + for the tapering. + + Example:: + + f_start, t_extra, f_min, original_f_min, f_isco = ( + get_conditioning_params_for_TEOB(parameters_gwsignal) + ) + # Generate modes at f_start (not the original f_min) + hlm_td, iota = generate_TD_modes_at_f_start(parameters, f_start) + condition_td_modes_for_TEOB_in_place( + hlm_td, t_extra, f_min, original_f_min, f_isco + ) + hlm_fd = td_modes_to_fd_modes(hlm_td, domain) + + Parameters + ---------- + hlm_td: dict + Dictionary with (l,m) keys and complex LAL COMPLEX16TimeSeries objects. + t_extra: float + Duration of extra time at the beginning for stage 1 cosine taper. + From get_conditioning_params_for_TEOB(). + f_min: float + Effective minimum frequency (possibly lowered to fisco_9). + Used for stage 2 beginning taper. + original_f_min: float + The originally requested starting frequency. + Used for stage 1 high-pass filter. + f_isco: float + ISCO frequency at 6M. Used for stage 2 end taper. + """ + for lm in list(hlm_td.keys()): + h = hlm_td[lm] + dt = h.dt.value + data = h.value.copy() + + # waveform_conditioning.py:106 — stage 1 uses original_f_min for high-pass + data = _condition_stage1_complex(data, dt, t_extra, original_f_min) + + # waveform_conditioning.py:109 — stage 2 uses adjusted f_min and f_isco + data = _condition_stage2_complex(data, dt, f_min, f_isco) + + # Replace with new TimeSeries (length may have changed due to trim) + hlm_td[lm] = TimeSeries( + data, t0=h.t0, dt=h.dt, unit=h.unit + ) diff --git a/dingo/pipe/parser.py b/dingo/pipe/parser.py index ce22535b8..e61a801bc 100644 --- a/dingo/pipe/parser.py +++ b/dingo/pipe/parser.py @@ -80,6 +80,23 @@ def create_parser(top_level=True, usage=None): ) calibration_parser.add( +<<<<<<< Updated upstream +======= + "--calibration-mode", + type=nonestr, + default=None, + choices=["marginalize", "sample", None], + help=( + "How to handle calibration uncertainty. 'marginalize' averages likelihoods " + "over multiple calibration draws. 'sample' treats calibration parameters as " + "part of the posterior (importance sampling over calibration). If None or " + "calibration-model is None, no calibration handling is performed. " + "Default is 'marginalize'." + ), + ) + + calibration_parser.add( +>>>>>>> Stashed changes "--calibration-correction-type", type=nonestr, default="data", From 59387b07dea4d7e688592b328950eab07ce27233 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 16 Feb 2026 14:55:39 +0100 Subject: [PATCH 38/48] implemented TEOB synthetic phase' --- .../waveform_generator/waveform_generator.py | 21 +- dingo/gw/waveform_generator/wfg_utils.py | 194 +++++++++++++++--- 2 files changed, 169 insertions(+), 46 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index c91a02039..eff471fb8 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -1033,7 +1033,7 @@ def _convert_parameters( f_min = self.domain.f_min # for SEOBNRv5EHM, the starting frequency must be the same as the reference frequency # note this is the orbit averaged reference frequency - if self.approximant_str == "SEOBNRv5EHM": + if self.approximant_str == "SEOBNRv5EHM" or self.approximant_str == "TEOBResumSDALI": f_min = self.f_ref # parameters needed for TD waveforms delta_t = 0.5 / self.domain.f_max @@ -1287,15 +1287,12 @@ def generate_hplus_hcross_m( # Step 1: generate waveform modes in L0 frame in native domain of # approximant (here: TD). # TEOB modes from GenerateTDModes are unconditioned (gwsignal's - # GenerateTDModes bypasses the conditioning pipeline). Generate - # at lower f_start for tapering buffer, then apply conditioning. + # GenerateTDModes bypasses the conditioning pipeline). Apply the + # same fallback conditioning (sigmoid start taper) that gwsignal + # uses for TEOB in generate_conditioned_td_waveform_from_td_fallback. parameters_gwsignal = self._convert_parameters( {**parameters, "f_ref": self.f_ref} ) - f_start, t_extra, f_min, original_f_min, f_isco = ( - wfg_utils.get_conditioning_params_for_TEOB(parameters_gwsignal) - ) - parameters_gwsignal["f22_start"] = f_start * u.Hz generator = new_interface_get_waveform_generator( self.approximant_str @@ -1311,12 +1308,12 @@ def generate_hplus_hcross_m( iota = parameters_gwsignal["inclination"].value # GenerateTDModes returns dimensionless strain modes; - # rescale by -nu * M_total / distance * (G/c^2) to get + # rescale by nu * M_total / distance * (G/c^2) to get # physical units consistent with GenerateFDWaveform. m1 = parameters_gwsignal["mass1"] m2 = parameters_gwsignal["mass2"] nu = m1 * m2 / (m1 + m2) ** 2 - distance_rescaling = -( + distance_rescaling = ( ( nu * (m1 + m2) @@ -1330,13 +1327,11 @@ def generate_hplus_hcross_m( for lm in hlm_td: hlm_td[lm] = hlm_td[lm] * distance_rescaling - wfg_utils.condition_td_modes_for_TEOB_in_place( - hlm_td, t_extra, f_min, original_f_min, f_isco - ) + wfg_utils.taper_td_gwpy_modes_in_place(hlm_td, iota) # Step 2: Transform modes to target domain. hlm_fd = wfg_utils.td_modes_to_fd_modes_gwpy( - hlm_td, self.domain + hlm_td, self.domain, iota ) # Step 3: One-sided modes -> polarizations organized by m. diff --git a/dingo/gw/waveform_generator/wfg_utils.py b/dingo/gw/waveform_generator/wfg_utils.py index 69f191f85..707449d08 100644 --- a/dingo/gw/waveform_generator/wfg_utils.py +++ b/dingo/gw/waveform_generator/wfg_utils.py @@ -1,7 +1,7 @@ import numpy as np import lal import lalsimulation as LS -from scipy.signal import butter, sosfiltfilt +from scipy.signal import butter, sosfiltfilt, find_peaks from lalsimulation.gwsignal.core import conditioning_subroutines as cond import astropy.units as u from gwpy.timeseries import TimeSeries @@ -54,7 +54,6 @@ def get_tapering_window_for_complex_time_series(h, tapering_flag: int = 1): # FIXME: using eps for numerical stability is not really robust here return window - def taper_td_modes_in_place(hlm_td, tapering_flag: int = 1): """ Taper the time domain modes in place. @@ -74,23 +73,128 @@ def taper_td_modes_in_place(hlm_td, tapering_flag: int = 1): window = get_tapering_window_for_complex_time_series(h, tapering_flag) h.data.data *= window -def taper_td_modes_in_place_gwsignal(hlm_td, tapering_flag: str = 'start'): +def taper_td_gwpy_modes_in_place(hlm_td, iota, phase=0.0): """ - Taper the time domain modes in place using gwsignal conditioning routines. - There are apparently slight differences in the tapering functions - between these two and since some wfs (EG TEOB) uses the gwsignal - routine, we use the gwsignal routine sometimes. + Apply sigmoid start taper to gwpy TimeSeries modes in place. + + This replicates the fallback conditioning used by gwsignal for + approximants with f_ref_spin=False (e.g. TEOBResumSDALI). + See generate_conditioned_td_waveform_from_td_fallback in + lalsimulation/gwsignal/core/waveform_conditioning.py:35-64. + + The gwsignal reference applies taper_gwpy_timeseries to the combined + polarizations hp and hc, where find_peaks determines the taper length. + Since hp and hc can yield different taper lengths, but we need a single + real-valued window to apply to each complex mode, we compute both windows + and use the one with shorter taper (preserving more signal). Parameters ---------- hlm_td: dict - Dictionary with (l,m) keys and the complex lal time series objects for the - corresponding modes. + Dictionary with (l,m) keys and gwpy TimeSeries objects. + iota: float + Inclination angle in radians. + phase: float + Reference phase for computing the combined polarizations used to + determine the taper window. Defaults to 0.0 since modes should be + phase-independent. + """ + max_len = max(len(h) for h in hlm_td.values()) + + # Sum modes to get combined hp and hc (unconditioned). + # These are needed to compute taper windows that match the reference path, + # where taper_gwpy_timeseries uses find_peaks on the combined polarizations. + h_complex = np.zeros(max_len, dtype=complex) + for (l, m), h in hlm_td.items(): + data = h.value.copy() + if len(data) < max_len: + data = np.pad(data, (0, max_len - len(data))) + ylm = lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, m) + h_complex += ylm * data + + hp_combined = h_complex.real + hc_combined = -h_complex.imag + + # Compute separate taper windows for hp and hc. The reference path + # (generate_conditioned_td_waveform_from_td_fallback) applies + # taper_gwpy_timeseries independently to hp and hc, which can yield + # different taper lengths because find_peaks finds different peaks. + # We apply W_hp to Re(h_lm) and W_hc to Im(h_lm) for each mode. + # For non-precessing systems the modes are purely real in the co-rotating + # frame, so Re(h_lm) and Im(h_lm) map cleanly to the two polarizations, + # giving exact agreement with the reference. + start_hp, n_hp = _compute_sigmoid_taper_params(hp_combined) + start_hc, n_hc = _compute_sigmoid_taper_params(hc_combined) + + if start_hp is None and start_hc is None: + return + + # Build windows, falling back to the other if one polarization is empty + if start_hp is None: + start_hp, n_hp = start_hc, n_hc + if start_hc is None: + start_hc, n_hc = start_hp, n_hp + + W_hp = _compute_sigmoid_window(max_len, start_hp, n_hp) + W_hc = _compute_sigmoid_window(max_len, start_hc, n_hc) + + for (l, m) in hlm_td: + data = hlm_td[(l, m)].value.copy() + if len(data) < max_len: + data = np.pad(data, (0, max_len - len(data))) + data = W_hp * data.real + 1j * W_hc * data.imag + hlm_td[(l, m)] = TimeSeries( + data, + t0=hlm_td[(l, m)].t0, + dt=hlm_td[(l, m)].dt, + ) +def _compute_sigmoid_taper_params(signal): + """Compute (start, n) for the sigmoid start taper. + + Replicates the peak-finding logic in taper_gwpy_timeseries for + taper_kind='start'. Returns (None, None) if the signal is empty + or too short. """ - lalseries_to_gwpy_timeseries_in_place(hlm_td) - for mode in hlm_td.keys(): - hlm_td[mode] = cond.taper_gwpy_timeseries(hlm_td[mode], tapering_flag) + LALSIMULATION_RINGING_EXTENT = 19 + + start = -1 + for idx, val in enumerate(signal): + if val != 0: + start = idx + break + if start == -1: + return None, None + + end = -1 + for idx, val in enumerate(signal[::-1]): + if val != 0: + end = len(signal) - 1 - idx + break + + if (end - start) <= 1: + return None, None + + mid = int((start + end) / 2) + pks, _ = find_peaks(abs(signal[start + 1 : mid])) + pks = pks[pks > LALSIMULATION_RINGING_EXTENT] + + if len(pks) < 2: + n = mid - start + else: + n = pks[1] + 1 + + return start, n + +def _compute_sigmoid_window(length, start, n): + """Compute the sigmoid taper window matching taper_gwpy_timeseries.""" + window = np.ones(length) + window[start] = 0.0 + realI = np.arange(1, n - 1) + z = (n - 1.0) / realI + (n - 1.0) / (realI - (n - 1.0)) + sigma = 1.0 / (np.exp(z) + 1.0) + window[start + 1 : start + n - 1] = sigma + return window def lalseries_to_gwpy_timeseries_in_place(hlm_td): """ @@ -199,16 +303,25 @@ def td_modes_to_fd_modes(hlm_td, domain): return hlm_fd -def td_modes_to_fd_modes_gwpy(hlm_td, domain): +def td_modes_to_fd_modes_gwpy(hlm_td, domain, iota, phase=0.0): """ Transform dict of gwpy TD modes to dict of one-sided FD modes via FFT. The td modes are expected to be conditioned (tapered/filtered). - Uses numpy's complex FFT (not gwpy's rfft-based .fft()) to correctly - handle complex-valued modes, then takes the positive-frequency half. - Applies a time shift correction to account for the epoch (physical start - time of the resized time series), matching the convention used by the - LAL-based td_modes_to_fd_modes and the reference FD waveform path. + Replicates the FFT pipeline in gwsignal's + generate_conditioned_fd_waveform_from_td (waveform_conditioning.py:510-522): + resize → fft → normalize by 1/(2*df) → set epoch + followed by the time shift from generate_FD_waveform + (waveform_generator.py:1178-1182): + dt = 1/df + epoch; h *= exp(-i*2π*dt*f) + + The epoch is computed by reconstructing the real h+ polarization from + the modes after resize and finding its peak with np.argmax, matching + the reference path where resize_gwpy_timeseries uses np.argmax(hp). + + Since the modes are complex-valued and gwpy's fft uses rfft (real input + only), each mode is split into real and imaginary parts which are FFT'd + separately and recombined: H_lm(f) = FFT(Re) + i*FFT(Im). Parameters ---------- @@ -217,6 +330,11 @@ def td_modes_to_fd_modes_gwpy(hlm_td, domain): corresponding conditioned modes. domain: dingo.gw.domains.UniformFrequencyDomain Target domain after FFT. + iota: float + Inclination angle in radians. Used to reconstruct h+ from modes + for peak-finding (epoch computation). + phase: float + Reference phase for reconstructing h+. Defaults to 0.0. Returns ------- @@ -232,31 +350,41 @@ def td_modes_to_fd_modes_gwpy(hlm_td, domain): chirplen = int(2 * f_nyquist / delta_f) frequency_array = domain() # [0, df, ..., f_max] - # Resize all modes and find consistent epoch from (2,2) mode + # Resize modes to chirplen (waveform_conditioning.py:510) resized_modes = {} - epoch = None for (l, m), h in hlm_td.items(): h_resized = cond.resize_gwpy_timeseries(h, len(h) - chirplen, chirplen) resized_modes[(l, m)] = h_resized - # Use (2,2) mode (or first mode) to determine epoch - if (l, m) == (2, 2) or epoch is None: - peak_idx = np.argmax(np.abs(np.asarray(h_resized))) - epoch = -peak_idx * delta_t + + # Reconstruct h+ from resized modes to find the peak, matching + # the reference path where resize_gwpy_timeseries uses np.argmax(hp) + # on the real h+ polarization. + hp_reconstructed = np.zeros(chirplen) + for (l, m), h in resized_modes.items(): + ylm = lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, m) + hp_reconstructed += (ylm * np.asarray(h)).real + epoch = -np.argmax(hp_reconstructed) * delta_t # Precompute time shift (same for all modes) dt_shift = 1.0 / delta_f + epoch time_shift = np.exp(-1j * 2 * np.pi * dt_shift * frequency_array) for (l, m), h in resized_modes.items(): - data = np.asarray(h, dtype=complex) - # Complex FFT -> take positive-frequency half - fft_full = np.fft.fft(data, n=chirplen) - fft_pos = fft_full[: chirplen // 2 + 1] - # Standard continuous FT normalization: H(f) = dt * DFT - hf = fft_pos * delta_t - # Time shift correction to place merger at t=0 - hf *= time_shift - hlm_fd[(l, m)] = hf + # Split complex mode into real and imaginary TimeSeries for gwpy fft + h_re = TimeSeries(np.asarray(h).real, dt=delta_t) + h_im = TimeSeries(np.asarray(h).imag, dt=delta_t) + + # FFT + normalize (waveform_conditioning.py:513-521) + hf_re = h_re.fft() + hf_re = hf_re / (2 * hf_re.df) + + hf_im = h_im.fft() + hf_im = hf_im / (2 * hf_im.df) + + # Combine and apply time shift + hf_combined = (hf_re.value[:len(frequency_array)] + + 1j * hf_im.value[:len(frequency_array)]) + hlm_fd[(l, m)] = hf_combined * time_shift return hlm_fd From 206da83bace4bbb6c84cebf80481b7bc6ac03394 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 16 Feb 2026 18:53:18 +0100 Subject: [PATCH 39/48] Refactor TEOB mode generation and remove dead code - Consolidate TEOB mode generation into generate_TD_modes_L0 with return_lal parameter (True by default, False for TEOB gwpy pipeline). Fixes incorrect negative distance rescaling sign for TEOB. - Fix typo in TEOB branch condition: "TEOBREsumSDALI" -> "TEOBResumSDALI" - Remove dead code from wfg_utils.py: lalseries_to_gwpy_timeseries_in_place, gwpyseries_to_lalseries_in_place, and the entire Stage 1/2 TEOB conditioning chain (get_conditioning_params_for_TEOB, condition_td_modes_for_TEOB_in_place, _high_pass_complex, _condition_stage1_complex, _condition_stage2_complex) which was superseded by the sigmoid taper approach. - Clean up unused imports (butter, sosfiltfilt, astropy.units). Co-Authored-By: Claude Opus 4.6 --- .../waveform_generator/waveform_generator.py | 170 ++++--- dingo/gw/waveform_generator/wfg_utils.py | 474 +++++------------- 2 files changed, 207 insertions(+), 437 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index eff471fb8..44275f3cc 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -1283,55 +1283,26 @@ def generate_hplus_hcross_m( hlm_fd, iota, parameters["phase"] ) - elif "TEOB" in self.approximant_str: - # Step 1: generate waveform modes in L0 frame in native domain of - # approximant (here: TD). - # TEOB modes from GenerateTDModes are unconditioned (gwsignal's - # GenerateTDModes bypasses the conditioning pipeline). Apply the - # same fallback conditioning (sigmoid start taper) that gwsignal - # uses for TEOB in generate_conditioned_td_waveform_from_td_fallback. - parameters_gwsignal = self._convert_parameters( - {**parameters, "f_ref": self.f_ref} - ) - - generator = new_interface_get_waveform_generator( - self.approximant_str - ) - hlm_td_gwpy = gws_wfm.GenerateTDModes( - parameters_gwsignal, generator - ) - hlm_td = { - key: value - for key, value in hlm_td_gwpy.items() - if type(key) != str - } - iota = parameters_gwsignal["inclination"].value - - # GenerateTDModes returns dimensionless strain modes; - # rescale by nu * M_total / distance * (G/c^2) to get - # physical units consistent with GenerateFDWaveform. - m1 = parameters_gwsignal["mass1"] - m2 = parameters_gwsignal["mass2"] - nu = m1 * m2 / (m1 + m2) ** 2 - distance_rescaling = ( - ( - nu - * (m1 + m2) - / parameters_gwsignal["distance"] - * ac.G - / ac.c ** 2 - ) - .to(u.dimensionless_unscaled) - .value + elif self.approximant_str in ["TEOBResumSDALI"]: + # Step 1: generate waveform modes in L0 frame in native domain + # of approximant (here: TD). Returns gwpy TimeSeries for TEOB. + # Then apply sigmoid start taper (TEOB fallback conditioning). + hlm_td, iota = self.generate_TD_modes_L0( + parameters, return_lal=False ) - for lm in hlm_td: - hlm_td[lm] = hlm_td[lm] * distance_rescaling - wfg_utils.taper_td_gwpy_modes_in_place(hlm_td, iota) # Step 2: Transform modes to target domain. - hlm_fd = wfg_utils.td_modes_to_fd_modes_gwpy( - hlm_td, self.domain, iota + # FFT without time shift, then apply phase-dependent time + # shift at spin_conversion_phase (peak at t≈0 for production). + hlm_fd_raw, resized_td_modes = ( + wfg_utils.td_modes_to_fd_modes_gwpy(hlm_td, self.domain) + ) + hlm_fd, self._deferred_timeshift_data = ( + wfg_utils.apply_time_shift_to_fd_modes( + hlm_fd_raw, resized_td_modes, iota, + self.spin_conversion_phase, self.domain, + ) ) # Step 3: One-sided modes -> polarizations organized by m. @@ -1339,7 +1310,7 @@ def generate_hplus_hcross_m( # use the non-precessing symmetry h_{l,-m}(t) = (-1)^l h*_{lm}(t) # to directly compute h+ and hx from the one-sided modes. pol_m = wfg_utils.get_polarizations_from_onesided_fd_modes_m( - hlm_fd, iota, parameters["phase"] + hlm_fd, iota, self.spin_conversion_phase ) elif generator.domain == "freq": @@ -1436,7 +1407,7 @@ def generate_FD_modes_L0(self, parameters): # Pending to adapt f"to test that this works as intended! Ideally, add some unit tests." ) - def generate_TD_modes_L0(self, parameters): + def generate_TD_modes_L0(self, parameters, return_lal=True): """ Generate TD modes in the L0 frame. @@ -1445,31 +1416,39 @@ def generate_TD_modes_L0(self, parameters): parameters: dict Dictionary of parameters for the waveform. For details see see self.generate_hplus_hcross. + return_lal: bool + If True (default), return modes as LAL COMPLEX16TimeSeries. + If False, return modes as gwpy TimeSeries. Returns ------- hlm_td: dict - Dictionary with (l,m) as keys and the corresponding TD modes in lal format as + Dictionary with (l,m) as keys and the corresponding TD modes as values. iota: float """ # TD approximants that are implemented in L0 frame. Currently tested for: - # 52: SEOBNRv4PHM # SEOBNRV5EHM + # TEOBResumSDALI parameters_gwsignal = self._convert_parameters( {**parameters, "f_ref": self.f_ref} ) - # TEOBREsumS returns the modes unscaled in the distance version when calling GenerateTDModes - # therefore, we rescale them here + generator = new_interface_get_waveform_generator(self.approximant_str) + hlm_td_gwpy = gws_wfm.GenerateTDModes(parameters_gwsignal, generator) + iota = parameters_gwsignal["inclination"].value + + # TEOB GenerateTDModes returns dimensionless strain modes; + # rescale by nu * M_total / distance * (G/c^2) to get physical units. if "TEOB" in self.approximant_str: - m1, m2 = parameters_gwsignal["mass1"], parameters_gwsignal["mass2"] - nu = m1 * m2 / (m1 + m2) ** 2 - distance_rescaling = -( + m1 = parameters_gwsignal["mass1"] + m2 = parameters_gwsignal["mass2"] + nu = m1 * m2 / (m1 + m2) ** 2 + distance_rescaling = ( ( nu - * (parameters_gwsignal["mass1"] + parameters_gwsignal["mass2"]) + * (m1 + m2) / parameters_gwsignal["distance"] * ac.G / ac.c ** 2 @@ -1477,29 +1456,27 @@ def generate_TD_modes_L0(self, parameters): .to(u.dimensionless_unscaled) .value ) - else: - distance_rescaling = 1.0 - generator = new_interface_get_waveform_generator(self.approximant_str) - hlm_td = gws_wfm.GenerateTDModes(parameters_gwsignal, generator) - hlms_lal = {} - - - for key, value in hlm_td.items(): + hlm_td = {} + for key, value in hlm_td_gwpy.items(): if type(key) != str: - hlm_lal = lal.CreateCOMPLEX16TimeSeries( - "hplus", - value.epoch.value, - 0, - value.dt.value, - lal.DimensionlessUnit, - len(value), - ) - hlm_lal.data.data = value.value * distance_rescaling - hlms_lal[key] = hlm_lal - + if "TEOB" in self.approximant_str: + value = value * distance_rescaling + if return_lal: + hlm_lal = lal.CreateCOMPLEX16TimeSeries( + "hplus", + value.epoch.value, + 0, + value.dt.value, + lal.DimensionlessUnit, + len(value), + ) + hlm_lal.data.data = value.value + hlm_td[key] = hlm_lal + else: + hlm_td[key] = value - return hlms_lal, parameters_gwsignal["inclination"].value + return hlm_td, iota def generate_TD_modes_L0_conditioned_extra_time(self, parameters): """ @@ -1680,15 +1657,56 @@ def generate_waveforms_parallel( return polarizations -def sum_contributions_m(x_m, phase_shift=0.0): +def sum_contributions_m(x_m, phase_shift=0.0, deferred_timeshift_data=None): """ Sum the contributions over m-components, optionally introducing a phase shift. + + When deferred_timeshift_data is provided and phase_shift != 0, a time shift + correction is applied to account for the phase-dependent epoch. The modes in + x_m carry a time shift computed at a reference phase; this correction adjusts + it to the target phase (phase_ref + phase_shift). + + Parameters + ---------- + x_m: dict + Dictionary with integer m keys and dict values containing arrays. + phase_shift: float + Phase shift to apply via exp(-1j * m * phase_shift). + deferred_timeshift_data: dict or None + If provided, must contain: resized_td_modes, iota, phase_ref, dt_ref, + delta_t, delta_f, frequency_array. Stored on the waveform generator + instance as _deferred_timeshift_data after calling + generate_hplus_hcross_m for TEOB approximants. + + Returns + ------- + result: dict + Dictionary with the same keys as x_m's inner dicts (e.g. "h_plus", + "h_cross"), containing the summed arrays. """ keys = next(iter(x_m.values())).keys() result = {key: 0.0 for key in keys} for key in keys: for m, x in x_m.items(): result[key] += x[key] * np.exp(-1j * m * phase_shift) + + if deferred_timeshift_data is not None and phase_shift != 0.0: + target_phase = deferred_timeshift_data["phase_ref"] + phase_shift + epoch_target = wfg_utils.compute_epoch_from_resized_td_modes( + deferred_timeshift_data["resized_td_modes"], + deferred_timeshift_data["iota"], + target_phase, + deferred_timeshift_data["delta_t"], + ) + dt_target = 1.0 / deferred_timeshift_data["delta_f"] + epoch_target + dt_correction = dt_target - deferred_timeshift_data["dt_ref"] + correction = np.exp( + -1j * 2 * np.pi * dt_correction + * deferred_timeshift_data["frequency_array"] + ) + for key in result: + result[key] *= correction + return result diff --git a/dingo/gw/waveform_generator/wfg_utils.py b/dingo/gw/waveform_generator/wfg_utils.py index 707449d08..04ee92f9c 100644 --- a/dingo/gw/waveform_generator/wfg_utils.py +++ b/dingo/gw/waveform_generator/wfg_utils.py @@ -1,9 +1,8 @@ import numpy as np import lal import lalsimulation as LS -from scipy.signal import butter, sosfiltfilt, find_peaks +from scipy.signal import find_peaks from lalsimulation.gwsignal.core import conditioning_subroutines as cond -import astropy.units as u from gwpy.timeseries import TimeSeries @@ -196,46 +195,6 @@ def _compute_sigmoid_window(length, start, n): window[start + 1 : start + n - 1] = sigma return window -def lalseries_to_gwpy_timeseries_in_place(hlm_td): - """ - Convert lal time series in place to gwpy time series. - - Parameters - ---------- - hlm_td: dict - Dictionary with (l,m) keys and the complex lal time series objects for the - corresponding modes. - - """ - from gwpy.timeseries import TimeSeries - for mode in hlm_td.keys(): - times = ( - hlm_td[mode].epoch.gpsSeconds - + hlm_td[mode].epoch.gpsNanoSeconds * 1e-9 - + np.arange(hlm_td[mode].data.length) * hlm_td[mode].deltaT - ) - - hlm_td[mode] = TimeSeries( - data=hlm_td[mode].data.data, - times=times, - name=f"h_{mode[0]}_{mode[1]}", - unit=u.dimensionless_unscaled, - ) - -def gwpyseries_to_lalseries_in_place(hlm_td): - """ - Convert gwpy time series in place to lal time series. - - Parameters - ---------- - hlm_td: dict - Dictionary with (l,m) keys and the complex gwpy time series objects for the - corresponding modes. - - """ - for mode in hlm_td.keys(): - hlm_td[mode] = hlm_td[mode].to_lal() - def td_modes_to_fd_modes(hlm_td, domain): """ Transform dict of td modes to dict of fd modes via FFT. The td modes are expected @@ -303,21 +262,23 @@ def td_modes_to_fd_modes(hlm_td, domain): return hlm_fd -def td_modes_to_fd_modes_gwpy(hlm_td, domain, iota, phase=0.0): +def td_modes_to_fd_modes_gwpy(hlm_td, domain): """ Transform dict of gwpy TD modes to dict of one-sided FD modes via FFT. The td modes are expected to be conditioned (tapered/filtered). Replicates the FFT pipeline in gwsignal's generate_conditioned_fd_waveform_from_td (waveform_conditioning.py:510-522): - resize → fft → normalize by 1/(2*df) → set epoch - followed by the time shift from generate_FD_waveform - (waveform_generator.py:1178-1182): - dt = 1/df + epoch; h *= exp(-i*2π*dt*f) + resize → fft → normalize by 1/(2*df) + + No time shift is applied to the FD modes. The caller is responsible for + computing and applying the time shift (which depends on the target phase + via np.argmax of the reconstructed hp). This enables the "deferred time + shift" approach where modes are FFT'd once (phase-independently) and the + phase-dependent time shift is applied at resummation time. - The epoch is computed by reconstructing the real h+ polarization from - the modes after resize and finding its peak with np.argmax, matching - the reference path where resize_gwpy_timeseries uses np.argmax(hp). + Also returns the resized TD mode data (numpy arrays) needed for the + deferred time shift computation. Since the modes are complex-valued and gwpy's fft uses rfft (real input only), each mode is split into real and imaginary parts which are FFT'd @@ -330,19 +291,18 @@ def td_modes_to_fd_modes_gwpy(hlm_td, domain, iota, phase=0.0): corresponding conditioned modes. domain: dingo.gw.domains.UniformFrequencyDomain Target domain after FFT. - iota: float - Inclination angle in radians. Used to reconstruct h+ from modes - for peak-finding (epoch computation). - phase: float - Reference phase for reconstructing h+. Defaults to 0.0. Returns ------- hlm_fd: dict Dictionary with (l,m) keys and numpy arrays with the corresponding - one-sided FD modes on [0, df, ..., f_max]. + one-sided FD modes on [0, df, ..., f_max]. No time shift applied. + resized_td_modes: dict + Dictionary with (l,m) keys and numpy complex arrays of the resized + TD mode data, needed for deferred time shift computation. """ hlm_fd = {} + resized_td_modes = {} delta_f = domain.delta_f delta_t = 0.5 / domain.f_max @@ -350,29 +310,15 @@ def td_modes_to_fd_modes_gwpy(hlm_td, domain, iota, phase=0.0): chirplen = int(2 * f_nyquist / delta_f) frequency_array = domain() # [0, df, ..., f_max] - # Resize modes to chirplen (waveform_conditioning.py:510) - resized_modes = {} for (l, m), h in hlm_td.items(): - h_resized = cond.resize_gwpy_timeseries(h, len(h) - chirplen, chirplen) - resized_modes[(l, m)] = h_resized - - # Reconstruct h+ from resized modes to find the peak, matching - # the reference path where resize_gwpy_timeseries uses np.argmax(hp) - # on the real h+ polarization. - hp_reconstructed = np.zeros(chirplen) - for (l, m), h in resized_modes.items(): - ylm = lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, m) - hp_reconstructed += (ylm * np.asarray(h)).real - epoch = -np.argmax(hp_reconstructed) * delta_t - - # Precompute time shift (same for all modes) - dt_shift = 1.0 / delta_f + epoch - time_shift = np.exp(-1j * 2 * np.pi * dt_shift * frequency_array) + # Resize to chirplen (waveform_conditioning.py:510) + start_id = len(h) - chirplen + h_resized = cond.resize_gwpy_timeseries(h, start_id, chirplen) + resized_td_modes[(l, m)] = np.asarray(h_resized).copy() - for (l, m), h in resized_modes.items(): # Split complex mode into real and imaginary TimeSeries for gwpy fft - h_re = TimeSeries(np.asarray(h).real, dt=delta_t) - h_im = TimeSeries(np.asarray(h).imag, dt=delta_t) + h_re = TimeSeries(np.asarray(h_resized).real, dt=delta_t) + h_im = TimeSeries(np.asarray(h_resized).imag, dt=delta_t) # FFT + normalize (waveform_conditioning.py:513-521) hf_re = h_re.fft() @@ -381,12 +327,98 @@ def td_modes_to_fd_modes_gwpy(hlm_td, domain, iota, phase=0.0): hf_im = h_im.fft() hf_im = hf_im / (2 * hf_im.df) - # Combine and apply time shift - hf_combined = (hf_re.value[:len(frequency_array)] - + 1j * hf_im.value[:len(frequency_array)]) - hlm_fd[(l, m)] = hf_combined * time_shift + hlm_fd[(l, m)] = (hf_re.value[:len(frequency_array)] + + 1j * hf_im.value[:len(frequency_array)]) + + return hlm_fd, resized_td_modes + + +def compute_epoch_from_resized_td_modes(resized_td_modes, iota, phase, delta_t): + """Compute the epoch (time of peak) from resized TD modes at a given phase. + + Reconstructs the real h+ polarization from the resized TD mode data using + spin-weighted spherical harmonics at the specified phase, then finds the + peak via np.argmax. This matches the epoch convention used by gwsignal's + resize_gwpy_timeseries (which calls np.argmax on the real hp). + + Parameters + ---------- + resized_td_modes: dict + Dictionary with (l,m) keys and numpy complex arrays of resized TD modes. + iota: float + Inclination angle in radians. + phase: float + Reference phase for spherical harmonic evaluation. + delta_t: float + Time step of the TD data. + + Returns + ------- + epoch: float + The epoch value (negative of peak index times delta_t). + """ + chirplen = len(next(iter(resized_td_modes.values()))) + hp = np.zeros(chirplen) + for (l, m), h_data in resized_td_modes.items(): + ylm = lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, m) + hp += (ylm * h_data).real + return -np.argmax(hp) * delta_t + + +def apply_time_shift_to_fd_modes(hlm_fd_raw, resized_td_modes, iota, phase, domain): + """Apply a phase-dependent time shift to raw (unshifted) FD modes. + + Computes the epoch by reconstructing hp from resized TD modes at the given + phase (matching gwsignal's resize_gwpy_timeseries convention), then applies + the time shift exp(-i 2pi dt f) to each FD mode. + + Also returns a deferred_timeshift_data dict that can be passed to + sum_contributions_m to correct the time shift when a phase_shift is applied. + + Parameters + ---------- + hlm_fd_raw : dict + Dictionary with (l,m) keys and numpy arrays of unshifted FD modes. + resized_td_modes : dict + Dictionary with (l,m) keys and numpy complex arrays of resized TD modes. + iota : float + Inclination angle in radians. + phase : float + Reference phase for epoch computation. + domain : dingo.gw.domains.UniformFrequencyDomain + Frequency domain. + + Returns + ------- + hlm_fd : dict + Dictionary with (l,m) keys and time-shifted FD mode arrays. + deferred_timeshift_data : dict + Data needed by sum_contributions_m for deferred time shift correction. + """ + delta_t = 0.5 / domain.f_max + delta_f = domain.delta_f + frequency_array = domain() + + epoch = compute_epoch_from_resized_td_modes( + resized_td_modes, iota, phase, delta_t + ) + dt = 1.0 / delta_f + epoch + time_shift = np.exp(-1j * 2 * np.pi * dt * frequency_array) + + hlm_fd = {lm: hf * time_shift for lm, hf in hlm_fd_raw.items()} + + deferred_timeshift_data = { + "resized_td_modes": resized_td_modes, + "iota": iota, + "phase_ref": phase, + "dt_ref": dt, + "delta_t": delta_t, + "delta_f": delta_f, + "frequency_array": frequency_array, + } + + return hlm_fd, deferred_timeshift_data - return hlm_fd def get_polarizations_from_fd_modes_m(hlm_fd, iota, phase): pol_m = {} @@ -603,283 +635,3 @@ def taper_td_modes_for_SEOBRNRv5_extra_time( return h_return -def get_conditioning_params_for_TEOB(parameters_gwsignal): - """ - Compute conditioning parameters for TEOBResumSDALI waveforms. - - Replicates the parameter computation from gwsignal's - generate_conditioned_td_waveform_from_td - (lalsimulation/gwsignal/core/waveform_conditioning.py:34-111). - - This is needed because gwsignal's GenerateTDModes - (lalsimulation/gwsignal/core/waveform.py:668) returns unconditioned modes, - while GenerateFDWaveform routes through the conditioning pipeline. To make - mode-decomposed waveforms match the reference, the same conditioning must be - applied to each mode individually. - - Parameters - ---------- - parameters_gwsignal: dict - Dictionary of gwsignal parameters (from _convert_parameters). - - Returns - ------- - f_start: float - Lower starting frequency for generating modes with buffer for tapering. - t_extra: float - Duration of the extra time to taper at the beginning (for stage 1). - f_min: float - Effective minimum frequency, possibly lowered to fisco_9 - (for stage 2 beginning taper). - original_f_min: float - The originally requested starting frequency (for stage 1 high-pass). - f_isco: float - ISCO frequency at 6M (for stage 2 end taper). - """ - # waveform_conditioning.py:51-52 - extra_time_fraction = 0.1 - extra_cycles = 3.0 - - # waveform_conditioning.py:55-61 - # Masses are stored as astropy Quantities in solar masses; .value gives - # solar masses, .si.value gives kg. LAL bound functions expect SI (kg). - f_min = parameters_gwsignal["f22_start"].value - m1_SI = parameters_gwsignal["mass1"].si.value - m2_SI = parameters_gwsignal["mass2"].si.value - s1z = parameters_gwsignal["spin1z"].value - s2z = parameters_gwsignal["spin2z"].value - original_f_min = f_min - - # waveform_conditioning.py:70-72 — clamp f_min to ISCO at 9^1.5 - fisco_9 = 1.0 / ( - np.power(9.0, 1.5) * np.pi * (m1_SI + m2_SI) * lal.MTSUN_SI / lal.MSUN_SI - ) - if f_min > fisco_9: - f_min = fisco_9 - - # waveform_conditioning.py:76 — upper bound on chirp time from f_min - tchirp = LS.SimInspiralChirpTimeBound(f_min, m1_SI, m2_SI, s1z, s2z) - - # waveform_conditioning.py:79 - s = LS.SimInspiralFinalBlackHoleSpinBound(s1z, s2z) - - # waveform_conditioning.py:82 — merger + ringdown time - tmerge = ( - LS.SimInspiralMergeTimeBound(m1_SI, m2_SI) - + LS.SimInspiralRingdownTimeBound(m1_SI + m2_SI, s) - ) - - # waveform_conditioning.py:87 — extra cycles near merger - textra = extra_cycles / f_min - - # waveform_conditioning.py:90 — lower starting frequency for tapering buffer - f_start = LS.SimInspiralChirpStartFrequencyBound( - (1.0 + extra_time_fraction) * tchirp + tmerge + textra, m1_SI, m2_SI - ) - - # waveform_conditioning.py:106 — t_extra for stage 1 cosine taper - t_extra = extra_time_fraction * tchirp + textra - - # waveform_conditioning.py:108 — ISCO at 6M for stage 2 end taper - f_isco = 1.0 / ( - np.power(6.0, 1.5) * np.pi * (m1_SI + m2_SI) * lal.MTSUN_SI / lal.MSUN_SI - ) - - return f_start, t_extra, f_min, original_f_min, f_isco - - -def _high_pass_complex(data, dt, f_min, attenuation=0.99, order=8): - """ - High-pass filter a complex time series using a Butterworth IIR filter. - - Replicates gwsignal's high_pass_time_series - (lalsimulation/gwsignal/core/conditioning_subroutines.py:10-46). - - Parameters - ---------- - data: np.ndarray (complex) - Time series data array. - dt: float - Sampling interval in seconds. - f_min: float - Minimum frequency for high-pass. - attenuation: float - Attenuation at the low-frequency cutoff (default 0.99). - order: int - Order of Butterworth filter (default 8). - - Returns - ------- - filtered: np.ndarray (complex) - High-pass filtered data. - """ - # conditioning_subroutines.py:33 — sampling frequency - fs = 1.0 / dt - - # conditioning_subroutines.py:37-39 — bilinear transform to compute cutoff - w1 = np.tan(np.pi * f_min * dt) - wc = w1 * (1.0 / attenuation**0.5 - 1) ** (1.0 / (2.0 * order)) - fc = fs * np.arctan(wc) / np.pi - - # conditioning_subroutines.py:42-43 — forward-backward Butterworth filter - sos = butter(order, fc, btype="highpass", output="sos", fs=fs) - filtered_re = sosfiltfilt(sos, data.real) - filtered_im = sosfiltfilt(sos, data.imag) - return filtered_re + 1j * filtered_im - - -def _condition_stage1_complex(data, dt, t_extra, f_min): - """ - Stage 1 conditioning: cosine taper at the beginning + high-pass filter. - - Replicates gwsignal's time_array_condition_stage1 - (lalsimulation/gwsignal/core/conditioning_subroutines.py:50-87). - - Parameters - ---------- - data: np.ndarray (complex) - Time series data array. - dt: float - Sampling interval in seconds. - t_extra: float - Duration of extra time at the beginning to taper. - f_min: float - Minimum frequency for high-pass filter. - - Returns - ------- - data: np.ndarray (complex) - Conditioned data. - """ - # conditioning_subroutines.py:71-77 — cosine (Hann) taper at beginning - Ntaper = int(np.round(t_extra / dt)) - if Ntaper > 0 and Ntaper < len(data): - taper_array = np.arange(Ntaper) - w = 0.5 - 0.5 * np.cos(taper_array * np.pi / Ntaper) - data[:Ntaper] *= w - - # conditioning_subroutines.py:80-81 — high-pass filter - data = _high_pass_complex(data, dt, f_min) - - # conditioning_subroutines.py:84-85 — trim trailing zeros - data = np.trim_zeros(data, trim='b') - - return data - - -def _condition_stage2_complex(data, dt, f_min, f_isco): - """ - Stage 2 conditioning: end taper (1 cycle at f_isco) + beginning taper - (1 cycle at f_min). - - Replicates gwsignal's time_array_condition_stage2 - (lalsimulation/gwsignal/core/conditioning_subroutines.py:90-144). - - Parameters - ---------- - data: np.ndarray (complex) - Time series data array. - dt: float - Sampling interval in seconds. - f_min: float - Minimum frequency for beginning taper. - f_isco: float - ISCO frequency for end taper. - - Returns - ------- - data: np.ndarray (complex) - Conditioned data. - """ - # conditioning_subroutines.py:111-114 - min_taper_samples = 4 - Nsize = len(data) - if Nsize < 2 * min_taper_samples: - return data - - # conditioning_subroutines.py:119-129 — end taper: 1 cycle at f_isco - ntaper_end = max(int(np.round(1.0 / (f_isco * dt))), min_taper_samples) - taper_array = np.arange(1, ntaper_end) - w_end = 0.5 - 0.5 * np.cos(taper_array * np.pi / ntaper_end) - data[Nsize - ntaper_end + 1 :] *= w_end[::-1] - - # conditioning_subroutines.py:133-142 — beginning taper: 1 cycle at f_min - ntaper_begin = max(int(np.round(1.0 / (f_min * dt))), min_taper_samples) - taper_array = np.arange(ntaper_begin) - w_begin = 0.5 - 0.5 * np.cos(taper_array * np.pi / ntaper_begin) - data[:ntaper_begin] *= w_begin - - return data - - -def condition_td_modes_for_TEOB_in_place(hlm_td, t_extra, f_min, original_f_min, f_isco): - """ - Apply gwsignal-style conditioning to TD modes in place. - - This replicates the conditioning that gwsignal's - generate_conditioned_td_waveform_from_td - (lalsimulation/gwsignal/core/waveform_conditioning.py:34-111) applies to the - combined TD polarizations, but applies it to individual complex modes. This - is necessary because GenerateTDModes (gwsignal/core/waveform.py:668) returns - unconditioned modes — it bypasses the conditioning pipeline entirely. - - The conditioning consists of two stages (following the C implementation at - XLALSimInspiralTDConditionStage1 / Stage2): - - Stage 1 (conditioning_subroutines.py:50-87): - - Cosine (Hann) taper over t_extra seconds at the beginning - - Order-8 Butterworth high-pass filter at original_f_min (attenuation 0.99) - - Stage 2 (conditioning_subroutines.py:90-144): - - End taper: 1 cycle at f_isco (ISCO at 6M) - - Beginning taper: 1 cycle at f_min - - Usage - ----- - The modes must be generated at the lower starting frequency f_start returned - by get_conditioning_params_for_TEOB(), so that there is sufficient buffer - for the tapering. - - Example:: - - f_start, t_extra, f_min, original_f_min, f_isco = ( - get_conditioning_params_for_TEOB(parameters_gwsignal) - ) - # Generate modes at f_start (not the original f_min) - hlm_td, iota = generate_TD_modes_at_f_start(parameters, f_start) - condition_td_modes_for_TEOB_in_place( - hlm_td, t_extra, f_min, original_f_min, f_isco - ) - hlm_fd = td_modes_to_fd_modes(hlm_td, domain) - - Parameters - ---------- - hlm_td: dict - Dictionary with (l,m) keys and complex LAL COMPLEX16TimeSeries objects. - t_extra: float - Duration of extra time at the beginning for stage 1 cosine taper. - From get_conditioning_params_for_TEOB(). - f_min: float - Effective minimum frequency (possibly lowered to fisco_9). - Used for stage 2 beginning taper. - original_f_min: float - The originally requested starting frequency. - Used for stage 1 high-pass filter. - f_isco: float - ISCO frequency at 6M. Used for stage 2 end taper. - """ - for lm in list(hlm_td.keys()): - h = hlm_td[lm] - dt = h.dt.value - data = h.value.copy() - - # waveform_conditioning.py:106 — stage 1 uses original_f_min for high-pass - data = _condition_stage1_complex(data, dt, t_extra, original_f_min) - - # waveform_conditioning.py:109 — stage 2 uses adjusted f_min and f_isco - data = _condition_stage2_complex(data, dt, f_min, f_isco) - - # Replace with new TimeSeries (length may have changed due to trim) - hlm_td[lm] = TimeSeries( - data, t0=h.t0, dt=h.dt, unit=h.unit - ) From 6a774830590390b5a77d0030b3cf6fe903fb8732 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 17 Feb 2026 09:14:12 +0100 Subject: [PATCH 40/48] Fix TEOB mode decomposition to use parameters["phase"] and update test Use parameters["phase"] instead of self.spin_conversion_phase for TEOB epoch computation and spherical harmonics, consistent with SEOBNRv5EHM. Pass deferred_timeshift_data in test for correct phase-dependent epoch. Co-Authored-By: Claude Opus 4.6 --- dingo/gw/waveform_generator/waveform_generator.py | 4 ++-- tests/gw/waveform_generator/test_wfg_m.py | 8 +++++++- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index 44275f3cc..9351850c5 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -1301,7 +1301,7 @@ def generate_hplus_hcross_m( hlm_fd, self._deferred_timeshift_data = ( wfg_utils.apply_time_shift_to_fd_modes( hlm_fd_raw, resized_td_modes, iota, - self.spin_conversion_phase, self.domain, + parameters["phase"], self.domain, ) ) @@ -1310,7 +1310,7 @@ def generate_hplus_hcross_m( # use the non-precessing symmetry h_{l,-m}(t) = (-1)^l h*_{lm}(t) # to directly compute h+ and hx from the one-sided modes. pol_m = wfg_utils.get_polarizations_from_onesided_fd_modes_m( - hlm_fd, iota, self.spin_conversion_phase + hlm_fd, iota, parameters["phase"] ) elif generator.domain == "freq": diff --git a/tests/gw/waveform_generator/test_wfg_m.py b/tests/gw/waveform_generator/test_wfg_m.py index a909b2fec..e98f25b59 100644 --- a/tests/gw/waveform_generator/test_wfg_m.py +++ b/tests/gw/waveform_generator/test_wfg_m.py @@ -171,7 +171,13 @@ def test_generate_hplus_hcross_m(intrinsic_prior, wfg, num_evaluations, toleranc phase_shift = np.random.uniform(high=2 * np.pi) pol_m = wfg.generate_hplus_hcross_m(p) - pol = sum_contributions_m(pol_m, phase_shift=phase_shift) + pol = sum_contributions_m( + pol_m, + phase_shift=phase_shift, + deferred_timeshift_data=getattr( + wfg, "_deferred_timeshift_data", None + ), + ) pol_ref = wfg.generate_hplus_hcross({**p, "phase": p["phase"] + phase_shift}) mismatches.append( From 15697904ae19f9ba6bfc53917c8e0847c64f337a Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Thu, 19 Feb 2026 12:02:56 +0100 Subject: [PATCH 41/48] updated pesummary to not require phi_jl --- dingo/gw/result.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/dingo/gw/result.py b/dingo/gw/result.py index 0dc3b86c9..052a85c1a 100644 --- a/dingo/gw/result.py +++ b/dingo/gw/result.py @@ -731,14 +731,16 @@ def get_pesummary_samples(self, num_processes=1): "waveform_generator" ].get("spin_conversion_phase") - # Redefine phase parameter to be consistent with Bilby. COMMENTED BECAUSE SLOW - samples = change_spin_conversion_phase( - samples, - self.f_ref, - spin_conversion_phase_old, - None, - num_processes=num_processes, - ) + # Redefine phase parameter to be consistent with Bilby. + # Only apply spin conversion when precessing spin parameters are present. + if "phi_jl" in samples.columns: + samples = change_spin_conversion_phase( + samples, + self.f_ref, + spin_conversion_phase_old, + None, + num_processes=num_processes, + ) self._pesummary_samples = samples From bf249d5029c23a521232235e6c1ae7f4ece3fa35 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Fri, 20 Mar 2026 11:38:35 +0100 Subject: [PATCH 42/48] allowed for more waveform failures --- .../gw/waveform_generator/waveform_generator.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index 9351850c5..fba4a5436 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -237,8 +237,12 @@ def generate_hplus_hcross( if not catch_waveform_errors: raise else: - EDOM = e.args[0] == "Internal function call failed: Input domain error" - if EDOM: + EDOM = ( + len(e.args) > 0 + and e.args[0] == "Internal function call failed: Input domain error" + ) + recoverable = EDOM or isinstance(e, RuntimeError) + if recoverable: warnings.warn( f"Evaluating the waveform failed with error: {e}\n" f"The parameters were {parameters_generator}\n" @@ -1163,14 +1167,17 @@ def generate_FD_waveform(self, parameters_gwsignal: Dict) -> Dict[str, np.ndarra frequency_array = self.domain() h_plus = np.zeros_like(frequency_array, dtype=complex) h_cross = np.zeros_like(frequency_array, dtype=complex) - # Ensure that length of wf agrees with length of domain. Enforce by truncating frequencies beyond f_max + # Ensure that length of wf agrees with length of domain. if len(hp) > len(frequency_array): warnings.warn( "GWSignal waveform longer than domain's `frequency_array`" f"({len(hp)} vs {len(frequency_array)}). Truncating gwsignal array." ) - h_plus = hp[: len(h_plus)].value - h_cross = hc[: len(h_cross)].value + h_plus[:] = hp[: len(h_plus)].value + h_cross[:] = hc[: len(h_cross)].value + elif len(hp) < len(frequency_array): + h_plus[: len(hp)] = hp.value + h_cross[: len(hc)] = hc.value else: h_plus = hp.value h_cross = hc.value From 5830821f20eec43bad2ffe51b8dd3188c2fc0c46 Mon Sep 17 00:00:00 2001 From: nihargupte-ph Date: Mon, 23 Mar 2026 10:18:59 +0100 Subject: [PATCH 43/48] Add comment explaining zero-padding for short TD waveforms --- dingo/gw/waveform_generator/waveform_generator.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/dingo/gw/waveform_generator/waveform_generator.py b/dingo/gw/waveform_generator/waveform_generator.py index fba4a5436..2a17aeca8 100644 --- a/dingo/gw/waveform_generator/waveform_generator.py +++ b/dingo/gw/waveform_generator/waveform_generator.py @@ -1176,6 +1176,11 @@ def generate_FD_waveform(self, parameters_gwsignal: Dict) -> Dict[str, np.ndarra h_plus[:] = hp[: len(h_plus)].value h_cross[:] = hc[: len(h_cross)].value elif len(hp) < len(frequency_array): + # For TD models (e.g. TEOBResumS), GenerateFDWaveform chooses the TD + # sampling rate based on the waveform's frequency content, not the + # requested f_max. High-mass systems merge at lower frequencies, so + # the FFT can have fewer bins than the domain. Zero-padding is correct + # since there is no signal power above that frequency. h_plus[: len(hp)] = hp.value h_cross[: len(hc)] = hc.value else: From 609ad08d2d497e056523c2c11fffa7c91a33f8fc Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 27 Apr 2026 15:04:15 +0200 Subject: [PATCH 44/48] added deffered phase reconstruction --- dingo/gw/likelihood.py | 64 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 63 insertions(+), 1 deletion(-) diff --git a/dingo/gw/likelihood.py b/dingo/gw/likelihood.py index 49b29941f..f15be844d 100644 --- a/dingo/gw/likelihood.py +++ b/dingo/gw/likelihood.py @@ -14,7 +14,7 @@ DecimateWaveformsAndASDS, create_mask_based_on_frequency_update, ) -from dingo.gw.waveform_generator import WaveformGenerator +from dingo.gw.waveform_generator import WaveformGenerator, wfg_utils from dingo.gw.domains import ( Domain, UniformFrequencyDomain, @@ -382,6 +382,18 @@ def _log_likelihood_phase_grid_mode_decomposed(self, theta, phases=None): pol_m = self.signal_m({**theta, "phase": 0}) pol_m = {k: pol["waveform"] for k, pol in pol_m.items()} + # For TEOB, generate_hplus_hcross_m applies a phase-dependent time shift + # bookkept in wfg._deferred_timeshift_data. The analytical precomputation + # below assumes mu(phi) = sum_m mu_m(0) * exp(-i m phi), which is only + # valid when that timeshift is absent. When it's present we fall back to + # a per-phase resum so the phase-grid likelihood stays consistent with + # the full direct likelihood. + dts = getattr(self.waveform_generator, "_deferred_timeshift_data", None) + if dts is not None: + return self._log_likelihood_phase_grid_mode_decomposed_with_dts( + pol_m, d, phases, dts + ) + # Step 2: Precompute complex inner products (mu, mu) and (d, mu) for the # individual modes m. min_idx = self.data_domain.min_idx @@ -475,6 +487,56 @@ def _log_likelihood_phase_grid_mode_decomposed(self, theta, phases=None): return log_likelihoods + def _log_likelihood_phase_grid_mode_decomposed_with_dts( + self, pol_m, d, phases, dts + ): + """Variant of _log_likelihood_phase_grid_mode_decomposed for approximants + (currently TEOB) whose mode decomposition carries a phase-dependent time + shift recorded in `dts` (= wfg._deferred_timeshift_data). + + We cannot use the analytical precomputation trick because the timeshift + depends on the argmax of the TD-resummed h+, which varies with phase. + Instead we resum modes per phase and apply the exp(-2 pi i f dt(phi)) + correction before computing inner products. The geocenter-level time + shift commutes with detector projection, so we can apply it directly to + the per-detector signals returned by signal_m. + """ + m_vals = sorted(pol_m.keys()) + ifos = list(d.keys()) + freqs = self.data_domain() + + log_likelihoods = np.empty(len(phases)) + for idx, phi in enumerate(phases): + # Analytical resum of per-detector mode contributions. + mu = { + ifo: sum( + pol_m[m][ifo] * np.exp(-1j * m * phi) for m in m_vals + ) + for ifo in ifos + } + # Deferred-timeshift correction (no-op when phi == 0.0). + if phi != 0.0: + target_phase = dts["phase_ref"] + phi + epoch_target = wfg_utils.compute_epoch_from_resized_td_modes( + dts["resized_td_modes"], dts["iota"], target_phase, + dts["delta_t"], + ) + dt_target = 1.0 / dts["delta_f"] + epoch_target + dt_correction = dt_target - dts["dt_ref"] + shift = np.exp(-2j * np.pi * dt_correction * freqs) + for ifo in ifos: + mu[ifo] = mu[ifo] * shift + + rho2opt = sum( + inner_product(mu_ifo, mu_ifo) for mu_ifo in mu.values() + ) + kappa2 = sum( + inner_product(d_ifo, mu_ifo) + for d_ifo, mu_ifo in zip(d.values(), mu.values()) + ) + log_likelihoods[idx] = self.log_Zn + kappa2 - 0.5 * rho2opt + return log_likelihoods + def _log_likelihood_phase_marginalized(self, theta): """ From fd6ef59d60530e00e494e5c080054762a721f04f Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Mon, 27 Apr 2026 15:18:49 +0200 Subject: [PATCH 45/48] Fix test_wfg_m.py CI collection failure Remove leftover debug raise that broke CI when EOBRun_module isn't installed. Restructure approximant detection so pyseobnr and EOBRun_module are checked independently and contribute their approximants additively. Co-Authored-By: Claude Opus 4.7 (1M context) --- tests/gw/waveform_generator/test_wfg_m.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/tests/gw/waveform_generator/test_wfg_m.py b/tests/gw/waveform_generator/test_wfg_m.py index e98f25b59..7dc900acf 100644 --- a/tests/gw/waveform_generator/test_wfg_m.py +++ b/tests/gw/waveform_generator/test_wfg_m.py @@ -151,16 +151,21 @@ def tolerances(approximant): return 1e-5, 1e-5 -# Uncomment to test only one approximant. +approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM"] + +try: + import pyseobnr # noqa: F401 + + approximant_list += ["SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM"] +except ImportError: + pass + try: - # import pyseobnr - import EOBRun_module + import EOBRun_module # noqa: F401 - # approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM", "SEOBNRv5PHM", "SEOBNRv5HM", "SEOBNRv5EHM", "TEOBResumSDALI"] - approximant_list = ["TEOBResumSDALI"] + approximant_list += ["TEOBResumSDALI"] except ImportError: - raise ValueError("TMP") - approximant_list = ["IMRPhenomXPHM", "SEOBNRv4PHM"] + pass @pytest.mark.parametrize("approximant", approximant_list) From 6545dccbd510d7d9d4584c2b25e95c5dda6633c8 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 28 Apr 2026 13:48:09 +0200 Subject: [PATCH 46/48] Disable debug plotting in test_generate_hplus_hcross_m The debug branch was enabled and tried to savefig to a hardcoded local path that does not exist on CI runners, causing the CI failure. Co-Authored-By: Claude Opus 4.7 (1M context) --- tests/gw/waveform_generator/test_wfg_m.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/gw/waveform_generator/test_wfg_m.py b/tests/gw/waveform_generator/test_wfg_m.py index 7dc900acf..ce3245430 100644 --- a/tests/gw/waveform_generator/test_wfg_m.py +++ b/tests/gw/waveform_generator/test_wfg_m.py @@ -197,7 +197,7 @@ def test_generate_hplus_hcross_m(intrinsic_prior, wfg, num_evaluations, toleranc ] ) - debug = True + debug = False if debug: maxval = max(mismatches[-1]) idx = mismatches[-1].index(maxval) @@ -212,7 +212,6 @@ def test_generate_hplus_hcross_m(intrinsic_prior, wfg, num_evaluations, toleranc plt.xscale("log") plt.xlim((5, 128)) plt.title(f"{p}, mismatch={maxval}") - plt.savefig("/work/nihargupte/src/dingo/tests/gw/mm.png") plt.show() mismatches = np.array(mismatches) From ee62fd4988a6e1f89ccc4a99cbc430b4f01f3f46 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 28 Apr 2026 14:45:30 +0200 Subject: [PATCH 47/48] Reuse gwsignal taper in TEOB mode conditioning Replace dingo's reimplementation of the sigmoid taper (_compute_sigmoid_taper_params, _compute_sigmoid_window) with a get_tapering_window_for_real_array helper that calls cond.taper_gwpy_timeseries on a copy and recovers the window via |tapered|/|original|, mirroring the LAL-side get_tapering_window_for_complex_time_series trick. TEOB mode-resum output is bit-identical to the previous implementation. Co-Authored-By: Claude Opus 4.7 (1M context) --- dingo/gw/waveform_generator/wfg_utils.py | 123 ++++++++--------------- 1 file changed, 42 insertions(+), 81 deletions(-) diff --git a/dingo/gw/waveform_generator/wfg_utils.py b/dingo/gw/waveform_generator/wfg_utils.py index 04ee92f9c..a57a03556 100644 --- a/dingo/gw/waveform_generator/wfg_utils.py +++ b/dingo/gw/waveform_generator/wfg_utils.py @@ -1,7 +1,6 @@ import numpy as np import lal import lalsimulation as LS -from scipy.signal import find_peaks from lalsimulation.gwsignal.core import conditioning_subroutines as cond from gwpy.timeseries import TimeSeries @@ -72,6 +71,36 @@ def taper_td_modes_in_place(hlm_td, tapering_flag: int = 1): window = get_tapering_window_for_complex_time_series(h, tapering_flag) h.data.data *= window +def get_tapering_window_for_real_array(signal, taper_kind="start"): + """Recover the multiplicative gwsignal taper window for a real signal. + + Mirrors `get_tapering_window_for_complex_time_series` (LAL path) for the + gwsignal/TEOB path: runs `cond.taper_gwpy_timeseries` on a copy and + derives the window as |tapered| / |original| with an eps regularizer. + Returning the window (instead of the tapered signal) lets the caller + apply it independently to the real and imaginary parts of complex modes. + + Parameters + ---------- + signal: np.ndarray + Real-valued 1D array to derive the window from. + taper_kind: str + Forwarded to `cond.taper_gwpy_timeseries`. Default "start". + + Returns + ------- + window: np.ndarray + Real-valued window of the same length as `signal`. All ones if the + signal is too short or all-zero (matching gwsignal's no-op behavior). + """ + ts = TimeSeries(signal.copy(), dt=1.0) + ts_tapered = cond.taper_gwpy_timeseries(ts, taper_kind) + abs_signal = np.abs(signal) + max_abs = np.max(abs_signal) if abs_signal.size else 0.0 + eps = 1e-20 * max_abs if max_abs > 0 else 1e-20 + return (np.abs(ts_tapered.value) + eps) / (abs_signal + eps) + + def taper_td_gwpy_modes_in_place(hlm_td, iota, phase=0.0): """ Apply sigmoid start taper to gwpy TimeSeries modes in place. @@ -81,11 +110,14 @@ def taper_td_gwpy_modes_in_place(hlm_td, iota, phase=0.0): See generate_conditioned_td_waveform_from_td_fallback in lalsimulation/gwsignal/core/waveform_conditioning.py:35-64. - The gwsignal reference applies taper_gwpy_timeseries to the combined - polarizations hp and hc, where find_peaks determines the taper length. - Since hp and hc can yield different taper lengths, but we need a single - real-valued window to apply to each complex mode, we compute both windows - and use the one with shorter taper (preserving more signal). + The gwsignal reference applies `taper_gwpy_timeseries` to the combined + polarizations hp and hc, which can yield different windows because + find_peaks finds different peaks. We therefore derive `W_hp` and `W_hc` + separately via `get_tapering_window_for_real_array` and apply `W_hp` to + Re(h_lm) and `W_hc` to Im(h_lm) for each mode. For non-precessing + systems the modes are purely real in the co-rotating frame, so Re(h_lm) + and Im(h_lm) map cleanly to the two polarizations, giving agreement with + the reference. Parameters ---------- @@ -100,9 +132,8 @@ def taper_td_gwpy_modes_in_place(hlm_td, iota, phase=0.0): """ max_len = max(len(h) for h in hlm_td.values()) - # Sum modes to get combined hp and hc (unconditioned). - # These are needed to compute taper windows that match the reference path, - # where taper_gwpy_timeseries uses find_peaks on the combined polarizations. + # Sum modes to get combined hp and hc (unconditioned). The reference + # path runs the taper on these summed polarizations. h_complex = np.zeros(max_len, dtype=complex) for (l, m), h in hlm_td.items(): data = h.value.copy() @@ -111,31 +142,8 @@ def taper_td_gwpy_modes_in_place(hlm_td, iota, phase=0.0): ylm = lal.SpinWeightedSphericalHarmonic(iota, np.pi / 2 - phase, -2, l, m) h_complex += ylm * data - hp_combined = h_complex.real - hc_combined = -h_complex.imag - - # Compute separate taper windows for hp and hc. The reference path - # (generate_conditioned_td_waveform_from_td_fallback) applies - # taper_gwpy_timeseries independently to hp and hc, which can yield - # different taper lengths because find_peaks finds different peaks. - # We apply W_hp to Re(h_lm) and W_hc to Im(h_lm) for each mode. - # For non-precessing systems the modes are purely real in the co-rotating - # frame, so Re(h_lm) and Im(h_lm) map cleanly to the two polarizations, - # giving exact agreement with the reference. - start_hp, n_hp = _compute_sigmoid_taper_params(hp_combined) - start_hc, n_hc = _compute_sigmoid_taper_params(hc_combined) - - if start_hp is None and start_hc is None: - return - - # Build windows, falling back to the other if one polarization is empty - if start_hp is None: - start_hp, n_hp = start_hc, n_hc - if start_hc is None: - start_hc, n_hc = start_hp, n_hp - - W_hp = _compute_sigmoid_window(max_len, start_hp, n_hp) - W_hc = _compute_sigmoid_window(max_len, start_hc, n_hc) + W_hp = get_tapering_window_for_real_array(h_complex.real) + W_hc = get_tapering_window_for_real_array(-h_complex.imag) for (l, m) in hlm_td: data = hlm_td[(l, m)].value.copy() @@ -148,53 +156,6 @@ def taper_td_gwpy_modes_in_place(hlm_td, iota, phase=0.0): dt=hlm_td[(l, m)].dt, ) -def _compute_sigmoid_taper_params(signal): - """Compute (start, n) for the sigmoid start taper. - - Replicates the peak-finding logic in taper_gwpy_timeseries for - taper_kind='start'. Returns (None, None) if the signal is empty - or too short. - """ - LALSIMULATION_RINGING_EXTENT = 19 - - start = -1 - for idx, val in enumerate(signal): - if val != 0: - start = idx - break - if start == -1: - return None, None - - end = -1 - for idx, val in enumerate(signal[::-1]): - if val != 0: - end = len(signal) - 1 - idx - break - - if (end - start) <= 1: - return None, None - - mid = int((start + end) / 2) - pks, _ = find_peaks(abs(signal[start + 1 : mid])) - pks = pks[pks > LALSIMULATION_RINGING_EXTENT] - - if len(pks) < 2: - n = mid - start - else: - n = pks[1] + 1 - - return start, n - -def _compute_sigmoid_window(length, start, n): - """Compute the sigmoid taper window matching taper_gwpy_timeseries.""" - window = np.ones(length) - window[start] = 0.0 - realI = np.arange(1, n - 1) - z = (n - 1.0) / realI + (n - 1.0) / (realI - (n - 1.0)) - sigma = 1.0 / (np.exp(z) + 1.0) - window[start + 1 : start + n - 1] = sigma - return window - def td_modes_to_fd_modes(hlm_td, domain): """ Transform dict of td modes to dict of fd modes via FFT. The td modes are expected From d5f0b7719e000b3653941728e2bf9a1125f19f6c Mon Sep 17 00:00:00 2001 From: nihargupte-ph Date: Wed, 6 May 2026 04:10:58 -0700 Subject: [PATCH 48/48] rename chi_1/chi_2 -> spin_1z/spin_2z for pesummary PESummary's standard_names dictionary maps spin_1z/spin_2z but not chi_1/chi_2, so without this rename the spin conversion path is skipped and a_1, chi_eff, chi_p, final_spin etc. collapse to zero in the generated webpages. Apply the rename in get_pesummary_samples and the matching pesummary_prior so derived spin quantities are computed from the actual posterior. Co-Authored-By: Claude Opus 4.7 (1M context) --- dingo/gw/result.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/dingo/gw/result.py b/dingo/gw/result.py index 052a85c1a..9a6d658ef 100644 --- a/dingo/gw/result.py +++ b/dingo/gw/result.py @@ -701,6 +701,10 @@ def get_pesummary_samples(self, num_processes=1): * The spin angles phi_jl and theta_jn are transformed to account for a difference in phase definition. * Some columns are dropped: delta_log_prob_target, log_prob + * Aligned-spin columns ``chi_1``/``chi_2`` are renamed to + ``spin_1z``/``spin_2z`` so that PESummary's standard-name + dictionary recognises them and derives ``a_1``, ``a_2``, + ``chi_eff``, etc. """ if hasattr(self, "_pesummary_samples"): return self._pesummary_samples @@ -727,6 +731,13 @@ def get_pesummary_samples(self, num_processes=1): if "time" in col: samples.loc[:, col] += self.t_ref + # Rename aligned-spin columns to match PESummary's standard names. + # PESummary's standard_names.py maps spin_1z/spin_2z but not chi_1/chi_2, + # so without this rename pesummary cannot derive a_1, chi_eff, chi_p, etc. + samples.rename( + columns={"chi_1": "spin_1z", "chi_2": "spin_2z"}, inplace=True + ) + spin_conversion_phase_old = self.base_metadata["dataset_settings"][ "waveform_generator" ].get("spin_conversion_phase") @@ -752,7 +763,9 @@ def pesummary_prior(self): By convention, Dingo stores all times *relative* to a reference time, typically the trigger time for an event. The prior returned here corrects for that offset to - be consistent with other codes. + be consistent with other codes. Aligned-spin keys ``chi_1``/``chi_2`` are + renamed to ``spin_1z``/``spin_2z`` to match the renaming applied in + :meth:`get_pesummary_samples`. """ prior = copy.deepcopy(self.prior) for p in prior: @@ -762,4 +775,7 @@ def pesummary_prior(self): prior[p].minimum += self.t_ref except AttributeError: continue + for old, new in (("chi_1", "spin_1z"), ("chi_2", "spin_2z")): + if old in prior and new not in prior: + prior[new] = prior.pop(old) return prior