示例03 纯 NumPy 手算 DFA

"""ex03_manual_dfa.py — 第 3 章 用纯 NumPy 手算 DFA，对照 fathon

目标：通过自己实现单窗口的 DFA 波动，理解 F(n) 到底在算什么，
然后与 fathon 的 C 加速结果对比，验证一致性。

对应教程：tutorials/03-dfa-principles.md 3.8、3.9 节。

运行：
python ex03_manual_dfa.py
"""

import numpy as np
import fathon
from fathon import fathonUtils as fu

def manual_F(profile, n, polOrd=1):
"""单尺度 n 的 DFA 波动（简化版，仅正向分段）。

DFA 五步：
  1. profile 已构造好（外部完成）
  2. 按 n 分段（不重叠）
  3. 每段拟合 polOrd 阶多项式，得残差
  4. 段内方差取平均
  5. （这里只算一个 n，多 n 的 log-log 拟合在外面做）
"""
L = len(profile)
N_n = L // n                      # 能分成多少段
F2_list = []
for nu in range(N_n):
    seg = profile[nu * n:(nu + 1) * n]
    t = np.arange(n, dtype=float)
    coeffs = np.polyfit(t, seg, polOrd)   # 最小二乘多项式拟合
    trend = np.polyval(coeffs, t)
    residual = seg - trend
    F2_list.append(np.mean(residual ** 2))
return np.sqrt(np.mean(F2_list))
​

def main():
np.random.seed(0)
L = 4096
profile = fu.toAggregated(np.random.randn(L))

win_sizes = np.array([16, 64, 256, 1024])

# ── 手算 F(n) ──
print("=== 手算 DFA ===")
F_manual = np.array([manual_F(profile, n) for n in win_sizes])
for n, F in zip(win_sizes, F_manual):
    print(f"  n={n:5d}  F(n)={F:.4f}")
H_manual = np.polyfit(np.log(win_sizes), np.log(F_manual), 1)[0]
print(f"  手算斜率 H = {H_manual:.4f}")

# ── fathon 验证 ──
print("\n=== fathon DFA ===")
dfa = fathon.DFA(profile)
n, F = dfa.computeFlucVec(win_sizes, revSeg=False, polOrd=1)
H, _ = dfa.fitFlucVec()
print("  窗口 n:", n)
print("  F(n):  ", np.round(F, 4))
print(f"  fathon H = {H:.4f}")

# ── 数值一致性检查 ──
diff = np.max(np.abs(F_manual - F))
print(f"\n|F_manual - F_fathon| max = {diff:.4e}")
print(f"|H_manual - H_fathon|     = {abs(H_manual - H):.4e}")
print("(差异源于 fathon 默认 revSeg=False 但内部实现细节略有不同；"
      "整体数量级应一致，斜率应高度吻合。)")
​

if name == "main":
main()