# ============================================
#  Excel → スペクトル読み込み
#  ピーク検出 → ローカル線形正規化
#  persistent homology → 可視化 → Excel保存
# ============================================

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.signal import find_peaks
from ripser import ripser
from persim import plot_diagrams

# -----------------------------
# 1. Excel からスペクトル読み込み
# -----------------------------
# Excel 形式：1列目 x、2列目 y と仮定
df = pd.read_excel("spectrum.xlsx")
x = df.iloc[:, 0].values
y = df.iloc[:, 1].values

# -----------------------------
# 2. ピーク検出（強いピークをアンカーにする）
# -----------------------------
# prominence や height は適宜調整
peaks, props = find_peaks(y, prominence=0.1, height=0.1)

# 強いピークだけを選ぶ（上位 N 個など）
N = 5
if len(peaks) > N:
    strong_idx = np.argsort(props["peak_heights"])[-N:]
    strong_peaks = peaks[strong_idx]
else:
    strong_peaks = peaks

# -----------------------------
# 3. ローカル線形正規化のスケール関数を作る
# -----------------------------
# 強いピークの (x, y) を取り出す
xp = x[strong_peaks]
yp = y[strong_peaks]

# x の範囲をカバーするために端点も追加
xp_full = np.concatenate([[x.min()], xp, [x.max()]])
yp_full = np.concatenate([[y.min()], yp, [y.min()]])

# 線形補間でスケール関数 s(x) を作る
s = np.interp(x, xp_full, yp_full)

# 正規化スペクトル
y_norm = y / s

# -----------------------------
# 4. persistent homology の計算
# -----------------------------
# 1次元関数 → (x, y_norm) の点群として扱う
# ただし x は距離スケールに影響するので正規化
X = np.vstack([x / (x.max() - x.min()), y_norm]).T

result = ripser(X, maxdim=1)
dgms = result["dgms"]

# -----------------------------
# 5. 可視化（元スペクトル・正規化スペクトル・PH）
# -----------------------------
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 元スペクトル
axes[0, 0].plot(x, y)
axes[0, 0].scatter(x[strong_peaks], y[strong_peaks], color="red")
axes[0, 0].set_title("Original Spectrum with Strong Peaks")

# 正規化スペクトル
axes[0, 1].plot(x, y_norm)
axes[0, 1].set_title("Locally Normalized Spectrum")

# persistence diagram
plot_diagrams(dgms, ax=axes[1, 0])
axes[1, 0].set_title("Persistence Diagram (H0, H1)")

# H1 のみ
plot_diagrams([dgms[1]], ax=axes[1, 1])
axes[1, 1].set_title("H1 Only")

plt.tight_layout()
plt.show()

# -----------------------------
# 6. Excel に保存（正規化スペクトルと PH 結果）
# -----------------------------
# PH の birth-death を DataFrame にまとめる
df_h0 = pd.DataFrame(dgms[0], columns=["birth", "death"])
df_h1 = pd.DataFrame(dgms[1], columns=["birth", "death"])

with pd.ExcelWriter("spectrum_PH_output.xlsx") as writer:
    pd.DataFrame({"x": x, "y": y, "y_norm": y_norm}).to_excel(writer, sheet_name="spectrum", index=False)
    df_h0.to_excel(writer, sheet_name="H0", index=False)
    df_h1.to_excel(writer, sheet_name="H1", index=False)

print("Excel に保存しました： spectrum_PH_output.xlsx")