From e62e8df0134453824485aa35314049218ed3bea0 Mon Sep 17 00:00:00 2001 From: mayge Date: Sat, 20 Sep 2025 12:02:24 -0400 Subject: [PATCH] =?UTF-8?q?feat:=20=E5=AE=8C=E6=88=90=E4=BA=86=E7=AE=97?= =?UTF-8?q?=E6=B3=95=E7=9A=84=E6=9C=80=E5=9F=BA=E7=A1=80=E9=83=A8=E5=88=86?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 3 +- README.md | 6 + core/basis.py | 241 ++++++++++++++++ core/freqency.py | 131 +++++++++ core/orthonormal_basis.py | 121 -------- core/robust.py | 191 ------------- core/sk_iter.py | 306 -------------------- img_formula_70.png | Bin 49695 -> 0 bytes main.py | 1 + param_vf.py | 583 -------------------------------------- 10 files changed, 381 insertions(+), 1202 deletions(-) create mode 100644 README.md create mode 100644 core/basis.py create mode 100644 core/freqency.py delete mode 100644 img_formula_70.png delete mode 100644 param_vf.py diff --git a/.gitignore b/.gitignore index 3117b89..2be33ae 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,5 @@ __pycache__/ *.pyd __pypackages__/ env/ -.venv/ \ No newline at end of file +.venv/ +.vscode/ \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..16c0227 --- /dev/null +++ b/README.md @@ -0,0 +1,6 @@ +1. 添加 A B C D 方法，对应不同的基 +2. 修改Levi部分和sk迭代部分，使用更通用的A B C D求解 +3. 添加正交基形式 +4. 加入多端口 +5. 引入QR，消除中间过程的残差求解，降低算法复杂度 +6. 使用最初始的表达式，获得无偏估计 \ No newline at end of file diff --git a/core/basis.py b/core/basis.py new file mode 100644 index 0000000..2f16635 --- /dev/null +++ b/core/basis.py @@ -0,0 +1,241 @@ +import numpy as np +from core.sk_iter import generate_starting_poles +from scipy.linalg import block_diag +import skrf as rf +from skrf import VectorFitting +from core.freqency import auto_select + +class BasicBasis: + def __init__(self,H,freqs,poles,weights=None): + self.least_squares_rms_error = None + self.least_squares_condition = None + self.eigenval_condition = None + self.eigenval_rms_error = None + + self.H = H + self.freqs = freqs + self.s = self.freqs * 2j * np.pi + self.P = len(poles) + self.poles = poles + self.Phi = self.generate_basis(self.s, self.poles) + self.A = self.matrix_A(self.poles) + self.B = self.vector_B(self.poles) + self.D = 1.0 + self.Cr,self.Cw = self.fit_denominator(self.H, d0=1.0, weights=weights) + + z = np.linalg.eigvals(self.A - self.B @ self.Cw) + p_next = -z + + # enforce LHP and pair ordering + p_next = np.where(np.real(p_next) < 0, p_next, -np.conj(p_next)) + p_next = np.sort_complex(p_next) + + self.next_poles = p_next + # z = np.where(np.real(z) < 0, z, -np.conj(z)) # enforce LHP + # self.next_poles = np.sort_complex(z) + self.eigenval_condition = np.linalg.cond(self.A - self.B @ self.Cw) + self.eigenval_rms_error = np.sqrt(np.mean(np.abs(np.real(z) - np.real(poles))**2 + np.abs(np.imag(z) - np.imag(poles))**2)) + + self.Dt = self.eval_Dt_state_space() + self.delta = self.Dt / weights if weights is not None else self.Dt + pass + + def eval_Dt_state_space(self): + """Return D(s_k)=C(s_k I - A)^(-1)B + D for all k (complex 1D array).""" + s = 1j * 2*np.pi * np.asarray(self.freqs, float).ravel() + A = np.asarray(self.A, np.complex128); n = A.shape[0] + B = np.asarray(self.B, np.complex128).reshape(n, 1) + C = np.asarray(self.Cw, float).reshape(1, n) + D = self.D + I = np.eye(n, dtype=np.complex128) + out = np.empty_like(s, dtype=np.complex128) + for k, sk in enumerate(s): + DS = D + (C @ np.linalg.inv(sk*I - A) @ B) + out[k] = DS[0, 0] + return out + + + def generate_basis(self,s, poles): + """Real basis of (15)-(16); returns Φ(s) and a layout for packing C.""" + cols = [] + i = 0 + while i < len(poles): + p = poles[i] + if p.real > 0: + raise ValueError("poles must be in the LHP") + if i+1 < len(poles) and np.isclose(poles[i+1], np.conj(p)): + pc = poles[i+1] + phi1 = 1/(s - p) + 1/(s - pc) # eq (15)generate_basis + phi2 = 1j*(1/(s - p) - 1/(s - pc)) # eq (16) (fixed sign) + cols += [phi1, phi2] + i += 2 + else: + cols.append(1/(s - p)) + i += 1 + Phi = np.column_stack(cols).astype(np.complex128) + return Phi + + def matrix_A(self, poles): + def A_block(p): + if abs(p.imag) < 1e-14: + return np.array([[p.real]], float) # A_p = [ p ] + return np.array([[p.real, p.imag], # A_p = [[Re p, Im p], + [-p.imag, p.real]], float) # [-Im p, Re p]] + A = None; i = 0 + while i < len(poles): + p = poles[i] + Ab = A_block(p) + if i+1 < len(poles) and np.isclose(poles[i+1], np.conj(p)): i += 2 + else: i += 1 + A = Ab if A is None else block_diag(A, Ab) + return A + + def vector_B(self, poles): + def B_block(p): + return np.array([[1.0]], float) if abs(p.imag)<1e-14 else np.array([[2.0],[0.0]], float) + B = None; i = 0 + while i < len(poles): + p = poles[i] + Bb = B_block(p) + if i+1 < len(poles) and np.isclose(poles[i+1], np.conj(p)): i += 2 + else: i += 1 + B = Bb if B is None else np.vstack([B, Bb]) + return B + + def fit_denominator(self, H, d0=1.0, weights=None): + """ + Solve formula (70) on the real basis Φ to obtain: + - d (real) → packs into C for this state's block structure + - gamma (complex) + Optional 'weights' (K,) apply row scaling: SK weighting if 1/|D_prev|. + """ + if weights is None: + weights = np.diag(np.ones(len(H), np.complex128)) + else: + weights = np.diag([1/res for res in weights]) + s = self.s + H = np.asarray(H, np.complex128).reshape(-1,1) + Phi = self.Phi + psi = weights @ Phi + psi = Phi + + HPhi = H * Phi + + A_re = np.hstack([np.real(-psi), np.real(-HPhi)]) + A_im = np.hstack([np.imag(-psi), np.imag(-HPhi)]) + + b_re = np.real(d0 * H) + b_im = np.imag(d0 * H) + + A = np.vstack([A_re, A_im]).astype(float) + # rown = np.linalg.norm(A, axis=1) + # rown = np.sqrt(rown) + # A = rown[:,None] * A + b = np.concatenate([b_re, b_im]).astype(float) + + x = np.linalg.inv(A.T @ A) @ A.T @ b + + self.least_squares_rms_error = np.sqrt(np.mean((A @ x - b)**2)) + self.least_squares_condition = np.linalg.cond(A) + + Cn,Cd = self.vector_C(x) + return Cn,Cd + + def vector_C(self,x): + Cn = np.asarray([x[:len(x)//2]], float).reshape(1,-1) + Cd = np.asarray([x[len(x)//2:]], float).reshape(1,-1) + return Cn, Cd + + def evaluate(self,freqs,poles,Cn,Cd,d0=1.0): + s = 1j * 2*np.pi * np.asarray(freqs, float).ravel() + phi = self.generate_basis(s, poles) + num = phi @ Cn.T + den = d0 + phi @ Cd.T + H = num / den + return H.ravel() + +if __name__ == "__main__": + network = rf.Network("/tmp/paramer/simulation/3000/3000.s2p") + K = 10 + H11,freqs = auto_select([network.y[i][0][0] for i in range(2,len(network.y))],network.f[2:],max_points=20) + poles = generate_starting_poles(2,beta_min=freqs[0]/1.1,beta_max=freqs[-1]*1.1) + + Dt_1 = np.ones((len(freqs),1),np.complex128) + # Levi step (no weighting): + basis = BasicBasis(H11,freqs,poles=poles) + Dt = basis.Dt + poles = basis.next_poles + + print("Levi step (no weighting):") + print("A:",basis.A) + print("B:",basis.B) + print("C:",basis.Cw) + print("D:",basis.D) + print("next_pozles:",basis.next_poles) + print("Dt:",Dt, "norm:",np.linalg.norm(Dt)) + + # SK weighting (optional, after first pass): + least_squares_condition = [] + least_squares_rms_error = [] + eigenval_condition = [] + eigenval_rms_error = [] + for i in range(K): + basis = BasicBasis(H11,freqs,poles=poles,weights=Dt) + Dt_1 = Dt + Dt = basis.Dt + poles = basis.next_poles + print(f"SK Iteration {i+1}/{K}") + print("A:",basis.A) + print("B:",basis.B) + print("C:",basis.Cw) + print("D:",basis.D) + print("z:",basis.next_poles) + print("Dt:",Dt) + print("Dt/Dt-1",np.linalg.norm(Dt) / np.linalg.norm(Dt_1)) + least_squares_condition.append(basis.least_squares_condition) + least_squares_rms_error.append(basis.least_squares_rms_error) + eigenval_condition.append(basis.eigenval_condition) + eigenval_rms_error.append(basis.eigenval_rms_error) + + # H11_evaluated = basis.evaluate_pole_residue(network.f[1:],poles,basis.C[0]) + H11_evaluated = basis.evaluate(network.f[2:], poles, basis.Cr[0],basis.Cw[0], d0=1.0) + import matplotlib.pyplot as plt + fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) + + ax0 = axes[0][0] + ax0.plot(network.f[2:], np.abs([network.y[i][0][0] for i in range(2,len(network.y))]), 'o', ms=4, color='red', label='Samples') + ax0.plot(network.f[2:], np.abs(H11_evaluated), '-', lw=2, color='k', label='Fit') + ax0.plot(freqs, np.abs(H11), 'x', ms=4, color='blue', label='Input Samples') + ax0.set_title("Response i=0, j=0") + ax0.set_ylabel("Magnitude") + ax0.legend(loc="best") + + ax1 = axes[1][0] + ax1.plot(least_squares_condition, label='Least Squares Condition') + ax1.set_title("least_squares_condition") + ax1.set_ylabel("Magnitude") + ax1.legend(loc="best") + + ax2 = axes[1][1] + ax2.plot(least_squares_rms_error, label='Least Squares RMS Error') + ax2.set_title("least_squares_rms_error") + ax2.set_ylabel("Magnitude") + ax2.legend(loc="best") + + ax3 = axes[2][0] + ax3.plot(eigenval_condition, label='Eigenvalue Condition') + ax3.set_title("eigenval_condition") + ax3.set_ylabel("Magnitude") + ax3.legend(loc="best") + + ax4 = axes[2][1] + ax4.plot(eigenval_rms_error, label='Eigenvalue RMS Error') + ax4.set_title("eigenval_rms_error") + ax4.set_ylabel("Magnitude") + ax4.legend(loc="best") + fig.tight_layout() + plt.savefig(f"basic_basis.png") + + + + diff --git a/core/freqency.py b/core/freqency.py new file mode 100644 index 0000000..617d296 --- /dev/null +++ b/core/freqency.py @@ -0,0 +1,131 @@ +import numpy as np +def auto_select(H, freq, + n_baseline=64, # log-spaced backbone points + peak_prominence=0.05, # fraction of |H| dB dynamic range for peak detection + peak_window=5, # take ±peak_window samples around each peak + topgrad_q=0.98, # keep top 2% largest slope/phase-change points + max_points=25, # final cap on selected samples (None = no cap) + ensure_ends=True): + """ + Select several significant sample points for vector fitting. + + Strategy: + 1) Always keep endpoints (optional). + 2) Add a log-spaced baseline over the band. + 3) Detect resonance peaks in |H| (on a log scale) and keep small windows around them. + 4) Add points with the largest magnitude slope and phase-change (w.r.t log-f). + 5) De-duplicate, sort, and optionally thin to 'max_points' with priority + to endpoints and detected peaks. + + Parameters + ---------- + H : (N,) complex array + Frequency response samples. + freq : (N,) float array + Frequency axis [Hz], strictly increasing. + n_baseline : int + Count of log-spaced baseline samples across the band. + peak_prominence : float + Peak prominence threshold as a fraction of the dynamic range in log|H|. + 0.05 ≈ keep peaks ≥ 5% of the range. + peak_window : int + Number of neighbor indices to include on each side of every detected peak. + topgrad_q : float in (0,1) + Quantile for selecting strong slope/phase points. + 0.98 ⇒ keep the top 2% largest derivatives. + max_points : int or None + If not None, cap the total number of selected indices to this value. + ensure_ends : bool + Always include the first and last samples. + + Returns + ------- + H_sel : (K,) complex array + freq_sel : (K,) float array + """ + H = np.asarray(H).reshape(-1) + f = np.asarray(freq).reshape(-1) + if H.size != f.size: + raise ValueError("H and freq must have the same length.") + N = f.size + if N < 4: + return H.copy(), f.copy() + + eps = 1e-16 + mag = np.abs(H) + logmag = np.log10(mag + eps) + phase = np.unwrap(np.angle(H)) + + # log-frequency axis (scale-invariant derivatives) + # keep it linear if any non-positive freq sneaks in + if np.all(f > 0): + lf = np.log(f) + else: + lf = f.copy() + + dlf = np.gradient(lf) + d_logmag = np.gradient(logmag) / (dlf + 1e-16) + d_phase = np.gradient(phase) / (dlf + 1e-16) + + idx = set() + if ensure_ends: + idx.update([0, N-1]) + + # 1) log-spaced baseline + if n_baseline > 0: + # map a log grid to nearest indices + grid = np.linspace(lf.min(), lf.max(), n_baseline) + base_idx = np.clip(np.searchsorted(lf, grid), 0, N-1) + idx.update(np.unique(base_idx).tolist()) + + # 2) peaks in |H| + try: + from scipy.signal import find_peaks + dyn = logmag.max() - logmag.min() + prom = peak_prominence * (dyn + 1e-12) + peaks, _ = find_peaks(logmag, prominence=prom) + except Exception: + # simple fallback: strict local maxima + peaks = np.where((mag[1:-1] > mag[:-2]) & (mag[1:-1] > mag[2:]))[0] + 1 + + for p in peaks: + lo = max(0, p - peak_window) + hi = min(N, p + peak_window + 1) + idx.update(range(lo, hi)) + + # 3) strongest slope / phase-change points + thr_slope = np.quantile(np.abs(d_logmag), topgrad_q) + thr_phase = np.quantile(np.abs(d_phase), topgrad_q) + idx.update(np.where(np.abs(d_logmag) >= thr_slope)[0].tolist()) + idx.update(np.where(np.abs(d_phase) >= thr_phase)[0].tolist()) + + # 4) finalize set + sel = np.array(sorted(idx), dtype=int) + + # 5) optional thinning with priority to endpoints and peaks + if max_points is not None and sel.size > max_points: + priority = np.zeros(sel.size, dtype=int) + if ensure_ends: + priority[(sel == 0) | (sel == N-1)] = 3 + if peaks.size: + priority[np.isin(sel, peaks)] = np.maximum(priority[np.isin(sel, peaks)], 2) + + keep = [] + budget = max_points + # keep highest-priority first + for lev in (3, 2, 1, 0): + cand = sel[priority == lev] + if cand.size == 0: + continue + if cand.size <= budget: + keep.extend(cand.tolist()) + budget -= cand.size + else: + step = max(1, int(np.ceil(cand.size / budget))) + keep.extend(cand[::step][:budget].tolist()) + budget = 0 + if budget == 0: + break + sel = np.array(sorted(set(keep)), dtype=int) + + return H[sel], f[sel] \ No newline at end of file diff --git a/core/orthonormal_basis.py b/core/orthonormal_basis.py index 096dea2..1f64ac7 100644 --- a/core/orthonormal_basis.py +++ b/core/orthonormal_basis.py @@ -1,36 +1,5 @@ import numpy as np -# ------------------------------------------------------------ -# 按论文公式 (9)(10)(11) 生成 Muntz–Laguerre 正交有理基 (解析形式)： -# -# 给定稳定极点集合 {p_k} (Re(p_k)<0)。论文记法中使用 -a_k，其中 Re(a_k)>0。 -# 对应关系：p_k = -a_k ⇒ a_k = -p_k, Re(a_k)= -Re(p_k) >0 -# -# 连续内积意义下（沿 jω 轴积分）这些 φ_k 解析正交。离散频率采样后数值上 -# 可能偏离，可再用加权 QR 做数值再正交（可选）。 -# -# 公式在“稳定极点 p 表达”下的改写： -# (原) 实极点: φ_p(s) = sqrt(2 Re(a_p)) / (s + a_p) * Π (s - a_i^*)/(s + a_i) -# 变换 a_p = -p ⇒ Re(a_p)= -Re(p) = σ >0 且 (s + a_p) = (s - p) -# 且乘积 (s - a_i^*)/(s + a_i) = (s + p_i^*)/(s - p_i) -# ⇒ φ_p(s) = sqrt(-2 Re(p)) / (s - p) * Π_{i0 的 p 作为首： -# (原) φ_p = sqrt(2 Re(a_p)) (s - |a_p|)/[(s + a_p)(s + a_p^*)] * Π(...) -# (原) φ_{p+1}= sqrt(2 Re(a_p)) (s + |a_p|)/[(s + a_p)(s + a_p^*)] * Π(...) -# 代入 a_p=-p： -# Re(a_p)= -Re(p)=σ>0, (s + a_p) = (s - p), (s + a_p^*)=(s - p^*) -# |a_p| = |p| -# 乘积同上 ⇒ Π_{i= self.N: -# raise StopIteration - -# p = self.poles[self.k] -# if np.real(p) >= 0: -# raise ValueError(f"极点必须在左半平面: {p}") - -# # 复对首 (正虚部) -# if np.iscomplex(p) and np.imag(p) > 0: -# if self.k + 1 >= self.N: -# raise ValueError("复极点缺少共轭") -# pc = self.poles[self.k + 1] -# if not np.isclose(pc, np.conj(p)): -# raise ValueError("复极点未按 (p, p*) 顺序排列 (正虚部在前)") -# sigma = -np.real(p) # >0 -# scale = np.sqrt(2 * sigma) -# r = np.abs(p) -# denom = (self.s - p) * (self.s - pc) - -# # 两个基函数 -# phi_p = scale * (self.s - r) / denom * self.product -# phi_pc = scale * (self.s + r) / denom * self.product - -# # product 先乘 (s + p^*)/(s - p)，再乘 (s + p)/(s - p^*) -# self.product = self.product * (self.s + pc) / (self.s - p) -# self.product = self.product * (self.s + p) / (self.s - pc) - -# self.k += 2 -# return [phi_p, phi_pc] - -# # 复对次 (负虚部) —— 应该被首元素处理，出现表示顺序错误 -# if np.iscomplex(p) and np.imag(p) < 0: -# raise ValueError("检测到负虚部复极点但其共轭尚未处理，请将正虚部成员放在前面。") - -# # 实极点 -# sigma = -np.real(p) -# if sigma <= 0: -# raise ValueError("实极点实部应为负 (稳定)。") -# scale = np.sqrt(2 * sigma) -# phi = scale / (self.s - p) * self.product -# # 更新乘积 -# self.product = self.product * (self.s + p) / (self.s - p) - -# self.k += 1 -# return [phi] - -# def generate_muntz_laguerre_basis(s: np.ndarray, init_poles: list | np.ndarray): -# """ -# 生成完整基函数列表: [φ_0=1, φ_1, φ_2, ...] -# """ -# basis = [np.ones_like(s, dtype=complex)] -# for block in MuntzLaguerreIterator(s, init_poles): -# basis.extend(block) -# return basis - -# if __name__ == "__main__": -# # 示例稳定极点 (复对正虚部在前) -# stable_poles = [ -# -0.8e9, -# -1.0e9 + 2.5e9j, -# -1.0e9 - 2.5e9j, -# -2.2e9 -# ] -# freqs = np.linspace(1e8, 8e9, 400) -# s = 1j * 2 * np.pi * freqs - -# basis = generate_muntz_laguerre_basis(s, stable_poles) -# print(f"生成 {len(basis)} 个基函数,{basis}") \ No newline at end of file diff --git a/core/robust.py b/core/robust.py index ec80240..e69de29 100644 --- a/core/robust.py +++ b/core/robust.py @@ -1,191 +0,0 @@ -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 - - \ No newline at end of file diff --git a/core/sk_iter.py b/core/sk_iter.py index f230398..47681f6 100644 --- a/core/sk_iter.py +++ b/core/sk_iter.py @@ -2,7 +2,6 @@ import numpy as np from dataclasses import dataclass from typing import List, Dict, Tuple -# ---------- 1. 生成初始极点（按论文 Re(-a_p) 小、Im 线性分布） ---------- def generate_starting_poles(n_pairs: int, beta_min: float, beta_max: float, alpha_scale: float = 0.01): """ 仅生成复共轭对: p = -alpha + j beta, p*。 @@ -19,308 +18,3 @@ def generate_starting_poles(n_pairs: int, beta_min: float, beta_max: float, alph poles += [p, np.conj(p)] print(f"生成 {len(poles)} 个初始极点 (复对) {poles}]") return poles - -# ---------- 2. Muntz–Laguerre 频率基 (与你现有 orthonormal_basis 中一致的核心) ---------- -def muntz_laguerre_basis(freqs_hz: np.ndarray, stable_poles: List[complex]) -> np.ndarray: - """ - 返回 Φ_ml (Nf, P) 列顺序: φ_0=1, 后续依次 (单极点 / 复对两列) - """ - s = 1j * 2 * np.pi * freqs_hz - basis = [np.ones_like(s, dtype=complex)] - product = np.ones_like(s, dtype=complex) - k = 0 - N = len(stable_poles) - while k < N: - p = stable_poles[k] - if np.real(p) >= 0: - raise ValueError("极点需在左半平面") - if np.imag(p) > 0: # 复对首 - if k + 1 >= N or not np.isclose(stable_poles[k+1], np.conj(p)): - raise ValueError("复对未配对") - pc = stable_poles[k+1] - sigma = -np.real(p) - scale = np.sqrt(2 * sigma) - r = np.abs(p) - denom = (s - p) * (s - pc) - phi_p = scale * (s - r) / denom * product - phi_pc = scale * (s + r) / denom * product - basis.extend([phi_p, phi_pc]) - # update product - product = product * (s + pc)/(s - p) * (s + p)/(s - pc) - k += 2 - elif np.imag(p) < 0: - raise ValueError("负虚部复极点必须跟在正虚部之后") - else: # 实极点 - sigma = -np.real(p) - scale = np.sqrt(2 * sigma) - phi = scale / (s - p) * product - basis.append(phi) - product = product * (s + p)/(s - p) - k += 1 - return np.column_stack(basis) # (Nf, P) - -# ---------- 3. 原始部分分式基 ---------- -def raw_partial_fraction_basis(freqs_hz: np.ndarray, stable_poles: List[complex]) -> np.ndarray: - s = 1j * 2 * np.pi * freqs_hz - cols = [np.ones_like(s, dtype=complex)] - for p in stable_poles: - cols.append(1.0 / (s - p)) - return np.column_stack(cols) # (Nf, P) - -# ---------- 4. 变换矩阵 Φ_ml T = G_raw ---------- -def compute_transform_matrix(Phi_ml: np.ndarray, G_raw: np.ndarray, weights: np.ndarray | None = None) -> np.ndarray: - # 加权最小二乘 T = (Φ^H W Φ)^{-1} Φ^H W G - if weights is None: - WPhi = Phi_ml - WG = G_raw - M = Phi_ml.conj().T @ WPhi - RHS = Phi_ml.conj().T @ WG - else: - w = weights[:, None] - M = Phi_ml.conj().T @ (w * Phi_ml) - RHS = Phi_ml.conj().T @ (w * G_raw) - T = np.linalg.solve(M, RHS) - return T # (P,P), Φ T = G - -# ---------- 5. 参数多项式正交基 ---------- -def generate_param_basis(param_matrix: np.ndarray, total_degree: int): - """ - param_matrix: (Ns, d) - 返回: - Qμ: (Ns, V) 正交列 - Rμ: (V, V) 上三角 - exps: 多指数列表 - """ - X = param_matrix - Ns, d = X.shape - # 生成多指数 - exps=[] - def rec(cur,i,rem): - if i==d: - exps.append(tuple(cur)); return - for k in range(rem+1): - cur.append(k); rec(cur,i+1,rem-k); cur.pop() - rec([],0,total_degree) - V = len(exps) - M = np.zeros((Ns, V)) - for idx,e in enumerate(exps): - v = np.ones(Ns) - for j,p in enumerate(e): - if p: v *= X[:,j]**p - M[:,idx]=v - Q,R = np.linalg.qr(M) - return Q,R,exps - -# ---------- 6. t=0 构建设计矩阵并解 C ---------- -def fit_t0(H_list: List[np.ndarray], Phi_f: np.ndarray, Qμ: np.ndarray) -> np.ndarray: - """ - H_list: 长度 Ns, 每项 (Nf,) 复数 - Phi_f: (Nf, P) - Qμ: (Ns, V) - 返回: C (P,V) - """ - Ns = len(H_list); Nf,P = Phi_f.shape; V = Qμ.shape[1] - cols = P*V - A = np.zeros((Ns*Nf, cols), dtype=complex) - b = np.zeros(Ns*Nf, dtype=complex) - r=0 - for n,Hn in enumerate(H_list): - blk = np.einsum('fp,v->fpv', Phi_f, Qμ[n]).reshape(Nf, cols) - A[r:r+Nf,:]=blk - b[r:r+Nf]=Hn - r+=Nf - x, *_ = np.linalg.lstsq(A, b, rcond=None) - C = x.reshape(P, V) - return C - -# ---------- 7. t=1 (首轮 SK) 解 C, Ct ---------- -def fit_t1_SK(H_list: List[np.ndarray], Phi_f: np.ndarray, Qμ: np.ndarray, - C_prev: np.ndarray, max_iter: int =1, tol: float =1e-3): - """ - 只做一轮 (或少量) SK，初始 D=1。 - 返回: C_new, Ct_new - """ - Ns = len(H_list); Nf,P = Phi_f.shape; V = Qμ.shape[1] - # 初始化分母系数 Ct: Ct[0,0]=1 - Ct = np.zeros((P,V), dtype=complex) - Ct[0,0]=1.0 - C = C_prev.copy() - for it in range(max_iter): - # 构建线性系统: (N - D*H)=0 - # 未知顺序: C(:), Ct(:) 去掉固定 Ct[0,0] - mask_fix = np.zeros((P,V), dtype=bool); mask_fix[0,0]=True - col_map={} - col=0 - for i in range(P): - for j in range(V): - col_map[('C',i,j)]=col; col+=1 - for i in range(P): - for j in range(V): - if mask_fix[i,j]: continue - col_map[('Ct',i,j)]=col; col+=1 - A = np.zeros((Ns*Nf, col), dtype=complex) - b = np.zeros(Ns*Nf, dtype=complex) - r=0 - # 评价当前 D_prev= Σ Ct φ ψ - for n,Hn in enumerate(H_list): - # 分子块 (Φ_f ⊗ ψ_n) - blk = np.einsum('fp,v->fpv', Phi_f, Qμ[n]).reshape(Nf, P*V) - # 填 C 部分 - for i in range(P): - for j in range(V): - A[r:r+Nf, col_map[('C',i,j)]] = blk[:, i*V + j] - # 分母块: - H * blk - for i in range(P): - for j in range(V): - if mask_fix[i,j]: continue - A[r:r+Nf, col_map[('Ct',i,j)]] = - Hn * blk[:, i*V + j] - r+=Nf - # 求解 - x, *_ = np.linalg.lstsq(A, b, rcond=None) - # 拆回 - idx=0 - for i in range(P): - for j in range(V): - C[i,j]=x[idx]; idx+=1 - Ct_new = Ct.copy() - for i in range(P): - for j in range(V): - if mask_fix[i,j]: continue - Ct_new[i,j]=x[idx]; idx+=1 - # 收敛性简单检查：分母变化 - diff = np.max(np.abs(Ct_new - Ct)) - Ct = Ct_new - if diff < tol: - break - return C, Ct - -# ---------- 8. 假极点检测 (在 raw 基上) ---------- -def detect_spurious_poles(C_ml: np.ndarray, Ct_ml: np.ndarray, T: np.ndarray, - Qμ: np.ndarray, - tau_cancel=1e-2, tau_small=1e-3, eps=1e-14): - """ - C_ml, Ct_ml: (P,V) (Muntz-Laguerre 基) - T: Φ_ml T = G_raw - Qμ: (Ns,V) - """ - # raw 系数: c_raw = T^{-1} c_ml - Tinv = np.linalg.inv(T) - C_raw = Tinv @ C_ml - Ct_raw = Tinv @ Ct_ml - P,V = C_ml.shape - Ns = Qμ.shape[0] - r_vals = C_raw @ Qμ.T - q_vals = Ct_raw @ Qμ.T - metrics={} - scales=[] - for k in range(1,P): # 跳过常数列 - rv = r_vals[k]; qv = q_vals[k] - diff = np.max(np.abs(rv - qv)) - scale = np.max([np.max(np.abs(rv)), np.max(np.abs(qv)), eps]) - eta = diff / scale - metrics[k] = {"diff": diff, "scale": scale, "eta": eta} - scales.append(scale) - S_max = max(scales) if scales else 1.0 - cancel_spurious = [k for k in range(1,P) - if metrics[k]["eta"] < tau_cancel and metrics[k]["scale"] > tau_small * S_max] - small_spurious = [k for k in range(1,P) - if metrics[k]["scale"] <= tau_small * S_max] - return { - "cancel_spurious": cancel_spurious, - "small_spurious": small_spurious, - "metrics": metrics, - "S_max": S_max, - "C_raw": C_raw, - "Ct_raw": Ct_raw - } - -# ---------- 9. 统一入口 ---------- -@dataclass -class OPVFResult: - C_ml: np.ndarray - Ct_ml: np.ndarray - C_raw: np.ndarray - Ct_raw: np.ndarray - T: np.ndarray - Qμ: np.ndarray - Rμ: np.ndarray - exps: List[Tuple[int]] - poles: List[complex] - spurious: Dict - -def opvf_from_H(H_list: List[np.ndarray], - freqs_hz: np.ndarray, - param_matrix: np.ndarray, - total_degree: int, - poles: List[complex], - do_t1: bool = True) -> OPVFResult: - """ - H_list: 长度 Ns; 每项 shape (Nf,) - freqs_hz: (Nf,) - param_matrix: (Ns,d) 归一化前参数值 (内部不自动标准化以保持可控) - poles: 初始稳定极点列表 - """ - # 频率基 - Phi_ml = muntz_laguerre_basis(freqs_hz, poles) # (Nf,P) - G_raw = raw_partial_fraction_basis(freqs_hz, poles) # (Nf,P) - # 简单权 (ω 权): w = Δf (因 (1/2π)∫ φφ* dω => ∑ Δf) - w = _trap_weights(freqs_hz) - T = compute_transform_matrix(Phi_ml, G_raw, w) - - # 参数基 - Qμ, Rμ, exps = generate_param_basis(param_matrix, total_degree) - - # t=0 - C_ml = fit_t0(H_list, Phi_ml, Qμ) - - if do_t1: - C_ml, Ct_ml = fit_t1_SK(H_list, Phi_ml, Qμ, C_ml, max_iter=1) - else: - # D=1 - P,V = C_ml.shape - Ct_ml = np.zeros((P,V), dtype=complex) - Ct_ml[0,0]=1.0 - - # 假极点检测 - spurious = detect_spurious_poles(C_ml, Ct_ml, T, Qμ) - - return OPVFResult( - C_ml=C_ml, - Ct_ml=Ct_ml, - C_raw=spurious["C_raw"], - Ct_raw=spurious["Ct_raw"], - T=T, - Qμ=Qμ, - Rμ=Rμ, - exps=exps, - poles=poles, - spurious=spurious - ) - -# ---------- 辅助: 梯形权 ---------- -def _trap_weights(f: np.ndarray): - if len(f)==1: return np.ones(1) - df = np.diff(f) - w = np.zeros_like(f) - w[0]=0.5*df[0]; w[-1]=0.5*df[-1] - if len(f)>2: - w[1:-1]=0.5*(df[:-1]+df[1:]) - return w - -# ---------- 简单测试占位 ---------- -if __name__ == "__main__": - # 虚构数据: 2 个参数样本 (Ns=2), 频率 200 点 - freqs = np.linspace(1e8, 5e9, 200) - # 真正模型 (示例): H = Σ_k R_k/(s - p_k) - true_poles = [-0.5e3 + 1.2e9j, -0.5e3 - 1.2e9j] - s = 1j*2*np.pi*freqs - def synth(Rs): - return Rs[0]/(s - true_poles[0]) + Rs[1]/(s - true_poles[1]) - H_list = [synth([0.8+0.1j, 1.2-0.2j]), synth([0.9+0.05j, 1.1+0.1j])] - params = np.array([[0.0],[1.0]]) # 1 维参数 - # 给一个冗余极点集合 (含真实 + 额外) - start_poles = generate_starting_poles(n_pairs=10, beta_min=5e8, beta_max=2.0e9, alpha_scale=0.000001) - res = opvf_from_H(H_list, freqs, params, total_degree=1, poles=start_poles) - print("C_ml shape:", res.C_ml.shape) - print("假极点(cancel):", res.spurious["cancel_spurious"]) - print("假极点(small):", res.spurious["small_spurious"]) \ No newline at end of file diff --git a/img_formula_70.png b/img_formula_70.png deleted file mode 100644 index 5b8eed476491c89904142c3168b9c62eb8fcacf1..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 49695 zcmb@uby$^Kv^|W0D1saV2?dS{N(d;@jfl8GIwcQC2uOEh;SmLNQ%ZM-(hZ7=bazT^ zx;uY!0q5Lv@AvP=^C+9W-&pU8ImZ}tto=q-TAY}Wk`NCMj~Mgdo;)7jv2pnK=gAZB z&S&xm0`QC1MoiU4!Q8;cUdKuwPfEwe;;FgKQzP9!?ewjzjm*vNu-v}I!g=%0CpI<~ z*1W8&rvDznVs2%~nv(je0%bl&-X>)tKXU-9Z;kO^Qu;W8R z@bG%8EBg*!;T`4Oe~0>+#~Jj(gW%7|W9Y@;k(1C5@Z<3BAs+w#9im~3s}Ex%62)LJ zvn+A4k1KrwbjW8}j^TM6h1+^as>#b^i7G65x3{+yvkat&DcJ5gug%)pJHF4hpVR8I zs$pi-E-n66QPGwlLo?I)rJ_%t$7ORO!MLw5=lS#JibZy&KPoGymM2>i`me$RcDWxR&INiA=xEZhOX11iy zRFJk?(o{y7Rv4dsd-=0tIXV?5q(V5Q*81!t;$GOM{jk=S+?YcS30yi+SAlmT4aQhOnY*%)0N4O zvJ9Li;)0g9mm0{NBi!8F7MHB++@FOz44uk((sWnbeJAH}#P8p~opv{HE5EM!Cd$R! zXK|#)E0242Wfl+bLNxR`Hg={bk1uU+*jI@CuW~%4R8=_<5fS5-cyGa-WrYMz9v&Wv z2M=T`Jc&s!T!@tjV4Rpt$!H(?@s`TelZYa^u~8N}UsS+pS@yw$zl~cHqb!GiPE3C_ zy&)ai2-i$cOr*y@ae|SPQ<=4NUH6%tr_Ph_f1g#$w}=YnFo`cO7v?bSR(yMr*~&QR z!-uF3*TflldDWghdzLI8=T&Gwm!4zN!JD(SP)XC4Whf(d?_Q!tfz@PdlAqi9@7JeG z?bG=QT)#fYH}1(zyA~=bA9LTkHOEB7b^a>>S&ocK*o#u5mg}b7*%;LG3WM(OzJ750 z#ha?KLpVG<4YhF*0=&pS`0Qu>X_XRXzkZp`8*pn;Qqhcq#!Y@Q43~q45?miVwldSB zLgTvp!K`rNI%XA*Tj$$x;!B(lXS#C~oSX{Xb~g;uG>eq`itO#DS-H4W1e{k}NmVqW zs{C%4G)0T~ECp41U%F8xsImR?l6h;UevM6oz}k3b4NHrPwrkvt2p4&n+if+hWz&jl zYXj~XBO`hO>w_-{PoHja-(8|vnH$h9vY)dfA}HJ2-OBmv>rrp=+jpgc*_(b>yk6N^ z8<>KVRHoX~=r98=OWsZ^lW}k@jfI7U#f62Y>8`Am4EG($0oS##Zs?rUwL45qk`xpa zV#ewPR{9uhok2@R*=};tTxpSTh}}&01k7CJ`}aM+ef#F6B&h>4-rvv9`l}tx3qoSz z0NOP5JZbFZYpjc{iE^Sj#%)&oXYhnSe1iV`X3!9BD!~?X^N}t!ug$OCEJNP{Xzd*~ zwq_zCBKEcEuHh#=(1a|GOOpZLrXT71m_07zL3I< ziIND6WM@Eg`iV}SyhTFAs}~$^#G;Yk=&$EPBh+uSb>!HIkDfXodvZ;!Q~mT#96!#T zwlx~0{SY%|l2J0IN_O_F&KVlP?a-XAr<~{>2ENn+4@2t4RNWUxFj;1Og$8T0eGOaP z@(K!YJLoB@sl|R=bnoJ2Xg=KJ5?6Y;^}`U(I)A-gn1i^qxLeSxV_kdrFb5Sr7+>Od z`20A8Ov~Y`zus(rNm03VwrNj$@sJmjm*yqB!RGrU5L~f^wo_x-m#L`I*7}`gpu@1V zR6JJq&`|yKzLfCp9Z^F6Ahm!Ky`Z2L%mbymU=F2H7ss8|9`o)(+Z0O9rw;u%s@u9R zhu!-BSZNU2(~u0j*<_m67Y|G7ZfR+$6Rd@>T<7hDN;EfIH%5b`x5nGkLMA&h=2Xp^ zqF+In3^Nd5g@Ea{y&w(?UWUWCm0q?5LRXIMl(ILqfZA&k8at;{la6#JSStEiT_eK^ z%XdbdmdAgye=nn#4&$Z65E2pHoLBAVnRPUkn2Zb3YHo@VXYFL!+p>+?Y;I}Mlt)i2 zB20iM+?;LHa^d{>nCxtJco@l=AXa-qEf^#}qhPPMZxu8&G@`b+9p)pA+f$PuUL?sx z2ztGFBbQ}R*9xn=&k+}946RmhrlFx}OHpYs(tP;vU97uj=UjJ=Bx?DscSUfTPvugI zW|7^qLPo{J!3%;nWWrOwz93*?WR!sLTaaEf%iCw$q3NTRYu3kzq2zb4_-!6IE&6SV zWUuV~WLruKtP?RaR?OG0U$N5_!wrF-K25a7$Fh`dFY20GM{Q~-{rfzt*1=2lbhABq z8R?~)#>ABD4??)im7hFGe99R@yOr7p;g8;i2Ht;9t7%=w(Ov3VB)Bnh-gRfC!)Yol zZ+WqX)#T?p3K&V%q*5Pi&Gy<_uBcw-yxwRZD!xewlkFYa?rAWML?HIdTUlj7(?%WA zSc~VMvlLBd-dzi_r?!&|kk1L8od2nZn&_r()RMy`PMi|GN~|Lr3U_` zcQN)|0|RB{7VF;>AqR-QKGW71C2G69Zi$|VlFu$xK{&7LY4?}MuVXqiu(?mWu2u4BNY`D zZh5lDQ1aS5`OrPK^Ko!^IA*LaWC|DI&eAwmN#m|q;%;iaq;MJyP+_+k zc>tMN3?<=9=EYI3shtwbEoR`+7}(h4u{P`TI+gDxuqAl~A}>yK!^C@r`^0OTyt@>! z_aKv?B}pNn8nQfGEPz(izRxz8$2vYl6U|}` zf?KI}v%N`y%<7YH<#tG4(U7-!rB13`fUJ_NloEJ?h&&QLW@KkqM4<{6+e)(}D<*SI zqim&J9sUdo<}^)?5`CR)P#4@?WS;>C$Is7O?993EF|{vg!SiExf``mYH{#RFcAZea zEUv9}_%kRd<`^pp?5+<>n?lk~H80!AEnoTlc1X8gBxpL}NUy`Ncj0QciEWQ*o)JK~ zU>HtM;!9&6Zvd{nYCUy%5f3l_QU7R}a)y@X;FsQBjV!mFO@p1ymAaM$nUHv+;GR5l zuDQY~mMpa#7Ip@>P zMwj2O$U<4383eD^q>|cgWb@iiaINCLkjJ%k%Z?kf^ArF=DRX zzu?3i;UWMBy{XD+-S}%No%-z(LlA&CS&OFG?%cZ7m%-h z*)3og^Rf{ARUiwX*{Jn-W_M>xslaNKEs+q0XpO_k$5AjY9Dk29Wj3t z`c5uQEyt@QOy&6U`1lh5dMV|bD5|x@OBD^S>gecTr_FU89P(?X_Q{e%+48cDn3xz3 zn}MaHn6P$!y|k5?p+FHMEuIbAa`YE zW%X&ExLFTj4q!xFu6qT$R95S%)=xlzaS+=Wcz7PUZZ7i}7JYMSjERo+4-EWCYl#{c zR8?QEY;2qi7b%i_Mn(_N0v5jVI78dS#Dv%1f9H=solp`O&oXRK4q#NWfkYUQr}@`k zkw=dmt4U-;&oK)GYjth_*p)W)TNQA%_}_OaZc8aZ9&d#_PKSZiZ@*#fAcgWyuz+*! z2@>j?LP8mikI?B7gNSW8F2*9x0Sm35=(4JYg~eyu35Z^ZJk-1{a;|F-NAYL2ugDfq z4a{3~CiPW7=-3N?{4qEtG042=HeOX#Ww}vW5}&fsZITf{7A8#sN0qtFlgqWm4|c*Z?)pMHq=DM`u9AbX@* zKJ;8%oycrD5u8n*3rmcNk57HBc*zF_KLG%>q+}3F8*qg_O-hMV0^{@4w{P=WT@RIf zRIMC&j`YTlpM^ziaH^`NRc4EZxP&xITLPcm^vBYzSw2V=7|lW(1U_2RHTk?LIV8CH z9b#zF8Z_=NE`SMjfs`~3EsaWF8vE&M3=EjjnxGED2-o<|PSpS=73pIqh#G+jt<>|) z(F3_7WmIg0#XbFTMhhDq6BF?9W8$-81dOz__c%=2TLC%NtIvfCINQ`P=OwyqEGDad ziQZhEh=DB5OD%)J$U*$QaPC|*goCN%l=S9w%_5Dhc*i5CQ+53cc%$xNEEiSYLx)1Q zp1agwX%8w3-sj8_y(D*cqlyy@cUpE0_jcnbn?*%MNhv9lU!5jnWMWEb7-tHt{JU+A-Gk=h)3%eoBdeH=hY8$uZP zr!R7wu6;Yrny6nBsFr6I7VbDkA+Y%Cnp&}gWt8p6$Vj5NFEyfdfLan@I7J{ZPj+P` zz(Ta$TsA^96FRyX(J|<`Nf>~~20M3eeW%j8#+CcV#q}p+2d5Sf_I}>_J!7j7tH&Yv zMeJzc9^HMeM=sy5A12jO%5N6hnxUP<9~Ti30a}L*#9drD0g(+wXK}(IBW8Xl6iH3)pXjeD4 zeQMvPq%r|C`2+2b@88E?)3`N4YDtCVMwVy`057xaDT2S3qad7iA-P6gUcNP6>QY=> z+yvxQi`J%m%f}t1I$F=EZ5+JU9TW<3k5^`z&ZzyUGySID4i6^Frk=ZfzJBRre72)? zlA!6{?oN*RfS@Y~J*eBz^$KlG=;-J$*ua~QqG2u|M1eXu%)URV;x)kKo z>(*5G`1nb{=z#OIK(d5b{Re1LCV;#Oe|uhPjTSrYy4)nLmSx~Q%)hX`$GW<|)4Ta$=I&91Oy1QqNI+B%`(E4XmPL&0!#%_4V zS}e@l<4lI7op8c_B4%mIh=3=Oe|u40_We%7h~7keh{Oy735{T?`_7cM6HLWrfQ*bB z97-_xqpkyEgRx8m_+4F9)dZA1o{GJ_MxTCSUEwnxV0v!axz?1H5qs{pW5~O2Y z=w1bGtTI|1fW$vn2<;M5Qc4mMUes>uItb;$pKTy^AqGobsW#j8`q#7fW!D~PGVZ)lF3qN4vSK$e|>y`6~XFJ5ryS6@KE z5FlSOM~6ntsq@#zfn!=#EeQc0MGzXH69BPNpV-6@Ai06TJ(n#NlKRL`FLH{KO@=%A$n?HZ#PfUc;&jMQMf?Q*SH%30-9&fo|TZ4 zQ~nB=%>KLBUQe!qWu!MqzZgV~C@0iYo$btAJ@H?P3jz9#Nt6HG_F z$30+}9NUGxGCNNh`U2^^<<%%_8IoAyF7Lw;w zMfppW;>iDbO?`0QeVBEFhu1)_|GH=ohvO-u*OB|P!Xw~p<8O`wsMeoujjN~Rv5GCK z5oWjotLN5RJ~+Dg$NMXydgS+?^}u?Cj4=L$4CE3OEyo|OufD}?_iEw|GVZsnJY;#g zO95boof;Vh1p|NwE~S@9W+^^_Lv1e)Iey|qGJtpKJO~E7p858(y^kPcG{^V~PLAum z&OuH2_yeKFAS}piZOKZKKonXaXXd!?x!DK3Jg56NU-1b6Bz@C+ahi6IgMO>OTMP%b zAn9dfHNtiE(_K+fsv01vR!+shpB7hElvGte0stDzTF*MsD?d#qmeI*uSU`cI4P=V+6HXz?5qG41f5q3#R0{ryZH7&o^0Q5f-a1g z-}(;3aUfWVkOfF7C?;mDWbbO-{_ffwi}CXELLrE1YO+0T>ZjM>nS)_br=_DKEfl1q z`(d=<#gyI%+PM5Ce$+$2s+&m5Mjbg8VvfLID}<4fjj<5i(};k8IDpLwmr{ML(Q3p~ z{RD6uLckM9_{N=|qe?d?nFQv`k0F5z!G{aoZr@&>0(fcB8yy{8ov-VE!oU)N`oHu6 zNSktoK+a7APOSyly&2*{g7aL7${mAR`s^JD2Y_D3@@F7zykoj=^v~1;unxjRA%w=E z>~s*7x=Wnxqh^jBJ4Q-PJ!|t~<)Az9n3lXe+#uT-0k4*x0`(&r)^C2&)%;`oBZQY) z%XD-T<`@FO6A%kJa!q^K4eKx0h4Lss8nU-Hq8I&V;lVKEnf55txNmUik6pIW=&ld{stMr$`kJ2k)TvWW!`@u%wv&peZ~=%Pg+JTn zaXI@S1Eb~03Z_OBV?3JK>{9(dju*4sZYnKJwSkpLBDDQE+-Z^ll^M{A>dYd%MI?&w z1`NW5%a>DubQpo2?^CHDFRw@L31o0gHP3CwiNZ4q03ahBofybBN@{AsL>2FK$nnFk zHK5sr&?nb`OBT$@6eLT8JI^+qm*M~YB zMkNgMO%aF#gw+(3XI3WLQ1t^=IPecvD$s$yX!egkfBEvb*!}yA{QMfAwYx_pl#BuC zYKEJjItBu-psK1s3L0Jj!6?aou0L5e>N$X80JI9O3qQ!nxEa&Mg2U?==0Ta(Ssi$Yi29 z4iyBXxcX>818-#ky*WnQw;hTCqc~gwWZ{XVgz)B+jMC(}(k*F#LMVAP>8_BLMmf|G6RtR z4p3m&jhbm+d3s))aJ8oM%&%a&!1+CFKY!%k)d<~bxcxYGVK{F3sZnpfdwWm-Ky0Li zV_QG#h$49s*ozX7Xev{)(hUPtcftAtY54J$Y$~Tis|5gs+t(d}iQOTe(BQhK#1t>Bs_Y;Rn;MHr_U;u}uQ#Lmtvr`%9ep+%CV5 z1-xcFSpE#o3J0y9l!}VEiu&VH&0Zl;M@EcKWN6;(bVpaFz9g`U`Y*jvq6AbzZW;i~ z8Dytmc>F=Y{P>X^Mu^u8*iYKyUpv6iQWu6R5yg^saVd5K2LIN@P{I=;O*zkga+G@W zX7ALp zf7n(~6HyfeX?&(a_ETPs_p1mnHB7nd@9K#Ix00FzWy5<rXQbcl@cJaR;P=wO(uA`@Ym^D8$0TyLSH_y~S?9KBHzf{?&*X< zSrtUvY?}#rR9}WQUvl5)GMdi}z#)JhwLxj~vd~k1I;Y0pqz9=gPJg=dIWDyFN5G7J zyBFh?DDJbk2?Ul+&>;~iLMR(JseSiuHB&|_Vy5cmB}yMY zdc|68<@)enu zMc33cXn25|6HU?!?&tM8rI-~j{Iu3OCyj-oWm?f}p(Gk8-N%pH^31im~{>|Fo);Bhipvuk*A{4;H z(P2>0qc zhAVr2tVKrtsc&E~38EdFl!B*D{ud-T<*N$|d%&3dx@GVGv8>EJX*m=Jq=0E;PK1P* z*!K7;sOB3R3%8*g@>@67*Xil$fA#6(;eMPVpC>dVSmqCXz0YMBfy)qXyw1mwjQ_89 zoEt^&-o3mwanO&Xs1K}(Ef2V=hteO=#fv2W{7DMTfOOHRac6t`0y+7(s>}ffOZgME z2Y(ry3L9HneC-DvSiXNp$0GvJpov1<0d4wBm~nFaz{nz5FyhzlHxHn}S zlyxgj8^Q(jHEwG#u!c(34(*RZAt1%&?dl-b49F#5QAum8L^zt9T=l0?RQ3q5dC=j6 z*3@cx;g$^d)xGc_>>gqI5@ z&Eg?c-+%p|TY+S>c+d}n#3#o;^B*kigVzUS=HT_g=sO%c%6s#W9UZ(rcxs1ZWv980 z@b`2Hj0+0xZO(7jDG#!w9n5LhL;puT>W4wX?ECWbFHGBiKmjBFJr$!w| zLjoJf*){*`t^Y0!IASdlS5`}T7-OEp!`7T0N1mgFg{Xt4raKE4aun@dV)?S~d%(kE zLJtdT>c#Wtv~Lddfqxj3k)HV=>C_JwmzTYka?}t$}#<@E+{nJot(#9oimecgf z#SK#Zf7nk%Ed8^VNX8F_bBL{F4PkDINV;5znJBIbG0Jd46 zk(HZ>Q=sbhWBt$1i~Fss9RKD!?1!H-M_eHf88=3XJX_%m@$Q$=klpq8XDM1c3!4L; z02#OF=oC~CBf)7#>g*00oWBk^QbA59zewOy5pI|gxv{Ujd?pmKC%a`(OAh63RS*sW z>I@VBIB}Zyr|#{pm2FOEHF)J@05O^J=IWP)EGm(t<`7)#D`eM8qE@OQAsx86V!a%B z=Q|>59!mKo1SNoaH{pR2ph!;-H0ANncTIpW()S>RP6B-Vux1T4{ywO&JHeV)in{ws zK9!vM8qcFgYXfttL=PqXHrIl%euqQCV>Txe6br1(cwy)VU3Vpv{%ipmrl{wcQ3|+=wMOUuZxHd$9lojd0rV18Q`OC(Ck<#Rax!>iObDjz{eZ%@}u0TJy&(hfKg zE`Y$84|Sjbt5@`->bVz)#^V`G`yU>Jv>HoE5~o*$q#zJ|=%(he8Y%Xj>&g;GIQ=54 z*0jy1v?U4-V}EMr>D!f7PENd6`n+^?{g~rsH+Pw3KrWfMmmE> zQ9fN{8Q>(JErj#H+Sdl4!gjn-1dB97AdD)Y*%l3!wx>RdL9h;tf0inrt zvFZkXyF2`!teAwo#P~d5B;q)WJAH#tY}ls*o$y(YbL07Q=a88ZT9oLzwTg`esY4DK z3qs6G9p$!aFLe{m1qr#}NJb&39P|*Z6dI^7Q*xQ5qCyLB%s^TWvtHE!_q`J6MI06# zK7JUUDDwzxFvbI=CH?Muu28c*&8u=KKzp?13q-A;05~{-o;^N{p zx87L`hKGcxd$yGuABJJ(Uh2k@xwj569twkmz zB=9kZ#3unuH<|AIf&{LRok zdAN-A%_3;ifmQQ@L+05Ec#0;QDXJU>N2OCo6TIB0juNb(w zI~P;|@ulYv5%z&JH_J$J(J~z~bO;ML%ub_%aUe3{PC@$ri<^(7y5Rbo!=&1l?|N}< zeci!62=Dap0Z{1w#(bm%MobM93NV4~hgV*VMmU1*I{Kd}|84NO#;xUv;g1*6LHU}c zKY|KQus-n6`2u^N?L6~iX#brIZfBgrAtzArJIHc|)ET)#1))#zZ0ete4b+vcWniR~ zc<8^ma@eR*7&Y^*Ofkfulv7oIXZN#9fvyC&l}qV_Z%}Az%Tznh#T_nFEQC?P&AO?#9D z(SgItp5b#av$_rAfNmu1C#a3xIY>u0E5X$bTFX86y`8+G8V}MR-^WKj6@BGh(%~dx z;0!!eS$%qrZN_5!mRl8M7(h{j5U31MriGY0Vcj8&+d3~)j|3Ob@T2k&>ss{+elLLt zkBvo=dYU1YTSK{t1$kh zW6}?Xe>QKf_^kZT5?%%(T~Ksf)gP#Z0tdUMJ1>ve`yeX!kz!zLC2o4jHkf?yv z2FxfDnZJp&@q4eBI}|Be{09#RjDxo~6cq_q#+xd>q=PyJdR*Yp(9kkiw;3P>Bc7?{ zvhd^Xnp}!9vs5Uz9B5k53r@H&`&KAuA%88<1m={W=WKyOhjL;0Bt$8$r=5|YJ6}-M zPLQTD?ffi)btT^hX%Je7lptmyn1%@O1I!Q;Uem>0`Ul{qSZ@8X1~bP9WL|*#cBwT` zZv$x%!q_!TGn}ggaQW&C@-9QzX?b)0Mk^5gyFGgP`V0aBeG@09cn1C@B;F9=6Qn}o zq2etEr$M!knYPQh&Dz>}rOPk^gSDIPi~>b=zIdzav(5|DSx|5TZ3!72l|jTJH01;T zG#xE%{_3m9+A#C~bl#C2p!|y9D?yRn6xJ0MF%c6$SIELci=3fky)sXqrlJu88p5pl zSs3zjI}io+Cl6@;0|E-6K%lYgC98Z;?D61&2^ACzPO3gw8+6^`hPfgn5THG$LL}0X~R58 zYCp5a3*8IKMJ4VK4#R?f)Jx+ya<(l^U4?*zS`liA%RT00MnJ5QJPnm^uWVm<^hu~e z%0q_~+RY>*TnGx1JZPnoc@Ft)8=m=+d&|*>%b2jY|3xW6H_)TW8U>jyE-pwELXs0Q zZXv-BLaC+&4Y=y;C=>~0tV^kFs@$lZ?CdxX!^_thfpI`VK3&~RD(R_^+A;h7PtaS6 zve+H+m%?7K9Hv_}xLJSdqdPlLP)%0POYhDz&wxe9^yi<#fUJ?fvsV)wov%FeadjjA zz-5Hh1)Q?4J%PEA2om~a=a(3y+QSkmDByL-JpC_s4}h?{2FuH`3|27Gjci7Kokb4* z+kSkqAOCwCb0qlQ`^1a^^ryP!QxRToqq~3c&R_qIz&i_h+5J95Qha_v`Ybr120Gf z6Txu?`fP7l+}ERvp7~^)BB6);&cxD0xS^r4QruJL7rPXQQ-xrA64oNwEI8_s=R|vd?!+ur+Hl=jIW9x!T1)Slun)w8pnGvtUg$~iJrMhs$?Z=+R z9qDPH0!_dyZQEauMc{W~OD9r(zehPv#!GNmKJvg>s|%M+r5C~aC<>>(7$g))rvEqJ zodN^|=WPH3pe;d8#tzSlDSG3g{~hv#rH_w#pAFH^{u75Bc!3%x(TnyTf*QKxHmdS+ z|43Zr5mbvq@*HeqA)a{Jy^hD$wFBYC-8O^s##*s z(dkbz(;0{Uq@P_A(O02G$6PpegeuVH>!aShr*dT8LH=|MvTg_Vt&p13`{)SzGw@gw zcSUbWOcc+-y^swr%X#NI`we&8DNS9*hPC3_PjsiO|zuwH2Rl)Lk9%VX+hY6#5pMY0Ue z)}x1aU;%Lev*Fx$T7W=PUQaJ}co*bW$&%~nbox`EcOg?DREkJ1UW|vOiYklUwNLxB zm{}tU@bRlE6Xty_<0JxLtZ5JSnKj@*+d)$Va3zvxn6-*kVUH1V?S?r$`mKrPyVIMT zJ80h*5Dv_NSu(tVsa`}o`^kqED_FF^J))szGahV`ZmFO|FXbpy8p;|eDG$ZvWvI?g zfdi0Ds--CgsZx01Kbf2w8v+9ZuVavvEQCqKsOt1Mj(3k=29Br2BuIzp-rL^_=Q~C! zFunoM0-RCus7MZ#J>@4Mu#EJ9)dse-wA3W?=X(H)q#R<|Ip0XPb!)A$u~9ya7~O{@ zuwx|>oGuz*`D`tAT&nLf_;0QH^Id;MJG)#kh0>9-M$$k)QJgvor+HxKlRX4=zn^qn z?J`sY!~|v@+j={YkL97tIh&ZlI(H~R%wI&6eki~r%L!POluDdzD^0=8HOoRA=YjT(fnDQ+%5r2#{KV`zSx@o< z@$V9+ z-1UPQif$GcTDow5X61|f(KEPjEHT)%?D|scLAR^ zlGKGW_u2TpWjYSWc)W#TygumqP?t!U3MRvpxfweTxcdBvNKv-cM~+cESzUPk=SIOk z9$H`^($zdBAOpWIuYpxzCbJtIQBDdy`VzSe0ZL=N9CsZBLH4%oJsbsgIRBeignPIQ zz#T+<@(+A@e3aN? z%}-HuU&E6I8=GnoQtTm0#lc1#1+*ClJm63v9p^-skuV=af$M^&_j~I z$v;a+&W1eTc7`Lt`hDvr&83u1a6cGeOlYaz zE2bSkY$j=OG{Fuii)v+VCG8QX%mGHCmr!_}>Hj9H3;Q}?H{3p@E;LA3SeSGlGfR)P zk2^>F356lBhj@AHLRBLL3`Nk{0kjieDybEq(##I(Fqeua7HJ;Cr8P&L%G@Yno6BQM z<-a0peFY2051*IN^HYa_<0(lQ8I^nYUcysf07KKP)q5F8E!mDs21t7Va*MEFS1BJD zgTb{ei3JTK>I@Zc{?h@}2o?dlxx}>ph2+|P z9QXn5(@MsG-!mITV7oy3BVrx13?}Z$zHiN>>kviHGXWP<{Cb6?G*5FYiws_h@LSSTKt% zPwF<~vs}2$tv0Ot!h$FBg(Nc@#>`NJXK`miCvnpp<>i+K9OfizSuGnK`!Gy7`p{y( zh%?^9MHst`?OI_j(K$a5vYtf_L_B{|IHe(wxw_c_^9Xx%H09h(b)*3 z!RkrKKTWNPIE{|2f9-L`*E-Tnlg=nRRv7cBS_y3QD!frVf?6aGlI;y{i#2Z zjV4RWg~hqcw4mxhalcE;w^&P9sb?@*1KZ+1&Y>%`ic*@I zR}#&k&d_-0ab}JCPN+*|M%mm^34px*`McE4|=q6IIHm#`WI#!48r>pyT(ls|^q`(VR%N1?hlTbrxI-JRi2 zkKg-U;q}hLC)4kJWFNM$iYJuSul2hooUc^o#s-O5R6V!V0JApc{KPj_oVM-5PW02B zCF);QS5nc^p=UeFD>3k?t6OMeMMd+vQQfDu57^QVt@jrH#{6a(uO*rbfNR{L+o4Ue zFMYiixEwkyFd3#lXS&BX_cb)syjRiikMP<@Q`)!ddjoH1md?!iDpv?zq%1$4|8~@g+_S=q*Jn)S2i%(j~e7&MgYZ|!xI<|t&A^1-og4JSmIt*z;FXPBo*7Jyn=fUg@ z+2U@~5fvksC{yznGkEU0>UQdln)O9}neaWk2vNb~&dG5WVPPTP5xp>KK4lXwQ+Wxq z^wI#==B(jY?8Jcr?}cbVVzil9`A8;yMMq$d2lrZ^WXds55mQ7rgUv3xP<^Uh=$U}V z>e65Zt!c-TG==DNxzO|$2iJ_I&4+I&6S;|79K?%cg@o``QsS*Vqat4r=zsqi?tF5$ zt7Ltn=j=MP;xD=Ve%`xw=Z->%w!7vfZAUSoxxSP(&XD%W+NZ;>=KWg5nyIj6&7Z33 z7T32VrGhoO1;)ez-059*|8n-352rKAnAUO4i!Gz2;$AY`?HRyoEaB5IGc45JApj|X zeYaTsxFn;>uLztv87;$L?VaPB>-3-R;^UKh9&Ks9pqdqbv?o_-FqngIeMtsiIa$x# ze*N>&tJg@MkQ4ITJa|TYi80bpX2WI2_P6LJk2&l($BTREF?%#fS6FCuwn#$iVe42J z-wJ&7#3nPJdYQg~afp%DyvsNv)=8pJYt_1@2Iqdu{3!R`f+>y2gpzVA!2~+0W}#<0 zf5*mb?D9ML*xoqZryf7u#3Q*%Khw2sbA5a&qH*pI;6*y>XD`GkQA)-8%f5|qW451pKq6qognV&l z%X~vQeY;`DX=S7*?yZZ0{78{&k;^6a_j9m!#AR#6FC!xZjEtkh6-w+dh}D!w`tXaH zBT<3QeP<6K;knN2V);rZ+nza8~Nm z7+FfJdbqO<>=X-SNzOM5i1t5tUVQerA*=I6?JDX{Q_hXxc!9x7o~~UT9i{aKepg^O zW~EzAU4U;!hO4Goy36+V$h9hNS2e+B&lcbO{a(hYa5c~1RvR{B&WXz)0@-jJ!jLHoV5?#+2NMrD7wjPQ3y1iIw}+pcW}%a7zd?UK(ZSq%mgwEfqa)IS>U+_{5xf`IXR zJTwoLY6pC%Hm5qjyrpcvOuMDu7#v=Z#X&ji^yoM|LV%G8^+l`>E(tmByf8-XJ7tuPVsQ zaAP2X=SQ-j%);v}=Nl$Fs&a)x`g=Aw_r4hYy*IR9nMnh}BsX#*(ro(Qo_^(u$^^eY zo$dpL%&WJ(ZHqB=r(pHB*NT{FjqepajkMUV=8Ox*R&O_!a8tk;T3{*!#3A#*(z(@62%qmQ{CoRq;`?v@28#to$Rmpt;-Hm?%RtoQX=Hm6C&n zkaM4qyN}8ZZ)KmOqZv2pOtyJ?#`+qVhw7By<)#;d3xRL0VLLCYV<#(=?Cv>zk}E-r zzN7`%;mYJJ_tbg2phXiQ^@(D#c%y1HSBCb@rl0Ri*n}!R`!X&rd@vXGaIv-;2?=2m z(uS?JvBe+6vD;SGqnF6(2j(p{4D^0+-BBn#a%9nOx(EMwXI@o(yK>pW;wfU~;g$t9 zADSQ;+4sA@9Muebk}i&3?-9q-H>keBuOY}`Qr(#+BM9N()R|o03)NiFzqlRc1xIUY zZt$x8OxIlGaZ(h#a%I)-hBWs;h9X9g_)>Leu9V)4qT2<;74JjojMHeaaIYTf06~ zez29jsW9iA?UdehA5#G9oI@?tLQ2bG#PwT~=``I6-f@`dZEqwjJtk%Bbaw&UH7vMy zb~Ae0%YHWTVV3b>i+xIek8#oywU2)exbuOrS7owO?x+0#CMijNw0iYJhGF+;AJ=4?A9!@Z z&-C+?9j+H@oouB$xr4j#c2hx?KlbWOr_KOw*jKVc6Sq7Zlro3B|0A{ZJmw=D0^?{O z1_;s|usH>NwUGkUV5ptOuLvyF-of>?w>US(#Kid(dM zsyDhJ`(3lxE~6-YWkINInmxr&p%_NOQ`@{2gS|8Lw&AXBV#JH-o#Iy(@BeX}Ub2hO z=H1Qw3!kvh{ga;JGYUVFP-x_N;rb}`PD5eRd$6z!|2^kNLQ_bh-12JXgBF%cTMy~mDwrX+^(pm?x=ed)z~EAHl@%>jS*+L)aO2C)MN{tBkfE z!uJ6om+=?3`N6Neu%#P)r^J%k+{BTR@GqK4$r)DN=~s*+dAj2tCuck`?I19Ye1Q3o zsXbiyKwuD8(P*0=!gS1CZcR?PVn=gf_}hw-l`0+tg0Ly5H8)r%h z8QC2x+}&$0u4Y$vfvD_;i>$)>3pI3eusRwSAHjWgAQU=h>y`q|Xp(*E{YN6%+ z>Ap$HIgZ&}DX);8*mowSw~iLZLP;?R_Cu+C(S5~hW1;a(T}EKVzjS9+9c>Sgooa z^7lFg)_!rZt6GxA^Omx{(eVo?vVHwCUckvxh|-}s@kKNn5nXUY_;96{%auM&p$GZD z34A0sD!owL^DWj4RS^wTt`WDtF&Pqjs^S zS}`3d2Lc(9@{T^u8p@@k*j^d=_@{*azM1Unwe}+#Gnxx7`nF~HJx`mSm?}rR_|t_P zC5_?c;9puZ*j-iUP*Y#>mZeO%?Pao&1SFGR>J6W}6zgY@Ulgzfx<3dLR%8S>V!ObJ zadZa0S7gZ@cR4RRJ3EuHz+rmt6+2zery|bE!mEATJ5pKtk}WQi3|kQ*BCA4vetp~D z@bbTl)9ftur7Kjj$HeSygyZ)5aqg+oR4wn`u}f8a9=LifjD~3%(EpFE>%qfUwQP*J z(zhiIip_^!MdT%}1c>puyKfEnd0NhviL%aja(i>qwOmeUif15IEgcl}+SH)#obAmJ zsJ}MJRAVQ$d0FFXUU3I?WUb_?z~Ihc$)GF!Ci9%@Nm+x1ZhZWJCiiTw8_djgH!->X zZ>KEZLZzjpqoZRR#yYHh0KU6v3|x+?k*?sM|Mbt^y?f@mNq)73-F9PHMkEuqIWyZd ztHrX?CkZ^3iayeC7cW$Ly?vrGn3MO^$fszN=}zgSjc;Pc%^sc-GDEj@iUjUES_^yI zJ=0z^C-}E3y{_)b{{Ed&p)^y!+cf!HsMlZLwV2Gz2p1sq^r9^dZbL(3ywdCW%+Rzu zfpaAGgiQwd5KGsMH}E}KVOq{J92hXyp}ObUYxh;%u=!y2Uh>Sr5!>4nS(BW}C2`mP z8(4IB)X9Xpe7!T?nf!6_Cue?^DyWpcK0{x>=2$tw_2z2>Q-{!l$65CaN=h@$xHyOO zH#k!h^4}#~o;lY-t!(xaoxX6{z%uVSR2k5wiAFG6+jQv%O$hC6_(BEW{p$?Bz`kJ*kDmr0 z?7{a#dSDria+eui1q0Ft{pt4EeYqNM9`B&+fyWOL|NciL=o^dRr~KF-|Mw7&h@#W; z#~p*7%gD&sThr|yp7Q@_d+&HI+dh7{t4nFprZO%?OQ?{op^}8`>`KX&oh_B7kwQo) z*(0)d3)x$C8QEmZ{JoDa*ZsSn=k+|lXT6?3+})S)JR9a_fvAP9yxLef^h^R zXMOMhRYAd+=~Lm||M};Gsl70QxepqL3LZ$bQ=a`?)8qu6EQQ7C;(`{WIpG~jdjN5GR6Cy#G}h73D5V;=Bmw48WxrN4mcO7IoIQmC zcShiTR{=h83gmJ23P+&&u#UZT0OM6T?$nItFW7(ZU}T;hGkE1SuEQSv+q>y4LDSwIJIrvnmm9S3nH0t40UV$W=KffO0$ZRrzPjW4ouox~!WQ8Ujv2dT^R1{hXJJRul^z+=>!^~a@2W1Y%F=b>|B-R1S#9wK1? z4-%P1aHteO^s|6VJlzz3IpOr*)HF3BQ_TeX8HomVA9UW18#njpasx3w-+f}S5+U?B zKqiFd31)m=ntQX`k_~*TJM5j*qtglw2naOvP8&9-$wGoB1H9&nH|+A1KXjS)VUu;ak@qeO7(m5{aZ!K{_TO3fYSDMW<}r}&r9 za8}1*??~8pp|U{0DFO8j@Uaf zuRpMNL8+zkG9ov?IIq!b0T57slHl;vtF+xeZSR`i7AEn0pU&yiF4_X3qOC+ytAE{H zFH)CQUq|4*am2?(%K6pTk|8irf?TT$bN)>M#r(1GT! z3_%|Pg>B@3ejvbv>KUlY5V^X9gaq3kcc_Uxoo}7v!L;uk1fh@4uQ)7VUv;<#u=;!j z+SonC+uK42(k8XVDgadDS9N2IW*+P3cXvM4rm0V8Ff1x8Jn-A3f(Ut_VF+}X41PP4 zbfCi4acTH!Sosl}DcLuWZn+iHE!3$(=o2^wQkH!+rmE%E5f-6HP~`69ei+SENZj;S zi`2rxtO;;OQRqoY_oE%9`a?Zb7|+IUGMjbc$p@KjL)YBP8(urk=mA^!xtHq^L7e9< z&d2+3YD8h_2|P#<3?mRG`>STMTd)d{2_hsC6Y))i*W|TKuDD^9|Y1G6eB`{W=sFI$aKpMf87fyn8UwKHnoLRUH#BN zbOZ^p5^9wV0!Hf_ZwiZ%^i6N)y8rWFh0gUXFmO+yzJ#0fXKm z+sF16xQ4!pkh|J{2)Yr(;R(pV9yofm(#<>M2+^>^$HYN$OCYveh)VrC;HNJjC4xBb zOYht-bYcaDQbD&Qi4DYo#y?Us_kiGTDjwv3RyusYaxojXh2%USD{%E?h9PN2PI=bSmAaUJrX1cz4nV@|&Ip zl9av5$!eO$JJ5W7H`F95_U7jWE& zb6VJG;ZHb3V`H#UGDC+XzkX#3LS@06XgC+xl}{BUUod+rI`>-)prBrNp`o-NW4{>> z3DM2)6>e{nd{NSzZq5s6H;{Mry`%TK&2V<{W&-hSgu}e_Y=mPCfe58g=0|^MG=!^p zcnJUq`mHn2$WVZbUTFIU%v9bB#29FJA4vT((F4iyj9neZj zVJfTbfBD^H7=-tyt}ZVx2eg>36vT;^ryd@B^YinifEV*OUVO!0c^KQ_&$scLF_p^Qbqx~h=GyZ5~ z{{n-t1ZE`wpUi`TS}{h z0}$fBfM3#s5gPu50O%svS>{Aw`8b{6uLBwFOhNE+ZR2BO+`uRQ8aJr`0EWSc$!X4% z8_FaXrKEVt8hwC`)`56XDd4j{fBMH(e-$}L4o>>5pQp`n@m~P%mWVUid>J1ZISMV~ zcHPg`g~WT0Hk%eO>YxL@zGXk$rVXrg$-{?h$Y=Ks$2{LJBiN`OMIsIzQc+Za1UdIV z|Ku4WS4aoatw0;6_ZYCe12qc(c^=WecC7?!{Q-)dzo1z%G@b@AW<8}^}*-kL~+nE_$*iMQcS`2JqGvIOCVKcQ2n!yF_! zT-6LsE6{~b3`JBE!RZI+|};oleMXr%GMWOfCQL(G%)apSs|WyxNZ`HNxQ%~g7M`% zapDBZB&jC@NWm9qbCBkrd3bn$-QMffIJTUw_*l1pSxf50L<2Hn5E0LID3-u=UM694 zk`xpef!ZR90r;gv)~+2zPK+3gjEaAR}C&kO}hC_ z3Znt?@SGAeKOZn#AP!zR`_SFJ4(H{N1eV#SSZi5rW<@+73`VSax*V~jL@f?f1r?$? zjzig^;SZ20NhW^O5}JR*hnxKXcT>XrTt^2Yrb_s+ya=bPeY^9PAcD@yO~C5I%z;tg zkt7P`4g}I>@`|gLvpLN{Kv|%3REzrn4Yz;%GDOq@q7FS=nh5Y)pOC$cjpGPG9Q449 z%va`nmNLb|!}v?0S*0z2D5$?aUYQ6NgSU@BfXM1XSfG~;NYD$AVylATu7v`C2HcTX zDA3RcAxsZRJ)hfLlr$<$ z7FsOOs*}4i4>YLvix;76g(XDipv)?f0twe8AOV3(evLtHP@c($JQc#zPL%e!nnnpbf@;$0JPu1{dMdLOypf3^yV|<1pjtoH+BgE7L z8WqY{995!B4&$=P#x1Jb*_8+*q6YebZj3@o5W?W9IF3yVd|eU};1FDhWW|ur;jAcw z1}p+TvURpW?h)1%QLiRSjAtKdnO`j@QaglXL8cx+bR3FIYQa(f8ABf890;j|utUx$ zHnBGZ!bef%cn$Q=iub@U=JyieN77$`QcyO;-H9lC$fp9koL^nH*EM!>Q5A;_x1eQJa{#6Z;o zHb(*354jK#b4FmNK$!d=8eTQhV8(`g*hPxQ5bj8kLdRv0l{9R#m=poJkG0I#UuQxQVMB!+Sq9Z44o z1+&2Rjt)2wF2AteDy%3|A=Qe!(@Vnx2)r_B0RYD1Y<@`GtySA1<@XaOsRNb^xWDJa z3lR3-AgabN)_iP!L`kqC0fd6Gr>J5gu~!XtMSN4Nhs3Jk#dYJRlv2HUL`m|xu{T(T z9iUQSsd@tD^8-F6tW<~ngDL=T*Nvg@kv9XY|90=*`T| z1`z$Pm=6txIw-W1LY>GjTqcO5!>DGLkekn&grHU!t2rpuoY8An{O}`G#M}2_MlX1! zU>EmvkOIZrpRqt&wZ*C2#a?LSWpw$dtxJNmaVAmHM*K{7I|lMD(#3FX3KkYg0BR0j zu7^`0L~^g16eBg@MWHitUJj&AXyP6Gd}JKzrTRn;Mk_Z5vb?qU9mt+E_NGJYN14o9 z!2u*D1P3MzBA{wG#G+FH&t-7N@rF4O&E}jR)F9fs`skYb%p@ty`u zdKgvqPAOLEGxO+r=DRaSBh7{u<)RmBsUNpE#JV^(rLjd~FOUKb(xV8NzQCz; z=zoU5nF__)!C7zLy*t{E`zBOiObMAu-kGa}Bx#^=IW5k%!JuN7D?yJt2|+4ULL?J| zxJkUakx@_qPgji;FARU)hwQr=t9hM1-zKO7#3*0Zc_>=qDQj(-F#jRxeD_0N=;h*R zGOuKQ=n7n~h26I&nJ~$*8NxtHB}p(OyOdD%K*HLsi}+m-(yIqhRnlw9I!@dn`XU|8 z;}~Z{=m!XBG9U&RS>v>yuBV&8QM^RV z92N-u_<2d}gKDIZAw-|8Vp2wUj`NJ)tY9#O!l}mN<9zzaF|}= z7ANsah!9W>viN9=!}ktvT;sh40L$15b|0<&+_-tOpGc@7hQkqP+G8|EP`Aj~-Mf%o z6Qx;FozZSVIs_`j#LA{$u~3i;(k2Ava3V}(-E8_mTLV8DN0b$Tblu&O7{le15#)La zFCqnH0DUepMqvQ~2wjh9xe3$&h*0GnW|?)-kw?f5La7;RW7%r8VnzR6L{En!k*P=* zL&=dSrz)V@Vs6Z>Tk#LuzI}b*##(QjEEzt7OtcCokMl(Cn|W)*ulq=nE&Yl=1fE&bGvYO~psi`5pk zZA9`K-iF*w287tkRIE+)1)K18;-O{Hg9YFkw^+Zz*_!CXAJ&5nB2rf%&0sy=0nALd z%Vm^QWyZg3mFYo2r)G9IH=YOpqQKZ~*amaQwtqhn$|j;-ghoT^iMx(maVFf-eNz7X zjsBc3ucE-&C|zuGbsku0cS{nB!h1p(EF8(BV6u#)BzqOqiETcH{Prr_EhyPe+9n9b zi7Xt6ZqK%I;%tt5zKFd_^dKkTv5}1_0{+d`;FMQbRCMUQQR~2h&52vDfCj=QB%GT- z&HxA!@DZ6gA~JxR3B|JYLL-W>m37Kqspx9h6#N`0t_KJnHZ!e|H;_qX%alRZS zmtY)*DCBfl)Jl*=CXkX4Z7{rFXqQN{%v|frJ)`FM=kp>7?C{+l7-tLh+$Nn5+ilYm z%mu6lE)?8uSL|v zTsO#)k3kN$^2F0u8uE?^1M7%}8I%Gwp$3|kC-m2f-vFbpaPImuGeeFH(LW%%9*{)) z^~#7E3xw-LqRR5MFT3>l)Z3GqX`>Kl!o7O?``em+rTdIRaWc{iFjaL+ULJLyz)9s8 z6(wq`0|NsqxB6@k!`g9*P|u@YQHd}QlH4Z{`^FyVdaaB^5uAzY!aKwYVkeMxRrMKD z=m@{>b!pmkZ6~tdrjEl6ByvZ_cmR{RAb!`~ z+=c=;ve538Vjpv%_#`p8wrJD>F67|BgKSQl={)N{J5Ju3Ap#GDvxP!!Iq!>Y0!zKS zgCSX=P72U*>A!YmMD@aW0N4cH!ESF;_fUAbFBtdNw|CTNhXSwGgHIyAhB%|e7lHH z4esk; zFuke=KEuK`3mKdH)piK}D55&di@paNku*nTc^mUY)q zFn0n}mPqnRNMLG1CAJb(Mwlp~QU+G_v$327cqu*nKLVM-#QQQ*-}y zfj!hZ!iWPwG$ZqD-y?JHni_?IBjadK(cKnCvuAU~`f8;5B<&%+1~Po2ng>vH2e+Nv z=TL3Ge&^XHQXbe2y;TToDnrmv5oO_VY?@v|HPxSjiR&Tbd_eC^g0d{sq$1E+xNz~J z#V!^|At=JM$w0#96FooDS682jmbm|KF1w<5<^yu+=BVl2PTjzm87biJ_ z0rD^}g&w_K@GSd+|91A}zNmarXgwg|oTyLZgb+n`$??x8^@va)G2*CL(b>k{LZL_V ze$Waq;~saRqfralPG00**1l%aP+O=%oVTVMRNeziLNaK~|KI2^{}Aj+)S|Omnw`R4 z-T#SjRaay^fLa{hHJY=S4>q@OYYDQ09!e2p`6$1UM_P!8RIZtbpY+|fW5>R43;U3( z35n=@Xn;9kLAzJ8hU*_^q2KBEw~zeSmr9p;!0Xp!UZ8mXVa}Bqo%D4~OiY@(C(BDp z6evyjP0n-E&7U6s_gBnhdiHnuaFIoCz-mEPSo{BEc4b$fnK#rn%hrny@dQYB`qk1_T~iBlcSvEDjisMO-NbVc?1EJIAj zmHG3wKeCghf2w=Ygx=`8GgYULZy84s`NtJME>Ky}`WSF~|IU>eJ;|J@N zE>oyP&LF0Xe3yw7UkB~On->~`d?JJ7JYHzXl0S8pJ9TpE3N?2bV;NT%rKFcKsc&#=z3O!G+}DFjvo%^C-%z2j|AM#S(_hAWgiNCY+tvO)DV6({x1OSxxPt4{$y8_jyni-) z@g3*UXYQ8O!)YK=kn3vRb2i?bBkI~%AK!Vo)uA(VcQoYH9b4zGOr2jyt&VCi6=<{IJNVhx$*iwv_98qd&2ZG()W;SI7e@);MiyK z((c?ZJx0-IgTM62cUQ!nKAB2iMq!k#7(RplAMWh9T!w6QaW+0G{nz}FrTS;qs=r#U z^(oR=HZMs&uE?UsH}qf7dgaE?s~RjsRRW=%;-o?-B?5a$v!#@nyKMJTfDq%m#LJ6bv9ZsaHY< zyG#7uA~ZS&C<_TUCeQ&lK)5-RyEvPBEH29s($CLlp<6(nMs9C0b7 zFOP~LtF4B7@;akjpx8(|74j>dxLZpCm_TksH(;5XF(=k+y@dh}${Q61kKpD_0~FsN z>m8sSh;JtQP#P`ObIpFp3omMD_-90w+)isZ2G2VG}P^P$Mq{I*}G2-)DzR}Kdj$54e##hh&`=#`! z*(XV~cz*^;4kC+lW28<;2CE@5lp^AFsO-VlDEIt3nd#01YT$>Cq`v}e6?qKQV>~>> zH}YSaxWpVQWbV{bccS)4;w!cOJKewkkFKS2{rU@CJCAie|0kVYQ&UsYPUVF%n*?bb zSfZVQc*23SBi4)D1~0Si`J^RvIVy`D>7Rd&CeBj(}qSK*@FpGBE=g6gAJDAD)D zeyM@fAge^S-BcZB$&o<>1)HYQ()fXuY1&+qt!UBZp zh;AZk)C3MkWL4cCbr2lIc=Q+Kg(!{_^<;ePAqcaPPlr@Bd=d(-QM z^i}IVckfX8vix=3yL|9Z;tko`L+oPDYGMtC&z>;bqgEiTlbKF2TX|pi1goiAbj^Kk z7NDhgjTo2^(}`sZOstnOpVrAXiIWsKUNQOYmO`hMGyk689e=g7zd7{cxyRMq@}>6B zd$^getET@Ph`J^qG3i{UNa0lH`LF5G{AXBr{CPvV;EQKn%Lbb6uJh2UQJHP9WRzu; zqO#z+xK4Su|38_wrvxvjp2c%DOt#kbZnn$f~HcmRxAX<0sBTY;Zo7&E3KjeYf{a zD>Hl}d>7>m9=I>RfxNHqdXIGD?o-^D&$*LRTmEn*``Av0E;}pr4dibab_Ksbc|t<2 zuWd=uVg{mGHJ(b=AkS`cqVmhz53`g^%=Ejl^ z#rwl6Gw34emX|wvPq5wQrBIAi(==HZlPj#3I*fU8BAzZ!yNH z^d+Z8ZmqZ`QlZQ=XvsugVZ;Ttb`D4Tcim-^hGM422 zy2SL~YhwvR=6qEp=V#NFS8dz&zgpG{Xy^}a>Rt*iAeNC$Tp-QXC3C)RPCfkt5+g=T zL_uby04*8JRq82fdHe9ygS3Bt=|c0}*Qlxzxiz#$h};@#*+dkwu<$%VKp{nm$Lz+z z$Evbtk`GZ8mlt?`S>B<<|hwm*IR-wh_C{Dx%+4VXHV)@zY!oMtvb{yda3YfO9@3C76; zF8~F6IuFX8v1EUsThO-ug`ak{)=s}Drp*IGUX-oTE&Ik<6y^3k$?Ku?1lOtM6@L1t zK;Q$%9;&4nMfg?mO;Z%30*FQ@N^BQk3OB@CJ;-?Whch{U$q@Y*QQuKPhA088L6p`Z z$R0?t7l{nao+Q)AGyC)|;U*w{MXq}JQ88r}$su&LAzG#e&ubn-25+KBvPW5D48235ClG6bR1*cGuBh-+n-w>%mwTTZFWuEyMG+KUvExS~ z6OohPpLo@cX3Rqr_%tSX4(wrO7OT9u?{O0%K{GQW<2Z`4$SqOM^*y8njVifTLE=?< zOr3&;dm{R>?fCSguw9EEJCjOlg2ZLgmu=UdG{c9pXmF6K@I3NJ5E<$K=X*znvZ4A` zi5k@p^0%G2HH~GtGW;AAOCtQF*tQ}71I39xKFXg;D7RU_y-i?#ED#+e?S_x(;FslWt^NPBwsV;bSQ$K-_SJd=C&w>_?7VI(qadnLKx4 z)GW(|DCvx67f2uE-e12QFP441@t24U6u;2gPQutQG9Qby-ykik+Bw)=9u!}|WOosv zoNA_3EEI<#(R%1Xll*>4IY=XVQck4;52dyYGuS?48XkKI1x zE4{_H8zE*D3jI7z3-<5A!~20=;v7vOJw$Xv``1jM*j9lmaU7vU{3J6oK`ug>?yyAw>Aez2F;u2V&=h=m7e}W2By>(It&upNchV); z*>IhUu11u)tjQ@=gGyBuM^RmkTThZ6$Nr=Gq}cmStxuh1a-&>jMoXhw_O4ZDO(WO1 zJeK^&?b?>p*Crjx{>&fsV-ef3CGdHPU+09q@M)c#o3kZ;COIWH%cf>U0CcPdE>oxl zJy8)`)Hf_E&(ifeVGb`TVfCVSq%0}<1g-Jf-s4hIQu_M(>qO(@;(Q$V8WlkJJZTs6w4wKrUH$RrNU>M>v7_s=S1JpXQ~LHdyx9a-I}?vpph9k4GgfXM> zJls`1;02@!S6f@V>ZYX@DnI9*(uZI_pZ9{?=`j+&6EqJQaMK6GL2bO_1Ut#Dt`<~j zJ(ZTJP_vuY^4!K0QT<}z_lv!p7!7reWH9v*D==W zJxi5~@O46`7ANAJIKDZpa)9k;NfjEqElg$%La+I=$c9fp+)Vo)vWy&xil}e>+cNQR zvCrChio>2GM)M92}HaMpp~E4i{ITwwlXpPtQ#Z)~_+C4Jy1p5`jdZ z{HC#Qh^$x#o5cODiBUhLdwVJ6!K9l_^y6El`A?CE3xFBRGJ)`rQ-*LctxvND5Z|Z?>@otU}z*KCwqkg9}=CCqD18M(D?HE z*uru;Kg75G84CAU_Ph91l$-^p`7u%keeYC*`u7q;EK{dX*cjSUy`_8{xgQIDt-{#4 zqL&9K8}t=x40{h6M9;3$i>TtP?|CKmarFburCq}>j&p$$gWcN|Ba zX!T$aYz^ARY!XD$d;&O}?>6-JSok<(y{@au$)0t%w{)9ScodV?EEQ6$MML>r_-ork zw~LbzIOzgCfLDLDh|%_Uo&M9MnpgLGj#IgqA39#FwID};@4`$n@8?i?u7|VMtG~J_ zta;pO>h%7lhbq8tx22nh=ghVYTOBAbEF`3^0l^KQw)M1gO4yxW+D}K}1At78AqTe$ zMd{ZtGFK@d6*7pey8GjK%ISb;Rbq8V&SRBUW5@^FiygW)~s&DL-^j;zn5`FnR~*q`_} zz6yL%agn^ADk&EXJu-Zk>57UE;QW+1wg_b#@i6)Z2A0V7^#PHzVpw^+NJ_FIb8+E< zm%?SuaBbJ{GyF^SLtGXE%m~U{7r%)}GDSPJhwTn~kJqA-0V^dr-DHzuVQQVOPCD~& zFN5D1ZB=He0@aq1RqpS{S7Qp3Nxh8Ij67-z#c>*Me3N7<(rYP+Hve4?jb=7I^g^QI zaT&XLpj6S6BVZ-);Q!w*|`s&zRwS?uXV_dj@Pnw7=P{Q*IaC) zMHcCws!tax3ZCwy04l$J{i+eCt>OcH>9c{5lq5}^&jJjF+q{B(mrl;f(EGJt9e}AJ zI>ux!xdG(viD)#5Fi{I3Vzp$@D#?#vCUoII5hh z6P?dxq9P*B6Ich~LuTVM^>uYG!K1p1OkO5^Xu;Wm4_#K7gO7I}U*q>JkaLLK^Cr`j zi(B{zuL3=C0x}!;&Ny?@Ni3eYvFXe=;qrw+-|J0Y-ybD zxKTN*SbuR_eMP&xo~8WKL0bOq&Ie^nY11ZI)}FsFb1tbFI%|LT3QkQ)>GM%Du0*%# zYQa#2__@?vx|utTYtB*3*L${89%y!xPH|nGl91H{QLljf?Mv2kT~kI{5TfzMESnD+ znPSX^7u}!p0``G~i}^$d#?W#g}0AtbA)C zBHn1mipn1~4xFd3<;2Maql#u?AbuX<;mgz<~iRC-dH_6qPmJox(K0r`%-~NO55;4D(}` ziqty|16YSrnv5wedkpM7e%Xpl#Uv07B$yMQnHMVJ#6F&jT+75(@HHYM!JCRwfZ9`- z!rp^%bCqC62))1n9vSpHllzan_q$5}F}i2=(C(V=pYA!5C)~6acLXa@R9-%OX`xDz z3(eQ2xX;0KR9h^Hy&fOx+sit|g<3m-Oz5oq<6#k#lB5?$(%i;|c7AqNZjXQcnAh&N*~Km4V=d1?4gSb&(8j z6KNa*4)gg{dIXURe1bzKUkD|&TpYTwY7MxWrN zfp>U=W~VN2LJwNv-oL*@W;r0@G?Zlc`>eF8{qLhoJNxAI@5d#Yciwyzg7Bd!F&_O> zeN1{(O}zeS-^G`~NP*db+W;isif_}jy{Lvo4u;LKWn-pHp5ED8~gL?(IK2X8H2d{)ocxDIO;4a5gku{x(O8J>pZ*ebn+;6OkZ z(J=eQW{XIX#_nlY*gjrf-c>=$%F2>nsAlW=v-JCi9V|{1d@;qBF4m!_R2}+vvQ$)* ziJY9eVW`+?`qdY*&n}++@A&u~;^QC~aty2G{N^+u7xn%N zT0mw2fEuz0Ao)M(N+hG9&BUmU<^i77^1*0g;xXRi$6a3lV-yv2NX3WEw>X7{+MWCG z&mRVM!_{0o#+N@`t#hs3_`>$d!hNJh{|4!2C9`8ABi9<(A$&UT@LX12{w;vDmuZE6 zJ3~{^HmUUScJ>#(zI67Q8XB=kN0A0&=%2uKjcOj0g^{H&YHa?*%aW<&v>2@JrJ!ev zka+XYhnw8*$eC2;URq9#kf6NPjZ8&_{;H9Qq)Tg6)f*f@va9|1O|M;^1>o*^)1VS2 zb)2`+(~CUD%`Ft8b?H*@r%$JG!1*rkm`?l$g&}k`65NGH3q6gN(-DVs^*HG^Lz3Cl zAsskzn``{j=#1mMD)r+#mrN=b|2Nl3UXSQqFy@4DGe5xb@&LloE|lOjw=GPjV#%!lrhC9xee#H4Y)f{9-%`vq^U6r0 znV6VUiYeVG7u?^o6)0RB&85e8~GUykGyNe`&$CzKilJ*>QSeS+lz+yVY&R&h$ zAA*%Pgw&2PwY1U7&nfxGE_H^%Ysc?WsQw-`nxXTK&59wM(l_p5gatzj z9Pa#~XLAx&_GIs05(}n1O9!mN+Hi2crix1BiG>jlpwRDQ)(h*O%#`7#)bg^=!M>d- z0SAm}J1B8kuuOy1s(F+R>oIcK8(d?&Zw`$N1yHbT`y#_hT0wY9abh7PD*1Q4eT?A=ze8Bjwdm_27)f3@PtwdVhLX#f9s=+7xs z4y{=>B)_Pss@g+n2Y}O)$kDRZrn>=mlGvPLJR}tx_BeX?zc!r#D~F6W6#*U@lHlsP zIuH-Ki(BU@26$OC8KqMGuCjM<-(zhCO5!-avq{HPMna;Jz}*4!(Gap4xMXA$3$#2x z&|pn*V!^bX8;pE|ed~MGw5#{fFF6l#p}+PAT6jqd_OfPLKM~zx7R?(+3G`0}CB!}v znt9;KIYh0$m!1q#TCw3UYtNcz&pH*e3OFPV9gg?7#Y4V*jZ02e6b+KmjOV+oq^cSP;~X6qCrwy9 zsW}s|w`lK*@?^}PF{SR~5fnMD2ofcdi4*3|T%>Z`1+{WdW% zaU5V~Rf0H0w)N`$-Lq|lBCJcX|9H;(%rqc=`vr%y#$3wDI0=LnWKT3EXg?So9koCu z`Mw?fcnHH%lGtN@b$<^j_aNcATRt{6MqU*e5#vSRE)KsxQvYOH?a))P6vR2V(pkAM z0W9O{>iT00_+EL8cG#ezsa8Cn^!D42+a9i83N9>Va`KBI$FyS-2&GDl>>r)Qs}nmd?v16_DX;tP4U<+ z#fZoB0xhZ}$UVi(^N)Wri+e+&(~~Ca?RHF|Xsk%^SqsC93o@X338Dcj<|Bh~!l+8F znL5&={lHy}Imnw~DJTeDTrY&tX28xEAP@-ToE`CSGAv-JC@`z0X`Ov#y7UTXKS@_$ zWJfJitTnQidv*uLS%#K6Ne=7-RP#yI(s1`)tnhFe1tgBU$g@Mrtg8L4$MOl9x*r%! z3ury?X5=wAbEdCh9U3;%?K#}K^=(NfzqbOP*U~OzUW12dZ@<<8&S7JQWi$bSN^+*3 zG07Tn1-=N@V9VU8@wFD_@eu|}iVl+ht{BGvRZGW+E^HsE*(BLb850~ZK^dZI<_I$! zot#|j&cFGu^xnbz5c6nCJE=zzFMg2d3ZMp5BgvH=`qo|W0%40sMMn<`Q_Wl&s+ipR zwf5W*N*l9h5&2?(zkd=qu>^8SSR@GfVbGio@bcE;%z4}|>C{{fo_nOO%50fJ8r4s5 zTWB$dp$|k1J>o9@M>Z!aa-PgGEhQGtn%k{h96znyuM4~x0sZnAHS$PQ0gNCR!-pLC z7VUje9}!kfV_%jay-3faOK&9_Z-p%qX$Vlw2EN7Vw6+$bU(zhTc%GWEw#lCKsr{B1 z^Fw5qQWKb5XSX3k&%Xr$w0L|T|9FN@P?5g>o~4*_Tq~ulz1180Ifs9P?Wx&E?aTF? z>~Q6+;|OP3;A_A{8Z1RT>SpX4Ot<{9;!!p-(bum+FhD5hfw$Wq*@BbHDnTUr2@;&Z z=(nK#spMOR6fdS4T zLTUNE(2dviqQ1WtC_qWGR+%()SF;B*@oyG0vouAkAKyAdEL+71?9^F%Q5)gcKgKZ5 z6E!#2Of~-V(yCdkkqtNdiZAGNlJ2o{>n1BY>*7S>Y;3YMCFEd}IJ3tIF7{r`6ly;A z?GBtdZ>;(zRF1enc<<#Rt>f_j_k|ZLJsh_UasE8=pjCgc=55v9B@;a^gDpvx zkI^*mqgnNG$~~U)z99JRV)4$RrGz9n^;z&+HN8hnY&P?O9Ft5DZcCQA;sp7=>s6}E z^x@g4Q&O=xYs4vCoq5`|g0`#`r7cOJ@oD0XTBo=UFPBv^$v4b;me=+%e%FhEE;{Lq2WtkUpMBEgpHt! zp?M0K~wWJ3jgmA;*ddg{@gd+(F<(rogt8XpwA-+;0I`AMI84s;jxPV3C4Tu z9d2v%5DO2W^lrWvKU+J4jYaUwR#^Bl>#=H>!fPe`gLCK8PkH>tA-(eRUpwJ=GDDL2pt*WZ}!r#BK zcl3Xe3ecW;E7x(?W9mh%nKZ-G{|dcH-{1`;{_of_j>2PM_n6-I_g7iEmb)W%Ngf2z zO~O8eM?xC^`|gBT5tMxt6}=JH6?YY)8KC2-SLAtJVkRjxs+fiPq40loLP)s(*vD2) zUVdxcq`PFrXG^>mihLsC-)(es^}T`s0AJG7)ckE@X_*Wtt~bDn7;xihDOC+XJkYh) z7H50dRql3--?;HsG~Q4SpzaNyA6lI6p<3>?{8vgbc!+xNU;$eX0(!H|mO~sI8xL}E za472N{JNStAHO)mMH$H1`_30bOFKF1$857lZL=ARG&q&RO@NEh;M}cyE}X5exagcB zUrJO|Wq1vh4Du&B(q{SZX2-5RqqGEpA@8Mz2FG=m)VzBpm@C<5mPr_Kgaz{2X~9M+ zNHm-4G^3}P*uKZnZdX_-79f|O6Alrb2vWtFsSa8QCFNtH8q(L-mxR%k5>=n>uNT|* zEVIL1$f8<`uNQ?j1*E_z)+r#~ZdBcJKJ)&*o{7y*DJO4+uD}w5jQ>}?hV%7z+UAF> z^^s$O=Qi&-AA?!<2rSICTw8~eWh%6$gD>x0Oay2XDCw~6MA9z=Ko zIh&w%Q#khN)2CM2HC{O4G@86lBcJDO@l|x;QFT=MR*SkXF>B6I>h6`L&or8UUMwFz zR2Q}QCs@ozLh&u!-M%)5LBZ^J+KL zpri_$h~+r8mrqgUps8SUX-%SBJI`LeMdwBi}b+G6SBhEz>p9 zQli)1n^kYwvyD>7IkP__ZLeXsl+fump)#l1thk17_Qeb-ihW&0p1JZq6=O=;MDwFv zfkV0aIR;wEb#uc*%f^@T!}GuMdls0P_RBoYy&y7oE<|%WnDx`rVvqlG>VZ-IO}759 z;frfT7aQ7I>vM`+*PNs5pUS)Vm#h(|T)9l#g@&*k zdD{YF_PTEHV$Gf%j!&BWDzk@SJGnD{=Cy)`zQ+zIP0j!2|C+WzqH6f~cHOc4H#WPk zBR_uyi=TZ>fAWfkQyoM-vfP@Ib2E;*Y$or``2*ftz4~aZy2EUbg@l}q#)r6uG!5b| z70BQIxCjTr%jp}E>2yT#vsboPDBrrGi6&*yC;#S6iwpG|G)Yp2r-pQn4_s0!>J{#FW+bR($BD0NBx0VBiAIX8* zmDK(%@ToXCdq;j?Rm#6+Xk<*eXnQlz?Uk4*D7$8KJM+Nu`@H@vc<&yx<}ACX_EQJ7 zjZV={$eLPGEp5`mJ*ULKkS)6YrkfWlZ?eB?%GH3r0JXhKuj|a}wSs&4^3C%NPAes{ zMh3@7vj@%p^NI30PF_de^#-wV#xRCXw08Kg&JVV5jn|}bEx+q=ku4TyntR62v`n@< z`q-U&TgP`@5t6yhWI@Oe`9E2Eh@9@05jyRq&^Az=xN$f#;Va$Jf&~wI27f-P6lC`_ z|2pdhFVRDBbw-QR2bM40#&YU`_7}Tr7Hza1-Y5ZcFguL3&2*QwFE5w_wyFCuPHvUx zUrbx=7e4H7t4EJ$^V2ot(ldA{hVD;08)nMCsd$C@=PPOash_H8=U#s3aZcQt`&D|? zv22N5e!I9N8W(KdE%e{_`0t;~nW@J(n*U6AY`NPjv@JJ7Jl;O*#PW&{le?4p#uU?E zzKhP@v86*uL9Ukm<3n=s;3+&;B`2h+%jkVsB%&e=>UejzJugpMUbOEXxN{wP*C9ct z$H$&`uu7N=MOK(-F7PaUw)6b?il4X19=7|LVD{(1+y2%$ndIbP;lxvvBT{<4V)J&L z&al4525||mi;@{L2~!H1P1}A(IQRVG<}l(oQ2(K6O3&bwz-RFl-(6Bxt;k7#d7gU3 zcuiLu)7x0PqKR{-($O#P&Dd|%bLUAo(%OX!vj)4opKS|AtJI8|Gx_hq>5G!EX2Eeu zHjc^(#3pY~4Y+-jolQFWmvn{7)(?B=cWx4~4||t7`J(>=*CC#;x>&!r)vGNs_f)2= zu=>F*DgJk4Ht*r$UB9&ql$4d(k=rUNDq`F)jf#$c7y_wylrYoPE%+ZM&DXjKcDL4B z1gEHYz1T+g)<3`Q$LM*7HMiu%`b_d@dOmO!@%(x??%3iy@uFJv$hnLh;oRuUzVkK# zG)8lOZ=M?xFF4z2Qh0FUc-1g%!WnUfBz8uDWTQvpHl}h6cP{Ii1nz$Fym3d_H$L;U z=EVKyE6%zw8EMXkUR=*?pILq5+4QC#q31qJB%CkP{T?`ZCn3IVFLIv59F^^s3qhle z`&xe2aKigrq)QbAU!K=rIP-3wt+)S=O3&%b8{MZ0w^j=o1wZfONZ=RRvH$B)&nEB5 z!^Lu9Ph`-(p2%*?VVYZ+%?(p!d*hGZlNg zwY7lr+ZXS6%2Fyv+IP3tG5PkSu9<2(k)0NDZ-0hfUO93m_p>O8wTnD{V##b3BpsFdrG4P?t1aU4Ph`SdCobgbDX+xC z7|BN;>Q_^!t$0&3vqSMX4|R}}Vtq_e%4k%|Xp!p<#oE+}hMSJDFR8y~@bZf)9lU#u z{BxJu#>GRv@%3zR)9V5|p4*pV+C_|ao+~ell<#{Y-~NwJj>lIXj@r~|h5DK(oX}d{ zKTi5hj&mLk9~1racg%OpU*k=&nxA^{!_?rLgun%+ZP8~;V!Ya~y`ZOi%NoCE{f6(H zMC)a?+bxexoyM}h9vt>4Pa2OE(xmxw*hW{_p`mqh+uDqc=PTY8*SVN!Oy>o(#1Xe- zn!uU9Q4Hnn!DF4rSbd8g2z~v>0%Pl@;gpUqX-(J6ljo^?(i>COb3Sk8Tr@CGXZ78nk#IGjs6Qje3?;1E zeCF~|=C?*|oS&W8g_N5NHVRwqyi?yYWZIM=QmOGdXUHW|(8e%hm#V$K;YegiLwJ}{ z&@O)|-jwN@?|obbhP<`k2Mm;Tq7rJmxz@*2Ut7)0m+WuU#Ax5THGxxXq?P{b*Nkm+ z{XTs8i$a|F942ZFAI*H*uQeOppst@viQP(P-ZQf>fv%pzRBC*5(i7SVFCM%W5@ zNQQMSR_r}E-7Gksde|?_DSi8G<77KyOZ(Qo4;H2QgHaK~GbXy{?PDA#S|kP+Xjaae z*wtCZ?7^k5yyp3ymT}nA`3x%+cO$e!)Ac<;ky1@wzQ^&P1nt$&BLo7-u-UCq6OfY%EPs2}hmLfg7H=gKJm^S=KtJ zRRtkc?UEVj8C^UhR`tG2uc6Yv-Jz^Ak+d!NS~A&3lp)P4#_F= zj1C#OlkX>R!@3B>Y_ywE$Jy3@d|SWOw+^dhk*USo*VIR;d{z}{tB+=AJa!XElQJv` zY1%cKpVoP7gU0U%e95^1xgl!ZjBjUlXq~^N(Nyy$BW&^p_158xeL;QxYDTjza{&*h z(!6^HU*=J7nGMYwJ9&+{TkH6Ny4c03aBGFSuL5i5=7p`hQfu^8aR8iJ*UcI0PdFVn z(4UBm9OPD=ZDZ-TN|d_Qq$FRvap8Ehkdub3iw>he%`fih(r3X<&Qu~9!kG?1agK&| zp`YpBuhg=j+B%y;)v7x?`q>l4s;YxFp@&hTmEP&ts~?pWDV_YhEK0m?v9<$&<(VxD z2o5f|KiR1HEC|ty#rCX~ zZOxy#joHLQ{TIbyQ2zyzUzonI~fKLQQZ=^hFL)`BgV8TZ{ulu6s6f?(J-9%2x`~?z>ws zIBLC@mbGc^!cbR8IW6aj=G&%0dU^<-4@K$^j5BVqpquTQlhClDTd#8Wif>L`Llx}CudO_GtJ%+lCZ#mZU#Hy}N}uZFu;}rj z(c78T6!8lAzM}8LIvkI z&=OI`{1!G(^aC<-4sGl{w!ta>L{Rm@Cc`WAaxMiYu_sTWMCL z;p!rCe}h<5ME}A%TN(4ax~E*+Be<9Hw-(1W)B8GPS850gNw)s2ac+KX&O7UM>=s{L zr!rc)Il7tNj;MTFaPT+M_|Z5?H?Z!e<-nnCZnhx36b|KEO?kc>VvT0URD5zGBI?`P zxQ34VJa=b^sySoQWS}=SQ{t|rK-K=f(J6iNz-qhq>qkPG6#0@EQlcI5?*05$#OVKz zN#;F4&Q4S18w+uPcM?w-xur!vpUMblrstaPD{}iQsnU9{)$^MXfI=4lkVlbsze50)D8)Z&q zUEi|2TN`##2kk!Q(tqD-tH^iWqz{8RxN7yVLxa&_I)({4INk(s>P{d14}J-rzWgBnwfeuWIG(hX|O8vWhBx@cW-Khgw;miz?Wx&w^f ze%mgFWgi{P{Kz!W+UK@bD(mjIw9K)T(Ck}H-aCS}8;bd>3R(>g#2Bqj*&iFyG%y>b zX#!SkoVxv|15%|-&2hWa_J3y-KK|4+YS5uNg$)hGl5MHSOF;|p-j3FbQ!*BX5aGTc%onu`AJPciB?+bO~nHQloq5TcD&Bbga; zfJ4;D`j6Oy&?M%vEtgx{g)B1jGL&u=UNF6JSCi|&@Nu8@H?JFf`n`?Kw&jB1?heND z8P6LXX&+l!H#sttz1X|jVDL$hMN`UQ5u5S)1$zVkg>MY=YaMj2X>FZuNr?%LQQ{ir z?6(S%{~-UuPdK7DM5s08`1huenEs!lFo4->=a?I-C#}vI_tgtBOLu>#ng@&mO{tc5 zXm3dQ9Xp%*r__W8_!Z82C+JU-DpMqo5FlOG0u~u!UM3JuyZ4+{!rqTPkFO@Vxea_w zxhYj5r&J`~a(owEDMM}lR%N6Ni!)QoE%~E`Y}|GtKCc%F?5drldnQIs|E$^m`N-U& z&6l)yrxJjetzNtyCp6Zh>Eu!7|MFeid`}w9^XVIt?W~TjR}To3oj7jMS#evU>Y#S@ z#JJsThRFWiPRF8U7F=7K!32eye&pY?p z@$P`ZG&p6GY8vy(?qYn<^IJ{ZWuCIgJeB-r>QEooEF3>o{-h?=NHV7BUDBEOsx8@! zJ`+hf!j=sZ*_tMHqp9|TPUDW5!`vL%KRt6pUTf+Km|koeY!l4+n$p|bR`Ju;yWt(* ziM9;&bH6WSTTF9(wzZXhJzWm=N=DbzQ^$?ZE2HLDkYLjtBX~wvY1?=&Db_=7o;VZSf&)pWFkiDsE9J5 z+lXp1SkxZVvI?n{LK%vRa_hpF+*-1UoX>aKv%9a|b6)4a^S7BX^E|)j_xzsU_x*W) zn=Ko2;=FEKo(XjeYmgo&&@q*n^M~5JEBZyUORLsM%+s~3uUSMAbivem@NlXm!PmnI;w_eXn5DRvpnoa-7w~2QR(E9rLG-A4fj|}((Th1 z@OL%p>Ic@CYESbiYfUg2Eq_YXkodSIuq!7?Pn(*>?k-&}*tuCSc^sbB4ERPCxhB0Z z+7xyk4*w!1bl0w2NXz>3iicO@sC<$H6qJ>im`LGCrk@c|HlRDC2^gOQ72;P_%YHes z`_-fB&c$|;Nrw7Mgq>RgtrcF1qn_FcHn<6H;Hgk0W__cUXW$NvQTA^I&Ej8Lr`Ve= zv_Gt?GNCgY2j^o^vZ_d&_f#63o3Dk?6spG9_nt6M%h>NEqK8kV#=ymsbKKUgtC0@M zGXQ2>|Kg60t#WH~%Mv=)s)^}{5Jxz4Y^5@nkCZh$f zO-tBqrSs%R$#(X%`fqV4r`R*~?b-fURd?*g#`(J&h>M8faZ!E&LGD&RKc#s_tJmWh z8<2H%(*s7QX+~q5eZorpTdda^LskR1k)F*v*A;#g68+oG|DSi^@7L5|A?>hfllr>N zsFqh&R@&#?sjM8@BqN#Lo;`a?+Q7wAual6F5RvPXi8sj#ePO|9$*jNcxyQy#=q`GN zuW0@|)Y3dX=>C$#-`Re!EeRXaCD^DbEZ(Zf~0{+9BBy@$&~Q?%wg3v%ht{` z_LC4q`;_}EmOxfWY-X$$m@HnmP`>9Jaj7&g&pBqyj@+pCQU9U75_)E7ThZsVL&O-X zV!J!Ad)q$yG5g+~%kw#%hGYqoW^<7|G}U6e*E&Q|Wn4&&&5j>|%Dy+0;c5$K(}sMl(aN z{b^4~AyRFL$Yuy_Muq&fyae9Nn=fw`&d`Wh&@Mj)CX5e}1ctKJ2b3~G;aF<$v%m9w z@yg@dYfz*$xVgwW95bt$^m2sCrJ})1I3&MEG^}Eg5_KKYuK28?1`HA=3!*C75b6$s z*49-ukfdD=xytQrRqpI|_DP<4%6*g_kf&r!1D_hJ@C$MI2xMg@t`JieLI!PVT$;2pl5~MOcR^YfgbIq_o7m zk@u{FV}qn8x5TXlYS!ZeDgZx1aV^FuCHCSag90~}svPDTKmj*e++U7Oq80MbSqs{^bp*t3Tq~k3H4+q*~5=Zfc!2QUm ziQM32!~Et*#g!9&)!gJiXJl$$d3m|AQ{5;LB0TGpCa^WN$IT(*og`N!Po89fU2;M=R5|G<7efLg z@Sa$Auj?(u+^AdZ&8L`nj9~!_BovLFACMKtRjnuNTn&P$@{|!fw1ta0?mB#(IEz7F zSBRMBFI;#uh-LK+3YmKIV9@3Fq1u{x@!}j|p(a!_P!!vRPAD@GH;x^8v5-euCUW`! z&rB~WvY6b`AB{cM<_(A?9;<32>AH3K!HF}|pw+RRJHNSyZnR`DhCN9f9TRgUqHv_MYV}1U_rg7kcYfQS7C7w4 zW7f2IL}0Qf4H98hn1a^wiwpa?%g~x8;C6vHFCh9}M>eE~2Q41})C2VsK*lIy!uoN# zWf`Kg(eC0m+RMic#^+X9Fn{GDZA8X1D3~rn@8^Q0iQO1og^2gr!#y!b0v$Se5RTY~ z2}A}Y))7{)+pj4zYHVyQEF>AD!-fqjLSqg|a~gMly2x}VkL@qT6@`1{Cb4fk?9w?i zjK1?lV3Jl)P(0~qs!Pyw^29cDxj_GpgrcF>t7rj(Xr9_;u7}V@Ardg2P0@=+q@rLW z&4E;jI3YU8;cqM4D1!karj|@unDY)+o}Acsr>NzJ|&ZG2<-VO7iyTVbv2&-SIQrDrA zT^;DRyU*C~Q7Tsdw$S3$lGhI?x}!+2GaC9Wjh6Q9o-LQde>u~%5gK$c7}lj*ZYjwy zQy;;FbYhypMt3YtoM#;^qa@oEv`CX-w8aE%v^s=;RLXzQ;se7&lDBgb`*aU+ znW-}~KQ{jpd7+P?e0{n+j)0rh2YQn+ZuSS%@A{|@Nd!g!;{(d~KVRbhJs1>_+lDnB zfbKc;_Cuaz%7f%=j-4oZ_h}g2PFpbFK+ly`vVUCgTpmtl-!rX|{Hj(IlZp==rtfmy z834$OUT#tE8yp7Xl?sgO*GnigF}SZZMr|9&Mq>B8rbi(%<>7HokB#REii9LWNeXsb ziXd|x5=oH_XLr}Lkw!j69yu^^Y><9@A+hX`$&+9dIE!KN&`w+}U-4n{6k<`8j=7zC zjr?uv$_%?OrILmK;tQr(b4oJ>beDLsHg#()jl#~M#5T7_(zI$AyXQzLbd!7(y5F^dKkkL148{S{S zfvYj|{CR5%S`)qs=7+dFZ5a$c?Z``4(+q>7XHPGH(5-Qm58z}Uc|2y=3W9SCZ9J2} zOFhG&7_4IuGEp5)#e)u)UlEFdBZhW=s+yb!YxzQ0U_0(HyOkdMyz6}hKz$H5iwuB% z8;XZJ{c--VTtl48IxaXJIYeT+Nk)fH(34`atX^EEaUr;Ik1ZzQYkEAL5#)mw%$o4f(nM0pYcFvZ{J|diMB0o68?^%2ohMQ|HkD;k%uWATWw@ ztho2NqFn{c^Hz- zV60#guPZr<5wR7MR{auH`vFLa-wvM;AgdjpW81kK;0G+I8`wd=yUJ{TMZe&ey3YZBk%Gl)Gv!1o3y4V#sXtg?suzio6%?{%8x7o zxau`OJ#~D@wKfG4IeDXL9;R3b1b*9~GyJg_xt{aM9?nTqrL*NL?CI_Kf2^$KYj#?y h|M$iB|LG+|ykXnBiYEiC7tmL+`_g_zj List[Tuple[int, ...]]: - """生成所有满足指数和 <= total_degree 的多元单项式指数。 - 对于 2 维 (W,L) 返回 (a,b) 列表。 - """ - exps: List[Tuple[int, ...]] = [] - def rec(prefix: Tuple[int, ...], remaining_dim: int, max_sum: int): - if remaining_dim == 1: - for i in range(max_sum + 1): - exps.append(prefix + (i,)) - return - for i in range(max_sum + 1): - rec(prefix + (i,), remaining_dim - 1, max_sum - i) - rec(tuple(), dim, total_degree) - return exps - - -class OrthonormalParamBasis: - """离散点集上 (W,L) 的多项式基 QR 正交化封装。""" - def __init__(self, W: np.ndarray, L: np.ndarray, total_degree: int): - assert W.shape == L.shape - self.W_mean = W.mean() - self.W_std = W.std() or 1.0 - self.L_mean = L.mean() - self.L_std = L.std() or 1.0 - self.exps = generate_monomial_exponents(2, total_degree) - # 构建单项式矩阵 M (Ns x M) - M = [] - Wn = (W - self.W_mean) / self.W_std - Ln = (L - self.L_mean) / self.L_std - for a, b in self.exps: - M.append((Wn ** a) * (Ln ** b)) - M = np.stack(M, axis=1) # (Ns, M) - # QR 分解获得离散正交基 Q, R 上三角 - Q, R = np.linalg.qr(M) - self.Q = Q # (Ns, Q) - self.R = R # (Q, Q) - self.R_inv = np.linalg.inv(R) - - @property - def n_basis(self) -> int: - return self.Q.shape[1] - - def phi_matrix(self) -> np.ndarray: - return self.Q # 对应样本点的基函数值 - - def eval(self, W: float, L: float) -> np.ndarray: - Wn = (W - self.W_mean) / self.W_std - Ln = (L - self.L_mean) / self.L_std - monomials = [] - for a, b in self.exps: - monomials.append((Wn ** a) * (Ln ** b)) - m = np.asarray(monomials) # (M,) - # M = Q R => 对新点 m = phi R => phi = m R^{-1} - phi = m @ self.R_inv - return phi # (Q,) - - -class ParametricVectorFitting: - """多参数 (W,L) 二端口/多端口 S 参数的简化 Parametric VF。 - - 模型形式 (单个端口对): - H(p; μ) ≈ Σ_{k=1..K} r_k(μ)/(p - a_k) + d(μ) + p e(μ) - 其中 μ = (W,L)。r_k, d, e 在参数域被正交基展开。 - """ - def __init__(self, config: ParametricVFConfig): - self.cfg = config - self.frequencies: np.ndarray | None = None # (Nf,) - self.poles: np.ndarray | None = None # (K,) - self.basis: OrthonormalParamBasis | None = None - # 残值系数: (nports, nports, K, Q) - self.residue_coeffs: np.ndarray | None = None - # d,e 系数: (nports, nports, Q) - self.d_coeffs: np.ndarray | None = None - self.e_coeffs: np.ndarray | None = None - self.meta: Dict[str, Any] = {} - - def _init_global_poles(self, freqs: np.ndarray) -> np.ndarray: - f_min, f_max = freqs.min(), freqs.max() - # 使用对数均匀分布 (排除 0) + 负实部保证稳定 - f_samples = np.logspace(np.log10(max(f_min, freqs[1] if len(freqs) > 1 else f_min*1.1)), np.log10(f_max), self.cfg.n_poles) - # 极点: -σ + j 2π f_k - sigma = self.cfg.damping * 2 * np.pi * f_max - poles = -sigma + 1j * 2 * np.pi * f_samples - return poles - - # ----------------- 内部工具函数 (前置以便类型检查) ----------------- - def _build_frequency_matrix(self, p_vec: np.ndarray, poles: np.ndarray) -> np.ndarray: - K = len(poles) - F_base = np.zeros((len(p_vec), K + 2), dtype=complex) - for k, a_k in enumerate(poles): - F_base[:, k] = 1.0 / (p_vec - a_k) - F_base[:, K] = 1.0 - F_base[:, K + 1] = p_vec - return F_base - - def _solve_residues_given_poles(self, aligned_S: List[np.ndarray], p_vec: np.ndarray): - assert self.poles is not None, "Poles not initialized" - nports = aligned_S[0].shape[1] - Ns = len(aligned_S) - K = len(self.poles) - residues_samples = np.zeros((nports, nports, K, Ns), dtype=complex) - d_samples = np.zeros((nports, nports, Ns), dtype=complex) - e_samples = np.zeros((nports, nports, Ns), dtype=complex) - F_base = self._build_frequency_matrix(p_vec, self.poles) - for s_idx, Sdata in enumerate(aligned_S): - for i in range(nports): - for j in range(nports): - y = Sdata[:, i, j] - x, *_ = np.linalg.lstsq(F_base, y, rcond=None) - residues_samples[i, j, :, s_idx] = x[:K] - d_samples[i, j, s_idx] = x[K] - e_samples[i, j, s_idx] = x[K+1] - return residues_samples, d_samples, e_samples - - def _select_characteristic_freqs(self, freqs: np.ndarray, S_list: List[np.ndarray], max_points: int, - method: str = 'curvature', threshold: Optional[float] = None) -> np.ndarray: - Nf = len(freqs) - if max_points >= Nf: - return np.arange(Nf) - if method == 'uniform': - idx = np.linspace(0, Nf-1, max_points, dtype=int) - return np.unique(idx) - mags = [] - for S in S_list: - mags.append(np.mean(np.abs(S), axis=(1,2))) # (Nf,) - agg = np.mean(np.stack(mags, axis=1), axis=1) # (Nf,) - curv = np.zeros_like(agg) - curv[1:-1] = np.abs(agg[2:] - 2*agg[1:-1] + agg[:-2]) - base = {0, Nf-1} - if threshold is not None: - cand = np.where(curv >= threshold)[0] - else: - order = np.argsort(-curv) - cand = order - selected = list(base) - for idx in cand: - if idx in base: - continue - selected.append(idx) - if len(selected) >= max_points: - break - return np.array(sorted(selected)) - - def _vf_iterate_poles(self, aligned_S: List[np.ndarray], p_vec: np.ndarray, - residues_samples: np.ndarray, d_samples: np.ndarray, e_samples: np.ndarray) -> np.ndarray: - assert self.poles is not None, "Poles not initialized" - nports = residues_samples.shape[0] - K = len(self.poles) - Ns = len(aligned_S) - Nf = len(p_vec) - n_pairs = nports * nports - n_r = n_pairs * K - n_c = K - n_d = n_pairs - n_e = n_pairs - total_unknowns = n_r + n_c + n_d + n_e - rows = Ns * Nf * n_pairs - A = np.zeros((rows, total_unknowns), dtype=complex) - b = np.zeros(rows, dtype=complex) - def pair_index(i,j): - return i*nports + j - row = 0 - for s_idx, Sdata in enumerate(aligned_S): - for f_idx, p in enumerate(p_vec): - denom = p - self.poles - basis_freq = 1.0/denom - for i in range(nports): - for j in range(nports): - H = Sdata[f_idx, i, j] - b[row] = H - pair = pair_index(i,j) - r_offset = pair*K - A[row, r_offset:r_offset+K] = basis_freq - A[row, n_r:n_r+K] = -H * basis_freq - A[row, n_r + n_c + pair] = 1.0 - A[row, n_r + n_c + n_d + pair] = p - row += 1 - x, *_ = np.linalg.lstsq(A, b, rcond=None) - c_k = x[n_r:n_r+n_c] - poly_total = np.poly(self.poles) - acc = poly_total.copy() - for k in range(K): - others = [self.poles[m] for m in range(K) if m != k] - term = c_k[k] * np.poly(others) - acc = np.polyadd(acc, term) - new_poles = np.roots(acc) - return new_poles - - def fit(self, dataset: List[ModelBasicDatasetUnit], network_loader=Capa.network): - # 载入所有 Network - networks: List[Network] = [] - W_list, L_list = [], [] - for unit in dataset: - net = network_loader(unit.result_dir, unit.id) - networks.append(net) - W_list.append(unit.parameters["W"]) - L_list.append(unit.parameters["L"]) - assert len(networks) > 0, "Empty dataset" - # 处理不同频率点数量: 取公共重叠频段并插值到统一网格 - freqs_list = [n.f for n in networks] - f_start = max(f[0] for f in freqs_list) - f_stop = min(f[-1] for f in freqs_list) - if f_stop <= f_start: - raise ValueError("No overlapping frequency range among networks") - # 选择基准网格: 使用第一个网络在重叠区间内的频点, 若其它网络缺该点则插值 - base_freqs = networks[0].f - mask = (base_freqs >= f_start) & (base_freqs <= f_stop) - freqs = base_freqs[mask] - # 如果其它网络点数差异较大, 可选: 统一稠密度 (这里只在差异>5%时使用最小长度网格) - for f in freqs_list[1:]: - # 判断差异比例 - if abs(len(f) - len(freqs)) / max(len(freqs),1) > 0.2: - # 使用所有网络在重叠区间内最少点数, 构造等间距网格 - n_common = min(len(ff[(ff>=f_start)&(ff<=f_stop)]) for ff in freqs_list) - freqs = np.linspace(f_start, f_stop, n_common) - warnings.warn("Frequency grids differ significantly; using uniform linspace for interpolation.") - break - # 对每个网络插值到 freqs (实部/虚部分开线性插值) - aligned_S = [] # List[(Nf, nports, nports)] - for net in networks: - f_src = net.f - # 获取在目标范围内原始数据 (若后面改用插值网格不同于原点, 仍需完整插值) - S_src = net.s # (Nf_src, nports, nports) - nports = net.nports - if not np.array_equal(f_src, freqs): - # 线性插值 (对每个端口对实部与虚部) - S_interp = np.zeros((len(freqs), nports, nports), dtype=complex) - for i in range(nports): - for j in range(nports): - y = S_src[:, i, j] - # numpy.interp 需要单调递增; 为类型检查器显式转为 np.ndarray - y_arr = np.asarray(y) - y_real = np.real(y_arr) - y_imag = np.imag(y_arr) - real_part = np.interp(freqs, f_src, y_real) - imag_part = np.interp(freqs, f_src, y_imag) - S_interp[:, i, j] = real_part + 1j * imag_part - aligned_S.append(S_interp) - else: - aligned_S.append(S_src[mask] if len(S_src)==len(base_freqs) and mask.sum()!=len(base_freqs) else S_src) - self.frequencies = freqs - nports = networks[0].nports - Ns = len(networks) - # 参数正交基 - self.basis = OrthonormalParamBasis(np.array(W_list), np.array(L_list), self.cfg.param_total_degree) - Phi = self.basis.phi_matrix() # (Ns, Q) - Qn = Phi.shape[1] - # (可选) 特征频率子采样 - if self.cfg.max_freq_points is not None and len(freqs) > self.cfg.max_freq_points: - sel_idx = self._select_characteristic_freqs(freqs, aligned_S, self.cfg.max_freq_points, - method=self.cfg.freq_selection_method, - threshold=self.cfg.curvature_threshold) - freqs = freqs[sel_idx] - aligned_S = [S[sel_idx, :, :] for S in aligned_S] - if self.cfg.verbose: - warnings.warn(f"Frequency reduced to {len(freqs)} points (method={self.cfg.freq_selection_method}).") - # 初始化极点 - self.poles = self._init_global_poles(freqs) - p_vec = 1j * 2 * np.pi * freqs - # 首次求解残值 - residues_samples, d_samples, e_samples = self._solve_residues_given_poles(aligned_S, p_vec) - # 极点迭代 (简化共享极点 VF) - if self.cfg.vf_iterations > 0: - for it in range(self.cfg.vf_iterations): - new_poles = self._vf_iterate_poles(aligned_S, p_vec, residues_samples, d_samples, e_samples) - # 反射不稳定极点 - new_poles = np.array([p if p.real < 0 else complex(-p.real, p.imag) for p in new_poles]) - # 去重: 近似重复 (距离 < 1e-3 * |p|) 过滤 - filtered = [] - for p in new_poles: - if all(abs(p - q) > 1e-3*max(1.0, abs(p)) for q in filtered): - filtered.append(p) - if self.cfg.verbose: - print(f"[VF-Iter {it+1}] poles -> {len(filtered)} kept") - self.poles = np.array(filtered, dtype=complex) - p_vec = 1j * 2 * np.pi * freqs - residues_samples, d_samples, e_samples = self._solve_residues_given_poles(aligned_S, p_vec) - # 在参数基上再拟合: 由于 Phi 正交, 系数 = Phi^H * values - PhiH = Phi.conj().T # (Q, Ns) - K = self.poles.shape[0] - self.residue_coeffs = np.zeros((nports, nports, K, Qn), dtype=complex) - self.d_coeffs = np.zeros((nports, nports, Qn), dtype=complex) - self.e_coeffs = np.zeros((nports, nports, Qn), dtype=complex) - for i in range(nports): - for j in range(nports): - for k in range(K): - vk = residues_samples[i, j, k, :] # (Ns,) - self.residue_coeffs[i, j, k, :] = PhiH @ vk - self.d_coeffs[i, j, :] = PhiH @ d_samples[i, j, :] - self.e_coeffs[i, j, :] = PhiH @ e_samples[i, j, :] - # 简单误差评估 (训练点重建 RMS) - F_base = self._build_frequency_matrix(p_vec, self.poles) - train_rms = self._compute_training_rms(residues_samples, d_samples, e_samples, F_base, aligned_S) - self.meta['train_rms'] = train_rms - self.meta['nports'] = nports - self.meta['Ns'] = Ns - self.meta['Q'] = Qn - return self - - def _compute_training_rms(self, residues_samples, d_samples, e_samples, F_base, original_S_list): - """计算每个样本的平均端口对 RMS 误差。""" - nports = residues_samples.shape[0] - K = residues_samples.shape[2] - Ns = residues_samples.shape[3] - Nf = F_base.shape[0] - rms = [] - for s in range(Ns): - err_sq_sum = 0.0 - count = 0 - S_ref = original_S_list[s] - for i in range(nports): - for j in range(nports): - x = np.concatenate([ - residues_samples[i, j, :, s], - [d_samples[i, j, s], e_samples[i, j, s]] - ]) # (K+2,) - y_hat = F_base @ x # (Nf,) - y_ref = S_ref[:, i, j] - err = y_ref - y_hat - err_sq_sum += np.mean(np.abs(err)**2) - count += 1 - rms.append(np.sqrt(err_sq_sum / max(count,1))) - return rms - - def evaluate(self, W: float, L: float, freqs: np.ndarray | None = None) -> np.ndarray: - """返回插值后的 S 参数 (Nf, nports, nports)。""" - assert self.frequencies is not None and self.poles is not None - assert self.basis is not None - assert self.residue_coeffs is not None and self.d_coeffs is not None and self.e_coeffs is not None, "Model not fitted" - if freqs is None: - freqs = self.frequencies - # 若请求频率网格不同, 需要重建频率矩阵 - p_vec = 1j * 2 * np.pi * freqs - K = len(self.poles) - nports = self.residue_coeffs.shape[0] - phi = self.basis.eval(W, L) # (Q,) - # 预分配 - S_eval = np.zeros((len(freqs), nports, nports), dtype=complex) - denom_cache = np.zeros((len(freqs), K), dtype=complex) - for k, a_k in enumerate(self.poles): - denom_cache[:, k] = 1.0 / (p_vec - a_k) - for i in range(nports): - for j in range(nports): - # 恢复参数依赖的残值/常数/斜率 - r_k = (self.residue_coeffs[i, j, :, :] @ phi) # (K,) - d = self.d_coeffs[i, j, :] @ phi - e = self.e_coeffs[i, j, :] @ phi - Hf = denom_cache @ r_k + d + p_vec * e # (Nf,) - S_eval[:, i, j] = Hf - return S_eval - - def export_s2p(self, W: float, L: float, filepath: str, freqs: np.ndarray | None = None, z0: float | complex = 50) -> str: - """在给定几何参数 (W,L) 下计算 S 并写入 Touchstone (.sNp) 文件。 - - 参数: - W, L: 几何参数 - filepath: 目标文件完整路径, 可包含或不包含 .s2p/.sNp 后缀 - freqs: 可选自定义频率数组 (需在训练频率范围内); 默认使用训练频率 - z0: 端口参考阻抗 (标量或与端口数相同长度可扩展) - - 返回: - 实际写入的文件路径 - """ - assert self.frequencies is not None, "Model not fitted" - S_eval = self.evaluate(W, L, freqs) - if freqs is None: - freqs = self.frequencies - nports = S_eval.shape[1] - freqs_arr = np.asarray(freqs) - # 构建 Network 对象 - ntw = rf.Network(f=freqs_arr, s=S_eval, z0=z0) - # 解析输出路径 - import os - base_dir = os.path.dirname(filepath) or '.' - os.makedirs(base_dir, exist_ok=True) - base_name = os.path.basename(filepath) - # 去掉可能已有的扩展名, 由 skrf 自动加 .sNp - if '.' in base_name: - base_name = '.'.join(base_name.split('.')[:-1]) or base_name - # 写文件 - ntw.write_touchstone(base_name, dir=base_dir) - # 构造最终文件名 (.sNp) - final_path = os.path.join(base_dir, f"{base_name}.s{nports}p") - return final_path - - def to_dict(self) -> Dict[str, Any]: - return { - 'config': self.cfg.__dict__, - 'frequencies': self.frequencies.tolist() if self.frequencies is not None else None, - 'poles': self.poles.tolist() if self.poles is not None else None, - 'basis': { - 'W_mean': getattr(self.basis, 'W_mean', None), - 'W_std': getattr(self.basis, 'W_std', None), - 'L_mean': getattr(self.basis, 'L_mean', None), - 'L_std': getattr(self.basis, 'L_std', None), - 'exps': getattr(self.basis, 'exps', None), - 'R': (lambda _R: _R.tolist() if _R is not None else None)(getattr(self.basis, 'R', None)), - }, - 'residue_coeffs': self.residue_coeffs.tolist() if self.residue_coeffs is not None else None, - 'd_coeffs': self.d_coeffs.tolist() if self.d_coeffs is not None else None, - 'e_coeffs': self.e_coeffs.tolist() if self.e_coeffs is not None else None, - 'meta': self.meta, - } - - @staticmethod - def from_dict(data: Dict[str, Any]) -> 'ParametricVectorFitting': - cfg = ParametricVFConfig(**data['config']) - obj = ParametricVectorFitting(cfg) - obj.frequencies = np.array(data['frequencies']) if data['frequencies'] is not None else None - obj.poles = np.array(data['poles']) if data['poles'] is not None else None - # 重建 basis (仅用于 eval); 重新构造 R_inv - if data.get('basis') and data['basis']['R'] is not None: - basis_stub = OrthonormalParamBasis(np.array([0.0, 1.0]), np.array([0.0, 1.0]), 1) # 占位 - basis_stub.W_mean = data['basis']['W_mean'] - basis_stub.W_std = data['basis']['W_std'] - basis_stub.L_mean = data['basis']['L_mean'] - basis_stub.L_std = data['basis']['L_std'] - basis_stub.exps = [tuple(e) for e in data['basis']['exps']] - R = np.array(data['basis']['R']) - basis_stub.R = R - basis_stub.R_inv = np.linalg.inv(R) - # Q 不需要 (仅用于训练点), 设为 None - basis_stub.Q = None # type: ignore - obj.basis = basis_stub - if data.get('residue_coeffs') is not None: - obj.residue_coeffs = np.array(data['residue_coeffs']) - if data.get('d_coeffs') is not None: - obj.d_coeffs = np.array(data['d_coeffs']) - if data.get('e_coeffs') is not None: - obj.e_coeffs = np.array(data['e_coeffs']) - obj.meta = data.get('meta', {}) - return obj - -# 划分 -def split_dataset(dataset: List[ModelBasicDatasetUnit], train_ratio: float = 0.7, shuffle: bool = True, seed: int = 0) -> Tuple[List[ModelBasicDatasetUnit], List[ModelBasicDatasetUnit]]: - """划分数据集 70% 训练 / 30% 测试 (默认可改)。""" - ds = list(dataset) - if shuffle: - rng = np.random.default_rng(seed) - rng.shuffle(ds) - n_train = int(len(ds) * train_ratio) - return ds[:n_train], ds[n_train:] - - -def _interp_s_to_freqs(S: np.ndarray, f_src: np.ndarray, f_tgt: np.ndarray) -> np.ndarray: - """将 S(f_src) 线性插值到 f_tgt (复数分离实虚部)。""" - if np.array_equal(f_src, f_tgt): - return S - nports = S.shape[1] - S_new = np.zeros((len(f_tgt), nports, nports), dtype=complex) - for i in range(nports): - for j in range(nports): - y = S[:, i, j] - real_part = np.interp(f_tgt, f_src, np.real(y)) - imag_part = np.interp(f_tgt, f_src, np.imag(y)) - S_new[:, i, j] = real_part + 1j * imag_part - return S_new - - -def evaluate_dataset(model: ParametricVectorFitting, - subset: List[ModelBasicDatasetUnit], - network_loader=Capa.network) -> Dict[str, Any]: - """对一个子集评估拟合质量，输出每样本与整体统计.""" - assert model.frequencies is not None - freqs = model.frequencies - results = [] - mag_err_db_all = [] - rms_all = [] - for unit in subset: - net = network_loader(unit.result_dir, unit.id) - S_ref = _interp_s_to_freqs(net.s, net.f, freqs) # (Nf, nports, nports) - S_pred = model.evaluate(unit.parameters["W"], unit.parameters["L"], freqs) - err = S_pred - S_ref - rms = float(np.sqrt(np.mean(np.abs(err)**2))) - # 幅度 dB 误差 - mag_ref_db = 20*np.log10(np.clip(np.abs(S_ref), 1e-15, None)) - mag_pred_db = 20*np.log10(np.clip(np.abs(S_pred), 1e-15, None)) - mag_err_db = float(np.mean(np.abs(mag_pred_db - mag_ref_db))) - rms_all.append(rms) - mag_err_db_all.append(mag_err_db) - results.append({ - 'id': unit.id, - 'W': unit.parameters['W'], - 'L': unit.parameters['L'], - 'rms': rms, - 'mag_err_db': mag_err_db - }) - summary = { - 'count': len(subset), - 'rms_mean': float(np.mean(rms_all)) if rms_all else None, - 'rms_max': float(np.max(rms_all)) if rms_all else None, - 'mag_err_db_mean': float(np.mean(mag_err_db_all)) if mag_err_db_all else None, - 'mag_err_db_max': float(np.max(mag_err_db_all)) if mag_err_db_all else None, - 'details': results - } - return summary - - -# 简要使用示例 (非执行): -if __name__ == "__main__": - - capa = Capa() - capa.sweep() - print("sweep finished") - capa.export('capa_results.json') - capa.clear() - capa.load('capa_results.json') - dataset = capa.results - - # 70% 训练 / 30% 测试 - train_set, test_set = split_dataset(dataset, train_ratio=0.7, shuffle=True, seed=42) - print(f"Train: {len(train_set)} Test: {len(test_set)}") - - cfg = ParametricVFConfig(n_poles=10, param_total_degree=3, vf_iterations=2, - max_freq_points=50) - pvf = ParametricVectorFitting(cfg).fit(train_set) - - train_metrics = evaluate_dataset(pvf, train_set) - test_metrics = evaluate_dataset(pvf, test_set) - - pvf.meta['train_eval'] = train_metrics - pvf.meta['test_eval'] = test_metrics - - def _print_metrics(title: str, m: Dict[str, Any]): - print(f"[{title}] samples={m['count']}" - f" RMS_mean={m['rms_mean']:.4e} RMS_max={m['rms_max']:.4e}" - f" |Δ|dB_mean={m['mag_err_db_mean']:.3f}dB |Δ|dB_max={m['mag_err_db_max']:.3f}dB") - - _print_metrics("TRAIN", train_metrics) - _print_metrics("TEST", test_metrics) - - # 示例：导出一个测试点的拟合结果为 s2p - if test_set: - sample = test_set[0] - out_path = pvf.export_s2p(sample.parameters["W"], sample.parameters["L"], filepath='outputs/pvf_test_sample') - print("Exported:", out_path) \ No newline at end of file