序列 PCA 定秩

#!/usr/bin/env python3
import argparse
import numpy as np
from larch.io import read_athena
from larch.xafs import pre_edge, find_e0
from larch.math import interp

def flat_window(g, emin, emax):
    e0 = find_e0(g.energy, g.mu)
    pre_edge(g, e0=e0, pre1=-150, pre2=-30, norm1=150, norm2=800, nvict=2)
    e, y = np.asarray(g.energy), np.asarray(g.flat)
    m = (e >= e0 + emin) & (e <= e0 + emax)
    return e[m] - e0, y[m]

def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("prj")
    ap.add_argument("--emin", type=float, default=-20)
    ap.add_argument("--emax", type=float, default=60)
    args = ap.parse_args()

    groups = list(read_athena(args.prj)._athena_groups.items())
    rel0, y0 = flat_window(groups[0][1], args.emin, args.emax)
    M = [y0]
    for _, g in groups[1:]:
        rel, y = flat_window(g, args.emin, args.emax)
        M.append(interp(rel, y, rel0))
    A = np.vstack(M)
    _, s, _ = np.linalg.svd(A - A.mean(0), full_matrices=False)
    tot = s.sum() or 1
    print(f"{'comp':>6} {'singular':>12} {'frac':>8} {'cum':>8}")
    cum = 0.0
    for i, si in enumerate(s, 1):
        cum += si / tot
        print(f"{i:6d} {si:12.4f} {si/tot:8.3f} {cum:8.3f}")

if __name__ == "__main__":
    main()