"""Inject a point-lensed leading-order chirp into synthetic Gaussian detector noise. Everything here is synthetic: the noise is drawn from the Ajith (2011) Eq. (4.7) fit to the Advanced LIGO zero-detuned high-power design PSD, and the source is the Newtonian stationary-phase chirp (Cutler & Flanagan 1994, Eqs. 41-42, leading term only). No GWOSC strain is read. Lensed templates use the analytic point-lens factor from point_lens_reference.py, tabulated on a fixed w grid and interpolated. For a template h with amplitude, arrival time and phase free, the maximized log likelihood ratio against pure noise is ln L = max_t |(t)|^2 / 2, where (t) = 4 df sum conj(a) b / S exp(-2 pi i f t). The lensing statistic is Delta ln L = ln L(best lensed template) - ln L(best unlensed template). Lensed templates keep the true component mass and search lens mass and y: a coarse grid, then a finer mass scan around the best coarse mass for each y, then a bounded polish. Unlensed templates search the component mass only. """ import json import numpy as np import mpmath as mp from point_lens_reference import analytic_amplification mp.mp.dps = 30 SOLAR_MASS_SECONDS = 1.3271244e20 / 299792458.0**3 SAMPLE_RATE, DURATION, COALESCENCE_TIME, LOW_FREQUENCY = 1024, 8.0, 6.6, 20.0 SAMPLES = int(SAMPLE_RATE * DURATION) FREQUENCY_STEP = 1 / DURATION FREQUENCIES = np.arange(SAMPLES) * FREQUENCY_STEP def design_psd(f): x = f / 245.4 return 1e-48 * (0.0152 * x**-4 + 0.2935 * x**2.25 + 2.7951 * x**1.5 - 6.5080 * x**0.75 + 17.7622) def chirp(component_mass): total = 2 * component_mass * SOLAR_MASS_SECONDS chirp_mass = (component_mass**2) ** 0.6 / (2 * component_mass) ** 0.2 * SOLAR_MASS_SECONDS upper = 1 / (6**1.5 * np.pi * total) f = FREQUENCIES band = (f >= LOW_FREQUENCY) & (f <= upper) rise = np.clip((f - LOW_FREQUENCY) / 2, 0, 1) fall = np.clip((upper - f) / (0.05 * upper), 0, 1) taper = np.sin(np.pi / 2 * rise) ** 2 * np.sin(np.pi / 2 * fall) ** 2 h = np.zeros(SAMPLES, dtype=complex) fb = f[band] phase = 2 * np.pi * fb * COALESCENCE_TIME - np.pi / 4 + 3 / 128 * (np.pi * chirp_mass * fb) ** (-5 / 3) h[band] = taper[band] * fb ** (-7 / 6) * np.exp(1j * phase) return h, band def inner(a, b, psd, band): return 4 * FREQUENCY_STEP * np.real(np.sum(np.conj(a[band]) * b[band] / psd[band])) def max_overlap(data, template, psd, band): weighted = np.zeros(SAMPLES, dtype=complex) weighted[band] = np.conj(template[band]) * data[band] / psd[band] return 4 * FREQUENCY_STEP * np.max(np.abs(np.fft.fft(weighted))) class AmplificationTable: """F(w, y) on a log w grid per y; linear interpolation of real and imaginary parts. 6000 points keep about 14 samples per fringe at the longest template delay (y = 3, w = 41). Measured worst relative interpolation error for y = 2 and 3, 5 <= w <= 41: 0.24 percent (1500 points gave 3.8 percent; see README). """ def __init__(self, y_values, w_min=0.01, w_max=50.0, points=6000): self.w = np.geomspace(w_min, w_max, points) self.values = {round(y, 6): np.array([analytic_amplification(w, y) for w in self.w]) for y in y_values} def __call__(self, w, y): table = self.values[round(y, 6)] return np.interp(w, self.w, table.real) + 1j * np.interp(w, self.w, table.imag) def lensed(h, band, lens_mass, y, table): out = h.copy() w = 8 * np.pi * SOLAR_MASS_SECONDS * lens_mass * FREQUENCIES[band] out[band] = h[band] * table(w, y) return out def noise(psd, band, seed): rng = np.random.default_rng(seed) sigma = np.sqrt(psd / (4 * FREQUENCY_STEP)) n = np.zeros(SAMPLES, dtype=complex) n[band] = sigma[band] * (rng.standard_normal(band.sum()) + 1j * rng.standard_normal(band.sum())) return n def best_log_likelihood(data, templates, psd, band): return max(0.5 * max_overlap(data, t, psd, band) ** 2 / inner(t, t, psd, band) for t in templates) def lensed_score(data, h0, band, psd, table, lens_mass, y): template = lensed(h0, band, lens_mass, y, table) return 0.5 * max_overlap(data, template, psd, band) ** 2 / inner(template, template, psd, band) def refined_lensed_best(data, h0, band, psd, table, mass_grid, y_grid, coarse_scores, refine_points): """Refine the lens mass between the neighbours of the best coarse mass, for every y. The coarse grid steps the mass by about 15 percent. At y = 2 the delay is 90 ms, so a 15 percent step moves the fringe phase at 200 Hz by several radians and a coarse-only search misses the signal. The scan below steps by about 0.6 percent, then a bounded scalar search polishes the best scan point. """ from scipy.optimize import minimize_scalar candidates = [] for column, y in enumerate(y_grid): index = int(np.argmax(coarse_scores[:, column])) low, high = mass_grid[max(index - 1, 0)], mass_grid[min(index + 1, len(mass_grid) - 1)] scan = np.geomspace(low, high, refine_points) values = [lensed_score(data, h0, band, psd, table, m, y) for m in scan] best = int(np.argmax(values)) candidates.append((values[best], scan[best], y, scan[max(best - 1, 0)], scan[min(best + 1, refine_points - 1)])) candidates.sort(reverse=True) polished = [] for value, mass, y, low, high in candidates[:2]: result = minimize_scalar(lambda log_m: -lensed_score(data, h0, band, psd, table, np.exp(log_m), y), bounds=(np.log(low), np.log(high)), method="bounded", options={"xatol": 1e-5}) polished.append(max((value, mass, y), (-result.fun, float(np.exp(result.x)), y))) return max(polished) def run(inputs): psd = np.ones(SAMPLES) source_mass = inputs["source_component_mass_solar"] h0, band = chirp(source_mass) psd[band] = design_psd(FREQUENCIES[band]) distance_scale = inputs["unlensed_snr"] / np.sqrt(inner(h0, h0, psd, band)) y_grid = inputs["template_grid"]["impact_parameters"] mass_grid = np.geomspace(*inputs["template_grid"]["lens_mass_range_solar"], inputs["template_grid"]["lens_mass_count"]) table = AmplificationTable(sorted(set(y_grid + [c["impact_parameter"] for c in inputs["cases"] if c["impact_parameter"] is not None]))) unlensed_bank = [chirp(m)[0] for m in np.linspace(*inputs["template_grid"]["unlensed_component_mass_range_solar"], inputs["template_grid"]["unlensed_component_mass_count"])] lensed_bank = [(m, y, lensed(h0, band, m, y, table)) for m in mass_grid for y in y_grid] refine_points = inputs["template_grid"]["lens_mass_refinement_points"] results = [] for case in inputs["cases"]: if case["impact_parameter"] is None: continue # Noise-free check: the search must recover the injected lens and the exact statistic. signal = lensed(h0, band, case["lens_mass_solar"], case["impact_parameter"], table) * distance_scale unlensed_best = best_log_likelihood(signal, unlensed_bank, psd, band) coarse = np.array([0.5 * max_overlap(signal, t, psd, band) ** 2 / inner(t, t, psd, band) for _, _, t in lensed_bank]).reshape(len(mass_grid), len(y_grid)) found, found_mass, found_y = refined_lensed_best(signal, h0, band, psd, table, mass_grid, y_grid, coarse, refine_points) results.append({ "id": case["id"] + "-noise-free", "label": case["label"] + ", no noise", "exact_template_delta_log_likelihood": float(0.5 * inner(signal, signal, psd, band) - unlensed_best), "coarse_grid_only_delta_log_likelihood": float(coarse.max() - unlensed_best), "search_delta_log_likelihood": float(found - unlensed_best), "search_lens_mass_solar": float(found_mass), "search_impact_parameter": found_y, }) for case in inputs["cases"]: if case["impact_parameter"] is None: signal = h0 * distance_scale else: signal = lensed(h0, band, case["lens_mass_solar"], case["impact_parameter"], table) * distance_scale rows = [] for seed in range(case["first_seed"], case["first_seed"] + case["realizations"]): data = signal + noise(psd, band, seed) unlensed_best = best_log_likelihood(data, unlensed_bank, psd, band) coarse = np.array([0.5 * max_overlap(data, t, psd, band) ** 2 / inner(t, t, psd, band) for _, _, t in lensed_bank]).reshape(len(mass_grid), len(y_grid)) lensed_best, best_mass, best_y = refined_lensed_best(data, h0, band, psd, table, mass_grid, y_grid, coarse, refine_points) rows.append({"seed": seed, "delta_log_likelihood": lensed_best - unlensed_best, "best_lens_mass_solar": best_mass, "best_impact_parameter": best_y}) deltas = np.array([r["delta_log_likelihood"] for r in rows]) results.append({ "id": case["id"], "label": case["label"], "realizations": len(rows), "optimal_snr_of_injection": float(np.sqrt(inner(signal, signal, psd, band))), "delta_log_likelihood_median": float(np.median(deltas)), "delta_log_likelihood_90th_percentile": float(np.percentile(deltas, 90)), "delta_log_likelihood_maximum": float(deltas.max()), "fraction_above_threshold": float(np.mean(deltas > inputs["reporting_threshold_delta_log_likelihood"])), "realization_with_largest_delta": max(rows, key=lambda r: r["delta_log_likelihood"]), }) return results if __name__ == "__main__": import sys inputs = json.load(open(sys.argv[1] if len(sys.argv) > 1 else "injection-inputs.json")) print(json.dumps(run(inputs), indent=2))