import numpy as np

import matplotlib.pyplot as plt

from gwpy.timeseries import TimeSeries

from pycbc.waveform import get_fd_waveform

from pycbc.filter import match, sigmasq

from scipy.optimize import minimize

# Step 1: Fetch strain data from GWOSC (run this locally with internet)

event = 'GW150914'

detector = 'H1' # Hanford; repeat for L1

url = f'https://gwosc.org/archive/links/{event}/R1/16KHZ/{detector}-{event}_R1-1126259446-32.gwf'

data = TimeSeries.read(url, start=1126259462.4 - 2, end=1126259462.4 + 2) # 4s window

data = data.whiten() # Whiten to normalize noise

# Step 2: Generate GR template (binary black hole example)

hp, hc = get_fd_waveform(approximant='IMRPhenomPv2', mass1=36, mass2=29, spin1z=0, spin2z=0,

delta_f=data.delta_f, f_lower=20)

# Step 3: Define modulation factor (Cycle Function projection)

def modulation(f, alpha, f_i): # For one i; product over multiple

return 1 + alpha * np.abs(np.sin(2 * np.pi * f / f_i))**2

# Custom template: h_custom(f) = h_GR(f) * prod(modulation)

def custom_template(params): # params = [alpha1, f1, alpha2, f2, ...]

mod = np.ones_like(hp)

for i in range(0, len(params), 2):

mod *= modulation(hp.sample_frequencies, params[i], params[i+1])

return hp * mod, hc * mod # Assuming same mod for plus/cross

# Step 4: Optimize custom template fit (maximize match with data)

def objective(params):

hp_custom, _ = custom_template(params)

m, _ = match(hp_custom, data.to_pycbc(), low_frequency_cutoff=20)

return -m # Minimize negative match

initial_guess = [0.01, 100, 0.01, 200] # e.g., 2 dimensions for simplicity

result = minimize(objective, initial_guess, method='Nelder-Mead')

best_params = result.x

hp_best, _ = custom_template(best_params)

# Step 5: Subtract best-fit and check residuals

residual = data.to_pycbc() - hp_best

residual_ts = residual.to_timeseries() # Time domain

# Plot residuals

plt.figure(figsize=(10, 4))

plt.plot(residual_ts.times, residual_ts, label='Residuals')

plt.xlabel('Time (s)')

plt.ylabel('Whitened Strain')

plt.title('Residuals After Custom Template Subtraction')

plt.legend()

plt.show()

# Statistical tests (e.g., Kolmogorov-Smirnov for Gaussianity)

from scipy.stats import kstest, anderson

ks_stat, ks_p = kstest(residual_ts.value, 'norm') # Test vs standard normal

ad_stat = anderson(residual_ts.value, dist='norm')

print(f'KS Test: stat={ks_stat:.3f}, p={ks_p:.3f} (p>0.05 consistent with noise)')

print(f'Anderson-Darling: stat={ad_stat.statistic:.3f} (compare to critical values {ad_stat.critical_values})')

# If p < 0.05 or AD stat > critical, possible structure/deviation