DOI: 10.5281/zenodo.17619811
title: TRMH

description: # TRMH Full Pipeline - Google Colab Version# Joel Emiliano Valdivia Muñoz

# -----------------------------# Instalación opcional de REBOUND# -----------------------------!pip install rebound --quiet

# -----------------------------# Librerías# -----------------------------import osimport jsonimport mathimport numpy as npimport pandas as pdimport matplotlib.pyplot as pltfrom itertools import productfrom datetime import datetime

# Mostrar gráficos inline%matplotlib inline

# -----------------------------# Configuración# -----------------------------OUTDIR = "trmh_output"os.makedirs(OUTDIR, exist_ok=True)

N_SYSTEMS = 2000ALPHA = 8.0SEED = 42np.random.seed(SEED)

CYCLES = {    "rotacion_planetaria": 1.0,    "lunacion": 29.53,    "traslacion_planetaria": 365.25,    "rotacion_estelar": 25.0,    "ciclo_magnetico_estelar": 4018,    "periodo_galactico": 230e6 * 365}

DEFAULT_WEIGHTS = {    "rotacion_planetaria": 0.20,    "lunacion": 0.10,    "traslacion_planetaria": 0.30,    "rotacion_estelar": 0.15,    "ciclo_magnetico_estelar": 0.18,    "periodo_galactico": 0.07}

Q_MEAN = 0.8Q_SIGMA = 0.12PHASE_SIGMA = 1.0WEIGHT_SAMPLES = 200

# -----------------------------# Funciones utilitarias# -----------------------------def sample_periods(P_ref):    out = {}    for k, v in P_ref.items():        if k == "rotacion_planetaria": sigma = 0.30        elif k == "lunacion": sigma = 0.40        elif k == "traslacion_planetaria": sigma = 0.25        elif k == "rotacion_estelar": sigma = 0.35        elif k == "ciclo_magnetico_estelar": sigma = 0.6        else: sigma = 0.5        factor = np.random.lognormal(mean=0.0, sigma=sigma)        out[k] = v * factor    return out

def sample_Qs(P_ref):    qs = {}    for k in P_ref.keys():        q = np.random.normal(Q_MEAN, Q_SIGMA)        qs[k] = float(np.clip(q, 0.0, 1.0))    return qs

def sample_phases(P_ref):    return {k: float(np.random.uniform(0, 2*math.pi)) for k in P_ref.keys()}

def normalize_weights_dict(wd):    s = sum(wd.values())    if s == 0:        return {k: 1.0/len(wd) for k in wd.keys()}    return {k: float(v)/s for k, v in wd.items()}

# -----------------------------# Funciones de R# -----------------------------def R_simplified(periods, P_ref, weights, alpha=ALPHA):    raw_terms = [weights[k]*math.exp(-alpha*((periods[k]-P_ref[k])/P_ref[k])**2) for k in P_ref.keys()]    raw = float(np.sum(raw_terms))    max_possible = float(np.sum(list(weights.values())))    R100 = 100.0 * raw / max_possible if max_possible>0 else 0.0    return R100, raw, raw_terms

def R_full(periods, freqs_ref, weights, Qs, phases, alpha=ALPHA, beta=1.0):    terms = []    for k in freqs_ref.keys():        P = periods[k]        f = 1.0/P        f_ref = 1.0/freqs_ref[k]        delta_f2 = (f-f_ref)**2        q = Qs.get(k, 1.0)        w = weights.get(k, 0.0)        phi = phases.get(k, 0.0)        delta_phi = phi        term = w*q*math.exp(-alpha*(delta_f2/(f_ref**2))) * math.cos(beta*delta_phi)        terms.append(term)    raw = float(np.sum(terms))    max_possible = float(np.sum([abs(weights[k]*1.0) for k in freqs_ref.keys()]))    R100 = 100.0 * raw / max_possible if max_possible>0 else 0.0    return R100, raw, terms

def sample_weight_vector(base_weights, scale=1.0):    keys = list(base_weights.keys())    base = np.array([base_weights[k] for k in keys], dtype=float)    noise = np.random.lognormal(0.0, 0.5*scale, size=len(keys))    w = base*noise    w /= np.sum(w)    return {k: float(wi) for k, wi in zip(keys, w)}

# -----------------------------# REBOUND (opcional)# -----------------------------try:    import rebound    USE_REBOUND = Trueexcept Exception:    USE_REBOUND = False

def run_rebound_test(system_periods):    if not USE_REBOUND: return {"rebound_available": False}    sim = rebound.Simulation()    sim.units = ('yr','AU','Msun')    sim.add(m=1.0)    P_days = system_periods.get("traslacion_planetaria", 365.25)    P_yr = P_days/365.25    a = P_yr**(2/3)    sim.add(m=3e-6, a=a)    sim.move_to_com()    try:        sim.integrate(1.0)        return {"rebound_available": True, "smoke_test": True, "a": a}    except Exception as e:        return {"rebound_available": True, "smoke_test": False, "error": str(e)}

# -----------------------------# Main pipeline# -----------------------------def main():    systems = []    for i in range(N_SYSTEMS):        P = sample_periods(CYCLES)        Qs = sample_Qs(CYCLES)        phases = sample_phases(CYCLES)        systems.append({"id": f"sys_{i:04d}", "periods": P, "Qs": Qs, "phases": phases})

    weight_samples = [sample_weight_vector(DEFAULT_WEIGHTS) for _ in range(WEIGHT_SAMPLES)]    results_rows = []

    for widx, wvec in enumerate(weight_samples):        wnorm = normalize_weights_dict(wvec)        for sys in systems:            periods = sys["periods"]            Qs = sys["Qs"]            phases = sys["phases"]            Rs_simp, raw_simp, terms_simp = R_simplified(periods, CYCLES, wnorm)            Rs_full, raw_full, terms_full = R_full(periods, CYCLES, wnorm, Qs, phases)            row = {"id": sys["id"], "weight_sample_idx": widx,                   "R_simp_0_100": Rs_simp, "R_full_0_100": Rs_full,                   "raw_simp": raw_simp, "raw_full": raw_full}            for k,v in periods.items(): row[f"P_{k}"]=v            for k,v in Qs.items(): row[f"Q_{k}"]=v            for k,v in phases.items(): row[f"phi_{k}"]=v            results_rows.append(row)

    df = pd.DataFrame(results_rows)    csv_all = os.path.join(OUTDIR, f"trmh_all_{N_SYSTEMS}systems.csv")    df.to_csv(csv_all, index=False)    print(f"Saved detailed results to {csv_all}")

    # Agregado por sistema    agg = df.groupby("id").agg({"R_simp_0_100":["mean","std"],"R_full_0_100":["mean","std"]})    agg.columns = ["".join(col).strip() for col in agg.columns.values]    agg = agg.reset_index()    rep_periods = {s["id"]: s["periods"] for s in systems}    df_rep = pd.DataFrame([{"id":sid, **rep_periods[sid]} for sid in rep_periods.keys()])    merged = agg.merge(df_rep, on="id", how="left")    csv_agg = os.path.join(OUTDIR, f"trmh_agg_{N_SYSTEMS}systems.csv")    merged.to_csv(csv_agg, index=False)    print(f"Saved aggregated per-system results to {csv_agg}")

    # Referencia Tierra    earth_periods = CYCLES.copy()    earth_Qs = {k: Q_MEAN for k in CYCLES.keys()}    earth_phases = {k:0.0 for k in CYCLES.keys()}    earth_Rsimp, _, _ = R_simplified(earth_periods, CYCLES, normalize_weights_dict(DEFAULT_WEIGHTS))    earth_Rfull, _, _ = R_full(earth_periods, CYCLES, normalize_weights_dict(DEFAULT_WEIGHTS), earth_Qs, earth_phases)

    # Gráficos    plt.figure(figsize=(8,5))    plt.hist(merged["R_simp_0_100mean"], bins=50)    plt.axvline(earth_Rsimp, color='k', linestyle='--', label="Earth")    plt.xlabel("R_simp (mean)")    plt.ylabel("N systems")    plt.legend()    plt.show()

    plt.figure(figsize=(8,5))    plt.hist(merged["R_full_0_100mean"], bins=50)    plt.axvline(earth_Rfull, color='k', linestyle='--', label="Earth")    plt.xlabel("R_full (mean)")    plt.ylabel("N systems")    plt.legend()    plt.show()

    plt.figure(figsize=(7,6))    plt.scatter(merged["R_simp_0_100mean"], merged["R_full_0_100mean"])    plt.xlabel("R_simp mean")    plt.ylabel("R_full mean")    plt.title("R_simp vs R_full (media por sistema)")    plt.show()

# Ejecutar pipelineif __name__ == "__main__":    main()