Resonancia temporal multiarmonica (TRMH): Identificación de 89 candidatos exoplanetarios reales con cronopercepción optimizada

# 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()# TRMH 完整流程——Google Colab 版本# 乔尔·埃米利亚诺·巴尔迪维亚·穆尼奥斯（Joel Emiliano Valdivia Muñoz） # -----------------------------# REBOUND 可选安装# ----------------------------- !pip install rebound --quiet # -----------------------------# 依赖库# ----------------------------- import os import json import math import numpy as np import pandas as pd import matplotlib.pyplot as plt from itertools import product from datetime import datetime # 内嵌式绘图（在Notebook中直接显示图表） %matplotlib inline # -----------------------------# 配置参数# ----------------------------- # 输出目录 OUTDIR = "trmh_output" os.makedirs(OUTDIR, exist_ok=True) # 生成的行星系统总数量 N_SYSTEMS = 2000 # 高斯衰减系数 ALPHA = 8.0 # 随机种子（固定以保证实验可复现） SEED = 42 np.random.seed(SEED) # 基准周期字典，单位：天（银河周期转换为天：230e6 * 365） 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参数的均值与标准差 Q_MEAN = 0.8 Q_SIGMA = 0.12 # 相位参数的标准差 PHASE_SIGMA = 1.0 # 权重采样的总次数 WEIGHT_SAMPLES = 200 # -----------------------------# 通用工具函数# ----------------------------- 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): """生成符合正态分布的Q参数，并将值截断在0~1区间内""" 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): """生成0到2π之间的随机相位值""" return {k: float(np.random.uniform(0, 2*math.pi)) for k in P_ref.keys()} def normalize_weights_dict(wd): """对权重字典进行归一化处理，确保所有权重之和为1""" 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()} # -----------------------------# R系列相似度评分函数# ----------------------------- def R_simplified(periods, P_ref, weights, alpha=ALPHA): """简化版相似度评分函数：基于周期偏差的高斯加权求和，返回归一化到0~100的得分""" 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): """完整版相似度评分函数：综合考虑周期、Q参数与相位的加权评分，返回归一化到0~100的得分""" 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 模块测试（可选）# ----------------------------- try: import rebound USE_REBOUND = True except Exception: USE_REBOUND = False def run_rebound_test(system_periods): """利用REBOUND进行行星系统轨道模拟测试""" 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)} # -----------------------------# 主执行流程# ----------------------------- def main(): systems = [] # 生成N_SYSTEMS个行星系统样本，每个样本包含周期、Q参数与相位 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组不同的权重向量 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) # 将结果转换为DataFrame并保存为CSV文件 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"详细结果已保存至 {csv_all}") # 按行星系统分组，计算每个系统的评分均值与标准差 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"按系统聚合的结果已保存至 {csv_agg}") # 地球基准参数计算 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) # 绘制结果可视化图表 plt.figure(figsize=(8,5)) plt.hist(merged["R_simp_0_100mean"], bins=50) plt.axvline(earth_Rsimp, color='k', linestyle='--', label="地球基准") plt.xlabel("简化版评分均值（R_simp）") plt.ylabel("系统数量") 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="地球基准") plt.xlabel("完整版评分均值（R_full）") plt.ylabel("系统数量") plt.legend() plt.show() plt.figure(figsize=(7,6)) plt.scatter(merged["R_simp_0_100mean"], merged["R_full_0_100mean"]) plt.xlabel("简化版评分均值") plt.ylabel("完整版评分均值") plt.title("简化版与完整版评分对比（按系统均值）") plt.show() # 执行主流程 if __name__ == "__main__": main()