From 05a9702b88f6acbde064377989cf42b66fe86137 Mon Sep 17 00:00:00 2001 From: mayge Date: Mon, 22 Sep 2025 06:37:28 -0400 Subject: [PATCH] =?UTF-8?q?feat:=20=E6=95=B4=E5=90=88=E4=BA=86basicQR?= =?UTF-8?q?=E5=92=8CrelaxedQR?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- core/relaxed_basisQR.py | 117 +++++++++++++++++++++++++++------------- 1 file changed, 80 insertions(+), 37 deletions(-) diff --git a/core/relaxed_basisQR.py b/core/relaxed_basisQR.py index 75dd73c..bb7b1b0 100644 --- a/core/relaxed_basisQR.py +++ b/core/relaxed_basisQR.py @@ -4,9 +4,10 @@ from scipy.linalg import block_diag import skrf as rf from skrf import VectorFitting from core.freqency 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): + 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 @@ -16,6 +17,7 @@ class RelaxedBasicBasisQR: self.dc_enforce = dc_enforce self.fit_constant = fit_constant + self.fit_proportional = fit_proportional # self.H = H # self.freqs = freqs @@ -27,7 +29,7 @@ class RelaxedBasicBasisQR: 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.fit_denominator(self.H, weights=weights) + self.Cw,self.w0,self.e = self.fit_denominator(self.H, weights=weights) self.D = self.w0 self.Cr = None @@ -119,7 +121,7 @@ class RelaxedBasicBasisQR: B = Bb if B is None else np.vstack([B, Bb]) return B - def fit_denominator(self, H, weights=None): + 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 @@ -129,6 +131,10 @@ class RelaxedBasicBasisQR: 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: @@ -136,17 +142,16 @@ class RelaxedBasicBasisQR: else: weights = np.diag([1/res for res in weights]) - Phi = self.Phi - Phi_w = np.hstack([one, Phi]) - M = np.hstack([Phi, -(H * Phi_w)]) # (K, 2N+1), complex - - dc_tol = 1e-18 - has_dc = self.dc_enforce and self.freqs[0] < dc_tol + 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 = np.ones(K, dtype=bool); keep[k0] = False + keep[k0] = False M_w = weights @ M A_re = np.real(M_w[keep, :]) A_im = np.imag(M_w[keep, :]) @@ -154,49 +159,78 @@ class RelaxedBasicBasisQR: # 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 - 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_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 ---- - A_blocks = [A_re, A_im] + # if A_dc_re is not None: # A_blocks += [A_dc_re, A_dc_im] - - A_blocks += [A_w0] - A = np.vstack(A_blocks).astype(float) - - m = A_re.shape[0] + A_im.shape[0] - b = np.zeros(m, float) - # if A_dc_re is not None: # b = np.concatenate([b, np.zeros(2, float)]) # DC rows → 0 - b = np.concatenate([b, b_w0]) # ---- 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") - x = np.linalg.solve(R, Q.T @ b) + + 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 = A @ 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)) # 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 Cw, w0 + # 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): @@ -215,12 +249,22 @@ class RelaxedBasicBasisQR: 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 = 5 - H11,freqs = auto_select([network.y[i][0][0] for i in range(start_point,len(network.y))],network.f[start_point:],max_points=20) + K = 10 + + full_freqences = network.f[start_point:] + noised_sampled_points = [(network.y[i][0][0]) for i in range(start_point,len(network.y))] + sampled_points = [network.y[i][0][0] 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) @@ -266,10 +310,9 @@ if __name__ == "__main__": fig, axes = plt.subplots(3, 2, figsize=(15, 16), sharex=False) ax00 = axes[0][0] - full_freqences = network.f[start_point:] - sliced_freqences = freqs - sampled_points = [network.y[i][0][0] for i in range(start_point,len(network.y))] 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')