algorithm/ovf
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 1 Wiki Activity

chore: 分离了sweep和vf部分,vf部分准备写为包

Browse Source
This commit is contained in:
mayge
2025-09-24 22:18:53 -04:00
parent 8e3472c839
commit ba87d6d784
24 changed files with 684 additions and 17745 deletions
Download Patch File Download Diff File Expand all files Collapse all files

362
core/basis/MultiPortOrthonormalBasis.py Normal file
View File

@@ -0,0 +1,362 @@
import numpy as np
import skrf as rf
from ..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):
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
self.dc_enforce = dc_enforce
self.fit_constant = fit_constant
self.fit_proportional = fit_proportional
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.C,self.w0,self.e = self.fit_denominator(self.H, weights=weights)
self.D = self.w0
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_enforce:
self.next_poles = self.passivity_enforce(z)
else:
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.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 = []
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, 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=float)
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."""
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]
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 = 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:
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_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]
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 = np.column_stack(cols).astype(float)
return A
def vector_B(self, poles):
return np.ones((len(poles), 1), float)
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(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_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))
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)
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)))
if self.fit_proportional and self.fit_constant:
d = C[1]
e = C[0]
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]
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
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.residuals.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 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 = self.w0 + phi @ self.residuals.T
if self.Cr is None:
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

440
core/basis/MultiplePortQR.py Normal file
View File

@@ -0,0 +1,440 @@
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")

0
core/basis/__init__.py Normal file
View File

241
core/basis/basis.py Normal file
View File

@@ -0,0 +1,241 @@
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")

259
core/basis/basisQR.py Normal file
View File

@@ -0,0 +1,259 @@
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")

361
core/basis/relaxed_basisQR.py Normal file
View File

@@ -0,0 +1,361 @@
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")