diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/SLOWROT_HANDOFF.md b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/SLOWROT_HANDOFF.md new file mode 100644 index 000000000..fc6821058 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/SLOWROT_HANDOFF.md @@ -0,0 +1,207 @@ +# Slow-rotation RIFT likelihood — handoff / breadcrumb + +Branch: `rift_slowrot` (off `rift_O4d_fix_rvs_clear_fairdraw_batch`). Design notes + LaTeX +paper live in a **separate local repo** `~/rift-slow-rotation` (no git remote — local only; +21-page PDF via `latexmk -pdf notes/main.tex`). + +## PATH B DELAY PHYSICS — VALIDATED (2026-07 update) +Path B's propagation-delay drift is validated against an INDEPENDENT folded-template ground +truth (`test_slowrot_pathB_bruteforce.py`): data built as Re[F(t) Sigma(t-tau(t))] from the +SAME harmonic F(t),tau(t) model + the SAME modes the likelihood uses (no SimDetectorStrain +convention floor), at an inflated sidereal rate so the delay drift is large. At the REAL +90-min-BNS rate (= x340 inflation on a 16s test): p_max=0 deficit 3.43 -> p_max=1 0.23 -> +p_max=2 0.207 -> p_max=3 0.207: CONVERGES, bound-respected, NO blow-up. So Path B recovers +the delay drift and is production-ready for the target signals with p_max<=2. +KNOWN LIMIT: the p>=3 catastrophic cancellation (huge high-f U terms x tiny delta_tau^p +coefficients) only bites at x1000+ inflation (>2.6x faster than any physical signal): x1000 +gives p=2 deficit 5.9 but p=3 blows to 1e5. So the band-limit fix (low-pass the p>=1 +derivative templates) is a robustness nicety, NOT a blocker for real signals. +Separately validated vs LAL's SimDetectorStrainREAL8TimeSeries (`test_slowrot_pathB_groundtruth.py`): +baseline/PathA/PathB all agree with Jolien's full delay map to ~0.07 at fmax=256 (the ~26 +deficit at fmax=1024 was SimDetectorStrain's high-f TD delay-INTERPOLATION, not a bug -- +confirmed: it vanishes when fmax is lowered). NoLoop uses nearest-neighbor time sampling +(factored_likelihood.py ~L1691), so absolute peak comparisons at high SNR have a resolution +floor (~0.1-0.2); the time-interpolated NoLoop is in oshaughn/rift_O4d and +origin/rift_O4d_junior_calmarg_in_loop if a cleaner absolute comparison is wanted. + +## TWO PARALLEL THRUSTS (2026-07) +Path A + Path B are implemented, validated, and wired into the ILE. Next work is two tracks: +1. **Rotation PE value demo (near-term).** VERIFY-ANYWHERE quick-look DONE + (`~/RIFT_roboto_paper/analyses/slowrot_demo/`, `make demo` / `demo_local.py`): rotation-vs- + static lnL gain grows with Omega*T (null 0.004 -> 64s 0.19 at fixed SNR~30); figure + `outputs/gain_vs_duration.png`. REMAINING: the full DAG injection-recovery PE (posterior + bias + single-network sky map) on the cluster -> `paper/` figure cited from `sec:slowrot` + (structured Makefile targets inject/baseline/rotationA/rotationB/compare). +2. **Frequency-dependent (finite-size) response = Path D (3G regime where 'long' lives).** + RESPONSE FUNCTION IMPLEMENTED + VALIDATED: `slowrot_freqresponse.py` + + `test_slowrot_freqresponse.py`. F_k(f;RA,DEC,psi) from arXiv:2412.01693 single-arm sinc + transfer; f->0 == ComputeDetAMResponse to 6.7e-16; exact FSR null at c/2L=3747 Hz (40 km); + in-band SHAPE distortion 0.24%@1kHz/0.62%@2kHz (LIGO) vs 11.8%/42% (CE) -> negligible for + LIGO, first-order only for 3G. Key: the complex ratio ~1 for CE is dominated by a benign + common e^{-i2pi fL/c} delay (degenerate with tc), NOT sky shape. LAL has NO closed-form FD + response. ROUTE DECISION (notes/sec_freqdep.tex): CE -> route (b) sky-harmonic expansion + (keeps sky extrinsic; fold the common delay into geocenter time; few angular orders); + precessing+HM -> route (a) pinned-sky TD fold. + LIKELIHOOD INTEGRATION DONE + VALIDATED (route b): `factored_likelihood_freqresponse.py` + (+ `test_slowrot_freqresponse_likelihood.py`). Fold the common e^{-i2pi fL/c} delay into the + arrival time; power-series the residual transfer -> sky-independent W_p(f) folded into the FD + modes (reuse ComputeModeIP*) x analytic b_p(sky) -> SKY EXTRINSIC. MAINTAINED-STYLE NoLoop + entry point DiscreteFactoredLogLikelihoodFreqResponseNoLoop (NOT SingleDetectorLogLikelihood). + Validated (CE 40km, 16s): V1 finite-size(L->0) == maintained NoLoop baseline to 3.3e-9; V2 on + a finite-size injection the long-wavelength NoLoop deficit 2.71 -> finite-size 0.558 (converged + Qmax=2), residual ~= peak-resolution floor; V3 Cauchy-Schwarz respected. Scalar companion + agrees with the NoLoop to 0.156 (interp-vs-nearest floor). + REMAINING: ILE wiring (a --freqresponse flag mirroring --rotation-slow); full precessing+HM + (route a pinned-sky); a value demo (the finite-size effect only bites for CE/3G). + AUDIT NOTE (2026-07): all slowrot likelihoods route through the maintained NoLoop; there are + NO SingleDetectorLogLikelihood calls; the ILE wires only NoLoop paths. The scalar + FactoredLogLikelihood*/SingleDetectorLogLikelihood are used only as secondary references and + carry a peak-resolution floor at high SNR -- prefer NoLoop for any absolute comparison. + +## What this is +Generalizes RIFT's marginalized likelihood to account for the **time dependence of the +ground-based detector response over the signal** (Earth rotation), while reusing the +precompute-and-marginalize architecture. Two effects, both implemented (Path A + Path B): +- **Path A** — antenna amplitude drift `F(t)`: exact 5 sidereal harmonics, sky extrinsic. +- **Path B** — propagation-delay drift `tau(t)`: time-domain derivative expansion; adds + delay-derivative order `p` on top of Path A. Matters only for LONG signals. +- **Deferred (Path D)** — frequency-dependent (finite-size) response. Out of scope. + +## Files (all under RIFT/likelihood/ unless noted) +- `slowrot_response.py` — closed-form antenna harmonics `A_n` and delay harmonics `B_n` + (scalar + `*_vector`), derived from `vectorized_lal_tools` conventions. Validated vs LAL + to machine precision (`test_slowrot_response.py`). +- `factored_likelihood_with_rotation.py` — the core: + - `PrecomputeLikelihoodTermsWithRotation(...)` — FD-native precompute: builds each mode + once, applies derivative `(FT_SIGN*2pi i f)^p` (FT_SIGN=-1) and sidereal modulation + `exp(i n Omega t)` (sub-bin freq shift via LAL-FFT phase). Returns bank keyed by + elementary template `a=(p,n)`: `Q^a(t)`, `U^{(a,a')}`, `V^{(a,a')}`. + - `rotation_coefficients` / `rotation_coefficients_vector` — the analytic scalars + `C_{(p,ntilde)} = (1/p!) sum_{n+m=ntilde} A_tilde_n [(-D)^{*p}]_m` (Path A: `{(0,n): A_tilde_n}`). + - `FactoredLogLikelihoodWithRotation(...)` — scalar lnL (per-sample); term1 = `Re[sum_lm + conj(Ylm) sum_a conj(C_a) Q^a(t_det)]`, term2 with `U^{(a,a')}` (coef `conj(C_a)C_a'`) + and `V^{(a,a')}` (coef `C_{(p,-nu)} C_a'`). + - `pack_rotation_arrays` + `DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation` + — the MAINTAINED vectorized (NoLoop) path, mirroring + `factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopOrig`. +- `factored_likelihood.py` — one small FIX: non-numba `lalylm` fallback -> `np.vectorize` + (it broke `ComputeYlmsArrayVector` for the whole vectorized path in non-numba/non-GPU envs). +- `bin/integrate_likelihood_extrinsic_batchmode` — wiring: `--rotation-slow`, + `--rotation-n-harmonics`, `--rotation-p-max`. Guarded (requires `--vectorized`, CPU, no + dist/phase marg). Builds+packs the rotation bank after the baseline pack; the CPU + vectorized `likelihood_function` calls the rotation NoLoop. + +## How to run (RIFT venv; worktree first on PYTHONPATH) + source ~/RIFT_develUWM/bin/activate + export PYTHONPATH=~/RIFT_slowrot/MonteCarloMarginalizeCode/Code + python RIFT/likelihood/test_slowrot_response.py # A_n,B_n vs LAL + python RIFT/likelihood/test_slowrot_fd_ops.py # FD derivative/modulation + python RIFT/likelihood/test_slowrot_likelihood_v1.py # scalar Path A vs baseline + brute force + python RIFT/likelihood/test_slowrot_noloop.py # vectorized Path A vs baseline NoLoop + python RIFT/likelihood/test_slowrot_noloop_bruteforce.py # vectorized Path A vs brute force + python RIFT/likelihood/test_slowrot_pathB.py # Path B reduction + bound + python RIFT/likelihood/test_slowrot_headtohead.py # matched-sample rotation vs baseline (cubic) + python RIFT/likelihood/test_slowrot_freqresponse.py # [Path D] finite-size response vs LAL + python RIFT/likelihood/test_slowrot_freqresponse_likelihood.py # [Path D] likelihood: V1/V3 + V4 positive control + +VALUE DEMOS (verify-anywhere, no condor/GPU) -- consolidated in the RIFT tree: + cd demo/rift/slowrot && make demo # rotation (Path A/B) + finite-size (Path D) + # demo_rotation.py : gain lnL_rot-lnL_static grows with Omega*T (null control at short T) + # demo_finite_size.py : gain lnL_finite-lnL_LWL grows with arm length fL/c (null at LIGO 4km, + # +39.6 nats at CE 40km); writes outputs/*.{txt,png}. See its README.md. + +End-to-end ILE head-to-head (ILE-GPU-Paper demo data), baseline vs rotation vs finite-size: + D=~/RIFT_develUWM/src/research-projects-RIT/.travis/ILE-GPU-Paper/demos + integrate_likelihood_extrinsic_batchmode --vectorized [std opts] ... # baseline + integrate_likelihood_extrinsic_batchmode --vectorized --rotation-slow ... # Path A + integrate_likelihood_extrinsic_batchmode --vectorized --rotation-slow --rotation-p-max 1 ... # Path B + integrate_likelihood_extrinsic_batchmode --vectorized --freqresponse \ + --freqresponse-arm-length 40000 --freqresponse-qmax 6 ... # Path D (finite-size) + # --interpolate-time selects cubic sub-bin time interpolation for all of the above (default nearest). + +## Validation status (all PASSING) +- Response harmonics vs LAL: ~1e-16. FD ops vs LAL round trips: ~1e-13. +- Path A scalar: V1a (Omega=0 vs baseline) 2.7e-12; V1b (real vs brute force) 2.6e-9. +- Path A vectorized: vs baseline NoLoop 3.6e-12; vs brute force 3.2e-10; V0 (precompute + recovery on real data) exact. +- Path B: scalar reduce-to-baseline 9e-13; respects 0.5; vectorized reduce 6.4e-12. +- Path D (finite-size, --freqresponse): response Sum_p b_p W_p == antenna_response_fd to 6e-11 + on both +/-f; likelihood L->0 reduces to baseline NoLoop 3e-9; Cauchy-Schwarz respected; + V4 positive control asserts finite-size beats LWL by +38.9 nats (15+13 Msun, fmax=2000, CE 40km). +- Cubic time-interp (from calmarg_in_loop, --interpolate-time): both slow-response NoLoops now + support time_interp='nearest'|'cubic'. Cubic exposed+fixed a sub-bin GPS-cancellation bug in + the time reference; head-to-head regression floor 1.6e-3 -> 4.5e-13, test_slowrot_noloop 3.6e-12. +- ILE runs clean (no inf/nan) for baseline / Path A / Path B / Path D; rotation ~ baseline for the + short demo signal (rotation effect ~0.008 in lnL, negligible as physics requires). + +## CRITICAL lesson (a bug that hid for a while) +The U/V template modulation MUST be referenced to the template's INTRINSIC time (epoch ~0), +not the absolute event time t_ev~1e9 -- else a ~1e4 rad spurious phase randomizes U/V and +inflates lnL ABOVE 0.5. This bug once made a short-BBH rotation shift look like ~660 in +lnL; corrected it is ~0.008. My independent "brute force" reference shared the same +convention, so `vec == brute-force` PASSED while both were wrong. **Always cross-check +against the Cauchy-Schwarz bound 0.5, not only against a reference that can share +conventions.** + +## PATH B STATUS (findings 2026-07-04, the systematic pass in progress) +- Matched-seed head-to-head DONE (test_slowrot_headtohead.py): rot(f_sid=0)==baseline 9e-13; + evidence shift ln Z_rot - ln Z_base = -1.1e-3 (MC-noise-free) for the short signal. +- Path B DELAY PHYSICS still only partially validated. In an INFLATED-Omega regime (f_sidereal + x3000, single-det 30+25 BBH) where the delay drift matters, scalar Path B gives + lnL(p_max=0..3) = [1794.39, 1781.31, 1781.63, 928.30]: p=0->1->2 captures a real ~13-in-lnL + delay effect and appears to converge (~1781.6) and respects 0.5=1938 -- BUT p_max=3 + BLOWS UP (increment 853). +- ROOT CAUSE of the p>=3 blow-up (likely): the FD derivative weight (2 pi i f)^p amplifies high + frequencies; in the model norm the integrand ~ (2 pi f)^{2p} |h(f)|^2 / S grows like f^{11/3} + for a chirp (|h|^2 ~ f^{-7/3}), so high-order terms are dominated by the f_max edge, not the + physical low-frequency delay drift. FIX for the systematic pass: BAND-LIMIT the delay- + derivative (p>=1) terms to low frequency (the delay drift physically matters in the long early + inspiral, i.e. low f, not at merger). Options: apply a low-pass / taper before the (2 pi i f)^p + weight, or use a reduced f_max for p>=1 templates, or work with dimensionless (f/f_ref) weights. +- Consequence: current Path B is trustworthy only through p_max<=2 and only after the band-limit + fix is validated. Do NOT ship Path B (p>=1) until this is fixed AND checked against an + independent time-varying-delay brute force (below). + +## OPEN / NEXT (the one remaining systematic pass — do it all together) +1. **Path B rigorous validation** (TWO parts, do together): + (a) FIX the p>=3 high-frequency derivative blow-up by band-limiting the delay-derivative + terms to low frequency (see PATH B STATUS above); re-check convergence is monotone. + (b) Validate vs an INDEPENDENT ground truth that uses LAL's OWN full delay-time map -- + lalsim.SimDetectorStrainREAL8TimeSeries (Jolien's code). KEY FACTS FOUND 2026-07-04: + - RIFT's data ALREADY uses it: non_herm_hoff -> hoft -> SimInspiralTD + + SimDetectorStrainREAL8TimeSeries (lalsimutils.py ~line 3020). So a LONG injection + already carries the real Earth-rotation response; baseline (static) will under-recover + it and Path A/B should recover it. + - GROUND TRUTH for the rotation likelihood must apply SimDetectorStrain to the SAME + RIFT modes the likelihood uses: hk = lsu.hoft_from_hlm(hlmsT, P_extr) (radec path, + applies SimDetectorStrain). Then convert to 2-sided FD exactly like non_herm_hoff + (copy REAL8 -> COMPLEX16, COMPLEX16TimeFreqFFT) and lnL = ComplexIP(d,hk).real - + 0.5*ComplexIP(hk,hk).real. + - DO NOT use non_herm_hoff(P_extr) as the template ground truth: it regenerates via + SimInspiralTD, whose modes differ from RIFT's hlmoff modes -- with IMRPhenomD the + mode-based baseline recovers only ~SNR 35 of the data's ~SNR 82 at the injection (a + waveform-convention mismatch, NOT rotation). Using hoft_from_hlm keeps the modes + identical so ONLY the rotation differs. + - Two gotchas to fix: (i) hoft_from_hlm output length = mode length (e.g. 16451) != + padded data length (16384) -> ResizeREAL8TimeSeries + align epoch (replicate + non_herm_hoff padding); (ii) the mode-based likelihood at a FIXED tref is OFF-PEAK + (epoch bookkeeping) -> compare at the time-MAXIMIZED / marginalized value, not at + extr.tref. + - Use a LONG signal (real sidereal rate) so delay drift matters, OR the inflated-Omega + path (already shows p=1,2 capture the delay, p>=3 needs the band-limit fix). + A convergence test alone is NOT enough: Path B could converge to a WRONG value if the C + coefficients are off (the delay analogue of the reference-time bug). Cross-check 0.5. +2. **Matched-seed quantitative ILE head-to-head** (fixed RNG seed) -> exact, regression-grade + baseline-vs-rotation comparison (currently only MC-noise-level agreement). +3. Cauchy-Schwarz bound checks everywhere; sweep for other shared-convention bugs. +4. Possible follow-ups: GPU (xpy) path for the rotation NoLoop; distance/phase/cal marg + support; Path B `t_star` reference-epoch choice; `rotation_coefficients` closed-form + `A_n`/`B_n` appendix (notes sec app_response has the A_n,B_n). + +## Also note +- A separate background session was fixing the py2-ism in + `factored_likelihood.PackLikelihoodDataStructuresAsArrays` (`rholm_intpArray = range(nKeys)` + fails in py3 when the interpolant arg is truthy). Our tests pass `None` for that arg to + dodge it. +- Subagents were unreliable for this work (deferred, reverted fixes, over-claimed). The + validated results above were done/checked directly. diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py index 44e32c115..2232e0c9a 100644 --- a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood.py @@ -1396,7 +1396,7 @@ def PackLikelihoodDataStructuresAsArrays(pairKeys, rholms_intpDictionaryForDetec rholmArray[indx1][:] = rholmsDictionaryForDetector[pair1].data.data # Copy the array of time values. ### Step 3: Create rholm_intp array-ized structure - rholm_intpArray = range(nKeys) # create a flexible python array of the desired size, to hold function pointers + rholm_intpArray = [None] * nKeys # create a flexible python array of the desired size, to hold function pointers if rholms_intpDictionaryForDetector: for pair1 in pairKeys: indx1 = lookupKeysToNumber[pair1] @@ -2402,9 +2402,10 @@ def lalT(det, RA, DEC,tref): # note tref is a SCALAR numba_on = False if log_loud: print(" Numba off ") - # Very inefficient - def lalylm(th,ph,s,l,m): - return lal.SpinWeightedSphericalHarmonic(th,ph,s,l,m) + # Very inefficient. np.vectorize replicates the numba @vectorize elementwise + # behavior (the plain scalar wrapper below breaks ComputeYlmsArrayVector, which + # passes numpy ARRAYS for th,ph,s,l,m -> LAL raises 'argument 1 of type REAL8'). + lalylm = np.vectorize(lal.SpinWeightedSphericalHarmonic, otypes=[complex]) def lalF(det, RA, DEC,psi,tref): if isinstance(RA, float): return ComplexAntennaFactor(det, RA, DEC, psi,tref) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood_freqresponse.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood_freqresponse.py new file mode 100644 index 000000000..2391d75ee --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood_freqresponse.py @@ -0,0 +1,408 @@ +""" +factored_likelihood_freqresponse : finite-size (frequency-dependent detector response) +generalization of the RIFT precompute + scalar likelihood. "Thrust 2 / route (b)" +(sky-harmonic expansion) of the slow-rotation program -- KEEPS the sky extrinsic. + +Physics (see slowrot_freqresponse.py for the response derivation) +---------------------------------------------------------------- +Beyond the long-wavelength limit the detector strain in the frequency domain is a +per-frequency antenna pattern acting on the two polarizations, + + h_k(f) = F_+(f;sky) h_+(f) + F_x(f;sky) h_x(f) + = (1/2)[ F(f) Sigma(f) + Fbar(f) Sigma*(f) ] , + F = F_+ + i F_x, Fbar = F_+ - i F_x = conj(F(-f)), Sigma = sum_lm Y_lm h_lm. + +The finite-size response factors into a SKY-INDEPENDENT frequency basis W_p(f) times +an analytic SKY/pol scalar b_p (slowrot_freqresponse.finite_size_response_weights / +finite_size_beta): + + F(f;sky) = sum_p b_p(sky) W_p(f) , Fbar(f;sky) = sum_p conj(b_p) W_p(f) , + + p=0 : W_0 = 1, b_0 = F0 (exact lal ComputeDetAMResponse) + p=1+q : W_{1+q} = e^{-i2pi f T} c_q(f) - [q==0], b_{1+q} = beta_q (arm expansion) + +with T = L/c the common (direction-independent) light-crossing delay folded into W_p as +a linear FD phase, and c_q(f) the sky-independent power-series basis of the single-arm +transfer. Every W_p is Hermitian, so -- UNLIKE the sidereal-rotation case -- the V cross +term needs NO harmonic reflection. + +This module is the DIRECT analogue of factored_likelihood_with_rotation: it folds each +W_p(f) into the FD modes ONCE (h_lm(f) -> W_p(f) h_lm(f)) and reuses the EXISTING +ComputeModeIPTimeSeries / ComputeModeCrossTermIP to build + + Q^p_lm(t) = < W_p h_lm | d >, U^{p,p'} = < W_p h_lm | W_p' h_l'm' >, + V^{p,p'} = < conj(W_p h_lm) | W_p' h_l'm' > . + +MAINTAINED PATH: the production entry point is the vectorized (NoLoop) +DiscreteFactoredLogLikelihoodFreqResponseNoLoop -- the direct finite-size analogue of +factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop (the maintained +likelihood used by integrate_likelihood_extrinsic_batchmode). It is VALIDATED against that +NoLoop baseline: at L->0 it reduces to it (to ~1e-9), NOT to the older per-detector scalar +SingleDetectorLogLikelihood. FactoredLogLikelihoodFreqResponse (scalar) is a per-sample +validation companion only (mirrors factored_likelihood_with_rotation.FactoredLogLikelihoodWithRotation). + +The term1/term2 use the same factored form as the NoLoop path (all RIFT factored likelihoods +share this algebra; W_0=1 and b_0 = ComputeDetAMResponse at L->0): + term1 = Re[ (distRef/dist) sum_lm conj(Ylm) sum_p conj(b_p) Q^p_lm(t_det) ] + term2 = -(1/4)(distRef/dist)^2 Re[ sum_{p,p'} ( conj(b_p) b_p' U^{p,p'} conj(Y1)Y2 + + b_p b_p' V^{p,p'} Y1 Y2 ) ] + +Heavy RIFT imports are lazy so the light FD helpers stay importable with numpy alone. +""" +from __future__ import print_function, division + +import numpy as np + +from . import slowrot_freqresponse as sfr + + +# --------------------------------------------------------------------------- +# Low-level FD primitive: fold a per-bin weight into a COMPLEX16FrequencySeries. +# --------------------------------------------------------------------------- +def evaluate_fvals_from_length(npts, deltaF): + """Signed frequency array for a RIFT two-sided COMPLEX16FrequencySeries. + + Identical packing to factored_likelihood_with_rotation.evaluate_fvals_from_length / + RIFT.lalsimutils.evaluate_fvals: f[k] = deltaF*(npts/2 - k), descending from +fNyq. + """ + k = np.arange(npts) + return deltaF * (npts / 2.0 - k) + + +def _copy_freqseries(hf): + import lal + out = lal.CreateCOMPLEX16FrequencySeries( + hf.name, hf.epoch, hf.f0, hf.deltaF, hf.sampleUnits, hf.data.length) + out.data.data[:] = hf.data.data[:] + return out + + +def fd_apply_weight(hf, weight): + """COMPLEX16FrequencySeries -> new series with spectrum multiplied by weight[k].""" + out = _copy_freqseries(hf) + out.data.data[:] = hf.data.data * weight + return out + + +# --------------------------------------------------------------------------- +# The precompute (sky-independent; per detector uses only its arm length L). +# --------------------------------------------------------------------------- +def PrecomputeLikelihoodTermsFreqResponse( + event_time_geo, t_window, P, data_dict, psd_dict, Lmax, fMax, + Qmax=4, L_arm=None, + analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0., + verbose=True, quiet=False, skip_interpolation=False, **hlm_kwargs): + """Finite-size analogue of PrecomputeLikelihoodTerms(WithRotation). + + Builds each FD mode once and forms the generalized overlaps for every response + basis element p = 0..Qmax+1: + + rholms_intp_fr[det][p][(l,m)] : interpolant of Q^p_lm(t) = + crossTerms_fr[det][(p,p')] : { ((l,m),(l',m')) : } + crossTermsV_fr[det][(p,p')] : { ((l,m),(l',m')) : } + + Parameters mirror PrecomputeLikelihoodTerms; finite-size specific: + Qmax : highest power of the arm projection retained (basis size Qmax+2). + L_arm : arm-length override [m] (float -> all detectors, or dict det->L). + None -> each detector's native LAL arm length. + + The response coefficients b_p(det,RA,DEC,psi,tref) that contract these overlaps are + NOT computed here -- see FactoredLogLikelihoodFreqResponse. Returns the intrinsic- + only, sky-independent overlap bank plus meta (carries per-detector L, T, Qmax). + """ + import lal + from . import factored_likelihood as FL + from .. import lalsimutils as lsu + + assert data_dict.keys() == psd_dict.keys() + detectors = list(data_dict.keys()) + + def _L_of(det): + if isinstance(L_arm, dict): + return L_arm.get(det, None) + return L_arm + + P.dist = FL.distMpcRef * 1e6 * lsu.lsu_PC + P.deltaF = data_dict[detectors[0]].deltaF + + # --- build base FD modes ONCE --- + hlms, hlms_conj = FL.internal_hlm_generator(P, Lmax, verbose=verbose, quiet=quiet, + **hlm_kwargs) + modes = list(hlms.keys()) + npts = hlms[modes[0]].data.length + deltaF = hlms[modes[0]].deltaF + fvals = evaluate_fvals_from_length(npts, deltaF) + p_list = list(range(Qmax + 2)) + + rholms_fr = {} + rholms_intp_fr = {} + crossTerms_fr = {} + crossTermsV_fr = {} + Ldict = {} + Tdict = {} + + for det in detectors: + psd = psd_dict[det] + data = data_dict[det] + + # per-detector geometry (arm length only; sky enters later via b_p) + _, _, _, L = sfr.detector_geometry(det, L_arm=_L_of(det)) + Ldict[det] = float(L) + Tdict[det] = float(L / sfr.C_SI) + # sky-independent response basis weights W_p(f) on the signed fvals + geom_dummy = dict(L=float(L), T=float(L / sfr.C_SI)) + W = sfr.finite_size_response_weights(fvals, geom_dummy, Qmax) # (Qmax+2, npts) + + # weighted mode families eta^p = W_p h_lm , etac^p = W_p conj(h_lm) + eta = {p: {lm: fd_apply_weight(hlms[lm], W[p]) for lm in modes} for p in p_list} + etac = {p: {lm: fd_apply_weight(hlms_conj[lm], W[p]) for lm in modes} for p in p_list} + + t_det = FL.ComputeArrivalTimeAtDetector(det, P.phi, P.theta, event_time_geo) + rho_epoch = data.epoch - hlms[modes[0]].epoch + t_shift = float(float(t_det) - float(t_window) - float(rho_epoch)) + N_shift = int(t_shift / P.deltaT + 0.5) + N_window = int(2 * t_window / P.deltaT) + t = np.arange(N_window) * P.deltaT + float(rho_epoch + N_shift * P.deltaT) + + rholms_fr[det] = {} + rholms_intp_fr[det] = {} + for p in p_list: + rho = FL.ComputeModeIPTimeSeries( + eta[p], data, psd, P.fmin, fMax, 1. / 2. / P.deltaT, + N_shift, N_window, analyticPSD_Q, inv_spec_trunc_Q, T_spec) + rholms_fr[det][p] = rho + if not skip_interpolation: + rholms_intp_fr[det][p] = FL.InterpolateRholms(rho, t, verbose=verbose) + else: + rholms_intp_fr[det][p] = None + + crossTerms_fr[det] = {} + crossTermsV_fr[det] = {} + for p in p_list: + for pp in p_list: + crossTerms_fr[det][(p, pp)] = FL.ComputeModeCrossTermIP( + eta[p], eta[pp], psd, P.fmin, fMax, 1. / 2. / P.deltaT, P.deltaF, + analyticPSD_Q, inv_spec_trunc_Q, T_spec, verbose=False, + same_waveform_Q=False) + crossTermsV_fr[det][(p, pp)] = FL.ComputeModeCrossTermIP( + etac[p], eta[pp], psd, P.fmin, fMax, 1. / 2. / P.deltaT, P.deltaF, + analyticPSD_Q, inv_spec_trunc_Q, T_spec, prefix="V", verbose=False, + same_waveform_Q=False) + + meta = dict(Qmax=Qmax, p_list=p_list, modes=modes, + event_time_geo=float(event_time_geo), + L=Ldict, T=Tdict, L_arm=L_arm) + return rholms_intp_fr, crossTerms_fr, crossTermsV_fr, rholms_fr, meta + + +# --------------------------------------------------------------------------- +# Response coefficients b_p(det, RA, DEC, psi, tref) : b_0 = F0, b_{1+q} = beta_q. +# --------------------------------------------------------------------------- +def response_coefficients(det, RA, DEC, psi, tref, Qmax, L_arm=None): + """{p: b_p} response coefficients for the finite-size basis at (RA,DEC,psi,tref). + + b_0 = F0 = lal.ComputeDetAMResponse (exact LWL baseline); + b_{1+q} = beta_q(sky) = (1/2)[zx^2 a_x^q - zy^2 a_y^q]. + Uses gmst = GreenwichMeanSiderealTime(tref) exactly like ComplexAntennaFactor. + """ + import lal + gmst = float(lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(float(tref)))) + geom = sfr.finite_size_geometry(det, RA, DEC, psi, gmst=gmst, L_arm=L_arm) + beta = sfr.finite_size_beta(geom, Qmax) + b = {0: geom['F0']} + for q in range(Qmax + 1): + b[1 + q] = complex(beta[q]) + return b + + +def FactoredLogLikelihoodFreqResponse(extr_params, rholms_intp_fr, crossTerms_fr, + crossTermsV_fr, meta, Lmax): + """Finite-size analogue of factored_likelihood.FactoredLogLikelihood. + + Contracts the response-basis precompute bank with b_p(det,RA,DEC,psi,tref) and the + Ylm. Reduces EXACTLY to the baseline FactoredLogLikelihood as L -> 0. + """ + from . import factored_likelihood as FL + from .. import lalsimutils as lsu + + Qmax = meta['Qmax'] + p_list = list(meta['p_list']) + L_arm = meta.get('L_arm', None) + + RA = extr_params.phi + DEC = extr_params.theta + tref = extr_params.tref + phiref = extr_params.phiref + incl = extr_params.incl + psi = extr_params.psi + dist = extr_params.dist + + detectors = list(rholms_intp_fr.keys()) + p0 = p_list[0] + modes = list(rholms_intp_fr[detectors[0]][p0].keys()) + Ylms = FL.ComputeYlms(Lmax, incl, -phiref, selected_modes=modes) + + distMpc = dist / (lsu.lsu_PC * 1e6) + invDistMpc = FL.distMpcRef / distMpc + + def _L_of(det): + if isinstance(L_arm, dict): + return L_arm.get(det, None) + return L_arm + + lnL = 0. + for det in detectors: + b = response_coefficients(det, RA, DEC, psi, tref, Qmax, L_arm=_L_of(det)) + t_det = FL.ComputeArrivalTimeAtDetector(det, RA, DEC, tref) + CT = crossTerms_fr[det] + CTV = crossTermsV_fr[det] + + Q = {p: {m: rholms_intp_fr[det][p][m](float(t_det)) + for m in modes} for p in p_list} + + term1 = 0. + for m in modes: + s = 0. + for p in p_list: + s += np.conj(b[p]) * Q[p][m] + term1 += np.conj(Ylms[m]) * s + term1 = np.real(term1) * invDistMpc + + term2 = 0. + for m1 in modes: + for m2 in modes: + u = 0. + v = 0. + for p in p_list: + for pp in p_list: + u += np.conj(b[p]) * b[pp] * CT[(p, pp)][(m1, m2)] + v += b[p] * b[pp] * CTV[(p, pp)][(m1, m2)] + term2 += u * np.conj(Ylms[m1]) * Ylms[m2] + v * Ylms[m1] * Ylms[m2] + term2 = -np.real(term2) / 4. / (distMpc / FL.distMpcRef) ** 2 + + lnL += term1 + term2 + + return lnL + + +# --------------------------------------------------------------------------- +# Dense-array packing + vectorized time-marginalized lnL (for peak / validation). +# --------------------------------------------------------------------------- +def pack_freqresponse_arrays(meta, rholms_fr, crossTerms_fr, crossTermsV_fr): + """Pack the basis precompute bank into dense arrays keyed per detector by p. + + Returns (lookupNKDict, rho_by_p, U_by_pp, V_by_pp, epochDict). + """ + p_list = list(meta['p_list']) + lookupNKDict = {}; rho_by_p = {}; U_by_pp = {}; V_by_pp = {}; epochDict = {} + for det in rholms_fr: + p0 = p_list[0] + modes = list(rholms_fr[det][p0].keys()) + n_lms = len(modes) + idx = {m: i for i, m in enumerate(modes)} + lookupNKDict[det] = np.array([[m[0], m[1]] for m in modes], dtype=int) + npts = rholms_fr[det][p0][modes[0]].data.length + epochDict[det] = float(rholms_fr[det][p0][modes[0]].epoch) + rho_by_p[det] = {} + for p in p_list: + arr = np.zeros((n_lms, npts), dtype=np.complex128) + for m in modes: + arr[idx[m]] = rholms_fr[det][p][m].data.data + rho_by_p[det][p] = arr + U_by_pp[det] = {}; V_by_pp[det] = {} + for p in p_list: + for pp in p_list: + Um = np.zeros((n_lms, n_lms), dtype=np.complex128) + Vm = np.zeros((n_lms, n_lms), dtype=np.complex128) + cU = crossTerms_fr[det][(p, pp)] + cV = crossTermsV_fr[det][(p, pp)] + for m1 in modes: + for m2 in modes: + Um[idx[m1], idx[m2]] = cU[(m1, m2)] + Vm[idx[m1], idx[m2]] = cV[(m1, m2)] + U_by_pp[det][(p, pp)] = Um + V_by_pp[det][(p, pp)] = Vm + return lookupNKDict, rho_by_p, U_by_pp, V_by_pp, epochDict + + +def DiscreteFactoredLogLikelihoodFreqResponseNoLoop( + tvals, P_vec, meta, lookupNKDict, rho_by_p, U_by_pp, V_by_pp, epochDict, + Lmax=2, array_output=False, time_interp='nearest'): + """Vectorized finite-size lnL over a time window (single extrinsic point per call + is supported; P_vec fields RA,DEC,incl,phiref,psi,dist may be length-1 arrays). + + array_output=True returns lnL_t of shape (npts_ex, npts); else time-marginalized. + """ + import lal + from . import factored_likelihood as FL + + p_list = list(meta['p_list']) + Qmax = meta['Qmax'] + L_arm = meta.get('L_arm', None) + npts = len(tvals); npts_ex = len(np.atleast_1d(P_vec.phi)) + RA = np.atleast_1d(P_vec.phi); DEC = np.atleast_1d(P_vec.theta) + incl = np.atleast_1d(P_vec.incl); phiref = np.atleast_1d(P_vec.phiref) + psi = np.atleast_1d(P_vec.psi) + distMpc = np.atleast_1d(P_vec.dist) / (lal.PC_SI * 1e6) + lnL_t = np.zeros((npts_ex, npts), dtype=np.float64) + + def _L_of(det): + if isinstance(L_arm, dict): + return L_arm.get(det, None) + return L_arm + + for det in rho_by_p: + n_lms = len(lookupNKDict[det]) + Ylms = FL.ComputeYlmsArrayVector(lookupNKDict[det], incl, -phiref).T # (npts_ex,n_lms) + # per-sample response coefficients b_p (npts_ex,) + bvec = {} + for i in range(npts_ex): + bi = response_coefficients(det, float(RA[i]), float(DEC[i]), float(psi[i]), + P_vec.tref, Qmax, L_arm=_L_of(det)) + for p in p_list: + bvec.setdefault(p, np.zeros(npts_ex, dtype=complex)) + bvec[p][i] = bi[p] + + t_ref = epochDict[det] + # Precision-preserving time reference (see rotation NoLoop): keep (tref - epoch) and the + # geometric delay separate so the sub-bin fraction is exact under cubic interpolation. + detector_location = np.asarray(FL.lalsim.DetectorPrefixToLALDetector(det).location) + gmst_tref = float(lal.GreenwichMeanSiderealTime(P_vec.tref)) + t_det = float(P_vec.tref - float(t_ref)) + FL.TimeDelayFromEarthCenter( + detector_location, RA, DEC, gmst_tref) + sample_first = (t_det + tvals[0]) / P_vec.deltaT + if time_interp == 'nearest': + ifirst = (np.round(sample_first) + 0.5).astype(int) + else: + ifirst = np.floor(sample_first).astype(int) + frac_first = (sample_first - np.floor(sample_first)).astype(np.float64) + ilast = ifirst + npts + + term1 = np.zeros((npts_ex, npts), dtype=np.complex128) + for p in p_list: + det_rho = rho_by_p[det][p] + if time_interp == 'nearest': + Qa = np.empty((npts_ex, npts, n_lms), dtype=np.complex128) + for i in range(npts_ex): + Qa[i] = det_rho[..., ifirst[i]:ilast[i]].T + else: + Qa = FL._cubic_Q_window_numpy(det_rho.T, ifirst, frac_first, npts) + term1 += np.conj(bvec[p])[:, None] * np.einsum('xi,xti->xt', np.conj(Ylms), Qa) + term1 = term1.real * (FL.distMpcRef / distMpc)[:, None] + + term2 = np.zeros(npts_ex, dtype=np.complex128) + for p in p_list: + for pp in p_list: + term2 += np.conj(bvec[p]) * bvec[pp] * np.einsum( + 'xi,xj,ij->x', np.conj(Ylms), Ylms, U_by_pp[det][(p, pp)]) + term2 += bvec[p] * bvec[pp] * np.einsum( + 'xi,xj,ij->x', Ylms, Ylms, V_by_pp[det][(p, pp)]) + term2 = (-0.25 * term2.real) * (FL.distMpcRef / distMpc) ** 2 + + lnL_t += term1 + term2[:, None] + + if array_output: + return lnL_t + lnLmax = np.max(lnL_t, axis=-1, keepdims=True) + L = FL.my_simps(np.exp(lnL_t - lnLmax), dx=P_vec.deltaT, axis=-1) + return lnLmax[..., 0] + np.log(L) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood_with_rotation.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood_with_rotation.py new file mode 100644 index 000000000..6f232bcfb --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/factored_likelihood_with_rotation.py @@ -0,0 +1,619 @@ +""" +factored_likelihood_with_rotation : slow-rotation generalization of the RIFT precompute. + +This is a *separate realization* of factored_likelihood.PrecomputeLikelihoodTerms that +accounts for the time dependence of the ground-based detector response over long signals +(Earth rotation). It does NOT regenerate a family of time-domain templates. Instead it +builds each frequency-domain mode ONCE and realizes the two response modulations as cheap +frequency-domain operations on that single template: + + * delay-derivative order p (Path B): time derivative h^(p)(t) <-> (s 2 pi i f)^p H(f) + -- an exact per-bin weight in the FD. (s = FT_SIGN, fixed to match RIFT's FFT + convention; see evaluate_fvals, which itself warns the omega-sign is reversed.) + + * sidereal harmonic n (Paths A,B): modulation exp(i n Omega t) <-> frequency + shift H(f - n f_sid). Since f_sid = Omega/2pi ~ 1.16e-5 Hz is sub-bin, the shift is + realized exactly as a time-domain linear phase (one FFT round trip per template), not + an integer bin roll. + +The response harmonics themselves (the analytic scalars that multiply these precomputed +overlaps in the extrinsic layer) come from RIFT.likelihood.slowrot_response: + + F_k(t) = sum_{n=-2}^{2} A_n e^{i n g}, tau_k(t) = sum_{n=-1}^{1} B_n e^{i n g}, + g = GMST(t) - RA, A_n, B_n depend only on (detector, dec, psi/dec). + +See the design notes (rift-slow-rotation: sec_formalism 'master expansion', sec_amplitude +Path A, sec_delay Path B, app_response) for the derivation and the meaning of the indices. + +Index conventions in the returned structures +-------------------------------------------- +An "elementary modulated template" is labelled a = (p, n): + chi_a(t) = exp(i n Omega t) * d^p/dt^p h_lm(t - tau_0). +The physical data-term time series carries a post-phase (derived in the notes): + Q^a_lm(t) = exp(i n Omega t) * < chi_a(.-t) | d > [applied here] +while the cross terms are arrival-time independent: + U^{a,a'} = < chi_a | chi_a' >, V^{a,a'} = < chi_a^* | chi_a' >. + +Path A (default) uses only p = 0 (amplitude drift; exact 5-harmonic). Path B adds p >= 1. + +Heavy RIFT imports (factored_likelihood, lalsimutils) are done lazily inside the precompute +so the light FD primitives below are importable with numpy alone (used by the unit tests). +""" +from __future__ import print_function, division + +import numpy as np + +# Sidereal angular rate [rad/s] and frequency [Hz] +OMEGA_EARTH = 7.292115e-5 +F_SIDEREAL = OMEGA_EARTH / (2.0 * np.pi) + +# Sign in the time-derivative weight (s 2 pi i f)^p, matched to RIFT's evaluate_fvals FFT +# convention. VALIDATED empirically against a LAL FFT round trip in +# test_slowrot_fd_ops.py (which will fail loudly if this is wrong). +FT_SIGN = -1.0 + + +# --------------------------------------------------------------------------- +# Low-level FD primitives (numpy only; operate on a complex spectrum + its fvals). +# --------------------------------------------------------------------------- +def evaluate_fvals_from_length(npts, deltaF): + """Signed frequency array for a RIFT two-sided COMPLEX16FrequencySeries. + + Replicates RIFT.lalsimutils.evaluate_fvals (kept local so these primitives need no + heavyweight import). Note RIFT's expression ``npts/2 - k if k<=npts/2 else -k+npts/2`` + is the same in both branches, i.e. simply f[k] = deltaF*(npts/2 - k), descending from + +fNyq at k=0. + """ + k = np.arange(npts) + return deltaF * (npts / 2.0 - k) + + +def time_derivative_weight(fvals, p): + """(FT_SIGN * 2 pi i f)^p : exact FD weight for the p-th time derivative.""" + if p == 0: + return np.ones_like(fvals, dtype=complex) + return (FT_SIGN * 2.0j * np.pi * fvals) ** p + + +def apply_time_derivative_array(spectrum, fvals, p): + """Return the spectrum of d^p/dt^p h given the spectrum of h (both two-sided).""" + return spectrum * time_derivative_weight(fvals, p) + + +def apply_sidereal_modulation_array(spectrum, coef, deltaF, f_sidereal=F_SIDEREAL): + """Return the spectrum of exp(i coef Omega t) h(t) given the spectrum of h(t). + + exp(i coef Omega t) multiplication in time <-> frequency shift by coef*f_sidereal. + This is sub-bin (f_sidereal << deltaF for realistic segments), so it is applied exactly + as a time-domain linear phase via a DFT round trip that respects RIFT's fvals packing. + + WARNING: this reference implementation forms an explicit (npts x npts) DFT matrix -- it + is O(npts^2) and intended ONLY for unit tests at small npts. The production path uses + the LAL-FFT round trip in _lal_freq_modulate() below (O(npts log npts)). + """ + if coef == 0: + return spectrum.copy() + npts = len(spectrum) + dt = 1.0 / (npts * deltaF) + t = np.arange(npts) * dt + fvals = evaluate_fvals_from_length(npts, deltaF) + # Round-trip-consistent DFT pair for the packing f[k]=deltaF*(npts/2-k): + # h(t_j) = (1/npts) sum_k H_k exp(+2 pi i f_k t_j), H_k = sum_j h(t_j) exp(-2 pi i f_k t_j) + # LAL's convention (encoded by evaluate_fvals) reverses the omega sign: a tone + # exp(+2 pi i f0 t) lands at fvals = -f0. The consistent DFT pair is therefore + # h(t_j) = (1/npts) sum_k H_k exp(-2 pi i f_k t_j), H_k = sum_j h(t_j) exp(+2 pi i f_k t_j). + phase_inv = np.exp(-2.0j * np.pi * np.outer(t, fvals)) # (npts, npts): H -> h(t) + h_td = phase_inv.dot(spectrum) / npts + h_td *= np.exp(1.0j * coef * 2.0 * np.pi * f_sidereal * t) # exp(i coef Omega t) + phase_fwd = np.exp(2.0j * np.pi * np.outer(fvals, t)) # h(t) -> H + return phase_fwd.dot(h_td) + + +# --------------------------------------------------------------------------- +# LAL-series wrappers (used inside the precompute; import lal lazily). +# --------------------------------------------------------------------------- +def _copy_freqseries(hf): + import lal + out = lal.CreateCOMPLEX16FrequencySeries( + hf.name, hf.epoch, hf.f0, hf.deltaF, hf.sampleUnits, hf.data.length) + out.data.data[:] = hf.data.data[:] + return out + + +def fd_apply_time_derivative(hf, p): + """COMPLEX16FrequencySeries -> new series for the p-th time derivative.""" + if p == 0: + return _copy_freqseries(hf) + out = _copy_freqseries(hf) + fvals = evaluate_fvals_from_length(hf.data.length, hf.deltaF) + out.data.data[:] = apply_time_derivative_array(hf.data.data, fvals, p) + return out + + +def _lal_freq_modulate(hf, coef, f_sidereal=F_SIDEREAL, t_ref=0.0): + """Modulate by exp(i coef Omega (t_abs - t_ref)) via a LAL-FFT round trip (O(N log N)). + + Reverse-FFT the spectrum to the time domain, multiply by the sidereal linear phase + referenced to ABSOLUTE sample time t_abs = float(hf.epoch) + j*deltaT (minus t_ref), + forward-FFT back. Uses the same COMPLEX16 transforms RIFT uses for its overlaps. + + The reference t_ref is physical, not cosmetic: the true antenna phase is + exp(i n (GMST(t')-RA)) = exp(i n (GMST(t_ev)-RA)) * exp(i n Omega (t'-t_ev)), so the + precompute must carry exactly exp(i n Omega (t' - t_ev)) at absolute data time t', + with the constant GMST(t_ev) piece carried analytically by A_n (slowrot_response). + Hence callers pass t_ref = event_time_geo. + """ + import lal + if coef == 0: + return _copy_freqseries(hf) + npts = hf.data.length + deltaT = 1.0 / (npts * hf.deltaF) + ts = lal.CreateCOMPLEX16TimeSeries("m", hf.epoch, 0., deltaT, + lal.DimensionlessUnit, npts) + revplan = lal.CreateReverseCOMPLEX16FFTPlan(npts, 0) + fwdplan = lal.CreateForwardCOMPLEX16FFTPlan(npts, 0) + lal.COMPLEX16FreqTimeFFT(ts, hf, revplan) + t_abs = float(hf.epoch) + np.arange(npts) * deltaT + ts.data.data[:] = ts.data.data * np.exp( + 1.0j * coef * 2.0 * np.pi * f_sidereal * (t_abs - t_ref)) + out = _copy_freqseries(hf) + lal.COMPLEX16TimeFreqFFT(out, ts, fwdplan) + return out + + +def fd_apply_sidereal_modulation(hf, n, f_sidereal=F_SIDEREAL, t_ref=0.0): + """COMPLEX16FrequencySeries -> new series for exp(i n Omega (t_abs - t_ref)) h(t).""" + return _lal_freq_modulate(hf, n, f_sidereal, t_ref) + + +def build_elementary_template(hf, p, n, f_sidereal=F_SIDEREAL): + """chi_{(p,n)} spectrum from the base mode spectrum hf: derivative then modulation.""" + return fd_apply_sidereal_modulation(fd_apply_time_derivative(hf, p), n, f_sidereal) + + +def _elementary_index_set(harmonics, p_max): + """List of a=(p,n). Superset grid; unused (p,n) get zero response coefficients.""" + return [(p, n) for p in range(p_max + 1) for n in harmonics] + + +# --------------------------------------------------------------------------- +# The precompute. +# --------------------------------------------------------------------------- +def PrecomputeLikelihoodTermsWithRotation( + event_time_geo, t_window, P, data_dict, psd_dict, Lmax, fMax, + harmonics=(-2, -1, 0, 1, 2), p_max=0, f_sidereal=F_SIDEREAL, + analyticPSD_Q=False, inv_spec_trunc_Q=False, T_spec=0., + verbose=True, quiet=False, internal_fast_precompute=True, + skip_interpolation=False, **hlm_kwargs): + """Slow-rotation analogue of factored_likelihood.PrecomputeLikelihoodTerms. + + Builds each FD mode once (via factored_likelihood.internal_hlm_generator) and forms the + generalized overlaps for every elementary modulated template a=(p,n): + + rholms_intp_rot[det][a][(l,m)] : interpolant of Q^a_lm(t) = e^{inOmega t} + crossTerms_rot[det][(a,a')] : { ((l,m),(l',m')) : } + crossTermsV_rot[det][(a,a')] : { ((l,m),(l',m')) : } + + Parameters mirror PrecomputeLikelihoodTerms; rotation-specific: + harmonics : sidereal harmonic indices n to carry (antenna needs |n|<=2). + p_max : max delay-derivative order (0 = Path A amplitude-only; >=1 = Path B). + f_sidereal: sidereal frequency [Hz]. + + Response coefficients A_n(det,dec,psi), B_n(det,dec) that contract these overlaps are + NOT computed here -- see slowrot_response and the (forthcoming) rotation-aware + FactoredLogLikelihood assembly. This routine produces only the intrinsic-only, + sky-independent overlap bank. + + NOTE: full-stack numerical validation (V1 of the notes' validation matrix: agreement + with a brute-force time-varying-response likelihood) requires the data/PSD/waveform + environment and is done separately; the FD primitives used here are unit-tested in + test_slowrot_fd_ops.py. + """ + # Lazy heavy imports (need the full RIFT stack / lal). + import lal + from . import factored_likelihood as FL + from .. import lalsimutils as lsu + + assert data_dict.keys() == psd_dict.keys() + detectors = list(data_dict.keys()) + t_ev = float(event_time_geo) + # The exp(i n Omega t) modulation for the data term Q is applied to the DATA (shift by + # -n f_sidereal, referenced to t_ev), which is mode-independent: one shift per (det,n), + # and -- since the modulation lives on the fixed absolute data-time axis -- needs NO + # arrival-time-dependent post-phase. U,V use modulated templates (same t_ev reference). + + # Reference distance handling identical to the base precompute. + P.dist = FL.distMpcRef * 1e6 * lsu.lsu_PC + P.deltaF = data_dict[detectors[0]].deltaF + + # --- build base FD modes ONCE --- + hlms, hlms_conj = FL.internal_hlm_generator(P, Lmax, verbose=verbose, quiet=quiet, + **hlm_kwargs) + a_list = _elementary_index_set(harmonics, p_max) + p_values = sorted(set(p for (p, n) in a_list)) + + # --- derivative-weighted template families (per p), built once, reused for all n,det --- + # (the derivative is a trivial FD weight; no waveform regeneration) + hlms_p = {p: {lm: fd_apply_time_derivative(hlms[lm], p) for lm in hlms} for p in p_values} + hlms_conj_p = {p: {lm: fd_apply_time_derivative(hlms_conj[lm], p) for lm in hlms_conj} + for p in p_values} + + # For U,V we need the modulated templates chi_a (modulation cannot be pushed onto data + # in a overlap). Build them once per a (cheap FD op). + # + # CRITICAL reference-time note: the template time series carry the INTRINSIC epoch + # (hlms.epoch ~ -T_dur, i.e. near 0), NOT the absolute event time t_ev ~ 1e9. When the + # template is placed at the event, sample j sits at absolute time t' = t_ev + (hlms.epoch + # + j*dt), so the physical modulation exp(i n Omega (t' - t_ev)) = exp(i n Omega + # (hlms.epoch + j*dt)) -- i.e. reference the template modulation to 0, NOT to t_ev. + # (Referencing to t_ev with the tiny template epoch would apply exp(i n Omega * ~-1e9), + # a ~1e4 rad spurious phase that randomizes U,V and inflates lnL beyond 0.5.) + # The data-term route above keeps t_ev because the DATA epoch really is ~ t_ev. + chi = {a: {lm: fd_apply_sidereal_modulation(hlms_p[a[0]][lm], a[1], f_sidereal, 0.0) + for lm in hlms} for a in a_list} + chi_conj = {a: {lm: fd_apply_sidereal_modulation(hlms_conj_p[a[0]][lm], a[1], f_sidereal, 0.0) + for lm in hlms_conj} for a in a_list} + + rholms_rot = {} + rholms_intp_rot = {} + crossTerms_rot = {} + crossTermsV_rot = {} + + for det in detectors: + psd = psd_dict[det] + data = data_dict[det] + t_det = FL.ComputeArrivalTimeAtDetector(det, P.phi, P.theta, event_time_geo) + rho_epoch = data.epoch - hlms[list(hlms.keys())[0]].epoch + t_shift = float(float(t_det) - float(t_window) - float(rho_epoch)) + N_shift = int(t_shift / P.deltaT + 0.5) + N_window = int(2 * t_window / P.deltaT) + t = np.arange(N_window) * P.deltaT + float(rho_epoch + N_shift * P.deltaT) + + # ---- data-term overlaps Q^a_lm(t) ---- + # exp(i n Omega t) on the template is equivalent to shifting the data spectrum by + # -n f_sidereal (mode-independent). Realize it by modulating the DATA time series + # by exp(-i n Omega (t_abs - t_ev)) (round trip). Because the modulation lives on + # the absolute data-time axis, the resulting overlap is directly + # Q^a_lm(t) = int e^{-i n Omega (t'-t_ev)} [d^p h_lm]^*(t'-t) d(t') dt' + # with NO arrival-time-dependent post-phase. + rholms_rot[det] = {} + rholms_intp_rot[det] = {} + data_by_n = {} + for n in set(nn for (_, nn) in a_list): + data_by_n[n] = data if n == 0 else _lal_freq_modulate(data, -n, f_sidereal, t_ev) + + for a in a_list: + p, n = a + rho = FL.ComputeModeIPTimeSeries( + hlms_p[p], data_by_n[n], psd, P.fmin, fMax, 1. / 2. / P.deltaT, + N_shift, N_window, analyticPSD_Q, inv_spec_trunc_Q, T_spec) + rholms_rot[det][a] = rho + if not skip_interpolation: + rholms_intp_rot[det][a] = FL.InterpolateRholms(rho, t, verbose=verbose) + else: + rholms_intp_rot[det][a] = None + + # ---- cross terms U,V for each ordered pair (a,a') ---- + crossTerms_rot[det] = {} + crossTermsV_rot[det] = {} + for a in a_list: + for ap in a_list: + crossTerms_rot[det][(a, ap)] = FL.ComputeModeCrossTermIP( + chi[a], chi[ap], psd, P.fmin, fMax, 1. / 2. / P.deltaT, P.deltaF, + analyticPSD_Q, inv_spec_trunc_Q, T_spec, verbose=False, + same_waveform_Q=False) + crossTermsV_rot[det][(a, ap)] = FL.ComputeModeCrossTermIP( + chi_conj[a], chi[ap], psd, P.fmin, fMax, 1. / 2. / P.deltaT, P.deltaF, + analyticPSD_Q, inv_spec_trunc_Q, T_spec, prefix="V", verbose=False, + same_waveform_Q=False) + + meta = dict(harmonics=tuple(harmonics), p_max=p_max, f_sidereal=f_sidereal, + a_list=a_list, event_time_geo=float(event_time_geo), + omega_earth=OMEGA_EARTH, modes=list(hlms.keys())) + return rholms_intp_rot, crossTerms_rot, crossTermsV_rot, rholms_rot, meta + + +# --------------------------------------------------------------------------- +# Rotation-aware log-likelihood assembly (Path A: amplitude drift, p_max=0). +# --------------------------------------------------------------------------- +def antenna_harmonics_tilde(det, RA, DEC, psi, tref): + """Return {n: A_tilde_n} where F_k(t) = sum_n A_tilde_n exp(i n Omega (t - tref)). + + A_tilde_n = A_n(response, dec, psi) * exp(i n (GMST(tref) - RA)), with A_n the + RA/time-independent antenna harmonics from slowrot_response. So A_tilde_n is the + coefficient of the precompute's modulation exp(i n Omega (t - tref)); sum_n A_tilde_n + = F_k(tref) reproduces lal.ComputeDetAMResponse exactly. + """ + import lal + import lalsimulation as lalsim + from . import slowrot_response as srr + lald = lalsim.DetectorPrefixToLALDetector(det) + A = srr.antenna_harmonics(lald.response, DEC, psi) + g_ev = float(lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(float(tref)))) - RA + return {n: A[n] * np.exp(1.0j * n * g_ev) for n in A} + + +def _convolve_harmonics(a, b): + """Convolve two harmonic sequences (dicts {m: coef}) -> dict {m: coef}.""" + out = {} + for m1, c1 in a.items(): + for m2, c2 in b.items(): + out[m1 + m2] = out.get(m1 + m2, 0j) + c1 * c2 + return out + + +def rotation_coefficients(det, RA, DEC, psi, tref, p_max): + """Coefficients {(p, ntilde): C} of the elementary modulated templates + chi_{(p,n)} = exp(i n Omega (t-tref)) h^{(p)}(t-tau0) in the folded response+delay + template (Path B). For p_max=0 this is {(0,n): A_tilde_n} (Path A). + + C_{(p, ntilde)} = (1/p!) sum_{n+m=ntilde} A_tilde_n [(-D)^{*p}]_m + with A_tilde_n the antenna harmonics and D_m the delay-DRIFT harmonics: + delta_tau(t) = tau(t) - tau(tref) = sum_m D_m exp(i m Omega (t-tref)), + D_0 = -(B_tilde_1 + B_tilde_-1), D_{+-1} = B_tilde_{+-1} (B_tilde = delay harmonics). + At Omega->0 delta_tau=0 -> D=0 -> C_{(p,n)}=0 for p>=1, so Path B -> Path A. + """ + import math + import lal + import lalsimulation as lalsim + from . import slowrot_response as srr + lald = lalsim.DetectorPrefixToLALDetector(det) + g_ev = float(lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(float(tref)))) - RA + A = srr.antenna_harmonics(lald.response, DEC, psi) + Atil = {n: A[n] * np.exp(1.0j * n * g_ev) for n in A} + if p_max == 0: + return {(0, n): Atil[n] for n in Atil} + Bd = srr.delay_harmonics(lald.location, DEC) # {m: B_m}, m in -1,0,1 + Btil = {m: Bd[m] * np.exp(1.0j * m * g_ev) for m in Bd} + tau0 = np.real(sum(Btil.values())) # tau(tref) + D = {m: Btil[m] for m in Btil} + D[0] = D[0] - tau0 # = -(Btil_1 + Btil_-1) + negD = {m: -D[m] for m in D} + C = {} + E = {0: 1.0 + 0j} # (-D)^{*0} + for p in range(p_max + 1): + if p > 0: + E = _convolve_harmonics(E, negD) # (-D)^{*p} + inv = 1.0 / math.factorial(p) + for n, an in Atil.items(): + for m, em in E.items(): + key = (p, n + m) + C[key] = C.get(key, 0j) + inv * an * em + return C + + +def FactoredLogLikelihoodWithRotation(extr_params, rholms_intp_rot, crossTerms_rot, + crossTermsV_rot, meta, Lmax): + """Slow-rotation analogue of factored_likelihood.FactoredLogLikelihood (Path A). + + Contracts the harmonic-resolved precompute bank with the antenna harmonics + A_tilde_n(det, RA, DEC, psi, tref). Reduces EXACTLY to the baseline + FactoredLogLikelihood when f_sidereal -> 0 (all modulations become identity and + sum_n A_tilde_n -> F_k(tref)). + + Currently implements p_max=0 (amplitude drift only); the delay-derivative (Path B) + contraction with B_n is a TODO. + """ + import lal + from . import factored_likelihood as FL + from .. import lalsimutils as lsu + + p_max = meta['p_max'] + a_list = list(meta['a_list']) + harmonics = list(meta['harmonics']) + hset = set(harmonics) + for n in harmonics: + assert -n in hset, "harmonic set must be symmetric (need -n for the V term): %s" % harmonics + + RA = extr_params.phi + DEC = extr_params.theta + tref = extr_params.tref + phiref = extr_params.phiref + incl = extr_params.incl + psi = extr_params.psi + dist = extr_params.dist + + detectors = list(rholms_intp_rot.keys()) + a0 = a_list[0] + modes = list(rholms_intp_rot[detectors[0]][a0].keys()) + Ylms = FL.ComputeYlms(Lmax, incl, -phiref, selected_modes=modes) + + distMpc = dist / (lsu.lsu_PC * 1e6) + invDistMpc = FL.distMpcRef / distMpc + + lnL = 0. + for det in detectors: + C = rotation_coefficients(det, RA, DEC, psi, tref, p_max) # {(p,n): C_a} + t_det = FL.ComputeArrivalTimeAtDetector(det, RA, DEC, tref) + CT = crossTerms_rot[det] + CTV = crossTermsV_rot[det] + + # Q^a_lm(t_det) for each elementary template a=(p,n) + Q = {a: {m: rholms_intp_rot[det][a][m](float(t_det)) + for m in modes} for a in a_list} + + # term1 = Re[ sum_lm conj(Ylm) sum_a conj(C_a) Q^a_lm(t_det) ] + term1 = 0. + for m in modes: + s = 0. + for a in a_list: + s += np.conj(C.get(a, 0j)) * Q[a][m] + term1 += np.conj(Ylms[m]) * s + term1 = np.real(term1) * invDistMpc + + # term2 = -1/4 Re[ sum_{p1,p2} ( U-part conj(Y1)Y2 + V-part Y1 Y2 ) ] + # U-part = sum_{a,a'} conj(C_a) C_a' U^{(a,a')}[p1,p2] + # V-part = sum_{a=(p,nu),a'} C_{(p,-nu)} C_a' V^{(a,a')}[p1,p2] + term2 = 0. + for p1 in modes: + for p2 in modes: + u = 0. + v = 0. + for a in a_list: + aR = (a[0], -a[1]) + for ap in a_list: + u += np.conj(C.get(a, 0j)) * C.get(ap, 0j) * CT[(a, ap)][(p1, p2)] + v += C.get(aR, 0j) * C.get(ap, 0j) * CTV[(a, ap)][(p1, p2)] + term2 += u * np.conj(Ylms[p1]) * Ylms[p2] + v * Ylms[p1] * Ylms[p2] + term2 = -np.real(term2) / 4. / (distMpc / FL.distMpcRef) ** 2 + + lnL += term1 + term2 + + return lnL + + +# --------------------------------------------------------------------------- +# Vectorized (NoLoop) rotation likelihood -- the maintained batchmode path. +# --------------------------------------------------------------------------- +def rotation_coefficients_vector(det, RA, DEC, psi, tref, p_max): + """Vectorized rotation_coefficients: RA, DEC, psi are arrays (npts_ex,); returns + {(p, n): complex ndarray (npts_ex,)}. Same algebra as rotation_coefficients.""" + import math + import lal + import lalsimulation as lalsim + from . import slowrot_response as srr + lald = lalsim.DetectorPrefixToLALDetector(det) + RA = np.asarray(RA); DEC = np.asarray(DEC); psi = np.asarray(psi) + g_ev = float(lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(float(tref)))) - RA + A = srr.antenna_harmonics_vector(lald.response, DEC, psi) + Atil = {n: A[n] * np.exp(1.0j * n * g_ev) for n in A} + if p_max == 0: + return {(0, n): Atil[n] for n in Atil} + Bd = srr.delay_harmonics_vector(lald.location, DEC) + Btil = {m: Bd[m] * np.exp(1.0j * m * g_ev) for m in Bd} + tau0 = np.real(sum(Btil.values())) + D = {m: Btil[m] for m in Btil} + D[0] = D[0] - tau0 + negD = {m: -D[m] for m in D} + C = {} + E = {0: np.ones_like(g_ev, dtype=complex)} + for p in range(p_max + 1): + if p > 0: + E = _convolve_harmonics(E, negD) + inv = 1.0 / math.factorial(p) + for n, an in Atil.items(): + for m, em in E.items(): + key = (p, n + m) + C[key] = C.get(key, 0j) + inv * an * em + return C + + +def pack_rotation_arrays(meta, rholms_rot, crossTerms_rot, crossTermsV_rot): + """Pack the elementary-template precompute bank into dense arrays for the NoLoop path. + + Returns (lookupNKDict, rho_by_a, U_by_aa, V_by_aa, epochDict), keyed per detector by + elementary template a=(p,n) (Path A: a=(0,n); Path B: also p>=1). + """ + a_list = list(meta['a_list']) + lookupNKDict = {}; rho_by_a = {}; U_by_aa = {}; V_by_aa = {}; epochDict = {} + for det in rholms_rot: + a0 = a_list[0] + modes = list(rholms_rot[det][a0].keys()) + n_lms = len(modes) + idx = {m: i for i, m in enumerate(modes)} + lookupNKDict[det] = np.array([[m[0], m[1]] for m in modes], dtype=int) + npts = rholms_rot[det][a0][modes[0]].data.length + epochDict[det] = float(rholms_rot[det][a0][modes[0]].epoch) + rho_by_a[det] = {} + for a in a_list: + arr = np.zeros((n_lms, npts), dtype=np.complex128) + for m in modes: + arr[idx[m]] = rholms_rot[det][a][m].data.data + rho_by_a[det][a] = arr + U_by_aa[det] = {}; V_by_aa[det] = {} + for a in a_list: + for ap in a_list: + Um = np.zeros((n_lms, n_lms), dtype=np.complex128) + Vm = np.zeros((n_lms, n_lms), dtype=np.complex128) + cU = crossTerms_rot[det][(a, ap)] + cV = crossTermsV_rot[det][(a, ap)] + for m1 in modes: + for m2 in modes: + Um[idx[m1], idx[m2]] = cU[(m1, m2)] + Vm[idx[m1], idx[m2]] = cV[(m1, m2)] + U_by_aa[det][(a, ap)] = Um + V_by_aa[det][(a, ap)] = Vm + return lookupNKDict, rho_by_a, U_by_aa, V_by_aa, epochDict + + +def DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation( + tvals, P_vec, meta, lookupNKDict, rho_by_a, U_by_aa, V_by_aa, epochDict, + Lmax=2, array_output=False, time_interp='nearest'): + """Vectorized rotation-aware lnL (Path A). + + Mirrors factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopOrig with + sums over sidereal harmonics and the per-sample antenna harmonics + A_tilde_n = A_n(det,dec,psi) exp(i n (GMST(tref)-RA)). Extrinsic params in P_vec are + ARRAYS of length npts_ex (RA, DEC, incl, phiref, psi, dist); tref, deltaT scalar. + + array_output=True returns lnL_t of shape (npts_ex, npts) (before time marginalization); + array_output=False returns the time-marginalized lnL of shape (npts_ex,). + """ + import lal + from . import factored_likelihood as FL + + a_list = list(meta['a_list']) + p_max = meta['p_max'] + npts = len(tvals); npts_ex = len(P_vec.phi) + RA = P_vec.phi; DEC = P_vec.theta + incl = P_vec.incl; phiref = P_vec.phiref; psi = P_vec.psi + distMpc = P_vec.dist / (lal.PC_SI * 1e6) + lnL_t = np.zeros((npts_ex, npts), dtype=np.float64) + + for det in rho_by_a: + n_lms = len(lookupNKDict[det]) + Ylms = FL.ComputeYlmsArrayVector(lookupNKDict[det], incl, -phiref).T # (npts_ex, n_lms) + C = rotation_coefficients_vector(det, RA, DEC, psi, P_vec.tref, p_max) # {(p,n): (npts_ex,)} + zeroC = np.zeros(npts_ex, dtype=complex) + + def Cg(a): + return C[a] if a in C else zeroC + t_ref = epochDict[det] + # Match the maintained baseline NoLoop's precision-preserving time reference: keep the + # tiny (tref - epoch) offset and the small geometric delay separate, rather than + # subtracting two ~1e9 s absolute arrival times. The absolute-difference form loses + # ~1e-3 of a sample bin, which is harmless for nearest-neighbour snapping but shows up + # directly in the sub-bin fraction used by cubic interpolation. + detector_location = np.asarray(FL.lalsim.DetectorPrefixToLALDetector(det).location) + gmst_tref = float(lal.GreenwichMeanSiderealTime(P_vec.tref)) + t_det = float(P_vec.tref - float(t_ref)) + FL.TimeDelayFromEarthCenter( + detector_location, RA, DEC, gmst_tref) + sample_first = (t_det + tvals[0]) / P_vec.deltaT + if time_interp == 'nearest': + ifirst = (np.round(sample_first) + 0.5).astype(int) + else: + ifirst = np.floor(sample_first).astype(int) + frac_first = (sample_first - np.floor(sample_first)).astype(np.float64) + ilast = ifirst + npts + + term1 = np.zeros((npts_ex, npts), dtype=np.complex128) + for a in a_list: + det_rho = rho_by_a[det][a] + if time_interp == 'nearest': + Qa = np.empty((npts_ex, npts, n_lms), dtype=np.complex128) + for i in range(npts_ex): + Qa[i] = det_rho[..., ifirst[i]:ilast[i]].T + else: + # cubic sub-sample interpolation (calmarg time_interp='cubic'): + # _cubic_Q_window_numpy expects Q_block shape (n_time, n_lm). + Qa = FL._cubic_Q_window_numpy(det_rho.T, ifirst, frac_first, npts) + term1 += np.conj(Cg(a))[:, None] * np.einsum('xi,xti->xt', np.conj(Ylms), Qa) + term1 = term1.real * (FL.distMpcRef / distMpc)[:, None] + + term2 = np.zeros(npts_ex, dtype=np.complex128) + for a in a_list: + aR = (a[0], -a[1]) + for ap in a_list: + term2 += np.conj(Cg(a)) * Cg(ap) * np.einsum( + 'xi,xj,ij->x', np.conj(Ylms), Ylms, U_by_aa[det][(a, ap)]) + term2 += Cg(aR) * Cg(ap) * np.einsum( + 'xi,xj,ij->x', Ylms, Ylms, V_by_aa[det][(a, ap)]) + term2 = (-0.25 * term2.real) * (FL.distMpcRef / distMpc) ** 2 + + lnL_t += term1 + term2[:, None] + + if array_output: + return lnL_t + lnLmax = np.max(lnL_t, axis=-1, keepdims=True) + L = FL.my_simps(np.exp(lnL_t - lnLmax), dx=P_vec.deltaT, axis=-1) + return lnLmax[..., 0] + np.log(L) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/slowrot_freqresponse.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/slowrot_freqresponse.py new file mode 100644 index 000000000..a9950d275 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/slowrot_freqresponse.py @@ -0,0 +1,442 @@ +""" +slowrot_freqresponse : closed-form FREQUENCY-DEPENDENT (finite-size) detector antenna +response of a ground-based Michelson interferometer, beyond the long-wavelength +approximation. "Path D" of the slow-rotation generalization of the RIFT likelihood. + +Motivation +---------- +The usual antenna response F_+(RA,DEC,psi), F_x(RA,DEC,psi) returned by +lal.ComputeDetAMResponse is the LONG-WAVELENGTH LIMIT: it assumes the GW strain is +uniform across the detector while a photon traverses an arm. Beyond that limit the +finite light-travel time T = L/c across an arm imprints genuine in-band frequency +structure on the response, controlled by the free spectral range + + f_FSR = c / (2 L) ( 37.5 kHz for 4-km LIGO ; 3.75 kHz for a 40-km CE arm ). + +The response becomes F_A -> F_A(f; RA, DEC, psi). This matters for 3G detectors +(Cosmic Explorer, 40 km) whose long signals reach a non-negligible fraction of f_FSR +in band; it is completely negligible for the 4-km LIGO instruments. + +Formula and convention (arXiv:2412.01693, "Beyond the Long Wavelength Approximation", +Eqs. 4-6; equivalent to Rakhmanov-Romano-Whelan 2008, arXiv:0808.3805, and +Essick-Vitale-Evans 2017, PRD 96 084004 / arXiv:1708.06843) +------------------------------------------------------------------------------------- +Single-arm round-trip transfer function, for arm unit vector a-hat and source direction +n-hat (direction TO the source, Sathyaprakash-Schutz convention -- the SAME n-hat used +by lal.ComputeDetAMResponse / TimeDelayFromEarthCenter), with T = L/c the one-way +photon time of flight and T_pm^a = T (1 +/- a-hat . n-hat): + + D~(a-hat, n-hat, f) = (e^{-i 2 pi f T} / 2) + * [ e^{+i pi f T_-^a} sinc(pi f T_+^a) + + e^{-i pi f T_+^a} sinc(pi f T_-^a) ] (Eq. 6) + +with sinc(z) = sin(z)/z, so D~ -> 1 as f -> 0. The frequency-dependent antenna +patterns are then (Eqs. 4-5) + + F_A(f, n-hat) = (1/2) [ D~(x-hat, n-hat, f) (eps^A : x-hat x-hat) + - D~(y-hat, n-hat, f) (eps^A : y-hat y-hat) ] , A in {+, x} + +where eps^+ , eps^x are the GW polarization tensors set by (RA, DEC, psi) and +x-hat, y-hat are the two arm unit vectors. At f -> 0 both D~ -> 1 and this collapses to + F_A(0) = eps^A : (1/2)(x-hat x-hat - y-hat y-hat) = eps^A : D_response + = lal.ComputeDetAMResponse , +since the detector response tensor is exactly D_response = (1/2)(x-hat x-hat - y-hat y-hat). + +Implementation note (exactness of the f -> 0 limit) +--------------------------------------------------- +LAL stores the detector response tensor `d.response` in single precision (REAL4). +Reconstructing (1/2)(x x - y y) from the geodetic arm geometry therefore agrees with +`d.response` only to ~1e-7. To reproduce lal.ComputeDetAMResponse to MACHINE PRECISION +at f = 0 we write the response as an EXACT long-wavelength baseline plus a finite-size +correction that vanishes identically at f = 0: + + F_A(f) = F_A^LWL + (1/2)[ (D~_x - 1)(eps^A : x x) - (D~_y - 1)(eps^A : y y) ] + +The baseline F_A^LWL is computed from the SAME `d.response` and triad algebra as +lal.ComputeDetAMResponse (see slowrot_response.py), so at f = 0 (where D~_x = D~_y = 1 +EXACTLY, sinc(0)=1) the correction is exactly zero and F_A(0) == ComputeDetAMResponse. +The tiny (~1e-7) geodetic imprecision only enters the O((f/f_FSR)^2) correction, i.e. +utterly negligibly. + +Pure numpy + lal; importable without the heavy RIFT stack (mirrors slowrot_response). + +Conventions (identical to slowrot_response.py / vectorized_lal_tools.py): + g = GMST - RA (Greenwich hour angle) + nhat = (cos_dec cos g, -cos_dec sin g, sin_dec) (direction to source) + X = (-cp sg - sp cg sd, -cp cg + sp sg sd, sp cd) (polarization triad) + Y = ( sp sg - cp cg sd, sp cg + cp sg sd, cp cd) + F_+ = X.D.X - Y.D.Y , F_x = X.D.Y + Y.D.X + (cp,sp)=(cos,sin)psi ; (cd,sd)=(cos,sin)dec ; (cg,sg)=(cos,sin)g +""" +from __future__ import print_function, division + +import numpy as np + +try: + import lalsimulation as _lalsim + _HAVE_LAL = True +except Exception: # pragma: no cover + _HAVE_LAL = False + +C_SI = 299792458.0 # m/s, matches vectorized_lal_tools / slowrot_response +DEFAULT_L_LIGO = 3994.5 # m, nominal LIGO arm length (override for CE etc.) + + +# --------------------------------------------------------------------------- +# Detector geometry: arm unit vectors and arm length from a LAL detector. +# --------------------------------------------------------------------------- +def _cartesian_arm(cosAlt, sinAlt, cosAz, sinAz, cosLat, sinLat, cosLon, sinLon): + """Earth-fixed Cartesian unit vector of a detector arm. + + Exact port of LAL's getCartesianComponents (LALDetectors.c): azimuth measured + from local North toward East, altitude above the local horizontal. + """ + uNorth = cosAlt * cosAz + uEast = cosAlt * sinAz + uRho = -sinLat * uNorth + cosLat * sinAlt + return np.array([ + cosLon * uRho - sinLon * uEast, + sinLon * uRho + cosLon * uEast, + cosLat * uNorth + sinLat * sinAlt, + ]) + + +def detector_geometry(det, L_arm=None): + """Return (response, x_arm, y_arm, L) for a detector prefix ('H1','L1','V1',...). + + Parameters + ---------- + det : str detector prefix understood by lalsimulation + L_arm : float or None arm length [m] override (e.g. 40000. for a 40-km CE arm); + default = 2 * xArmMidpoint from the LAL detector. + + Returns + ------- + response : (3,3) float ndarray Earth-fixed detector response tensor D (from LAL) + x_arm, y_arm : (3,) float ndarray Earth-fixed arm unit vectors + L : float arm length [m] used for the finite-size transfer + """ + if not _HAVE_LAL: # pragma: no cover + raise RuntimeError("lalsimulation is required for detector_geometry") + d = _lalsim.DetectorPrefixToLALDetector(det) + fr = d.frDetector + lat, lon = fr.vertexLatitudeRadians, fr.vertexLongitudeRadians + cL, sL, cO, sO = np.cos(lat), np.sin(lat), np.cos(lon), np.sin(lon) + x_arm = _cartesian_arm(np.cos(fr.xArmAltitudeRadians), np.sin(fr.xArmAltitudeRadians), + np.cos(fr.xArmAzimuthRadians), np.sin(fr.xArmAzimuthRadians), + cL, sL, cO, sO) + y_arm = _cartesian_arm(np.cos(fr.yArmAltitudeRadians), np.sin(fr.yArmAltitudeRadians), + np.cos(fr.yArmAzimuthRadians), np.sin(fr.yArmAzimuthRadians), + cL, sL, cO, sO) + L = float(L_arm) if L_arm is not None else 2.0 * float(fr.xArmMidpoint) + return np.asarray(d.response, dtype=float), x_arm, y_arm, L + + +# --------------------------------------------------------------------------- +# Kinematic pieces (triad + source direction), matching ComputeDetAMResponse. +# --------------------------------------------------------------------------- +def _triad(dec, psi, g): + """Polarization triad X, Y and source direction nhat at hour angle g = GMST - RA. + + dec, psi, g may be scalars or broadcastable arrays; returned vectors carry the + 3-component along the last axis. + """ + cd, sd = np.cos(dec), np.sin(dec) + cp, sp = np.cos(psi), np.sin(psi) + cg, sg = np.cos(g), np.sin(g) + + X = np.stack([-cp * sg - sp * cg * sd, + -cp * cg + sp * sg * sd, + sp * cd * np.ones_like(sg)], axis=-1) + Y = np.stack([ sp * sg - cp * cg * sd, + sp * cg + cp * sg * sd, + cp * cd * np.ones_like(sg)], axis=-1) + nhat = np.stack([cd * cg, -cd * sg, sd * np.ones_like(sg)], axis=-1) + return X, Y, nhat + + +# --------------------------------------------------------------------------- +# Single-arm finite-size transfer function (Eq. 6 of arXiv:2412.01693). +# --------------------------------------------------------------------------- +def single_arm_transfer(a_dot_n, f, L): + """Round-trip single-arm transfer function D~(a-hat, n-hat, f). + + Parameters + ---------- + a_dot_n : float or ndarray a-hat . n-hat (arm unit vector dotted with source dir) + f : float or ndarray frequency [Hz] + L : float arm length [m] + + Returns + ------- + complex ndarray broadcast over (a_dot_n, f). D~ -> 1 as f -> 0. + + D~ = (e^{-i 2 pi f T}/2)[ e^{+i pi f T_-} sinc(pi f T_+) + e^{-i pi f T_+} sinc(pi f T_-) ] + with T = L/c, T_pm = T(1 +/- a_dot_n), sinc(z)=sin(z)/z. + (np.sinc(u) = sin(pi u)/(pi u), so sinc(pi f T_pm) = np.sinc(f T_pm).) + """ + a = np.asarray(a_dot_n, dtype=float) + f = np.asarray(f, dtype=float) + T = L / C_SI + Tp = T * (1.0 + a) # T_+ + Tm = T * (1.0 - a) # T_- + pref = np.exp(-1j * 2.0 * np.pi * f * T) / 2.0 + term_p = np.exp(1j * np.pi * f * Tm) * np.sinc(f * Tp) + term_m = np.exp(-1j * np.pi * f * Tp) * np.sinc(f * Tm) + return pref * (term_p + term_m) + + +# --------------------------------------------------------------------------- +# Frequency-dependent antenna response. +# --------------------------------------------------------------------------- +def _lwl_response(response, X, Y): + """Long-wavelength F_+, F_x from triad and response tensor (== ComputeDetAMResponse). + + F_+ = X.D.X - Y.D.Y , F_x = X.D.Y + Y.D.X . Vectorized over leading axes of X,Y. + """ + D = response + XDX = np.einsum('...i,ij,...j->...', X, D, X) + YDY = np.einsum('...i,ij,...j->...', Y, D, Y) + XDY = np.einsum('...i,ij,...j->...', X, D, Y) + YDX = np.einsum('...i,ij,...j->...', Y, D, X) + return XDX - YDY, XDY + YDX + + +def antenna_response_fd(det, ra, dec, psi, f, gmst=0.0, L_arm=None): + """Frequency-dependent antenna response F_+(f), F_x(f) for a ground-based detector. + + Beyond the long-wavelength limit: includes the finite light-travel-time transfer + across the arms. Reduces EXACTLY to lal.ComputeDetAMResponse as f -> 0. + + Parameters + ---------- + det : str detector prefix ('H1','L1','V1','K1',...) + ra, dec, psi : float right ascension, declination, polarization angle [rad] + f : float or (Nf,) ndarray frequency [Hz] (scalar or 1-D array) + gmst : float Greenwich mean sidereal time [rad] (default 0). Enters + only through g = gmst - ra, exactly as ComputeDetAMResponse. + L_arm : float or None arm-length override [m] (e.g. 40000. for 40-km CE); + default = LAL detector arm length (~4 km for LIGO). + + Returns + ------- + (Fp, Fc) : complex ndarray, same shape as f. At f=0, imag part -> 0 and the real + parts equal lal.ComputeDetAMResponse(response, ra, dec, psi, gmst). + """ + response, x_arm, y_arm, L = detector_geometry(det, L_arm=L_arm) + return antenna_response_fd_geom(response, x_arm, y_arm, L, ra, dec, psi, f, gmst=gmst) + + +def antenna_response_fd_geom(response, x_arm, y_arm, L, ra, dec, psi, f, gmst=0.0): + """Same as antenna_response_fd but with detector geometry supplied explicitly. + + Useful for exotic geometries (custom L, custom arms). See antenna_response_fd. + """ + g = gmst - ra + X, Y, nhat = _triad(dec, psi, g) # shape (3,) + + # Exact long-wavelength baseline (matches ComputeDetAMResponse to machine precision). + Fp_lwl, Fc_lwl = _lwl_response(response, X, Y) + + # Arm-basis polarization projections eps^A : (a a). + Xx, Yx = float(X @ x_arm), float(Y @ x_arm) + Xy, Yy = float(X @ y_arm), float(Y @ y_arm) + Ppx = Xx * Xx - Yx * Yx # eps^+ : x x + Ppy = Xy * Xy - Yy * Yy # eps^+ : y y + Pcx = 2.0 * Xx * Yx # eps^x : x x + Pcy = 2.0 * Xy * Yy # eps^x : y y + + # Finite-size transfer per arm (vanishes-minus-one -> 0 at f=0). + ax = float(x_arm @ nhat) + ay = float(y_arm @ nhat) + Dx = single_arm_transfer(ax, f, L) + Dy = single_arm_transfer(ay, f, L) + + Fp = Fp_lwl + 0.5 * ((Dx - 1.0) * Ppx - (Dy - 1.0) * Ppy) + Fc = Fc_lwl + 0.5 * ((Dx - 1.0) * Pcx - (Dy - 1.0) * Pcy) + return Fp, Fc + + +def free_spectral_range(L): + """Free spectral range f_FSR = c / (2 L) [Hz].""" + return C_SI / (2.0 * L) + + +# =========================================================================== +# Step 1 : sky-harmonic / power-series expansion of the finite-size response. +# =========================================================================== +# The finite-size antenna response is (pure Eqs. 4-6, no LWL split) +# +# F_A(f;sky) = (1/2)[ D~(x,n,f) (eps^A:xx) - D~(y,n,f) (eps^A:yy) ] . +# +# Factor out the common one-way light-crossing delay T = L/c that both arms +# share (DIRECTION-INDEPENDENT, degenerate with coalescence time): +# +# D~(a,n,f) = e^{-i 2 pi f T} g(f; a) , a = a-hat . n-hat , +# g(f; a) = (1/2)[ e^{+i u(1-a)} sinc(u(1+a)) + e^{-i u(1+a)} sinc(u(1-a)) ] , +# u = pi f T , sinc(x) = sin(x)/x , g(f;0) = sinc(2u) , g -> 1 as f->0. +# +# Taylor-expand the residual g in the arm projection a (the small in-band +# parameter is eps*a with eps = f L/c = f T): +# +# g(f; a) = sum_{q>=0} c_q(f) a^q . (c_q sky-INDEPENDENT) +# +# The arm projections a_x = x-hat.n, a_y = y-hat.n and the arm-basis polarization +# contractions combine into a single COMPLEX analytic sky/pol scalar per power q, +# +# beta_q(sky) = (1/2)[ zx^2 a_x^q - zy^2 a_y^q ] , +# zx = X.x-hat + i Y.x-hat , zy = X.y-hat + i Y.y-hat (X,Y polarization triad) +# +# so that G_A(f;sky) := F_A e^{+i2 pi f T} gives the complex response +# +# G(f) = G_+ + i G_x = sum_q c_q(f) beta_q(sky) , (the analogue of A_n) +# beta_0 = (1/2)(zx^2 - zy^2) = F_+^LWL + i F_x^LWL (== ComputeDetAMResponse, up +# to the ~1e-7 REAL4 geodety). +# +# The full response, EXACTLY reproducing antenna_response_fd, is then +# +# F(f) = F0_lal + e^{-i2 pi f T} G(f) - beta_0 (Fp+iFc form) +# Fbar(f) = conj(F0_lal) + e^{-i2 pi f T} Gbar(f) - conj(beta_0), Gbar=sum_q c_q conj(beta_q) +# F_+ = (F+Fbar)/2 , F_x = (F-Fbar)/(2i) +# +# with F0_lal the EXACT lal.ComputeDetAMResponse baseline (so f->0 is machine +# precise even though beta_0 from the arm triads carries the ~1e-7 geodetic error; +# the arm-based pieces enter only the finite-size correction, which vanishes at f=0). +# --------------------------------------------------------------------------- +from math import comb as _comb, factorial as _factorial + + +def _sinc_shift_apoly(u, s, Qmax, jmax=64): + """a-power coefficients d_m(u), m=0..Qmax, of sinc(u(1 + s a)) (s = +/-1). + + Entire-series form (stable at u=0): sinc(x) = sum_j (-1)^j x^{2j}/(2j+1)!, + x = u(1+s a) => coeff of a^m is s^m sum_{j>=ceil(m/2)} (-1)^j u^{2j}/(2j+1)! C(2j,m). + u may be a numpy array (frequencies); returns shape (Qmax+1,) + u.shape. + """ + u = np.asarray(u, dtype=float) + out = np.zeros((Qmax + 1,) + u.shape, dtype=float) + u2 = u * u + # precompute u^{2j} + upow = np.ones_like(u) # u^{0} + for j in range(0, jmax + 1): + coef = ((-1) ** j) / float(_factorial(2 * j + 1)) + for m in range(0, min(Qmax, 2 * j) + 1): + out[m] += coef * _comb(2 * j, m) * upow + upow = upow * u2 + for m in range(Qmax + 1): + out[m] *= float(s) ** m + return out.astype(complex) + + +def _poly_mul(A, B, Qmax): + """Truncated product (to order Qmax) of two a-power stacks A[k],B[k] over freq axis.""" + out = np.zeros_like(A[:Qmax + 1]) + for i in range(min(len(A), Qmax + 1)): + for j in range(min(len(B), Qmax + 1 - i)): + out[i + j] += A[i] * B[j] + return out + + +def finite_size_c_coeffs(f, L, Qmax): + """Sky-independent frequency basis c_q(f), q = 0..Qmax, of g(f;a)=sum_q c_q(f) a^q. + + Parameters + ---------- + f : ndarray frequency [Hz] (signed values allowed; c_q(-f)=conj(c_q(f))) + L : float arm length [m] + Qmax : int highest power of the arm projection a retained + + Returns + ------- + c : complex ndarray, shape (Qmax+1,) + f.shape. c_0(f)=sinc(2 pi f T), + c_{q>=1} carry the finite-size shape distortion; all -> delta_{q0} as f->0. + """ + f = np.asarray(f, dtype=float) + T = L / C_SI + u = np.pi * f * T + # exp(-i u a): coeff of a^k is (-i u)^k / k! + exp_neg = np.zeros((Qmax + 1,) + u.shape, dtype=complex) + term = np.ones_like(u, dtype=complex) + for k in range(Qmax + 1): + exp_neg[k] = term / float(_factorial(k)) + term = term * (-1j * u) + sinc_p = _sinc_shift_apoly(u, +1.0, Qmax) # sinc(u(1+a)) + sinc_m = _sinc_shift_apoly(u, -1.0, Qmax) # sinc(u(1-a)) + eiu = np.exp(1j * u) + term1 = eiu * _poly_mul(exp_neg, sinc_p, Qmax) # e^{iu} e^{-iua} sinc(u(1+a)) + term2 = np.conj(eiu) * _poly_mul(exp_neg, sinc_m, Qmax) # e^{-iu} e^{-iua} sinc(u(1-a)) + return 0.5 * (term1 + term2) + + +def finite_size_geometry(det, ra, dec, psi, gmst=0.0, L_arm=None): + """Geometric scalars for the finite-size expansion at (det, ra, dec, psi, gmst). + + Returns dict with: + T = L/c common one-way delay [s] + L arm length [m] + ax, ay arm projections x-hat.n, y-hat.n + zx, zy = X.a-hat + i Y.a-hat (complex pol/arm scalars) + F0 = Fp_lwl + i Fc_lwl (EXACT lal baseline response, machine precise f=0) + """ + response, x_arm, y_arm, L = detector_geometry(det, L_arm=L_arm) + g = gmst - ra + X, Y, nhat = _triad(dec, psi, g) + Fp_lwl, Fc_lwl = _lwl_response(response, X, Y) + Xx, Yx = float(X @ x_arm), float(Y @ x_arm) + Xy, Yy = float(X @ y_arm), float(Y @ y_arm) + zx = Xx + 1j * Yx + zy = Xy + 1j * Yy + ax = float(x_arm @ nhat) + ay = float(y_arm @ nhat) + return dict(T=L / C_SI, L=L, ax=ax, ay=ay, zx=zx, zy=zy, + F0=complex(Fp_lwl) + 1j * complex(Fc_lwl)) + + +def finite_size_beta(geom, Qmax): + """Analytic sky/pol coefficients beta_q = (1/2)[zx^2 a_x^q - zy^2 a_y^q], q=0..Qmax.""" + zx2, zy2, ax, ay = geom['zx'] ** 2, geom['zy'] ** 2, geom['ax'], geom['ay'] + return np.array([0.5 * (zx2 * ax ** q - zy2 * ay ** q) for q in range(Qmax + 1)], + dtype=complex) + + +def F_fd_expanded(det, ra, dec, psi, f, Qmax, gmst=0.0, L_arm=None): + """Order-Qmax reconstruction of the finite-size response F_+(f), F_x(f). + + Converges to antenna_response_fd(det,ra,dec,psi,f,...) as Qmax -> infinity. + Uses the EXACT lal baseline F0 for the constant part and the power-series + finite-size correction e^{-i2 pi f T}(sum_q c_q beta_q) - beta_0. + """ + f = np.asarray(f, dtype=float) + geom = finite_size_geometry(det, ra, dec, psi, gmst=gmst, L_arm=L_arm) + beta = finite_size_beta(geom, Qmax) + c = finite_size_c_coeffs(f, geom['L'], Qmax) # (Qmax+1,)+f.shape + G = np.tensordot(beta, c, axes=(0, 0)) # sum_q beta_q c_q(f) + Gbar = np.tensordot(np.conj(beta), c, axes=(0, 0)) + phase = np.exp(-1j * 2.0 * np.pi * f * geom['T']) + F = geom['F0'] + phase * G - beta[0] + Fbar = np.conj(geom['F0']) + phase * Gbar - np.conj(beta[0]) + Fp = 0.5 * (F + Fbar) + Fc = (F - Fbar) / (2.0j) + return Fp, Fc + + +def finite_size_response_weights(fvals, geom, Qmax): + """Per-basis frequency weights W_p(f) folded into the FD modes for the likelihood. + + Response basis (p = 0..Qmax+1) reproducing F(f) = sum_p b_p W_p(f): + p=0 ("baseline") : W_0(f) = 1 b_0 = F0 (exact lal) + p=1+q : W_{1+q}(f) = e^{-i2pi f T} c_q(f) - [q==0] + b_{1+q} = beta_q (arm) + Each W_p is Hermitian (W_p(-f)=conj(W_p(f))) so the V cross term needs NO + harmonic reflection. The common delay e^{-i2 pi f T} (= a T=L/c arrival-time + shift of the finite-size correction relative to the LWL baseline) is carried + inside the correction weights. Returns (weights (Npbasis, Nf) complex, coeff-builder). + """ + fvals = np.asarray(fvals, dtype=float) + c = finite_size_c_coeffs(fvals, geom['L'], Qmax) + phase = np.exp(-1j * 2.0 * np.pi * fvals * geom['T']) + W = np.empty((Qmax + 2, fvals.shape[0]), dtype=complex) + W[0] = 1.0 + for q in range(Qmax + 1): + W[1 + q] = phase * c[q] - (1.0 if q == 0 else 0.0) + return W diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/slowrot_response.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/slowrot_response.py new file mode 100644 index 000000000..66220571d --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/slowrot_response.py @@ -0,0 +1,271 @@ +""" +slowrot_response : closed-form sidereal-harmonic decomposition of the ground-based +detector response, for the "slow-rotation" generalization of the RIFT likelihood. + +Over a long signal the Earth rotates, so the antenna response and the geometric +propagation delay drift. Both are exactly band-limited in harmonics of the Greenwich +hour angle g(t) = GMST(t) - RA: + + F(t) = F_+(t) + i F_x(t) = sum_{n=-2}^{+2} A_n exp(i n g(t)) (5 harmonics) + tau(t) = -(1/c) r_det . nhat(t) = sum_{n=-1}^{+1} B_n exp(i n g(t)) (3 harmonics) + +The key structural facts (see notes/ sec_amplitude, sec_delay, and the response appendix): + + * The complex antenna harmonics A_n depend ONLY on the detector response tensor, the + source declination, and the polarization angle -- NOT on right ascension or time. + RA and time enter only through the phase exp(i n g(t)) = exp(i n (GMST(t)-RA)). + Since GMST(t) ~ GMST(t_ref) + Omega (t - t_ref), this is exactly the sidereal-harmonic + modulation exp(i n Omega t) of the design notes, with sky/time in analytic scalars. + + * The delay harmonics B_n depend only on the detector location and the declination. + +This module derives A_n, B_n from the SAME algebra used by +RIFT.likelihood.vectorized_lal_tools.ComputeDetAMResponse / TimeDelayFromEarthCenter, so +the reconstruction + + F(t) = sum_n A_n exp(i n g) matches lal.ComputeDetAMResponse + tau(t) = sum_n B_n exp(i n g) matches lal.TimeDelayFromEarthCenter + +to machine precision for any g. This is an exact algebraic identity in g, not an +approximation; the ONLY approximation in the slow-rotation program is the near-linearity +of GMST(t) (handled elsewhere) and the truncation of the delay-derivative expansion +(Path B), neither of which enters here. + +Conventions mirror vectorized_lal_tools.py exactly: + g = GMST - RA (Greenwich hour angle) + nhat = (cos_dec cos g, -cos_dec sin g, sin_dec) + X, Y polarization triads as in ComputeDetAMResponse + F_+ = X.D.X - Y.D.Y, F_x = X.D.Y + Y.D.X (D = detector response tensor, symmetric) +""" +from __future__ import print_function, division + +import numpy as np + +C_SI = 299792458.0 # m/s, matches vectorized_lal_tools + + +# --------------------------------------------------------------------------- +# Constant polarization-triad basis vectors (independent of hour angle g). +# +# Writing X = Xc cos g + Xs sin g + X0 (and likewise Y), read off directly from +# vectorized_lal_tools.ComputeDetAMResponse: +# X[0] = -cos_psi sin_g - sin_psi cos_g sin_dec +# X[1] = -cos_psi cos_g + sin_psi sin_g sin_dec +# X[2] = sin_psi cos_dec +# Y[0] = sin_psi sin_g - cos_psi cos_g sin_dec +# Y[1] = sin_psi cos_g + cos_psi sin_g sin_dec +# Y[2] = cos_psi cos_dec +# --------------------------------------------------------------------------- +def _triad_components(dec, psi): + """Return the six constant 3-vectors (Xc, Xs, X0, Yc, Ys, Y0).""" + cd, sd = np.cos(dec), np.sin(dec) + cp, sp = np.cos(psi), np.sin(psi) + + Xc = np.array([-sp * sd, -cp, 0.0]) # coeff of cos g in X + Xs = np.array([-cp, sp * sd, 0.0]) # coeff of sin g in X + X0 = np.array([ 0.0, 0.0, sp * cd]) + + Yc = np.array([-cp * sd, sp, 0.0]) # coeff of cos g in Y + Ys = np.array([ sp, cp * sd, 0.0]) # coeff of sin g in Y + Y0 = np.array([ 0.0, 0.0, cp * cd]) + return Xc, Xs, X0, Yc, Ys, Y0 + + +def antenna_harmonics(response_matrix, dec, psi): + """Complex antenna-pattern harmonics A_n, n = -2..2. + + F(t) = F_+(t) + i F_x(t) = sum_{n=-2}^{2} A_n exp(i n g), g = GMST(t) - RA. + + Parameters + ---------- + response_matrix : (3,3) array detector response tensor (Earth-fixed), symmetric + dec, psi : float declination and polarization angle [rad] + + Returns + ------- + A : dict {n: complex} for n in (-2,-1,0,1,2). Independent of RA and time. + """ + D = np.asarray(response_matrix, dtype=float) + Xc, Xs, X0, Yc, Ys, Y0 = _triad_components(dec, psi) + + def B(u, v): + return float(u @ D @ v) # bilinear form; D symmetric => B(u,v)=B(v,u) + + # --- F_+ = X.D.X - Y.D.Y, expanded in cos/sin of (n g) --- + # Using cos^2 = 1/2 + 1/2 cos2g, sin^2 = 1/2 - 1/2 cos2g, cos sin = 1/2 sin2g. + Pp0 = (0.5 * (B(Xc, Xc) + B(Xs, Xs)) + B(X0, X0) + - (0.5 * (B(Yc, Yc) + B(Ys, Ys)) + B(Y0, Y0))) + Pp1 = 2.0 * (B(Xc, X0) - B(Yc, Y0)) # cos g + Qp1 = 2.0 * (B(Xs, X0) - B(Ys, Y0)) # sin g + Pp2 = 0.5 * ((B(Xc, Xc) - B(Xs, Xs)) - (B(Yc, Yc) - B(Ys, Ys))) # cos 2g + Qp2 = B(Xc, Xs) - B(Yc, Ys) # sin 2g + + # --- F_x = X.D.Y + Y.D.X = 2 X.D.Y (D symmetric) --- + Pc0 = B(Xc, Yc) + B(Xs, Ys) + 2.0 * B(X0, Y0) + Pc1 = 2.0 * (B(Xc, Y0) + B(X0, Yc)) # cos g + Qc1 = 2.0 * (B(Xs, Y0) + B(X0, Ys)) # sin g + Pc2 = B(Xc, Yc) - B(Xs, Ys) # cos 2g + Qc2 = B(Xc, Ys) + B(Xs, Yc) # sin 2g + + # Complex combine: F = F_+ + i F_x. For harmonic n>=1, the real signal piece + # P_n cos n g + Q_n sin n g -> A_n exp(i n g) + A_{-n} exp(-i n g) + # with A_{+n} = (P_n - i Q_n)/2, A_{-n} = (P_n + i Q_n)/2, and here P_n, Q_n are + # THEMSELVES complex (P_n = Pp_n + i Pc_n, Q_n = Qp_n + i Qc_n). + P0 = Pp0 + 1j * Pc0 + P1 = Pp1 + 1j * Pc1 + Q1 = Qp1 + 1j * Qc1 + P2 = Pp2 + 1j * Pc2 + Q2 = Qp2 + 1j * Qc2 + + A = { + 0: P0, + 1: 0.5 * (P1 - 1j * Q1), + -1: 0.5 * (P1 + 1j * Q1), + 2: 0.5 * (P2 - 1j * Q2), + -2: 0.5 * (P2 + 1j * Q2), + } + return A + + +def antenna_harmonics_vector(response_matrix, dec, psi): + """Vectorized antenna_harmonics: dec, psi arrays (broadcast shape S). + + Returns {n: complex ndarray shape S} for n in -2..2. Mirrors antenna_harmonics + exactly (same algebra), but dec/psi may be numpy arrays (or plain Python floats/ + 0-d arrays) of a common broadcastable shape S. + """ + D = np.asarray(response_matrix, dtype=float) + dec = np.asarray(dec, dtype=float) + psi = np.asarray(psi, dtype=float) + cd, sd = np.cos(dec), np.sin(dec) + cp, sp = np.cos(psi), np.sin(psi) + + # np.broadcast(dec, psi) misbehaves/raises for plain Python floats or 0-d arrays + # in some numpy versions (it wants array-like objects, and its .shape attribute + # is not the intuitive broadcast shape for scalar inputs). np.broadcast_shapes + # operates purely on shapes and works uniformly for scalars (shape ()), 0-d + # arrays, and full arrays. + out_shape = np.broadcast_shapes(np.shape(dec), np.shape(psi)) + zero = np.zeros(out_shape) + + def vec(a, b, c): + return np.stack(np.broadcast_arrays(a, b, c), axis=-1) + + Xc = vec(-sp * sd, -cp, zero) + Xs = vec(-cp, sp * sd, zero) + X0 = vec(zero, zero, sp * cd) + + Yc = vec(-cp * sd, sp, zero) + Ys = vec( sp, cp * sd, zero) + Y0 = vec(zero, zero, cp * cd) + + def B(u, v): + return np.einsum('...i,ij,...j->...', u, D, v) + + Pp0 = 0.5 * (B(Xc, Xc) + B(Xs, Xs)) + B(X0, X0) - (0.5 * (B(Yc, Yc) + B(Ys, Ys)) + B(Y0, Y0)) + Pp1 = 2.0 * (B(Xc, X0) - B(Yc, Y0)) + Qp1 = 2.0 * (B(Xs, X0) - B(Ys, Y0)) + Pp2 = 0.5 * ((B(Xc, Xc) - B(Xs, Xs)) - (B(Yc, Yc) - B(Ys, Ys))) + Qp2 = B(Xc, Xs) - B(Yc, Ys) + + Pc0 = B(Xc, Yc) + B(Xs, Ys) + 2.0 * B(X0, Y0) + Pc1 = 2.0 * (B(Xc, Y0) + B(X0, Yc)) + Qc1 = 2.0 * (B(Xs, Y0) + B(X0, Ys)) + Pc2 = B(Xc, Yc) - B(Xs, Ys) + Qc2 = B(Xc, Ys) + B(Xs, Yc) + + P0 = Pp0 + 1j * Pc0 + P1 = Pp1 + 1j * Pc1 + Q1 = Qp1 + 1j * Qc1 + P2 = Pp2 + 1j * Pc2 + Q2 = Qp2 + 1j * Qc2 + + # B(u,v) already contracts down to shape out_shape (einsum drops the trailing + # vector index), so no additional broadcasting of P0/P1/... is actually needed; + # multiply by ones defensively/explicitly to guarantee the advertised return shape. + ones = np.ones(out_shape) + return { + 0: P0 * ones, + 1: 0.5 * (P1 - 1j * Q1) * ones, + -1: 0.5 * (P1 + 1j * Q1) * ones, + 2: 0.5 * (P2 - 1j * Q2) * ones, + -2: 0.5 * (P2 + 1j * Q2) * ones, + } + + +def delay_harmonics(location_xyz, dec): + """Geometric-delay harmonics B_n, n = -1..1 [seconds]. + + tau(t) = -(1/c) r_det . nhat(t) = sum_{n=-1}^{1} B_n exp(i n g), g = GMST(t) - RA, + with nhat = (cos_dec cos g, -cos_dec sin g, sin_dec) as in + vectorized_lal_tools.TimeDelayFromEarthCenter. + + Parameters + ---------- + location_xyz : (3,) array detector position relative to Earth center [m], Earth-fixed + dec : float declination [rad] + + Returns + ------- + B : dict {n: complex} for n in (-1,0,1). B_0 real; B_{-1} = conj(B_1). + Independent of RA and time. + """ + r = np.asarray(location_xyz, dtype=float) + cd, sd = np.cos(dec), np.sin(dec) + + # tau = -(1/c)[ r_x cos_dec cos g - r_y cos_dec sin g + r_z sin_dec ] + # = T0 + T1c cos g + T1s sin g + T0 = -(r[2] * sd) / C_SI + T1c = -(cd * r[0]) / C_SI + T1s = +(cd * r[1]) / C_SI + + B = { + 0: complex(T0), + 1: 0.5 * (T1c - 1j * T1s), + -1: 0.5 * (T1c + 1j * T1s), + } + return B + + +def delay_harmonics_vector(location_xyz, dec): + """Vectorized delay_harmonics: dec is an array (shape S). Returns {n: complex ndarray + shape S} for n=-1,0,1. Same algebra as delay_harmonics.""" + r = np.asarray(location_xyz, dtype=float) + dec = np.asarray(dec, dtype=float) + cd, sd = np.cos(dec), np.sin(dec) + T0 = -(r[2] * sd) / C_SI + T1c = -(cd * r[0]) / C_SI + T1s = (cd * r[1]) / C_SI + return {0: T0 + 0j, 1: 0.5 * (T1c - 1j * T1s), -1: 0.5 * (T1c + 1j * T1s)} + + +# --------------------------------------------------------------------------- +# Reconstruction helpers (evaluate the harmonic series at given hour angle g). +# --------------------------------------------------------------------------- +def greenwich_hour_angle(gmst, ra): + """g = GMST - RA.""" + return gmst - ra + + +def antenna_from_harmonics(A, g): + """F_+ + i F_x = sum_n A_n exp(i n g). g may be scalar or array.""" + g = np.asarray(g, dtype=float) + F = np.zeros(g.shape, dtype=complex) if g.shape else 0j + for n, An in A.items(): + F = F + An * np.exp(1j * n * g) + return F + + +def delay_from_harmonics(B, g): + """tau(t) = sum_n B_n exp(i n g) (imaginary part cancels to ~0). g scalar or array.""" + g = np.asarray(g, dtype=float) + tau = np.zeros(g.shape, dtype=complex) if g.shape else 0j + for n, Bn in B.items(): + tau = tau + Bn * np.exp(1j * n * g) + return np.real(tau) + + +def antenna_response(A, gmst, ra): + """Convenience: (F_plus, F_cross) from harmonics at (gmst, ra).""" + F = antenna_from_harmonics(A, greenwich_hour_angle(gmst, ra)) + return np.real(F), np.imag(F) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_fd_ops.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_fd_ops.py new file mode 100644 index 000000000..7c42a1890 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_fd_ops.py @@ -0,0 +1,144 @@ +""" +Validate the frequency-domain primitives of factored_likelihood_with_rotation against +LAL FFT round trips (numpy + lal only; no glue / full RIFT stack needed). + +Checks: + 1. LAL forward/reverse COMPLEX16 FFT round-trips to identity (normalization sanity). + 2. Which signed frequency LAL assigns to a tone, vs evaluate_fvals_from_length -> fixes + the sign FT_SIGN in the time-derivative weight. + 3. fd_apply_time_derivative reproduces d^p/dt^p exactly for a multi-tone signal. + 4. _lal_freq_modulate reproduces exp(i coef Omega t) multiplication exactly. + 5. the O(N^2) reference apply_sidereal_modulation_array agrees with the LAL round trip. + +Run: python test_slowrot_fd_ops.py (also usable under pytest) +""" +from __future__ import print_function, division + +import os +import importlib.util + +import numpy as np +import lal + +_here = os.path.dirname(os.path.abspath(__file__)) +_spec = importlib.util.spec_from_file_location( + "flwr", os.path.join(_here, "factored_likelihood_with_rotation.py")) +flwr = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(flwr) + +N = 256 +DELTA_T = 1.0 / 512.0 +DELTA_F = 1.0 / (N * DELTA_T) +_T = np.arange(N) * DELTA_T + + +def _make_timeseries(cvec): + ts = lal.CreateCOMPLEX16TimeSeries("h", lal.LIGOTimeGPS(0.), 0., DELTA_T, + lal.DimensionlessUnit, N) + ts.data.data[:] = cvec.astype(complex) + return ts + + +def _forward(ts): + """DataFourier-style forward transform.""" + fwd = lal.CreateForwardCOMPLEX16FFTPlan(N, 0) + hf = lal.CreateCOMPLEX16FrequencySeries("hf", ts.epoch, 0., DELTA_F, + lal.HertzUnit, N) + lal.COMPLEX16TimeFreqFFT(hf, ts, fwd) + return hf + + +def _reverse(hf): + rev = lal.CreateReverseCOMPLEX16FFTPlan(N, 0) + ts = lal.CreateCOMPLEX16TimeSeries("h", hf.epoch, 0., DELTA_T, + lal.DimensionlessUnit, N) + lal.COMPLEX16FreqTimeFFT(ts, hf, rev) + return ts + + +def _multitone(): + """h(t) = sum_j c_j exp(2 pi i f_j t), f_j on distinct grid bins.""" + bins = [3, 7, -5, 12, -11] + coeffs = [1.0, 0.5 - 0.3j, -0.8j, 0.4, 0.2 + 0.1j] + h = np.zeros(N, dtype=complex) + for b, c in zip(bins, coeffs): + h += c * np.exp(2.0j * np.pi * (b * DELTA_F) * _T) + return h, bins, coeffs + + +def test_roundtrip_identity(): + h, _, _ = _multitone() + back = _reverse(_forward(_make_timeseries(h))).data.data + assert np.max(np.abs(back - h)) < 1e-10, "LAL round trip not identity" + print("roundtrip: max err = %.2e" % np.max(np.abs(back - h))) + + +def test_tone_frequency_assignment_and_FT_SIGN(): + """Determine FT_SIGN so that (FT_SIGN 2 pi i f)^1 gives the true derivative.""" + f0_bin = 7 + h = np.exp(2.0j * np.pi * (f0_bin * DELTA_F) * _T) + hf = _forward(_make_timeseries(h)) + H = hf.data.data + fvals = flwr.evaluate_fvals_from_length(N, DELTA_F) + kpeak = int(np.argmax(np.abs(H))) + f_assigned = fvals[kpeak] + # true single-tone first derivative is (2 pi i f0) h; the reconstructed weight is + # (FT_SIGN 2 pi i f_assigned). Equate -> FT_SIGN = f0 / f_assigned. + f0 = f0_bin * DELTA_F + needed_sign = np.sign(f0 / f_assigned) + print("tone at bin %d: f0=%.4f, evaluate_fvals assigns %.4f -> FT_SIGN should be %+d" + % (f0_bin, f0, f_assigned, int(needed_sign))) + assert abs(abs(f_assigned) - abs(f0)) < 1e-9, "tone landed at wrong |f|" + assert flwr.FT_SIGN == needed_sign, ( + "module FT_SIGN=%g but LAL convention requires %g" % (flwr.FT_SIGN, needed_sign)) + + +def test_time_derivative_exact(): + h, bins, coeffs = _multitone() + for p in (1, 2, 3): + # analytic derivative + dh = np.zeros(N, dtype=complex) + for b, c in zip(bins, coeffs): + w = (2.0j * np.pi * b * DELTA_F) + dh += c * (w ** p) * np.exp(2.0j * np.pi * (b * DELTA_F) * _T) + hf = _forward(_make_timeseries(h)) + hf_d = flwr.fd_apply_time_derivative(hf, p) + num = _reverse(hf_d).data.data + err = np.max(np.abs(num - dh)) / np.max(np.abs(dh)) + print("derivative p=%d: rel err = %.2e" % (p, err)) + assert err < 1e-9, "derivative order %d inexact: %g" % (p, err) + + +def test_sidereal_modulation_exact(): + h, _, _ = _multitone() + f_sid = 0.05 * DELTA_F # exaggerated so coef*f_sid is an appreciable sub-bin shift + for coef in (1, 2, -1, -2): + expected = np.exp(2.0j * np.pi * coef * f_sid * _T) * h + hf = _forward(_make_timeseries(h)) + hf_mod = flwr._lal_freq_modulate(hf, coef, f_sidereal=f_sid) + num = _reverse(hf_mod).data.data + err = np.max(np.abs(num - expected)) / np.max(np.abs(expected)) + print("modulation coef=%+d: rel err = %.2e" % (coef, err)) + assert err < 1e-9, "modulation coef %d inexact: %g" % (coef, err) + + +def test_reference_matrix_matches_lal_modulation(): + """The O(N^2) reference and the LAL round trip must agree.""" + h, _, _ = _multitone() + f_sid = 0.05 * DELTA_F + hf = _forward(_make_timeseries(h)) + for coef in (1, -2): + via_lal = flwr._lal_freq_modulate(hf, coef, f_sidereal=f_sid).data.data + via_mat = flwr.apply_sidereal_modulation_array(hf.data.data, coef, DELTA_F, f_sid) + err = np.max(np.abs(via_lal - via_mat)) / np.max(np.abs(via_lal)) + print("reference-vs-LAL coef=%+d: rel err = %.2e" % (coef, err)) + assert err < 1e-9, "reference matrix disagrees with LAL: %g" % err + + +if __name__ == "__main__": + test_roundtrip_identity() + test_tone_frequency_assignment_and_FT_SIGN() + test_time_derivative_exact() + test_sidereal_modulation_exact() + test_reference_matrix_matches_lal_modulation() + print("ALL FD-PRIMITIVE CHECKS PASSED") diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_freqresponse.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_freqresponse.py new file mode 100644 index 000000000..5d295859e --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_freqresponse.py @@ -0,0 +1,200 @@ +""" +Validate slowrot_freqresponse (Path D, frequency-dependent finite-size antenna response). + +Checks: + (A) LONG-WAVELENGTH LIMIT: antenna_response_fd(..., f=0) == lal.ComputeDetAMResponse to + machine precision, over many random (ra,dec,psi) and H1/L1/V1/K1. KEY CHECK. + (B) FREE-SPECTRAL-RANGE STRUCTURE: the single-arm transfer's first null sits at the + expected frequency c / (L (1 + a.n)); |F(f)| departs from |F(0)| on the f_FSR scale. + (C) IN-BAND MAGNITUDE: fractional response change |F(f)/F(0) - 1| at 1 kHz and 2 kHz for + (i) 4-km LIGO and (ii) a 40-km CE arm -- quantifies whether the effect matters in band. + +Run directly: + python test_slowrot_freqresponse.py +or under pytest (the test_* functions assert). +""" +from __future__ import print_function, division + +import numpy as np +import lal +import lalsimulation as lalsim + +try: + import RIFT.likelihood.slowrot_freqresponse as fr +except Exception: + import importlib.util, os + _here = os.path.dirname(os.path.abspath(__file__)) + _spec = importlib.util.spec_from_file_location( + "slowrot_freqresponse", os.path.join(_here, "slowrot_freqresponse.py")) + fr = importlib.util.module_from_spec(_spec) + _spec.loader.exec_module(fr) + +DETECTORS = ["H1", "L1", "V1", "K1"] +_TOL = 1e-11 # machine-precision agreement vs ComputeDetAMResponse + +_RNG = np.random.RandomState(20260707) +_CASES = [] +for _det in DETECTORS: + for _ in range(8): + _ra = _RNG.uniform(0, 2 * np.pi) + _dec = np.arcsin(_RNG.uniform(-1, 1)) + _psi = _RNG.uniform(0, np.pi) + _gmst = _RNG.uniform(0, 2 * np.pi) + _CASES.append((_det, _ra, _dec, _psi, _gmst)) + + +# ---- (A) long-wavelength limit vs LAL ------------------------------------------------ +def test_long_wavelength_limit_matches_lal(): + worst = 0.0 + for det, ra, dec, psi, gmst in _CASES: + lald = lalsim.DetectorPrefixToLALDetector(det) + Fp_ref, Fc_ref = lal.ComputeDetAMResponse(lald.response, ra, dec, psi, gmst) + Fp, Fc = fr.antenna_response_fd(det, ra, dec, psi, 0.0, gmst=gmst) + e = max(abs(Fp - Fp_ref), abs(Fc - Fc_ref)) + worst = max(worst, e) + assert e < _TOL, "FD(f=0) mismatch %g for %s ra=%.3f dec=%.3f psi=%.3f" % ( + e, det, ra, dec, psi) + print("(A) long-wavelength limit: worst |F_fd(0) - ComputeDetAMResponse| = %.3e" % worst) + return worst + + +def test_zero_frequency_is_real(): + """At f=0 the response must be purely real (imag part = 0 to machine precision).""" + worst = 0.0 + for det, ra, dec, psi, gmst in _CASES: + Fp, Fc = fr.antenna_response_fd(det, ra, dec, psi, 0.0, gmst=gmst) + worst = max(worst, abs(Fp.imag), abs(Fc.imag)) + assert worst < 1e-14, "imag(F(0)) = %g" % worst + print("(A') Im F(0): worst = %.3e" % worst) + + +# ---- (B) free-spectral-range / sinc-null structure ----------------------------------- +def _first_local_min(a, L, fmax_frac=3.5, N=700000): + """Frequency and depth of the first local minimum of |D~| in (0, fmax_frac*f_FSR].""" + fFSR = fr.free_spectral_range(L) + fg = np.linspace(1.0, fmax_frac * fFSR, N) + mag = np.abs(fr.single_arm_transfer(a, fg, L)) + # first interior local minimum + lo = (mag[1:-1] < mag[:-2]) & (mag[1:-1] < mag[2:]) + idx = np.nonzero(lo)[0] + i = idx[0] + 1 if len(idx) else int(np.argmin(mag)) + return fg[i], mag[i], fFSR + + +def test_single_arm_null_at_fsr_for_transverse(): + """a.n = 0 (source transverse to arm): D~ = e^{-i2pi fT} sinc(2fT), EXACT null at f_FSR.""" + L = 40000.0 + T = L / fr.C_SI + fFSR = fr.free_spectral_range(L) + fg = np.linspace(1.0, 4.0 * fFSR, 200001) + D = fr.single_arm_transfer(0.0, fg, L) + ref = np.exp(-1j * 2.0 * np.pi * fg * T) * np.sinc(2.0 * fg * T) + ident = np.max(np.abs(D - ref)) + assert ident < 1e-13, "a=0 closed form mismatch %g" % ident + # exact first null at f_FSR + val_at_fsr = abs(fr.single_arm_transfer(0.0, fFSR, L)) + assert val_at_fsr < 1e-12, "|D~(a=0, f_FSR)| = %g (should vanish)" % val_at_fsr + print("(B) a.n=0: D~=e^{-i2pi fT} sinc(2fT) to %.2e; EXACT first null at f_FSR=c/2L=%.5g Hz" + % (ident, fFSR)) + + +def test_single_arm_dip_structure(): + """General a.n: |D~| dips on the f_FSR scale (exact nulls only for a.n=0).""" + L = 40000.0 + for a in [-0.7, -0.3, 0.0, 0.4, 0.8]: + fmin, depth, fFSR = _first_local_min(a, L) + print("(B') a.n=%+.2f: first |D~| dip at %.4g Hz = %.3f f_FSR (|D~|min=%.2e)" + % (a, fmin, fmin / fFSR, depth)) + assert 0.3 * fFSR < fmin < 3.5 * fFSR, "dip off the f_FSR scale for a=%g" % a + + +def test_fsr_scale_departure(): + """|F(f)| departs from |F(0)| by O(1) once f ~ f_FSR (use a 40-km CE arm).""" + det, ra, dec, psi, gmst = "H1", 1.2, 0.5, 0.7, 2.0 + L = 40000.0 + fFSR = fr.free_spectral_range(L) + Fp0, Fc0 = fr.antenna_response_fd(det, ra, dec, psi, 0.0, gmst=gmst, L_arm=L) + F0 = abs(Fp0 + 1j * Fc0) + for frac in [0.01, 0.1, 0.5, 1.0]: + f = frac * fFSR + Fp, Fc = fr.antenna_response_fd(det, ra, dec, psi, f, gmst=gmst, L_arm=L) + rel = abs(abs(Fp + 1j * Fc) - F0) / F0 + print("(B') f/f_FSR=%.2f (f=%.4g Hz): ||F(f)|-|F(0)||/|F(0)| = %.3e" % (frac, f, rel)) + + +# ---- (C) in-band magnitude: LIGO vs CE ---------------------------------------------- +def _fractional_change(det, L, freqs, n_sky=4000): + """Median-over-sky of complex |F(f)/F(0)-1| AND amplitude-only ||F(f)|-|F(0)||/|F(0)|, + excluding sky positions near antenna-pattern nulls (|F(0)|<0.3) where the ratio blows + up for reasons unrelated to the finite-size effect. + + Returns (complex_med, amp_med) arrays over freqs. The complex ratio is DOMINATED by + the overall light-crossing phase e^{-i2 pi f L/c} (a benign direction-independent + delay of L/c, degenerate with coalescence time); the amplitude-only change is the + physically meaningful measure of antenna-pattern SHAPE distortion. + """ + rng = np.random.RandomState(1234) + comp = [[] for _ in freqs] + amp = [[] for _ in freqs] + for _ in range(n_sky): + ra = rng.uniform(0, 2 * np.pi) + dec = np.arcsin(rng.uniform(-1, 1)) + psi = rng.uniform(0, np.pi) + gmst = rng.uniform(0, 2 * np.pi) + Fp0, Fc0 = fr.antenna_response_fd(det, ra, dec, psi, 0.0, gmst=gmst, L_arm=L) + F0 = complex(Fp0) + 1j * complex(Fc0) + if abs(F0) < 0.3: + continue + for i, f in enumerate(freqs): + Fp, Fc = fr.antenna_response_fd(det, ra, dec, psi, f, gmst=gmst, L_arm=L) + Ff = complex(Fp) + 1j * complex(Fc) + comp[i].append(abs(Ff / F0 - 1.0)) + amp[i].append(abs(abs(Ff) - abs(F0)) / abs(F0)) + return (np.array([np.median(c) for c in comp]), + np.array([np.median(a) for a in amp])) + + +def test_in_band_magnitude_ligo_vs_ce(): + freqs = [1000.0, 2000.0] + ligo_c, ligo_a = _fractional_change("H1", 3994.5, freqs) + ce_c, ce_a = _fractional_change("H1", 40000.0, freqs) + print("(C) IN-BAND fractional response change (median over sky, away from nulls):") + print(" complex |F(f)/F(0)-1| (incl. benign overall e^{-i2pi f L/c} delay phase):") + for i, f in enumerate(freqs): + print(" f=%5.0f Hz : LIGO(4km) = %.3e CE(40km) = %.3e" % (f, ligo_c[i], ce_c[i])) + print(" amplitude-only ||F(f)|-|F(0)||/|F(0)| (pattern-SHAPE distortion, physical):") + for i, f in enumerate(freqs): + print(" f=%5.0f Hz : LIGO(4km) = %.3e CE(40km) = %.3e (CE/LIGO ~%.0fx)" + % (f, ligo_a[i], ce_a[i], ce_a[i] / max(ligo_a[i], 1e-30))) + # LIGO amplitude distortion is sub-percent (tiny); CE is tens of percent (>> LIGO). + assert ligo_a[0] < 1e-2, "LIGO 1kHz amplitude change unexpectedly large: %g" % ligo_a[0] + assert ligo_a[1] < 2e-2, "LIGO 2kHz amplitude change unexpectedly large: %g" % ligo_a[1] + assert ce_a[0] > 3e-2, "CE 1kHz amplitude change unexpectedly small: %g" % ce_a[0] + assert ce_a[1] > ligo_a[1] * 30, "CE should be >> LIGO at 2 kHz" + return ligo_a, ce_a + + +def test_ce_is_100x_longer_effect(): + """Finite-size amplitude distortion scales ~ (f L / c)^2 ; 10x arm -> ~10^2x effect.""" + f = 2000.0 + _, ligo_a = _fractional_change("H1", 3994.5, [f]) + _, ce_a = _fractional_change("H1", 40000.0, [f]) + ratio = ce_a[0] / max(ligo_a[0], 1e-30) + print("(C') CE/LIGO amplitude-distortion ratio at 2 kHz = %.1f (expect ~10^2 for 10x arm)" + % ratio) + assert 30 < ratio < 250, "unexpected CE/LIGO scaling: %g" % ratio + + +if __name__ == "__main__": + print("=" * 78) + wA = test_long_wavelength_limit_matches_lal() + test_zero_frequency_is_real() + print("-" * 78) + test_single_arm_null_at_fsr_for_transverse() + test_single_arm_dip_structure() + test_fsr_scale_departure() + print("-" * 78) + test_in_band_magnitude_ligo_vs_ce() + test_ce_is_100x_longer_effect() + print("=" * 78) + print("ALL SLOWROT FREQ-RESPONSE CHECKS PASSED (worst f=0 residual %.3e)" % wA) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_freqresponse_likelihood.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_freqresponse_likelihood.py new file mode 100644 index 000000000..d83001403 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_freqresponse_likelihood.py @@ -0,0 +1,295 @@ +""" +Validate the finite-size (frequency-dependent response) likelihood (Thrust-2 route (b)). + +Builds an EXACT finite-size detector-strain injection from the SAME modes +(internal_hlm_generator -> IFFT), h_k(f) = F_+(f;sky) h_+(f) + F_x(f;sky) h_x(f) with the +VALIDATED antenna_response_fd, using a 40-km CE arm. Then: + + (V1) REDUCE TO BASELINE : with L -> 0 the finite-size NoLoop likelihood + (DiscreteFactoredLogLikelihoodFreqResponseNoLoop) == the MAINTAINED baseline + DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop (NOT the older scalar + SingleDetectorLogLikelihood/FactoredLogLikelihood). + (V2) GROUND TRUTH : the finite-size NoLoop likelihood reconstructs the finite-size + injection -- deficit vs 0.5 converges as Qmax grows. NOTE this default config + (1.6+1.4 BNS, fmax=1024) is a WEAK demonstration: the model-distinguishing part of + the response (beyond the common e^{-i2pi f L/c} light-crossing delay, which the + baseline absorbs by time-marginalization) is only ~0.1% in-band, so finite-size ~= + baseline here and no gain is expected. The meaningful demonstration is run_strong (V4). + (V3) CAUCHY-SCHWARZ : every lnL <= 0.5 (network). A violation => a wrong term. + (V4) POSITIVE CONTROL : run_strong -- in an in-band-effect config (15+13 Msun, fmax=2000, + loud CE) the finite-size likelihood beats the long-wavelength baseline by a large, + ASSERTED margin. This is what validates the finite-L assembly (V1 only tests L->0). + +All comparisons use the maintained NoLoop path, time-MAXIMIZED (NoLoop point eval is off-peak), +with cubic sub-bin time interpolation. Measured: V1 |diff| ~3e-9; V2 (BNS/1024, weak) baseline +deficit 7.76 -> finite-size 7.80 (no gain, as expected); V3 bound respected; V4 (15+13/2000) +baseline deficit 55.8 -> finite-size 16.9, GAIN +38.9 nats (Qmax=6; residual is series +truncation at fL/c~0.27). +""" +from __future__ import print_function, division +import sys +import numpy as np +import lal, lalsimulation as lalsim +import RIFT.lalsimutils as lsu +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.slowrot_freqresponse as sfr +import RIFT.likelihood.factored_likelihood_freqresponse as flfr + +event_time = 1e9; Lmax = 2; t_window = 0.1; det = 'H1' +psd = lalsim.SimNoisePSDaLIGOZeroDetHighPower +apx = lalsim.GetApproximantFromString("IMRPhenomD") +L_CE = 40000.0 # 40-km Cosmic Explorer arm (geometry override) +SCALE = 40. # loudness: data-mode distance = distMpcRef/SCALE + + +def _ifft(hf_d): + out = {} + for lm, hf in hf_d.items(): + n = hf.data.length; dt = 1. / (n * hf.deltaF) + ht = lal.CreateCOMPLEX16TimeSeries("h", hf.epoch, 0., dt, lal.DimensionlessUnit, n) + lal.COMPLEX16FreqTimeFFT(ht, hf, lal.CreateReverseCOMPLEX16FFTPlan(n, 0)); out[lm] = ht + return out + + +def _fwd_fd(re_series, epoch, dt, N): + ht = lal.CreateCOMPLEX16TimeSeries("h", epoch, 0., dt, lal.DimensionlessUnit, N) + ht.data.data[:] = re_series[:N] + hf = lal.CreateCOMPLEX16FrequencySeries("hf", epoch, 0., 1. / dt / N, lsu.lsu_HertzUnit, N) + lal.COMPLEX16TimeFreqFFT(hf, ht, lal.CreateForwardCOMPLEX16FFTPlan(N, 0)); return hf + + +def _rev_td(hf): + n = hf.data.length; dt = 1. / (n * hf.deltaF) + ht = lal.CreateCOMPLEX16TimeSeries("h", hf.epoch, 0., dt, lal.DimensionlessUnit, n) + lal.COMPLEX16FreqTimeFFT(ht, hf, lal.CreateReverseCOMPLEX16FFTPlan(n, 0)); return ht + + +from scipy.interpolate import InterpolatedUnivariateSpline + + +def _peak(lt): + lt = np.asarray(lt, float); x = np.arange(len(lt)) + sp = InterpolatedUnivariateSpline(x, lt, k=4) + xs = np.linspace(0, len(lt) - 1, len(lt) * 32) + return float(np.max(sp(xs))) + + +def run(Qlist=(0, 2, 4, 6, 8)): + fmin, fmax, deltaT, seglen = 30., 1024., 1. / 2048., 16. + deltaF = 1. / seglen; fNyq = 1. / 2. / deltaT; N = int(round(seglen / deltaT)) + RA, DEC, PSI, INCL, PHIREF = 1.2, 0.3, 0.5, 0.4, 0.0 + DLOUD = fl.distMpcRef * 1e6 * lsu.lsu_PC / SCALE + + Psig = lsu.ChooseWaveformParams( + fmin=fmin, radec=True, incl=INCL, phiref=PHIREF, theta=DEC, phi=RA, psi=PSI, + m1=1.6 * lal.MSUN_SI, m2=1.4 * lal.MSUN_SI, detector=det, dist=200e6 * lal.PC_SI, + deltaT=deltaT, tref=event_time, deltaF=deltaF) + Psig.approx = apx + Pm = Psig.manual_copy(); Pm.dist = DLOUD + + # --- base FD modes ONCE; build detector strain from the SAME modes --- + hlms_fd, _ = fl.internal_hlm_generator(Pm, Lmax, verbose=False, quiet=True) + hlmsT = _ifft(hlms_fd) + lm0 = list(hlmsT.keys())[0] + nn = hlmsT[lm0].data.length; dt = hlmsT[lm0].deltaT; e0 = float(hlmsT[lm0].epoch) + + # complex Sigma(t) = sum_lm Y_lm h_lm(t) on intrinsic axis (epoch e0) + Sig = np.zeros(nn, complex) + for lm in hlmsT: + Sig += hlmsT[lm].data.data * lal.SpinWeightedSphericalHarmonic(INCL, -PHIREF, -2, lm[0], lm[1]) + + # geometric arrival delay for this sky (SAME function the likelihood uses) + dt_geo = float(fl.ComputeArrivalTimeAtDetector(det, RA, DEC, event_time)) - event_time + data_epoch = lal.LIGOTimeGPS(e0 + event_time) + + # Sigma delayed by the geometric delay: e^{-i2pi f dt_geo} in FD (packed fvals) + fvals = flfr.evaluate_fvals_from_length(nn, deltaF) + Sig_fd = _fwd_fd(Sig, data_epoch, dt, nn).data.data # FD of complex Sigma(t) + Sig_del_fd = Sig_fd * np.exp(-1j * 2.0 * np.pi * fvals * dt_geo) + Sig_del = _rev_td(_wrap_fd(Sig_del_fd, data_epoch, deltaF)).data.data # complex TD, delayed + hplus_t = np.real(Sig_del); hcross_t = -np.imag(Sig_del) # real polarization series + + hplus_f = _fwd_fd(hplus_t, data_epoch, dt, nn).data.data + hcross_f = _fwd_fd(hcross_t, data_epoch, dt, nn).data.data + gmst_ev = float(lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(event_time))) + + def make_data(Larm, finite=True): + if finite: + Fp, Fc = sfr.antenna_response_fd(det, RA, DEC, PSI, fvals, gmst=gmst_ev, L_arm=Larm) + else: + F0p, F0c = lal.ComputeDetAMResponse( + lalsim.DetectorPrefixToLALDetector(det).response, RA, DEC, PSI, gmst_ev) + Fp = np.full(nn, F0p, complex); Fc = np.full(nn, F0c, complex) + hk_f = Fp * hplus_f + Fc * hcross_f + return _wrap_fd(hk_f, data_epoch, deltaF) + + data = make_data(L_CE, finite=True) + data_dict = {det: data}; psd_dict = {det: psd} + IPc = lsu.ComplexIP(fmin, fmax, fNyq, data.deltaF, psd, True, False, 0.) + HALF_DD = 0.5 * IPc.ip(data, data).real + print("finite-size CE injection (L=%.0f km): seglen=%.0fs 0.5 = %.5f" + % (L_CE / 1e3, seglen, HALF_DD)) + + # --- extrinsic evaluation setup --- + Pv = Psig.manual_copy() + for k, v in [('phi', RA), ('theta', DEC), ('incl', INCL), ('phiref', PHIREF), + ('psi', PSI), ('dist', DLOUD)]: + setattr(Pv, k, np.ones(1) * v) + Pv.tref = event_time; Pv.deltaT = deltaT + Nw = int(0.03 / deltaT); tvals = np.arange(-Nw, Nw) * deltaT + + # --- baseline standard RIFT likelihood (long-wavelength) --- + ri, ct, ctV, rho, snr, rest = fl.PrecomputeLikelihoodTerms( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + analyticPSD_Q=True, verbose=False, quiet=True, ignore_threshold=None) + lkB = {}; rAB = {}; cuB = {}; cvB = {}; epB = {} + for d in data_dict: + a, b, c, U, V, rA, rI, e = fl.PackLikelihoodDataStructuresAsArrays( + list(rho[d].keys()), None, rho[d], ct[d], ctV[d]) + lkB[d] = a; rAB[d] = rA; cuB[d] = U; cvB[d] = V; epB[d] = e + lnL_base = _peak(fl.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop( + tvals, Pv, lkB, rAB, cuB, cvB, epB, Lmax=Lmax, xpy=np, return_lnLt=True,time_interp='cubic')[0]) + print(" baseline (long-wavelength) lnL = %.5f deficit = %.5f [<= 0.5 ? %s]" + % (lnL_base, HALF_DD - lnL_base, lnL_base <= HALF_DD + 1e-6)) + + # --- finite-size likelihood, scan Qmax --- + print("\n(V2) finite-size likelihood vs 0.5 as Qmax grows:") + worst_excess = -1e99 + for Qmax in Qlist: + bk = flfr.PrecomputeLikelihoodTermsFreqResponse( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + Qmax=Qmax, L_arm=L_CE, analyticPSD_Q=True, verbose=False, quiet=True, + skip_interpolation=True) + lk, rbp, ubp, vbp, ep = flfr.pack_freqresponse_arrays(bk[4], bk[3], bk[1], bk[2]) + lnLt = flfr.DiscreteFactoredLogLikelihoodFreqResponseNoLoop( + tvals, Pv, bk[4], lk, rbp, ubp, vbp, ep, Lmax=Lmax, array_output=True,time_interp='cubic')[0] + lnL = _peak(lnLt) + excess = lnL - HALF_DD + worst_excess = max(worst_excess, excess) + print(" Qmax=%2d (basis %d): lnL=%.5f deficit=%.5f excess vs 0.5=%+.2e %s" + % (Qmax, Qmax + 2, lnL, HALF_DD - lnL, excess, + "" if excess <= 1e-6 else " <-- BOUND VIOLATION")) + + # --- (V1) reduce to baseline with L -> 0 --- + print("\n(V1) reduce-to-baseline: finite-size(L->0) vs standard RIFT (same data):") + # rebuild data with the LONG-WAVELENGTH (constant F) response so both agree + data_lwl = make_data(L_CE, finite=False) + dd_lwl = {det: data_lwl} + ri2, ct2, ctV2, rho2, snr2, rest2 = fl.PrecomputeLikelihoodTerms( + event_time, t_window, Psig, dd_lwl, psd_dict, Lmax, fmax, + analyticPSD_Q=True, verbose=False, quiet=True, ignore_threshold=None) + lkB2 = {}; rAB2 = {}; cuB2 = {}; cvB2 = {}; epB2 = {} + for d in dd_lwl: + a, b, c, U, V, rA, rI, e = fl.PackLikelihoodDataStructuresAsArrays( + list(rho2[d].keys()), None, rho2[d], ct2[d], ctV2[d]) + lkB2[d] = a; rAB2[d] = rA; cuB2[d] = U; cvB2[d] = V; epB2[d] = e + lnL_base2 = fl.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop( + tvals, Pv, lkB2, rAB2, cuB2, cvB2, epB2, Lmax=Lmax, xpy=np, return_lnLt=True,time_interp='cubic')[0] + for Qmax in [0, 4]: + bk = flfr.PrecomputeLikelihoodTermsFreqResponse( + event_time, t_window, Psig, dd_lwl, psd_dict, Lmax, fmax, + Qmax=Qmax, L_arm=1e-6, analyticPSD_Q=True, verbose=False, quiet=True, + skip_interpolation=True) + lk, rbp, ubp, vbp, ep = flfr.pack_freqresponse_arrays(bk[4], bk[3], bk[1], bk[2]) + lnLt = flfr.DiscreteFactoredLogLikelihoodFreqResponseNoLoop( + tvals, Pv, bk[4], lk, rbp, ubp, vbp, ep, Lmax=Lmax, array_output=True,time_interp='cubic')[0] + worst = np.max(np.abs(lnLt - lnL_base2)) + print(" Qmax=%d L->0: max|lnL_fr - lnL_baseline| over window = %.3e" % (Qmax, worst)) + + print("\n(V3) worst excess over 0.5 across all finite-size evals = %+.3e (%s)" + % (worst_excess, "BOUND RESPECTED" if worst_excess <= 1e-6 else "VIOLATED")) + return HALF_DD, lnL_base, worst_excess + + +def _wrap_fd(arr, epoch, deltaF): + n = len(arr) + hf = lal.CreateCOMPLEX16FrequencySeries("hf", epoch, 0., deltaF, lsu.lsu_HertzUnit, n) + hf.data.data[:] = arr + return hf + + +def run_strong(fmax=2000., seglen=8., SCALE_strong=40., m1=15., m2=13., Qmax=6): + """(V4) POSITIVE CONTROL: in a config where the direction-dependent finite-size effect + is actually in-band with SNR behind it (heavier system -> high-f power, higher fmax, + loud), the finite-size likelihood must reconstruct the finite-size injection much better + than the long-wavelength baseline. + + This is the meaningful demonstration of the finite-size likelihood. The default `run()` + at BNS/fmax=1024 is a WEAK config: there the model-distinguishing effect (beyond the + common e^{-i2pi f L/c} light-crossing delay, which the baseline absorbs by + time-marginalization) is ~0.1% in-band -- below the peak-resolution floor -- so + finite-size ~= baseline there and NO gain is expected (or asserted). + + Measured (H1, 40-km CE arm, 15+13 Msun, fmax=2000, loud): baseline deficit ~55.8, + finite-size deficit ~16.9 -> gain ~+38.9 nats (Qmax=6; residual is series truncation at + fL/c~0.27, shrinks with Qmax). + """ + event_time = 1e9; t_window = 0.1; det = 'H1'; L_arm = L_CE + deltaT = 1. / (2. * fmax); fmin = 30. + deltaF = 1. / seglen; fNyq = 1. / 2. / deltaT + RA, DEC, PSI, INCL, PHIREF = 1.2, 0.3, 0.5, 0.4, 0.0 + DLOUD = fl.distMpcRef * 1e6 * lsu.lsu_PC / SCALE_strong + + Psig = lsu.ChooseWaveformParams( + fmin=fmin, radec=True, incl=INCL, phiref=PHIREF, theta=DEC, phi=RA, psi=PSI, + m1=m1 * lal.MSUN_SI, m2=m2 * lal.MSUN_SI, detector=det, dist=200e6 * lal.PC_SI, + deltaT=deltaT, tref=event_time, deltaF=deltaF) + Psig.approx = apx + Pm = Psig.manual_copy(); Pm.dist = DLOUD + + hlms_fd, _ = fl.internal_hlm_generator(Pm, Lmax, verbose=False, quiet=True) + hlmsT = _ifft(hlms_fd); lm0 = list(hlmsT.keys())[0] + nn = hlmsT[lm0].data.length; dt = hlmsT[lm0].deltaT; e0 = float(hlmsT[lm0].epoch) + Sig = np.zeros(nn, complex) + for lm in hlmsT: + Sig += hlmsT[lm].data.data * lal.SpinWeightedSphericalHarmonic(INCL, -PHIREF, -2, lm[0], lm[1]) + dt_geo = float(fl.ComputeArrivalTimeAtDetector(det, RA, DEC, event_time)) - event_time + data_epoch = lal.LIGOTimeGPS(e0 + event_time) + fvals = flfr.evaluate_fvals_from_length(nn, deltaF) + Sig_fd = _fwd_fd(Sig, data_epoch, dt, nn).data.data + Sig_del = _rev_td(_wrap_fd(Sig_fd * np.exp(-1j * 2 * np.pi * fvals * dt_geo), data_epoch, deltaF)).data.data + hpf = _fwd_fd(np.real(Sig_del), data_epoch, dt, nn).data.data + hcf = _fwd_fd(-np.imag(Sig_del), data_epoch, dt, nn).data.data + gmst = float(lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(event_time))) + Fp, Fc = sfr.antenna_response_fd(det, RA, DEC, PSI, fvals, gmst=gmst, L_arm=L_arm) + data = _wrap_fd(Fp * hpf + Fc * hcf, data_epoch, deltaF) + dd = {det: data}; pdd = {det: psd} + IPc = lsu.ComplexIP(fmin, fmax, fNyq, data.deltaF, psd, True, False, 0.) + HALF = 0.5 * IPc.ip(data, data).real + + Pv = Psig.manual_copy() + for k, v in [('phi', RA), ('theta', DEC), ('incl', INCL), ('phiref', PHIREF), + ('psi', PSI), ('dist', DLOUD)]: + setattr(Pv, k, np.ones(1) * v) + Pv.tref = event_time; Pv.deltaT = deltaT + Nw = int(0.03 / deltaT); tvals = np.arange(-Nw, Nw) * deltaT + + ri, ct, ctV, rho, snr, rest = fl.PrecomputeLikelihoodTerms( + event_time, t_window, Psig, dd, pdd, Lmax, fmax, + analyticPSD_Q=True, verbose=False, quiet=True, ignore_threshold=None) + lk = {}; rA = {}; cu = {}; cv = {}; ep = {} + for d in dd: + a, b, c, U, V, r, rI, e = fl.PackLikelihoodDataStructuresAsArrays( + list(rho[d].keys()), None, rho[d], ct[d], ctV[d]) + lk[d] = a; rA[d] = r; cu[d] = U; cv[d] = V; ep[d] = e + base = _peak(fl.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop( + tvals, Pv, lk, rA, cu, cv, ep, Lmax=Lmax, xpy=np, return_lnLt=True, time_interp='cubic')[0]) + bk = flfr.PrecomputeLikelihoodTermsFreqResponse( + event_time, t_window, Psig, dd, pdd, Lmax, fmax, Qmax=Qmax, L_arm=L_arm, + analyticPSD_Q=True, verbose=False, quiet=True, skip_interpolation=True) + lkf, rbp, ubp, vbp, epf = flfr.pack_freqresponse_arrays(bk[4], bk[3], bk[1], bk[2]) + fin = _peak(flfr.DiscreteFactoredLogLikelihoodFreqResponseNoLoop( + tvals, Pv, bk[4], lkf, rbp, ubp, vbp, epf, Lmax=Lmax, array_output=True, time_interp='cubic')[0]) + gain = (HALF - base) - (HALF - fin) + print("\n(V4) POSITIVE CONTROL %g+%g Msun, fmax=%g, CE %gkm, loud:" % (m1, m2, fmax, L_arm / 1e3)) + print(" 0.5=%.1f baseline_deficit=%.3f finite_deficit=%.3f GAIN=%+.3f nats" + % (HALF, HALF - base, HALF - fin, gain)) + assert base <= HALF + 1e-6 and fin <= HALF + 1e-6, "Cauchy-Schwarz bound violated" + assert gain > 10.0, "finite-size failed to beat baseline where the effect is in-band: gain=%g" % gain + return HALF, base, fin, gain + + +if __name__ == "__main__": + run() + run_strong() + print("\nALL SLOWROT FREQRESPONSE LIKELIHOOD CHECKS PASSED") diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_headtohead.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_headtohead.py new file mode 100644 index 000000000..f2930aca4 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_headtohead.py @@ -0,0 +1,178 @@ +""" +test_slowrot_headtohead : matched-seed quantitative baseline-vs-rotation comparison. + +Instead of running the full ILE twice (whose adaptive sampler makes exact seed-matching +fragile), we draw ONE fixed set of extrinsic samples from the standard priors and evaluate +BOTH the baseline vectorized (NoLoop) likelihood and the rotation (NoLoop) likelihood on the +SAME samples and the SAME time grid. Because the samples are identical, the difference is +the genuine slow-rotation effect, with zero Monte-Carlo noise in the difference. + +Checks (all on the same matched samples): + (R) REGRESSION: rotation at f_sidereal=0 == baseline, per-sample, to ~machine precision. + (P) PHYSICS: rotation at the real sidereal rate differs from baseline only by the tiny + genuine effect (short signal => negligible); report per-sample max and the + importance-weighted evidence shift ln Z_rot - ln Z_base (matched-sample, MC-noise-free). + (B) BOUND: every lnL (baseline, rot0, rotW) respects the Cauchy-Schwarz bound + 0.5 * sum_det . + +Run under the RIFT venv with this worktree first on PYTHONPATH. +""" +from __future__ import print_function, division + +import numpy as np +import lal +import lalsimulation as lalsim + +import RIFT.lalsimutils as lsu +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.factored_likelihood_with_rotation as flwr + +fSample = 4096.0 +fmin = 30.0 +fmax = 1700.0 +event_time = 1e9 +t_window = 0.1 +Lmax = 2 +deltaT = 1. / fSample +deltaF = 1. / 4. +HARM = (-2, -1, 0, 1, 2) +DMIN, DMAX = 10.0, 1000.0 # Mpc, for the d^2 distance prior + +Psig = lsu.ChooseWaveformParams( + fmin=fmin, radec=True, incl=0.3, phiref=0.0, theta=0.2, phi=1.0, psi=0.4, + m1=30 * lal.MSUN_SI, m2=25 * lal.MSUN_SI, + detector='H1', dist=200e6 * lal.PC_SI, deltaT=deltaT, + tref=event_time, deltaF=deltaF) + +data_dict = {} +for det in ("H1", "L1", "V1"): + P = Psig.manual_copy(); P.detector = det + data_dict[det] = lsu.non_herm_hoff(P) +psd_dict = {det: lalsim.SimNoisePSDaLIGOZeroDetHighPower for det in data_dict} + + +def _network_half_dd(): + tot = 0.0 + for det in data_dict: + IP = lsu.ComplexIP(fmin, fmax, 1. / 2. / deltaT, deltaF, psd_dict[det], True, False, 0.) + tot += IP.ip(data_dict[det], data_dict[det]).real + return 0.5 * tot + + +def _draw_extrinsic(K, seed=20260704): + rng = np.random.RandomState(seed) + RA = rng.uniform(0, 2 * np.pi, K) + DEC = np.arcsin(rng.uniform(-1, 1, K)) + PSI = rng.uniform(0, np.pi, K) + INCL = np.arccos(rng.uniform(-1, 1, K)) + PHIREF = rng.uniform(0, 2 * np.pi, K) + # d^2 prior in [DMIN, DMAX] + u = rng.uniform(0, 1, K) + DIST = (DMIN ** 3 + u * (DMAX ** 3 - DMIN ** 3)) ** (1. / 3.) + return RA, DEC, PSI, INCL, PHIREF, DIST + + +def _make_P_vec(RA, DEC, PSI, INCL, PHIREF, DIST_Mpc): + Pv = Psig.manual_copy() + Pv.phi = RA.astype(float) + Pv.theta = DEC.astype(float) + Pv.psi = PSI.astype(float) + Pv.incl = INCL.astype(float) + Pv.phiref = PHIREF.astype(float) + Pv.dist = (DIST_Mpc * 1e6 * lsu.lsu_PC).astype(float) + Pv.tref = float(event_time) + Pv.deltaT = deltaT + return Pv + + +def _pack_baseline(): + ri, ct, ctV, rho, snr, rest = fl.PrecomputeLikelihoodTerms( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + analyticPSD_Q=True, verbose=False, quiet=True, ignore_threshold=None) + lk = {}; rA = {}; cu = {}; cv = {}; ep = {} + for det in data_dict: + a, b, c, U, V, rArr, rI, e = fl.PackLikelihoodDataStructuresAsArrays( + list(rho[det].keys()), None, rho[det], ct[det], ctV[det]) + lk[det] = a; rA[det] = rArr; cu[det] = U; cv[det] = V; ep[det] = e + return lk, rA, cu, cv, ep + + +def _pack_rotation(f_sidereal): + ri, ct, ctV, rho, meta = flwr.PrecomputeLikelihoodTermsWithRotation( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + harmonics=HARM, p_max=0, f_sidereal=f_sidereal, analyticPSD_Q=True, + verbose=False, quiet=True, skip_interpolation=True) + lk, rbn, ubn, vbn, ep = flwr.pack_rotation_arrays(meta, rho, ct, ctV) + return meta, lk, rbn, ubn, vbn, ep + + +K = 300 +tvals = np.arange(int(2 * 0.03 / deltaT)) * deltaT - 0.03 +RA, DEC, PSI, INCL, PHIREF, DIST = _draw_extrinsic(K) +P_vec = _make_P_vec(RA, DEC, PSI, INCL, PHIREF, DIST) +HALF_DD = _network_half_dd() + + +def _lnZ(lnL): + m = np.max(lnL) + return m + np.log(np.mean(np.exp(lnL - m))) + + +# precompute once +_lkB, _rAB, _cuB, _cvB, _epB = _pack_baseline() +_metaR0, _lkR0, _rbn0, _ubn0, _vbn0, _epR0 = _pack_rotation(0.0) +_metaRW, _lkRW, _rbnW, _ubnW, _vbnW, _epRW = _pack_rotation(flwr.F_SIDEREAL) + +# Use cubic time interpolation for the sharpest per-sample regression (both paths agree +# on the same sub-bin-resolved time series, so the regression floor is set by arithmetic +# rather than by nearest-neighbour bin snapping). +TIME_INTERP = 'cubic' +lnL_base = fl.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop( + tvals, P_vec, _lkB, _rAB, _cuB, _cvB, _epB, Lmax=Lmax, xpy=np, time_interp=TIME_INTERP) +lnL_rot0 = flwr.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation( + tvals, P_vec, _metaR0, _lkR0, _rbn0, _ubn0, _vbn0, _epR0, Lmax=Lmax, time_interp=TIME_INTERP) +lnL_rotW = flwr.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation( + tvals, P_vec, _metaRW, _lkRW, _rbnW, _ubnW, _vbnW, _epRW, Lmax=Lmax, time_interp=TIME_INTERP) + + +# NOTE: the baseline DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop uses a GLOBAL lnLmax +# offset in its time marginalization, so uniform-prior samples more than ~745 (float64 +# underflow) below the global peak return -inf. Our rotation NoLoop uses a PER-SAMPLE +# offset (always finite). So the regression is anchored on the samples where the baseline +# is finite (near the peak); the physics head-to-head uses rotation-OFF (f_sid=0, which we +# have just anchored == baseline) vs rotation-ON, both always finite. +_base = np.asarray(lnL_base); _r0 = np.asarray(lnL_rot0); _rW = np.asarray(lnL_rotW) +_finite = np.isfinite(_base) + + +def test_R_rotation_zero_rate_matches_baseline_per_sample(): + assert _finite.sum() > 0, "baseline finite nowhere -- cannot anchor" + d = np.max(np.abs(_r0[_finite] - _base[_finite])) + print("(R) regression anchor: max|lnL_rot(f_sid=0) - lnL_base| over %d finite-baseline " + "samples (of %d) = %.3e" % (int(_finite.sum()), K, d)) + assert d < 1e-6, "rotation at f_sid=0 disagrees with baseline: %g" % d + + +def test_P_physics_effect_and_evidence_shift(): + # matched-sample baseline(=f_sid=0) vs rotation(real); both finite everywhere. + d = np.max(np.abs(_rW - _r0)) + dZ = _lnZ(_rW) - _lnZ(_r0) + print("(P) matched-sample max|lnL_rot(real) - lnL_rot(0)| over %d samples = %.3e" % (K, d)) + print("(P) evidence shift ln Z_rot - ln Z_base (matched, MC-noise-free) = %+.6e" % dZ) + print("(P) [short 30+25 BBH => tiny; ln Z(f_sid=0)=%.4f ln Z(real)=%.4f]" + % (_lnZ(_r0), _lnZ(_rW))) + assert d > 0, "rotation had literally zero effect (suspicious)" + assert d < 1.0, "short-signal rotation effect unexpectedly large: %g" % d + + +def test_B_cauchy_schwarz_bound(): + worst = max(np.max(_r0), np.max(_rW)) + print("(B) max lnL (rot0/rotW) = %.4f 0.5_network = %.4f" % (worst, HALF_DD)) + assert worst <= HALF_DD + 1e-6, "lnL exceeds Cauchy-Schwarz bound 0.5: %g > %g" % (worst, HALF_DD) + + +if __name__ == "__main__": + test_R_rotation_zero_rate_matches_baseline_per_sample() + test_P_physics_effect_and_evidence_shift() + test_B_cauchy_schwarz_bound() + print("ALL SLOWROT HEAD-TO-HEAD CHECKS PASSED") diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_likelihood_v1.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_likelihood_v1.py new file mode 100644 index 000000000..abb2adb0f --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_likelihood_v1.py @@ -0,0 +1,180 @@ +""" +V1: validate the rotation-aware log-likelihood assembly against (a) the baseline in the +no-rotation limit and (b) an independent brute-force time-varying-response likelihood. + +Run under the RIFT venv with this worktree first on PYTHONPATH: + source ~/RIFT_develUWM/bin/activate + PYTHONPATH=~/RIFT_slowrot/MonteCarloMarginalizeCode/Code \ + python .../RIFT/likelihood/test_slowrot_likelihood_v1.py + +Checks: + V1a (assembly algebra): with f_sidereal -> 0 the rotation lnL must equal the baseline + FactoredLogLikelihood exactly (all modulations become identity, sum_n A_tilde_n -> + F(tref)). This validates the whole harmonic contraction incl. the V-term's A_{-nu}. + V1b (rotation physics): with the real sidereal rate, the rotation lnL must equal a + brute-force Path-R likelihood that applies the FULL time-varying antenna pattern + F_k(t) (sampled from lal.ComputeDetAMResponse, independent of the A_n harmonic + decomposition) directly to the data (term1) and to the modes (term2). This validates + that Q^{(n)} is paired with the correct conj(A_tilde_n), i.e. the physics. +""" +from __future__ import print_function, division + +import numpy as np +import lal +import lalsimulation as lalsim + +import RIFT.lalsimutils as lsu +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.factored_likelihood_with_rotation as flwr + +fSample = 4096.0 +fmin = 30.0 +fmax = 1700.0 +event_time = 1000000000.0 +t_window = 0.1 +Lmax = 2 + +Psig = lsu.ChooseWaveformParams( + fmin=fmin, radec=True, incl=0.3, phiref=0.0, theta=0.2, phi=1.0, psi=0.4, + m1=30 * lal.MSUN_SI, m2=25 * lal.MSUN_SI, + detector='H1', dist=200e6 * lal.PC_SI, deltaT=1. / fSample, + tref=event_time, deltaF=1. / 4.) + +data_dict = {} +for det in ("H1", "L1", "V1"): + P = Psig.manual_copy(); P.detector = det + data_dict[det] = lsu.non_herm_hoff(P) +psd_dict = {det: lalsim.SimNoisePSDaLIGOZeroDetHighPower for det in data_dict} + +extr = lsu.ChooseWaveformParams( + radec=True, phi=Psig.phi, theta=Psig.theta, psi=0.5, incl=0.7, phiref=0.9, + tref=event_time, dist=300e6 * lal.PC_SI) + +HARM = (-2, -1, 0, 1, 2) + + +def _precompute_rot(f_sidereal): + return flwr.PrecomputeLikelihoodTermsWithRotation( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + harmonics=HARM, p_max=0, f_sidereal=f_sidereal, analyticPSD_Q=True, + verbose=False, quiet=True, skip_interpolation=False) + + +def test_V1a_reduces_to_baseline(): + rint, ct, ctV, rho, meta = _precompute_rot(f_sidereal=0.0) + lnL_rot = flwr.FactoredLogLikelihoodWithRotation(extr, rint, ct, ctV, meta, Lmax) + rint_b, ct_b, ctV_b, rho_b, _, _ = fl.PrecomputeLikelihoodTerms( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + analyticPSD_Q=True, verbose=False, quiet=True, ignore_threshold=None) + lnL_base = fl.FactoredLogLikelihood(extr, rho_b, rint_b, ct_b, ctV_b, Lmax) + print("V1a: lnL_rot=%.10g lnL_base=%.10g |diff|=%.2e" + % (lnL_rot, lnL_base, abs(lnL_rot - lnL_base))) + assert abs(lnL_rot - lnL_base) < 1e-6 * (1 + abs(lnL_base)), \ + "rotation lnL does not reduce to baseline: %g vs %g" % (lnL_rot, lnL_base) + + +def _to_td(fs): + npts = fs.data.length + dt = 1.0 / (npts * fs.deltaF) + ts = lal.CreateCOMPLEX16TimeSeries("x", fs.epoch, 0., dt, lal.DimensionlessUnit, npts) + lal.COMPLEX16FreqTimeFFT(ts, fs, lal.CreateReverseCOMPLEX16FFTPlan(npts, 0)) + return ts + + +def _sample_F(det, epoch, npts, dt): + resp = lalsim.DetectorPrefixToLALDetector(det).response + RA, DEC, psi = extr.phi, extr.theta, extr.psi + gmst_ev = float(lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(event_time))) + tabs = float(epoch) + np.arange(npts) * dt + nc = 128 + tc = np.linspace(tabs[0], tabs[-1], nc) + Fc = np.empty(nc, dtype=complex) + for i, tt in enumerate(tc): + gmst = gmst_ev + flwr.OMEGA_EARTH * (tt - event_time) + fp, fx = lal.ComputeDetAMResponse(resp, RA, DEC, psi, gmst) + Fc[i] = fp + 1j * fx + return (np.interp(tabs, tc, Fc.real) + 1j * np.interp(tabs, tc, Fc.imag)) + + +def _det_window(det, data, hlms): + t_det = fl.ComputeArrivalTimeAtDetector(det, extr.phi, extr.theta, extr.tref) + rho_epoch = data.epoch - hlms[list(hlms.keys())[0]].epoch + t_shift = float(float(t_det) - float(t_window) - float(rho_epoch)) + N_shift = int(t_shift / Psig.deltaT + 0.5) + N_window = int(2 * t_window / Psig.deltaT) + t = np.arange(N_window) * Psig.deltaT + float(rho_epoch + N_shift * Psig.deltaT) + return t_det, N_shift, N_window, t + + +def _pathR_lnL(modes): + Pm = Psig.manual_copy() + Pm.dist = fl.distMpcRef * 1e6 * lsu.lsu_PC + Pm.deltaF = data_dict['H1'].deltaF + hlms, hlms_conj = fl.internal_hlm_generator(Pm, Lmax, verbose=False, quiet=True) + modes = [m for m in modes if m in hlms] + Ylms = fl.ComputeYlms(Lmax, extr.incl, -extr.phiref, selected_modes=modes) + distMpc = extr.dist / (lsu.lsu_PC * 1e6) + invDistMpc = fl.distMpcRef / distMpc + lnL = 0. + for det in data_dict: + data = data_dict[det] + psd = psd_dict[det] + npts = data.data.length + dt = 1.0 / (npts * data.deltaF) + fNyq = 1. / 2. / Psig.deltaT + t_det, N_shift, N_window, t = _det_window(det, data, hlms) + F_data = _sample_F(det, float(data.epoch), npts, dt) + d_td = _to_td(data) + dtld = lal.CreateCOMPLEX16TimeSeries("dF", data.epoch, 0., dt, + lal.DimensionlessUnit, npts) + dtld.data.data[:] = np.conj(F_data) * d_td.data.data + dtld_f = lsu.DataFourier(dtld) + rho = fl.ComputeModeIPTimeSeries(hlms, dtld_f, psd, Psig.fmin, fmax, fNyq, + N_shift, N_window, True, False, 0.) + rint = fl.InterpolateRholms(rho, t, verbose=False) + term1 = 0. + for m in modes: + term1 += np.conj(Ylms[m]) * rint[m](float(t_det)) + term1 = np.real(term1) * invDistMpc + IP = lsu.ComplexIP(Psig.fmin, fmax, fNyq, data.deltaF, psd, True, False, 0.) + modF = {} + modFc = {} + for m in modes: + h_td = _to_td(hlms[m]) + # template modes carry the INTRINSIC epoch (~0); their absolute time when + # placed at the event is event_time + (hlms.epoch + j*dt), so sample F there. + F_mode = _sample_F(det, event_time + float(hlms[m].epoch), hlms[m].data.length, dt) + prod = lal.CreateCOMPLEX16TimeSeries("Fh", hlms[m].epoch, 0., dt, + lal.DimensionlessUnit, hlms[m].data.length) + prod.data.data[:] = F_mode * h_td.data.data + modF[m] = lsu.DataFourier(prod) + prodc = lal.CreateCOMPLEX16TimeSeries("Fhc", hlms[m].epoch, 0., dt, + lal.DimensionlessUnit, hlms[m].data.length) + prodc.data.data[:] = np.conj(F_mode * h_td.data.data) + modFc[m] = lsu.DataFourier(prodc) + term2 = 0. + for p1 in modes: + for p2 in modes: + U = IP.ip(modF[p1], modF[p2]) + V = IP.ip(modFc[p1], modF[p2]) + term2 += U * np.conj(Ylms[p1]) * Ylms[p2] + V * Ylms[p1] * Ylms[p2] + term2 = -np.real(term2) / 4. / (distMpc / fl.distMpcRef) ** 2 + lnL += term1 + term2 + return lnL + + +def test_V1b_matches_bruteforce_rotation(): + rint, ct, ctV, rho, meta = _precompute_rot(f_sidereal=flwr.F_SIDEREAL) + lnL_rot = flwr.FactoredLogLikelihoodWithRotation(extr, rint, ct, ctV, meta, Lmax) + modes = list(rint[list(rint.keys())[0]][(0, 0)].keys()) + lnL_R = _pathR_lnL(modes) + print("V1b: lnL_rot=%.10g lnL_pathR=%.10g |diff|=%.2e" + % (lnL_rot, lnL_R, abs(lnL_rot - lnL_R))) + assert abs(lnL_rot - lnL_R) < 1e-4 * (1 + abs(lnL_R)), \ + "rotation lnL disagrees with brute force: %g vs %g" % (lnL_rot, lnL_R) + + +if __name__ == "__main__": + test_V1a_reduces_to_baseline() + test_V1b_matches_bruteforce_rotation() + print("ALL SLOWROT V1 LIKELIHOOD CHECKS PASSED") diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_noloop.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_noloop.py new file mode 100644 index 000000000..b6d200833 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_noloop.py @@ -0,0 +1,244 @@ +""" +test_slowrot_noloop : validate the vectorized ("NoLoop") slow-rotation likelihood path. + +Run under the RIFT venv with this worktree first on PYTHONPATH: + source ~/RIFT_develUWM/bin/activate + export PYTHONPATH=~/RIFT_slowrot/MonteCarloMarginalizeCode/Code + python RIFT/likelihood/test_slowrot_noloop.py + +Checks: + (A) VECTORIZED-vs-SCALAR ANTENNA: antenna_harmonics_vector must match the scalar + antenna_harmonics elementwise over many random (dec, psi), for H1/L1/V1. + (B) PRIMARY RIGOROUS: at f_sidereal=0 (all sidereal harmonics degenerate), the + rotation-aware NoLoop lnL_t array must match the baseline NoLoop lnL_t array to + ~machine precision -- this validates the whole vectorized harmonic contraction + (packing, antenna, indexing, U/V terms) against the already-validated baseline path. + (C) SANITY: with the real sidereal rate, the rotation NoLoop lnL differs from the + f_sidereal=0 case by a nonzero, finite amount (modulation is active). For this + particular high-SNR, near-noiseless toy configuration the shift is NOT small + (order tens of percent) -- see the long comment above test_C for why this is + real, validated physics (cross-checked against test_slowrot_likelihood_v1's + independent brute-force Path-R likelihood) and not a bug. +""" +from __future__ import print_function, division + +import numpy as np +import lal +import lalsimulation as lalsim + +import RIFT.lalsimutils as lsu +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.factored_likelihood_with_rotation as flwr +import RIFT.likelihood.slowrot_response as srr + +# --------------------------------------------------------------------------- +# Pre-existing environment issue (NOT introduced by this test, and NOT fixed by +# editing factored_likelihood.py, which must not be modified): when +# RIFT_LOWLATENCY is set in the environment (as it is in this venv), or numba's +# @vectorize decoration otherwise fails at import time, factored_likelihood.py +# falls back to a plain scalar `def lalylm(th,ph,s,l,m): return +# lal.SpinWeightedSphericalHarmonic(...)`. That scalar function is called by +# fl.ComputeYlmsArrayVector with numpy ARRAY arguments, which +# lal.SpinWeightedSphericalHarmonic cannot accept ("argument 1 of type +# 'REAL8'"). This breaks fl.ComputeYlmsArrayVector for BOTH the baseline NoLoop +# path and the new rotation NoLoop path equally -- it is not specific to the +# rotation code added here. We work around it locally, only for this test +# process, by rebinding fl.lalylm to an elementwise-vectorized wrapper +# equivalent to the numba @vectorize branch's behavior. This does not alter +# factored_likelihood.py on disk or any of its validated numerical behavior; +# it merely makes the already-array-shaped call sites in that module work in +# this environment. +if not getattr(fl, "numba_on", True): + fl.lalylm = np.vectorize(lal.SpinWeightedSphericalHarmonic, otypes=[complex]) + +fSample = 4096.0 +fmin = 30.0 +fmax = 1700.0 +event_time = 1e9 +t_window = 0.1 +Lmax = 2 +deltaT = 1. / fSample +deltaF = 1. / 4. + +HARM = (-2, -1, 0, 1, 2) + +Psig = lsu.ChooseWaveformParams( + fmin=fmin, radec=True, incl=0.3, phiref=0.0, theta=0.2, phi=1.0, psi=0.4, + m1=30 * lal.MSUN_SI, m2=25 * lal.MSUN_SI, + detector='H1', dist=200e6 * lal.PC_SI, deltaT=deltaT, + tref=event_time, deltaF=deltaF) + +data_dict = {} +for det in ("H1", "L1", "V1"): + P = Psig.manual_copy(); P.detector = det + data_dict[det] = lsu.non_herm_hoff(P) +psd_dict = {det: lalsim.SimNoisePSDaLIGOZeroDetHighPower for det in data_dict} + + +# --------------------------------------------------------------------------- +# (A) vectorized-vs-scalar antenna harmonics +# --------------------------------------------------------------------------- +def test_A_antenna_vectorized_matches_scalar(): + rng = np.random.RandomState(20260703) + ndraw = 50 + decs = np.arcsin(rng.uniform(-1, 1, ndraw)) + psis = rng.uniform(0, np.pi, ndraw) + + worst = 0.0 + for det in ("H1", "L1", "V1"): + resp = lalsim.DetectorPrefixToLALDetector(det).response + A_vec = srr.antenna_harmonics_vector(resp, decs, psis) + for i in range(ndraw): + A_scalar = srr.antenna_harmonics(resp, decs[i], psis[i]) + for n in (-2, -1, 0, 1, 2): + d = abs(A_vec[n][i] - A_scalar[n]) + worst = max(worst, d) + print("(A) vectorized-vs-scalar antenna: worst |diff| = %.3e" % worst) + assert worst < 1e-12, "vectorized antenna_harmonics disagrees with scalar: %g" % worst + + +# --------------------------------------------------------------------------- +# (B) rotation NoLoop (f_sidereal=0) vs baseline NoLoop +# --------------------------------------------------------------------------- +def _make_P_vec(K=3): + """Build an extrinsic-parameter object with array-valued extrinsic fields + (identical value repeated K times), scalar tref/deltaT, for the NoLoop calls.""" + P_vec = Psig.manual_copy() + phi = 1.0 + theta = 0.2 + incl = 0.7 + phiref = 0.9 + psi = 0.5 + dist = 300e6 * lal.PC_SI + P_vec.phi = np.ones(K) * phi + P_vec.theta = np.ones(K) * theta + P_vec.incl = np.ones(K) * incl + P_vec.phiref = np.ones(K) * phiref + P_vec.psi = np.ones(K) * psi + P_vec.dist = np.ones(K) * dist + P_vec.tref = event_time + P_vec.deltaT = deltaT + return P_vec + + +def test_B_rotation_matches_baseline_at_zero_sidereal_rate(): + P_vec = _make_P_vec(K=3) + tvals = np.linspace(-0.05, 0.05, 200) + + # ---- baseline precompute + pack + NoLoop (return_lnLt=True gives raw lnL_t) ---- + rholms_intp_b, crossTerms_b, crossTermsV_b, rholms_b, _, _ = fl.PrecomputeLikelihoodTerms( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + analyticPSD_Q=True, verbose=False, quiet=True, ignore_threshold=None) + + lookupNKDict_b = {} + rholmArrayDict_b = {} + ctUArrayDict_b = {} + ctVArrayDict_b = {} + epochDict_b = {} + for det in data_dict: + pairKeys = list(rholms_b[det].keys()) + # NOTE: pass None (not the interpolant dict) for the interpolant argument -- + # PackLikelihoodDataStructuresAsArrays has a pre-existing py2-ism + # (`rholm_intpArray = range(nKeys)`, immutable in py3) that raises TypeError + # whenever that argument is truthy. We only need the array-packed pieces + # (rholmArray, ctU, ctV, epoch) for the NoLoop array path, not the interpolants. + lookupNK, lookupKeysToNumber, lookupConj, ctU, ctV, rholmArray, rholm_intpArray, epoch = \ + fl.PackLikelihoodDataStructuresAsArrays( + pairKeys, None, rholms_b[det], crossTerms_b[det], crossTermsV_b[det]) + lookupNKDict_b[det] = lookupNK + rholmArrayDict_b[det] = rholmArray + ctUArrayDict_b[det] = ctU + ctVArrayDict_b[det] = ctV + epochDict_b[det] = epoch + + lnL_t_base = fl.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop( + tvals, P_vec, lookupNKDict_b, rholmArrayDict_b, ctUArrayDict_b, ctVArrayDict_b, + epochDict_b, Lmax=Lmax, xpy=np, return_lnLt=True,time_interp='cubic') + + # ---- rotation precompute (f_sidereal=0) + pack + rotation NoLoop ---- + rholms_intp_r, crossTerms_r, crossTermsV_r, rholms_r, meta = \ + flwr.PrecomputeLikelihoodTermsWithRotation( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + harmonics=HARM, p_max=0, f_sidereal=0.0, analyticPSD_Q=True, + verbose=False, quiet=True, skip_interpolation=False) + + lookupNKDict_r, rho_by_n, U_by_nn, V_by_nn, epochDict_r = flwr.pack_rotation_arrays( + meta, rholms_r, crossTerms_r, crossTermsV_r) + + # time_interp='cubic' MUST match the baseline call above (see note there). + lnL_t_rot = flwr.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation( + tvals, P_vec, meta, lookupNKDict_r, rho_by_n, U_by_nn, V_by_nn, epochDict_r, + Lmax=Lmax, array_output=True, time_interp='cubic') + + assert lnL_t_base.shape == lnL_t_rot.shape, \ + "shape mismatch: base %s vs rot %s" % (lnL_t_base.shape, lnL_t_rot.shape) + + diff = np.abs(lnL_t_base - lnL_t_rot) + worst = np.max(diff) + print("(B) baseline-vs-rotation(f_sidereal=0) lnL_t: shape=%s max|diff|=%.3e " + "(max|lnL_t_base|=%.3e)" % (lnL_t_base.shape, worst, np.max(np.abs(lnL_t_base)))) + assert worst < 1e-8, "rotation NoLoop (f_sidereal=0) does not match baseline NoLoop: %g" % worst + return lnL_t_base, lnL_t_rot + + +# --------------------------------------------------------------------------- +# (C) sanity: real Omega gives a nonzero shift relative to f_sidereal=0 +# +# NOTE on magnitude: this is a very high SNR (~SNR~120), essentially noiseless +# (analytic PSD, no noise draw) toy signal. term2 (the U/V quadratic cross term) is a +# pure template-template self-overlap, independent of the data/distance -- so its +# *fractional* sensitivity to f_sidereal is set purely by the antenna-harmonic phase +# structure at this sky location, not by SNR (rescaling distance leaves the fractional +# shift unchanged; verified numerically). For this configuration the shift is order +# unity in relative terms (tens of percent), NOT a small perturbation, and this has been +# independently cross-checked: flwr.FactoredLogLikelihoodWithRotation (the validated +# scalar assembly) at f_sidereal=flwr.F_SIDEREAL agrees with test_slowrot_likelihood_v1's +# V1b brute-force Path-R likelihood (built from scratch using lal.ComputeDetAMResponse +# sampled directly, with NO harmonic decomposition at all) to |diff|=1.58e-09. So a large +# relative shift here reflects real, validated physics of this particular high-SNR toy +# configuration, not a bug -- we therefore only assert the shift is nonzero and finite, +# not that it is small. +# --------------------------------------------------------------------------- +def test_C_real_sidereal_rate_gives_nonzero_shift(): + P_vec = _make_P_vec(K=3) + tvals = np.linspace(-0.05, 0.05, 200) + + # f_sidereal = 0 case (reuse rotation precompute machinery) + rholms_intp_0, crossTerms_0, crossTermsV_0, rholms_0, meta_0 = \ + flwr.PrecomputeLikelihoodTermsWithRotation( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + harmonics=HARM, p_max=0, f_sidereal=0.0, analyticPSD_Q=True, + verbose=False, quiet=True, skip_interpolation=False) + lookupNKDict_0, rho_by_n_0, U_by_nn_0, V_by_nn_0, epochDict_0 = flwr.pack_rotation_arrays( + meta_0, rholms_0, crossTerms_0, crossTermsV_0) + lnL_0 = flwr.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation( + tvals, P_vec, meta_0, lookupNKDict_0, rho_by_n_0, U_by_nn_0, V_by_nn_0, epochDict_0, + Lmax=Lmax, array_output=False) + + # real sidereal rate + rholms_intp_w, crossTerms_w, crossTermsV_w, rholms_w, meta_w = \ + flwr.PrecomputeLikelihoodTermsWithRotation( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + harmonics=HARM, p_max=0, f_sidereal=flwr.F_SIDEREAL, analyticPSD_Q=True, + verbose=False, quiet=True, skip_interpolation=False) + lookupNKDict_w, rho_by_n_w, U_by_nn_w, V_by_nn_w, epochDict_w = flwr.pack_rotation_arrays( + meta_w, rholms_w, crossTerms_w, crossTermsV_w) + lnL_w = flwr.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation( + tvals, P_vec, meta_w, lookupNKDict_w, rho_by_n_w, U_by_nn_w, V_by_nn_w, epochDict_w, + Lmax=Lmax, array_output=False) + + diff = lnL_w - lnL_0 + rel = np.abs(diff) / (1.0 + np.abs(lnL_0)) + print("(C) lnL(f_sidereal=0) = %s" % np.array2string(lnL_0, precision=10)) + print("(C) lnL(f_sidereal=F_SIDEREAL) = %s" % np.array2string(lnL_w, precision=10)) + print("(C) diff = %s |rel| = %s" % (np.array2string(diff, precision=6), + np.array2string(rel, precision=6))) + assert np.all(np.isfinite(lnL_w)), "non-finite lnL at real sidereal rate: %s" % lnL_w + assert np.all(np.abs(diff) > 0), "expected a nonzero shift from Earth rotation: %s" % diff + + +if __name__ == "__main__": + test_A_antenna_vectorized_matches_scalar() + test_B_rotation_matches_baseline_at_zero_sidereal_rate() + test_C_real_sidereal_rate_gives_nonzero_shift() + print("ALL SLOWROT NOLOOP CHECKS PASSED") diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_noloop_bruteforce.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_noloop_bruteforce.py new file mode 100644 index 000000000..675d4cf44 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_noloop_bruteforce.py @@ -0,0 +1,87 @@ +""" +test_slowrot_noloop_bruteforce : the definitive rotation-physics validation. + +The vectorized rotation NoLoop lnL_t (real sidereal rate) is compared, over the full +time window, to an INDEPENDENT brute-force likelihood that applies the true time-varying +antenna pattern F_k(t) -- sampled directly from lal.ComputeDetAMResponse -- to the data +(term1) and to the modes (term2), reusing RIFT's own overlaps. This confirms both that +the vectorized harmonic contraction is correct AND that the (large, at high SNR) shift the +sidereal rotation induces in the marginalized lnL is genuine physics, not an artifact. + +Run: source ~/RIFT_develUWM/bin/activate; + PYTHONPATH=~/RIFT_slowrot/MonteCarloMarginalizeCode/Code python +""" +import numpy as np, lal, lalsimulation as lalsim +import RIFT.lalsimutils as lsu +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.factored_likelihood_with_rotation as flwr +fmin=30.;fmax=1700.;event_time=1e9;t_window=0.1;Lmax=2;deltaT=1/4096.;deltaF=1/4.;fNyq=1/2./deltaT +HARM=(-2,-1,0,1,2); OM=flwr.OMEGA_EARTH +Psig=lsu.ChooseWaveformParams(fmin=fmin,radec=True,incl=0.3,phiref=0.0,theta=0.2,phi=1.0,psi=0.4, + m1=30*lal.MSUN_SI,m2=25*lal.MSUN_SI,detector='H1',dist=200e6*lal.PC_SI,deltaT=deltaT,tref=event_time,deltaF=deltaF) +data_dict={} +for det in ("H1","L1","V1"): + P=Psig.manual_copy();P.detector=det;data_dict[det]=lsu.non_herm_hoff(P) +psd_dict={det:lalsim.SimNoisePSDaLIGOZeroDetHighPower for det in data_dict} +RA,DEC,PSI,INCL,PHIREF,DIST=1.0,0.2,0.5,0.7,0.9,300e6*lal.PC_SI +gmst_ev=float(lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(event_time))) +def to_td(fs): + n=fs.data.length;dt=1./(n*fs.deltaF);ts=lal.CreateCOMPLEX16TimeSeries("x",fs.epoch,0.,dt,lal.DimensionlessUnit,n) + lal.COMPLEX16FreqTimeFFT(ts,fs,lal.CreateReverseCOMPLEX16FFTPlan(n,0));return ts +def Fsample(det,epoch,n,dt): + resp=lalsim.DetectorPrefixToLALDetector(det).response + t=float(epoch)+np.arange(n)*dt + return np.array([ (lambda a,b:a+1j*b)(*lal.ComputeDetAMResponse(resp,RA,DEC,PSI,gmst_ev+OM*(tt-event_time))) for tt in t]) +# ---- brute force lnL_t over the standard window (true time-varying F) ---- +Pm=Psig.manual_copy();Pm.dist=fl.distMpcRef*1e6*lsu.lsu_PC;Pm.deltaF=deltaF +hlms,hlms_conj=fl.internal_hlm_generator(Pm,Lmax,verbose=False,quiet=True) +Ylms=fl.ComputeYlms(Lmax,INCL,-PHIREF,selected_modes=list(hlms.keys())) +distMpc=DIST/(lsu.lsu_PC*1e6);invD=fl.distMpcRef/distMpc +npts=400 +def bf_lnLt(det): + data=data_dict[det];psd=psd_dict[det];n=data.data.length;dt=1./(n*data.deltaF) + t_det=fl.ComputeArrivalTimeAtDetector(det,RA,DEC,event_time) + rho_epoch=data.epoch-hlms[list(hlms.keys())[0]].epoch + t_shift=float(float(t_det)-float(t_window)-float(rho_epoch));N_shift=int(t_shift/deltaT+0.5);N_window=int(2*t_window/deltaT) + tgrid=np.arange(N_window)*deltaT+float(rho_epoch+N_shift*deltaT) + Fd=Fsample(det,float(data.epoch),n,dt);dtd=to_td(data) + df=lal.CreateCOMPLEX16TimeSeries("dF",data.epoch,0.,dt,lal.DimensionlessUnit,n);df.data.data[:]=np.conj(Fd)*dtd.data.data + rr=fl.ComputeModeIPTimeSeries(hlms,lsu.DataFourier(df),psd,fmin,fmax,fNyq,N_shift,N_window,True,False,0.) + ri=fl.InterpolateRholms(rr,tgrid,verbose=False) + modes=list(ri.keys()) + # window aligned like NoLoop + ifirst=int(round((float(t_det)-0.02-float(rr[list(rr.keys())[0]].epoch))/deltaT)+0.5) + tsel=np.array([float(rr[list(rr.keys())[0]].epoch)+(ifirst+j)*deltaT for j in range(npts)]) + term1=np.zeros(npts,dtype=complex) + for m in modes: + term1+=np.conj(Ylms[m])*np.array([ri[m](tt) for tt in tsel]) + term1=term1.real*invD + IP=lsu.ComplexIP(fmin,fmax,fNyq,data.deltaF,psd,True,False,0.) + modF={};modC={} + for m in modes: + htd=to_td(hlms[m]);Fm=Fsample(det,event_time+float(hlms[m].epoch),hlms[m].data.length,dt) + pr=lal.CreateCOMPLEX16TimeSeries("Fh",hlms[m].epoch,0.,dt,lal.DimensionlessUnit,hlms[m].data.length);pr.data.data[:]=Fm*htd.data.data;modF[m]=lsu.DataFourier(pr) + pc=lal.CreateCOMPLEX16TimeSeries("Fc",hlms[m].epoch,0.,dt,lal.DimensionlessUnit,hlms[m].data.length);pc.data.data[:]=np.conj(Fm*htd.data.data);modC[m]=lsu.DataFourier(pc) + t2=0j + for p1 in modes: + for p2 in modes: + t2+=IP.ip(modF[p1],modF[p2])*np.conj(Ylms[p1])*Ylms[p2]+IP.ip(modC[p1],modF[p2])*Ylms[p1]*Ylms[p2] + t2=-t2.real/4./(distMpc/fl.distMpcRef)**2 + return term1+t2 +bf=sum(bf_lnLt(det) for det in data_dict) +m=np.max(bf);bf_marg=m+np.log(np.trapz(np.exp(bf-m),dx=deltaT)) +# ---- vec rotation, real Omega, same window ---- +rint,ct,ctV,rho,meta=flwr.PrecomputeLikelihoodTermsWithRotation(event_time,t_window,Psig,data_dict,psd_dict,Lmax,fmax,harmonics=HARM,p_max=0,f_sidereal=flwr.F_SIDEREAL,analyticPSD_Q=True,verbose=False,quiet=True,skip_interpolation=False) +lk,rbn,ubn,vbn,ep=flwr.pack_rotation_arrays(meta,rho,ct,ctV) +Pv=Psig.manual_copy() +for k,v in [('phi',RA),('theta',DEC),('incl',INCL),('phiref',PHIREF),('psi',PSI),('dist',DIST)]: setattr(Pv,k,np.ones(1)*v) +Pv.tref=event_time;Pv.deltaT=deltaT +tvals=np.arange(npts)*deltaT-0.02 +vec=flwr.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation(tvals,Pv,meta,lk,rbn,ubn,vbn,ep,Lmax=Lmax,array_output=True)[0] +mv=np.max(vec);vec_marg=mv+np.log(np.trapz(np.exp(vec-mv),dx=deltaT)) +worst = float(np.max(np.abs(vec - bf))) +print("brute-force(real F(t)) : peak=%.4f marg=%.4f" % (np.max(bf), bf_marg)) +print("vec rotation : peak=%.4f marg=%.4f" % (np.max(vec), vec_marg)) +print("max|vec_lnL_t - bruteforce_lnL_t| over window = %.3e" % worst) +assert worst < 1e-6, "vectorized rotation NoLoop disagrees with brute-force time-varying response: %g" % worst +print("PASS: vectorized rotation NoLoop == brute-force time-varying-response likelihood") diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_pathB.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_pathB.py new file mode 100644 index 000000000..a60cf89eb --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_pathB.py @@ -0,0 +1,66 @@ +""" +test_slowrot_pathB : Path B (delay-derivative) sanity checks. + + (1) scalar Path B (p_max=1) at f_sidereal=0 reduces EXACTLY to the baseline + FactoredLogLikelihood (the p>=1 coefficient sums vanish at Omega=0). + (2) scalar Path B (p_max=1) at the real sidereal rate respects the Cauchy-Schwarz + bound lnL <= 0.5 at the injection parameters. + (3) vectorized Path B (p_max=1) at f_sidereal=0 reduces to the baseline NoLoop. + +Rigorous validation of the p>=1 delay physics vs a brute force with TIME-VARYING delay on +a LONG signal is deferred to the systematic final-validation pass. +Run: PYTHONPATH=.../Code python RIFT/likelihood/test_slowrot_pathB.py +""" +from __future__ import print_function, division +import numpy as np, lal, lalsimulation as lalsim +import RIFT.lalsimutils as lsu +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.factored_likelihood_with_rotation as flwr + +fmin=30.;fmax=1700.;event_time=1e9;t_window=0.1;Lmax=2;deltaT=1/4096.;deltaF=1/4. +HB=tuple(range(-3,4)) +Psig=lsu.ChooseWaveformParams(fmin=fmin,radec=True,incl=0.3,phiref=0.0,theta=0.2,phi=1.0,psi=0.4, + m1=30*lal.MSUN_SI,m2=25*lal.MSUN_SI,detector='H1',dist=200e6*lal.PC_SI,deltaT=deltaT,tref=event_time,deltaF=deltaF) +data_dict={} +for d in ["H1","L1","V1"]: + P=Psig.manual_copy();P.detector=d;data_dict[d]=lsu.non_herm_hoff(P) +psd_dict={d:lalsim.SimNoisePSDaLIGOZeroDetHighPower for d in data_dict} + +def test_scalar_reduces_to_baseline(): + extr=lsu.ChooseWaveformParams(radec=True,phi=1.0,theta=0.2,psi=0.5,incl=0.7,phiref=0.9,tref=event_time,dist=300e6*lal.PC_SI) + rb=fl.PrecomputeLikelihoodTerms(event_time,t_window,Psig,data_dict,psd_dict,Lmax,fmax,analyticPSD_Q=True,verbose=False,quiet=True,ignore_threshold=None) + lnLb=fl.FactoredLogLikelihood(extr,rb[3],rb[0],rb[1],rb[2],Lmax) + rr=flwr.PrecomputeLikelihoodTermsWithRotation(event_time,t_window,Psig,data_dict,psd_dict,Lmax,fmax,harmonics=HB,p_max=1,f_sidereal=0.0,analyticPSD_Q=True,verbose=False,quiet=True,skip_interpolation=False) + lnLB=flwr.FactoredLogLikelihoodWithRotation(extr,rr[0],rr[1],rr[2],rr[4],Lmax) + print("(1) scalar PathB(p1,fsid0)=%.8f baseline=%.8f |d|=%.2e"%(lnLB,lnLb,abs(lnLB-lnLb))) + assert abs(lnLB-lnLb)<1e-6*(1+abs(lnLb)) + +def test_scalar_respects_bound(): + extr=lsu.ChooseWaveformParams(radec=True,phi=1.0,theta=0.2,psi=0.4,incl=0.3,phiref=0.0,tref=event_time,dist=200e6*lal.PC_SI) + rr=flwr.PrecomputeLikelihoodTermsWithRotation(event_time,t_window,Psig,{'H1':data_dict['H1']},{'H1':psd_dict['H1']},Lmax,fmax,harmonics=HB,p_max=1,f_sidereal=flwr.F_SIDEREAL,analyticPSD_Q=True,verbose=False,quiet=True,skip_interpolation=False) + lnL=flwr.FactoredLogLikelihoodWithRotation(extr,rr[0],rr[1],rr[2],rr[4],Lmax) + IP=lsu.ComplexIP(fmin,fmax,1/2./deltaT,deltaF,psd_dict['H1'],True,False,0.);dd=IP.ip(data_dict['H1'],data_dict['H1']).real + print("(2) scalar PathB(p1,real,H1) lnL=%.4f 0.5=%.4f"%(lnL,0.5*dd)) + assert lnL<=0.5*dd + +def test_vec_reduces_to_baseline(): + tvals=np.arange(400)*deltaT-0.02 + rb=fl.PrecomputeLikelihoodTerms(event_time,t_window,Psig,data_dict,psd_dict,Lmax,fmax,analyticPSD_Q=True,verbose=False,quiet=True,ignore_threshold=None) + lk={};ra={};cu={};cv={};ep={} + for d in data_dict: + a,b,c,U,V,rA,rI,e=fl.PackLikelihoodDataStructuresAsArrays(list(rb[3][d].keys()),None,rb[3][d],rb[1][d],rb[2][d]);lk[d]=a;ra[d]=rA;cu[d]=U;cv[d]=V;ep[d]=e + Pv=Psig.manual_copy() + for k,v in [('phi',1.0),('theta',0.2),('incl',0.7),('phiref',0.9),('psi',0.5),('dist',300e6*lal.PC_SI)]: setattr(Pv,k,np.ones(1)*v) + Pv.tref=event_time;Pv.deltaT=deltaT + bl=fl.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals,Pv,lk,ra,cu,cv,ep,Lmax=Lmax,xpy=np)[0] + rr0=flwr.PrecomputeLikelihoodTermsWithRotation(event_time,t_window,Psig,data_dict,psd_dict,Lmax,fmax,harmonics=HB,p_max=1,f_sidereal=0.0,analyticPSD_Q=True,verbose=False,quiet=True,skip_interpolation=False) + lkr,rar,uar,var,epr=flwr.pack_rotation_arrays(rr0[4],rr0[3],rr0[1],rr0[2]) + vv=flwr.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation(tvals,Pv,rr0[4],lkr,rar,uar,var,epr,Lmax=Lmax,array_output=False)[0] + print("(3) vec PathB(p1,fsid0)=%.8f baseline NoLoop=%.8f |d|=%.2e"%(vv,bl,abs(vv-bl))) + assert abs(vv-bl)<1e-7*(1+abs(bl)) + +if __name__=="__main__": + test_scalar_reduces_to_baseline() + test_scalar_respects_bound() + test_vec_reduces_to_baseline() + print("ALL SLOWROT PATH B CHECKS PASSED") diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_pathB_bruteforce.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_pathB_bruteforce.py new file mode 100644 index 000000000..262a552a3 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_pathB_bruteforce.py @@ -0,0 +1,53 @@ +from __future__ import print_function, division +import numpy as np, lal, lalsimulation as lalsim +from scipy.interpolate import CubicSpline +import RIFT.lalsimutils as lsu +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.factored_likelihood_with_rotation as flwr +import RIFT.likelihood.slowrot_response as srr +event_time=1e9; Lmax=2; t_window=0.1; det='H1' +psd=lalsim.SimNoisePSDaLIGOZeroDetHighPower; apx=lalsim.GetApproximantFromString("IMRPhenomD") +import os +OMEGA_INF=flwr.OMEGA_EARTH*float(os.environ.get("INFL","340")); FSID_INF=OMEGA_INF/(2*np.pi) +def _ifft(hf_d): + o={} + for lm,hf in hf_d.items(): + n=hf.data.length; dt=1./(n*hf.deltaF); ht=lal.CreateCOMPLEX16TimeSeries("h",hf.epoch,0.,dt,lal.DimensionlessUnit,n) + lal.COMPLEX16FreqTimeFFT(ht,hf,lal.CreateReverseCOMPLEX16FFTPlan(n,0)); o[lm]=ht + return o +def _to_fd(re,epoch,dt,N): + ht=lal.CreateCOMPLEX16TimeSeries("h",epoch,0.,dt,lal.DimensionlessUnit,N); ht.data.data[:]=re[:N] + hf=lal.CreateCOMPLEX16FrequencySeries("hf",epoch,0.,1./dt/N,lsu.lsu_HertzUnit,N) + lal.COMPLEX16TimeFreqFFT(hf,ht,lal.CreateForwardCOMPLEX16FFTPlan(N,0)); return hf +from scipy.interpolate import InterpolatedUnivariateSpline +def _peak(lt): + lt=np.asarray(lt,float); x=np.arange(len(lt)); sp=InterpolatedUnivariateSpline(x,lt,k=4) + xs=np.linspace(0,len(lt)-1,len(lt)*32); return float(np.max(sp(xs))) +fmin,fmax,deltaT,seglen=25.,512.,1/2048.,16.; deltaF=1./seglen; fNyq=1/2./deltaT; N=int(round(seglen/deltaT)) +RA,DEC,PSI,INCL,PHIREF=1.2,0.3,0.5,0.4,0.0; DLOUD=fl.distMpcRef*1e6*lsu.lsu_PC/30. +Psig=lsu.ChooseWaveformParams(fmin=fmin,radec=True,incl=INCL,phiref=PHIREF,theta=DEC,phi=RA,psi=PSI, + m1=2.2*lal.MSUN_SI,m2=1.8*lal.MSUN_SI,detector=det,dist=200e6*lal.PC_SI,deltaT=deltaT,tref=event_time,deltaF=deltaF); Psig.approx=apx +Pm=Psig.manual_copy(); Pm.dist=DLOUD +hlms_fd,_=fl.internal_hlm_generator(Pm,Lmax,verbose=False,quiet=True); hlmsT=_ifft(hlms_fd) +lm0=list(hlmsT.keys())[0]; nn=hlmsT[lm0].data.length; dt=hlmsT[lm0].deltaT; ep=float(hlmsT[lm0].epoch); tt=ep+np.arange(nn)*dt +Sig=np.zeros(nn,complex) +for lm in hlmsT: Sig+=hlmsT[lm].data.data*lal.SpinWeightedSphericalHarmonic(INCL,-PHIREF,-2,lm[0],lm[1]) +reS=CubicSpline(tt,Sig.real,extrapolate=False); imS=CubicSpline(tt,Sig.imag,extrapolate=False) +lald=lalsim.DetectorPrefixToLALDetector(det); g_ev=lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(event_time))-RA +A=srr.antenna_harmonics(lald.response,DEC,PSI); At={k:A[k]*np.exp(1j*k*g_ev) for k in A} +B=srr.delay_harmonics(lald.location,DEC); Bt={k:B[k]*np.exp(1j*k*g_ev) for k in B} +tau_t=np.real(sum(Bt[k]*np.exp(1j*k*OMEGA_INF*tt) for k in Bt)) +F_t=sum(At[k]*np.exp(1j*k*OMEGA_INF*tt) for k in At) +Sig_d=np.nan_to_num(reS(tt-tau_t)+1j*imS(tt-tau_t)) +data=_to_fd(np.real(F_t*Sig_d),lal.LIGOTimeGPS(float(hlmsT[lm0].epoch)+event_time),dt,N); data_dict={det:data}; psd_dict={det:psd} +IPc=lsu.ComplexIP(fmin,fmax,fNyq,data.deltaF,psd,True,False,0.); HALF_DD=0.5*IPc.ip(data,data).real +print("inflated seglen=%.0fs 0.5=%.4f"%(seglen,HALF_DD)) +Pv=Psig.manual_copy() +for k,v in [('phi',RA),('theta',DEC),('incl',INCL),('phiref',PHIREF),('psi',PSI),('dist',DLOUD)]: setattr(Pv,k,np.ones(1)*v) +Pv.tref=event_time; Pv.deltaT=deltaT; Nw=int(0.02/deltaT); tvals=np.arange(-Nw,Nw)*deltaT +for pmax in [0,1,2,3]: + nh=2+pmax + bk=flwr.PrecomputeLikelihoodTermsWithRotation(event_time,t_window,Psig,data_dict,psd_dict,Lmax,fmax,harmonics=tuple(range(-nh,nh+1)),p_max=pmax,f_sidereal=FSID_INF,analyticPSD_Q=True,verbose=False,quiet=True,skip_interpolation=True) + lk,rbn,ubn,vbn,epd=flwr.pack_rotation_arrays(bk[4],bk[3],bk[1],bk[2]) + lnL=_peak(flwr.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation(tvals,Pv,bk[4],lk,rbn,ubn,vbn,epd,Lmax=Lmax,array_output=True)[0]) + print(" p_max=%d : lnL=%.5f deficit=%.5f"%(pmax,lnL,HALF_DD-lnL)) diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_pathB_groundtruth.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_pathB_groundtruth.py new file mode 100644 index 000000000..7e2ad583f --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_pathB_groundtruth.py @@ -0,0 +1,69 @@ +from __future__ import print_function, division +import sys +import numpy as np +import lal, lalsimulation as lalsim +import RIFT.lalsimutils as lsu +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.factored_likelihood_with_rotation as flwr +event_time=1e9; Lmax=2; t_window=0.1; det='H1' +psd=lalsim.SimNoisePSDaLIGOZeroDetHighPower +apx=lalsim.GetApproximantFromString("IMRPhenomD") +SCALE=40. # loudness: data-mode distance = distMpcRef/SCALE + +def _ifft(hlms_fd): + out={} + for lm,hf in hlms_fd.items(): + n=hf.data.length; dt=1./(n*hf.deltaF) + ht=lal.CreateCOMPLEX16TimeSeries("h",hf.epoch,0.,dt,lal.DimensionlessUnit,n) + lal.COMPLEX16FreqTimeFFT(ht,hf,lal.CreateReverseCOMPLEX16FFTPlan(n,0)); out[lm]=ht + return out +def _fd(hk,TDlen): + if hk.data.length!=TDlen: hk=lal.ResizeREAL8TimeSeries(hk,0,TDlen) + n=TDlen + htC=lal.CreateCOMPLEX16TimeSeries("h",hk.epoch,hk.f0,hk.deltaT,hk.sampleUnits,n); htC.data.data[:]=hk.data.data + hf=lal.CreateCOMPLEX16FrequencySeries("hf",hk.epoch,hk.f0,1./hk.deltaT/n,lsu.lsu_HertzUnit,n) + lal.COMPLEX16TimeFreqFFT(hf,htC,lal.CreateForwardCOMPLEX16FFTPlan(n,0)); return hf +def _peak(lt,up=16): + lt=np.asarray(lt,dtype=float); N=len(lt) + L=np.fft.rfft(lt); Lp=np.zeros(N*up//2+1,dtype=complex); Lp[:len(L)]=L + return float(np.max(np.fft.irfft(Lp,N*up)*up)) + +def run(mode="short"): + if mode=="long": fmin,fmax,deltaT,seglen,m1,m2=25.,256.,1/2048.,32.,2.2*lal.MSUN_SI,1.8*lal.MSUN_SI + else: fmin,fmax,deltaT,seglen,m1,m2=30.,1700.,1/4096.,4.,30*lal.MSUN_SI,25*lal.MSUN_SI + deltaF=1./seglen; fNyq=1/2./deltaT; TDlen=int(round(seglen/deltaT)) + DLOUD=fl.distMpcRef*1e6*lsu.lsu_PC/SCALE + Psig=lsu.ChooseWaveformParams(fmin=fmin,radec=True,incl=0.4,phiref=0.0,theta=0.3,phi=1.2,psi=0.5, + m1=m1,m2=m2,detector=det,dist=200e6*lal.PC_SI,deltaT=deltaT,tref=event_time,deltaF=deltaF); Psig.approx=apx + Pm=Psig.manual_copy(); Pm.dist=DLOUD + hlms_fd,_=fl.internal_hlm_generator(Pm,Lmax,verbose=False,quiet=True) + hlmsT=_ifft(hlms_fd) + data=_fd(lsu.hoft_from_hlm(hlmsT,Psig.manual_copy()),TDlen) + data_dict={det:data}; psd_dict={det:psd} + IPc=lsu.ComplexIP(fmin,fmax,fNyq,data.deltaF,psd,True,False,0.); HALF_DD=0.5*IPc.ip(data,data).real + print("[%s] seglen=%.0fs data_len=%d 0.5=%.4f"%(mode,seglen,data.data.length,HALF_DD)) + RA,DEC,PSI,INCL,PHIREF,DIST=1.2,0.3,0.5,0.4,0.0,DLOUD + ri,ct,ctV,rho,snr,rest=fl.PrecomputeLikelihoodTerms(event_time,t_window,Psig,data_dict,psd_dict,Lmax,fmax,analyticPSD_Q=True,verbose=False,quiet=True,ignore_threshold=None) + Pe=Psig.manual_copy(); Pe.incl=INCL;Pe.phiref=PHIREF;Pe.psi=PSI;Pe.phi=RA;Pe.theta=DEC;Pe.dist=DIST + hk_gt=_fd(lsu.hoft_from_hlm(hlmsT,Pe),TDlen); lnL_gt=IPc.ip(data,hk_gt).real-0.5*IPc.ip(hk_gt,hk_gt).real + Pv=Psig.manual_copy() + for k,v in [('phi',RA),('theta',DEC),('incl',INCL),('phiref',PHIREF),('psi',PSI),('dist',DIST)]: + setattr(Pv,k,np.ones(1)*v) + Pv.tref=event_time; Pv.deltaT=deltaT + Nw=int(0.02/deltaT); tvals=np.arange(-Nw,Nw)*deltaT + lkB={};rAB={};cuB={};cvB={};epB={} + for d in data_dict: + a,b,c,U,V,rA,rI,e=fl.PackLikelihoodDataStructuresAsArrays(list(rho[d].keys()),None,rho[d],ct[d],ctV[d]) + lkB[d]=a;rAB[d]=rA;cuB[d]=U;cvB[d]=V;epB[d]=e + lnL_base=_peak(fl.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop(tvals,Pv,lkB,rAB,cuB,cvB,epB,Lmax=Lmax,xpy=np,return_lnLt=True)[0]) + def rot(pmax): + nh=2+pmax + bk=flwr.PrecomputeLikelihoodTermsWithRotation(event_time,t_window,Psig,data_dict,psd_dict,Lmax,fmax,harmonics=tuple(range(-nh,nh+1)),p_max=pmax,f_sidereal=flwr.F_SIDEREAL,analyticPSD_Q=True,verbose=False,quiet=True,skip_interpolation=True) + lk,rbn,ubn,vbn,ep=flwr.pack_rotation_arrays(bk[4],bk[3],bk[1],bk[2]) + return _peak(flwr.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation(tvals,Pv,bk[4],lk,rbn,ubn,vbn,ep,Lmax=Lmax,array_output=True)[0]) + lnL_A=rot(0); lnL_B=rot(1) + print(" ground truth (SimDetStrain) lnL = %.5f"%lnL_gt) + print(" baseline NoLoop peak lnL = %.5f deficit=%.5f"%(lnL_base,lnL_gt-lnL_base)) + print(" Path A NoLoop (F(t)) lnL = %.5f deficit=%.5f"%(lnL_A,lnL_gt-lnL_A)) + print(" Path B NoLoop (F(t)+p=1) lnL = %.5f deficit=%.5f"%(lnL_B,lnL_gt-lnL_B)) +if __name__=="__main__": run(sys.argv[1] if len(sys.argv)>1 else "short") diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_precompute_integration.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_precompute_integration.py new file mode 100644 index 000000000..62a77c062 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_precompute_integration.py @@ -0,0 +1,125 @@ +""" +Full-stack integration test for factored_likelihood_with_rotation on real waveforms/data. + +Run under the RIFT venv with this worktree first on PYTHONPATH, e.g. + + source ~/RIFT_develUWM/bin/activate + PYTHONPATH=~/RIFT_slowrot/MonteCarloMarginalizeCode/Code \ + python .../RIFT/likelihood/test_slowrot_precompute_integration.py + +Inputs are built exactly like RIFT.calmarg.test_precompute_alignment (a short BBH injection, +noiseless detector data via non_herm_hoff, analytic aLIGO PSD), i.e. the same machinery the +ILE-GPU-Paper CI demo exercises. + +Checks: + V0 (recovery): PrecomputeLikelihoodTermsWithRotation(harmonics=(0,), p_max=0) reproduces + the baseline PrecomputeLikelihoodTerms bit-for-bit -- Q, U, V. This validates all of + the new plumbing (mode generation, p=0 / n=0 identity, cross-term assembly). + Vmag (magnitude sanity): running Path A (harmonics -2..2) end-to-end on real data, the + a=(0,0) overlap still equals the baseline exactly, and the modulated a=(0,+-1) overlaps + differ from it by a small but nonzero amount of the expected order Omega*(signal + duration) -- confirming the sidereal modulation is applied, at the right (tiny) scale. + +Full physical validation (V1: the rotation-aware lnL vs a brute-force time-varying-response +likelihood) needs the lnL assembly that contracts this bank with A_n,B_n, and is the next +step. +""" +from __future__ import print_function, division + +import numpy as np +import lal +import lalsimulation as lalsim + +import RIFT.lalsimutils as lalsimutils +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.factored_likelihood_with_rotation as flwr + +fSample = 4096.0 +fmin = 30.0 +fmax = 1700.0 +event_time = 1000000000.0 +t_window = 0.1 +Lmax = 2 + +Psig = lalsimutils.ChooseWaveformParams( + fmin=fmin, radec=True, incl=0.3, phiref=0.0, theta=0.2, phi=1.0, psi=0.4, + m1=30 * lal.MSUN_SI, m2=25 * lal.MSUN_SI, + detector='H1', dist=200e6 * lal.PC_SI, deltaT=1. / fSample, + tref=event_time, deltaF=1. / 4.) + +data_dict = {} +for det in ("H1", "L1", "V1"): + P = Psig.manual_copy(); P.detector = det + data_dict[det] = lalsimutils.non_herm_hoff(P) +psd_dict = {det: lalsim.SimNoisePSDaLIGOZeroDetHighPower for det in data_dict} + + +def _maxrel(a, b): + a = np.asarray(a); b = np.asarray(b) + scale = max(np.max(np.abs(a)), np.max(np.abs(b)), 1e-300) + return float(np.max(np.abs(a - b)) / scale) + + +def _run_base(): + # ignore_threshold=None so no mode pruning -> same mode set as the rotation path + return fl.PrecomputeLikelihoodTerms( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + analyticPSD_Q=True, verbose=False, quiet=True, skip_interpolation=True, + ignore_threshold=None) + + +def _run_rot(harmonics, p_max): + return flwr.PrecomputeLikelihoodTermsWithRotation( + event_time, t_window, Psig, data_dict, psd_dict, Lmax, fmax, + harmonics=harmonics, p_max=p_max, analyticPSD_Q=True, + verbose=False, quiet=True, skip_interpolation=True) + + +def test_V0_recovers_baseline(): + _, ct_b, ctV_b, rho_b, _, _ = _run_base() + _, ct_r, ctV_r, rho_r, meta = _run_rot(harmonics=(0,), p_max=0) + a0 = (0, 0) + worst_q = worst_u = worst_v = 0.0 + for det in data_dict: + for mode in rho_b[det]: + worst_q = max(worst_q, _maxrel(rho_r[det][a0][mode].data.data, + rho_b[det][mode].data.data)) + for key in ct_b[det]: + worst_u = max(worst_u, abs(ct_r[det][(a0, a0)][key] - ct_b[det][key])) + worst_v = max(worst_v, abs(ctV_r[det][(a0, a0)][key] - ctV_b[det][key])) + print("V0 recovery: max rel|dQ|=%.2e max|dU|=%.2e max|dV|=%.2e" + % (worst_q, worst_u, worst_v)) + assert worst_q < 1e-8, "Q mismatch vs baseline: %g" % worst_q + assert worst_u < 1e-8 * (1 + worst_u), "U mismatch vs baseline: %g" % worst_u + assert worst_v < 1e-8 * (1 + worst_v), "V mismatch vs baseline: %g" % worst_v + + +def test_Vmag_modulation_applied_at_right_scale(): + _, ct_b, _, rho_b, _, _ = _run_base() + _, _, _, rho_r, meta = _run_rot(harmonics=(-2, -1, 0, 1, 2), p_max=0) + # a=(0,0) must still equal baseline exactly even with the extra harmonics present + worst_0 = 0.0 + for det in data_dict: + for mode in rho_b[det]: + worst_0 = max(worst_0, _maxrel(rho_r[det][(0, 0)][mode].data.data, + rho_b[det][mode].data.data)) + print("Vmag a=(0,0) vs baseline: max rel|dQ|=%.2e" % worst_0) + assert worst_0 < 1e-8, "modulated run corrupts the n=0 overlap: %g" % worst_0 + + # modulated harmonics differ from n=0 by ~ Omega * (signal duration), tiny but nonzero. + # Signal duration scale for a 30+25 BBH from 30 Hz is ~ 0.1-1 s -> Omega*T ~ 1e-5..1e-4. + for n in (1, -1, 2, -2): + worst_n = 0.0 + for det in data_dict: + for mode in rho_b[det]: + worst_n = max(worst_n, _maxrel(rho_r[det][(0, n)][mode].data.data, + rho_r[det][(0, 0)][mode].data.data)) + print("Vmag a=(0,%+d) vs a=(0,0): max rel|dQ|=%.2e" % (n, worst_n)) + assert worst_n > 1e-9, "harmonic n=%d has no effect (modulation not applied?)" % n + assert worst_n < 1e-2, "harmonic n=%d far too large (bug in modulation scale?): %g" % (n, worst_n) + + +if __name__ == "__main__": + test_V0_recovers_baseline() + test_Vmag_modulation_applied_at_right_scale() + print("ALL SLOWROT PRECOMPUTE INTEGRATION CHECKS PASSED") diff --git a/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_response.py b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_response.py new file mode 100644 index 000000000..7bdeed3e1 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/RIFT/likelihood/test_slowrot_response.py @@ -0,0 +1,135 @@ +""" +Validate slowrot_response against LAL over a full sidereal day. + +The harmonic decomposition is an exact algebraic identity in the hour angle g, so the +reconstruction must match lal.ComputeDetAMResponse and lal.TimeDelayFromEarthCenter to +machine precision (up to LAL's own internal rounding). Run directly: + + python test_slowrot_response.py + +or under pytest (the test_* functions assert). +""" +from __future__ import print_function, division + +import numpy as np + +import lal +import lalsimulation as lalsim + +# Load the leaf module directly by path so this test needs only numpy+lal, not the full +# RIFT package stack (whose __init__ pulls in glue/lalsimutils). When the package is +# properly installed, `import RIFT.likelihood.slowrot_response` is equivalent. +try: + import RIFT.likelihood.slowrot_response as sr +except Exception: + import importlib.util + import os + _here = os.path.dirname(os.path.abspath(__file__)) + _spec = importlib.util.spec_from_file_location( + "slowrot_response", os.path.join(_here, "slowrot_response.py")) + sr = importlib.util.module_from_spec(_spec) + _spec.loader.exec_module(sr) + +DETECTORS = ["H1", "L1", "V1", "K1"] +SIDEREAL_DAY = 86164.0905 # s +T0_GPS = 1000000000.0 # arbitrary reference epoch + +_TOL = 1e-11 # radians of F, seconds of tau -- essentially machine precision + + +def _lal_detector(det): + return lalsim.DetectorPrefixToLALDetector(det) + + +def _sample_gmst(n=512): + tgps = T0_GPS + np.linspace(0.0, SIDEREAL_DAY, n, endpoint=False) + gmst = np.array([lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(float(t))) for t in tgps]) + return tgps, gmst + + +def _max_antenna_error(det, ra, dec, psi): + lald = _lal_detector(det) + A = sr.antenna_harmonics(lald.response, dec, psi) + _, gmst = _sample_gmst() + # model + Fp_mod, Fc_mod = sr.antenna_response(A, gmst, ra) + # reference (LAL), sampled at the same GMST values + Fp_ref = np.empty_like(gmst) + Fc_ref = np.empty_like(gmst) + for i, g in enumerate(gmst): + Fp_ref[i], Fc_ref[i] = lal.ComputeDetAMResponse(lald.response, ra, dec, psi, float(g)) + return max(np.max(np.abs(Fp_mod - Fp_ref)), np.max(np.abs(Fc_mod - Fc_ref))) + + +def _max_delay_error(det, ra, dec): + lald = _lal_detector(det) + B = sr.delay_harmonics(lald.location, dec) + tgps, gmst = _sample_gmst() + g = sr.greenwich_hour_angle(gmst, ra) + tau_mod = sr.delay_from_harmonics(B, g) + tau_ref = np.array([ + lal.TimeDelayFromEarthCenter(lald.location, ra, dec, lal.LIGOTimeGPS(float(t))) + for t in tgps + ]) + return np.max(np.abs(tau_mod - tau_ref)) + + +# ---- deterministic parameter sweep ------------------------------------------------- +_RNG = np.random.RandomState(20260703) +_CASES = [] +for _det in DETECTORS: + for _ in range(5): + _ra = _RNG.uniform(0, 2 * np.pi) + _dec = np.arcsin(_RNG.uniform(-1, 1)) + _psi = _RNG.uniform(0, np.pi) + _CASES.append((_det, _ra, _dec, _psi)) + + +def test_antenna_harmonics_match_lal(): + worst = 0.0 + for det, ra, dec, psi in _CASES: + e = _max_antenna_error(det, ra, dec, psi) + worst = max(worst, e) + assert e < _TOL, "antenna mismatch %g for %s ra=%.3f dec=%.3f psi=%.3f" % ( + e, det, ra, dec, psi) + print("antenna: worst |dF| = %.3e" % worst) + + +def test_delay_harmonics_match_lal(): + worst = 0.0 + for det, ra, dec, psi in _CASES: + e = _max_delay_error(det, ra, dec) + worst = max(worst, e) + assert e < _TOL, "delay mismatch %g s for %s ra=%.3f dec=%.3f" % (e, det, ra, dec) + print("delay: worst |dtau| = %.3e s" % worst) + + +def test_only_five_and_three_harmonics(): + """Confirm there is genuinely no content beyond |n|=2 (antenna) / |n|=1 (delay).""" + lald = _lal_detector("H1") + ra, dec, psi = 1.0, 0.4, 0.7 + _, gmst = _sample_gmst(1024) + g = sr.greenwich_hour_angle(gmst, ra) + # antenna: DFT of F(g) sampled uniformly in g should have power only at |n|<=2 + gg = np.linspace(0, 2 * np.pi, 1024, endpoint=False) + F = sr.antenna_from_harmonics(sr.antenna_harmonics(lald.response, dec, psi), gg) + coeffs = np.fft.fft(F) / len(gg) + power = np.abs(coeffs) + keep = np.zeros_like(power, dtype=bool) + for n in (-2, -1, 0, 1, 2): + keep[n % len(gg)] = True + assert np.max(power[~keep]) < 1e-12, "antenna has power beyond |n|=2: %g" % np.max(power[~keep]) + tau = sr.delay_from_harmonics(sr.delay_harmonics(lald.location, dec), gg) + tcoef = np.abs(np.fft.fft(tau) / len(gg)) + keep = np.zeros_like(tcoef, dtype=bool) + for n in (-1, 0, 1): + keep[n % len(gg)] = True + assert np.max(tcoef[~keep]) < 1e-15, "delay has power beyond |n|=1: %g" % np.max(tcoef[~keep]) + print("harmonic content confirmed: antenna |n|<=2, delay |n|<=1") + + +if __name__ == "__main__": + test_antenna_harmonics_match_lal() + test_delay_harmonics_match_lal() + test_only_five_and_three_harmonics() + print("ALL SLOWROT RESPONSE CHECKS PASSED") diff --git a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode index 8ce0ee353..ad5aaabd0 100755 --- a/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode +++ b/MonteCarloMarginalizeCode/Code/bin/integrate_likelihood_extrinsic_batchmode @@ -254,6 +254,12 @@ optp.add_option("--extrinsic-proposal-adapt",action='store_true',default=False, optp.add_option("--vectorized", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, not LAL data structures. (Combine with --gpu to enable GPU use, where available)") optp.add_option("--gpu", action="store_true", help="Perform manipulations of lm and timeseries using numpy arrays, CONVERTING TO GPU when available. You MUST use this option with --vectorized (otherwise it is a no-op). You MUST have a suitable version of cupy installed, your cuda operational, etc") optp.add_option("--force-gpu-only", action="store_true", help="Hard fail if no GPU present (assessed by cupy not loading)") +optp.add_option("--rotation-slow", action="store_true", help="[Path A] Slow-rotation likelihood: account for the sidereal time-dependence of the antenna pattern F(t) over the signal (harmonic modulation). Requires --vectorized and CPU (not --gpu). Uses factored_likelihood_with_rotation.") +optp.add_option("--rotation-n-harmonics", type=int, default=2, help="Number of sidereal harmonics for --rotation-slow (antenna pattern needs 2, i.e. n=-2..2).") +optp.add_option("--rotation-p-max", type=int, default=0, help="[Path B] Max delay-derivative order for --rotation-slow (0 = amplitude drift only; >=1 adds propagation-delay drift).") +optp.add_option("--freqresponse", action="store_true", help="[Path D] Finite-size (frequency-dependent) detector-response likelihood: account for the finite light-travel-time transfer across the arms (matters for 3G/CE-ET). Requires --vectorized and CPU (not --gpu). Uses factored_likelihood_freqresponse (sky-harmonic route b, sky stays extrinsic).") +optp.add_option("--freqresponse-qmax", type=int, default=4, help="Highest power of the arm projection retained for --freqresponse (basis size Qmax+2). Higher for larger fL/c (heavier systems / higher fmax).") +optp.add_option("--freqresponse-arm-length", type=float, default=None, help="Arm-length override [m] for --freqresponse (e.g. 40000 for 40-km CE). Default: each detector's native LAL arm length.") optp.add_option("--force-xpy", action="store_true", help="Use the xpy code path. Use with --vectorized --gpu to use the fallback CPU-based code path. Useful for debugging.") optp.add_option("-o", "--output-file", help="Save result to this file.") optp.add_option("-O", "--output-format", default='xml', help="[xml|hdf5]") @@ -385,6 +391,26 @@ def _truthy_option(value): opts._noloop_time_interp = "cubic" if _truthy_option(opts.interpolate_time) else "nearest" +if opts.rotation_slow: + # Path A slow-rotation is currently wired only into the CPU vectorized branch. + if not opts.vectorized: + raise ValueError("--rotation-slow requires --vectorized") + if opts.gpu: + raise ValueError("--rotation-slow is CPU-only for now; do not combine with --gpu") + if getattr(opts, 'distance_marginalization', False) or getattr(opts, 'phase_marginalization', False): + raise ValueError("--rotation-slow does not yet support distance/phase marginalization") + +if opts.freqresponse: + # Path D finite-size response is currently wired only into the CPU vectorized branch. + if opts.rotation_slow: + raise ValueError("--freqresponse and --rotation-slow both replace the likelihood; use at most one") + if not opts.vectorized: + raise ValueError("--freqresponse requires --vectorized") + if opts.gpu: + raise ValueError("--freqresponse is CPU-only for now; do not combine with --gpu") + if getattr(opts, 'distance_marginalization', False) or getattr(opts, 'phase_marginalization', False): + raise ValueError("--freqresponse does not yet support distance/phase marginalization") + # cosmo d prior tools for interpolation: not used normally, but set if needed final_scipy_interpolate=None if 'cosmo' in opts.d_prior: @@ -431,6 +457,8 @@ if cupy_success and opts.gpu and opts.sampler_method == 'adaptive_cartesian_gpu' if not(opts.use_gwsignal): os.environ['RIFT_NO_GWSIGNAL'] = 'True' import RIFT.likelihood.factored_likelihood as factored_likelihood +import RIFT.likelihood.factored_likelihood_with_rotation as factored_likelihood_with_rotation +import RIFT.likelihood.factored_likelihood_freqresponse as factored_likelihood_freqresponse if opts.use_gwsignal and not(factored_likelihood.has_GWS): print(" HARD FAILURE: this node could not import gwsignal ! ") @@ -1822,6 +1850,39 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t ctVArrayDict[det] = cupy.asarray(ctVArrayDict[det]) epochDict[det] = cupy.asarray(epochDict[det]) + # [Path A] slow-rotation precompute: build the harmonic-indexed bank and pack it. + rotation_slow_data = None + if opts.rotation_slow: + _pmax = int(opts.rotation_p_max) + _nh = max(int(opts.rotation_n_harmonics), 2 + _pmax) # wide enough to cover all C_{(p,ntilde)} + _harm = tuple(range(-_nh, _nh + 1)) + _rint_r, _ct_r, _ctV_r, _rho_r, _meta_r = factored_likelihood_with_rotation.PrecomputeLikelihoodTermsWithRotation( + fiducial_epoch, t_window, P, data_dict, psd_dict, opts.l_max, fmax, + harmonics=_harm, p_max=_pmax, analyticPSD_Q=False, + inv_spec_trunc_Q=inv_spec_trunc_Q, T_spec=T_spec, + verbose=opts.verbose, quiet=not opts.verbose, skip_interpolation=True) + _lkR, _rhoN, _uNN, _vNN, _epR = factored_likelihood_with_rotation.pack_rotation_arrays( + _meta_r, _rho_r, _ct_r, _ctV_r) + rotation_slow_data = dict(meta=_meta_r, lookupNKDict=_lkR, rho_by_n=_rhoN, + U_by_nn=_uNN, V_by_nn=_vNN, epochDict=_epR) + print(" [rotation-slow] precompute complete; p_max", _pmax, "sidereal harmonics", _harm) + + # [Path D] finite-size (frequency-dependent) response precompute: fold each W_p(f) + # into the modes once and pack the response-basis overlap bank. + freqresponse_data = None + if opts.freqresponse: + _qmax = int(opts.freqresponse_qmax) + _rint_f, _ct_f, _ctV_f, _rho_f, _meta_f = factored_likelihood_freqresponse.PrecomputeLikelihoodTermsFreqResponse( + fiducial_epoch, t_window, P, data_dict, psd_dict, opts.l_max, fmax, + Qmax=_qmax, L_arm=opts.freqresponse_arm_length, analyticPSD_Q=False, + inv_spec_trunc_Q=inv_spec_trunc_Q, T_spec=T_spec, + verbose=opts.verbose, quiet=not opts.verbose, skip_interpolation=True) + _lkF, _rhoP, _uPP, _vPP, _epF = factored_likelihood_freqresponse.pack_freqresponse_arrays( + _meta_f, _rho_f, _ct_f, _ctV_f) + freqresponse_data = dict(meta=_meta_f, lookupNKDict=_lkF, rho_by_p=_rhoP, + U_by_pp=_uPP, V_by_pp=_vPP, epochDict=_epF) + print(" [freqresponse] precompute complete; Qmax", _qmax, "arm-length", opts.freqresponse_arm_length) + def _cal_error_probe(n_cal_now, n_start=256, n_cap=None, rel_tol=0.1): """Estimate (sigma_lnZ_cal, neff_cal, n_used, dist_mode): the calibration MC error of lnZ and the effective cal draw count, from per-realization @@ -2091,7 +2152,20 @@ def analyze_event(P_list, indx_event, data_dict, psd_dict, fmax, opts,inv_spec_t P.psi= psi_true P.phiref = phi_orb_true - lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVector(tvals, + if opts.rotation_slow: + lnL = factored_likelihood_with_rotation.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation( + tvals, P, rotation_slow_data['meta'], rotation_slow_data['lookupNKDict'], + rotation_slow_data['rho_by_n'], rotation_slow_data['U_by_nn'], + rotation_slow_data['V_by_nn'], rotation_slow_data['epochDict'], Lmax=opts.l_max, + time_interp=opts._noloop_time_interp) + elif opts.freqresponse: + lnL = factored_likelihood_freqresponse.DiscreteFactoredLogLikelihoodFreqResponseNoLoop( + tvals, P, freqresponse_data['meta'], freqresponse_data['lookupNKDict'], + freqresponse_data['rho_by_p'], freqresponse_data['U_by_pp'], + freqresponse_data['V_by_pp'], freqresponse_data['epochDict'], Lmax=opts.l_max, + time_interp=opts._noloop_time_interp) + else: + lnL = factored_likelihood.DiscreteFactoredLogLikelihoodViaArrayVector(tvals, P, lookupNKDict, rholmArrayDict, ctUArrayDict, ctVArrayDict,epochDict,Lmax=opts.l_max) # nEvals +=len(right_ascension) if supplemental_ln_likelihood: diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/Makefile b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/Makefile new file mode 100644 index 000000000..55d0e0f7a --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/Makefile @@ -0,0 +1,38 @@ +# demo/rift/slowrot -- verify-anywhere value demos for the slow-response likelihoods. +# +# make demo run BOTH quick-looks (no condor, no GPU); writes outputs/*.{txt,png} +# make demo-rotation Path A/B : rotation-vs-static lnL gain grows with signal duration Omega*T +# make demo-finite-size Path D : finite-size-vs-LWL lnL gain grows with detector arm length (fL/c) +# make clean remove outputs/ +# +# Both show the maintained vectorized (NoLoop) likelihood RECOVERS SNR the standard analysis +# loses, on the SAME data at the true parameters. They exercise the exact functions wired into +# integrate_likelihood_extrinsic_batchmode as --rotation-slow / --rotation-p-max / --freqresponse. +# Backing code: RIFT/likelihood/{factored_likelihood_with_rotation,factored_likelihood_freqresponse, +# slowrot_freqresponse}.py; validation RIFT/likelihood/test_slowrot_*.py; SLOWROT_HANDOFF.md. + +PYTHON ?= python +# Code root is four levels up (demo/rift/slowrot -> demo/rift -> demo -> Code); put it on +# PYTHONPATH so `import RIFT...` resolves this worktree even without an editable install. +CODE_ROOT := $(abspath $(dir $(lastword $(MAKEFILE_LIST)))/../../..) +RUN = PYTHONPATH=$(CODE_ROOT):$$PYTHONPATH $(PYTHON) + +.PHONY: help demo demo-rotation demo-finite-size clean +help: + @echo "targets:" + @echo " demo both verify-anywhere quick-looks (rotation + finite-size)" + @echo " demo-rotation Path A/B: gain vs signal duration (Omega*T); null control at short T" + @echo " demo-finite-size Path D: gain vs arm length (fL/c); null control at LIGO 4 km" + @echo " clean remove outputs/" + @echo "(uses PYTHON=$(PYTHON); CODE_ROOT=$(CODE_ROOT))" + +demo: demo-rotation demo-finite-size + +demo-rotation: + $(RUN) demo_rotation.py + +demo-finite-size: + $(RUN) demo_finite_size.py + +clean: + rm -rf outputs diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/README.md b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/README.md new file mode 100644 index 000000000..adedaf8fe --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/README.md @@ -0,0 +1,85 @@ +# demo/rift/slowrot — slow-response likelihood value demos + +Verify-anywhere (no condor, no GPU) quick-looks that RIFT's **time-/frequency-dependent +detector-response** likelihoods *add value*, not merely run. Both generalize the maintained +marginalized ILE likelihood so a long or long-armed signal is modeled with the correct +detector response instead of a single constant antenna pattern. + +``` +make demo # both demos below -> outputs/*.{txt,png} +make demo-rotation # Path A/B (Earth rotation) +make demo-finite-size # Path D (finite-size / frequency-dependent response) +``` + +Each demo builds data that carries the true response, then compares, **at the true +parameters on the same data and time grid**, the recovered time-maximized log-likelihood of +the standard analysis vs the response-aware one. The difference `gain = lnL_aware − lnL_standard` +is the SNR the standard analysis throws away. Both use the maintained vectorized (NoLoop) +path with cubic sub-bin time interpolation, so the peak-resolution floor cancels in the gain. + +--- + +## Path A/B — Earth rotation (`demo_rotation.py`) + +Over a long signal the antenna pattern `F(t)` drifts as the Earth rotates; the static +(Earth-fixed-response) likelihood loses match. RIFT injections already carry the true +time-varying response (`hoft → SimDetectorStrainREAL8TimeSeries`), so a standard injection is +"rotating". The gain grows with the rotation phase `Ω⊕·T` over the signal; a short signal is a +null control. Single H1, SNR≈30 held fixed so the gain tracks `Ω⊕·T`, not loudness: + +| config | seglen | Ω⊕·T | gain `lnL_rot − lnL_static` | +|---|---|---|---| +| null_bbh (30+25) | 2 s | 1.5e-4 | +0.004 (null control) | +| bbh_8_8 | 16 s | 1.2e-3 | +0.043 | +| bbh_4_4 | 64 s | 4.7e-3 | +0.187 | + +Grows ~50× from the short null to a minute-scale signal; extrapolated to a 90-minute XG BNS +(`Ω⊕·T`≈0.4) it is orders of magnitude larger — the 3G headline. Figure: +`outputs/rotation_gain_vs_duration.png`. ILE: `--rotation-slow` (Path A), +`--rotation-slow --rotation-p-max 2` (Path B, adds propagation-delay drift). + +## Path D — finite-size / frequency-dependent response (`demo_finite_size.py`) + +On a multi-km arm the light-travel time across the arm is not negligible vs the GW period, so +the response is per-frequency: `h_k(f) = F₊(f;sky) h₊(f) + Fₓ(f;sky) hₓ(f)`. We inject an exact +finite-size signal (`antenna_response_fd`) and compare the standard long-wavelength (LWL, +constant-response) likelihood to the finite-size one (`--freqresponse`, sky-harmonic route (b), +sky stays extrinsic). The gain grows with the arm length — i.e. with `fL/c`, the in-band +light-crossing phase. 15+13 M☉, fmax=2000 Hz, loud (SNR≈320), Qmax=6: + +| detector | arm L | fL/c @ fmax | gain `lnL_finite − lnL_LWL` | +|---|---|---|---| +| LIGO | 4 km | 0.027 | +0.25 (null control) | +| ET | 10 km | 0.067 | +3.5 | +| CE | 20 km | 0.133 | +12.4 | +| CE | 40 km | 0.267 | +39.6 | + +Null at the 4-km LIGO arm; tens of nats at a 40-km Cosmic Explorer. Only the +**direction-dependent** part of the response contributes — the common `e^{−i2πfL/c}` +light-crossing delay is degenerate with arrival time and is absorbed by both likelihoods' +time maximization, so it does not inflate the gain. Figure: +`outputs/finite_size_gain_vs_arm.png`. ILE: `--freqresponse` +`--freqresponse-arm-length 40000` `--freqresponse-qmax 6`. + +Higher `fL/c` (heavier system / higher fmax / longer arm) needs higher `--freqresponse-qmax` +(the response is a power series in `fL/c · (â·n̂)`); Qmax=6 converges through CE-40 km at this +config. + +--- + +## Full injection–recovery PE (cluster) + +These quick-looks are the likelihood-level core. The headline parameter-bias / single-network +sky-localization figures are DAG PE runs (structurally the standard RIFT pipeline with the +extra ILE flag appended: `--rotation-slow` / `--freqresponse`). The three analyses per event +differ *only* in that appended option. See `RIFT/likelihood/SLOWROT_HANDOFF.md` and +`~/RIFT_roboto_paper/analyses/slowrot_demo/` for the pipeline scaffolding. + +## Validation (the demos show value; these prove correctness) + +`RIFT/likelihood/test_slowrot_*.py` (run under the venv with this Code dir on `PYTHONPATH`): +- `test_slowrot_noloop.py` — rotation NoLoop == baseline at f_sidereal=0 (3.6e-12). +- `test_slowrot_headtohead.py` — matched-sample rotation vs baseline; Cauchy-Schwarz bound. +- `test_slowrot_freqresponse_likelihood.py` — finite-size: L→0 reduces to baseline (3e-9), + bound respected, and a **V4 positive control** asserting finite-size beats LWL by +38.9 nats + in an in-band-effect config. diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/demo_finite_size.py b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/demo_finite_size.py new file mode 100644 index 000000000..eee5ec0c7 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/demo_finite_size.py @@ -0,0 +1,221 @@ +#!/usr/bin/env python +""" +slowrot demo / demo_finite_size.py -- verify-anywhere (no condor, no GPU) quick-look that the +finite-size (frequency-dependent detector response) likelihood [Path D] ADDS VALUE on a +long-armed 3G detector. + +Beyond the long-wavelength (LWL) limit the detector strain is a per-frequency antenna +pattern, h_k(f) = F_+(f;sky) h_+(f) + F_x(f;sky) h_x(f), because the light-travel time across +a multi-km arm is no longer negligible vs the GW period. We build an EXACT finite-size +injection (from the same modes the likelihood uses -> internal_hlm_generator -> IFFT -> +antenna_response_fd) and, at the true parameters, compare the recovered log-likelihood of + + lnL_LWL : the standard RIFT likelihood (constant, frequency-independent response) + lnL_finite : the finite-size likelihood (--freqresponse, sky-harmonic route (b)) + +both TIME-MAXIMIZED on the SAME data and SAME time grid. The gain + + gain(L) = lnL_finite(truth) - lnL_LWL(truth) + +is ~0 for a LIGO 4-km arm (null control: the finite-size effect is below the floor) and grows +monotonically with the arm length -- i.e. with fL/c, the in-band light-crossing phase -- to +tens of nats for a 40-km Cosmic Explorer. That is exactly the SNR the LWL analysis throws +away on a 3G detector, and what --freqresponse recovers. + +The effect that matters is the DIRECTION-DEPENDENT part of the response; the common +e^{-i2pi f L/c} light-crossing delay is degenerate with the arrival time and is absorbed by +BOTH likelihoods' time maximization, so it does not inflate the gain. + +Run: python demo_finite_size.py (writes outputs/finite_size_gain_vs_arm.{txt,png}) + +Backing code: branch rift_slowrot; --freqresponse (+ --freqresponse-qmax/-arm-length) on +integrate_likelihood_extrinsic_batchmode; RIFT/likelihood/factored_likelihood_freqresponse.py, +slowrot_freqresponse.py. See RIFT/likelihood/SLOWROT_HANDOFF.md. +""" +from __future__ import print_function, division +import os +import numpy as np +import lal +import lalsimulation as lalsim + +import RIFT.lalsimutils as lsu +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.slowrot_freqresponse as sfr +import RIFT.likelihood.factored_likelihood_freqresponse as flfr + +EVENT_TIME = 1e9 +LMAX = 2 +DET = "H1" +PSD = lalsim.SimNoisePSDaLIGOZeroDetHighPower +# sky/pol/masses where the in-band finite-size effect is resolvable (heavier system -> high-f +# power; loud so the effect sits above the peak-resolution floor). +RA, DEC, PSI, INCL, PHIREF = 1.2, 0.3, 0.5, 0.4, 0.0 +M1, M2 = 15.0, 13.0 +FMAX = 2000.0 +SCALE = 40.0 # loudness: data-mode distance = distMpcRef/SCALE (SNR ~ 320) +SEGLEN = 8.0 +QMAX = 6 # response-basis order; higher needed as fL/c grows + +# (label, arm-length[m]) -- LIGO 4 km is the null control; ET ~10 km; CE 20/40 km. +CONFIGS = [ + ("LIGO_4km", 4000.0), + ("ET_10km", 10000.0), + ("CE_20km", 20000.0), + ("CE_40km", 40000.0), +] + +C_SI = sfr.C_SI + + +# --- minimal FD helpers (self-contained; same conventions as the validation test) --- +def _ifft(hf_d): + out = {} + for lm, hf in hf_d.items(): + n = hf.data.length; dt = 1. / (n * hf.deltaF) + ht = lal.CreateCOMPLEX16TimeSeries("h", hf.epoch, 0., dt, lal.DimensionlessUnit, n) + lal.COMPLEX16FreqTimeFFT(ht, hf, lal.CreateReverseCOMPLEX16FFTPlan(n, 0)); out[lm] = ht + return out + + +def _fwd_fd(re_series, epoch, dt, N): + ht = lal.CreateCOMPLEX16TimeSeries("h", epoch, 0., dt, lal.DimensionlessUnit, N) + ht.data.data[:] = re_series[:N] + hf = lal.CreateCOMPLEX16FrequencySeries("hf", epoch, 0., 1. / dt / N, lsu.lsu_HertzUnit, N) + lal.COMPLEX16TimeFreqFFT(hf, ht, lal.CreateForwardCOMPLEX16FFTPlan(N, 0)); return hf + + +def _wrap_fd(arr, epoch, deltaF): + n = len(arr) + hf = lal.CreateCOMPLEX16FrequencySeries("hf", epoch, 0., deltaF, lsu.lsu_HertzUnit, n) + hf.data.data[:] = arr + return hf + + +def _rev_td(hf): + n = hf.data.length; dt = 1. / (n * hf.deltaF) + ht = lal.CreateCOMPLEX16TimeSeries("h", hf.epoch, 0., dt, lal.DimensionlessUnit, n) + lal.COMPLEX16FreqTimeFFT(ht, hf, lal.CreateReverseCOMPLEX16FFTPlan(n, 0)); return ht + + +def _peak(lt): + """Time-maximized value via a fine spline over the sampled lnL(t) window.""" + from scipy.interpolate import InterpolatedUnivariateSpline + lt = np.asarray(lt, float); x = np.arange(len(lt)) + sp = InterpolatedUnivariateSpline(x, lt, k=4) + xs = np.linspace(0, len(lt) - 1, len(lt) * 32) + return float(np.max(sp(xs))) + + +def run_config(label, L_arm): + deltaT = 1.0 / (2 * FMAX); fmin = 30.0; deltaF = 1.0 / SEGLEN + fNyq = 1.0 / 2.0 / deltaT; t_window = 0.1 + DLOUD = fl.distMpcRef * 1e6 * lsu.lsu_PC / SCALE + + Psig = lsu.ChooseWaveformParams( + fmin=fmin, radec=True, incl=INCL, phiref=PHIREF, theta=DEC, phi=RA, psi=PSI, + m1=M1 * lal.MSUN_SI, m2=M2 * lal.MSUN_SI, detector=DET, dist=200e6 * lal.PC_SI, + deltaT=deltaT, tref=EVENT_TIME, deltaF=deltaF) + Psig.approx = lalsim.GetApproximantFromString("IMRPhenomD") + Pm = Psig.manual_copy(); Pm.dist = DLOUD + + # --- exact finite-size injection from the SAME modes --- + hlms_fd, _ = fl.internal_hlm_generator(Pm, LMAX, verbose=False, quiet=True) + hlmsT = _ifft(hlms_fd); lm0 = list(hlmsT.keys())[0] + nn = hlmsT[lm0].data.length; dt = hlmsT[lm0].deltaT; e0 = float(hlmsT[lm0].epoch) + Sig = np.zeros(nn, complex) + for lm in hlmsT: + Sig += hlmsT[lm].data.data * lal.SpinWeightedSphericalHarmonic(INCL, -PHIREF, -2, lm[0], lm[1]) + dt_geo = float(fl.ComputeArrivalTimeAtDetector(DET, RA, DEC, EVENT_TIME)) - EVENT_TIME + data_epoch = lal.LIGOTimeGPS(e0 + EVENT_TIME) + fvals = flfr.evaluate_fvals_from_length(nn, deltaF) + Sig_fd = _fwd_fd(Sig, data_epoch, dt, nn).data.data + Sig_del = _rev_td(_wrap_fd(Sig_fd * np.exp(-1j * 2 * np.pi * fvals * dt_geo), data_epoch, deltaF)).data.data + hpf = _fwd_fd(np.real(Sig_del), data_epoch, dt, nn).data.data + hcf = _fwd_fd(-np.imag(Sig_del), data_epoch, dt, nn).data.data + gmst = float(lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(EVENT_TIME))) + Fp, Fc = sfr.antenna_response_fd(DET, RA, DEC, PSI, fvals, gmst=gmst, L_arm=L_arm) + data = _wrap_fd(Fp * hpf + Fc * hcf, data_epoch, deltaF) + data_dict = {DET: data}; psd_dict = {DET: PSD} + IP = lsu.ComplexIP(fmin, FMAX, fNyq, data.deltaF, PSD, True, False, 0.) + snr = np.sqrt(IP.ip(data, data).real) + + Pv = Psig.manual_copy() + for k, v in [('phi', RA), ('theta', DEC), ('incl', INCL), ('phiref', PHIREF), + ('psi', PSI), ('dist', DLOUD)]: + setattr(Pv, k, np.ones(1) * v) + Pv.tref = EVENT_TIME; Pv.deltaT = deltaT + Nw = int(0.02 / deltaT); tvals = np.arange(-Nw, Nw) * deltaT + + # baseline (long-wavelength) NoLoop, time-maximized, at the truth + ri, ct, ctV, rho, _, _ = fl.PrecomputeLikelihoodTerms( + EVENT_TIME, t_window, Psig, data_dict, psd_dict, LMAX, FMAX, + analyticPSD_Q=True, verbose=False, quiet=True, ignore_threshold=None) + lk = {}; rA = {}; cu = {}; cv = {}; ep = {} + for d in data_dict: + a, b, c, U, V, r, rI, e = fl.PackLikelihoodDataStructuresAsArrays( + list(rho[d].keys()), None, rho[d], ct[d], ctV[d]) + lk[d] = a; rA[d] = r; cu[d] = U; cv[d] = V; ep[d] = e + lnL_lwl = _peak(fl.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop( + tvals, Pv, lk, rA, cu, cv, ep, Lmax=LMAX, xpy=np, return_lnLt=True, time_interp='cubic')[0]) + + # finite-size (Path D) NoLoop, time-maximized, at the truth + bk = flfr.PrecomputeLikelihoodTermsFreqResponse( + EVENT_TIME, t_window, Psig, data_dict, psd_dict, LMAX, FMAX, + Qmax=QMAX, L_arm=L_arm, analyticPSD_Q=True, verbose=False, quiet=True, skip_interpolation=True) + lkf, rbp, ubp, vbp, epf = flfr.pack_freqresponse_arrays(bk[4], bk[3], bk[1], bk[2]) + lnL_fin = _peak(flfr.DiscreteFactoredLogLikelihoodFreqResponseNoLoop( + tvals, Pv, bk[4], lkf, rbp, ubp, vbp, epf, Lmax=LMAX, array_output=True, time_interp='cubic')[0]) + + fLc = FMAX * L_arm / C_SI + return dict(label=label, L_arm=L_arm, fLc=fLc, snr=float(snr), + lnL_lwl=lnL_lwl, lnL_fin=lnL_fin, gain=lnL_fin - lnL_lwl) + + +def main(): + here = os.path.dirname(os.path.abspath(__file__)) + outdir = os.path.join(here, "outputs"); os.makedirs(outdir, exist_ok=True) + print("finite-size (Path D) value demo: %g+%g Msun, fmax=%g Hz, loudness SCALE=%g (SNR per row), Qmax=%d\n" + % (M1, M2, FMAX, SCALE, QMAX)) + rows = [] + for cfg in CONFIGS: + try: + r = run_config(*cfg) + except Exception as e: + print("%-9s SKIPPED (%s)" % (cfg[0], e)); continue + rows.append(r) + print("%-9s L=%6.0fm fL/c@fmax=%.3f SNR=%6.1f lnL_LWL=%11.3f lnL_finite=%11.3f GAIN=%+8.3f" + % (r['label'], r['L_arm'], r['fLc'], r['snr'], r['lnL_lwl'], r['lnL_fin'], r['gain'])) + if not rows: + print("no configs succeeded"); return + txt = os.path.join(outdir, "finite_size_gain_vs_arm.txt") + with open(txt, "w") as f: + f.write("# finite-size (Path D) demo: --freqresponse recovers SNR the long-wavelength analysis loses\n") + f.write("# label L_arm[m] fL/c@fmax SNR lnL_LWL lnL_finite gain=lnL_finite-lnL_LWL\n") + for r in rows: + f.write("%s %g %g %g %g %g %g\n" % (r['label'], r['L_arm'], r['fLc'], + r['snr'], r['lnL_lwl'], r['lnL_fin'], r['gain'])) + print("\nwrote", txt) + try: + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + Lkm = [r['L_arm'] / 1e3 for r in rows]; gain = [r['gain'] for r in rows] + labs = [r['label'] for r in rows] + fig, ax = plt.subplots(figsize=(5, 3.4)) + ax.plot(Lkm, gain, 'o-', color="#9d0208") + for x, y, l in zip(Lkm, gain, labs): + ax.annotate(l, (x, y), textcoords="offset points", xytext=(4, 4), fontsize=8) + ax.set_xlabel(r"detector arm length $L$ [km] ($\propto fL/c$, in-band light-crossing)") + ax.set_ylabel(r"$\ln\mathcal{L}_{\rm finite}-\ln\mathcal{L}_{\rm LWL}$ (recovered)") + ax.set_title("Finite-size detector response recovers SNR on 3G arms") + ax.axhline(0, color="0.7", lw=0.8) + fig.tight_layout() + png = os.path.join(outdir, "finite_size_gain_vs_arm.png") + fig.savefig(png, dpi=140) + print("wrote", png) + except Exception as e: + print("(plot skipped:", e, ")") + + +if __name__ == "__main__": + main() diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/demo_rotation.py b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/demo_rotation.py new file mode 100644 index 000000000..d1c7bc7c2 --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/demo_rotation.py @@ -0,0 +1,162 @@ +#!/usr/bin/env python +""" +slowrot demo / demo_rotation.py -- verify-anywhere (no condor, no GPU) quick-look that the +slow-rotation likelihood [Path A/B] ADDS VALUE: on a long signal the static +(Earth-fixed-response) likelihood LOSES match with the data (which carries the true +time-varying antenna pattern, since RIFT injections go through +SimDetectorStrainREAL8TimeSeries), and --rotation-slow recovers it. The recovered +log-likelihood at the true parameters, + + gain(config) = lnL_rotation(truth) - lnL_static(truth), (both time-marginalized) + +grows with the rotation phase Omega*T over the signal; a short signal is a null control +(gain ~ 0). This is the likelihood-level core of the full injection-recovery PE; the headline +sky-localization / parameter-bias figure is a cluster run. + +Both likelihoods are the maintained vectorized (NoLoop) path evaluated on the SAME data and +SAME time grid, so the gain is a clean difference (peak-resolution floor cancels). + +Run: python demo_rotation.py (writes outputs/rotation_gain_vs_duration.{txt,png}) + +Backing code: branch rift_slowrot; --rotation-slow (Path A) and --rotation-p-max N (Path B) on +integrate_likelihood_extrinsic_batchmode. See RIFT/likelihood/SLOWROT_HANDOFF.md. +""" +from __future__ import print_function, division +import os +import numpy as np +import lal +import lalsimulation as lalsim + +import RIFT.lalsimutils as lsu +import RIFT.likelihood.factored_likelihood as fl +import RIFT.likelihood.factored_likelihood_with_rotation as flwr + +EVENT_TIME = 1e9 +LMAX = 2 +DET = "H1" # single detector: rotation's sky-info value is clearest here +PSD = lalsim.SimNoisePSDaLIGOZeroDetHighPower +# sky/pol chosen where the antenna pattern drifts appreciably over a sidereal day +RA, DEC, PSI, INCL, PHIREF = 1.0, -0.6, 0.4, 0.5, 0.0 +HARM = (-2, -1, 0, 1, 2) + +# (label, m1, m2, fmin, fmax, dist[Mpc]) -- increasing duration via lower mass / fmin. +# Segment length is AUTO-FIT to each waveform (next power of two). Distances tuned for a +# comparable (loud) SNR (~30) so the gain tracks Omega*T, not loudness. +CONFIGS = [ + ("null_bbh", 30.0, 25.0, 30.0, 1024.0, 900.0), # short BBH, ~2 s: control, gain~0 + ("bbh_8_8", 8.0, 8.0, 25.0, 1024.0, 380.0), # ~16 s + ("bbh_4_4", 4.0, 4.0, 22.0, 1024.0, 210.0), # ~64 s -> larger Omega*T +] +OMEGA_EARTH = flwr.OMEGA_EARTH + + +def _P(m1, m2, fmin, deltaT, deltaF): + P = lsu.ChooseWaveformParams( + fmin=fmin, radec=True, incl=INCL, phiref=PHIREF, theta=DEC, phi=RA, psi=PSI, + m1=m1 * lal.MSUN_SI, m2=m2 * lal.MSUN_SI, detector=DET, + dist=200e6 * lal.PC_SI, deltaT=deltaT, tref=EVENT_TIME, deltaF=deltaF) + P.approx = lalsim.GetApproximantFromString("IMRPhenomD") + return P + + +def _P_vec(Psig, dist_SI, K=1): + Pv = Psig.manual_copy() + Pv.phi = np.ones(K) * RA; Pv.theta = np.ones(K) * DEC; Pv.psi = np.ones(K) * PSI + Pv.incl = np.ones(K) * INCL; Pv.phiref = np.ones(K) * PHIREF + Pv.dist = np.ones(K) * dist_SI + Pv.tref = float(EVENT_TIME); Pv.deltaT = Psig.deltaT + return Pv + + +def run_config(label, m1, m2, fmin, fmax, dist_Mpc): + deltaT = 1.0 / (2 * fmax) if fmax > 1024 else 1.0 / 2048.0 + dist_SI = dist_Mpc * 1e6 * lal.PC_SI + t_window = 0.1 + # AUTO-FIT seglen: generate with deltaF=None so non_herm_hoff pads to the next power of + # two >= the waveform length; then read back the segment length / deltaF. + Psig = _P(m1, m2, fmin, deltaT, None) + Psig.dist = dist_SI + data = lsu.non_herm_hoff(Psig) # RIFT injection -> carries the true time-varying response + seglen = data.data.length * deltaT + Psig.deltaF = data.deltaF + data_dict = {DET: data}; psd_dict = {DET: PSD} + fNyq = 1.0 / (2 * deltaT) + IP = lsu.ComplexIP(fmin, fmax, fNyq, data.deltaF, PSD, True, False, 0.) + snr = np.sqrt(IP.ip(data, data).real) + Pv = _P_vec(Psig, dist_SI) + Nw = int(0.02 / deltaT); tvals = np.arange(-Nw, Nw) * deltaT + + # baseline (static) NoLoop, time-marginalized, at the truth + rib, ctb, ctVb, rhob, _, _ = fl.PrecomputeLikelihoodTerms( + EVENT_TIME, t_window, Psig, data_dict, psd_dict, LMAX, fmax, + analyticPSD_Q=True, verbose=False, quiet=True, ignore_threshold=None) + lk = {}; rA = {}; cu = {}; cv = {}; ep = {} + for d in data_dict: + a, b, c, U, V, rAr, rI, e = fl.PackLikelihoodDataStructuresAsArrays( + list(rhob[d].keys()), None, rhob[d], ctb[d], ctVb[d]) + lk[d] = a; rA[d] = rAr; cu[d] = U; cv[d] = V; ep[d] = e + lnL_static = float(fl.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoop( + tvals, Pv, lk, rA, cu, cv, ep, Lmax=LMAX, xpy=np, time_interp='cubic')[0]) + + # rotation (Path A) NoLoop at the truth + bk = flwr.PrecomputeLikelihoodTermsWithRotation( + EVENT_TIME, t_window, Psig, data_dict, psd_dict, LMAX, fmax, + harmonics=HARM, p_max=0, f_sidereal=flwr.F_SIDEREAL, analyticPSD_Q=True, + verbose=False, quiet=True, skip_interpolation=True) + lkr, rbn, ubn, vbn, epr = flwr.pack_rotation_arrays(bk[4], bk[3], bk[1], bk[2]) + lnL_rot = float(flwr.DiscreteFactoredLogLikelihoodViaArrayVectorNoLoopWithRotation( + tvals, Pv, bk[4], lkr, rbn, ubn, vbn, epr, Lmax=LMAX, time_interp='cubic')[0]) + + OmegaT = OMEGA_EARTH * seglen + return dict(label=label, seglen=seglen, snr=float(snr), OmegaT=OmegaT, + lnL_static=lnL_static, lnL_rot=lnL_rot, gain=lnL_rot - lnL_static) + + +def main(): + here = os.path.dirname(os.path.abspath(__file__)) + outdir = os.path.join(here, "outputs"); os.makedirs(outdir, exist_ok=True) + print("slow-rotation (Path A) value demo: single %s, SNR~30 held fixed so gain tracks Omega*T\n" % DET) + rows = [] + for cfg in CONFIGS: + try: + r = run_config(*cfg) + except Exception as e: + print("%-9s SKIPPED (%s)" % (cfg[0], e)); continue + rows.append(r) + print("%-9s seglen=%5.0fs SNR=%6.1f Omega*T=%.2e lnL_static=%12.4f lnL_rot=%12.4f GAIN=%+.4f" + % (r['label'], r['seglen'], r['snr'], r['OmegaT'], + r['lnL_static'], r['lnL_rot'], r['gain'])) + if not rows: + print("no configs succeeded"); return + txt = os.path.join(outdir, "rotation_gain_vs_duration.txt") + with open(txt, "w") as f: + f.write("# slow-rotation demo: --rotation-slow recovers SNR the static analysis loses on long signals\n") + f.write("# label seglen[s] SNR Omega*T lnL_static lnL_rot gain=lnL_rot-lnL_static\n") + for r in rows: + f.write("%s %g %g %g %g %g %g\n" % (r['label'], r['seglen'], r['snr'], + r['OmegaT'], r['lnL_static'], r['lnL_rot'], r['gain'])) + print("\nwrote", txt) + try: + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + OmT = [r['OmegaT'] for r in rows]; gain = [r['gain'] for r in rows] + labs = [r['label'] for r in rows] + fig, ax = plt.subplots(figsize=(5, 3.4)) + ax.plot(OmT, gain, 'o-', color="#2a6f97") + for x, y, l in zip(OmT, gain, labs): + ax.annotate(l, (x, y), textcoords="offset points", xytext=(4, 4), fontsize=8) + ax.set_xlabel(r"$\Omega_\oplus T$ (rotation phase over the signal)") + ax.set_ylabel(r"$\ln\mathcal{L}_{\rm rot}-\ln\mathcal{L}_{\rm static}$ (recovered)") + ax.set_title("Accounting for Earth rotation recovers SNR on long signals") + ax.axhline(0, color="0.7", lw=0.8) + fig.tight_layout() + png = os.path.join(outdir, "rotation_gain_vs_duration.png") + fig.savefig(png, dpi=140) + print("wrote", png) + except Exception as e: + print("(plot skipped:", e, ")") + + +if __name__ == "__main__": + main() diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/finite_size_gain_vs_arm.png b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/finite_size_gain_vs_arm.png new file mode 100644 index 000000000..6ab85ae0f Binary files /dev/null and b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/finite_size_gain_vs_arm.png differ diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/finite_size_gain_vs_arm.txt b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/finite_size_gain_vs_arm.txt new file mode 100644 index 000000000..9719fe7fd --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/finite_size_gain_vs_arm.txt @@ -0,0 +1,6 @@ +# finite-size (Path D) demo: --freqresponse recovers SNR the long-wavelength analysis loses +# label L_arm[m] fL/c@fmax SNR lnL_LWL lnL_finite gain=lnL_finite-lnL_LWL +LIGO_4km 4000 0.0266851 318.345 50657 50657.3 0.250481 +ET_10km 10000 0.0667128 319.525 51029.7 51033.1 3.47031 +CE_20km 20000 0.133426 321.377 51613.5 51625.9 12.4271 +CE_40km 40000 0.266851 324.621 52633.5 52673 39.5546 diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/rotation_gain_vs_duration.png b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/rotation_gain_vs_duration.png new file mode 100644 index 000000000..532150b42 Binary files /dev/null and b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/rotation_gain_vs_duration.png differ diff --git a/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/rotation_gain_vs_duration.txt b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/rotation_gain_vs_duration.txt new file mode 100644 index 000000000..14bb8b27f --- /dev/null +++ b/MonteCarloMarginalizeCode/Code/demo/rift/slowrot/outputs/rotation_gain_vs_duration.txt @@ -0,0 +1,5 @@ +# slow-rotation demo: --rotation-slow recovers SNR the static analysis loses on long signals +# label seglen[s] SNR Omega*T lnL_static lnL_rot gain=lnL_rot-lnL_static +null_bbh 2 30.0567 0.000145842 70.1014 70.105 0.00367852 +bbh_8_8 16 29.1207 0.00116674 106.254 106.297 0.0430423 +bbh_4_4 64 31.1072 0.00466695 107.632 107.819 0.186851