feat: 打包以及添加数据集

.gitignore

@@ -6,4 +6,9 @@ __pypackages__/
 env/
 .venv/
 .vscode/
-outputs/
+.cache/
+test/
+outputs/
+ovf.egg-info/
+dist/
+build/

core/VFManager.py

@@ -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

core/basis/MultiplePortQR.py

@@ -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")
-
-    
-
-

core/basis/basis.py

@@ -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")
-
-    
-
-

core/basis/basisQR.py

@@ -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")
-
-    
-
-

core/basis/relaxed_basisQR.py

@@ -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")
-
-    
-
-

env

@@ -1,2 +0,0 @@
-export PYTHONPATH=$(pwd)/core:$PYTHONPATH
-export PYTHONPATH=$(pwd)/models:$PYTHONPATH

main.py

@@ -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")

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"]

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

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):
+    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_condition = None
-        self.least_squares_row_condition = None
+        self._least_squares_row_condition = None
-        self.least_squares_col_condition = None
+        self._least_squares_col_condition = None
-        self.least_squares_rms_error  = None
+        self._least_squares_rms_error  = None
-        self.eigenval_condition = None
+        self._eigenval_condition = None
-        self.eigenval_row_condition = None
+        self._eigenval_row_condition = None
-        self.eigenval_col_condition = None
+        self._eigenval_col_condition = None
-        self.eigenval_rms_error = 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.pre_poles = pre_poles
-        self.poles = poles
 
-        self.Phi = self.generate_basis(self.s, self.poles)
+        self.Phi = self.generate_basis(self.s, self.pre_poles)
-        self.A = self.matrix_A(self.poles)
+        self.A = self.matrix_A(self.pre_poles)
-        self.B = self.vector_B(self.poles)
+        self.B = self.vector_B(self.pre_poles)
-        self.C,self.w0,self.e = self.fit_denominator(self.H, weights=weights)
+        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:
+        if self.passivity_enforce:
-            self.next_poles = self.passivity_enforce(z)
+            poles = self._passivity_enforce(z)
         else:
-            self.next_poles = z
+            poles = z
 
-        self.eigenval_condition,\
+        return poles
-            self.eigenval_row_condition,\
+    
-            self.eigenval_col_condition,\
+    @property
-            self.eigenval_rms_error = self.eigen_metric()
+    def least_squares_metric(self):
-
+        return self._least_squares_condition,\
-        self.Dt = self.eval_Dt_state_space()
+            self._least_squares_row_condition,\
-        self.Dt_Dt_1 = np.linalg.norm(self.Dt) / np.linalg.norm(weights) if weights is not None else np.linalg.norm(self.Dt)
+            self._least_squares_col_condition,\
-        pass
+            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)
+        cond = np.linalg.cond(R)
-        rms = np.sqrt(np.mean((A @ np.linalg.pinv(A) @ b - b)**2))
+        residual = Q @ R @ x -b
+        rms = np.sqrt(np.mean(residual**2))
 
-        row_cond = cond_row_inf(A)
+        row_cond = cond_row_inf(R)
-        col_cond = cond_col_one(A)
+        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_kp)
-            A_re  = np.real(M_w_kp)
+            A_im  = np.imag(M_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_kp)
-            A_re = np.real(M_w_kp)
+            A_im = np.imag(M_kp)
-            A_im = np.imag(M_w_kp)
-            # A_dc_re = A_dc_im = None
 
         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_condition,\
-            self.least_squares_row_condition,\
+            self._least_squares_row_condition,\
-            self.least_squares_col_condition,\
+            self._least_squares_col_condition,\
-            self.least_squares_rms_error = self.least_squares_metric(A, b)
+            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)

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

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"]

test/test_numpy.py

@@ -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)

test/test_orthonormal_basis.py

@@ -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)

test/test_vector_fitting.py

@@ -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')
-