From de7d1fb473590e295ea6b759d811504bee7e1055 Mon Sep 17 00:00:00 2001 From: mayge Date: Sat, 27 Sep 2025 05:33:12 -0400 Subject: [PATCH] =?UTF-8?q?feat:=20=E6=89=93=E5=8C=85=E4=BB=A5=E5=8F=8A?= =?UTF-8?q?=E6=B7=BB=E5=8A=A0=E6=95=B0=E6=8D=AE=E9=9B=86?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .gitignore | 7 +- core/VFManager.py | 176 ------- core/basis/MultiplePortQR.py | 440 ------------------ core/basis/basis.py | 241 ---------- core/basis/basisQR.py | 259 ----------- core/basis/relaxed_basisQR.py | 361 -------------- env | 2 - main.py | 136 ------ ovf/__init__.py | 6 + ovf/core/VFManager.py | 317 +++++++++++++ {core => ovf/core}/__init__.py | 0 .../core}/basis/MultiPortOrthonormalBasis.py | 162 ++++--- {core => ovf/core}/basis/__init__.py | 0 {core => ovf/core}/sample.py | 0 {core => ovf/core}/utils.py | 5 +- pyproject.toml | 67 +++ test/test_numpy.py | 9 - test/test_orthonormal_basis.py | 102 ---- test/test_vector_fitting.py | 60 --- 19 files changed, 492 insertions(+), 1858 deletions(-) delete mode 100644 core/VFManager.py delete mode 100644 core/basis/MultiplePortQR.py delete mode 100644 core/basis/basis.py delete mode 100644 core/basis/basisQR.py delete mode 100644 core/basis/relaxed_basisQR.py delete mode 100644 env delete mode 100644 main.py create mode 100644 ovf/__init__.py create mode 100644 ovf/core/VFManager.py rename {core => ovf/core}/__init__.py (100%) rename {core => ovf/core}/basis/MultiPortOrthonormalBasis.py (75%) rename {core => ovf/core}/basis/__init__.py (100%) rename {core => ovf/core}/sample.py (100%) rename {core => ovf/core}/utils.py (88%) create mode 100644 pyproject.toml delete mode 100644 test/test_numpy.py delete mode 100644 test/test_orthonormal_basis.py delete mode 100644 test/test_vector_fitting.py diff --git a/.gitignore b/.gitignore index b25ab98..6c39187 100644 --- a/.gitignore +++ b/.gitignore @@ -6,4 +6,9 @@ __pypackages__/ env/ .venv/ .vscode/ -outputs/ \ No newline at end of file +.cache/ +test/ +outputs/ +ovf.egg-info/ +dist/ +build/ \ No newline at end of file diff --git a/core/VFManager.py b/core/VFManager.py deleted file mode 100644 index d45c22b..0000000 --- a/core/VFManager.py +++ /dev/null @@ -1,176 +0,0 @@ -import matplotlib.pyplot as plt -import numpy as np -from .basis.MultiPortOrthonormalBasis import MultiPortOrthonormalBasis -from .utils import generate_starting_poles - -class VFManager(): - def __init__( - self, - npoles_cplx, - freqs, - H, - model=MultiPortOrthonormalBasis, - iterations:int=5, - fit_constant:bool=True, - fit_proportional:bool=False, - dc_enforce:bool=False, - passivity_enforce:bool=True, - verbose:bool=True - ): - - - self.freqs=freqs - self.H=H - self.iterations=iterations - self.fit_constant=fit_constant - self.fit_proportional=fit_proportional - self.dc_enforce=dc_enforce - self.passivity_enforce=passivity_enforce - self.verbose=verbose - - self.nports = H.shape[1] - self.npoles_cplx = npoles_cplx - - self.least_squares_condition = [] - self.least_squares_row_condition = [] - self.least_squares_col_condition = [] - self.least_squares_rms_error = [] - self.eigenval_condition = [] - self.eigenval_row_condition = [] - self.eigenval_col_condition = [] - self.eigenval_rms_error = [] - - self.model_instance = None - self.model_responses_freqs = None - self.model_responses_H = None - - self.model=model - - def fit(self): - self.levi() - self.model_instance = self.sk_iteration() - return self.model - - def levi(self): - self.poles = generate_starting_poles(self.npoles_cplx,beta_min=1e4,beta_max=self.freqs[-1]*1.1) - self.model_instance=self.model( - H=self.H, - freqs=self.freqs, - poles=self.poles, - fit_constant=self.fit_constant, - fit_proportional=self.fit_proportional, - dc_enforce=self.dc_enforce, - passivity_enforce=self.passivity_enforce - ) - return self.model_instance - - def sk_iteration(self): - for i in range(self.iterations): - assert self.model_instance is not None ,"Please run levi() first." - self.poles = self.model_instance.next_poles - self.weights = self.model_instance.Dt - self.model_instance = self.model( - H=self.H, - freqs=self.freqs, - poles=self.poles, - weights=self.weights, - fit_constant=self.fit_constant, - fit_proportional=self.fit_proportional, - dc_enforce=self.dc_enforce, - passivity_enforce=self.passivity_enforce - ) - if self.verbose: - print(f"Iteration {i+1}/{self.iterations}") - print("A:",self.model_instance.A) - print("B:",self.model_instance.B) - print("C:",self.model_instance.C) - print("D:",self.model_instance.D) - print("next_pozles:",self.model_instance.next_poles) - print("Dt:",self.model_instance.Dt) - print("Dt/Dt_1:",np.linalg.norm(self.model_instance.Dt_Dt_1)) - self.least_squares_condition.append(self.model_instance.least_squares_condition) - self.least_squares_row_condition.append(self.model_instance.least_squares_row_condition) - self.least_squares_col_condition.append(self.model_instance.least_squares_col_condition) - self.least_squares_rms_error.append(self.model_instance.least_squares_rms_error) - self.eigenval_condition.append(self.model_instance.eigenval_condition) - self.eigenval_row_condition.append(self.model_instance.eigenval_row_condition) - self.eigenval_col_condition.append(self.model_instance.eigenval_col_condition) - self.eigenval_rms_error.append(self.model_instance.eigenval_rms_error) - return self.model_instance - - def plot_metrics(self,show:bool=True,save_path=None): - plt.figure(figsize=(16, 12)) - plt.subplot(4, 2, 1) - plt.plot(self.least_squares_condition, label='Least Squares Condition') - plt.legend() - plt.subplot(4, 2, 2) - plt.plot(self.least_squares_row_condition, label='Least Squares Row Condition') - plt.legend() - plt.subplot(4, 2, 3) - plt.plot(self.least_squares_col_condition, label='Least Squares Col Condition') - plt.legend() - plt.subplot(4, 2, 4) - plt.plot(self.least_squares_rms_error, label='Least Squares RMS Error') - plt.legend() - plt.subplot(4, 2, 5) - plt.plot(self.eigenval_condition, label='Eigenvalue Condition') - plt.legend() - plt.subplot(4, 2, 6) - plt.plot(self.eigenval_row_condition, label='Eigenvalue Row Condition') - plt.legend() - plt.subplot(4, 2, 7) - plt.plot(self.eigenval_col_condition, label='Eigenvalue Col Condition') - plt.legend() - plt.subplot(4, 2, 8) - plt.plot(self.eigenval_rms_error, label='Eigenvalue RMS Error') - plt.legend() - if show: - plt.show() - if save_path is not None: - if self.verbose: - print(f"Saving metrics plot to {save_path}/fitting_metrics.png") - plt.savefig(f"{save_path}/fitting_metrics.png") - - def plot_model_responses(self,show:bool=True,save_path=None): - assert self.model_responses_freqs is not None and self.model_responses_H is not None, "Please run get_model_responses() first." - for i in range(self.nports): - for j in range(self.nports): - plt.figure(figsize=(12, 6)) - plt.subplot(2, 2, 1) - plt.plot(self.freqs, np.abs(self.H[:,i,j]), 'o', ms=4, color='red', label='Input Samples') - plt.plot(self.model_responses_freqs, np.abs(self.model_responses_H[:,i,j]), '-', lw=2, color='k', label='Fit') - plt.title(f"Response i={i+1}, j={j+1}") - plt.ylabel("Magnitude") - plt.legend(loc="best") - plt.subplot(2, 2, 2) - plt.plot(self.freqs, np.angle(self.H[:,i,j],deg=True), 'o', ms=4, color='red', label='Input Samples') - plt.plot(self.model_responses_freqs, np.angle(self.model_responses_H[:,i,j],deg=True), '-', lw=2, color='k', label='Fit') - plt.title(f"Response i={i+1}, j={j+1}") - plt.ylabel("Phase (deg)") - plt.legend(loc="best") - plt.tight_layout() - plt.subplot(2, 2, 3) - plt.plot(self.freqs, np.real(self.H[:,i,j]), 'o', ms=4, color='red', label='Input Samples') - plt.plot(self.model_responses_freqs, np.real(self.model_responses_H[:,i,j]), '-', lw=2, color='k', label='Fit') - plt.title(f"Response i={i+1}, j={j+1}") - plt.ylabel("Real Part") - plt.legend(loc="best") - plt.subplot(2, 2, 4) - plt.plot(self.freqs, np.imag(self.H[:,i,j]), 'o', ms=4, color='red', label='Input Samples') - plt.plot(self.model_responses_freqs, np.imag(self.model_responses_H[:,i,j]), '-', lw=2, color='k', label='Fit') - plt.title(f"Response i={i+1}, j={j+1}") - plt.ylabel("Imag Part") - plt.legend(loc="best") - plt.tight_layout() - if show: - plt.show() - if save_path is not None: - if self.verbose: - print(f"Saving response plot for port {i+1},{j+1} to {save_path}/response_{i+1}_{j+1}.png") - plt.savefig(f"{save_path}/response_{i+1}_{j+1}.png") - - def get_model_responses(self,freqs): - assert self.model_instance is not None ,"Please run levi() and sk_iteration() first." - self.model_responses_freqs = freqs - self.model_responses_H = self.model_instance.get_model_responses(freqs) - return self.model_responses_H \ No newline at end of file diff --git a/core/basis/MultiplePortQR.py b/core/basis/MultiplePortQR.py deleted file mode 100644 index 2ea21c8..0000000 --- a/core/basis/MultiplePortQR.py +++ /dev/null @@ -1,440 +0,0 @@ -import numpy as np -from core.utils import generate_starting_poles -from scipy.linalg import block_diag -import skrf as rf -from skrf import VectorFitting -from core.sample import auto_select_multple_ports -import matplotlib.pyplot as plt -import random as rnd - -class MultiplePortQR: - def __init__(self,H,freqs,poles,weights=None,passivity=True,dc_enforce=False,fit_constant=True,fit_proportional=False): - self.least_squares_rms_error = None - self.least_squares_condition = None - self.eigenval_condition = None - self.eigenval_rms_error = None - - self.dc_tol = 1e-18 - - self.dc_enforce = dc_enforce - self.fit_constant = fit_constant - self.fit_proportional = fit_proportional - - # self.H = H - # self.freqs = freqs - self.freqs = freqs - self.H = H - self.ports = H.shape[1] - 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.Cw,self.w0,self.e = self.fit_denominator(self.H, weights=weights) - self.D = self.w0 - - self.Cr = None - - z = np.linalg.eigvals(self.A - self.B @ self.Cw) - p_next = z - - if passivity: - self.next_poles = self.passivity_enforce(p_next) - else: - 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 passivity_enforce(self,poles): - """enforce poles' real parts to be negative""" - enforced_poles = [] - for pole in poles: - if pole.real > 0: - pole = -np.conj(pole) - enforced_poles.append(pole) - return enforced_poles - - 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, weights=None, d0 = 1.0): - """ - 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|. - """ - K, N = self.Phi.shape - one = np.ones((K, 1), np.complex128) - Phi = self.Phi - dc_tol = 1e-18 - has_dc = self.dc_enforce and self.freqs[0] < dc_tol - keep = np.ones(K, dtype=bool) - - # SK weighting (applied only to the (73) rows we keep in LS) - - if has_dc: - # Enforce DC response exactly: - k0 = int(np.argmin(np.abs(self.freqs))) - keep[k0] = False - - if self.fit_constant: - Phi_w = np.hstack([one, Phi]) - index = 0 - M_kp = None - for i in range(self.ports): - for j in range(self.ports): - M0 = np.zeros((K,N*self.ports**2),dtype=complex) - M0[:,index*N:(index+1)*N] = Phi - M0 = np.hstack([M0, -(H[:,i,j].reshape(-1,1) * Phi_w)]).reshape((K, -1))[keep,:] # (K, 2N), complex - index+=1 - M_kp = M0 if M_kp is None else np.vstack([M_kp, M0]) - assert M_kp is not None - else: - index = 0 - M_kp = None - for i in range(self.ports): - for j in range(self.ports): - M0 = np.zeros((K,N*self.ports**2),dtype=complex) - M0[:,index*N:(index+1)*N] = Phi - M0 = np.hstack([M0, -(H[:,i,j].reshape(-1,1) * Phi)]).reshape((K, -1))[keep,:] # (K, 2N), complex - index+=1 - M_kp = M0 if M_kp is None else np.vstack([M_kp, M0]) - assert M_kp is not None - - if weights is None: - weights_kp = np.diag(np.ones(len(freqs[keep]) * self.ports**2, np.complex128)) - else: - weights_kp0 = weights[keep] - weights0 = [] - for i in range(self.ports **2 ): - for res in weights_kp0: - weights0.append(1/res) - weights_kp = np.diag(np.array(weights0)) - - - if has_dc: - M_w_kp = weights_kp @ M_kp - A_re = np.real(M_w_kp) - A_im = np.imag(M_w_kp) - mask = np.ones(K, dtype=bool); mask[k0] = False - # exact (unweighted) DC rows: - # A_dc_re = np.real(M_kp).reshape(1, -1) - # A_dc_im = np.imag(M_kp).reshape(1, -1) - else: - M_w_kp = weights_kp @ M_kp - A_re = np.real(M_w_kp) - A_im = np.imag(M_w_kp) - # A_dc_re = A_dc_im = None - - A_blocks = [A_re, A_im] - - if self.fit_constant: - Hk_sum = [] - for i in range(self.ports): - Hk_sum.append([]) - for j in range(self.ports): - Hk_kp0 = H[:,i,j][keep] - Hk_sum[i].append(np.sum(np.abs(Hk_kp0)**2)) - # Hk_kp = Hk_kp0 if Hk_kp is None else np.hstack([Hk_kp, Hk_kp0]) - K_keep = int(np.count_nonzero(keep)) - A_w0 = [] - b_w0 = [] - # Hk_sum = np.sum(np.abs(Hk_kp)**2) - for i in range(self.ports): - for j in range(self.ports): - beta_ij = float(np.sqrt(Hk_sum[i][j])) - mean_row = (beta_ij / K_keep) * np.sum(Phi_w[keep, :], axis=0) - A_w0.append(np.concatenate([np.zeros(N*self.ports**2, float), - np.real(mean_row).astype(float)] - ).reshape(1, -1)) - b_w0.append(np.array([beta_ij], float)) - b_w0 = np.asarray(b_w0).ravel() - - A_blocks += A_w0 - m = A_re.shape[0] + A_im.shape[0] - b = np.zeros(m, float) - b = np.concatenate([b, b_w0]) - else: - H_kp = None - for i in range(self.ports): - for j in range(self.ports): - H_0 = H[:,i,j][keep] - H_kp = H_0 if H_kp is None else np.hstack([H_kp, H_0]) - assert H_kp is not None - H_kp = weights_kp @ H_kp.reshape(-1,1) - - b_re = np.real(d0 * H_kp) - b_im = np.imag(d0 * H_kp) - b = np.concatenate([b_re.ravel(), b_im.ravel()]).astype(float) - - # ---- build final stacked-real system ---- - - - # if A_dc_re is not None: - # A_blocks += [A_dc_re, A_dc_im] - # b = np.concatenate([b, np.zeros(2, float)]) # DC rows → 0 - - # ---- QR solve for x = [c_H (N); c_w (N+1)] ---- - A = np.vstack(A_blocks).astype(float) - Q, R = np.linalg.qr(A, mode="reduced") - - if self.fit_constant: - Q2 = Q[:,Phi.shape[1] * self.ports**2:] - R22 = R[Phi.shape[1] * self.ports**2:,Phi.shape[1] * self.ports**2:] - else: - Q2 = Q[:,Phi.shape[1] * self.ports**2:] - R22 = R[Phi.shape[1] * self.ports**2:,Phi.shape[1] * self.ports**2:] - - x = np.linalg.solve(R22, Q2.T @ b) - - # diagnostics - resid = Q2 @ R22 @ x - b - self.least_squares_rms_error = float(np.sqrt(np.mean(resid**2))) - self.least_squares_condition = float(np.linalg.cond(R)) - - # split cw and return - # cw = x[N:] # last (N+1) entries = [w0, w_1..w_N] - # w0 = float(cw[0]) - # Cw = cw[1:].reshape(1, N) # row vector (1, N) - return self.extract_Cw_d_e(x,N,d0) - - def extract_Cw_d_e(self,C,N,d0=1.0): - if self.fit_proportional and self.fit_constant: - d = C[1] - e = C[0] - return C[2:].reshape(1, -1), d, e - elif self.fit_proportional and not self.fit_constant: - d = 0.0 - e = C[0] - return C[1:].reshape(1, -1), d, e - elif not self.fit_proportional and self.fit_constant: - d = C[0] - e = 0.0 - return C[1:].reshape(1, -1), d, e - else: - return C.reshape(1, -1), d0, 0.0 - - - def non_bias_Cr(self,w0): - A = np.asarray(self.Phi) - den = np.diag((w0 + self.Phi @ self.Cw.T).ravel()) - Cr = [] - for i in range(self.ports): - Cr.append([]) - for j in range(self.ports): - b = np.asarray(den) @ self.H[:,i,j].reshape(-1,1) - Cr_ij, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None) - Cr[i].append(Cr_ij) - return Cr - - def evaluate(self,freqs, w0): - H = np.zeros((len(freqs),self.ports,self.ports),dtype=complex) - s = 1j * 2*np.pi * np.asarray(freqs, float).ravel() - phi = self.generate_basis(s, self.poles) - den = w0 + phi @ self.Cw.T - if self.Cr is None: - self.Cr = self.non_bias_Cr(w0=w0) - for i in range(self.ports): - for j in range(self.ports): - num = phi @ self.Cr[i][j] - H[:,i,j] = (num / den).reshape(1,-1) - return H - -def noise(n:complex,coeff:float=0.05): - noise_r = rnd.gauss(-coeff * n.real, coeff * n.real) - noise_i = rnd.gauss(-coeff * n.imag, coeff * n.imag) - return complex(n.real + noise_r, n.imag + noise_i) - -if __name__ == "__main__": - start_point = 0 - id = 3000 - network = rf.Network(f"/tmp/paramer/simulation/{id}/{id}.s2p") - # network = rf.data.ring_slot - ports = network.nports - K = 10 - - full_freqences = network.f[start_point:] - noised_sampled_points = network.y[start_point:,:,:].reshape(-1,ports,ports) - sampled_points = network.y[start_point:,:,:].reshape(-1,ports,ports) - - # noised_sampled_points = network.y[start_point:,1,0].reshape(-1,1,1) - # sampled_points = network.y[start_point:,1,0].reshape(-1,1,1) - - H,freqs = auto_select_multple_ports(noised_sampled_points,full_freqences,max_points=20) - poles = generate_starting_poles(2,beta_min=1e4,beta_max=freqs[-1]*1.1) - - Dt_1 = np.ones((len(freqs),1),np.complex128) - # Levi step (no weighting): - basis = MultiplePortQR(H,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 = MultiplePortQR(H,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]) - H_evaluated = basis.evaluate(full_freqences, w0=basis.w0) - fitted_points = H_evaluated - sliced_freqences = freqs - - input_points = H - for i in range(ports): - for j in range(ports): - fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) - ax00 = axes[0][0] - ax00.plot(full_freqences, np.abs(sampled_points[:,i,j]), 'o', ms=4, color='red', label='Samples') - ax00.plot(full_freqences, np.abs(fitted_points[:,i,j]), '-', lw=2, color='k', label='Fit') - ax00.plot(sliced_freqences, np.abs(input_points[:,i,j]), 'x', ms=4, color='blue', label='Input Samples') - ax00.set_title(f"Response i={i+1}, j={j+1}") - ax00.set_ylabel("Magnitude") - ax00.legend(loc="best") - - ax01 = axes[0][1] - ax01.set_title(f"Response i={i+1}, j={j+1}") - ax01.set_ylabel("Phase (deg)") - ax01.plot(full_freqences, np.angle(sampled_points[:,i,j],deg=True), 'o', ms=4, color='red', label='Samples') - ax01.plot(full_freqences, np.angle(fitted_points[:,i,j],deg=True), '-', lw=2, color='k', label='Fit') - ax01.plot(sliced_freqences, np.angle(input_points[:,i,j],deg=True), 'x', ms=4, color='blue', label='Input Samples') - ax01.legend(loc="best") - - # ax00 = axes[0][0] - # ax00.plot(full_freqences, np.real(sampled_points[:,i,j]), 'o', ms=4, color='red', label='Samples') - # ax00.plot(full_freqences, np.real(fitted_points[:,i,j]), '-', lw=2, color='k', label='Fit') - # ax00.plot(sliced_freqences, np.real(input_points[:,i,j]), 'x', ms=4, color='blue', label='Input Samples') - # ax00.set_title(f"Response i={i+1}, j={j+1}") - # ax00.set_ylabel("Real Part") - # ax00.legend(loc="best") - - # ax01 = axes[0][1] - # ax01.set_title(f"Response i={i+1}, j={j+1}") - # ax01.set_ylabel("Imag Part") - # ax01.plot(full_freqences, np.imag(sampled_points[:,i,j]), 'o', ms=4, color='red', label='Samples') - # ax01.plot(full_freqences, np.imag(fitted_points[:,i,j]), '-', lw=2, color='k', label='Fit') - # ax01.plot(sliced_freqences, np.imag(input_points[:,i,j]), 'x', ms=4, color='blue', label='Input Samples') - # ax01.legend(loc="best") - - ax10 = axes[1][0] - ax10.plot(least_squares_condition, label='Least Squares Condition') - ax10.set_title("least_squares_condition") - ax10.set_ylabel("Magnitude") - ax10.legend(loc="best") - - ax11 = axes[1][1] - ax11.plot(least_squares_rms_error, label='Least Squares RMS Error') - ax11.set_title("least_squares_rms_error") - ax11.set_ylabel("Magnitude") - ax11.legend(loc="best") - - ax20 = axes[2][0] - ax20.plot(eigenval_condition, label='Eigenvalue Condition') - ax20.set_title("eigenval_condition") - ax20.set_ylabel("Magnitude") - ax20.legend(loc="best") - - ax21 = axes[2][1] - ax21.plot(eigenval_rms_error, label='Eigenvalue RMS Error') - ax21.set_title("eigenval_rms_error") - ax21.set_ylabel("Magnitude") - ax21.legend(loc="best") - fig.tight_layout() - plt.savefig(f"MultiplePortQR_port_{i+1}{j+1}.png") - print(f"Saved MultiplePortQR_port_{i+1}{j+1}.png") - - - - diff --git a/core/basis/basis.py b/core/basis/basis.py deleted file mode 100644 index 70bf970..0000000 --- a/core/basis/basis.py +++ /dev/null @@ -1,241 +0,0 @@ -import numpy as np -from core.utils import generate_starting_poles -from scipy.linalg import block_diag -import skrf as rf -from skrf import VectorFitting -from core.sample 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/basis/basisQR.py b/core/basis/basisQR.py deleted file mode 100644 index 6a00b1e..0000000 --- a/core/basis/basisQR.py +++ /dev/null @@ -1,259 +0,0 @@ -import numpy as np -from core.utils import generate_starting_poles -from scipy.linalg import block_diag -import skrf as rf -from skrf import VectorFitting -from core.sample import auto_select - -class BasicBasisQR: - 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.Cw = self.fit_denominator(self.H, d0=1.0, weights=weights) - self.Cr = None - - 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) - psi = weights @ self.Phi - - HPhi = H * psi - - A_re = np.hstack([np.real(psi), np.real(-HPhi)]) - A_im = np.hstack([np.imag(psi), np.imag(-HPhi)]) - - b_re = np.real(weights @ (d0 * H)) - b_im = np.imag(weights @ (d0 * H)) - - A = np.vstack([A_re, A_im]).astype(float) - - Q,R = np.linalg.qr(A) - R22 = R[A_re.shape[1]//2:, A_re.shape[1]//2:] - Q2 = Q[:, A_re.shape[1]//2:] - # 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(R22) @ (Q2.T @ b) - - self.least_squares_rms_error = np.sqrt(np.mean((Q2 @ R22 @ x - b)**2)) - self.least_squares_condition = np.linalg.cond(Q2 @ R22) - - Cw = self.vector_Cw(x) - return Cw - - def vector_Cw(self,x): - return x.T - - def non_bias_Cr(self,d0 = 1.0): - A = np.asarray(self.Phi) - den = np.diag((d0 + self.Phi @ self.Cw.T).ravel()) - b = np.asarray(den) @ self.H.reshape(-1,1) - Cr, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None) - return Cr - - def evaluate(self,freqs, d0=1.0): - s = 1j * 2*np.pi * np.asarray(freqs, float).ravel() - phi = self.generate_basis(s, self.poles) - den = d0 + phi @ self.Cw.T - if self.Cr is None: - self.Cr = self.non_bias_Cr(d0=d0) - num = phi @ self.Cr - 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 = BasicBasisQR(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 = BasicBasisQR(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:], d0=1.0) - import matplotlib.pyplot as plt - fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) - - ax00 = axes[0][0] - ax00.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') - ax00.plot(network.f[2:], np.abs(H11_evaluated), '-', lw=2, color='k', label='Fit') - ax00.plot(freqs, np.abs(H11), 'x', ms=4, color='blue', label='Input Samples') - ax00.set_title("Response i=0, j=0") - ax00.set_ylabel("Magnitude") - ax00.legend(loc="best") - - ax01 = axes[0][1] - ax01.set_title("Response i=0, j=0") - ax01.set_ylabel("Phase (deg)") - ax01.plot(network.f[2:], np.angle([network.y[i][0][0] for i in range(2,len(network.y))],deg=True), 'o', ms=4, color='red', label='Samples') - ax01.plot(network.f[2:], np.angle(H11_evaluated,deg=True), '-', lw=2, color='k', label='Fit') - ax01.plot(freqs, np.angle(H11,deg=True), 'x', ms=4, color='blue', label='Input Samples') - ax01.legend(loc="best") - - ax10 = axes[1][0] - ax10.plot(least_squares_condition, label='Least Squares Condition') - ax10.set_title("least_squares_condition") - ax10.set_ylabel("Magnitude") - ax10.legend(loc="best") - - ax11 = axes[1][1] - ax11.plot(least_squares_rms_error, label='Least Squares RMS Error') - ax11.set_title("least_squares_rms_error") - ax11.set_ylabel("Magnitude") - ax11.legend(loc="best") - - ax20 = axes[2][0] - ax20.plot(eigenval_condition, label='Eigenvalue Condition') - ax20.set_title("eigenval_condition") - ax20.set_ylabel("Magnitude") - ax20.legend(loc="best") - - ax21 = axes[2][1] - ax21.plot(eigenval_rms_error, label='Eigenvalue RMS Error') - ax21.set_title("eigenval_rms_error") - ax21.set_ylabel("Magnitude") - ax21.legend(loc="best") - fig.tight_layout() - plt.savefig(f"basic_basis_QR.png") - - - - diff --git a/core/basis/relaxed_basisQR.py b/core/basis/relaxed_basisQR.py deleted file mode 100644 index bf957ae..0000000 --- a/core/basis/relaxed_basisQR.py +++ /dev/null @@ -1,361 +0,0 @@ -import numpy as np -from core.utils import generate_starting_poles -from scipy.linalg import block_diag -import skrf as rf -from skrf import VectorFitting -from core.sample import auto_select -import random as rnd - -class RelaxedBasicBasisQR: - def __init__(self,H,freqs,poles,weights=None,passivity=True,dc_enforce=True,fit_constant=True,fit_proportional=False): - self.least_squares_rms_error = None - self.least_squares_condition = None - self.eigenval_condition = None - self.eigenval_rms_error = None - - self.dc_tol = 1e-18 - - self.dc_enforce = dc_enforce - self.fit_constant = fit_constant - self.fit_proportional = fit_proportional - - # self.H = H - # self.freqs = freqs - self.freqs = freqs - self.H = H - 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.Cw,self.w0,self.e = self.fit_denominator(self.H, weights=weights) - self.D = self.w0 - - self.Cr = None - - z = np.linalg.eigvals(self.A - self.B @ self.Cw) - p_next = -z - - if passivity: - self.next_poles = self.passivity_enforce(p_next) - else: - 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 passivity_enforce(self,poles): - """enforce poles' real parts to be negative""" - enforced_poles = [] - for pole in poles: - if pole.real > 0: - pole = -np.conj(pole) - enforced_poles.append(pole) - return enforced_poles - - 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, weights=None, d0 = 1.0): - """ - 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|. - """ - H = np.asarray(H, np.complex128).reshape(-1,1) - K, N = self.Phi.shape - one = np.ones((K, 1), np.complex128) - Phi = self.Phi - dc_tol = 1e-18 - has_dc = self.dc_enforce and self.freqs[0] < dc_tol - keep = np.ones(K, dtype=bool) - - # SK weighting (applied only to the (73) rows we keep in LS) - if weights is None: - weights = np.diag(np.ones(len(H), np.complex128)) - else: - weights = np.diag([1/res for res in weights]) - - if self.fit_constant: - Phi_w = np.hstack([one, Phi]) - M = np.hstack([Phi, -(H * Phi_w)]) # (K, 2N+1), complex - else: - M = np.hstack([Phi, -(H * Phi)]) # (K, 2N), complex - - if has_dc: - # Enforce DC response exactly: - k0 = int(np.argmin(np.abs(self.freqs))) - keep[k0] = False - M_w = weights @ M - A_re = np.real(M_w[keep, :]) - A_im = np.imag(M_w[keep, :]) - mask = np.ones(K, dtype=bool); mask[k0] = False - # exact (unweighted) DC rows: - A_dc_re = np.real(M[k0, :]).reshape(1, -1) - A_dc_im = np.imag(M[k0, :]).reshape(1, -1) - else: - M_w = weights @ M - A_re = np.real(M_w) - A_im = np.imag(M_w) - A_dc_re = A_dc_im = None - - A_blocks = [A_re, A_im] - - if self.fit_constant: - beta = float(np.sqrt(np.sum(np.abs(H)**2))) - mean_row = (beta / K) * np.sum(Phi_w, axis=0) - A_w0 = np.concatenate([np.zeros(N, float), - np.real(mean_row).astype(float)] - ).reshape(1, -1) - b_w0 = np.array([beta], float) - - A_blocks += [A_w0] - m = A_re.shape[0] + A_im.shape[0] - b = np.zeros(m, float) - b = np.concatenate([b, b_w0]) - else: - H_kp = (weights @ H)[keep,:] - b_re = np.real(d0 * H_kp) - b_im = np.imag(d0 * H_kp) - b = np.concatenate([b_re.ravel(), b_im.ravel()]).astype(float) - - # ---- build final stacked-real system ---- - - - # if A_dc_re is not None: - # A_blocks += [A_dc_re, A_dc_im] - # b = np.concatenate([b, np.zeros(2, float)]) # DC rows → 0 - - # ---- QR solve for x = [c_H (N); c_w (N+1)] ---- - A = np.vstack(A_blocks).astype(float) - Q, R = np.linalg.qr(A, mode="reduced") - - if self.fit_constant: - Q2 = Q[:,A.shape[1]//2:] - R22 = R[A.shape[1]//2:,A.shape[1]//2:] - else: - Q2 = Q[:,A.shape[1]//2:] - R22 = R[A.shape[1]//2:,A.shape[1]//2:] - - x = np.linalg.solve(R22, Q2.T @ b) - - # diagnostics - resid = Q2 @ R22 @ x - b - self.least_squares_rms_error = float(np.sqrt(np.mean(resid**2))) - self.least_squares_condition = float(np.linalg.cond(R)) - - # split cw and return - # cw = x[N:] # last (N+1) entries = [w0, w_1..w_N] - # w0 = float(cw[0]) - # Cw = cw[1:].reshape(1, N) # row vector (1, N) - return self.extract_Cw_d_e(x,N,d0) - - def extract_Cw_d_e(self,C,N,d0=1.0): - if self.fit_proportional and self.fit_constant: - d = C[1] - e = C[0] - return C[2:].reshape(1, -1), d, e - elif self.fit_proportional and not self.fit_constant: - d = 0.0 - e = C[0] - return C[1:].reshape(1, -1), d, e - elif not self.fit_proportional and self.fit_constant: - d = C[0] - e = 0.0 - return C[1:].reshape(1, -1), d, e - else: - return C.reshape(1, -1), d0, 0.0 - - - def non_bias_Cr(self,w0): - A = np.asarray(self.Phi) - den = np.diag((w0 + self.Phi @ self.Cw.T).ravel()) - b = np.asarray(den) @ self.H.reshape(-1,1) - Cr, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None) - return Cr - - def evaluate(self,freqs, w0): - s = 1j * 2*np.pi * np.asarray(freqs, float).ravel() - phi = self.generate_basis(s, self.poles) - den = w0 + phi @ self.Cw.T - if self.Cr is None: - self.Cr = self.non_bias_Cr(w0=w0) - num = phi @ self.Cr - H = num / den - return H.ravel() - -def noise(n:complex,coeff:float=0.05): - noise_r = rnd.gauss(-coeff * n.real, coeff * n.real) - noise_i = rnd.gauss(-coeff * n.imag, coeff * n.imag) - return complex(n.real + noise_r, n.imag + noise_i) - -if __name__ == "__main__": - start_point = 0 - network = rf.Network("/tmp/paramer/simulation/3000/3000.s2p") - K = 10 - - full_freqences = network.f[start_point:] - noised_sampled_points = [(network.y[i][1][1]) for i in range(start_point,len(network.y))] - sampled_points = [network.y[i][1][1] for i in range(start_point,len(network.y))] - - H11,freqs = auto_select(noised_sampled_points,full_freqences,max_points=20) - poles = generate_starting_poles(2,beta_min=1e4,beta_max=freqs[-1]*1.1) - - Dt_1 = np.ones((len(freqs),1),np.complex128) - # Levi step (no weighting): - basis = RelaxedBasicBasisQR(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 = RelaxedBasicBasisQR(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[start_point:], w0=basis.w0) - import matplotlib.pyplot as plt - fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) - - ax00 = axes[0][0] - fitted_points = H11_evaluated - sliced_freqences = freqs - - input_points = H11 - ax00.plot(full_freqences, np.abs(sampled_points), 'o', ms=4, color='red', label='Samples') - ax00.plot(full_freqences, np.abs(fitted_points), '-', lw=2, color='k', label='Fit') - ax00.plot(sliced_freqences, np.abs(input_points), 'x', ms=4, color='blue', label='Input Samples') - ax00.set_title("Response i=0, j=0") - ax00.set_ylabel("Magnitude") - ax00.legend(loc="best") - - ax01 = axes[0][1] - ax01.set_title("Response i=0, j=0") - ax01.set_ylabel("Phase (deg)") - ax01.plot(full_freqences, np.angle(sampled_points,deg=True), 'o', ms=4, color='red', label='Samples') - ax01.plot(full_freqences, np.angle(fitted_points,deg=True), '-', lw=2, color='k', label='Fit') - ax01.plot(sliced_freqences, np.angle(input_points,deg=True), 'x', ms=4, color='blue', label='Input Samples') - ax01.legend(loc="best") - - ax10 = axes[1][0] - ax10.plot(least_squares_condition, label='Least Squares Condition') - ax10.set_title("least_squares_condition") - ax10.set_ylabel("Magnitude") - ax10.legend(loc="best") - - ax11 = axes[1][1] - ax11.plot(least_squares_rms_error, label='Least Squares RMS Error') - ax11.set_title("least_squares_rms_error") - ax11.set_ylabel("Magnitude") - ax11.legend(loc="best") - - ax20 = axes[2][0] - ax20.plot(eigenval_condition, label='Eigenvalue Condition') - ax20.set_title("eigenval_condition") - ax20.set_ylabel("Magnitude") - ax20.legend(loc="best") - - ax21 = axes[2][1] - ax21.plot(eigenval_rms_error, label='Eigenvalue RMS Error') - ax21.set_title("eigenval_rms_error") - ax21.set_ylabel("Magnitude") - ax21.legend(loc="best") - fig.tight_layout() - plt.savefig(f"relaxed_basic_basis_QR.png") - print("Saved relaxed_basic_basis_QR.png") - - - - diff --git a/env b/env deleted file mode 100644 index 27df2ed..0000000 --- a/env +++ /dev/null @@ -1,2 +0,0 @@ -export PYTHONPATH=$(pwd)/core:$PYTHONPATH -export PYTHONPATH=$(pwd)/models:$PYTHONPATH \ No newline at end of file diff --git a/main.py b/main.py deleted file mode 100644 index b8e565e..0000000 --- a/main.py +++ /dev/null @@ -1,136 +0,0 @@ -import skrf as rf -from core.VFManager import VFManager -from core.utils import generate_starting_poles -from core.sample import auto_select_multple_ports -from core.basis.MultiPortOrthonormalBasis import MultiPortOrthonormalBasis - - -if __name__ == "__main__": - start_point = 0 - id = 3000 - network = rf.Network(f"/tmp/paramer/simulation/{id}/{id}.s2p") - # network = rf.data.ring_slot - ports = network.nports - K = 100 - - full_freqences = network.f[start_point:] - noised_sampled_points = network.y[start_point:,:,:].reshape(-1,ports,ports) - sampled_points = network.y[start_point:,:,:].reshape(-1,ports,ports) - - # noised_sampled_points = network.y[start_point:,0,0].reshape(-1,1,1) - # sampled_points = network.y[start_point:,0,0].reshape(-1,1,1) - - H,freqs = auto_select_multple_ports(noised_sampled_points,full_freqences,max_points=20) - - vf = VFManager(npoles_cplx=2,freqs=freqs,H=H,model=MultiPortOrthonormalBasis,iterations=K,verbose=False) - model = vf.fit() - vf.plot_metrics(show=False,save_path="outputs") - model_responses = vf.get_model_responses(full_freqences) - vf.plot_model_responses(show=False,save_path="outputs") - - # # Original plot functions - - # Dt_1 = np.ones((len(freqs),1),np.complex128) - # # Levi step (no weighting): - # basis = MultiPortOrthonormalBasis(H,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.C) - # 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 = MultiPortOrthonormalBasis(H,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.C) - # 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]) - # H_evaluated = basis.get_model_responses(full_freqences) - # fitted_points = H_evaluated - # sliced_freqences = freqs - - # input_points = H - # for i in range(ports): - # for j in range(ports): - # fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) - # ax00 = axes[0][0] - # ax00.plot(full_freqences, np.abs(sampled_points[:,i,j]), 'o', ms=4, color='red', label='Samples') - # ax00.plot(full_freqences, np.abs(fitted_points[:,i,j]), '-', lw=2, color='k', label='Fit') - # ax00.plot(sliced_freqences, np.abs(input_points[:,i,j]), 'x', ms=4, color='blue', label='Input Samples') - # ax00.set_title(f"Response i={i+1}, j={j+1}") - # ax00.set_ylabel("Magnitude") - # ax00.legend(loc="best") - - # ax01 = axes[0][1] - # ax01.set_title(f"Response i={i+1}, j={j+1}") - # ax01.set_ylabel("Phase (deg)") - # ax01.plot(full_freqences, np.angle(sampled_points[:,i,j],deg=True), 'o', ms=4, color='red', label='Samples') - # ax01.plot(full_freqences, np.angle(fitted_points[:,i,j],deg=True), '-', lw=2, color='k', label='Fit') - # ax01.plot(sliced_freqences, np.angle(input_points[:,i,j],deg=True), 'x', ms=4, color='blue', label='Input Samples') - # ax01.legend(loc="best") - - # # ax00 = axes[0][0] - # # ax00.plot(full_freqences, np.real(sampled_points[:,i,j]), 'o', ms=4, color='red', label='Samples') - # # ax00.plot(full_freqences, np.real(fitted_points[:,i,j]), '-', lw=2, color='k', label='Fit') - # # ax00.plot(sliced_freqences, np.real(input_points[:,i,j]), 'x', ms=4, color='blue', label='Input Samples') - # # ax00.set_title(f"Response i={i+1}, j={j+1}") - # # ax00.set_ylabel("Real Part") - # # ax00.legend(loc="best") - - # # ax01 = axes[0][1] - # # ax01.set_title(f"Response i={i+1}, j={j+1}") - # # ax01.set_ylabel("Imag Part") - # # ax01.plot(full_freqences, np.imag(sampled_points[:,i,j]), 'o', ms=4, color='red', label='Samples') - # # ax01.plot(full_freqences, np.imag(fitted_points[:,i,j]), '-', lw=2, color='k', label='Fit') - # # ax01.plot(sliced_freqences, np.imag(input_points[:,i,j]), 'x', ms=4, color='blue', label='Input Samples') - # # ax01.legend(loc="best") - - # ax10 = axes[1][0] - # ax10.plot(least_squares_condition, label='Least Squares Condition') - # ax10.set_title("least_squares_condition") - # ax10.set_ylabel("Magnitude") - # ax10.legend(loc="best") - - # ax11 = axes[1][1] - # ax11.plot(least_squares_rms_error, label='Least Squares RMS Error') - # ax11.set_title("least_squares_rms_error") - # ax11.set_ylabel("Magnitude") - # ax11.legend(loc="best") - - # ax20 = axes[2][0] - # ax20.plot(eigenval_condition, label='Eigenvalue Condition') - # ax20.set_title("eigenval_condition") - # ax20.set_ylabel("Magnitude") - # ax20.legend(loc="best") - - # ax21 = axes[2][1] - # ax21.plot(eigenval_rms_error, label='Eigenvalue RMS Error') - # ax21.set_title("eigenval_rms_error") - # ax21.set_ylabel("Magnitude") - # ax21.legend(loc="best") - # fig.tight_layout() - # plt.savefig(f"MultiplePortQR_port_{i+1}{j+1}.png") - # print(f"Saved MultiplePortQR_port_{i+1}{j+1}.png") \ No newline at end of file diff --git a/ovf/__init__.py b/ovf/__init__.py new file mode 100644 index 0000000..ec39e38 --- /dev/null +++ b/ovf/__init__.py @@ -0,0 +1,6 @@ +from ovf.core.VFManager import VFManager +from ovf.core.sample import auto_select_multple_ports +from ovf.core.basis.MultiPortOrthonormalBasis import MultiPortOrthonormalBasis + + +__all__ = ["VFManager","auto_select_multple_ports","MultiPortOrthonormalBasis"] \ No newline at end of file diff --git a/ovf/core/VFManager.py b/ovf/core/VFManager.py new file mode 100644 index 0000000..692cf11 --- /dev/null +++ b/ovf/core/VFManager.py @@ -0,0 +1,317 @@ +import matplotlib.pyplot as plt +import numpy as np +import os +from ovf.core.basis.MultiPortOrthonormalBasis import MultiPortOrthonormalBasis +from ovf.core.utils import generate_starting_poles +import json +import pickle + +class VFManager(): + def __init__( + self, + npoles_cplx, + freqs, + H, + model=MultiPortOrthonormalBasis, + iterations:int=5, + fit_constant:bool=True, + fit_proportional:bool=False, + dc_enforce:bool=False, + passivity_enforce:bool=True, + verbose:bool=True + ): + + + self.freqs=freqs + self.H=H + self.iterations=iterations + self.fit_constant=fit_constant + self.fit_proportional=fit_proportional + self.dc_enforce=dc_enforce + self.passivity_enforce=passivity_enforce + self.verbose=verbose + + self.nports = H.shape[1] + self.npoles_cplx = npoles_cplx + + self.least_squares_condition = [] + self.least_squares_row_condition = [] + self.least_squares_col_condition = [] + self.least_squares_rms_error = [] + self.eigenval_condition = [] + self.eigenval_row_condition = [] + self.eigenval_col_condition = [] + self.eigenval_rms_error = [] + + self.model_instance:MultiPortOrthonormalBasis|None = None + self.model_responses_freqs = None + self.model_responses_H = None + + self.model=model + + @classmethod + def load(cls,dirname): + instance = cls._load_npz(f"{dirname}/model.npz") + instance._load_model_instance(f"{dirname}/model_instance.pkl") + return instance + + def write(self,dirname): + os.makedirs(dirname, exist_ok=True) + self._write_npz(f"{dirname}/model.npz") + self._write_model_instance(f"{dirname}/model_instance.pkl") + + def _write_npz(self,filename): + assert self.model_instance is not None ,"Please run levi() and sk_iteration() first." + np.savez(filename, + model_name=self.model.__name__, + poles=self.model_instance.poles, + freqs=self.freqs, + H=self.H, + iterations=self.iterations, + verbose=self.verbose, + nports=self.nports, + npoles_cplx=self.npoles_cplx, + least_squares_condition=self.least_squares_condition, + least_squares_row_condition=self.least_squares_row_condition, + least_squares_col_condition=self.least_squares_col_condition, + least_squares_rms_error=self.least_squares_rms_error, + eigenval_condition=self.eigenval_condition, + eigenval_row_condition=self.eigenval_row_condition, + eigenval_col_condition=self.eigenval_col_condition, + eigenval_rms_error=self.eigenval_rms_error, + fit_constant=self.fit_constant, + fit_proportional=self.fit_proportional, + dc_enforce=self.dc_enforce, + passivity_enforce=self.passivity_enforce, + denominator=self.model_instance.denominator, + allow_pickle=True + ) + if self.verbose: + print(f"Model parameters saved to {filename}") + + def _write_model_instance(self,filename): + assert self.model_instance is not None ,"Please run levi() and sk_iteration() first." + with open(filename,"wb") as f: + pickle.dump(self.model_instance,f) + if self.verbose: + print(f"Model instance saved to {filename}") + + def _load_model_instance(self,filename): + with open(filename,"rb") as f: + self.model_instance = pickle.load(f) + if self.verbose: + print(f"Model instance loaded from {filename}") + return self.model_instance + + @classmethod + def _load_npz(cls,filename): + instance = cls( + npoles_cplx=1, # 临时值，稍后会被覆盖 + freqs=np.array([]), # 临时值，稍后会被覆盖 + H=np.array([[]]), # 临时值，稍后会被覆盖 + model=MultiPortOrthonormalBasis, # 临时值，稍后会被覆盖 + iterations=1, # 临时值，稍后会被覆盖 + verbose=False # 临时值，稍后会被覆盖 + ) + data = np.load(filename,allow_pickle=True) + instance.freqs = data["freqs"] + instance.H = data["H"] + instance.iterations = int(data["iterations"]) + instance.verbose = bool(data["verbose"]) + instance.nports = int(data["nports"]) + instance.npoles_cplx = int(data["npoles_cplx"]) + instance.least_squares_condition = data["least_squares_condition"].tolist() + instance.least_squares_row_condition = data["least_squares_row_condition"].tolist() + instance.least_squares_col_condition = data["least_squares_col_condition"].tolist() + instance.least_squares_rms_error = data["least_squares_rms_error"].tolist() + instance.eigenval_condition = data["eigenval_condition"].tolist() + instance.eigenval_row_condition = data["eigenval_row_condition"].tolist() + instance.eigenval_col_condition = data["eigenval_col_condition"].tolist() + instance.eigenval_rms_error = data["eigenval_rms_error"].tolist() + instance.fit_constant = bool(data["fit_constant"]) + instance.fit_proportional = bool(data["fit_proportional"]) + instance.dc_enforce = bool(data["dc_enforce"]) + instance.passivity_enforce = bool(data["passivity_enforce"]) + poles = data["poles"] + denominator = data["denominator"] + instance.model = globals()[data["model_name"].item()] + if instance.verbose: + print(f"Model parameters loaded from {filename}") + return instance + + + def fit(self): + self.levi() + self.model_instance = self.sk_iteration() + return self.model + + def levi(self): + self.poles = generate_starting_poles(self.npoles_cplx,beta_min=max(self.freqs[0],1e4),beta_max=self.freqs[-1]*1.1,verbose=self.verbose) + self.model_instance=self.model( + H=self.H, + freqs=self.freqs, + pre_poles=self.poles, + fit_constant=self.fit_constant, + fit_proportional=self.fit_proportional, + dc_enforce=self.dc_enforce, + passivity_enforce=self.passivity_enforce + ) + return self.model_instance + + def sk_iteration(self): + for i in range(self.iterations): + assert self.model_instance is not None ,"Please run levi() first." + self.poles = self.model_instance.poles + self.model_instance = self.model( + H=self.H, + freqs=self.freqs, + pre_poles=self.poles, + fit_constant=self.fit_constant, + fit_proportional=self.fit_proportional, + dc_enforce=self.dc_enforce, + passivity_enforce=self.passivity_enforce + ) + if self.verbose: + print(f"Iteration {i+1}/{self.iterations}") + print("A:",self.model_instance.A) + print("B:",self.model_instance.B) + print("C:",self.model_instance.C) + print("D:",self.model_instance.D) + print("poles:",self.model_instance.poles) + print("denominator:",self.model_instance.denominator) + + least_squares_condition,\ + least_squares_row_condition,\ + least_squares_col_condition,\ + least_squares_rms_error = self.model_instance.least_squares_metric + eigenval_condition,\ + eigenval_row_condition,\ + eigenval_col_condition,\ + eigenval_rms_error = self.model_instance.eigen_metric + self.least_squares_condition.append(least_squares_condition) + self.least_squares_row_condition.append(least_squares_row_condition) + self.least_squares_col_condition.append(least_squares_col_condition) + self.least_squares_rms_error.append(least_squares_rms_error) + self.eigenval_condition.append(eigenval_condition) + self.eigenval_row_condition.append(eigenval_row_condition) + self.eigenval_col_condition.append(eigenval_col_condition) + self.eigenval_rms_error.append(eigenval_rms_error) + return self.model_instance + + def plot_metrics(self,show:bool=True,save_path=None): + plt.figure(figsize=(16, 12)) + plt.subplot(4, 2, 1) + plt.plot( + range(1, len(self.least_squares_condition) + 1), + self.least_squares_condition, + label='Least Squares Condition' + ) + plt.legend() + plt.subplot(4, 2, 2) + plt.plot( + range(1, len(self.least_squares_row_condition) + 1), + self.least_squares_row_condition, + label='Least Squares Row Condition' + ) + plt.legend() + + plt.subplot(4, 2, 3) + plt.plot( + range(1, len(self.least_squares_col_condition) + 1), + self.least_squares_col_condition, + label='Least Squares Col Condition' + ) + plt.legend() + + plt.subplot(4, 2, 4) + plt.plot( + range(1, len(self.least_squares_rms_error) + 1), + self.least_squares_rms_error, + label='Least Squares RMS Error' + ) + plt.legend() + + plt.subplot(4, 2, 5) + plt.plot( + range(1, len(self.eigenval_condition) + 1), + self.eigenval_condition, + label='Eigenvalue Condition' + ) + plt.legend() + + plt.subplot(4, 2, 6) + plt.plot( + range(1, len(self.eigenval_row_condition) + 1), + self.eigenval_row_condition, + label='Eigenvalue Row Condition' + ) + plt.legend() + + plt.subplot(4, 2, 7) + plt.plot( + range(1, len(self.eigenval_col_condition) + 1), + self.eigenval_col_condition, + label='Eigenvalue Col Condition' + ) + plt.legend() + + plt.subplot(4, 2, 8) + plt.plot( + range(1, len(self.eigenval_rms_error) + 1), + self.eigenval_rms_error, + label='Eigenvalue RMS Error' + ) + plt.legend() + + if show: + plt.show() + if save_path is not None: + if self.verbose: + print(f"Saving metrics plot to {save_path}/fitting_metrics.png") + os.makedirs(save_path, exist_ok=True) + plt.savefig(f"{save_path}/fitting_metrics.png") + + def plot_model_responses(self,show:bool=True,save_path=None): + assert self.model_responses_freqs is not None and self.model_responses_H is not None, "Please run get_model_responses() first." + for i in range(self.nports): + for j in range(self.nports): + plt.figure(figsize=(12, 6)) + plt.subplot(2, 2, 1) + plt.plot(self.freqs, np.abs(self.H[:,i,j]), 'o', ms=4, color='red', label='Input Samples') + plt.plot(self.model_responses_freqs, np.abs(self.model_responses_H[:,i,j]), '-', lw=2, color='k', label='Fit') + plt.title(f"Response i={i+1}, j={j+1}") + plt.ylabel("Magnitude") + plt.legend(loc="best") + plt.subplot(2, 2, 2) + plt.plot(self.freqs, np.angle(self.H[:,i,j],deg=True), 'o', ms=4, color='red', label='Input Samples') + plt.plot(self.model_responses_freqs, np.angle(self.model_responses_H[:,i,j],deg=True), '-', lw=2, color='k', label='Fit') + plt.title(f"Response i={i+1}, j={j+1}") + plt.ylabel("Phase (deg)") + plt.legend(loc="best") + plt.tight_layout() + plt.subplot(2, 2, 3) + plt.plot(self.freqs, np.real(self.H[:,i,j]), 'o', ms=4, color='red', label='Input Samples') + plt.plot(self.model_responses_freqs, np.real(self.model_responses_H[:,i,j]), '-', lw=2, color='k', label='Fit') + plt.title(f"Response i={i+1}, j={j+1}") + plt.ylabel("Real Part") + plt.legend(loc="best") + plt.subplot(2, 2, 4) + plt.plot(self.freqs, np.imag(self.H[:,i,j]), 'o', ms=4, color='red', label='Input Samples') + plt.plot(self.model_responses_freqs, np.imag(self.model_responses_H[:,i,j]), '-', lw=2, color='k', label='Fit') + plt.title(f"Response i={i+1}, j={j+1}") + plt.ylabel("Imag Part") + plt.legend(loc="best") + plt.tight_layout() + if show: + plt.show() + if save_path is not None: + if self.verbose: + print(f"Saving response plot for port {i+1},{j+1} to {save_path}/response_{i+1}_{j+1}.png") + os.makedirs(save_path, exist_ok=True) + plt.savefig(f"{save_path}/response_{i+1}_{j+1}.png") + + def get_model_responses(self,freqs): + assert self.model_instance is not None ,"Please run levi() and sk_iteration() first." + self.model_responses_freqs = freqs + self.model_responses_H = self.model_instance.get_model_responses(freqs) + return self.model_responses_H \ No newline at end of file diff --git a/core/__init__.py b/ovf/core/__init__.py similarity index 100% rename from core/__init__.py rename to ovf/core/__init__.py diff --git a/core/basis/MultiPortOrthonormalBasis.py b/ovf/core/basis/MultiPortOrthonormalBasis.py similarity index 75% rename from core/basis/MultiPortOrthonormalBasis.py rename to ovf/core/basis/MultiPortOrthonormalBasis.py index e35144a..3dc8b58 100644 --- a/core/basis/MultiPortOrthonormalBasis.py +++ b/ovf/core/basis/MultiPortOrthonormalBasis.py @@ -1,18 +1,18 @@ import numpy as np import skrf as rf -from ..utils import cond_row_inf, cond_col_one, generate_starting_poles +from ovf.core.utils import cond_row_inf, cond_col_one, generate_starting_poles class MultiPortOrthonormalBasis: - def __init__(self,H,freqs,poles,weights=None,passivity_enforce=True,dc_enforce=True,fit_constant=True,fit_proportional=False): - self.least_squares_condition = None - self.least_squares_row_condition = None - self.least_squares_col_condition = None - self.least_squares_rms_error = None - self.eigenval_condition = None - self.eigenval_row_condition = None - self.eigenval_col_condition = None - self.eigenval_rms_error = None + def __init__(self,H,freqs,pre_poles,passivity_enforce=True,dc_enforce=True,fit_constant=True,fit_proportional=False): + self._least_squares_condition = None + self._least_squares_row_condition = None + self._least_squares_col_condition = None + self._least_squares_rms_error = None + self._eigenval_condition = None + self._eigenval_row_condition = None + self._eigenval_col_condition = None + self._eigenval_rms_error = None self.Cr = None self.dc_tol = 1e-18 @@ -20,69 +20,109 @@ class MultiPortOrthonormalBasis: self.dc_enforce = dc_enforce self.fit_constant = fit_constant self.fit_proportional = fit_proportional + self.passivity_enforce = passivity_enforce self.freqs = freqs self.H = H self.ports = H.shape[1] self.s = self.freqs * 2j * np.pi - self.P = len(poles) - self.poles = poles + self.pre_poles = pre_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.C,self.w0,self.e = self.fit_denominator(self.H, weights=weights) + self.Phi = self.generate_basis(self.s, self.pre_poles) + self.A = self.matrix_A(self.pre_poles) + self.B = self.vector_B(self.pre_poles) + self.C,self.w0,self.e = self.fit_denominator(self.H) self.D = self.w0 + + @property + def constant(self): + return self.w0 + + @property + def proportional(self): + return self.e - self.residuals = self.C / np.sqrt(2 * np.real(-np.array(self.poles))) - + @property + def poles(self): z = np.linalg.eigvals(self.A - self.B @ self.C) - if passivity_enforce: - self.next_poles = self.passivity_enforce(z) + if self.passivity_enforce: + poles = self._passivity_enforce(z) else: - self.next_poles = z + poles = z - self.eigenval_condition,\ - self.eigenval_row_condition,\ - self.eigenval_col_condition,\ - self.eigenval_rms_error = self.eigen_metric() - - self.Dt = self.eval_Dt_state_space() - self.Dt_Dt_1 = np.linalg.norm(self.Dt) / np.linalg.norm(weights) if weights is not None else np.linalg.norm(self.Dt) - pass + return poles + + @property + def least_squares_metric(self): + return self._least_squares_condition,\ + self._least_squares_row_condition,\ + self._least_squares_col_condition,\ + self._least_squares_rms_error + @property def eigen_metric(self): + self._eigenval_condition,\ + self._eigenval_row_condition,\ + self._eigenval_col_condition,\ + self._eigenval_rms_error = self._eigen_metric() + return self._eigenval_condition, self._eigenval_row_condition, self._eigenval_col_condition, self._eigenval_rms_error + + + @property + def residuals(self): + C = self.C / np.sqrt(2 * np.real(-np.array(self.pre_poles))) + return C + + @property + def numerator(self): + if self.Cr is None: + self.Cr = self.non_bias_Cr(w0=self.w0) + phi = self.Phi + num = np.zeros((len(self.freqs),self.ports,self.ports),dtype=complex) + for i in range(self.ports): + for j in range(self.ports): + num[:,i,j] = (phi @ self.Cr[i][j]).reshape(-1) + return num + + @property + def denominator(self): + phi = self.Phi + den = self.w0 + phi @ self.residuals.T + return den + + def _eigen_metric(self): """Return condition number and RMS error of eigenvalues of A-BC.""" - z = np.linalg.eigvals(self.A - self.B @ self.C) + z = self.poles cond = np.linalg.cond(self.A - self.B @ self.C) - rms = np.sqrt(np.mean(np.abs(np.real(z) - np.real(self.poles))**2 + np.abs(np.imag(z) - np.imag(self.poles))**2)) + rms = np.sqrt(np.mean(np.abs(np.real(z) - np.real(self.pre_poles))**2 + np.abs(np.imag(z) - np.imag(self.pre_poles))**2)) row_cond = cond_row_inf(self.A - self.B @ self.C) col_cond = cond_col_one(self.A - self.B @ self.C) return cond,row_cond,col_cond,rms - def least_squares_metric(self,A,b): + def _least_squares_metric(self,Q,R,x,b): """Return condition number and RMS error of least-squares matrix A and rhs b.""" - cond = np.linalg.cond(A) - rms = np.sqrt(np.mean((A @ np.linalg.pinv(A) @ b - b)**2)) + cond = np.linalg.cond(R) + residual = Q @ R @ x -b + rms = np.sqrt(np.mean(residual**2)) - row_cond = cond_row_inf(A) - col_cond = cond_col_one(A) + row_cond = cond_row_inf(R) + col_cond = cond_col_one(R) return cond,row_cond,col_cond,rms - def passivity_enforce(self,poles): + def _passivity_enforce(self,poles): """enforce poles' real parts to be negative""" enforced_poles = [] for pole in poles: if pole.real > 0: pole = -np.conj(pole) enforced_poles.append(pole) - return enforced_poles + return np.array(enforced_poles).reshape(-1) def eval_Dt_state_space(self): """Return D(s_k)=C(s_k I - A)^(-1)B + D for all k (complex 1D array).""" @@ -167,7 +207,7 @@ class MultiPortOrthonormalBasis: def vector_B(self, poles): return np.ones((len(poles), 1), float) - def fit_denominator(self, H, weights=None, d0 = 1.0): + def fit_denominator(self, H, d0 = 1.0): """ Solve formula (70) on the real basis Φ to obtain: - d (real) → packs into C for this state's block structure @@ -175,7 +215,6 @@ class MultiPortOrthonormalBasis: Optional 'weights' (K,) apply row scaling: SK weighting if 1/|D_prev|. """ K, N = self.Phi.shape - one = np.ones((K, 1), np.complex128) Phi = self.Phi dc_tol = 1e-18 has_dc = self.dc_enforce and self.freqs[0] < dc_tol @@ -189,6 +228,7 @@ class MultiPortOrthonormalBasis: keep[k0] = False if self.fit_constant: + one = np.ones((K, 1), np.complex128) Phi_w = np.hstack([one, Phi]) index = 0 M_kp = None @@ -212,31 +252,15 @@ class MultiPortOrthonormalBasis: M_kp = M0 if M_kp is None else np.vstack([M_kp, M0]) assert M_kp is not None - if weights is None: - weights_kp = np.diag(np.ones(len(self.freqs[keep]) * self.ports**2, np.complex128)) - else: - weights_kp = np.diag(np.ones(len(self.freqs[keep]) * self.ports**2, np.complex128)) - # weights_kp0 = weights[keep] - # weights0 = [] - # for i in range(self.ports **2 ): - # for res in weights_kp0: - # weights0.append(1/res) - # weights_kp = np.diag(np.array(weights0)) - if has_dc: - M_w_kp = weights_kp @ M_kp - A_re = np.real(M_w_kp) - A_im = np.imag(M_w_kp) + A_re = np.real(M_kp) + A_im = np.imag(M_kp) mask = np.ones(K, dtype=bool); mask[k0] = False - # exact (unweighted) DC rows: - # A_dc_re = np.real(M_kp).reshape(1, -1) - # A_dc_im = np.imag(M_kp).reshape(1, -1) + else: - M_w_kp = weights_kp @ M_kp - A_re = np.real(M_w_kp) - A_im = np.imag(M_w_kp) - # A_dc_re = A_dc_im = None + A_re = np.real(M_kp) + A_im = np.imag(M_kp) A_blocks = [A_re, A_im] @@ -273,7 +297,7 @@ class MultiPortOrthonormalBasis: H_0 = H[:,i,j][keep] H_kp = H_0 if H_kp is None else np.hstack([H_kp, H_0]) assert H_kp is not None - H_kp = weights_kp @ H_kp.reshape(-1,1) + H_kp = H_kp.reshape(-1,1) b_re = np.real(d0 * H_kp) b_im = np.imag(d0 * H_kp) @@ -300,18 +324,18 @@ class MultiPortOrthonormalBasis: x = np.linalg.solve(R22, Q2.T @ b) # diagnostics - resid = Q2 @ R22 @ x - b + # resid = Q2 @ R22 @ x - b # self.least_squares_rms_error = float(np.sqrt(np.mean(resid**2))) # self.least_squares_condition = float(np.linalg.cond(R)) - self.least_squares_condition,\ - self.least_squares_row_condition,\ - self.least_squares_col_condition,\ - self.least_squares_rms_error = self.least_squares_metric(A, b) + self._least_squares_condition,\ + self._least_squares_row_condition,\ + self._least_squares_col_condition,\ + self._least_squares_rms_error = self._least_squares_metric(Q2,R22,x,b) return self.extract_C_d_e(x,N,d0) def extract_C_d_e(self,C,N,d0=1.0): - a = np.sqrt(2 * np.real(-np.array(self.poles))) + a = np.sqrt(2 * np.real(-np.array(self.pre_poles))) if self.fit_proportional and self.fit_constant: d = C[1] e = C[0] @@ -347,7 +371,7 @@ class MultiPortOrthonormalBasis: def get_model_responses(self,freqs): H = np.zeros((len(freqs),self.ports,self.ports),dtype=complex) s = 1j * 2*np.pi * np.asarray(freqs, float).ravel() - phi = self.generate_basis(s, self.poles) + phi = self.generate_basis(s, self.pre_poles) den = self.w0 + phi @ self.residuals.T if self.Cr is None: self.Cr = self.non_bias_Cr(w0=self.w0) diff --git a/core/basis/__init__.py b/ovf/core/basis/__init__.py similarity index 100% rename from core/basis/__init__.py rename to ovf/core/basis/__init__.py diff --git a/core/sample.py b/ovf/core/sample.py similarity index 100% rename from core/sample.py rename to ovf/core/sample.py diff --git a/core/utils.py b/ovf/core/utils.py similarity index 88% rename from core/utils.py rename to ovf/core/utils.py index dbeaabe..4859462 100644 --- a/core/utils.py +++ b/ovf/core/utils.py @@ -14,7 +14,7 @@ def cond_col_one(A, use_pinv=True): Ainv = np.linalg.pinv(A) if (use_pinv or A.shape[0] != A.shape[1]) else np.linalg.inv(A) return np.linalg.norm(A, ord=1) * np.linalg.norm(Ainv, ord=1) -def generate_starting_poles(n_pairs: int, beta_min: float, beta_max: float, alpha_scale: float = 0.01): +def generate_starting_poles(n_pairs: int, beta_min: float, beta_max: float, alpha_scale: float = 0.01,verbose: bool = False) -> List[complex]: """ 仅生成复共轭对: p = -alpha + j beta, p*。 n_pairs: 复对数量 (总极点数 = 2*n_pairs) @@ -28,5 +28,6 @@ def generate_starting_poles(n_pairs: int, beta_min: float, beta_max: float, alph alpha = alpha_scale * b p = -alpha + 1j * b poles += [p, np.conj(p)] - print(f"生成 {len(poles)} 个初始极点 (复对) {poles}]") + if verbose: + print(f"生成 {len(poles)} 个初始极点 (复对) {poles}]") return poles diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..ec3563e --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,67 @@ +[build-system] +requires = ["setuptools>=61.0", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "ovf" +version = "0.1.5" +description = "A package for orthonormal basis and rational function fitting" +authors = [ + {name = "M4yGem1ni", email = "M4yGem1ni@outlook.com"} +] +readme = "README.md" +requires-python = ">=3.12" + +classifiers = [ + "Development Status :: 3 - Alpha", + "Intended Audience :: Developers", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.12", + "Programming Language :: Python :: 3.13", +] + +dependencies = [ + "matplotlib", + "requests", + "pandas", + "numpy", + "pydantic", + "sympy", + "scikit-rf", + "setuptools" +] + +[project.urls] +Homepage = "https://github.com/M4yGem1ni/your-repo" + +[tool.black] +line-length = 88 +target-version = ['py312'] +include = '\.pyi?$' +extend-exclude = ''' +/( + # directories + \.eggs + | \.git + | \.hg + | \.mypy_cache + | \.tox + | \.venv + | build + | dist +)/ +''' + +[tool.setuptools] +include-package-data = true + +[tool.setuptools.packages.find] +include = ["ovf*"] + +# [tool.setuptools.package-data] +# 包含 json 数据文件、md 文档等到最终包中 +# "ovf.examples.data" = ["datasets/**/*"] +# "docs" = ["*.md"] + +# 如果 test 不是包的一部分可以不用包含 +# "test" = ["*.py"] \ No newline at end of file diff --git a/test/test_numpy.py b/test/test_numpy.py deleted file mode 100644 index 5e03142..0000000 --- a/test/test_numpy.py +++ /dev/null @@ -1,9 +0,0 @@ -import numpy as np - -# 创建一个复数矩阵 -A = np.array([[1+2j, 2-1j], [3+4j, 4+0j]]) - -Q, R = np.linalg.qr(A) - -print("Q =\n", Q) -print("R =\n", R) \ No newline at end of file diff --git a/test/test_orthonormal_basis.py b/test/test_orthonormal_basis.py deleted file mode 100644 index c5116b4..0000000 --- a/test/test_orthonormal_basis.py +++ /dev/null @@ -1,102 +0,0 @@ -import numpy as np -from core.orthonormal_basis import generate_muntz_laguerre_basis - -# ---------- 可选：离散再正交 (加权 QR) ---------- -# def trapezoid_weights(freqs: np.ndarray): -# if len(freqs) == 1: -# return np.ones(1) -# df = np.diff(freqs) -# w = np.zeros_like(freqs, dtype=float) -# w[0] = 0.5 * df[0] -# w[-1] = 0.5 * df[-1] -# if len(freqs) > 2: -# w[1:-1] = 0.5 * (df[:-1] + df[1:]) -# return w - -def weighted_qr_from_basis(basis_cols: list[np.ndarray], weights: np.ndarray | None = None): - A = np.column_stack(basis_cols) # (Nf, M) - if weights is None: - sw = np.ones(A.shape[0]) - else: - sw = np.sqrt(weights.real) - Aw = sw[:, None] * A - Qw, R = np.linalg.qr(Aw) - Phi = Qw / (sw[:, None] + 1e-30) - return Phi, R # Raw = Phi R - -def check_discrete_orthogonality(Phi: np.ndarray, w: np.ndarray): - G = Phi.conj().T @ (w[:, None] * Phi) - off = np.max(np.abs(G - np.eye(G.shape[0]))) - return G, off - -def verify_orthonormal(Phi: np.ndarray, w: np.ndarray, atol=1e-10, rtol=1e-8): - """ - 返回: - G : Gram 矩阵 (Φ^H W Φ) - max_off : 最大非对角幅值 - diag_err : max |diag(G)-1| - passed : 是否满足阈值 - """ - G = Phi.conj().T @ (w[:, None] * Phi) - I = np.eye(G.shape[0]) - diag_err = np.max(np.abs(np.diag(G) - 1.0)) - max_off = np.max(np.abs(G - I + np.diag(np.diag(G) - 1.0))) - passed = (diag_err <= atol) and (max_off <= rtol) - return G, max_off, diag_err, passed - -def omega_weights(freqs_hz: np.ndarray): - """ - 基于 ω=2πf 的梯形法得到 w_ω = Δω/(2π)，使得 - (1/2π) ∫_{-∞}^{∞} → Σ w_ω,k - """ - f = freqs_hz - if len(f) == 1: - return np.ones(1) - df = np.diff(f) - w_f = np.zeros_like(f) - w_f[0] = 0.5 * df[0] - w_f[-1] = 0.5 * df[-1] - if len(f) > 2: - w_f[1:-1] = 0.5 * (df[:-1] + df[1:]) - # dω = 2π df, (1/2π) * dω = df ⇒ 直接 w_f 就是 w_ω - return w_f # 已等价于 Δω/(2π) - -def evaluate(basis,freqs): - w = omega_weights(freqs) - Phi_num, R = weighted_qr_from_basis(basis, w) - Gram, off = check_discrete_orthogonality(Phi_num, w) - print("离散 Gram 最大非对角元素 =", off) - print("R 形状:", R.shape) - # 验证 Raw ≈ Phi R - raw = np.column_stack(basis) - err = np.max(np.abs(raw - Phi_num @ R)) - print("重构误差 ||Raw - Phi R||_∞ =", err) - - # 验证正交性 - print("离散 Gram 矩阵 (前5x5):") - print(Gram[:5, :5]) - - Gcheck, max_off, diag_err, ok = verify_orthonormal(Phi_num, w) - print(f"Diag 误差={diag_err:.3e}, Max off={max_off:.3e}, Orthonormal={ok}") - # 额外: 验证 R - # raw = Φ R => R ≈ Φ^H W raw (因为 Φ^H W Φ = I) - R_alt = Phi_num.conj().T @ (w[:,None] * raw) - print("R 差异 ||R - R_alt||_max =", np.max(np.abs(R - R_alt))) - -# ------------------ 示例 ------------------ -if __name__ == "__main__": - # 示例稳定极点 (复对正虚部在前) - init_poles = [ - -1.0e3 + 2.5e9j, - -1.0e3 - 2.5e9j, - ] - freqs = np.linspace(1e8, 8e9, 40) - s = 1j * 2 * np.pi * freqs - - basis = generate_muntz_laguerre_basis(s, init_poles) - print("解析基函数数量 =", len(basis)) - print("基函数:") - for i in range(len(basis)): - print(f"φ_{i}:", basis[i][:]) - - evaluate(basis, freqs) diff --git a/test/test_vector_fitting.py b/test/test_vector_fitting.py deleted file mode 100644 index 5c64569..0000000 --- a/test/test_vector_fitting.py +++ /dev/null @@ -1,60 +0,0 @@ -import matplotlib.pyplot as plt -from skrf import Network -from core.freqency import auto_select_multple_ports -from core.freqency import auto_select -from core.vectorFitting import VectorFitting -import numpy as np -from typing import Literal -import skrf as rf - -def vector_fitting_info(vf:VectorFitting,sampled_freqs,sampled_responses,title:str='vf',parameter_type:Literal['y','s','z']='y'): - for i in range(vf.network.nports): - for j in range(vf.network.nports): - rms_error = vf.get_rms_error(i,j,parameter_type=parameter_type) - print(f'RMS Error Port {i+1} to Port {j+1}: {rms_error}') - fitted_points = vf.get_model_response(i,j,sampled_freqs) - if parameter_type=='y': - input_points = vf.network.y[:,i,j] - elif parameter_type=='s': - input_points = vf.network.s[:,i,j] - elif parameter_type=='z': - input_points = vf.network.z[:,i,j] - plt.figure(figsize=(20,12)) - plt.suptitle(f'{title} Port {i+1} to Port {j+1}') - plt.subplot(2,2,1) - plt.title('Magnitude') - plt.plot(sampled_freqs,np.abs(sampled_responses[:,i,j]),'o', ms=4, color='red', label='Samples') - plt.plot(vf.network.f,np.abs(input_points),'x', ms=4, color='blue', label='Input Samples') - plt.plot(sampled_freqs,np.abs(fitted_points),'-', lw=2, color='k', label='Fit') - plt.ylabel('Magnitude (dB)') - plt.grid() - plt.legend() - plt.subplot(2,2,2) - plt.title('Phase') - plt.plot(sampled_freqs,np.angle(sampled_responses[:,i,j]),'o', ms=4, color='red', label='Samples') - plt.plot(vf.network.f,np.angle(input_points),'x', ms=4, color='blue', label='Input Samples') - plt.plot(sampled_freqs,np.angle(fitted_points),'-', lw=2, color='k', label='Fit') - plt.ylabel('Phase (deg)') - plt.xlabel('Frequency (GHz)') - plt.grid() - plt.legend() - plt.subplot(2,2,3) - plt.title('RMS_Error') - plt.plot(rms_error, 'b', label='RMS Error') - plt.ylabel('RMS Error (dB)') - plt.grid() - plt.legend() - plt.savefig(f"outputs/{title.replace(' ','_')}_{i+1}_{j+1}.png") - - -original_network = Network('/tmp/paramer/simulation/3000/3000.s2p') -# original_network = rf.data.ring_slot - -H,freqs = auto_select_multple_ports(original_network.y,original_network.f,max_points=20) -fitted_network = Network(y=H,f=freqs) -vf = VectorFitting(fitted_network) -vf.vector_fit(n_poles_cmplx=2,n_poles_real=0,parameter_type='y') -# vf.auto_fit(parameter_type='y') - -vector_fitting_info(vf,original_network.f,original_network.y,parameter_type='y') -