diff --git a/core/MultiPortOrthonormalBasis.py b/core/MultiPortOrthonormalBasis.py index f33e19c..ff162d4 100644 --- a/core/MultiPortOrthonormalBasis.py +++ b/core/MultiPortOrthonormalBasis.py @@ -7,12 +7,29 @@ from core.freqency import auto_select_multple_ports import matplotlib.pyplot as plt import random as rnd +def cond_row_inf(A, use_pinv=True): + """行条件数 κ∞(A) = ||A||∞ * ||A^{-1}||∞；矩形阵用广义逆。""" + A = np.asarray(A) + 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=np.inf) * np.linalg.norm(Ainv, ord=np.inf) + +def cond_col_one(A, use_pinv=True): + """列条件数 κ1(A) = ||A||1 * ||A^{-1}||1；矩形阵用广义逆。""" + A = np.asarray(A) + 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) + class MultiPortOrthonormalBasis: - 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 + def __init__(self,H,freqs,poles,weights=None,passivity=True,dc_enforce=True,fit_constant=True,fit_proportional=False): self.least_squares_condition = None + self.least_squares_row_condition = None + self.least_squares_col_condition = None + self.least_squares_rms_error = None self.eigenval_condition = None + self.eigenval_row_condition = None + self.eigenval_col_condition = None self.eigenval_rms_error = None + self.Cr = None self.dc_tol = 1e-18 @@ -20,38 +37,60 @@ class MultiPortOrthonormalBasis: 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.C,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 + self.residuals = self.C / np.sqrt(2 * np.real(-np.array(self.poles))) + + z = np.linalg.eigvals(self.A - self.B @ self.C) if passivity: - self.next_poles = self.passivity_enforce(p_next) + self.next_poles = self.passivity_enforce(z) 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.next_poles = z + + self.eigenval_condition,\ + self.eigenval_row_condition,\ + self.eigenval_col_condition,\ + self.eigenval_rms_error = self.eigen_metric() self.Dt = self.eval_Dt_state_space() - self.delta = self.Dt / weights if weights is not None else self.Dt + self.Dt_Dt_1 = np.linalg.norm(self.Dt) / np.linalg.norm(weights) if weights is not None else np.linalg.norm(self.Dt) pass + 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) + 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)) + + 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): + """Return condition number and RMS error of least-squares matrix A and rhs b.""" + cond = np.linalg.cond(A) + rms = np.sqrt(np.mean((A @ np.linalg.pinv(A) @ b - b)**2)) + + row_cond = cond_row_inf(A) + col_cond = cond_col_one(A) + + return cond,row_cond,col_cond,rms + + def passivity_enforce(self,poles): """enforce poles' real parts to be negative""" enforced_poles = [] @@ -64,12 +103,12 @@ class MultiPortOrthonormalBasis: 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) + A = np.asarray(self.A, float); n = A.shape[0] + B = np.asarray(self.B, float).reshape(n, 1) + C = np.asarray(self.C, float).reshape(1, n) D = self.D - I = np.eye(n, dtype=np.complex128) - out = np.empty_like(s, dtype=np.complex128) + I = np.eye(n, dtype=float) + out = np.empty_like(s, dtype=float) for k, sk in enumerate(s): DS = D + (C @ np.linalg.inv(sk*I - A) @ B) out[k] = DS[0, 0] @@ -78,7 +117,14 @@ class MultiPortOrthonormalBasis: def generate_basis(self,s, poles): """Real basis of (15)-(16); returns Φ(s) and a layout for packing C.""" + def all_pass(s,ap_list): + res = 1.0 +0.0j + for ap in ap_list: + res *= (s - np.conj(ap)) / (s + ap) + return res + cols = [] + ap_list = [] i = 0 while i < len(poles): ap = -poles[i] @@ -86,42 +132,56 @@ class MultiPortOrthonormalBasis: raise ValueError("poles must be in the LHP") if i+1 < len(poles) and np.isclose(poles[i+1], np.conj(-ap)): ap1 = - poles[i+1] - phi1 = 1/(s + ap) + 1/(s + ap1) # eq (15)generate_basis - phi2 = 1j*(1/(s + ap) - 1/(s + ap1)) # eq (16) (fixed sign) + phi1 = np.sqrt(2 * ap.real) * all_pass(s, ap_list) * ((s - np.abs(ap))/((s + ap)*(s + ap1))) + phi2 = np.sqrt(2 * ap.real) * all_pass(s, ap_list) * ((s + np.abs(ap))/((s + ap)*(s + ap1))) cols += [phi1, phi2] i += 2 + ap_list.append(ap) + ap_list.append(ap1) else: - cols.append(1/(s + ap)) + basis = np.sqrt(2 * ap.real) * all_pass(s, ap_list) * (1/(s + ap)) + cols.append(basis) i += 1 + ap_list.append(ap) 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 + def A_col(p:np.complex128,index:int): + ap = -p + if abs(ap.imag) < 1e-14: + col = [] + for i in range(index): + col.append(0.0) + col.append(-ap.real) + for i in range(len(poles)-index-1): + col.append(2*(-ap).real) + return np.array([col], float) + else: + col1 = [] + col2 = [] + for i in range(index): + col1.append(0.0) + col2.append(0.0) + col1.append(-ap.real); col2.append(-ap.real - np.abs(ap)) + col1.append(-ap.real + np.abs(ap)); col2.append(-ap.real) + for i in range(len(poles)-index-2): + col1.append(2*(-ap).real) + col2.append(2*(-ap).real) + return np.array([col1, col2], float) + + i = 0 + cols = [] while i < len(poles): p = poles[i] - Ab = A_block(p) + cols.extend(A_col(p,i)) 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) + A = np.column_stack(cols).astype(float) 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 + return np.ones((len(poles), 1), float) def fit_denominator(self, H, weights=None, d0 = 1.0): """ @@ -161,7 +221,7 @@ class MultiPortOrthonormalBasis: M_kp = None for i in range(self.ports): for j in range(self.ports): - M0 = np.zeros((K,N*ports**2),dtype=complex) + 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 @@ -169,7 +229,7 @@ class MultiPortOrthonormalBasis: assert M_kp is not None if weights is None: - weights_kp = np.diag(np.ones(len(freqs[keep]) * self.ports**2, np.complex128)) + weights_kp = np.diag(np.ones(len(self.freqs[keep]) * self.ports**2, np.complex128)) else: weights_kp0 = weights[keep] weights0 = [] @@ -225,10 +285,10 @@ class MultiPortOrthonormalBasis: H_kp = None for i in range(self.ports): for j in range(self.ports): - H_kp0 = weights_kp @ (H[:,i,j]).reshape(1,-1)[keep,:] - H_kp = H_kp0 if H_kp is None else np.hstack([H_kp, H_kp0]) + 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 = H_kp.reshape(-1,1) + H_kp = weights_kp @ H_kp.reshape(-1,1) b_re = np.real(d0 * H_kp) b_im = np.imag(d0 * H_kp) @@ -256,35 +316,40 @@ class MultiPortOrthonormalBasis: # 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)) + # self.least_squares_rms_error = float(np.sqrt(np.mean(resid**2))) + # self.least_squares_condition = float(np.linalg.cond(R)) + self.least_squares_condition,\ + self.least_squares_row_condition,\ + self.least_squares_col_condition,\ + self.least_squares_rms_error = self.least_squares_metric(A, b) - # 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) + return self.extract_C_d_e(x,N,d0) - def extract_Cw_d_e(self,C,N,d0=1.0): + def extract_C_d_e(self,C,N,d0=1.0): + a = np.sqrt(2 * np.real(-np.array(self.poles))) if self.fit_proportional and self.fit_constant: d = C[1] e = C[0] - return C[2:].reshape(1, -1), d, e + C = a * C[2:] + return C.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 + C = a * C[1:] + return C.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 + C = a * C[1:] + return C.reshape(1, -1), d, e else: + C = a * C 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()) + den = np.diag((w0 + self.Phi @ self.residuals.T).ravel()) Cr = [] for i in range(self.ports): Cr.append([]) @@ -294,19 +359,132 @@ class MultiPortOrthonormalBasis: Cr[i].append(Cr_ij) return Cr - def evaluate(self,freqs, w0): + 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) - den = w0 + phi @ self.Cw.T + den = self.w0 + phi @ self.residuals.T if self.Cr is None: - self.Cr = self.non_bias_Cr(w0=w0) + self.Cr = self.non_bias_Cr(w0=self.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 +class VFUtils(): + def __init__(self,npoles_cplx,freqs,H,model=MultiPortOrthonormalBasis,iterations:int=5): + poles = generate_starting_poles(npoles_cplx,beta_min=1e4,beta_max=freqs[-1]*1.1) + self.model=model(H=H,freqs=freqs,poles=poles) + self.freqs=freqs + self.H=H + self.iterations=iterations + + self.nports = H.shape[1] + + 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_responses_freqs = None + self.model_responses_H = None + + def fit(self): + for i in range(self.iterations): + print(f"Iteration {i+1}/{self.iterations}") + poles = self.model.next_poles + weights = self.model.Dt + self.model = self.model.__class__(H=self.H,freqs=self.freqs,poles=poles,weights=weights) + print("A:",self.model.A) + print("B:",self.model.B) + print("C:",self.model.C) + print("D:",self.model.D) + print("next_pozles:",self.model.next_poles) + print("Dt:",self.model.Dt) + print("Dt/Dt_1:",np.linalg.norm(self.model.Dt_Dt_1)) + self.least_squares_condition.append(self.model.least_squares_condition) + self.least_squares_row_condition.append(self.model.least_squares_row_condition) + self.least_squares_col_condition.append(self.model.least_squares_col_condition) + self.least_squares_rms_error.append(self.model.least_squares_rms_error) + self.eigenval_condition.append(self.model.eigenval_condition) + self.eigenval_row_condition.append(self.model.eigenval_row_condition) + self.eigenval_col_condition.append(self.model.eigenval_col_condition) + self.eigenval_rms_error.append(self.model.eigenval_rms_error) + return self.model + + def plot_metrics(self): + 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() + plt.savefig("fit_metrics.png") + + def plot_model_responses(self): + 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() + plt.savefig(f"model_response_{i+1}{j+1}.png") + print(f"Saved model_response_{i+1}{j+1}.png") + + def get_model_responses(self,freqs): + self.model_responses_freqs = freqs + self.model_responses_H = self.model.get_model_responses(freqs) + return self.model_responses_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) @@ -318,7 +496,7 @@ if __name__ == "__main__": network = rf.Network(f"/tmp/paramer/simulation/{id}/{id}.s2p") # network = rf.data.ring_slot ports = network.nports - K = 5 + K = 100 full_freqences = network.f[start_point:] noised_sampled_points = network.y[start_point:,:,:].reshape(-1,ports,ports) @@ -329,111 +507,119 @@ if __name__ == "__main__": 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) + + vf = VFUtils(npoles_cplx=2,freqs=freqs,H=H,model=MultiPortOrthonormalBasis,iterations=K) + model = vf.fit() + vf.plot_metrics() + model_responses = vf.get_model_responses(full_freqences) + vf.plot_model_responses() + + # # 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 + # 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.Cw) - print("D:",basis.D) - print("next_pozles:",basis.next_poles) - print("Dt:",Dt, "norm:",np.linalg.norm(Dt)) + # 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.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) + # # 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.evaluate(full_freqences, w0=basis.w0) - fitted_points = H_evaluated - sliced_freqences = freqs + # # 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") + # 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") + # 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") + # # 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") + # # 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") + # 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") + # 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") + # 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") + # ax21 = axes[2][1] + # ax21.plot(eigenval_rms_error, label='Eigenvalue RMS Error') + # ax21.set_title("eigenval_rms_error") + # ax21.set_ylabel("Magnitude") + # ax21.legend(loc="best") + # fig.tight_layout() + # plt.savefig(f"MultiplePortQR_port_{i+1}{j+1}.png") + # print(f"Saved MultiplePortQR_port_{i+1}{j+1}.png") diff --git a/core/orthonormal_basis.py b/core/orthonormal_basis.py deleted file mode 100644 index 1f64ac7..0000000 --- a/core/orthonormal_basis.py +++ /dev/null @@ -1,53 +0,0 @@ -import numpy as np - -def generate_laguerre_basis(poles:list|np.ndarray, s: np.ndarray): - basis = np.zeros((len(poles)+1,len(s)),dtype=complex) - - product = np.ones(len(s),dtype=complex) - - basis[0] = np.ones(len(s),dtype=complex) # φ_0 = 1 - i = 0 - while i < len(poles): - if np.real(poles[i]) >= 0: - raise ValueError(f"极点必须在左半平面: {poles[i]}") - - # 复对首 (正虚部) - if np.iscomplex(poles[i]) and np.imag(poles[i]) > 0: - if i + 1 >= len(poles): - raise ValueError("复极点缺少共轭") - pc = poles[i + 1] - if not np.isclose(pc, np.conj(poles[i])): - raise ValueError("复极点未按 (p, p*) 顺序排列 (正虚部在前)") - sigma = -np.real(poles[i]) # >0 - scale = np.sqrt(2 * sigma) - r = np.abs(poles[i]) - denom = (s - poles[i]) * (s - pc) - - # 两个基函数 - phi_p = scale * (s - r) / denom * product - phi_pc = scale * (s + r) / denom * product - - # product 先乘 (s + p^*)/(s - p)，再乘 (s + p)/(s - p^*) - product = product * (s + pc) / (s - poles[i]) - product = product * (s + poles[i]) / (s - pc) - - basis[i + 1] = phi_p - basis[i + 2] = phi_pc - i += 2 - continue - - # 复对次 (负虚部) —— 应该被首元素处理，出现表示顺序错误 - if np.iscomplex(poles[i]) and np.imag(poles[i]) < 0: - raise ValueError("检测到负虚部复极点但其共轭尚未处理，请将正虚部成员放在前面。") - - # 实极点 - sigma = -np.real(poles[i]) - if sigma <= 0: - raise ValueError("实极点实部应为负 (稳定)。") - scale = np.sqrt(2 * sigma) - phi = scale / (s - poles[i]) * product - # 更新乘积 - product = product * (s + poles[i]) / (s - poles[i]) - i += 1 - basis[i + 1] = phi - return basis diff --git a/models/basic.py b/models/basic.py index 10f81ad..722209c 100644 --- a/models/basic.py +++ b/models/basic.py @@ -16,7 +16,7 @@ class ModelBasicParametersUnit(BaseModel): range: List[Union[float,int,str,bool]] class ModelBasicInfo(BaseModel): - base_url = "http://localhost:8105/api/v1" + base_url: str = Field(default="http://localhost:8105/api/v1") nports: int = Field(default=2) cell_name: str template_name: str @@ -82,6 +82,7 @@ class ModelBasic: parameters_list.append({parameters_name[i]:res[i] for i in range(len(res))}) for res in parameters_list: + print(f"Simulating with parameters: {res}") request_unit = SimulationRequestUnit( user_id=self.info.user_id, template_name=self.info.template_name, @@ -125,21 +126,26 @@ class ModelBasic: simulation_response_model:SimulationResponseUnit = SimulationResponseUnit(**(response[0])) - time.sleep(0.01) + time.sleep(0.05) response = send_get_request(f"{self.info.base_url}/simulations/input_hash/{simulation_response_model.input_hash}") + assert isinstance(response, dict), "Response is not a dictionary." uuid_model = UuidResponseUnit(**response) - time.sleep(0.01) # Wait for 2 seconds before checking the status again + time.sleep(0.05) # Wait for 2 seconds before checking the status again status = uuid_model.status while status != "completed" and status != "failed": - time.sleep(0.01) # Wait for 2 seconds before checking again + time.sleep(0.05) # Wait for 2 seconds before checking again response = send_get_request(f"{self.info.base_url}/simulations/input_hash/{simulation_response_model.input_hash}") assert isinstance(response, dict), "Response is not a dictionary." - uuid_model = UuidResponseUnit(**response) + try: + uuid_model = UuidResponseUnit(**response) + except Exception as e: + print(f"Error parsing response: {e}") + continue assert response is not None, "No response received from the server." status = uuid_model.status diff --git a/models/mlin.py b/models/mlin.py new file mode 100644 index 0000000..c080f67 --- /dev/null +++ b/models/mlin.py @@ -0,0 +1,45 @@ +from typing import List,Optional +from pydantic import BaseModel, Field +from schemas.paramer import SimulationRequestUnit +from skrf import Network +import re +from models.basic import ModelBasic, ModelBasicParametersUnit, ModelBasicInfo, ModelBasicDatasetUnit + +W = [] +L = [] +i = 15.52 +while i <= 100: + W.append(i) + L.append(i) + i = int(i*1.05*100 + 0.5) / 100.0 + +class Mlin(ModelBasic): + def __init__(self): + super().__init__() + + @property + def info(self) -> ModelBasicInfo: + return ModelBasicInfo( + nports=2, + cell_name="mlin", + template_name="em_interface_compound", + user_id=0, + template_version="" + ) + + @property + def settings(self) -> dict: + return {} + + @property + def parameters(self) -> List[ModelBasicParametersUnit]: + return [ + ModelBasicParametersUnit(name="W",type="number",range=W), + ModelBasicParametersUnit(name="L",type="number",range=L) + ] + + + + + + \ No newline at end of file diff --git a/schemas/paramer.py b/schemas/paramer.py index be0fc37..f0f4563 100644 --- a/schemas/paramer.py +++ b/schemas/paramer.py @@ -21,8 +21,7 @@ class UuidResponseUnit(BaseModel): input_hash:str storage_type: str result_path:str - started_at:str - completed_at:str + completed_at:Optional[str] result_json:Optional[dict] error_message:Optional[str] failure_count:Optional[int] diff --git a/sweep.py b/sweep.py new file mode 100644 index 0000000..7ae768d --- /dev/null +++ b/sweep.py @@ -0,0 +1,13 @@ +from models.mlin import Mlin + +W = [] +L = [] +i = 15.52 +while i <= 100: + W.append(i) + L.append(i) + i = int(i*1.05*100 + 0.5) / 100.0 + +model = Mlin() +model.sweep() +model.export("mlin_sweep") \ No newline at end of file