init: 验证robust算法

core/robust.py

@@ -0,0 +1,191 @@
+from models.basic import ModelBasic
+from typing import List,Literal,Dict
+import numpy as np
+import random
+from pydantic import BaseModel
+from skrf import Network
+from models.basic import ModelBasicDatasetUnit
+
+class RobustParametricConfig(BaseModel):
+    n_poles: int
+    max_iter: int = 10
+    parameter_type: Literal["s","y","z"]
+
+class RobustParametricModel:
+    def __init__(self,model:ModelBasic,config:RobustParametricConfig):
+        self.model = model
+        self.config = config
+
+    # 区分训练集和测试集
+    def _train_test_split(self,train_ratio:float=0.8,random_state:int=42):
+        random.seed(random_state)
+        dataset = self.model.results
+        random.shuffle(dataset)
+        train_size = int(len(dataset)*train_ratio)
+        self.train_set = dataset[:train_size]
+        self.test_set = dataset[train_size:]
+        return self.train_set, self.test_set
+    
+    def first_iteration_step(self, datasets: List[ModelBasicDatasetUnit],
+                             param_degree: int = 1):
+        '''
+        SK 迭代的第一步 (t=0) – 对应论文公式 (3).
+
+        目标:
+            最小化  sum_{样本 n, 频点 f} | N^{(0)}(s_f, g_n) - H(s_f, g_n) |^2
+        在 t=0 阶段我们令 D^{(0)}(s,g)=1, 仅拟合分子:
+            N^{(0)}(s,g) = Σ_{p∈P} Σ_{v∈V} c_{p,v} * φ_freq_p(s) * φ_param_v(g)
+
+        因此是一个纯线性最小二乘:
+            H ≈  Σ_{p,v} c_{p,v}  F[:,p] ⊗ G[n,v]
+
+        记:
+            F: (Nf, P)  频率正交(或原始)基列
+            G: (Ns, V)  参数多项式（含常数）基列
+            设计矩阵 A 大小 (Ns*Nf, P*V), A[(n*Nf + f), (p*V + v)] = F[f,p] * G[n,v]
+
+        输出:
+            self.first_iter_coeffs[(i,j)] = C_{p,v}  (P × V)  每个端口对一套系数
+            self.freq_basis_F = F
+            self.param_basis_G = G
+            self.poles 初始极点 (供后续构造有理正交基 / SK 加权使用)
+
+        参数:
+            datasets: 训练用样本列表
+            param_degree: 参数多项式最大总次数 (默认 1 → 常数 + 线性)
+
+        注意:
+            这里使用简单频率基 [1, 1/(s-a_k)]，未做 QR 正交化;
+            后续 SK 迭代 / 正交化可在第二步再进行。
+        '''
+        assert len(datasets) > 0, "空数据集"
+        # ---------------- 收集频率与端口信息 ----------------
+        Ns = len(datasets)
+        freqs = datasets[0].freqs  # 频率需要对齐
+        Nf = freqs.shape[0]
+        nports = self.model.info.nports
+
+        # ---------------- 构造初始极点 (对数分布 + 阻尼) ----------------
+        def _init_poles(n_poles: int):
+            fmin, fmax = freqs[0], freqs[-1]
+            if n_poles <= 0:
+                return np.array([], dtype=complex)
+            # 避免 0 Hz，用第二个点或微小偏移
+            start_f = max(fmin if fmin > 0 else freqs[1] if Nf > 1 else 1.0e-3, 1e-12)
+            f_samples = np.logspace(np.log10(start_f), np.log10(fmax), n_poles)
+            sigma = 0.02 * 2 * np.pi * fmax  # 固定阻尼，可放入 config
+            return -sigma + 1j * 2 * np.pi * f_samples
+
+        self.poles = _init_poles(self.config.n_poles)
+        s_vec = 1j * 2 * np.pi * freqs  # (Nf,)
+
+        # ---------------- 频率基 F ----------------
+        # 列0: 常数 1, 后续列: 1/(s - a_k)
+        P = 1 + len(self.poles)
+        F = np.zeros((Nf, P), dtype=complex)
+        F[:, 0] = 1.0
+        for k, a in enumerate(self.poles, start=1):
+            F[:, k] = 1.0 / (s_vec - a)
+        self.freq_basis_F = F  # (Nf, P)
+
+        # ---------------- 参数基 G ----------------
+        # 取出参数向量 g = [param1, param2, ...]
+        # 假设 datasets[i].parameters 为 dict
+        param_keys = list(datasets[0].parameters.keys())
+        d = len(param_keys)
+        # 构造参数矩阵 (Ns,d)
+        Pmat = np.zeros((Ns, d), dtype=float)
+        for i, unit in enumerate(datasets):
+            for j, k in enumerate(param_keys):
+                Pmat[i, j] = float(unit.parameters[k])
+        # 标准化
+        mean = Pmat.mean(axis=0)
+        std = Pmat.std(axis=0) + 1e-15
+        Xn = (Pmat - mean) / std
+
+        # 生成多项式指数 (总次数 <= param_degree)
+        def _gen_param_exps(dim, D):
+            exps = []
+            def rec(cur, idx, rem):
+                if idx == dim:
+                    exps.append(tuple(cur)); return
+                for t in range(rem + 1):
+                    cur.append(t); rec(cur, idx + 1, rem - t); cur.pop()
+            rec([], 0, D)
+            return exps
+
+        exps = _gen_param_exps(d, param_degree)  # V_exps
+        V = len(exps)
+        G = np.zeros((Ns, V), dtype=float)
+        for vidx, e in enumerate(exps):
+            val = np.ones(Ns)
+            for j, p in enumerate(e):
+                if p:
+                    val *= Xn[:, j] ** p
+            G[:, vidx] = val
+        self.param_basis_G = G  # (Ns, V)
+        self.param_basis_meta = {
+            "keys": param_keys,
+            "mean": mean.tolist(),
+            "std": std.tolist(),
+            "exponents": exps,
+            "degree": param_degree
+        }
+
+        # ---------------- 构造设计矩阵 A (Kronecker) ----------------
+        # A shape: (Ns*Nf, P*V)
+        # 利用广播：A_block[n,f,p,v] = F[f,p] * G[n,v]
+        A4 = np.einsum('fp,nv->nfpv', F, G)   # (Ns, Nf, P, V)
+        A = A4.reshape(Ns * Nf, P * V)
+
+        # ---------------- 采集目标 H 数据 ----------------
+        # 对每个端口对 (i,j) 分别解一个向量 c_{p,v}
+        # 选择参数类型
+        param_type = self.config.parameter_type.lower()
+        valid_types = {"s", "y", "z"}
+        if param_type not in valid_types:
+            raise ValueError(f"parameter_type 必须在 {valid_types}")
+
+        # 预先载入全部 Network (避免重复 IO)
+        # H_data[(i,j)] -> (Ns,Nf) 复数
+        H_data = {}
+        for i_port in range(nports):
+            for j_port in range(nports):
+                H_mat = np.zeros((Ns, Nf), dtype=complex)
+                for n, unit in enumerate(datasets):
+                    net: Network = unit.network
+                    if param_type == "s":
+                        sij = net.s[:, i_port, j_port]
+                    elif param_type == "y":
+                        sij = net.y[:, i_port, j_port]
+                    else:
+                        sij = net.z[:, i_port, j_port]
+                    # 插值或对齐假设已完成
+                    H_mat[n, :] = sij
+                H_data[(i_port, j_port)] = H_mat
+
+        # ---------------- 最小二乘求系数 ----------------
+        self.first_iter_coeffs = {}
+        # 可选: 预计算 A^+ (伪逆) 如果数据不巨大
+        # pinv = np.linalg.pinv(A)
+        for key, Hmat in H_data.items():
+            b = Hmat.reshape(Ns * Nf)  # (Ns*Nf,)
+            # 解 x  (长度 P*V)
+            x, *_ = np.linalg.lstsq(A, b, rcond=None)
+            C = x.reshape(P, V)
+            self.first_iter_coeffs[key] = C
+
+        # ---------------- 存储便于后续 SK 迭代使用 ----------------
+        self.meta_first_iter = {
+            "A_shape": A.shape,
+            "num_samples": Ns,
+            "num_freqs": Nf,
+            "freq_basis_dim": P,
+            "param_basis_dim": V
+        }
+        if getattr(self.config, "verbose", True):
+            print(f"[t=0] 线性最小二乘完成: A={A.shape}, 频率基P={P}, 参数基V={V}, 端口对={nports*nports}")
+
+        return self.first_iter_coeffs
+
+        

core/sk_iter.py

@@ -0,0 +1,250 @@
+import numpy as np
+from dataclasses import dataclass
+from typing import List, Sequence, Tuple, Optional
+
+@dataclass
+class OPVFConfig:
+    n_poles: int
+    max_iter: int = 5
+    tol: float = 1e-3
+    add_constraint: bool = True
+    real_split: bool = True
+    lambda_reg: float = 0.0
+    verbose: bool = True
+
+# ---------- 参数正交基 ----------
+class ParamOrthoBasis:
+    def __init__(self, params: np.ndarray, total_degree: int):
+        # params: (Ns, d)
+        self.mean = params.mean(0)
+        self.std = params.std(0) + 1e-15
+        self.exps = self._gen_exps(params.shape[1], total_degree)
+        M = self._build_monomial_matrix(params)  # (Ns,M)
+        Q,R = np.linalg.qr(M)
+        self.Q = Q
+        self.R = R
+    def _gen_exps(self, dim, D):
+        exps=[]
+        def rec(cur, i, rem):
+            if i==dim:
+                exps.append(tuple(cur)); return
+            for k in range(rem+1):
+                cur.append(k); rec(cur, i+1, rem-k); cur.pop()
+        rec([],0,D)
+        return exps
+    def _build_monomial_matrix(self, params):
+        X = (params - self.mean)/self.std
+        Ns = X.shape[0]
+        M = np.zeros((Ns, len(self.exps)))
+        for idx,e in enumerate(self.exps):
+            v = np.ones(Ns)
+            for j,p in enumerate(e):
+                if p: v *= X[:,j]**p
+            M[:,idx]=v
+        return M
+    def eval(self, g: np.ndarray) -> np.ndarray:
+        x = (g - self.mean)/self.std
+        mono=[]
+        for e in self.exps:
+            val=1.0
+            for j,p in enumerate(e):
+                if p: val*= x[j]**p
+            mono.append(val)
+        mono=np.array(mono)
+        # φ = mono @ R^{-1}
+        return mono @ np.linalg.inv(self.R)
+
+# ---------- 频率正交(有理)基 ----------
+class RationalFreqBasis:
+    def __init__(self, freqs: np.ndarray, poles: np.ndarray):
+        # 构造原始矩阵 G: 列0 常数 1, 后面 1/(s-a_k)
+        s = 1j*2*np.pi*freqs
+        cols=[np.ones_like(s)]
+        for a in poles:
+            cols.append(1.0/(s - a))
+        G = np.vstack(cols).T  # (Nf, K+1)
+        # 加权可加 w_f，这里统一 1
+        Q, R = np.linalg.qr(G)
+        self.Q = Q            # Φ(f)  (Nf, K+1) 频率正交基取样
+        self.R = R            # 变换: G = Q R
+        self.freqs = freqs
+    def eval(self) -> np.ndarray:
+        return self.Q   # 直接返回已正交基取样 (Nf, K+1)
+
+# ---------- 主模型 ----------
+class OrthonormalParametricVF:
+    """
+    H(s, μ) = N(s,μ)/D(s,μ)
+    N(s,μ)= Σ_k c_k(μ) φ_k(s);  D(s,μ)= Σ_k ĉ_k(μ) φ_k(s)
+    其中 φ_k(s) 是频率正交有理基; c_k(μ), ĉ_k(μ) 在参数正交基上展开:
+    c_k(μ)= Σ_q C[k,q] ψ_q(μ);  ĉ_k(μ)= Σ_q Ctilde[k,q] ψ_q(μ)
+    """
+    def __init__(self, cfg: OPVFConfig, freq_basis: RationalFreqBasis, param_basis: ParamOrthoBasis):
+        self.cfg = cfg
+        self.fb = freq_basis
+        self.pb = param_basis
+        Kp1 = self.fb.Q.shape[1]
+        Qp = self.pb.Q.shape[1]
+        self.C = np.zeros((Kp1, Qp), dtype=complex)        # 分子
+        self.Ct = np.zeros((Kp1, Qp), dtype=complex)       # 分母
+        # 初始化分母：ĉ_0 ≈ 1，其余 0
+        self.Ct[0,0] = 1.0
+
+    def _assemble_phi_param(self, params: np.ndarray) -> np.ndarray:
+        # 返回 (Ns, Qp)
+        return self.pb.Q
+
+    def fit(self, H_samples: List[np.ndarray], params: np.ndarray):
+        """
+        H_samples: list 长度 Ns, 每个 (Nf,) 复
+        params: (Ns,d)
+        """
+        Φf = self.fb.eval()           # (Nf, K+1)
+        Ns = len(H_samples)
+        Nf, Kp1 = Φf.shape
+        Φμ = self._assemble_phi_param(params)  # (Ns, Qp)
+        Qp = Φμ.shape[1]
+
+        # t=0: 只解 C (式3, D=1)
+        self._solve_t0(H_samples, Φf, Φμ)
+
+        D_prev = self._eval_D(Φf, Φμ)  # (Nf, Ns)
+
+        for it in range(1, self.cfg.max_iter+1):
+            A, b = self._build_iter_system(H_samples, Φf, Φμ, D_prev)
+            if self.cfg.lambda_reg>0:
+                lam = self.cfg.lambda_reg
+                A = np.vstack([A, lam*np.eye(A.shape[1])])
+                b = np.concatenate([b, np.zeros(A.shape[1])])
+            x, *_ = np.linalg.lstsq(A, b, rcond=None)
+            self._unpack_iter_solution(x, Kp1, Qp)
+            D_new = self._eval_D(Φf, Φμ)
+            rel = np.max(np.abs(D_new-D_prev)/(np.abs(D_prev)+1e-12))
+            if self.cfg.verbose:
+                print(f"[OPVF-Iter {it}] rel_change={rel:.3e}")
+            D_prev = D_new
+            if rel < self.cfg.tol:
+                if self.cfg.verbose:
+                    print(f"[OPVF] converged at {it}")
+                break
+        return self
+
+    def _solve_t0(self, H_samples, Φf, Φμ):
+        Ns = len(H_samples); Nf,Kp1 = Φf.shape; Qp = Φμ.shape[1]
+        # 设计矩阵行数 Ns*Nf；未知数 Kp1*Qp
+        cols = Kp1*Qp
+        A = np.zeros((Ns*Nf, cols), dtype=complex)
+        b = np.zeros(Ns*Nf, dtype=complex)
+        r = 0
+        for n,Hn in enumerate(H_samples):
+            phiμ = Φμ[n]  # (Qp,)
+            # Φf (Nf,Kp1) 外积 phiμ -> (Nf, Kp1*Qp)
+            blk = np.einsum('fk,q->fkq', Φf, phiμ).reshape(Nf, cols)
+            A[r:r+Nf,:] = blk
+            b[r:r+Nf] = Hn
+            r += Nf
+        x, *_ = np.linalg.lstsq(A, b, rcond=None)
+        self.C = x.reshape(Kp1, Qp)
+        # 分母保持初始 (ĉ_0=1)
+
+    def _build_iter_system(self, H_samples, Φf, Φμ, D_prev):
+        Ns = len(H_samples); Nf,Kp1 = Φf.shape; Qp=Φμ.shape[1]
+        # 未知：C (Kp1*Qp) + Ct (Kp1*Qp) 但固定 Ct[0,0]=1，可去掉该变量
+        mask_fix = np.zeros((Kp1,Qp), dtype=bool)
+        mask_fix[0,0]=True
+        idx_map={}
+        col=0
+        for i in range(Kp1):
+            for j in range(Qp):
+                idx_map[('C',i,j)] = col; col+=1
+        for i in range(Kp1):
+            for j in range(Qp):
+                if mask_fix[i,j]: continue
+                idx_map[('Ct',i,j)] = col; col+=1
+        cols = col
+        rows = Ns*Nf
+        A = np.zeros((rows, cols), dtype=complex)
+        b = np.zeros(rows, dtype=complex)
+        r=0
+        for n,Hn in enumerate(H_samples):
+            phiμ = Φμ[n]
+            Dp = D_prev[:,n]  # (Nf,)
+            invDp = 1.0/(Dp + 1e-12)
+            # 分子块: (Φf φμ^T) / D_prev
+            Num_blk = np.einsum('fk,q->fkq', Φf, phiμ).reshape(Nf, Kp1*Qp)
+            Num_blk = (Num_blk.T * invDp).T
+            # 填 C 部分
+            for i in range(Kp1):
+                for j in range(Qp):
+                    A[r:r+Nf, idx_map[('C',i,j)]] = Num_blk[:, i*Qp + j]
+            # 分母块: -(H * Φf φμ^T)/D_prev
+            Den_blk = (Num_blk.T * Hn).T * (-1.0)
+            for i in range(Kp1):
+                for j in range(Qp):
+                    if mask_fix[i,j]: continue
+                    A[r:r+Nf, idx_map[('Ct',i,j)]] = Den_blk[:, i*Qp + j]
+            # 右端 0
+            r += Nf
+
+        # 约束行 (式5)(8) 可选：Re( Σ_k D^{(t)}/D^{(t-1)} ) = K+1
+        if self.cfg.add_constraint:
+            row_c = np.zeros(cols, dtype=complex)
+            # D^{(t)} ≈ Σ Ct_k φ_k(s) ; 用代表样本 n=0, 对所有频率平均
+            n0=0
+            phiμ0=Φμ[n0]
+            mean_phiμ = phiμ0  # 可换平均
+            # φ_k(s) 平均
+            mean_phif = Φf.mean(0)  # (Kp1,)
+            for i in range(Kp1):
+                for j in range(Qp):
+                    if mask_fix[i,j]: continue
+                    row_c[idx_map[('Ct',i,j)]] = mean_phif[i]*mean_phiμ[j]
+            rhs = (Kp1)  # K+1
+            A = np.vstack([A, np.real(row_c) if self.cfg.real_split else row_c])
+            b = np.concatenate([b, [rhs]])
+
+        # 拆实虚
+        if self.cfg.real_split:
+            A_real = np.vstack([np.real(A), np.imag(A)])
+            b_real = np.concatenate([np.real(b), np.imag(b)])
+        else:
+            A_real, b_real = A, b
+        return A_real, b_real
+
+    def _unpack_iter_solution(self, x, Kp1, Qp):
+        # 重新填回 C, Ct (保持 Ct[0,0]=1)
+        # 构建与 _build_iter_system 相同的 idx_map
+        mask_fix = np.zeros((Kp1,Qp), dtype=bool); mask_fix[0,0]=True
+        idx_C = Kp1*Qp
+        # 注意：我们当时 C 索引从 0..(Kp1*Qp-1)
+        self.C = x[:idx_C].reshape(Kp1, Qp)
+        Ct_new = self.Ct.copy()
+        pos = idx_C
+        for i in range(Kp1):
+            for j in range(Qp):
+                if mask_fix[i,j]: continue
+                Ct_new[i,j]=x[pos]; pos+=1
+        self.Ct = Ct_new
+
+    def _eval_C_mu(self, phiμ):
+        # 返回 c_k(μ) (Kp1,)
+        return self.C @ phiμ
+    def _eval_Ct_mu(self, phiμ):
+        return self.Ct @ phiμ
+    def _eval_D(self, Φf, Φμ):
+        # D(f, sample)= Σ_k ĉ_k(μ_n) φ_k(f)
+        Kp1 = Φf.shape[1]
+        Ns = Φμ.shape[0]
+        D = np.zeros((Φf.shape[0], Ns), dtype=complex)
+        for n in range(Ns):
+            ct_mu = self._eval_Ct_mu(Φμ[n])
+            D[:,n] = Φf @ ct_mu
+        return D
+
+    def evaluate(self, freqs: np.ndarray, g: np.ndarray):
+        assert np.allclose(freqs, self.fb.freqs), "频率需在训练网格上 (示例简化)"
+        Φf = self.fb.Q            # (Nf,K+1)
+        phiμ = self.pb.eval(g)    # (Qp,)
+        num = Φf @ (self.C @ phiμ)
+        den = Φf @ (self.Ct @ phiμ)
+        return num / (den + 1e-15)