feat: Add QR part

core/basisQR.py

@@ -5,7 +5,7 @@ import skrf as rf
 from skrf import VectorFitting
 from core.freqency import auto_select
 
-class BasicBasis:
+class BasicBasisQR:
     def __init__(self,H,freqs,poles,weights=None):
         self.least_squares_rms_error  = None
         self.least_squares_condition = None
@@ -21,7 +21,8 @@ class BasicBasis:
         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)
+        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
@@ -128,29 +129,40 @@ class BasicBasis:
         b_im = np.imag(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(A.T @ A) @ A.T @ b
+        x = np.linalg.inv(R22) @ (Q2.T @ b)
 
-        self.least_squares_rms_error = np.sqrt(np.mean((A @ x - b)**2))
-        self.least_squares_condition = np.linalg.cond(A)
+        self.least_squares_rms_error = np.sqrt(np.mean((Q2 @ R22 @ x - b)**2))
+        self.least_squares_condition = np.linalg.cond(Q2 @ R22)
 
-        Cn,Cd = self.vector_C(x)
-        return Cn,Cd
+        Cw = self.vector_Cw(x)
+        return Cw
     
-    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 vector_Cw(self,x):
+        return x.T
     
-    def evaluate(self,freqs,poles,Cn,Cd,d0=1.0):
+    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, poles)
-        num = phi @ Cn.T
-        den = d0 + phi @ Cd.T
+        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()
 
@@ -162,7 +174,7 @@ if __name__ == "__main__":
     
     Dt_1 = np.ones((len(freqs),1),np.complex128)
     # Levi step (no weighting):
-    basis = BasicBasis(H11,freqs,poles=poles)
+    basis = BasicBasisQR(H11,freqs,poles=poles)
     Dt = basis.Dt
     poles = basis.next_poles
 
@@ -180,7 +192,7 @@ if __name__ == "__main__":
     eigenval_condition = []
     eigenval_rms_error = []
     for i in range(K):
-        basis = BasicBasis(H11,freqs,poles=poles,weights=Dt)
+        basis = BasicBasisQR(H11,freqs,poles=poles,weights=Dt)
         Dt_1 = Dt
         Dt = basis.Dt
         poles = basis.next_poles
@@ -198,7 +210,7 @@ if __name__ == "__main__":
         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)
+    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)
 
@@ -234,7 +246,7 @@ if __name__ == "__main__":
     ax4.set_ylabel("Magnitude")
     ax4.legend(loc="best")
     fig.tight_layout()
-    plt.savefig(f"basic_basis.png")
+    plt.savefig(f"basic_basis_QR.png")