diff --git a/pycbc/inference/models/__init__.py b/pycbc/inference/models/__init__.py index 05f6c2b8766..ae0eba1a4e1 100644 --- a/pycbc/inference/models/__init__.py +++ b/pycbc/inference/models/__init__.py @@ -34,7 +34,8 @@ from .marginalized_gaussian_noise import MarginalizedTime from .brute_marg import BruteParallelGaussianMarginalize from .brute_marg import BruteLISASkyModesMarginalize -from .gated_gaussian_noise import (GatedGaussianNoise, GatedGaussianMargPol) +from .gated_gaussian_noise import (GatedGaussianNoise, GatedGaussianMargPol, + GatedGaussianMargPhase) from .single_template import SingleTemplate from .relbin import Relative, RelativeTime, RelativeTimeDom from .hierarchical import (HierarchicalModel, MultiSignalModel, @@ -202,6 +203,7 @@ def read_from_config(cp, **kwargs): BruteLISASkyModesMarginalize, GatedGaussianNoise, GatedGaussianMargPol, + GatedGaussianMargPhase, SingleTemplate, Relative, RelativeTime, diff --git a/pycbc/inference/models/gated_gaussian_noise.py b/pycbc/inference/models/gated_gaussian_noise.py index 7a09b9ce9c4..e08e72d5733 100644 --- a/pycbc/inference/models/gated_gaussian_noise.py +++ b/pycbc/inference/models/gated_gaussian_noise.py @@ -53,7 +53,7 @@ def __init__(self, variable_params, data, low_frequency_cutoff, psds=None, self._overwhitened_data = {} # cache the current gated data self._gated_data = {} - # cache terms related to normalization and gating + # cache terms related to normalization and gating self._invasds = {} self._Rss = {} self._lognorm = {} @@ -77,7 +77,7 @@ def __init__(self, variable_params, data, low_frequency_cutoff, psds=None, def from_config(cls, cp, data_section='data', data=None, psds=None, **kwargs): """Adds addiotional keyword arguments based on config file. - + Additional keyword arguments are: * ``highpass_waveforms`` : waveforms will be highpassed. @@ -142,10 +142,10 @@ def psds(self, psds): if self.normalize: self._set_covfit(det) self._overwhitened_data = self.whiten(self.data, 2, inplace=False) - + def _set_covfit(self, det): """Sets the fit function for estimating the covariance determinant. - + This must be called after the PSDs have been set, otherwise a ValueError will be raised. """ @@ -163,7 +163,7 @@ def _set_covfit(self, det): def logdet_fit(self, cov, p): """Construct a linear regression from a sample of truncated covariance matrices. - + Returns the sample points used for linear fit generation as well as the linear fit parameters. """ @@ -197,7 +197,7 @@ def logdet_fit(self, cov, p): x = numpy.vstack([sample_sizes, numpy.ones(len(sample_sizes))]).T m, b = numpy.linalg.lstsq(x, sample_dets, rcond=None)[0] return (sample_sizes, sample_dets), (m, b) - + @BaseGaussianNoise.normalize.setter def normalize(self, normalize): """Clears the current stats if the normalization state is changed. @@ -226,11 +226,11 @@ def gate_indices(self, det): gt = gate_start + window lindex, rindex = ts.get_gate_indices(gt, window) return lindex, rindex - + def det_lognorm(self, det, start_index=None, end_index=None): """Calculate the normalization term from the truncated covariance matrix. - + Determinant is estimated using a linear fit to logdet vs truncated matrix size. """ @@ -414,14 +414,14 @@ def _get_gate_times(self, gatestart, gateend, ra, dec): Parameters ---------- gatestart : float - Geocentric start time of the gate. + Start time of the gate. gateend : float - Geocentric end time of the gate. + End time of the gate. ra : float Right ascension of the signal. dec : float Declination of the signal. - + Returns ------- dict : @@ -432,10 +432,9 @@ def _get_gate_times(self, gatestart, gateend, ra, dec): thisdet = Detector(det) # account for the time delay between the waveforms of the # different detectors - gatestartdelay = gatestart + thisdet.time_delay_from_earth_center( - ra, dec, gatestart) - gateenddelay = gateend + thisdet.time_delay_from_earth_center( - ra, dec, gateend) + refdet = self.current_params.get('tc_ref_frame', 'geocentric') + gatestartdelay = thisdet.arrival_time(gatestart, ra, dec, refdet) + gateenddelay = thisdet.arrival_time(gateend, ra, dec, refdet) dgatedelay = gateenddelay - gatestartdelay gatetimes[det] = (gatestartdelay, dgatedelay) return gatetimes @@ -672,14 +671,31 @@ def _loglikelihood(self): rr = 4 * invpsd.delta_f * rtilde[slc].inner(gated_rtilde[slc]).real logl += norm - 0.5*rr return float(logl) - + + @property + def _extra_stats(self): + """Adds ``loglr``, plus ``cplx_loglr`` and ``optimal_snrsq`` in each + detector.""" + return ['loglr', 'maxl_phase'] + ['{}_optimal_snrsq'.format(det) for det in self._data] + + def _nowaveform_loglr(self): + """Convenience function to set loglr values if no waveform generated. + """ + setattr(self._current_stats, 'loglikelihood', -numpy.inf) + # maxl phase doesn't exist, so set it to nan + setattr(self._current_stats, 'maxl_phase', numpy.nan) + for det in self._data: + # snr can't be < 0 by definition, so return 0 + setattr(self._current_stats, '{}_optimal_snrsq'.format(det), 0.) + return -numpy.inf + @property def multi_signal_support(self): """ The list of classes that this model supports in a multi-signal likelihood """ return [type(self)] - + def multi_loglikelihood(self, models): """ Calculate a multi-model (signal) likelihood """ @@ -797,7 +813,7 @@ def get_gated_waveforms(self): # the waveforms are a dictionary of (hp, hc) pols = [] for h in wfs[det]: - ht = h.to_timeseries() + ht = h.to_timeseries() ht = ht.gate(gatestartdelay + dgatedelay/2, window=dgatedelay/2, copy=False, invpsd=invpsd, method='paint') @@ -805,7 +821,7 @@ def get_gated_waveforms(self): pols.append(h) out[det] = tuple(pols) return out - + def get_gate_times_hmeco(self): """Gets the time to apply a gate based on the current sky position. Returns @@ -931,14 +947,14 @@ def _loglikelihood(self): # compute the marginalized log likelihood marglogl = special.logsumexp(loglr) + lognl - numpy.log(len(self.pol)) return float(marglogl) - + @property def multi_signal_support(self): """ The list of classes that this model supports in a multi-signal likelihood """ return [type(self)] - + @catch_waveform_error def multi_loglikelihood(self, models): """ Calculate a multi-model (signal) likelihood @@ -965,3 +981,211 @@ def multi_loglikelihood(self, models): self._current_wfs = combine return self._loglikelihood() + + +class GatedGaussianMargPhase(BaseGatedGaussian): + r"""Gated Gaussian noise model that analytically marginalizes over the + phase of a signal. + + The phase to be marginalized over is specified by the user using the + `ref_phase` argument. If a model consists of multiple modes each with their + own phase, only the reference phase is marginalized over. All phases must + be specified with the `phase_names` argument. This can be passed as a list + or a string delimited by spaces (e.g. 'phase1 phase2 phase3'). + + Marginalization is done using explicit numerical integration over 500 + thousand integration points by default. This method assumes that the + waveform h can be written in terms of an overall phase phi as + + h = h_c * cos(phi) + h_s * sin(phi), + + where h_c and h_s are the waveform with phi set to zero and pi/2 + respectively. The number of integration points can be controlled via the + `phase_samples` argument. + """ + name = 'gated_gaussian_margphase' + + def __init__(self, variable_params, data, low_frequency_cutoff, psds=None, + high_frequency_cutoff=None, normalize=False, + static_params=None, phase_samples=500000, phase_names=None, + ref_phase=None, **kwargs): + # set up the boiler-plate attributes + super().__init__( + variable_params, data, low_frequency_cutoff, psds=psds, + high_frequency_cutoff=high_frequency_cutoff, normalize=normalize, + static_params=static_params, **kwargs) + self.det_names = list(self.data.keys()) + self.dets = {} + # phase marginalization parameters + self.phase_samples = int(phase_samples) + self.phases = numpy.linspace(0, 2*numpy.pi, self.phase_samples) + if ref_phase is None: + raise KeyError('ref_phase is set to None. Please specify the ' + 'name of the phase parameter to marginalize ' + 'over') + self.ref_phase = ref_phase + if phase_names is None: + logging.warning('No phase_names provided. Assuming single mode ' + f'specified by ref_phase {ref_phase}') + self.phase_names = [ref_phase] + elif type(phase_names) == list: + self.phase_names = phase_names + elif type(phase_names) == str: + self.phase_names = phase_names.split(' ') + else: + raise TypeError('Unrecognized format for phase_names arg. Accepts ' + 'string, list, or None') + # create the waveform generator + self.waveform_generator = create_waveform_generator( + self.variable_params, self.data, + waveform_transforms=self.waveform_transforms, + recalibration=self.recalibration, + generator_class=generator.FDomainDetFrameTwoPhaseGenerator, + **self.static_params) + + def get_waveforms(self): + r"""Generate the waveforms. + """ + if self._current_wfs is None: + params = self.current_params + # generate the cosine and sine terms + wfs = self.waveform_generator.generate(phases=self.phase_names, + ref_phase=self.ref_phase, + **params) + for det, (hc, hs) in wfs.items(): + # make the same length as the data + hc.resize(len(self.data[det])) + hs.resize(len(self.data[det])) + # apply high pass + if self.highpass_waveforms: + hc = highpass( + hc.to_timeseries(), + frequency=self.highpass_waveforms).to_frequencyseries() + hs = highpass( + hs.to_timeseries(), + frequency=self.highpass_waveforms).to_frequencyseries() + wfs[det] = (hc, hs) + self._current_wfs = wfs + return self._current_wfs + + def get_gated_waveforms(self): + r"""Generate the gated waveforms. + """ + wfs = self.get_waveforms() + out = {} + # apply the gate + for det, (hc, hs) in wfs.items(): + hct = hc.to_timeseries() + hst = hs.to_timeseries() + invpsd = self._invpsds[det] + gate_times = self.get_gate_times() + gatestartdelay, dgatedelay = gate_times[det] + hct = hct.gate(gatestartdelay + dgatedelay/2, + window=dgatedelay/2, copy=False, + invpsd=invpsd, method='paint') + hst = hst.gate(gatestartdelay + dgatedelay/2, + window=dgatedelay/2, copy=False, + invpsd=invpsd, method='paint') + hc = hct.to_frequencyseries() + hs = hst.to_frequencyseries() + out[det] = (hc, hs) + return out + + @property + def _extra_stats(self): + """Adds the maxL phase and corresponding likelihood.""" + return ['maxl_phase', 'maxl_logl'] + + @catch_waveform_error + def _loglikelihood(self): + r"""Computes the log likelihood. + """ + # get waveforms + wfs = self.get_waveforms() + gated_wfs = self.get_gated_waveforms() + # get data + data = self.get_data() + gated_data = self.get_gated_data() + # cycle over all detectors + norm = 0. + hchc = 0. + hchs = 0. + hshc = 0. + hshs = 0. + dhc = 0. + dhs = 0. + hcd = 0. + hsd = 0. + dd = 0. + for det in self.det_names: + if det not in self.dets: + self.dets[det] = Detector(det) + # we always filter the entire segment starting from kmin, since the + # gated series may have high frequency components + slc = slice(self._kmin[det], self._kmax[det]) + invpsd = self._invpsds[det] + d = data[det].copy() + hc, hs = wfs[det] + gated_d = gated_data[det].copy() + gated_hc, gated_hs = gated_wfs[det] + # overwhiten gated waveforms and data + gated_hc *= 2 * invpsd.delta_f * invpsd + gated_hs *= 2 * invpsd.delta_f * invpsd + gated_d *= 2 * invpsd.delta_f * invpsd + # evaluate the inner products + hchc += hc[slc].inner(gated_hc[slc]).real + hchs += hc[slc].inner(gated_hs[slc]).real + hshc += hs[slc].inner(gated_hc[slc]).real + hshs += hs[slc].inner(gated_hs[slc]).real + dhc += d[slc].inner(gated_hc[slc]).real + dhs += d[slc].inner(gated_hs[slc]).real + hcd += hc[slc].inner(gated_d[slc]).real + hsd += hs[slc].inner(gated_d[slc]).real + dd += d[slc].inner(gated_d[slc]).real + # get the normalization in this detector + if self.normalize: + start_index, end_index = self.gate_indices(det) + else: + start_index = end_index = None + norm += self.det_lognorm(det, start_index, end_index) + # numerical marginalization over phases + cphi = numpy.cos(self.phases) + sphi = numpy.sin(self.phases) + hh = cphi*cphi*hchc + sphi*sphi*hshs + cphi*sphi*(hchs+hshc) + dh = cphi*dhc + sphi*dhs + hd = cphi*hcd + sphi*hsd + loglr = -(hh-dh-hd) + lognl = -dd + # get the maxL phase + maxlidx = loglr.argmax() + setattr(self._current_stats, 'maxl_phase', self.phases[maxlidx]) + setattr(self._current_stats, 'maxl_logl', loglr[maxlidx] + lognl + norm) + # get the marginalized log likelihood ratio + marglogl = special.logsumexp(loglr) + lognl + norm - numpy.log(self.phase_samples) + return marglogl + + @property + def multi_signal_support(self): + """ The list of classes that this model supports in a multi-signal + likelihood + """ + return [type(self)] + + @catch_waveform_error + def multi_loglikelihood(self, models): + """ Calculate a multi-model (signal) likelihood + """ + # Generate the waveforms for each submodel + wfs = [] + for m in models + [self]: + wf = m.get_waveforms() + wfs.append(wf) + # combine into a single waveform + combine = {} + for det in self.data: + # get max waveform length + mlen = max([len(x[det]) for x in wfs]) + [x[det].resize(mlen) for x in wfs] + combine[det] = sum([x[det] for x in wfs]) + self._current_wfs = combine + return self._loglikelihood() diff --git a/pycbc/waveform/generator.py b/pycbc/waveform/generator.py index 435b70ab3be..b736c749a6e 100644 --- a/pycbc/waveform/generator.py +++ b/pycbc/waveform/generator.py @@ -42,6 +42,7 @@ from pycbc.pool import use_mpi import lal as _lal from pycbc import strain +from numpy import pi # utility functions/class @@ -965,6 +966,175 @@ def select_rframe_generator(approximant, domain): return select_waveform_generator(approximant, domain) +class FDomainDetFrameTwoPhaseGenerator(BaseFDomainDetFrameGenerator): + """Generates frequency-domain waveform in a specific frame. + + This class assumes that the radiation-frame waveform can be decomposed in + terms of a phase phi such that + + h = h_c * cos(phi) + h_s * sin(phi), + + where h_c and h_s are the waveform evaluated at phi = 0 and phi = pi/2 + respectively. The output waveforms have the full detector response applied. + + Parameters + ---------- + rFrameGeneratorClass : class + The class to use for generating the waveform in the radiation frame, + e.g., FDomainCBCGenerator. This should be the class, not an + instance of the class (the class will be initialized with the + appropriate arguments internally). + detectors : {None, list of strings} + The names of the detectors to use. If provided, all location parameters + must be included in either the variable args or the frozen params. If + None, the generate function will just return the plus polarization + returned by the rFrameGeneratorClass shifted by any desired time shift. + epoch : {float, lal.LIGOTimeGPS} + The epoch start time to set the waveform to. A time shift = tc - epoch is + applied to waveforms before returning. + variable_args : {(), list or tuple} + A list or tuple of strings giving the names and order of parameters + that will be passed to the generate function. + ref_phase : str + The phase used as a reference for generating h_c and h_s. If multiple + phases are in the signal, the relative difference between phases is + maintained during generation. For example, if a waveform has two phases + phi1 = pi/3 and phi2 = pi/2, with phi1 set as the reference, h_c will + be generated with phi1 = 0 and phi2 = pi/6, and h_s will be generated + with phi1 = pi/2 and phi2 = 2pi/3. + phases : {list, None}, optional + The names of the phase parameters taken in by the waveform approximant. + If None, it is assumed that the ref_phase argument is the only phase in + the signal. + \**frozen_params + Keyword arguments setting the parameters that will not be changed from + call-to-call of the generate function. + + Attributes + ---------- + detectors : dict + The dictionary of detectors that antenna patterns are calculated for + on each call of generate. If no detectors were provided, will be + ``{'RF': None}``, where "RF" means "radiation frame". + detector_names : list + The list of detector names. If no detectors were provided, then this + will be ['RF'] for "radiation frame". + epoch : lal.LIGOTimeGPS + The GPS start time of the frequency series returned by the generate function. + A time shift is applied to the waveform equal to tc-epoch. Update by using + ``set_epoch``. + current_params : dict + A dictionary of name, value pairs of the arguments that were last + used by the generate function. + rframe_generator : instance of rFrameGeneratorClass + The instance of the radiation-frame generator that is used for waveform + generation. All parameters in current_params except for the + location params are passed to this class's generate function. + frozen_location_args : dict + Any location parameters that were included in the frozen_params. + variable_args : tuple + The list of names of arguments that are passed to the generate + function. + """ + + location_args = set(['tc', 'ra', 'dec', 'polarization']) + """set(['tc', 'ra', 'dec', 'polarization']): + The set of location parameters. These are not passed to the rFrame + generator class; instead, they are used to apply the detector response + function and/or shift the waveform in time. The parameters are: + + * tc: The GPS time of coalescence. + * ra: Right ascension. + * dec: declination + * polarization: polarization. + * tc_ref_frame (optional): reference frame in which tc is defined. + Must be one of: 'geocentric', for geocentric time, or one of the + detector names. Default 'geocentric.' + + All of these must be provided in either the variable args or the + frozen params if detectors is not None. If detectors + is None, tc may optionally be provided. + """ + + def generate(self, phases=None, ref_phase=None, **kwargs): + """Generates a waveform, applies a time shift and the detector response + function from the given kwargs. + """ + self.current_params.update(kwargs) + rfparams = {param: self.current_params[param] + for param in kwargs if param not in self.location_args} + # generate the cosine term: ref_phase = 0 + if rfparams[ref_phase] != 0.: + raise ValueError(f'Reference phase {ref_phase}={rfparams[ref_phase]} is ' + 'not zero') + hpc, hcc = self.rframe_generator.generate(**rfparams) + # generate the sine term: shift all phases by pi/2 + sin_params = rfparams.copy() + for i in phases: + sin_params[i] = rfparams[i] + pi/2 + hps, hcs = self.rframe_generator.generate(**sin_params) + if isinstance(hpc, TimeSeries): + df = self.current_params['delta_f'] + hpc = hpc.to_frequencyseries(delta_f=df) + hcc = hcc.to_frequencyseries(delta_f=df) + hps = hps.to_frequencyseries(delta_f=df) + hcs = hcs.to_frequencyseries(delta_f=df) + # time-domain waveforms will not be shifted so that the peak amp + # happens at the end of the time series (as they are for f-domain), + # so we add an additional shift to account for it + tshift = 1./df - abs(hpc._epoch) + else: + tshift = 0. + hpc._epoch = hcc._epoch = hps._epoch = hcs._epoch = self._epoch + h = {} + if self.detector_names != ['RF']: + ra = self.current_params['ra'] + dec = self.current_params['dec'] + ref_tc = self.current_params['tc'] + pol = self.current_params['polarization'] + refframe = self.current_params.get('tc_ref_frame', 'geocentric') + for detname, det in self.detectors.items(): + tc = det.arrival_time(ref_tc, ra, dec, refframe) + # apply response function + fp, fc = det.antenna_pattern(ra, dec, pol, tc) + thishc = fp*hpc + fc*hcc + thishs = fp*hps + fc*hcs + # apply time shift + hc = apply_fd_time_shift(thishc, tc+tshift, copy=False) + hs = apply_fd_time_shift(thishs, tc+tshift, copy=False) + if self.recalib: + # recalibrate with given calibration model + hc = self.recalib[detname].map_to_adjust(hc, + **self.current_params) + hs = self.recalib[detname].map_to_adjust(hs, + **self.current_params) + h[detname] = (hc, hs) + else: + # no detector response, just use the + polarization + if 'tc' in self.current_params: + hpc = apply_fd_time_shift(hpc, self.current_params['tc']+tshift, + copy=False) + hps = apply_fd_time_shift(hps, self.current_params['tc']+tshift, + copy=False) + h['RF'] = (hpc, hps) + if self.gates is not None: + # resize all to nearest power of 2 + for ifo, (hc, hs) in h.items(): + hc.resize(ceilpow2(len(hc)-1) + 1) + hs.resize(ceilpow2(len(hs)-1) + 1) + # apply gates to wfs + h[ifo] = (strain.gate_data(hc, self.gates[ifo]), + strain.gate_data(hs, self.gates[ifo])) + return h + + @staticmethod + def select_rframe_generator(approximant, domain): + """Returns a radiation frame generator class based on the approximant + string. + """ + return select_waveform_generator(approximant, domain) + + class FDomainDetFrameModesGenerator(BaseFDomainDetFrameGenerator): r"""Generates frequency-domain waveform modes in a specific frame. @@ -1230,8 +1400,8 @@ def get_td_generator(approximant, modes=False): if approximant in supernovae.supernovae_td_approximants: return TDomainSupernovaeGenerator - raise ValueError(f"No time-domain generator found for " - "approximant: {approximant}") + raise ValueError(f"No time-domain generator found for " + "approximant: {approximant}") def get_fd_generator(approximant, modes=False): """Returns the frequency-domain generator for the given approximant.""" @@ -1246,7 +1416,7 @@ def get_fd_generator(approximant, modes=False): return FDomainFreqTauRingdownGenerator raise ValueError(f"No frequency-domain generator found for " - "approximant: {approximant}") + "approximant: {approximant}") def select_waveform_generator(approximant, domain=None): """Returns the single-IFO generator for the approximant. diff --git a/pycbc/waveform/ringdown.py b/pycbc/waveform/ringdown.py index 6f4feca6dab..f0cd0da3af9 100644 --- a/pycbc/waveform/ringdown.py +++ b/pycbc/waveform/ringdown.py @@ -716,7 +716,7 @@ def multimode_base(input_params, domain, freq_tau_approximant=False): be generated with :py:func:`td_damped_sinusoid` (:py:func:`fd_damped_sinusoid`). freq_tau_approximant : {False, bool}, optional - Choose choose the waveform approximant to use. Either based on + Choose the waveform approximant to use. Either based on mass/spin (set to False, default), or on frequencies/damping times of the modes (set to True). @@ -770,9 +770,10 @@ def multimode_base(input_params, domain, freq_tau_approximant=False): input_params['t_final']) sample_times = outplus.sample_times.numpy() elif domain == 'fd': - outplus, outcross = fd_output_vector(freqs, taus, - input_params['delta_f'], input_params['f_final']) kmin = int(input_params['f_lower'] / input_params['delta_f']) + outplus, outcross = fd_output_vector(freqs, taus, + input_params['delta_f'], + input_params['f_final']) sample_freqs = outplus.sample_frequencies.numpy()[kmin:] else: raise ValueError('unrecognised domain argument {}; ' @@ -921,7 +922,6 @@ def get_td_from_final_mass_spin(template=None, **kwargs): input_params = props(template, mass_spin_required_args, td_args, **kwargs) return multimode_base(input_params, domain='td') - def get_fd_from_final_mass_spin(template=None, **kwargs): """Return frequency domain ringdown with all the modes specified. @@ -1017,7 +1017,6 @@ def get_fd_from_final_mass_spin(template=None, **kwargs): input_params = props(template, mass_spin_required_args, fd_args, **kwargs) return multimode_base(input_params, domain='fd') - def get_td_from_freqtau(template=None, **kwargs): """Return time domain ringdown with all the modes specified. @@ -1121,7 +1120,6 @@ def get_td_from_freqtau(template=None, **kwargs): input_params = props(template, freqtau_required_args, td_args, **kwargs) return multimode_base(input_params, domain='td', freq_tau_approximant=True) - def get_fd_from_freqtau(template=None, **kwargs): """Return frequency domain ringdown with all the modes specified. @@ -1222,7 +1220,6 @@ def get_fd_from_freqtau(template=None, **kwargs): input_params = props(template, freqtau_required_args, fd_args, **kwargs) return multimode_base(input_params, domain='fd', freq_tau_approximant=True) - # Approximant names ########################################################### ringdown_fd_approximants = { 'FdQNMfromFinalMassSpin': get_fd_from_final_mass_spin,