Factored Likelihood

Code to compute the log likelihood of parameters of a gravitational waveform. Precomputes terms that depend only on intrinsic parameters and computes the log likelihood for given values of extrinsic parameters

Requires python SWIG bindings of the LIGO Algorithms Library (LAL)

RIFT.likelihood.factored_likelihood.AssertCalibrationSelfTermInnerProductSupported(IP_factory)

Raise ValueError if IP_factory cannot supply the calibration self-term cross terms.

Called from PrecomputeLikelihoodTerms BEFORE any waveform is generated, so a misconfigured job dies at once instead of producing a likelihood whose U/V cross terms and whose per-realization calibration self terms come from two DIFFERENT inner products.

RIFT.likelihood.factored_likelihood.BuildCalibrationSelfTermBasis(calibration_realizations, base_weights2side, rank_tol=1e-10, rank_max=None, verbose=False, use_cache=True)

Low-rank basis for the per-realization |C_c(f)|^2 profiles used by the fused- calmarg self-term fix. The calibration factors are random SPLINE draws made in-flight (per iteration / proposal round), so |C_c|^2 must be handled PER DRAW at runtime – but NOT by a per-draw band integral (that is what this avoids).

The SVD is the ONE-TIME “expand on a fixed basis” step; it depends only on the draws + PSD band (not the template), so it is memoized on a content fingerprint and reused across the intrinsic grid (a multi-point ILE job pays the SVD once, then only rank band integrals per point + a cheap per-draw linear combo).

|C_c(f)|^2 = amp_c(f)^2 with amp_c a cubic spline on ~n_spline (~10) log-f nodes, so the SET {|C_c|^2}_c spans a small fixed subspace (dim <= n_spline*(n_spline+1)/2, independent of n_cal). An SVD of the in-band |C_c(f)|^2 matrix recovers an orthonormal basis {b_k(f)} of that subspace and the per-draw coefficients:

|C_c(f)|^2 = sum_k alpha[k,c] b_k(f) (exact to rank_tol).

Then the per-realization weighted cross terms are a cheap linear combination of a FIXED, once-computed set of basis-weighted cross-term tensors (see CalibrationSelfTermCrossTermsFromBasis); no band integral is redone per draw, and a freshly drawn C_c only needs its coefficients alpha_c = B^+ |C_c|^2 (a matmul against the cached basis), not a new integral.

base_weights2sideComplexIP.weights2side (two-sided 1/S, zero out of band).

The band (nonzero-weight) support defines where |C_c|^2 must be represented; outside it the integrand h_lm* h_l’m’/S vanishes, so |C_c|^2 there is irrelevant.

Returns a dict with:

weights2side(rank, len2side) float64 – base_weights2side * b_k(f) (0 out of band),: ready to drop into ComplexIP.weights2side.

alpha : (rank, n_cal) float64 – per-draw expansion coefficients. rank, resid, n_cal.

RIFT.likelihood.factored_likelihood.CalibrationSelfTermCrossTermsFromBasis(IP, hlmsA, hlmsB, cal_basis, prefix='U', same_waveform_Q=False)

Per-realization |C_c|^2-weighted mode cross terms (fused-calmarg self-term fix), formed from the pre-built low-rank basis (BuildCalibrationSelfTermBasis) WITHOUT a per-draw band integral.

For each of the rank basis functions we compute the basis-weighted cross-term tensor U^(k)_{lm,l’m’} = 2 dF sum_f b_k(f)/S(f) h_lm*(f) h_l’m’(f) (rank band integrals, ONCE), then each realization’s cross terms are the cheap linear combo

U_c = sum_k alpha[k,c] U^(k).

Returns a LIST of n_cal cross-term dicts, keyed exactly like ComputeModeCrossTermIP.

Because |C_c(f)|^2 = sum_k alpha[k,c] b_k(f) exactly (to the basis rank_tol), this reproduces the per-realization <C_c h|C_c h> to that tolerance. alpha is real, so the U-hermitian / V-symmetric structure of each U^(k) (real weight) carries to U_c.

IP MUST be a frequency-domain-diagonal inner product (lsu.ComplexIP): the basis weight is applied by assigning IP.weights2side, which only a diagonal kernel reads. Callers reach this through PrecomputeLikelihoodTerms, which refuses any other inner product up front (CalibrationSelfTermInnerProductSupported).

RIFT.likelihood.factored_likelihood.CalibrationSelfTermInnerProductSupported(IP_factory): True if IP_factory may be used for the |C_c|^2-weighted calibration self-term cross terms. None (the default frequency-domain-diagonal lsu path) always is; any other factory must advertise supports_calibration_self_term = True.

RIFT.likelihood.factored_likelihood.ComplexAntennaFactor(det, RA, DEC, psi, tref)

Function to compute the complex-valued antenna pattern function: F+ + i Fx

‘det’ is a detector prefix string (e.g. ‘H1’) ‘RA’ and ‘DEC’ are right ascension and declination (in radians) ‘psi’ is the polarization angle ‘tref’ is the reference GPS time

RIFT.likelihood.factored_likelihood.ComputeArrivalTimeAtDetector(det, RA, DEC, tref)

Function to compute the time of arrival at a detector from the time of arrival at the geocenter.

‘det’ is a detector prefix string (e.g. ‘H1’) ‘RA’ and ‘DEC’ are right ascension and declination (in radians) ‘tref’ is the reference time at the geocenter. It can be either a float (in which case the return is a float) or a GPSTime object (in which case it returns a GPSTime)

RIFT.likelihood.factored_likelihood.ComputeArrivalTimeAtDetectorWithoutShift(det, RA, DEC, tref)

Function to compute the time of arrival at a detector from the time of arrival at the geocenter.

‘det’ is a detector prefix string (e.g. ‘H1’) ‘RA’ and ‘DEC’ are right ascension and declination (in radians) ‘tref’ is the reference time at the geocenter. It can be either a float (in which case the return is a float) or a GPSTime object (in which case it returns a GPSTime)

RIFT.likelihood.factored_likelihood.ComputeModeCrossTermIP(hlmsA, hlmsB, psd, fmin, fMax, fNyq, deltaF, analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0.0, verbose=True, prefix='U', same_waveform_Q=False, batched=None, IP_factory=None)

Compute the ‘cross terms’ between waveform modes, i.e. < h_lm | h_l’m’ >. The inner product is weighted by power spectral density ‘psd’ and integrated over the interval [-fNyq, -fmin] union [fmin, fNyq]

Returns a dictionary of inner product values keyed by tuples of mode indices i.e. ((l,m),(l’,m’))

same_waveform_Q: if true, it will do a loop that assumes (a) the same keys in both loops and (b) the same waveform in both (i.e., conjugate transpose): Note we WILL use the ‘prefix’ quantity to identify if we are computing ‘V’ and therefore using a different symmetry!
IP_factory: optional inner-product factory (see RIFT.likelihood.td_inner_product,: task T3 of TD_LIKELIHOOD.md). When None (the default), the inner product is exactly the current lsu.ComplexIP construction (byte-identical default path). When provided, IP_factory.ip_scalar(…) returns a drop-in object whose .ip(h1, h2) yields the same scalar <h1|h2>.

RIFT.likelihood.factored_likelihood.ComputeModeCrossTermIPCal(hlmsA, hlmsB, psd, fmin, fMax, fNyq, deltaF, calibration_realizations, analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0.0, verbose=False, prefix='U', same_waveform_Q=False, cal_basis=None)

Per-realization |C_c(f)|^2-weighted mode cross terms for the fused-calmarg self-term fix (rho_sq_c = <C_c h | C_c h>; the calmarg self-term-bias analysis note).

Uses the low-rank basis expansion (BuildCalibrationSelfTermBasis + CalibrationSelfTermCrossTermsFromBasis): rank (<= n_spline*(n_spline+1)/2, NOT n_cal) band integrals ONCE, then each realization’s U_c is a cheap linear combination – no per-draw band integral. Pass a pre-built cal_basis (dict) to share the SVD across the U and V calls; otherwise it is built here from the local IP weights.

Returns a LIST of n_cal cross-term dicts, keyed like ComputeModeCrossTermIP.

RIFT.likelihood.factored_likelihood.ComputeModeIPTimeSeries(hlms, data, psd, fmin, fMax, fNyq, N_shift, N_window, analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0.0, calibration_realizations=None, calibration_conjugate=False, IP_factory=None)

Compute the complex-valued overlap between each member of a SphHarmFrequencySeries ‘hlms’ and the interferometer data COMPLEX16FrequencySeries ‘data’, weighted the power spectral density REAL8FrequencySeries ‘psd’.

The integrand is non-zero in the range: [-fNyq, -fmin] union [fmin, fNyq]. This integrand is then inverse-FFT’d to get the inner product at a discrete series of time shifts.

Returns a SphHarmTimeSeries object containing the complex inner product for discrete values of the reference time tref. The epoch of the SphHarmTimeSeries object is set to account for the transformation

IP_factory: optional inner-product factory (see RIFT.likelihood.td_inner_product,: task T3 of TD_LIKELIHOOD.md). When None (the default), the overlap object is exactly the current lsu.ComplexOverlap construction (byte-identical default path). When provided, IP_factory.overlap(…) returns a drop-in object whose .ip(hlm, data) yields the same (rho, rhoTS, rhoIdx, rhoPhase) contract.

RIFT.likelihood.factored_likelihood.ComputeYlms(Lmax, theta, phi, selected_modes=None): Return a dictionary keyed by tuples (l,m) that contains the values of all -2Y_lm(theta,phi) with l <= Lmax -l <= m <= l

RIFT.likelihood.factored_likelihood.ComputeYlmsArray(lookupNK, theta, phi)

Returns an array Ylm[k] where lookup(k) = l,m. Only computes the LM values needed. theta, phi arguments are scalars

SHOULD BE DEPRECATED

RIFT.likelihood.factored_likelihood.ComputeYlmsArrayVector(lookupNK, theta, phi)

Returns an array Ylm[k] where lookup(k) = l,m. Only computes the LM values needed. theta, phi arguments are vectors. Shape is (len(th),len(lookup(NK)))

Should be combined with the previous routine ComputeYlmsArray (redundant)

Example:: th = np.linspace(0,np.pi, 5); lookupNK =[[2,-2], [2,2]] factored_likelihood.ComputeYlmsArrayVector(lookupNK,th,th)

RIFT.likelihood.factored_likelihood.DiscreteFactoredLogLikelihoodViaArray(tvals, P, lookupNKDict, rholmsArrayDict, ctUArrayDict, ctVArrayDict, epochDict, Lmax=2, array_output=False)

DiscreteFactoredLogLikelihoodViaArray uses the array-ized data structures to compute the log likelihood, either as an array vs time or marginalized in time. This generally is marginally faster, particularly if Lmax is large.

The timeseries quantities are computed via discrete shifts of an existing grid

Note ‘P’ must have the sampling rate set to correctly interpret the event time.: Note arguments passed are STILL SCALARS

RIFT.likelihood.factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVector(tvals, P_vec, lookupNKDict, rholmsArrayDict, ctUArrayDict, ctVArrayDict, epochDict, Lmax=2, array_output=False, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

DiscreteFactoredLogLikelihoodViaArray uses the array-ized data structures to compute the log likelihood, either as an array vs time or marginalized in time. This generally is marginally faster, particularly if Lmax is large.

The timeseries quantities are computed via discrete shifts of an existing grid

Note ‘P’ must have the sampling rate set to correctly interpret the event time.: Note arguments passed are NOW ARRAYS, in contrast to similar function which does not have ‘Vector’ postfix

RIFT.likelihood.factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals, P_vec, lookupNKDict, rholmsArrayDict, ctUArrayDict, ctVArrayDict, epochDict, Lmax=2, array_output=False, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>, loglikelihood=<function _factored_lnL_helper>, return_lnLt=False, phase_marginalization=False, n_cal=1, cal_method='loop', cal_distmarg=None, cal_log_weights=None, return_cal_components=False, time_interp='nearest', ctUArrayDict_cal=None, ctVArrayDict_cal=None, time_quadrature=None, explicit_time_values=False, return_time_draw=False, time_draw_uniforms=None, q_deltaT=None, td_subsample=1)

DiscreteFactoredLogLikelihoodViaArray uses the array-ized data structures to compute the log likelihood, either as an array vs time or marginalized in time. This generally is marginally faster, particularly if Lmax is large. The timeseries quantities are computed via discrete shifts of an existing grid Note ‘P’ must have the sampling rate set to correctly interpret the event time.

Note arguments passed are NOW ARRAYS, in contrast to similar function which does not have ‘Vector’ postfix

Calibration marginalization (n_cal>=1)

When n_cal>1, the rholm timeseries are assumed to hold n_cal contiguous calibration realizations, each of length N_window = npts_full/n_cal (built by ComputeModeIPTimeSeries with the calibration applied to the data). The data term kappa = <h|C_c d> changes per realization, selected by shifting the window offset ifirst -> ifirst + c*N_window.

The template self-term is NOT calibration-independent once the calibration is applied to the data: the statistically complete term is the per-realization rho_sq_c = <C_c h | C_c h> = <h| |C_c|^2 |h>. When the |C_c|^2-weighted cross terms are supplied (ctUArrayDict_cal/ctVArrayDict_cal), this per-realization rho_sq_c is used in place of the shared, cal-independent rho_sq=<h|h> – for every n_cal, INCLUDING n_cal==1 (a single applied draw still carries its own self-term). Omitting them (ctUArrayDict_cal=None) reproduces the historical, cal-independent-<h|h> behavior byte-for-byte. (The bare <h|h> shortcut broke the C->lambda*C distance-degeneracy invariance and biased lnZ; see ComputeModeCrossTermIPCal / BuildCalibrationSelfTermBasis.)

We then Monte-Carlo marginalize:: Z_cal(theta) = (1/n_cal) sum_c integral dt exp( lnL_t(theta, c) )

via a streaming log-sum-exp over the n_cal realizations (memory unchanged vs the n_cal==1 path; one extra GPU kernel launch per realization).

cal_method selects the n_cal>1 reduction:

‘loop’ (default, Option B): Python loop over realizations reusing the: existing Q_inner_product kernel; works on CPU and GPU and with any loglikelihood callback (distance/phase marginalization).
‘fused’ (Option C): a single fused kernel (RIFT.likelihood.Q_fused_calmarg): does the Q-product, the loglikelihood, and the cal+time log-sum-exp in one launch – a CUDA kernel on GPU, or pure numpy on CPU. Supports phase marginalization (|kappa|). With cal_distmarg=None it uses the default (distance-unmarginalized) helper kernel; with cal_distmarg=<table dict> it uses the separate distance-marginalization kernel (Q_fused_calmarg_distmarg) reproducing distmarg_loglikelihood. Raises NotImplementedError for unsupported combinations; the ‘loop’ path is the fallback for everything (incl. distmarg, on CPU and GPU).

cal_distmargdict or None

Distance-marginalization table+params for the fused distmarg kernel; see RIFT.likelihood.Q_fused_calmarg.Q_fused_calmarg_distmarg_cupy.

loglikelihoodcallable(kappa_sq, rho_sq) -> lnL(t)

MUST NOT WRITE INTO rho_sq. rho_sq is time-independent, so it is passed as a stride-0 broadcast_to view over an (npts_extrinsic,) vector rather than as a materialized (npts_extrinsic, npts) array. Every npts column therefore aliases one address: on CPU an in-place write raises ValueError: output array is read-only, but on GPU cupy’s broadcast is WRITABLE and an in-place write races, giving a wrong and irreproducible answer with no error. Every in-tree callback allocates (_factored_lnL_helper and the driver’s distmarg_loglikelihood), so nothing is broken today; a caller supplying its own must allocate too, or call xpy.ascontiguousarray first.

time_interp{‘nearest’, ‘cubic’, ‘sinc’}

Detector-time sampling convention for the data term. ‘nearest’ preserves the historical NoLoop integer-bin gather. ‘cubic’ evaluates the precomputed Q_lm time series at the fractional detector arrival time using a four-sample cubic Lagrange stencil, with zero extension outside the precomputed buffer.

CHOOSING BETWEEN ‘cubic’ AND ‘sinc’: neither is uniformly better. ‘cubic’ (4-point Lagrange) has O(h^4) error, so it improves fast with oversampling and is poor near Nyquist; ‘sinc’ (Lanczos) is window-limited, so its error is flat in oversampling. The operative quantity is NOT fNyq/fmax: Q^a_lm(t) is band-limited by whichever is lower, fmax or the TEMPLATE’s own cutoff, so the right choice depends on the masses, on fmin AND on srate. The recommendation is NOT restated here – it has been superseded twice and docstring copies went stale both times. Live value: RIFT.likelihood.time_interp_choice.CROSSOVER_GUIDANCE; measured tables: RIFT/likelihood/DESIGN_q_window_stencil.md. THE DEFAULT IS ‘nearest’, NOT ‘cubic’: this argument defaults to ‘nearest’, and the batch-mode CLI’s –interpolate-time defaults to off, which also resolves to ‘nearest’. Omitting either therefore keeps the historical nearest-bin behavior, whose errors the measured guidance calls scientifically significant (200-443 nats at SNR 100, reaching 1 nat by SNR 2-6); ‘cubic’ is only what a legacy truthy –interpolate-time value maps to. Ask for a stencil explicitly if you want one. All three stencils have CPU and GPU implementations. See _sinc_Q_window_numpy and RIFT/likelihood/DESIGN_q_window_stencil.md for the measured tables.

time_quadrature{‘simpson’, ‘trapezoid’, ‘bandlimited’, ‘peak-local’} or None

Rule used for the time integral. None (the default) defers to the module-level TIME_QUADRATURE_DEFAULT, which is ‘simpson’. Each name is the same rule on CPU and GPU (time_quadrature_rule).

‘trapezoid’ is the trapezoidal rule at the same fixed spacing. For a peak well inside the window it converges exponentially in sigma_t/deltaT and beats Simpson; it does not fix under-resolution.

‘simpson’ is the historical rule: Simpson’s rule at the FIXED spacing deltaT=1/srate. The integrand’s width, however, is set by the signal – sigma_t = 1/(2 pi rho sigma_f) – not by the data’s sample rate, so this under-resolves its own integrand at high SNR, and Simpson’s (4T_h-T_2h)/3 form makes an under-resolved peak worse than trapezoid rather than better.

‘bandlimited’ refines the grid to the integrand using ONLY the samples already computed: kappa(t) is band-limited below Nyquist by construction and rho_sq is time-independent on this path, so the continuous lnL(t) is recovered exactly by a zero-padded FFT per row. The refinement factor is DERIVED from the measured peak width and re-asserted on the refined grid; it is not a settable accuracy knob. Restricted to n_cal==1 and to the integrated outputs and continuous posterior draws (not return_lnLt / return_cal_components); anything else raises rather than quietly falling back. A time draw returns (time_offset, lnL_at_draw) and uses the same validated dense representation as the integral. Rationale, measured before/after and the exclusions: RIFT.likelihood.time_marginalization_quadrature.

‘peak-local’ is the same argument with the refined grid placed only where the integrand has support: kappa’s extrema are ENUMERATED on a small, SNR-independent upsample, an interval of a few sigma_t is built around each, overlapping intervals are MERGED into disjoint ones, and each is integrated at its own derived spacing. The mass left outside is bounded per row and checked, not assumed; a row whose bound is not small enough, or whose local grid would cost more than the dense one, is given the ‘bandlimited’ value. Same exclusions as ‘bandlimited’, PLUS phase marginalization, which it refuses. RIFT.likelihood.time_marginalization_peak_local.

Time-indexed cross terms (ctUArrayDict[det].ndim == 3)

For the time-domain fixed-segment code path, U and V depend on the arrival time and are supplied as (npts_full, n_lms, n_lms) arrays on the SAME time grid and epoch as rholmsArrayDict[det] (axis 0 of U <-> axis 1 of rholm). The self term is then evaluated per sample and per time at the same window index (ifirst + j, with the same time_interp handling) as the data term, and rho_sq becomes genuinely (npts_extrinsic, npts). Detectors may mix (some scalar, some time-indexed). All n_cal==1 reductions are supported (time marginalization, return_lnLt, and any elementwise loglikelihood such as the distance-marginalization table). Analytic orbital-phase marginalization is refused: the finite-segment norm depends on phase, even for the (2,+/-2) pair. Sample phase explicitly. The self term is gathered with the ‘nearest’ and ‘cubic’ stencils on the likelihood’s own time grid only: time_interp=’sinc’, an independently spaced Q pregrid (q_deltaT != deltaT), explicit_time_values, and the ‘bandlimited’ and ‘peak-local’ time quadratures (which assume a time-independent rho_sq) are refused with NotImplementedError when any detector is time-indexed.

Calibration marginalization with time-indexed terms (the td-gs code path) runs on cal_method=’loop’ only. The rholm buffers hold n_cal contiguous blocks as on the FD code path; a time-indexed base U(t) may then be ONE block long (N_window = npts_full/n_blocks; it is gathered block-locally) or the full length. The per-realization self term may be time-indexed: ctUArrayDict_cal[det] of shape (n_blocks, N_window, n_lms, n_lms) – U_c(t) of the template filtered with realization c, on the time grid of block c – gathered per realization at the same window index as its Q_c block, exactly as the time-indexed n_cal==1 self term (with td_subsample=K the block length is K*N_window and the gathers step by K). The calmarg reduction (streaming log-sum-exp over realizations, importance weights, return_lnLt, return_cal_components) is unchanged. cal_method=’fused’ and the SVD-compressed base U refuse these inputs with NotImplementedError.

SVD-compressed cross terms (ctUArrayDict[det] a TimeIndexedCrossTermsSVD)

compress_time_indexed_cross_terms(U3, V3, tol) replaces the pair (U(t), V(t)) by one object holding r component time series and their pair-space right vectors; pass it as ctUArrayDict[det] (ctVArrayDict[det] the same object or None). The self term is then r Q-style gathers per sample instead of n_lm (n_lm+1); exact up to the dropped singular components. Same restrictions as the 3-D arrays.

Sub-sample time grid (td_subsample = K > 1; td-gs –td-subsample K)

The rholm series (and the time-indexed U, V on the same grid) hold K entries per sample of P_vec.deltaT: entry r sits at epochDict[det] + r*deltaT/K (td_gs_precompute, “Sub-sample placement”). Each extrinsic sample’s first arrival time is then located on that fine grid (nearest fine entry, or the four-point cubic stencil on the fine grid for time_interp=’cubic’), and its time bins – still tvals, spaced deltaT – are gathered with stride K, so one gather per bin as before and the Simpson integral keeps dx = deltaT. The residual placement step of ‘nearest’ is deltaT/K. K = 1 is the historical code path, unchanged. With calibration marginalization the realization blocks of the rholms and the time-indexed per-realization self terms carry the same fine grid (K*N_window entries per block) and the loop reduction gathers them with the same stride; cal_method=’fused’ is refused with K > 1, as are the options refused for time-indexed cross terms above.

RIFT.likelihood.factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopOrig(tvals, P_vec, lookupNKDict, rholmsArrayDict, ctUArrayDict, ctVArrayDict, epochDict, Lmax=2, array_output=False): DiscreteFactoredLogLikelihoodViaArray uses the array-ized data structures to compute the log likelihood, either as an array vs time or marginalized in time. This generally is marginally faster, particularly if Lmax is large. The timeseries quantities are computed via discrete shifts of an existing grid Note ‘P’ must have the sampling rate set to correctly interpret the event time.

Note arguments passed are NOW ARRAYS, in contrast to similar function which does not have ‘Vector’ postfix

RIFT.likelihood.factored_likelihood.DiscreteSingleDetectorLogLikelihoodData(epoch, rholmsDictionary, tStart, nBins, RA, DEC, thS, phiS, psi, dist, Lmax, det)

Returns lnLdata array, evaluated at the geocenter, based on a DISCRETE timeshift.

At low sampling rates, this procedure will be considerably time offset (~ ms). It will also be undersampled for use in integration.

Return value is

a RAW numpy array
associated with nBins following tStart

Uses ‘tStart’ (a GPSTime) to identify the current detector orientations and hence time of flight delay.

the assumption is that nBins will be very small

Does NOT

resample the Q array : it is nearest-neighbor timeshifted before computing lnL

RIFT.likelihood.factored_likelihood.DiscreteSingleDetectorLogLikelihoodDataViaArray(tvals, extr_params, lookupNK, rholmsArrayDict, Lmax=2, det='H1'): SingleDetectorLogLikelihoodDataViaArray evaluates everything using arrays for each (l,m) pair. Uses discrete arrays. Compare to FactoredLogLikelihoodTimeMarginalized

RIFT.likelihood.factored_likelihood.FactoredLogLikelihood(extr_params, rholms, rholms_intp, crossTerms, crossTermsV, Lmax, interpolate=True)

Compute the log-likelihood = -1/2 < d - h | d - h > from:

extr_params is an object containing values of all extrinsic parameters
rholms_intp is a dictionary of interpolating functions < h_lm(t) | d >
crossTerms is a dictionary of < h_lm | h_l’m’ >
Lmax is the largest l-index of any h_lm mode considered

N.B. rholms_intp and crossTerms are the first two outputs of the function ‘PrecomputeLikelihoodTerms’

RIFT.likelihood.factored_likelihood.FactoredLogLikelihoodTimeMarginalized(tvals, extr_params, rholms_intp, rholms, crossTerms, crossTermsV, Lmax, interpolate=False)

Compute the log-likelihood = -1/2 < d - h | d - h > from:

extr_params is an object containing values of all extrinsic parameters
rholms_intp is a dictionary of interpolating functions < h_lm(t) | d >
crossTerms is a dictionary of < h_lm | h_l’m’ >
Lmax is the largest l-index of any h_lm mode considered

tvals is an array of timeshifts relative to the detector, used to compute the marginalized integral. It provides both the time prior and the sample points used for the integral.

N.B. rholms_intp and crossTerms are the first two outputs of the function ‘PrecomputeLikelihoodTerms’

RIFT.likelihood.factored_likelihood.FactoredLogLikelihoodViaArray(epoch, P, lookupNKDict, rholms_intpArrayDict, ctUArrayDict, ctVArrayDict)

FactoredLogLikelihoodViaArray uses the array-ized data structures to compute the log likelihood, a single scalar value. This generally is marginally faster, particularly if Lmax is large. The timeseries quantities are computed via interpolation onto the desired grid

Speed-wise, because we extract a single scalar value, this code has the same efficiency as FactoredLogLikelihoood

Note ‘P’ must have the sampling rate set to correctly interpret the event time. Note arguments passed are STILL SCALARS

RIFT.likelihood.factored_likelihood.IdentifyEffectiveModesForDetector(crossTermsOneDetector, fac, det)

RIFT.likelihood.factored_likelihood.InterpolateRholm(rholm, t, verbose=False)

RIFT.likelihood.factored_likelihood.InterpolateRholms(rholms, t, verbose=False)

Return a dictionary keyed on mode index tuples, (l,m) where each value is an interpolating function of the overlap against data as a function of time shift: rholm_intp(t) = < h_lm(t) | d >

‘rholms’ is a dictionary keyed on (l,m) containing discrete time series of < h_lm(t_i) | d > ‘t’ is an array of the discrete times: [t_0, t_1, …, t_N]

RIFT.likelihood.factored_likelihood.NetworkLogLikelihoodPolarizationMarginalized(epoch, rholmsDictionary, crossTerms, crossTermsV, tref, RA, DEC, thS, phiS, psi, dist, Lmax, detList): DOCUMENT ME!!!

RIFT.likelihood.factored_likelihood.NetworkLogLikelihoodTimeMarginalized(epoch, rholmsDictionary, crossTerms, crossTermsV, tref, RA, DEC, thS, phiS, psi, dist, Lmax, detList): DOCUMENT ME!!!

RIFT.likelihood.factored_likelihood.NetworkLogLikelihoodTimeMarginalizedDiscrete(epoch, rholmsDictionary, crossTerms, tref, deltaTWindow, RA, DEC, thS, phiS, psi, dist, Lmax, detList, array_output=False)

Uses DiscreteSingleDetectorLogLikelihoodData to calculate the lnL(t) on a discrete grid. - Mode 1: array_output = False

Computes int L dt/T over tref+deltaTWindow (GPSTime + pair, with |deltaTWindow|=T)

Mode 2: array_output = True Returns lnL(t) as a raw numpy array.

RIFT.likelihood.factored_likelihood.PackCalCrossTermsAsArrays(pairKeys, lookupKeysToNumber, crossTermsCalList, crossTermsCalListV)

Pack the per-realization |C_c|^2-weighted cross terms (fused-calmarg self-term fix) into arrays, using the SAME (l,m)->index mapping (lookupKeysToNumber) that PackLikelihoodDataStructuresAsArrays produced for the cal-independent U,V. This guarantees ctUArrayDict_cal[det][c] is index-aligned with ctUArrayDict[det].

crossTermsCalList, crossTermsCalListVlists of n_cal cross-term dicts (the: ComputeModeCrossTermIPCal output for one detector; or, on the td-gs code path, dicts of time-indexed 1-D arrays of length N_window, the per-realization U_c(t), V_c(t)).

Returns (U_cal, V_cal), each (n_cal, nKeys, nKeys) complex128 for scalar entries, or (n_cal, N_window, nKeys, nKeys) for time-indexed ones (axis 1 shares the time grid of the realization’s rholm block).

RIFT.likelihood.factored_likelihood.PackLikelihoodDataStructuresAsArrays(pairKeys, rholms_intpDictionaryForDetector, rholmsDictionaryForDetector, crossTermsForDetector, crossTermsForDetectorV)

Accepts list of LM pairs, dictionary for rholms against keys, and cross terms (a dictionary)

PROBLEM: Different detectors may have different time zeros. User must use the returned arrays with great care.

Cross-term values are either scalars (every frequency-domain code path) or 1-D arrays of length npts, the length of rholms[mode].data.data (the time-domain fixed-segment code path, where <h|h> depends on the arrival time because the template slides relative to a fixed data segment). Scalars pack to crossTermsArrayU/V of shape (nKeys, nKeys); time-indexed values pack to shape (npts, nKeys, nKeys), with index j along axis 0 sharing the time grid and epoch of rholmArray[:, j]. All pairs of one detector must be of one kind. With calibration marginalization the rholms hold n_cal contiguous realization blocks; the (calibration-independent) time-indexed cross terms may then be one block long (npts / n_cal, which must divide npts) and pack to (N_window, nKeys, nKeys).

RIFT.likelihood.factored_likelihood.PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict, psd_dict, Lmax, fMax, analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0.0, verbose=True, quiet=False, internal_fast_precompute=True, NR_group=None, NR_param=None, ignore_threshold=0.0001, extra_waveform_kwargs={}, use_gwsignal=False, use_gwsignal_approx=None, use_external_EOB=False, nr_lookup=False, nr_lookup_valid_groups=None, no_memory=True, perturbative_extraction=False, perturbative_extraction_full=False, hybrid_use=False, hybrid_method='taper_add', use_provided_strain=False, ROM_group=None, ROM_param=None, ROM_use_basis=False, ROM_limit_basis_size=None, skip_interpolation=False, force_22_mode=False, calibration_realizations=None, calibration_conjugate=False, return_calibration_crossterms=False, calibration_self_term=True, IP_factory=None, td_gs_contexts=None)

Compute < h_lm(t) | d > and < h_lm | h_l’m’ >

Returns:

Dictionary of interpolating functions, keyed on detector, then (l,m) e.g. rholms_intp[‘H1’][(2,2)]
Dictionary of “cross terms” <h_lm | h_l’m’ > keyed on (l,m),(l’,m’) e.g. crossTerms[((2,2),(2,1))]
Dictionary of discrete time series of < h_lm(t) | d >, keyed the same as the interpolating functions. Their main use is to validate the interpolating functions

IP_factory: optional inner-product factory threaded to ComputeModeCrossTermIP and

ComputeModeIPTimeSeries (task T3 of TD_LIKELIHOOD.md). Default None keeps the frequency-domain-diagonal lsu inner products (byte-identical default); a factory (e.g. RIFT.likelihood.td_inner_product.TDInnerProductFactory) swaps in the time-domain Toeplitz inner products at these two Precompute call sites only.

td_gs_contexts: optional {det: RIFT.likelihood.td_gs_precompute.TDGSDetectorContext}

(ILE –internal-ile-inner-product td-gs). When given, data_dict holds TIME-domain strain and the whole computation is handed to the time-domain code path (precompute_likelihood_terms_td_gs): TD modes with conditioning off, Q/U/V by direct quadratic forms against a Toeplitz covariance on the analysis pieces, no FFT of data or templates. On that code path crossTerms[det][(p1,p2)] and crossTermsV[det][(p1,p2)] are ARRAYS of length N_window (U(t), V(t) aligned with rholms), not scalars. Default None: this function is unchanged.

RIFT.likelihood.factored_likelihood.ReconstructPrecomputedLikelihoodTermsROM(P, acat_rom, rho_intp_rom, crossTerms_rom, crossTermsV_rom, rho_rom, verbose=True): Using a set of ROM coefficients for hlm[lm] = coef[l,m,basis] w[basis], reconstructs <h[lm]|data>, <h[lm]|h[l’m’]> Requires ROM also be loaded in top level, for simplicity

RIFT.likelihood.factored_likelihood.SingleDetectorLogLikelihood(rholm_vals, crossTerms, crossTermsV, Ylms, F, dist)

Compute the value of the log-likelihood at a single detector from several intermediate pieces of data.

Inputs:

rholm_vals: A dictionary of values of inner product between data
and h_lm modes, < h_lm(t*) | d >, at a single time of interest t*
crossTerms: A dictionary of inner products between h_lm modes:
< h_lm | h_l’m’ >
Ylms: Dictionary of values of -2-spin-weighted spherical harmonic modes
for a certain inclination and ref. phase, Y_lm(incl, - phiref)
F: Complex-valued antenna pattern depending on sky location and
polarization angle, F = F_+ + i F_x
dist: The distance from the source to detector in meters

Outputs: The value of ln L for a single detector given the inputs.

RIFT.likelihood.factored_likelihood.SingleDetectorLogLikelihoodData(epoch, rholmsDictionary, tref, RA, DEC, thS, phiS, psi, dist, Lmax, det): DOCUMENT ME!!!

RIFT.likelihood.factored_likelihood.SingleDetectorLogLikelihoodDataViaArray(epoch, lookupNK, rholms_intpArrayDict, tref, RA, DEC, thS, phiS, psi, dist, det)

SingleDetectorLogLikelihoodDataViaArray evaluates everything using arrays for each (l,m) pair Note arguments passed are STILL SCALARS

DEPRECATED: use DiscreteFactoredLogLikelihoodViaArray for end-to-end uuse USED IN : FactoredLogLikelihoodViaArray

RIFT.likelihood.factored_likelihood.SingleDetectorLogLikelihoodModel(crossTermsDictionary, crossTermsVDictionary, tref, RA, DEC, thS, phiS, psi, dist, Lmax, det): DOCUMENT ME!!!

RIFT.likelihood.factored_likelihood.SingleDetectorLogLikelihoodModelViaArray(lookupNKDict, ctUArrayDict, ctVArrayDict, tref, RA, DEC, thS, phiS, psi, dist, det): DOCUMENT ME!!!

RIFT.likelihood.factored_likelihood.TIME_QUADRATURE_DEFAULT = 'simpson': Time-marginalization quadrature used when a caller does not pass time_quadrature explicitly. ‘simpson’ is the historical fixed-deltaT Simpson rule and remains the DEFAULT; ‘bandlimited’ refines the time grid to the integrand actually present, using only the samples already computed (see RIFT.likelihood.time_marginalization_quadrature). Drivers set this once from their CLI so every call site inherits it; tests pass the kwarg directly.

class RIFT.likelihood.factored_likelihood.TimeIndexedCrossTermsSVD(c, right, iu, ju, n_lm, npts_full, singular_values, tol, rel_err_fro, rel_err_max)

Bases: object

Time-indexed cross terms U(t), V(t) compressed by an SVD over time.

Built by compress_time_indexed_cross_terms. Stack the pseudo-mode time series of _time_indexed_self_term_block as M (npts_full x n_cols), with n_cols = n_lm (n_lm+1), and take the thin SVD M = W diag(s) Z^H. Keeping the r components with s_k > tol * s_0,

M[t, col] = sum_k c_k(t) Zh[k, col], c_k(t) = W[t, k] s_k,

so the self term with the pseudo-mode weights w_e[col] of _time_indexed_self_term_weights becomes

rho_sq(e, t) = Re sum_col w_e[col] M[t_e, col]
= Re sum_k g_{e,k} c_k(t_e), g_{e,k} = sum_col w_e[col] Zh[k, col].

The weights are the same complex numbers as on the uncompressed path, so the only approximation is the dropped singular components (each bounded by tol * s_0 in spectral norm); with tol at rounding level the two agree to rounding. g is formed once per extrinsic sample (r contractions), and the time gather over the r series c_k uses the same Q_inner_product kernels / numpy windows as the data term, for nearest and cubic time_interp alike (cubic interpolation is linear in the series, so interpolating c_k is identical to interpolating M).

Attributes

c : (npts_full, r) complex128, C-contiguous – the component time series right : (n_cols, r) complex128, C-contiguous – Zh^T, so g = w @ right iu, ju : numpy int arrays, the unordered pair index (as the uncompressed path) n_lm, npts_full, n_cols, r : ints singular_values : numpy float64, all min(npts_full, n_cols) singular values tol : the relative cutoff used rel_err_fro : ||M - c right^T||_F / ||M||_F (0 when nothing is dropped) rel_err_max : max |M - c right^T| / max |M|

component_matrices(k, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

(U_k, V_k): upper-triangular n_lm x n_lm matrices such that, for any F, Y,: |F|^2 conj(Y)^T U_k Y + F^2 Y^T V_k Y = sum_col w[col] right[col, k] = g_k

exactly as complex numbers (the k-th contraction). The pseudo-mode of the pair (i<j) carries U_ij + conj(U_ji), so U_k is not Hermitian; the self term is Re sum_k g_k c_k(t), and g_k must be kept complex.

property ndim

property shape: The shape of the (uncompressed) time-indexed array this stands for.

to_xpy(xpy): Copy with c and right as xpy arrays (e.g. cupy); the index arrays and diagnostics stay numpy.

RIFT.likelihood.factored_likelihood.build_reflected_q_pregrid(rholms, factor=8, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

Build a one-time finer Q grid without changing the likelihood time grid.

The finite cut Q window is reflected before FFT interpolation so its unlike endpoints are never identified. Only the forward interval is retained; consequently the epoch is unchanged and every factor-th sample must reproduce the input. This helper is intentionally opt-in at the driver.

RIFT.likelihood.factored_likelihood.compress_time_indexed_cross_terms(U3, V3, tol, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

SVD-over-time compression of time-indexed cross terms; see TimeIndexedCrossTermsSVD for the formulation.

U3, V3 : (npts_full, n_lm, n_lm) arrays (numpy, or xpy arrays). tol : relative singular-value cutoff, 0 < tol < 1; components with

s_k > tol * s_0 are kept. 1e-13 keeps the result at rounding level for well-scaled U(t) (see docs/td_gs_route.md).

The SVD itself is done with xpy.linalg.svd (LAPACK on numpy, cuSOLVER on cupy). Returns a TimeIndexedCrossTermsSVD whose arrays live on xpy.

RIFT.likelihood.factored_likelihood.constructLMIterator(Lmax)

RIFT.likelihood.factored_likelihood.estimateEventTimeRelative(theEpochFiducial, rholms, rholms_intp)

RIFT.likelihood.factored_likelihood.estimateUpperDistanceBoundInMpc(rholms, crossTerms)

RIFT.likelihood.factored_likelihood.evaluateFast1dInterpolator(x, y, xlow, xhigh)

RIFT.likelihood.factored_likelihood.internal_hlm_generator(P, Lmax, verbose=True, quiet=False, NR_group=None, NR_param=None, extra_waveform_kwargs={}, use_gwsignal=False, use_gwsignal_approx=None, use_external_EOB=False, nr_lookup=False, nr_lookup_valid_groups=None, no_memory=True, perturbative_extraction=False, perturbative_extraction_full=False, hybrid_use=False, hybrid_method='taper_add', use_provided_strain=False, ROM_group=None, ROM_param=None, ROM_use_basis=False, ROM_limit_basis_size=None, skip_interpolation=False, force_22_mode=False, **kwargs): internal_hlm_generator: top-level front end to all waveform generators used. Needs to be restructured so it works on a ‘hook’ basis, so we are not constantly changing the source code

RIFT.likelihood.factored_likelihood.lalF(det, RA, DEC, psi, tref)

RIFT.likelihood.factored_likelihood.lalT(det, RA, DEC, tref)

RIFT.likelihood.factored_likelihood.makeFast1dInterpolator(y, xlow, xhigh)

RIFT.likelihood.factored_likelihood.marginalization_time_grid(integration_window_half, deltaT, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

The time-marginalization window grid, shared by every ILE likelihood path.

THE ONE CONSTRUCTOR. Both extrinsic drivers (bin/integrate_likelihood_extrinsic_batchmode and bin/integrate_likelihood_extrinsic_jax, the latter via RIFT.likelihood.jax_ile.wrapper) must obtain their grid here, or they silently evaluate different likelihoods – see issue #146, which measured up to 67.8 nats per sample between the two conventions that preceded this.

Convention:

npts  = int(2*integration_window_half/deltaT)
tvals = (arange(npts) - npts//2) * deltaT

Two properties, and both matter:

Spacing is EXACTLY deltaT. Every consumer of this grid (DiscreteFactoredLogLikelihoodViaArrayVector*, the JAX fused_log_likelihood*) reads only tvals[0] and len(tvals): it gathers rho[ifirst:ifirst+npts], i.e. steps by one sample, and integrates with dx=deltaT. So tvals[k] is a LABEL for a sample the code reaches by stepping deltaT from tvals[0], and the label is only truthful if the grid is deltaT-spaced. The former batchmode linspace(-iwh, iwh, npts) is spaced 2*iwh/(npts-1), which mislabelled its own samples by up to 1.4 samples (3.4e-4 s at iwh=0.075 s, srate=4096) at the window edge – visible wherever tvals is used as a time label, e.g. the time-resampling export in batchmode.
npts comes from ``int(2*iwh/deltaT)``, not ``2*int(iwh/deltaT)``. These differ whenever 2*iwh/deltaT has a fractional part below 0.5: at iwh=0.075 s they disagree at srate 1024, 2048 and 16384 (the low-mass production rate) and agree at 4096 and 8192. Taking the int(2*iwh/deltaT) form preserves batchmode’s window LENGTH at every rate; it lengthens the former JAX default by one sample at the rates above, which is the deliberate choice made here – a window that is a sample too SHORT truncates the time marginalization, and matching the longer-standing production length is the lower-risk direction.

With npts even the grid is [-npts/2, npts/2) samples, exactly reproducing the former JAX arange(-Nw, Nw)*deltaT; with npts odd it is symmetric, [-(npts//2), +(npts//2)]. Either way t=0 – the fiducial epoch – is on the grid exactly, which the old linspace only achieved for odd npts.

Parameters

integration_window_halffloat: Half-width of the marginalization window in seconds (the drivers’ t_ref_wind / --data-integration-window-half).
deltaTfloat: Sample spacing in seconds.
xpymodule: numpy or cupy; the array is built with xpy.arange.

Returns

array of shape (npts,), spacing exactly deltaT, containing 0.0.

RIFT.likelihood.factored_likelihood.non_herm_hoff(P)

RIFT.likelihood.factored_likelihood.prepare_reflected_q_pregrid(rholms_by_detector, factor=8, transfer=None, cleanup=None)

Transactionally build and optionally transfer a detector Q pregrid.

A backend OOM after one detector transfer cannot leave a mixed host/device, coarse/fine dictionary. Partial temporaries are dropped, cleanup is invoked (normally CuPy’s memory-pool release), and the original coarse Q dictionary is transferred instead. The caller can then restore its prior stencil and continue with an explicit fallback telemetry record.

RIFT.likelihood.factored_likelihood.profile(fn)

RIFT.likelihood.factored_likelihood.simpson_function(xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>): Time-integration rule for the xpy likelihoods. numpy uses scipy directly; other backends apply scipy’s weights, so CPU and GPU integrate with the same rule for odd and even npts (optimized_gpu_tools.simps handles an even npts differently from current scipy).

RIFT.likelihood.factored_likelihood.simpson_weights(npts, dx, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>): Weights w with my_simps(y, dx=dx, axis=-1) == y @ w, for any npts (scipy’s rule, including its treatment of an even number of samples). Computed once per (npts, dx) with scipy and moved to the array backend.

RIFT.likelihood.factored_likelihood.time_quadrature_rule(time_quadrature, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

The fixed-grid time rule for time_quadrature, the same on every backend.

‘trapezoid’ gives the trapezoidal rule. Every other name gives scipy’s Simpson rule, which is also what ‘bandlimited’ and ‘peak-local’ apply to the rows they do not refine. On numpy, Simpson is scipy.integrate.simpson itself, so CPU results are unchanged bit for bit. The callable takes (y, dx=..., axis=-1).

RIFT.likelihood.factored_likelihood.validate_time_interp(time_interp, on_gpu=False)

Reject unknown stencils loudly.

All three stencils now have both a CPU and a GPU implementation (‘sinc’ via the Q_inner_sinc kernel added alongside Q_inner and Q_inner_cubic), so on_gpu no longer restricts the choice. It is kept in the signature because the callers pass it and because it documents, at each call site, that the stencil has to be legal on the backend actually in use.

Overview

The factored_likelihood module provides tools to compute the log-likelihood of gravitational waveform parameters. It employs a “factored” approach where terms depending only on intrinsic parameters are precomputed, allowing for efficient evaluation over extrinsic parameters.

Main Driver Functions

RIFT.likelihood.factored_likelihood.PrecomputeLikelihoodTerms(event_time_geo, t_window, P, data_dict, psd_dict, Lmax, fMax, analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0.0, verbose=True, quiet=False, internal_fast_precompute=True, NR_group=None, NR_param=None, ignore_threshold=0.0001, extra_waveform_kwargs={}, use_gwsignal=False, use_gwsignal_approx=None, use_external_EOB=False, nr_lookup=False, nr_lookup_valid_groups=None, no_memory=True, perturbative_extraction=False, perturbative_extraction_full=False, hybrid_use=False, hybrid_method='taper_add', use_provided_strain=False, ROM_group=None, ROM_param=None, ROM_use_basis=False, ROM_limit_basis_size=None, skip_interpolation=False, force_22_mode=False, calibration_realizations=None, calibration_conjugate=False, return_calibration_crossterms=False, calibration_self_term=True, IP_factory=None, td_gs_contexts=None)

Compute < h_lm(t) | d > and < h_lm | h_l’m’ >

Returns:

Dictionary of interpolating functions, keyed on detector, then (l,m) e.g. rholms_intp[‘H1’][(2,2)]
Dictionary of “cross terms” <h_lm | h_l’m’ > keyed on (l,m),(l’,m’) e.g. crossTerms[((2,2),(2,1))]
Dictionary of discrete time series of < h_lm(t) | d >, keyed the same as the interpolating functions. Their main use is to validate the interpolating functions

IP_factory: optional inner-product factory threaded to ComputeModeCrossTermIP and

ComputeModeIPTimeSeries (task T3 of TD_LIKELIHOOD.md). Default None keeps the frequency-domain-diagonal lsu inner products (byte-identical default); a factory (e.g. RIFT.likelihood.td_inner_product.TDInnerProductFactory) swaps in the time-domain Toeplitz inner products at these two Precompute call sites only.

td_gs_contexts: optional {det: RIFT.likelihood.td_gs_precompute.TDGSDetectorContext}

(ILE –internal-ile-inner-product td-gs). When given, data_dict holds TIME-domain strain and the whole computation is handed to the time-domain code path (precompute_likelihood_terms_td_gs): TD modes with conditioning off, Q/U/V by direct quadratic forms against a Toeplitz covariance on the analysis pieces, no FFT of data or templates. On that code path crossTerms[det][(p1,p2)] and crossTermsV[det][(p1,p2)] are ARRAYS of length N_window (U(t), V(t) aligned with rholms), not scalars. Default None: this function is unchanged.

Computes the inner products between waveform modes and data ($langle h_{lm}(t) | d rangle$) and the cross terms between modes ($langle h_{lm} | h_{l’m’} rangle$).

RIFT.likelihood.factored_likelihood.FactoredLogLikelihood(extr_params, rholms, rholms_intp, crossTerms, crossTermsV, Lmax, interpolate=True)

Compute the log-likelihood = -1/2 < d - h | d - h > from:

extr_params is an object containing values of all extrinsic parameters
rholms_intp is a dictionary of interpolating functions < h_lm(t) | d >
crossTerms is a dictionary of < h_lm | h_l’m’ >
Lmax is the largest l-index of any h_lm mode considered

N.B. rholms_intp and crossTerms are the first two outputs of the function ‘PrecomputeLikelihoodTerms’

Computes the log-likelihood for a given set of extrinsic parameters using precomputed terms.

RIFT.likelihood.factored_likelihood.FactoredLogLikelihoodTimeMarginalized(tvals, extr_params, rholms_intp, rholms, crossTerms, crossTermsV, Lmax, interpolate=False)

Compute the log-likelihood = -1/2 < d - h | d - h > from:

extr_params is an object containing values of all extrinsic parameters
rholms_intp is a dictionary of interpolating functions < h_lm(t) | d >
crossTerms is a dictionary of < h_lm | h_l’m’ >
Lmax is the largest l-index of any h_lm mode considered

tvals is an array of timeshifts relative to the detector, used to compute the marginalized integral. It provides both the time prior and the sample points used for the integral.

N.B. rholms_intp and crossTerms are the first two outputs of the function ‘PrecomputeLikelihoodTerms’

Computes the log-likelihood marginalized over the arrival time.

Helper Functions

RIFT.likelihood.factored_likelihood.ComputeModeIPTimeSeries(hlms, data, psd, fmin, fMax, fNyq, N_shift, N_window, analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0.0, calibration_realizations=None, calibration_conjugate=False, IP_factory=None)

Compute the complex-valued overlap between each member of a SphHarmFrequencySeries ‘hlms’ and the interferometer data COMPLEX16FrequencySeries ‘data’, weighted the power spectral density REAL8FrequencySeries ‘psd’.

The integrand is non-zero in the range: [-fNyq, -fmin] union [fmin, fNyq]. This integrand is then inverse-FFT’d to get the inner product at a discrete series of time shifts.

Returns a SphHarmTimeSeries object containing the complex inner product for discrete values of the reference time tref. The epoch of the SphHarmTimeSeries object is set to account for the transformation

IP_factory: optional inner-product factory (see RIFT.likelihood.td_inner_product,: task T3 of TD_LIKELIHOOD.md). When None (the default), the overlap object is exactly the current lsu.ComplexOverlap construction (byte-identical default path). When provided, IP_factory.overlap(…) returns a drop-in object whose .ip(hlm, data) yields the same (rho, rhoTS, rhoIdx, rhoPhase) contract.

Computes the complex-valued overlap between each waveform mode and the interferometer data, weighted by the power spectral density (PSD).

RIFT.likelihood.factored_likelihood.InterpolateRholms(rholms, t, verbose=False)

Return a dictionary keyed on mode index tuples, (l,m) where each value is an interpolating function of the overlap against data as a function of time shift: rholm_intp(t) = < h_lm(t) | d >

‘rholms’ is a dictionary keyed on (l,m) containing discrete time series of < h_lm(t_i) | d > ‘t’ is an array of the discrete times: [t_0, t_1, …, t_N]

Returns a dictionary of interpolating functions for the overlap $langle h_{lm}(t) | d rangle$.

RIFT.likelihood.factored_likelihood.ComputeModeCrossTermIP(hlmsA, hlmsB, psd, fmin, fMax, fNyq, deltaF, analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0.0, verbose=True, prefix='U', same_waveform_Q=False, batched=None, IP_factory=None)

Compute the ‘cross terms’ between waveform modes, i.e. < h_lm | h_l’m’ >. The inner product is weighted by power spectral density ‘psd’ and integrated over the interval [-fNyq, -fmin] union [fmin, fNyq]

Returns a dictionary of inner product values keyed by tuples of mode indices i.e. ((l,m),(l’,m’))

same_waveform_Q: if true, it will do a loop that assumes (a) the same keys in both loops and (b) the same waveform in both (i.e., conjugate transpose): Note we WILL use the ‘prefix’ quantity to identify if we are computing ‘V’ and therefore using a different symmetry!
IP_factory: optional inner-product factory (see RIFT.likelihood.td_inner_product,: task T3 of TD_LIKELIHOOD.md). When None (the default), the inner product is exactly the current lsu.ComplexIP construction (byte-identical default path). When provided, IP_factory.ip_scalar(…) returns a drop-in object whose .ip(h1, h2) yields the same scalar <h1|h2>.

Computes the inner products (cross terms) between waveform modes $langle h_{lm} | h_{l’m’} rangle$.

RIFT.likelihood.factored_likelihood.ComplexAntennaFactor(det, RA, DEC, psi, tref)

Function to compute the complex-valued antenna pattern function: F+ + i Fx

‘det’ is a detector prefix string (e.g. ‘H1’) ‘RA’ and ‘DEC’ are right ascension and declination (in radians) ‘psi’ is the polarization angle ‘tref’ is the reference GPS time

Computes the complex-valued antenna pattern function $F_+ + i F_times$.

RIFT.likelihood.factored_likelihood.ComputeYlms(Lmax, theta, phi, selected_modes=None)

Return a dictionary keyed by tuples (l,m) that contains the values of all -2Y_lm(theta,phi) with l <= Lmax -l <= m <= l

Computes the spin-weighted spherical harmonics $-2Y_{lm}(theta, phi)$.

RIFT.likelihood.factored_likelihood.ComputeArrivalTimeAtDetector(det, RA, DEC, tref)

Function to compute the time of arrival at a detector from the time of arrival at the geocenter.

‘det’ is a detector prefix string (e.g. ‘H1’) ‘RA’ and ‘DEC’ are right ascension and declination (in radians) ‘tref’ is the reference time at the geocenter. It can be either a float (in which case the return is a float) or a GPSTime object (in which case it returns a GPSTime)

Computes the time of arrival at a detector from the geocenter reference time.