From ae98a7acb1515c21943a6858a01a2234864ade68 Mon Sep 17 00:00:00 2001 From: mayge Date: Tue, 23 Sep 2025 23:02:51 -0400 Subject: [PATCH] feat: prepare to add orthonormal basis --- core/MultiPortOrthonormalBasis.py | 440 ++++++++++++++++++++++++++++++ 1 file changed, 440 insertions(+) create mode 100644 core/MultiPortOrthonormalBasis.py diff --git a/core/MultiPortOrthonormalBasis.py b/core/MultiPortOrthonormalBasis.py new file mode 100644 index 0000000..f33e19c --- /dev/null +++ b/core/MultiPortOrthonormalBasis.py @@ -0,0 +1,440 @@ +import numpy as np +from core.sk_iter import generate_starting_poles +from scipy.linalg import block_diag +import skrf as rf +from skrf import VectorFitting +from core.freqency import auto_select_multple_ports +import matplotlib.pyplot as plt +import random as rnd + +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 + 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): + ap = -poles[i] + if ap.real < 0: + 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) + cols += [phi1, phi2] + i += 2 + else: + cols.append(1/(s + ap)) + 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*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_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]) + assert H_kp is not None + H_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 = 5 + + 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) + 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 = 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)) + + # 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) + + # 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") + + + +