TRMH

# 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()