"""
================================================================================
 CONTROL DE HIPÓTESIS NULA – FASE ALEATORIA (MÉTODO DE 8 CANALES VIRTUALES)
 ================================================================================
 Descarga un segmento oficial de LIGO (enero 2024), aplica el protocolo de
 8 canales virtuales y mide Γ. Luego repite el cálculo tras multiplicar el strain
 por una fase aleatoria. Si Γ se mantiene alto, el método es insensible a la
 fase real y, por tanto, espurio.
 ================================================================================
"""

import numpy as np
import matplotlib.pyplot as plt
import gc, sys, warnings
warnings.filterwarnings('ignore')

# 1. Dependencias
try:
    from gwpy.timeseries import TimeSeries
except ImportError:
    import subprocess, sys
    subprocess.check_call([sys.executable, "-m", "pip", "install", "gwpy"])
    from gwpy.timeseries import TimeSeries

# 2. Parámetros fijos (los usados en aquel análisis)
TAU_BASE = 0.3697
K_EARLY = 0.967

def calcular_gamma(matriz_8, central, t, tau):
    """
    Protocolo 33 original: suma constructiva de 8 canales con compensación de fase.
    """
    phi = 2 * np.pi * tau * np.log(t + 1e-9)
    torsion_8 = [matriz_8[k] * np.exp(-1j*(phi + (k*np.pi/4))) for k in range(8)]
    suma_8 = np.sum(torsion_8, axis=0)
    central_demod = central * np.exp(-1j*phi)
    if np.std(suma_8)==0 or np.std(central_demod)==0:
        return 0.0
    corr = np.corrcoef(np.real(suma_8), np.real(central_demod))[0,1]
    return np.abs(corr) * K_EARLY * np.log(36)

# 3. Descarga del segmento de referencia (enero 2024, GPS 1389379584)
print("Descargando segmento oficial (GPS 1389379584, 32 s)...")
data = TimeSeries.fetch_open_data('H1', 1389379584, 1389379584+32, verbose=False)
fs = data.sample_rate.value
strain = data.whiten().bandpass(20, 500).value.astype(np.float64)
t = np.arange(len(strain)) / fs

# 4. Construcción de 8 canales virtuales (el método original)
central = strain
mat = np.array([central * np.cos(k*np.pi/4) + np.random.randn(len(strain))*1e-4 for k in range(8)])

# 5. Medición de Γ sobre el strain real
tau_max_real, gamma_max_real = 0, 0
for tau_centro in [TAU_BASE]:
    taus = np.linspace(tau_centro*0.8, tau_centro*1.2, 50)
    gammas = [calcular_gamma(mat, central, t, tau) for tau in taus]
    idx = np.argmax(gammas)
    tau_max_real, gamma_max_real = taus[idx], gammas[idx]
print(f"\nΓ máximo sobre strain real = {gamma_max_real:.4f} (en τ≈{tau_max_real:.4f})")

# 6. Prueba de fase aleatoria: destruimos cualquier estructura de fase
rng = np.random.default_rng(42)
fase_lineal_aleatoria = 2 * np.pi * rng.uniform(0, 1, size=len(t)) * np.arange(len(t)) / len(t)
strain_modulado = central * np.exp(1j * fase_lineal_aleatoria)  # modulación compleja
central_mod = np.real(strain_modulado)
mat_mod = np.array([central_mod * np.cos(k*np.pi/4) + np.random.randn(len(central_mod))*1e-4 for k in range(8)])

# Medición de Γ sobre la señal con fase aleatoria
tau_max_mod, gamma_max_mod = 0, 0
for tau_centro in [TAU_BASE]:
    taus_mod = np.linspace(tau_centro*0.8, tau_centro*1.2, 50)
    gammas_mod = [calcular_gamma(mat_mod, central_mod, t, tau) for tau in taus_mod]
    idx_mod = np.argmax(gammas_mod)
    tau_max_mod, gamma_max_mod = taus_mod[idx_mod], gammas_mod[idx_mod]
print(f"Γ máximo tras fase aleatoria = {gamma_max_mod:.4f} (en τ≈{tau_max_mod:.4f})")

# 7. Interpretación
print("\n" + "="*60)
if abs(gamma_max_real - gamma_max_mod) < 0.1:
    print(" CONCLUSIÓN: Γ es prácticamente idéntico en ambos casos.")
    print(" La fase logarítmica real no influye → el método genera un artefacto.")
else:
    print(" Γ cambió significativamente; el método podría tener sensibilidad real.")
print("="*60)