From c1612cb8a97f56036aa6653277473e2f6409b1ce Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Mon, 28 Sep 2026 19:46:53 +0200 Subject: [PATCH 01/18] Update pyproject.toml Updated env to add sympy and optional python-mumps, not available in Pyodide --- pyproject.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 104e1ff..0ef14f9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,6 +30,8 @@ dependencies = [ "numpy>=1.15", "scipy>=1.2", "scikit-rf>=0.30", + "sympy>=1.14", + "python-mumps>=0.0.6; sys_platform != 'emscripten'", # If a dependency ships native code that can't be installed in Pyodide, # tag it with `; sys_platform != 'emscripten'` so a normal pip install # still pulls it eagerly but micropip in the browser skips it. Pair From 3e56ca8b8375b1c90600cc9295607018796e0d5d Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Mon, 28 Sep 2026 19:48:59 +0200 Subject: [PATCH 02/18] added netlist_to_statespace Added NetlistStateSpace class in pathsim_rf for netlist integration --- src/pathsim_rf/netlist_to_statespace.py | 1068 +++++++++++++++++++++++ 1 file changed, 1068 insertions(+) create mode 100644 src/pathsim_rf/netlist_to_statespace.py diff --git a/src/pathsim_rf/netlist_to_statespace.py b/src/pathsim_rf/netlist_to_statespace.py new file mode 100644 index 0000000..14fefaa --- /dev/null +++ b/src/pathsim_rf/netlist_to_statespace.py @@ -0,0 +1,1068 @@ +""" +netlist_to_statespace.py + +Turns a small SPICE-like netlist of R, L, C, K (mutual inductance / coupling +coefficient), V (independent voltage source) and I (independent current +source) elements into an explicit LTI state-space model + + xdot = A x + B u + y = C x + D u + +via Modified Nodal Analysis (MNA) + elimination of the algebraic +(non-storage) unknowns. This is the standard "state-variable method" used +to derive circuit ODEs by hand, automated in code so it scales to +messy topologies and handles coupled inductors correctly. + +Only LINEAR elements are supported (no diodes/switches/etc - PathSim itself +handles those fine via its algebraic-loop solver and Function blocks, see +the accompanying explanation). +""" + +from __future__ import annotations +import sympy as sp +import numpy as np +from dataclasses import dataclass +from pathlib import Path +import importlib.util +import logging +import re +from typing import Literal +import warnings + +from pathsim.blocks.lti import StateSpace + + +# -------------------------------------------------------------------------- +# Netlist parsing +# -------------------------------------------------------------------------- + +@dataclass +class Element: + kind: str # 'R','L','C','V','I','K' + name: str + n1: str = None + n2: str = None + value: float = None + # for K elements: + l1: str = None + l2: str = None + + +SI_MULTIPLIERS = { + "f": 1e-15, + "p": 1e-12, + "n": 1e-9, + "u": 1e-6, + "µ": 1e-6, + "m": 1e-3, + "": 1.0, + "k": 1e3, + "K": 1e3, + "meg": 1e6, + "Meg": 1e6, + "M": 1e6, + "g": 1e9, + "G": 1e9, + "t": 1e12, + "T": 1e12, +} + +_PARAM_PATTERN = re.compile(r"^\.param\s+([A-Za-z_]\w*)\s*=\s*(.+)$") +_VALUE_PATTERN = re.compile( + r"^\s*([+-]?\d*\.?\d+(?:[eE][+-]?\d+)?)([A-Za-zµ]*)\s*$" +) +_VALID_KINDS = {"R", "L", "C", "V", "I", "K"} +_IGNORED_LINE_STARTS = (".", "*", '"') +ReductionMode = Literal["symbolic", "fast"] +_MUMPS_AVAILABLE = importlib.util.find_spec("mumps") is not None +_LOGGER = logging.getLogger("pathsim.netlist_to_statespace") + + +def _normalize_node(node: str) -> str: + """Normalize ground aliases to canonical node '0'.""" + return "0" if node.strip().lower() in {"0", "gnd"} else node + + +def _get_si_multiplier(suffix: str) -> float: + """Return numeric multiplier for an SI prefix/suffix string.""" + if suffix in SI_MULTIPLIERS: + return SI_MULTIPLIERS[suffix] + if suffix.lower() in SI_MULTIPLIERS: + return SI_MULTIPLIERS[suffix.lower()] + + candidates = sorted(SI_MULTIPLIERS.keys(), key=len, reverse=True) + for key in candidates: + if not key: + continue + if suffix.startswith(key): + return SI_MULTIPLIERS[key] + low_key = key.lower() + if suffix.lower().startswith(low_key): + return SI_MULTIPLIERS[low_key] + raise ValueError(f"Unknown SI prefix: '{suffix}'") + + +def parse_value_with_units(value_str: str, params: dict[str, str] | None = None) -> float: + """ + Parse SPICE-like numeric tokens with SI prefixes and optional unit tails. + + Supported forms include: + - plain/scientific floats: ``1``, ``100.e-3``, ``2.2e6`` + - SI-prefixed tokens: ``4u``, ``10k``, ``3Meg``, ``2.2kOhm``, ``1mH`` + - parameter references: ``{RMAIN}`` when ``params`` is provided + """ + if value_str is None: + raise ValueError("Missing value") + + value = value_str.strip() + try: + return float(value) + except ValueError: + pass + + if value.startswith("{") and value.endswith("}"): + if params is None: + raise ValueError(f"Unknown parameter reference: '{value}'") + name = value[1:-1].strip() + if name not in params: + raise ValueError(f"Unknown parameter: '{name}'") + return parse_value_with_units(params[name], params) + + match = _VALUE_PATTERN.fullmatch(value) + if match is None: + raise ValueError(f"Invalid value format: '{value_str}'") + + number, suffix = match.groups() + return float(number) * _get_si_multiplier(suffix) + + +def _parse_optional_source_value(raw: str, params: dict[str, str]) -> float | None: + """Parse source literal values, returning None for waveform expressions.""" + try: + return parse_value_with_units(raw, params) + except ValueError: + # Many source statements use waveforms (PWL, SIN, EXP, etc.); those + # are runtime waveforms and do not affect linearization here. + return None + + +def _parse_behavioral_source(parts: list[str], line_number: int, line: str) -> Element: + """ + Parse LTspice behavioral source shorthand: + Bx n+ n- I= -> modeled as current source input placeholder + Bx n+ n- V= -> modeled as voltage source input placeholder + """ + if len(parts) < 4: + raise ValueError(f"Malformed behavioral source at line {line_number}: '{line}'") + + name = parts[0] + n1, n2 = _normalize_node(parts[1]), _normalize_node(parts[2]) + expr = " ".join(parts[3:]).strip() + expr_upper = expr.upper() + if expr_upper.startswith("I="): + return Element(kind="I", name=name, n1=n1, n2=n2, value=None) + if expr_upper.startswith("V="): + return Element(kind="V", name=name, n1=n1, n2=n2, value=None) + raise ValueError( + f"Unsupported behavioral source expression at line {line_number}: '{line}'" + ) + + +def parse_netlist(text: str) -> list[Element]: + """ + Parse netlist text into linear element records. + + One element per line, whitespace separated; '#' or ';' starts a comment. + R n1 n2 value + L n1 n2 value + C n1 n2 value + V n+ n- value (value is a placeholder; real waveform comes + from the PathSim Source block at sim time) + I n+ n- value (same) + B n+ n- I= (behavioral current source placeholder) + B n+ n- V= (behavioral voltage source placeholder) + K Lname1 Lname2 k (coupling coefficient, -1<=k<=1) + Node '0' (or 'gnd') is ground. + Deck/meta lines beginning with '.', '*', or a quoted header are ignored. + """ + params: dict[str, str] = {} + raw_element_lines: list[tuple[int, str]] = [] + + for line_number, raw in enumerate(text.splitlines(), start=1): + line = raw.split("#")[0].split(";")[0].strip() + if not line: + continue + match = _PARAM_PATTERN.fullmatch(line) + if match is not None: + params[match.group(1)] = match.group(2).strip() + continue + raw_element_lines.append((line_number, line)) + + elements: list[Element] = [] + for line_number, line in raw_element_lines: + parts = line.split() + if not parts: + continue + if parts[0].startswith(_IGNORED_LINE_STARTS): + continue + + name = parts[0] + kind = name[0].upper() + if kind == "B": + elements.append(_parse_behavioral_source(parts, line_number, line)) + continue + + if kind not in _VALID_KINDS: + raise ValueError( + f"Unsupported element '{name}' at line {line_number}: '{line}'" + ) + + if kind == 'K': + if len(parts) < 4: + raise ValueError(f"Malformed K element at line {line_number}: '{line}'") + elements.append(Element(kind='K', name=name, l1=parts[1], l2=parts[2], + value=parse_value_with_units(parts[3], params))) + else: + if len(parts) < 3: + raise ValueError(f"Malformed element at line {line_number}: '{line}'") + n1, n2 = _normalize_node(parts[1]), _normalize_node(parts[2]) + val_token = " ".join(parts[3:]) if len(parts) > 3 else None + if kind in {"R", "L", "C"}: + if val_token is None: + raise ValueError(f"Missing value for element '{name}' at line {line_number}") + val = parse_value_with_units(val_token, params) + else: + val = _parse_optional_source_value(val_token, params) if val_token else None + elements.append(Element(kind=kind, name=name, n1=n1, n2=n2, value=val)) + return elements + + +def parse_netlist_file(path: str | Path) -> list[Element]: + """Load and parse a netlist file.""" + with open(path, "r", encoding="utf-8") as handle: + return parse_netlist(handle.read()) + + +GND = {'0', 'gnd', 'GND'} + + +# -------------------------------------------------------------------------- +# MNA + symbolic DAE -> explicit state space reduction +# -------------------------------------------------------------------------- + +class CircuitModel: + """ + Build an explicit state-space model from a linear circuit netlist. + + The circuit is first stamped as the modified nodal analysis (MNA) system + + ``E * dz/dt + G * z = B * u``, + + where ``z`` contains node voltages, inductor currents, and ideal voltage + source currents. The MNA unknowns are partitioned into differential states + ``x_d`` and algebraic unknowns ``x_a``. Eliminating ``x_a`` gives + + ``x_a = Phi * x_d + Psi * u`` + + and the final explicit model + + ``dx_d/dt = A * x_d + B_ss * u``. + + Both reduction modes produce the same matrices and state ordering. They + differ only in the arithmetic and linear solver used during elimination. + + Parameters + ---------- + elements: + Output of :func:`parse_netlist` / :func:`parse_netlist_file`. + reduction_mode: + Method used to eliminate the algebraic MNA variables: + + ``"symbolic"`` + Stamp and reduce with SymPy matrices. Matrix inverses and products + retain symbolic arithmetic until the final conversion to floating + point. This is the default and is useful for small circuits, + reference results, and diagnosing rank deficiencies. Its runtime + and memory use can grow quickly on large or densely coupled + circuits. + + ``"fast"`` + Stamp directly into floating-point arrays. Eliminate the algebraic + subsystem with a sparse MUMPS factorization, then solve the + differential mass system numerically. This mode is intended for + large RLC/K netlists and repeated output construction. A singular + algebraic solve is retried with progressively larger diagonal + regularization and emits :class:`RuntimeWarning` when regularization + is used. When ``python-mumps`` is unavailable, construction falls + back to symbolic reduction and logs a warning through PathSim. + + The mode affects model construction only + + Attributes + ---------- + A: + State matrix of the reduced explicit system. + B_ss: + Input matrix of the reduced explicit system. + Phi: + Map from differential states to eliminated algebraic variables. + Psi: + Map from inputs to eliminated algebraic variables. + state_labels: + Labels matching the row and column ordering of the state matrices. + """ + + def __init__( + self, + elements: list[Element], + reduction_mode: ReductionMode = "symbolic", + ): + if reduction_mode not in {"symbolic", "fast"}: + raise ValueError( + f"Unsupported reduction_mode '{reduction_mode}'. " + "Expected 'symbolic' or 'fast'." + ) + if reduction_mode == "fast" and not _MUMPS_AVAILABLE: + _LOGGER.warning( + "Fast netlist reduction requires python-mumps; " + "falling back to symbolic reduction." + ) + reduction_mode = "symbolic" + self.reduction_mode = reduction_mode + self.elements = elements + + self.R = [e for e in elements if e.kind == 'R'] + self.L = [e for e in elements if e.kind == 'L'] + self.C = [e for e in elements if e.kind == 'C'] + self.V = [e for e in elements if e.kind == 'V'] + self.I = [e for e in elements if e.kind == 'I'] + self.K = [e for e in elements if e.kind == 'K'] + self.dipoles_by_name = { + e.name: e for e in elements if e.kind in {"R", "L", "C", "V", "I"} + } + + # ---- nodes ---- + nodes = [] + for e in elements: + if e.kind == 'K': + continue + for n in (e.n1, e.n2): + if n not in GND and n not in nodes: + nodes.append(n) + self.nodes = nodes # ordered list of non-ground node names + self.node_idx = {n: i for i, n in enumerate(nodes)} + self.n_nodes = len(nodes) + + self.L_names = [e.name for e in self.L] + self.L_idx = {name: i for i, name in enumerate(self.L_names)} + self.n_L = len(self.L) + + self.V_names = [e.name for e in self.V] + self.V_idx = {name: i for i, name in enumerate(self.V_names)} + self.n_V = len(self.V) + + # unknown ordering: [node voltages] + [inductor currents] + [V-source currents] + self.n_total = self.n_nodes + self.n_L + self.n_V + + def col_node(n): + return None if n in GND else self.node_idx[n] + + def col_iL(name): + return self.n_nodes + self.L_idx[name] + + def col_iV(name): + return self.n_nodes + self.n_L + self.V_idx[name] + + self._col_node, self._col_iL, self._col_iV = col_node, col_iL, col_iV + numeric = self.reduction_mode == "fast" + + # ---- inductance matrix (with mutual terms) ---- + Lmat = ( + np.zeros((self.n_L, self.n_L), dtype=float) + if numeric + else sp.zeros(self.n_L, self.n_L) + ) + l_values = {e.name: e.value for e in self.L} + for e in self.L: + i = self.L_idx[e.name] + Lmat[i, i] = float(e.value) if numeric else sp.nsimplify(e.value) + for k in self.K: + i, j = self.L_idx[k.l1], self.L_idx[k.l2] + Li = l_values[k.l1] + Lj = l_values[k.l2] + if numeric: + M = float(k.value) * np.sqrt(float(Li) * float(Lj)) + else: + M = sp.nsimplify(k.value) * sp.sqrt(sp.nsimplify(Li) * sp.nsimplify(Lj)) + Lmat[i, j] += M + Lmat[j, i] += M + self.Lmat = Lmat + + # ---- inputs: one column per V source, then one column per I source ---- + self.input_names = self.V_names + [e.name for e in self.I] + self.n_u = len(self.input_names) + + # ---- build E (dynamic) and G (algebraic) matrices, and B (input map) ---- + n = self.n_total + if numeric: + E = np.zeros((n, n), dtype=float) + G = np.zeros((n, n), dtype=float) + B = np.zeros((n, self.n_u), dtype=float) + else: + E = sp.zeros(n, n) + G = sp.zeros(n, n) + B = sp.zeros(n, self.n_u) + + def stamp_G(row, col, val): + if row is not None and col is not None: + G[row, col] += val + + # Resistors: contribute to node KCL rows only + for e in self.R: + a, b = col_node(e.n1), col_node(e.n2) + g = (1.0 / float(e.value)) if numeric else (1 / sp.nsimplify(e.value)) + stamp_G(a, a, g); stamp_G(b, b, g) + stamp_G(a, b, -g); stamp_G(b, a, -g) + + # Capacitors: contribute dv/dt terms to node KCL rows (the E matrix) + for e in self.C: + a, b = col_node(e.n1), col_node(e.n2) + c = float(e.value) if numeric else sp.nsimplify(e.value) + if a is not None: E[a, a] += c + if b is not None: E[b, b] += c + if a is not None and b is not None: + E[a, b] -= c + E[b, a] -= c + + # Inductors: KCL stamp (current unknown enters/leaves nodes) + + # dedicated branch row v_na - v_nb - sum_j Lij * d(iLj)/dt = 0 + for e in self.L: + a, b = col_node(e.n1), col_node(e.n2) + iL_col = col_iL(e.name) + row = iL_col # branch row shares index with its current unknown + stamp_G(a, iL_col, 1) + stamp_G(b, iL_col, -1) + if a is not None: G[row, a] += 1 + if b is not None: G[row, b] -= 1 + # branch eqn: v_na - v_nb - L*d(iL)/dt - sum_j M_ij*d(iLj)/dt = 0 + # => E[row, iLj] = -Lij (note the minus sign!) + i = self.L_idx[e.name] + for j, name_j in enumerate(self.L_names): + Lij = self.Lmat[i, j] + if Lij != 0: + E[row, col_iL(name_j)] -= Lij + + # Voltage sources: KCL stamp + branch row v_na - v_nb = u(t) + for e in self.V: + a, b = col_node(e.n1), col_node(e.n2) + iV_col = col_iV(e.name) + row = iV_col + stamp_G(a, iV_col, 1) + stamp_G(b, iV_col, -1) + if a is not None: G[row, a] += 1 + if b is not None: G[row, b] -= 1 + u_col = self.input_names.index(e.name) + B[row, u_col] = 1 + + # Current sources: pure RHS injection into node KCL rows. + # Convention: positive I flows from n1 -> n2 *through the source*, + # i.e. it delivers current INTO n2 and draws it OUT of n1 from the + # external circuit's point of view. + for e in self.I: + a, b = col_node(e.n1), col_node(e.n2) + u_col = self.input_names.index(e.name) + if a is not None: B[a, u_col] -= 1 + if b is not None: B[b, u_col] += 1 + + self.E, self.G, self.B = E, G, B + + # ---- differential / algebraic partition ---- + diff_rows = [r for r in range(n) if any(E[r, c] != 0 for c in range(n))] + alg_rows = [r for r in range(n) if r not in diff_rows] + self.diff_rows, self.alg_rows = diff_rows, alg_rows + + # sanity: the state columns should be exactly node-voltages that own a + # nonzero E column plus all inductor currents; algebraic columns = rest + diff_cols = sorted(set(c for r in diff_rows for c in range(n) if E[r, c] != 0)) + alg_cols = [c for c in range(n) if c not in diff_cols] + if len(diff_cols) != len(diff_rows) or len(alg_cols) != len(alg_rows): + raise ValueError( + "Circuit is degenerate for this reduction (e.g. an all-capacitor " + "loop, an all-inductor cutset, or a floating node). " + f"diff_rows={len(diff_rows)} diff_cols={len(diff_cols)} " + f"alg_rows={len(alg_rows)} alg_cols={len(alg_cols)}" + ) + self.diff_cols, self.alg_cols = diff_cols, alg_cols + + self._reduce() + self._selected_output_rows: list[tuple[np.ndarray, np.ndarray]] = [] + self.output_labels: list[str] = [] + + def _reduce(self): + """ + Reduce the stamped MNA differential-algebraic system. + + This dispatcher selects the arithmetic backend requested by + :attr:`reduction_mode`. Both implementations populate ``A``, ``B_ss``, + ``Phi``, and ``Psi`` and preserve the same differential-state ordering. + State labels are assigned only after a successful reduction. + """ + if self.reduction_mode == "fast": + self._reduce_fast() + else: + self._reduce_symbolic() + self._set_state_labels() + + def _set_state_labels(self): + """Set labels for the differential state vector in `self.diff_cols` order.""" + dc = self.diff_cols + labels = [] + for c in dc: + if c < self.n_nodes: + labels.append(f"v_{self.nodes[c]}") + else: + labels.append(f"i_{self.L_names[c - self.n_nodes]}") + self.state_labels = labels + self.state_label_idx = {label: i for i, label in enumerate(labels)} + + def _raise_algebraic_singular(self): + raise ValueError( + "Algebraic subsystem is singular. Usually means a loop made " + "purely of ideal voltage sources (and/or 0-ohm shorts), or a " + "node with no DC path to ground." + ) + + def _raise_state_singular(self): + raise ValueError( + "State/mass matrix is singular: this circuit's capacitor " + "voltages (or inductor currents) aren't independent, so the " + "naive 'one state per capacitor-touched node' selection " + "over-counted states. Classic cause: a capacitor whose *both* " + "terminals only reach the rest of the circuit through that " + "same capacitor (no other cap ties either node down " + "independently) -- e.g. 'R1 in a / C1 a b / R2 b 0' with " + "nothing else at a or b. Only one of v_a, v_b is really an " + "independent state there; the other is pinned by KCL. This " + "reduction doesn't do full tree/cotree state selection, so it " + "can't detect that automatically yet. Workarounds: (1) add a " + "TINY STRAY CAPACITANCE FROM ONE OF THE FLOATING NODES TO " + "GROUND (not a resistor -- the degeneracy lives in the " + "capacitor/mass matrix, a parallel resistor doesn't touch it " + "at all). e.g. 'Cstray b 0 1e-15' -- this gives that node's " + "own row independent rank; it adds one extra, extremely fast " + "eigenvalue (a numerical artifact, orders of magnitude faster " + "than your real dynamics) alongside the correct physical " + "pole(s), or (2) hand-pick the true independent capacitor " + "voltage as the state and eliminate the redundant node " + "yourself before building the netlist." + ) + + def _solve_with_mumps(self, mat: np.ndarray, rhs: np.ndarray, context: str) -> np.ndarray: + """ + Solve one or more right-hand sides with one MUMPS factorization. + + Parameters + ---------- + mat: + Square coefficient matrix. + rhs: + Vector or matrix of right-hand sides. Each matrix column is solved + independently while reusing the factorization of ``mat``. + context: + Subsystem name included in errors and regularization warnings. + + Returns + ------- + numpy.ndarray + Solution with the same one- or two-dimensional convention as + ``rhs``. + + Notes + ----- + The unmodified matrix is attempted first. If factorization or solution + fails, the method retries ``mat + eps * I`` for increasing ``eps``. + Regularization permits reduction of numerically singular algebraic + systems, but it perturbs their constraints; every successful + regularized solve therefore emits :class:`RuntimeWarning`. + + MUMPS deliberately remains in its default unsymmetric mode. General + MNA systems may contain voltage-source saddle-point blocks, so they + are not positive definite, and the reduced transient matrix is not + generally symmetric. Selecting MUMPS ``sym=1`` (SPD) or passing + ``symmetric=True`` would therefore impose an invalid matrix property. + """ + try: + import mumps + import scipy.sparse as sps + except ImportError as exc: + raise ImportError( + "Fast reduction mode requires `python-mumps` (module `mumps`) " + "and scipy.sparse to be available." + ) from exc + + rhs_2d = rhs if rhs.ndim == 2 else rhs.reshape((-1, 1)) + if mat.shape[0] != mat.shape[1]: + raise ValueError(f"{context} matrix must be square, got {mat.shape}") + if rhs_2d.shape[0] != mat.shape[0]: + raise ValueError( + f"{context} rhs has incompatible shape {rhs_2d.shape} for matrix {mat.shape}" + ) + + last_error: Exception | None = None + for eps in (0.0, 1e-15, 1e-12, 1e-9, 1e-6): + try: + if eps > 0.0: + mat_reg = mat + np.eye(mat.shape[0], dtype=float) * eps + else: + mat_reg = mat + + ctx = mumps.Context() + ctx.set_matrix(sps.csc_matrix(mat_reg)) + ctx.factor() + + cols = [] + for j in range(rhs_2d.shape[1]): + bj = np.array(rhs_2d[:, j], dtype=float, order="F") + xj = np.array(ctx.solve(bj), dtype=float) + cols.append(xj) + x = np.column_stack(cols) + + if eps > 0.0: + warnings.warn( + f"Fast reduction regularized singular {context} with eps={eps:.1e}", + RuntimeWarning, + ) + return x if rhs.ndim == 2 else x[:, 0] + except (mumps.MUMPSError, ValueError, RuntimeError, TypeError) as exc: + last_error = exc + continue + + raise ValueError( + f"Failed to solve {context} even after regularization attempts." + ) from last_error + + @staticmethod + def _safe_matmul(a: np.ndarray, b: np.ndarray) -> np.ndarray: + """ + Matrix multiply without relying on NumPy BLAS-backed matmul. + This avoids environment-specific native linear algebra crashes. + """ + if a.shape[1] != b.shape[0]: + raise ValueError(f"Incompatible shapes for matmul: {a.shape} and {b.shape}") + out = np.zeros((a.shape[0], b.shape[1]), dtype=float) + for i in range(a.shape[0]): + for k in range(a.shape[1]): + aik = a[i, k] + if aik != 0.0: + out[i, :] += aik * b[k, :] + return out + + def _solve_with_sympy_numeric(self, mat: np.ndarray, rhs: np.ndarray, context: str) -> np.ndarray: + """Solve linear system using numeric SymPy LU; supports regularization ladder.""" + rhs_2d = rhs if rhs.ndim == 2 else rhs.reshape((-1, 1)) + m = sp.Matrix(mat) + r = sp.Matrix(rhs_2d) + + try: + sol = m.LUsolve(r) + x = np.array(sol, dtype=float) + return x if rhs.ndim == 2 else x[:, 0] + except (ValueError, sp.matrices.exceptions.NonInvertibleMatrixError): + pass + + for eps in (1e-15, 1e-12, 1e-9, 1e-6): + try: + reg = m + sp.eye(m.rows) * eps + sol = reg.LUsolve(r) + warnings.warn( + f"Fast reduction regularized singular {context} with eps={eps:.1e}", + RuntimeWarning, + ) + x = np.array(sol, dtype=float) + return x if rhs.ndim == 2 else x[:, 0] + except (ValueError, sp.matrices.exceptions.NonInvertibleMatrixError): + continue + + raise ValueError(f"Failed to solve {context} with numeric LU.") + + def _reduce_symbolic(self): + """ + Eliminate algebraic variables using symbolic SymPy arithmetic. + + With differential and algebraic variables denoted by ``x_d`` and + ``x_a``, the algebraic MNA rows are + + ``G_alg_d*x_d + G_alg_a*x_a = B_alg*u``. + + The method computes + + ``Phi = -inv(G_alg_a)*G_alg_d`` and + ``Psi = inv(G_alg_a)*B_alg``, + + so that ``x_a = Phi*x_d + Psi*u``. Substitution into the differential + rows yields + + ``A = -inv(E_dd)*(G_diff_d + G_diff_a*Phi)`` and + ``B_ss = inv(E_dd)*(B_diff - G_diff_a*Psi)``. + + SymPy retains symbolic values through elimination and converts only the + final state matrices to ``float`` arrays. This makes the mode useful as + a correctness reference for small circuits, but explicit symbolic + inverses can become expensive in time and memory as the netlist grows. + Singular algebraic and state matrices are reported without numerical + regularization. + """ + E, G, B = self.E, self.G, self.B + dr, ar = self.diff_rows, self.alg_rows + dc, ac = self.diff_cols, self.alg_cols + + E_dd = E[dr, dc] # square, the "mass"/coupling matrix + G_alg_d = G[ar, dc] + G_alg_a = G[ar, ac] + G_diff_d = G[dr, dc] + G_diff_a = G[dr, ac] + B_alg = B[ar, :] + B_diff = B[dr, :] + + try: + G_alg_a_inv = G_alg_a.inv() + except sp.matrices.exceptions.NonInvertibleMatrixError: + self._raise_algebraic_singular() + Phi = -G_alg_a_inv * G_alg_d # x_a = Phi * x_d + Psi * u + Psi = G_alg_a_inv * B_alg + + try: + E_dd_inv = E_dd.inv() + except sp.matrices.exceptions.NonInvertibleMatrixError: + self._raise_state_singular() + A = -E_dd_inv * (G_diff_d + G_diff_a * Phi) + Bmat = E_dd_inv * (B_diff - G_diff_a * Psi) + + self.A = np.array(A.evalf(), dtype=float) + self.B_ss = np.array(Bmat.evalf(), dtype=float) + self.Phi = Phi + self.Psi = Psi + + def _reduce_fast(self): + """ + Eliminate algebraic variables using floating-point linear solves. + + This method implements the same block elimination and equations as + :meth:`_reduce_symbolic`, but avoids symbolic matrix inversion: + + 1. Extract ``E_dd`` and the differential/algebraic blocks of ``G`` and + ``B`` as ``float`` arrays. + 2. Factor ``G_alg_a`` with MUMPS and solve for ``Phi`` and ``Psi``. + 3. Form the Schur-complement terms involving ``G_diff_a``. + 4. Solve the differential mass system ``E_dd`` for ``A`` and ``B_ss``. + + The algebraic MUMPS factorization is reused across right-hand sides. + A small diagonal regularization ladder is available for numerically + singular systems and is always announced with a warning. The state + mass solve uses numeric SymPy LU in this implementation to avoid + environment-specific native dense-linear-algebra failures. + + ``"fast"`` changes only construction cost and numerical precision; it + does not simplify the circuit, discard states, or alter output + equations. It is generally preferred for large transmission-line and + densely coupled RLC/K models. + """ + E, G, B = self.E, self.G, self.B + dr, ar = self.diff_rows, self.alg_rows + dc, ac = self.diff_cols, self.alg_cols + + if isinstance(E, np.ndarray): + E_dd = E[np.ix_(dr, dc)] + G_alg_d = G[np.ix_(ar, dc)] + G_alg_a = G[np.ix_(ar, ac)] + G_diff_d = G[np.ix_(dr, dc)] + G_diff_a = G[np.ix_(dr, ac)] + B_alg = B[np.ix_(ar, range(self.n_u))] + B_diff = B[np.ix_(dr, range(self.n_u))] + else: + E_dd = np.array(E[dr, dc].evalf(), dtype=float) + G_alg_d = np.array(G[ar, dc].evalf(), dtype=float) + G_alg_a = np.array(G[ar, ac].evalf(), dtype=float) + G_diff_d = np.array(G[dr, dc].evalf(), dtype=float) + G_diff_a = np.array(G[dr, ac].evalf(), dtype=float) + B_alg = np.array(B[ar, :].evalf(), dtype=float) + B_diff = np.array(B[dr, :].evalf(), dtype=float) + + try: + Phi = -self._solve_with_mumps(G_alg_a, G_alg_d, "algebraic subsystem") + Psi = self._solve_with_mumps(G_alg_a, B_alg, "algebraic subsystem") + except (ImportError, ValueError): + self._raise_algebraic_singular() + + try: + gphi = self._safe_matmul(G_diff_a, Phi) + gpsi = self._safe_matmul(G_diff_a, Psi) + A = -self._solve_with_sympy_numeric(E_dd, G_diff_d + gphi, "state/mass subsystem") + Bmat = self._solve_with_sympy_numeric(E_dd, B_diff - gpsi, "state/mass subsystem") + except (ImportError, ValueError): + self._raise_state_singular() + + self.A = A + self.B_ss = Bmat + self.Phi = Phi + self.Psi = Psi + + def _build_output_rows( + self, + coefficient_rows: list[dict[str, float]], + ) -> tuple[np.ndarray, np.ndarray]: + """ + Project raw MNA output expressions to state-space rows. + + Each coefficient mapping uses keys of the form ``node:name``, + ``iL:name``, or ``iV:name``. Algebraic unknowns are substituted through + ``Phi`` and ``Psi`` so each result satisfies ``y = C*x + D*u``. + """ + if not coefficient_rows: + return ( + np.empty((0, len(self.state_labels)), dtype=float), + np.empty((0, self.n_u), dtype=float), + ) + + vec = np.zeros((len(coefficient_rows), self.n_total), dtype=float) + for row, coeffs in enumerate(coefficient_rows): + for key, val in coeffs.items(): + try: + typ, name = key.split(":", 1) + except ValueError as exc: + raise ValueError(f"Invalid output coefficient key '{key}'") from exc + + if typ == "node": + col = self._col_node(name) + elif typ == "iL": + col = self._col_iL(name) + elif typ == "iV": + col = self._col_iV(name) + else: + raise ValueError(f"Unknown output coefficient type '{typ}'") + if col is not None: + vec[row, col] += float(val) + + vec_d = vec[:, self.diff_cols] + vec_a = vec[:, self.alg_cols] + if self.reduction_mode == "fast": + C = vec_d + self._safe_matmul(vec_a, self.Phi) + D = self._safe_matmul(vec_a, self.Psi) + return C, D + + vec_sym = sp.Matrix(vec.tolist()) + vec_d_sym = vec_sym[:, self.diff_cols] + vec_a_sym = vec_sym[:, self.alg_cols] + C = vec_d_sym + vec_a_sym * self.Phi + D = vec_a_sym * self.Psi + return np.array(C.evalf(), dtype=float), np.array(D.evalf(), dtype=float) + + def _build_output_row( + self, + coefficients: dict[str, float], + ) -> tuple[np.ndarray, np.ndarray]: + """Project one raw MNA expression to one state-space output row.""" + C, D = self._build_output_rows([coefficients]) + return C[0], D[0] + + def _state_derivative_row( + self, + node: str, + ) -> tuple[np.ndarray, np.ndarray]: + """Return rows representing ``dV(node)/dt = C*x + D*u``.""" + n_states = len(self.state_labels) + if node in GND: + return np.zeros(n_states), np.zeros(self.n_u) + label = f"v_{node}" + if label not in self.state_labels: + raise ValueError( + f"node '{node}' has no capacitor attached, so its voltage isn't " + f"a state and dv/dt isn't defined this way. States: {self.state_labels}" + ) + idx = self.state_label_idx[label] + return self.A[idx, :], self.B_ss[idx, :] + + def add_node_voltage_output(self, node_name: str) -> None: + """ + Add a node-to-ground voltage to the selected system outputs. + + Parameters + ---------- + node_name: + Netlist node whose voltage is measured relative to ground. + + Raises + ------ + ValueError + If ``node_name`` is not present in the circuit. + """ + node = _normalize_node(node_name) + if node not in GND and node not in self.node_idx: + raise ValueError(f"Unknown node '{node_name}'") + + coefficients = {} if node in GND else {f"node:{node}": 1.0} + self._selected_output_rows.append(self._build_output_row(coefficients)) + self.output_labels.append(f"V({node})") + + def add_dipole_current_output(self, dipole_name: str) -> None: + """ + Add the current through a two-terminal netlist element. + + Parameters + ---------- + dipole_name: + Name of a resistor, inductor, capacitor, voltage source, or current + source. Positive current follows the element declaration from + ``n1`` to ``n2``. + + Raises + ------ + ValueError + If no supported dipole has the requested name. + + Notes + ----- + Resistor current is derived from Ohm's law. Inductor and voltage-source + currents are MNA branch unknowns. Capacitor current is computed from its + voltage derivative. Current-source current is its corresponding input + signal directly. + """ + try: + element = self.dipoles_by_name[dipole_name] + except KeyError as exc: + raise ValueError(f"Unknown dipole '{dipole_name}'") from exc + + if element.kind == "R": + coefficients = {} + if element.n1 not in GND: + coefficients[f"node:{element.n1}"] = 1.0 / element.value + if element.n2 not in GND: + key = f"node:{element.n2}" + coefficients[key] = coefficients.get(key, 0.0) - 1.0 / element.value + row = self._build_output_row(coefficients) + elif element.kind == "L": + row = self._build_output_row({f"iL:{element.name}": 1.0}) + elif element.kind == "C": + A_n1, B_n1 = self._state_derivative_row(element.n1) + A_n2, B_n2 = self._state_derivative_row(element.n2) + row = ( + element.value * (A_n1 - A_n2), + element.value * (B_n1 - B_n2), + ) + elif element.kind == "V": + row = self._build_output_row({f"iV:{element.name}": 1.0}) + elif element.kind == "I": + C_row = np.zeros(len(self.state_labels), dtype=float) + D_row = np.zeros(self.n_u, dtype=float) + D_row[self.input_names.index(element.name)] = 1.0 + row = C_row, D_row + else: + raise ValueError( + f"Element '{dipole_name}' of type '{element.kind}' is not a dipole" + ) + + self._selected_output_rows.append(row) + self.output_labels.append(f"I({element.name})") + + def get_system( + self, + ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """ + Return the reduced system with all selected outputs. + + Returns + ------- + A, B, C, D: + Matrices satisfying ``dx/dt = A*x + B*u`` and + ``y = C*x + D*u``. Rows of ``C`` and ``D`` follow the order in + which outputs were added. With no selected outputs, ``C`` and ``D`` + have zero rows and retain the correct state/input column counts. + """ + if not self._selected_output_rows: + return ( + self.A, + self.B_ss, + np.empty((0, len(self.state_labels)), dtype=float), + np.empty((0, self.n_u), dtype=float), + ) + C = np.vstack([row[0] for row in self._selected_output_rows]) + D = np.vstack([row[1] for row in self._selected_output_rows]) + return self.A, self.B_ss, C, D + + +class NetlistStateSpace(StateSpace): + """ + PathSim state-space block constructed directly from a linear netlist. + + Parameters + ---------- + netlist: + Existing netlist path or inline netlist text. A :class:`str` is treated + as a path when it names an existing file and as netlist text otherwise. + A :class:`Path` is always treated as an explicit file path. + output_voltages: + Node names exposed as node-to-ground voltage outputs. + output_currents: + Names of R, L, C, voltage-source, or current-source dipoles exposed as + current outputs. Positive current follows each netlist ``n1 -> n2`` + declaration. + reduction_mode: + Circuit reduction backend passed to :class:`CircuitModel`. + initial_value: + Initial differential state. Defaults to zero for every state. + + Notes + ----- + Output ports are ordered with all requested voltages first, followed by all + requested currents. Their labels are ``V(node)`` and ``I(dipole)``. + Input ports retain the voltage-source-then-current-source ordering of the + netlist model. The underlying :class:`CircuitModel` is available as + :attr:`model`. + + Examples + -------- + >>> block = NetlistStateSpace( + ... "filter.net", + ... output_voltages=["n1"], + ... output_currents=["Rload"], + ... reduction_mode="fast", + ... ) + """ + + def __init__( + self, + netlist: str | Path, + output_voltages: list[str] | None = None, + output_currents: list[str] | None = None, + reduction_mode: ReductionMode = "symbolic", + initial_value: np.ndarray | None = None, + ): + if isinstance(netlist, Path): + elements = parse_netlist_file(netlist) + elif isinstance(netlist, str): + candidate = Path(netlist) + try: + is_file = candidate.is_file() + except OSError: + is_file = False + elements = parse_netlist_file(candidate) if is_file else parse_netlist(netlist) + else: + raise TypeError("netlist must be a string or pathlib.Path") + + self.model = CircuitModel(elements, reduction_mode=reduction_mode) + for node_name in output_voltages or []: + self.model.add_node_voltage_output(node_name) + for dipole_name in output_currents or []: + self.model.add_dipole_current_output(dipole_name) + + A, B, C, D = self.model.get_system() + super().__init__( + A=A, + B=B, + C=C, + D=D, + initial_value=initial_value, + state_labels=list(self.model.state_labels), + input_labels=list(self.model.input_names), + output_labels=list(self.model.output_labels), + ) \ No newline at end of file From ab625b3e06e3772a5b6fe92f1a00f2405335374e Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Mon, 28 Sep 2026 19:50:50 +0200 Subject: [PATCH 03/18] added tests for new NetlistStateSpace cover the new class as well as associated subclass --- tests/filter.net | 16 ++++ tests/filter_without_load.net | 16 ++++ tests/test_netlist_parser.py | 155 +++++++++++++++++++++++++++++++ tests/test_netlist_statespace.py | 117 +++++++++++++++++++++++ 4 files changed, 304 insertions(+) create mode 100644 tests/filter.net create mode 100644 tests/filter_without_load.net create mode 100644 tests/test_netlist_parser.py create mode 100644 tests/test_netlist_statespace.py diff --git a/tests/filter.net b/tests/filter.net new file mode 100644 index 0000000..a7ce135 --- /dev/null +++ b/tests/filter.net @@ -0,0 +1,16 @@ +*validation netlist for pathsim +.param FREQ = 9k +.param DUTY = 0.5 +.param V_LOW = 0 +.param V_HIGH = 1 +.param T_PER = {1/FREQ} +.param T_ON = {T_PER * DUTY} +.param T_RISE = 10n +.param T_FALL = 10n + +V1 in 0 PULSE({V_LOW} {V_HIGH} 0 {T_RISE} {T_FALL} {T_ON} {T_PER}) +Lfilter in n1 0.5e-3 +Cfilter n1 0 50e-6 +Rload n1 0 50 +.tran 0.00001 0.2 0 0.00001 +.end diff --git a/tests/filter_without_load.net b/tests/filter_without_load.net new file mode 100644 index 0000000..19c47ac --- /dev/null +++ b/tests/filter_without_load.net @@ -0,0 +1,16 @@ +*validation netlist for pathsim +.param FREQ = 9k +.param DUTY = 0.5 +.param V_LOW = 0 +.param V_HIGH = 1 +.param T_PER = {1/FREQ} +.param T_ON = {T_PER * DUTY} +.param T_RISE = 10n +.param T_FALL = 10n + +V1 in 0 PULSE({V_LOW} {V_HIGH} 0 {T_RISE} {T_FALL} {T_ON} {T_PER}) +Bload n1 0 I=V(n1)/50 +Lfilter in n1 0.5e-3 +Cfilter n1 0 50e-6 +.tran 0.00001 0.2 0 0.00001 +.end diff --git a/tests/test_netlist_parser.py b/tests/test_netlist_parser.py new file mode 100644 index 0000000..6304b91 --- /dev/null +++ b/tests/test_netlist_parser.py @@ -0,0 +1,155 @@ +######################################################################################## +## +## TESTS FOR +## 'netlist_to_statespace.py' +## +######################################################################################## + +# IMPORTS ============================================================================== +import unittest +from pathlib import Path +from unittest.mock import patch +import numpy as np + +import pathsim_rf.netlist_to_statespace +from pathsim_rf.netlist_to_statespace import CircuitModel, parse_netlist, parse_value_with_units, parse_netlist_file + +# TESTS ================================================================================ + +class TestNetlistParser(unittest.TestCase): + """Test the netlist parser functions and CircuitModel used in NetlistStateSpace block.""" + + def test_parse_value_with_si_suffixes(self): + """Test conversion between spice si suffixes and power of ten.""" + cases = { + "10f": 10e-15, + "10p": 10e-12, + "10n": 10e-9, + "10u": 10e-6, + "10µ": 10e-6, + "10m": 10e-3, + "10k": 10e3, + "10K": 10e3, + "10meg": 10e6, + "10Meg": 10e6, + "10M": 10e6, + "10g": 10e9, + "10G": 10e9, + "10t": 10e12, + "10T": 10e12, + "100.e-3": 100e-3, + "2.2kOhm": 2200.0, + "3mH": 3e-3, + "4uF": 4e-6, + } + for raw, expected in cases.items(): + self.assertAlmostEqual(parse_value_with_units(raw), expected) + + def test_parse_param_reference(self): + """Test .param in spice list are taken into account""" + elements = parse_netlist( + """ + .param RMAIN = 12k + R1 n1 0 {RMAIN} + L1 n1 0 1mH + """ + ) + self.assertEqual(len(elements), 2) + self.assertAlmostEqual(elements[0].value, 12000.0) + self.assertAlmostEqual(elements[1].value, 1e-3) + + def test_parse_current_source_waveform_keeps_placeholder(self): + """Test if current sources with args are correctly parsed""" + elements = parse_netlist('I1 0 N3 PWL file="SIGNAL.TXT"') + self.assertEqual(len(elements), 1) + self.assertEqual(elements[0].kind, "I") + self.assertIsNone(elements[0].value) + + def test_parse_behavioral_source_placeholder(self): + """Test behavioral sources with are correctly parsed + and segragated between current and voltage sources""" + elements = parse_netlist("B1 0 N3 I=10*(exp(-t)-exp(-2*t)) \n B2 0 N4 V=10*(exp(-t)-exp(-2*t))") + self.assertEqual(len(elements), 2) + self.assertEqual(elements[0].kind, "I") + self.assertIsNone(elements[0].value) + self.assertEqual(elements[1].kind, "V") + self.assertIsNone(elements[0].value) + + def test_unsupported_element_raises(self): + """Test unsupported netlist elements raising an error""" + with self.assertRaises(ValueError): + parse_netlist("X1 n1 0 some_subckt") + + def test_get_system_without_outputs_has_zero_rows(self): + """Test system with zero output correctly formed""" + model = CircuitModel( + parse_netlist( + """ + V1 n1 0 1 + R1 n1 n2 1k + C1 n2 0 1u + """ + ) + ) + A, B, C, D = model.get_system() + self.assertEqual(C.shape, (0, A.shape[0])) + self.assertEqual(D.shape, (0, B.shape[1])) + self.assertEqual(model.output_labels, []) + + def test_fast_mode_falls_back_to_symbolic_without_mumps(self): + """Test class falling back to sympy if mumps not installed""" + elements = parse_netlist("V1 in 0 1\nR1 in out 10\nC1 out 0 1u") + with patch.object(pathsim_rf.netlist_to_statespace, "_MUMPS_AVAILABLE", False): + model = CircuitModel(elements, reduction_mode="fast") + + self.assertEqual(model.reduction_mode, "symbolic") + + def test_invalid_reduction_mode_raises(self): + """Test class raising an error for unknown reduction type""" + with self.assertRaises(ValueError): + CircuitModel(parse_netlist("R1 n1 0 1k"), reduction_mode="unknown") + + def test_outputs_accumulate_for_all_supported_dipoles(self): + """Test output currents can be selected as StateSpace output for all dipole types""" + model = CircuitModel( + parse_netlist( + """ + V1 n1 0 1 + R1 n1 n2 10 + C1 n2 0 1u + L1 n2 0 2m + I1 n2 0 1 + """ + ), + reduction_mode="symbolic", + ) + + model.add_node_voltage_output("n2") + for name in ["R1", "L1", "C1", "V1", "I1"]: + model.add_dipole_current_output(name) + A, B, C, D = model.get_system() + + self.assertEqual(C.shape, (6, A.shape[0])) + self.assertEqual(D.shape, (6, B.shape[1])) + self.assertEqual( + model.output_labels, + ["V(n2)", "I(R1)", "I(L1)", "I(C1)", "I(V1)", "I(I1)"], + ) + np.testing.assert_allclose(C[-1], 0.0) + np.testing.assert_allclose(D[-1], [0.0, 1.0]) + + def test_unknown_node_and_dipole_raise(self): + """Test error in case of unknown dipole or node selected as output""" + model = CircuitModel( + parse_netlist("V1 n1 0 1\nR1 n1 n2 1k\nC1 n2 0 1u"), + reduction_mode="symbolic", + ) + with self.assertRaisesRegex(ValueError, "Unknown node"): + model.add_node_voltage_output("missing") + with self.assertRaisesRegex(ValueError, "Unknown dipole"): + model.add_dipole_current_output("R404") + +# RUN TESTS LOCALLY ==================================================================== + +if __name__ == '__main__': + unittest.main(verbosity=2) \ No newline at end of file diff --git a/tests/test_netlist_statespace.py b/tests/test_netlist_statespace.py new file mode 100644 index 0000000..7a239f7 --- /dev/null +++ b/tests/test_netlist_statespace.py @@ -0,0 +1,117 @@ +######################################################################################## +## +## TESTS FOR +## 'netlist_to_statespace.py' +## +######################################################################################## + +# IMPORTS ============================================================================== +import unittest +from pathlib import Path +import numpy as np +import importlib.util + +from pathsim_rf.netlist_to_statespace import CircuitModel, NetlistStateSpace, parse_netlist_file, parse_netlist +from pathsim.blocks.lti import StateSpace + +# TESTS ================================================================================ + +class TestNetlistStateSpace(unittest.TestCase): + """Test the NetlistStateSpace block""" + + def test_path_string_matches_manual_circuit_model(self): + """Test CircuitModel is correctly used within NetlistStateSpace""" + netlist_path = "filter.net" + model = CircuitModel(parse_netlist_file(netlist_path), reduction_mode="symbolic") + model.add_node_voltage_output("n1") + model.add_dipole_current_output("Rload") + A, B, C, D = model.get_system() + + block = NetlistStateSpace( + str(netlist_path), + output_voltages=["n1"], + output_currents=["Rload"], + reduction_mode="symbolic", + ) + + self.assertIsInstance(block, StateSpace) + np.testing.assert_allclose(block.A, A) + np.testing.assert_allclose(block.B, B) + np.testing.assert_allclose(block.C, C) + np.testing.assert_allclose(block.D, D) + self.assertEqual(block.state_labels, model.state_labels) + self.assertEqual(block.input_labels, ["V1"]) + self.assertEqual(block.output_labels, ["V(n1)", "I(Rload)"]) + + def test_path_object_preserves_feedback_source_input(self): + """Test current sources feedback used to chain models""" + block = NetlistStateSpace("filter_without_load.net", + output_voltages=["n1"], + output_currents=["Lfilter"], + reduction_mode="symbolic", + ) + + self.assertEqual(block.input_labels, ["V1", "Bload"]) + self.assertEqual(block.output_labels, ["V(n1)", "I(Lfilter)"]) + self.assertEqual(len(block.inputs), 2) + self.assertEqual(len(block.outputs), 2) + + def test_inline_netlist_string(self): + """Test direct netlist use to class for instance generation""" + block = NetlistStateSpace( + """ + V1 in 0 1 + R1 in out 10 + C1 out 0 1u + """, + output_voltages=["out"], + output_currents=["R1", "V1"], + reduction_mode="symbolic", + ) + + self.assertEqual(block.input_labels, ["V1"]) + self.assertEqual(block.output_labels, ["V(out)", "I(R1)", "I(V1)"]) + self.assertEqual(block.C.shape[0], 3) + self.assertEqual(block.D.shape[0], 3) + + def test_explicit_missing_path_raises(self): + """Test missing .net file provided""" + with self.assertRaises(FileNotFoundError): + NetlistStateSpace(Path("missing.net")) + + def test_invalid_output_names_raise(self): + """Test invalid output names provided for state space generated""" + netlist_path = "filter.net" + with self.assertRaisesRegex(ValueError, "Unknown node"): + NetlistStateSpace(netlist_path, output_voltages=["missing"]) + with self.assertRaisesRegex(ValueError, "Unknown dipole"): + NetlistStateSpace(netlist_path, output_currents=["R404"]) + +@unittest.skipUnless(importlib.util.find_spec("mumps"), "python-mumps is not installed") +class TestFastReduction(unittest.TestCase): + """Class validation the fast implementation for reduction mode using mumps + using python-mumps""" + def test_fast_mode_matches_symbolic_on_small_case(self): + """Test if fast reduction mode is equivalent to symbolic + skipped if mumps not installed""" + elements = parse_netlist( + """ + V1 n1 0 1 + R1 n1 n2 2k + C1 n2 0 2u + L1 n2 n3 3m + L2 n3 0 5m + R2 n3 0 4k + K1 L1 L2 0.1 + """ + ) + symbolic = CircuitModel(elements, reduction_mode="symbolic") + fast = CircuitModel(elements, reduction_mode="fast") + + np.testing.assert_allclose(fast.A, symbolic.A, rtol=1e-9, atol=1e-12) + np.testing.assert_allclose(fast.B_ss, symbolic.B_ss, rtol=1e-9, atol=1e-12) + + +# RUN TESTS LOCALLY ==================================================================== +if __name__ == '__main__': + unittest.main(verbosity=2) \ No newline at end of file From 34d0e3a66a781229c0fa017ee3a9180dd72820c9 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Mon, 28 Sep 2026 19:52:46 +0200 Subject: [PATCH 04/18] Import NetlistStateSpace from netlist_to_statespace import only NetlistStateSpace as to keep netlist parser not accessible directly --- src/pathsim_rf/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/pathsim_rf/__init__.py b/src/pathsim_rf/__init__.py index 0531686..c3fe14b 100644 --- a/src/pathsim_rf/__init__.py +++ b/src/pathsim_rf/__init__.py @@ -15,6 +15,7 @@ from .transmission_line import * from .amplifier import * from .mixer import * +from .netlist_to_statespace import NetlistStateSpace try: from .network import * From 55f88888902fc6fa94210b54f8bd49d1d92705b4 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Mon, 28 Sep 2026 20:24:00 +0200 Subject: [PATCH 05/18] Add files via upload --- docs/source/examples/netlist_import.ipynb | 303 ++++++++++++++++++++++ 1 file changed, 303 insertions(+) create mode 100644 docs/source/examples/netlist_import.ipynb diff --git a/docs/source/examples/netlist_import.ipynb b/docs/source/examples/netlist_import.ipynb new file mode 100644 index 0000000..4e79ce6 --- /dev/null +++ b/docs/source/examples/netlist_import.ipynb @@ -0,0 +1,303 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Netlist State-Space Validation\n", + "\n", + "This example validates the `NetlistStateSpace` block against an equivalent PathSim implementation of an LC filter with a resistive load.\n", + "\n", + "Two circuit descriptions are used. The first netlist contains the load explicitly as a resistor. The second represents the same load using a behavioral current source. Both netlists are passed directly to `NetlistStateSpace` without requiring external `.net` files.\n", + "\n", + "The simulation compares the PathSim implementation and the imported netlist models using the recorded capacitor voltage and load current signals.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "from pathsim import Simulation, Connection\n", + "from pathsim.blocks import (\n", + " Adder,\n", + " Amplifier,\n", + " Integrator,\n", + " PulseSource,\n", + " Scope\n", + ")\n", + "from pathsim.solvers import EUB\n", + "from pathsim_rf.netlist_to_statespace import NetlistStateSpace\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Rload = 50\n", + "Cfilter = 50e-6\n", + "Lfilter = 0.5e-3\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Netlists\n", + "\n", + "The first netlist contains the 50 Ω load explicitly as `Rload`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "filter_netlist = \"\"\"\n", + "*validation netlist for pathsim\n", + ".param FREQ = 9k\n", + ".param DUTY = 0.5\n", + ".param V_LOW = 0\n", + ".param V_HIGH = 1\n", + ".param T_PER = {1/FREQ}\n", + ".param T_ON = {T_PER * DUTY}\n", + ".param T_RISE = 10n\n", + ".param T_FALL = 10n\n", + "\n", + "V1 in 0 PULSE({V_LOW} {V_HIGH} 0 {T_RISE} {T_FALL} {T_ON} {T_PER})\n", + "Lfilter in n1 0.5e-3\n", + "Cfilter n1 0 50e-6\n", + "Rload n1 0 50\n", + ".tran 0.00001 0.2 0 0.000001\n", + ".end\n", + "\"\"\"\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The second netlist represents the load using a behavioral current source.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "filter_without_load_netlist = \"\"\"\n", + "*validation netlist for pathsim\n", + ".param FREQ = 9k\n", + ".param DUTY = 0.5\n", + ".param V_LOW = 0\n", + ".param V_HIGH = 1\n", + ".param T_PER = {1/FREQ}\n", + ".param T_ON = {T_PER * DUTY}\n", + ".param T_RISE = 10n\n", + ".param T_FALL = 10n\n", + "\n", + "V1 in 0 PULSE({V_LOW} {V_HIGH} 0 {T_RISE} {T_FALL} {T_ON} {T_PER})\n", + "Bload n1 0 I=V(n1)/50\n", + "Lfilter in n1 0.5e-3\n", + "Cfilter n1 0 50e-6\n", + ".tran 0.00001 0.2 0 0.00001\n", + ".end\n", + "\"\"\"\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Blocks\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Sources\n", + "pulsesource = PulseSource(\n", + " T=1/9000,\n", + " amplitude=1,\n", + " duty=0.5,\n", + " t_rise=10e-9, \n", + " t_fall=10e-9\n", + ")\n", + "\n", + "state_space = NetlistStateSpace(\n", + " filter_netlist,\n", + " output_voltages=[\"n1\"],\n", + " output_currents=[\"Rload\"],\n", + " reduction_mode=\"fast\",\n", + ")\n", + "\n", + "state_space_bis = NetlistStateSpace(\n", + " filter_without_load_netlist,\n", + " output_voltages=[\"n1\"],\n", + " output_currents=[\"Lfilter\"],\n", + " reduction_mode=\"symbolic\",\n", + ")\n", + "\n", + "# Dynamic\n", + "integrator = Integrator(initial_value=0)\n", + "block_2 = Integrator()\n", + "\n", + "# Algebraic\n", + "adder = Adder(operations=\"+-\")\n", + "amplifier = Amplifier(gain=1/Lfilter)\n", + "block_5 = Adder(operations=\"+-\")\n", + "block_6 = Amplifier(gain=1/Cfilter)\n", + "block_7 = Amplifier(gain=1/Rload)\n", + "block_8 = Amplifier(gain=1/Rload)\n", + "\n", + "# Recording\n", + "scope = Scope(\n", + " sampling_period=1/100000,\n", + " labels=[\n", + " 'V(pwm)',\n", + " 'V(Cfull_pathsim)',\n", + " 'V(C_net_imported)',\n", + " 'V(Bload)'\n", + " ]\n", + ")\n", + "\n", + "blocks = [\n", + " pulsesource,\n", + " integrator,\n", + " block_2,\n", + " adder,\n", + " amplifier,\n", + " block_5,\n", + " block_6,\n", + " block_7,\n", + " scope,\n", + " state_space,\n", + " state_space_bis,\n", + " block_8\n", + "]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Connections\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "conn_0 = Connection(pulsesource[0], adder[0], state_space[0], state_space_bis[0])\n", + "conn_1 = Connection(adder[0], integrator[0])\n", + "conn_2 = Connection(integrator[0], amplifier[0])\n", + "conn_3 = Connection(amplifier[0], block_5[0])\n", + "conn_4 = Connection(block_5[0], block_2[0])\n", + "conn_5 = Connection(block_2[0], block_6[0])\n", + "conn_6 = Connection(block_6[0], block_7[0])\n", + "conn_7 = Connection(block_7[0], block_5[1])\n", + "conn_8 = Connection(block_6[0], adder[1])\n", + "conn_9 = Connection(pulsesource[0], scope[0])\n", + "conn_10 = Connection(block_6[0], scope[1])\n", + "conn_11 = Connection(state_space[0], scope[2])\n", + "conn_12 = Connection(state_space_bis[0], block_8[0], scope[3])\n", + "conn_13 = Connection(block_8[0], state_space_bis[1])\n", + "\n", + "connections = [\n", + " conn_0,\n", + " conn_1,\n", + " conn_2,\n", + " conn_3,\n", + " conn_4,\n", + " conn_5,\n", + " conn_6,\n", + " conn_7,\n", + " conn_8,\n", + " conn_9,\n", + " conn_10,\n", + " conn_11,\n", + " conn_12,\n", + " conn_13\n", + "]\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulation\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim = Simulation(\n", + " blocks,\n", + " connections,\n", + " Solver=EUB,\n", + " dt=0.5e-6,\n", + " dt_min=1e-16,\n", + " dt_max=1e-6,\n", + " tolerance_lte_rel=1e-6,\n", + " tolerance_lte_abs=1e-10,\n", + " tolerance_fpi=1e-10,\n", + ")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Results\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sim.run(duration=0.05)\n", + "\n", + "sim.plot()\n", + "plt.show()\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.x", + "mimetype": "text/x-python", + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "pygments_lexer": "ipython3", + "nbconvert_exporter": "python", + "file_extension": ".py" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From bc9c676bcd2cceb2582605d00dc84cbd4418fac9 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Tue, 29 Sep 2026 09:43:31 +0200 Subject: [PATCH 06/18] Add NetlistStateSpace entry to README --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index 5abdd88..638adbe 100644 --- a/README.md +++ b/README.md @@ -27,6 +27,7 @@ PathSim-RF extends the [PathSim](https://github.com/pathsim/pathsim) simulation | Block | Description | Key Parameters | |-------|-------------|----------------| | `RFNetwork` | N-port network from S-parameter data (Touchstone) via vector fitting | `ntwk`, `auto_fit` | +| `NetlistStateSpace` | Conversion of spice netlist into pathsim StateSpace | `netlist`, `output_voltages`, `output_currents` | | `TransmissionLine` | Lossy delay-based transmission line (scattering domain) | `length`, `er`, `attenuation`, `Z0` | | `RFAmplifier` | Amplifier with optional IP3 nonlinearity | `gain` [dB], `IIP3` [dBm], `P1dB` [dBm], `Z0` | | `RFMixer` | Ideal frequency converter (time-domain multiplication) | `conversion_gain` [dB], `Z0` | From 5660a8cc94f745a4ba337f9930e2ccdb0107d689 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Tue, 29 Sep 2026 10:52:30 +0200 Subject: [PATCH 07/18] Downgrade python-mumps dependency version --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 0ef14f9..24cc698 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,7 @@ dependencies = [ "scipy>=1.2", "scikit-rf>=0.30", "sympy>=1.14", - "python-mumps>=0.0.6; sys_platform != 'emscripten'", + "python-mumps>=0.0.4; sys_platform != 'emscripten'", # If a dependency ships native code that can't be installed in Pyodide, # tag it with `; sys_platform != 'emscripten'` so a normal pip install # still pulls it eagerly but micropip in the browser skips it. Pair From e06f7785a61c7c071fc683454e498f01b98a0bb4 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Tue, 29 Sep 2026 13:31:41 +0200 Subject: [PATCH 08/18] Update python-mumps dependency condition --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 24cc698..0ab1505 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -31,7 +31,7 @@ dependencies = [ "scipy>=1.2", "scikit-rf>=0.30", "sympy>=1.14", - "python-mumps>=0.0.4; sys_platform != 'emscripten'", + "python-mumps>=0.0.4; python_version >= '3.11' and sys_platform != 'emscripten'", # If a dependency ships native code that can't be installed in Pyodide, # tag it with `; sys_platform != 'emscripten'` so a normal pip install # still pulls it eagerly but micropip in the browser skips it. Pair From 6b5304ccf13a002c7bcbb5162352e868e93bf781 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Tue, 29 Sep 2026 14:40:22 +0200 Subject: [PATCH 09/18] Implemented and tested the revised "symbolic" to use scipy improve significantly the time for large netlist with k couplings --- src/pathsim_rf/netlist_to_statespace.py | 129 ++++++++++++++---------- 1 file changed, 73 insertions(+), 56 deletions(-) diff --git a/src/pathsim_rf/netlist_to_statespace.py b/src/pathsim_rf/netlist_to_statespace.py index 14fefaa..6c94e3d 100644 --- a/src/pathsim_rf/netlist_to_statespace.py +++ b/src/pathsim_rf/netlist_to_statespace.py @@ -279,12 +279,11 @@ class CircuitModel: Method used to eliminate the algebraic MNA variables: ``"symbolic"`` - Stamp and reduce with SymPy matrices. Matrix inverses and products - retain symbolic arithmetic until the final conversion to floating - point. This is the default and is useful for small circuits, - reference results, and diagnosing rank deficiencies. Its runtime - and memory use can grow quickly on large or densely coupled - circuits. + Stamp the MNA system with SymPy, convert the selected matrix blocks + to floating point, and reduce them with SciPy sparse LU solves. + This is the default and does not require ``python-mumps``. Keeping + symbolic stamping preserves the reference assembly path while + avoiding expensive symbolic matrix inverses on large circuits. ``"fast"`` Stamp directly into floating-point arrays. Eliminate the algebraic @@ -685,61 +684,87 @@ def _solve_with_sympy_numeric(self, mat: np.ndarray, rhs: np.ndarray, context: s raise ValueError(f"Failed to solve {context} with numeric LU.") - def _reduce_symbolic(self): - """ - Eliminate algebraic variables using symbolic SymPy arithmetic. - - With differential and algebraic variables denoted by ``x_d`` and - ``x_a``, the algebraic MNA rows are - - ``G_alg_d*x_d + G_alg_a*x_a = B_alg*u``. - - The method computes + @staticmethod + def _solve_with_scipy(mat: np.ndarray, rhs: np.ndarray, context: str) -> np.ndarray: + """Solve right-hand sides with SciPy sparse LU and bounded regularization.""" + import scipy.sparse as sps + from scipy.sparse.linalg import splu - ``Phi = -inv(G_alg_a)*G_alg_d`` and - ``Psi = inv(G_alg_a)*B_alg``, + rhs_2d = rhs if rhs.ndim == 2 else rhs.reshape((-1, 1)) + if mat.shape[0] != mat.shape[1]: + raise ValueError(f"{context} matrix must be square, got {mat.shape}") + if rhs_2d.shape[0] != mat.shape[0]: + raise ValueError( + f"{context} rhs has incompatible shape {rhs_2d.shape} for matrix {mat.shape}" + ) - so that ``x_a = Phi*x_d + Psi*u``. Substitution into the differential - rows yields + last_error: Exception | None = None + for eps in (0.0, 1e-15, 1e-12, 1e-9, 1e-6): + try: + matrix = mat if eps == 0.0 else mat + eps * np.eye(mat.shape[0]) + factor = splu(sps.csc_matrix(matrix)) + columns = [ + np.asarray(factor.solve(np.asarray(column, dtype=float)), dtype=float) + for column in rhs_2d.T + ] + solution = np.column_stack(columns) + if eps > 0.0: + warnings.warn( + f"Symbolic reduction regularized singular {context} with eps={eps:.1e}", + RuntimeWarning, + ) + return solution if rhs.ndim == 2 else solution[:, 0] + except (RuntimeError, ValueError) as exc: + last_error = exc - ``A = -inv(E_dd)*(G_diff_d + G_diff_a*Phi)`` and - ``B_ss = inv(E_dd)*(B_diff - G_diff_a*Psi)``. + raise ValueError(f"Failed to solve {context} with SciPy sparse LU.") from last_error - SymPy retains symbolic values through elimination and converts only the - final state matrices to ``float`` arrays. This makes the mode useful as - a correctness reference for small circuits, but explicit symbolic - inverses can become expensive in time and memory as the netlist grows. - Singular algebraic and state matrices are reported without numerical + def _reduce_symbolic(self): + """ + Stamp with SymPy and eliminate algebraic variables using SciPy. + + The symbolic MNA blocks are converted to floating-point arrays before + factorization. SciPy sparse LU then solves for ``Phi`` and ``Psi`` and + the differential mass system without forming explicit inverses. This + retains the default SymPy stamping implementation while making the + reduction practical for substantially larger circuits. Unlike fast + mode, this path requires no ``python-mumps`` and applies no diagonal regularization. """ E, G, B = self.E, self.G, self.B dr, ar = self.diff_rows, self.alg_rows dc, ac = self.diff_cols, self.alg_cols - E_dd = E[dr, dc] # square, the "mass"/coupling matrix - G_alg_d = G[ar, dc] - G_alg_a = G[ar, ac] - G_diff_d = G[dr, dc] - G_diff_a = G[dr, ac] - B_alg = B[ar, :] - B_diff = B[dr, :] + E_dd = np.array(E[dr, dc], dtype=float) + G_alg_d = np.array(G[ar, dc], dtype=float) + G_alg_a = np.array(G[ar, ac], dtype=float) + G_diff_d = np.array(G[dr, dc], dtype=float) + G_diff_a = np.array(G[dr, ac], dtype=float) + B_alg = np.array(B[ar, :], dtype=float) + B_diff = np.array(B[dr, :], dtype=float) try: - G_alg_a_inv = G_alg_a.inv() - except sp.matrices.exceptions.NonInvertibleMatrixError: + Phi = -self._solve_with_scipy(G_alg_a, G_alg_d, "algebraic subsystem") + Psi = self._solve_with_scipy(G_alg_a, B_alg, "algebraic subsystem") + except ValueError: self._raise_algebraic_singular() - Phi = -G_alg_a_inv * G_alg_d # x_a = Phi * x_d + Psi * u - Psi = G_alg_a_inv * B_alg try: - E_dd_inv = E_dd.inv() - except sp.matrices.exceptions.NonInvertibleMatrixError: + A = -self._solve_with_scipy( + E_dd, + G_diff_d + self._safe_matmul(G_diff_a, Phi), + "state/mass subsystem", + ) + Bmat = self._solve_with_scipy( + E_dd, + B_diff - self._safe_matmul(G_diff_a, Psi), + "state/mass subsystem", + ) + except ValueError: self._raise_state_singular() - A = -E_dd_inv * (G_diff_d + G_diff_a * Phi) - Bmat = E_dd_inv * (B_diff - G_diff_a * Psi) - self.A = np.array(A.evalf(), dtype=float) - self.B_ss = np.array(Bmat.evalf(), dtype=float) + self.A = A + self.B_ss = Bmat self.Phi = Phi self.Psi = Psi @@ -845,17 +870,9 @@ def _build_output_rows( vec_d = vec[:, self.diff_cols] vec_a = vec[:, self.alg_cols] - if self.reduction_mode == "fast": - C = vec_d + self._safe_matmul(vec_a, self.Phi) - D = self._safe_matmul(vec_a, self.Psi) - return C, D - - vec_sym = sp.Matrix(vec.tolist()) - vec_d_sym = vec_sym[:, self.diff_cols] - vec_a_sym = vec_sym[:, self.alg_cols] - C = vec_d_sym + vec_a_sym * self.Phi - D = vec_a_sym * self.Psi - return np.array(C.evalf(), dtype=float), np.array(D.evalf(), dtype=float) + C = vec_d + self._safe_matmul(vec_a, self.Phi) + D = self._safe_matmul(vec_a, self.Psi) + return C, D def _build_output_row( self, @@ -1065,4 +1082,4 @@ def __init__( state_labels=list(self.model.state_labels), input_labels=list(self.model.input_names), output_labels=list(self.model.output_labels), - ) \ No newline at end of file + ) From f6c2c4b0327907cd5621598c7f5e5c4ad8204575 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Wed, 30 Sep 2026 10:47:40 +0200 Subject: [PATCH 10/18] passed to full mumps solve in case of "fast" option selected --- src/pathsim_rf/netlist_to_statespace.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pathsim_rf/netlist_to_statespace.py b/src/pathsim_rf/netlist_to_statespace.py index 6c94e3d..9e6da8c 100644 --- a/src/pathsim_rf/netlist_to_statespace.py +++ b/src/pathsim_rf/netlist_to_statespace.py @@ -822,8 +822,8 @@ def _reduce_fast(self): try: gphi = self._safe_matmul(G_diff_a, Phi) gpsi = self._safe_matmul(G_diff_a, Psi) - A = -self._solve_with_sympy_numeric(E_dd, G_diff_d + gphi, "state/mass subsystem") - Bmat = self._solve_with_sympy_numeric(E_dd, B_diff - gpsi, "state/mass subsystem") + A = -self._solve_with_mumps(E_dd, G_diff_d + gphi, "state/mass subsystem") + Bmat = self._solve_with_mumps(E_dd, B_diff - gpsi, "state/mass subsystem") except (ImportError, ValueError): self._raise_state_singular() From 1ad3ae19316e97ccb679726cbec8a7d16f71367b Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Wed, 30 Sep 2026 15:06:15 +0200 Subject: [PATCH 11/18] Remove sympy and python-mumps dependencies Removed sympy and python-mumps dependencies for compatibility. --- pyproject.toml | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 0ab1505..104e1ff 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,8 +30,6 @@ dependencies = [ "numpy>=1.15", "scipy>=1.2", "scikit-rf>=0.30", - "sympy>=1.14", - "python-mumps>=0.0.4; python_version >= '3.11' and sys_platform != 'emscripten'", # If a dependency ships native code that can't be installed in Pyodide, # tag it with `; sys_platform != 'emscripten'` so a normal pip install # still pulls it eagerly but micropip in the browser skips it. Pair From 988e1b0ccfeaba973a82e845fd3a3888ceeab6fc Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Wed, 30 Sep 2026 15:08:04 +0200 Subject: [PATCH 12/18] Removed mumps and sympy dependancies --- src/pathsim_rf/netlist_to_statespace.py | 298 ++---------------------- 1 file changed, 20 insertions(+), 278 deletions(-) diff --git a/src/pathsim_rf/netlist_to_statespace.py b/src/pathsim_rf/netlist_to_statespace.py index 9e6da8c..5667b4a 100644 --- a/src/pathsim_rf/netlist_to_statespace.py +++ b/src/pathsim_rf/netlist_to_statespace.py @@ -19,14 +19,10 @@ """ from __future__ import annotations -import sympy as sp import numpy as np from dataclasses import dataclass from pathlib import Path -import importlib.util -import logging import re -from typing import Literal import warnings from pathsim.blocks.lti import StateSpace @@ -73,9 +69,6 @@ class Element: ) _VALID_KINDS = {"R", "L", "C", "V", "I", "K"} _IGNORED_LINE_STARTS = (".", "*", '"') -ReductionMode = Literal["symbolic", "fast"] -_MUMPS_AVAILABLE = importlib.util.find_spec("mumps") is not None -_LOGGER = logging.getLogger("pathsim.netlist_to_statespace") def _normalize_node(node: str) -> str: @@ -247,7 +240,7 @@ def parse_netlist_file(path: str | Path) -> list[Element]: # -------------------------------------------------------------------------- -# MNA + symbolic DAE -> explicit state space reduction +# MNA DAE -> explicit state space reduction # -------------------------------------------------------------------------- class CircuitModel: @@ -268,34 +261,15 @@ class CircuitModel: ``dx_d/dt = A * x_d + B_ss * u``. - Both reduction modes produce the same matrices and state ordering. They - differ only in the arithmetic and linear solver used during elimination. + The MNA matrices are stamped directly into floating-point arrays and + reduced using SciPy sparse LU solves. Singular solves are retried with + bounded diagonal regularization and emit :class:`RuntimeWarning` when + regularization is used. Parameters ---------- elements: Output of :func:`parse_netlist` / :func:`parse_netlist_file`. - reduction_mode: - Method used to eliminate the algebraic MNA variables: - - ``"symbolic"`` - Stamp the MNA system with SymPy, convert the selected matrix blocks - to floating point, and reduce them with SciPy sparse LU solves. - This is the default and does not require ``python-mumps``. Keeping - symbolic stamping preserves the reference assembly path while - avoiding expensive symbolic matrix inverses on large circuits. - - ``"fast"`` - Stamp directly into floating-point arrays. Eliminate the algebraic - subsystem with a sparse MUMPS factorization, then solve the - differential mass system numerically. This mode is intended for - large RLC/K netlists and repeated output construction. A singular - algebraic solve is retried with progressively larger diagonal - regularization and emits :class:`RuntimeWarning` when regularization - is used. When ``python-mumps`` is unavailable, construction falls - back to symbolic reduction and logs a warning through PathSim. - - The mode affects model construction only Attributes ---------- @@ -311,23 +285,8 @@ class CircuitModel: Labels matching the row and column ordering of the state matrices. """ - def __init__( - self, - elements: list[Element], - reduction_mode: ReductionMode = "symbolic", - ): - if reduction_mode not in {"symbolic", "fast"}: - raise ValueError( - f"Unsupported reduction_mode '{reduction_mode}'. " - "Expected 'symbolic' or 'fast'." - ) - if reduction_mode == "fast" and not _MUMPS_AVAILABLE: - _LOGGER.warning( - "Fast netlist reduction requires python-mumps; " - "falling back to symbolic reduction." - ) - reduction_mode = "symbolic" - self.reduction_mode = reduction_mode + def __init__(self, elements: list[Element]): + self.elements = elements self.R = [e for e in elements if e.kind == 'R'] @@ -373,26 +332,18 @@ def col_iV(name): return self.n_nodes + self.n_L + self.V_idx[name] self._col_node, self._col_iL, self._col_iV = col_node, col_iL, col_iV - numeric = self.reduction_mode == "fast" # ---- inductance matrix (with mutual terms) ---- - Lmat = ( - np.zeros((self.n_L, self.n_L), dtype=float) - if numeric - else sp.zeros(self.n_L, self.n_L) - ) + Lmat = np.zeros((self.n_L, self.n_L), dtype=float) l_values = {e.name: e.value for e in self.L} for e in self.L: i = self.L_idx[e.name] - Lmat[i, i] = float(e.value) if numeric else sp.nsimplify(e.value) + Lmat[i, i] = float(e.value) for k in self.K: i, j = self.L_idx[k.l1], self.L_idx[k.l2] Li = l_values[k.l1] Lj = l_values[k.l2] - if numeric: - M = float(k.value) * np.sqrt(float(Li) * float(Lj)) - else: - M = sp.nsimplify(k.value) * sp.sqrt(sp.nsimplify(Li) * sp.nsimplify(Lj)) + M = float(k.value) * np.sqrt(float(Li) * float(Lj)) Lmat[i, j] += M Lmat[j, i] += M self.Lmat = Lmat @@ -403,14 +354,9 @@ def col_iV(name): # ---- build E (dynamic) and G (algebraic) matrices, and B (input map) ---- n = self.n_total - if numeric: - E = np.zeros((n, n), dtype=float) - G = np.zeros((n, n), dtype=float) - B = np.zeros((n, self.n_u), dtype=float) - else: - E = sp.zeros(n, n) - G = sp.zeros(n, n) - B = sp.zeros(n, self.n_u) + E = np.zeros((n, n), dtype=float) + G = np.zeros((n, n), dtype=float) + B = np.zeros((n, self.n_u), dtype=float) def stamp_G(row, col, val): if row is not None and col is not None: @@ -419,14 +365,14 @@ def stamp_G(row, col, val): # Resistors: contribute to node KCL rows only for e in self.R: a, b = col_node(e.n1), col_node(e.n2) - g = (1.0 / float(e.value)) if numeric else (1 / sp.nsimplify(e.value)) + g = 1.0 / float(e.value) stamp_G(a, a, g); stamp_G(b, b, g) stamp_G(a, b, -g); stamp_G(b, a, -g) # Capacitors: contribute dv/dt terms to node KCL rows (the E matrix) for e in self.C: a, b = col_node(e.n1), col_node(e.n2) - c = float(e.value) if numeric else sp.nsimplify(e.value) + c = float(e.value) if a is not None: E[a, a] += c if b is not None: E[b, b] += c if a is not None and b is not None: @@ -493,25 +439,11 @@ def stamp_G(row, col, val): ) self.diff_cols, self.alg_cols = diff_cols, alg_cols + self._set_state_labels() self._reduce() self._selected_output_rows: list[tuple[np.ndarray, np.ndarray]] = [] self.output_labels: list[str] = [] - def _reduce(self): - """ - Reduce the stamped MNA differential-algebraic system. - - This dispatcher selects the arithmetic backend requested by - :attr:`reduction_mode`. Both implementations populate ``A``, ``B_ss``, - ``Phi``, and ``Psi`` and preserve the same differential-state ordering. - State labels are assigned only after a successful reduction. - """ - if self.reduction_mode == "fast": - self._reduce_fast() - else: - self._reduce_symbolic() - self._set_state_labels() - def _set_state_labels(self): """Set labels for the differential state vector in `self.diff_cols` order.""" dc = self.diff_cols @@ -556,90 +488,6 @@ def _raise_state_singular(self): "yourself before building the netlist." ) - def _solve_with_mumps(self, mat: np.ndarray, rhs: np.ndarray, context: str) -> np.ndarray: - """ - Solve one or more right-hand sides with one MUMPS factorization. - - Parameters - ---------- - mat: - Square coefficient matrix. - rhs: - Vector or matrix of right-hand sides. Each matrix column is solved - independently while reusing the factorization of ``mat``. - context: - Subsystem name included in errors and regularization warnings. - - Returns - ------- - numpy.ndarray - Solution with the same one- or two-dimensional convention as - ``rhs``. - - Notes - ----- - The unmodified matrix is attempted first. If factorization or solution - fails, the method retries ``mat + eps * I`` for increasing ``eps``. - Regularization permits reduction of numerically singular algebraic - systems, but it perturbs their constraints; every successful - regularized solve therefore emits :class:`RuntimeWarning`. - - MUMPS deliberately remains in its default unsymmetric mode. General - MNA systems may contain voltage-source saddle-point blocks, so they - are not positive definite, and the reduced transient matrix is not - generally symmetric. Selecting MUMPS ``sym=1`` (SPD) or passing - ``symmetric=True`` would therefore impose an invalid matrix property. - """ - try: - import mumps - import scipy.sparse as sps - except ImportError as exc: - raise ImportError( - "Fast reduction mode requires `python-mumps` (module `mumps`) " - "and scipy.sparse to be available." - ) from exc - - rhs_2d = rhs if rhs.ndim == 2 else rhs.reshape((-1, 1)) - if mat.shape[0] != mat.shape[1]: - raise ValueError(f"{context} matrix must be square, got {mat.shape}") - if rhs_2d.shape[0] != mat.shape[0]: - raise ValueError( - f"{context} rhs has incompatible shape {rhs_2d.shape} for matrix {mat.shape}" - ) - - last_error: Exception | None = None - for eps in (0.0, 1e-15, 1e-12, 1e-9, 1e-6): - try: - if eps > 0.0: - mat_reg = mat + np.eye(mat.shape[0], dtype=float) * eps - else: - mat_reg = mat - - ctx = mumps.Context() - ctx.set_matrix(sps.csc_matrix(mat_reg)) - ctx.factor() - - cols = [] - for j in range(rhs_2d.shape[1]): - bj = np.array(rhs_2d[:, j], dtype=float, order="F") - xj = np.array(ctx.solve(bj), dtype=float) - cols.append(xj) - x = np.column_stack(cols) - - if eps > 0.0: - warnings.warn( - f"Fast reduction regularized singular {context} with eps={eps:.1e}", - RuntimeWarning, - ) - return x if rhs.ndim == 2 else x[:, 0] - except (mumps.MUMPSError, ValueError, RuntimeError, TypeError) as exc: - last_error = exc - continue - - raise ValueError( - f"Failed to solve {context} even after regularization attempts." - ) from last_error - @staticmethod def _safe_matmul(a: np.ndarray, b: np.ndarray) -> np.ndarray: """ @@ -656,34 +504,6 @@ def _safe_matmul(a: np.ndarray, b: np.ndarray) -> np.ndarray: out[i, :] += aik * b[k, :] return out - def _solve_with_sympy_numeric(self, mat: np.ndarray, rhs: np.ndarray, context: str) -> np.ndarray: - """Solve linear system using numeric SymPy LU; supports regularization ladder.""" - rhs_2d = rhs if rhs.ndim == 2 else rhs.reshape((-1, 1)) - m = sp.Matrix(mat) - r = sp.Matrix(rhs_2d) - - try: - sol = m.LUsolve(r) - x = np.array(sol, dtype=float) - return x if rhs.ndim == 2 else x[:, 0] - except (ValueError, sp.matrices.exceptions.NonInvertibleMatrixError): - pass - - for eps in (1e-15, 1e-12, 1e-9, 1e-6): - try: - reg = m + sp.eye(m.rows) * eps - sol = reg.LUsolve(r) - warnings.warn( - f"Fast reduction regularized singular {context} with eps={eps:.1e}", - RuntimeWarning, - ) - x = np.array(sol, dtype=float) - return x if rhs.ndim == 2 else x[:, 0] - except (ValueError, sp.matrices.exceptions.NonInvertibleMatrixError): - continue - - raise ValueError(f"Failed to solve {context} with numeric LU.") - @staticmethod def _solve_with_scipy(mat: np.ndarray, rhs: np.ndarray, context: str) -> np.ndarray: """Solve right-hand sides with SciPy sparse LU and bounded regularization.""" @@ -710,7 +530,7 @@ def _solve_with_scipy(mat: np.ndarray, rhs: np.ndarray, context: str) -> np.ndar solution = np.column_stack(columns) if eps > 0.0: warnings.warn( - f"Symbolic reduction regularized singular {context} with eps={eps:.1e}", + f"Circuit reduction regularized singular {context} with eps={eps:.1e}", RuntimeWarning, ) return solution if rhs.ndim == 2 else solution[:, 0] @@ -719,18 +539,8 @@ def _solve_with_scipy(mat: np.ndarray, rhs: np.ndarray, context: str) -> np.ndar raise ValueError(f"Failed to solve {context} with SciPy sparse LU.") from last_error - def _reduce_symbolic(self): - """ - Stamp with SymPy and eliminate algebraic variables using SciPy. - - The symbolic MNA blocks are converted to floating-point arrays before - factorization. SciPy sparse LU then solves for ``Phi`` and ``Psi`` and - the differential mass system without forming explicit inverses. This - retains the default SymPy stamping implementation while making the - reduction practical for substantially larger circuits. Unlike fast - mode, this path requires no ``python-mumps`` and applies no diagonal - regularization. - """ + def _reduce(self): + """Eliminate algebraic variables and solve the mass system with SciPy.""" E, G, B = self.E, self.G, self.B dr, ar = self.diff_rows, self.alg_rows dc, ac = self.diff_cols, self.alg_cols @@ -768,70 +578,6 @@ def _reduce_symbolic(self): self.Phi = Phi self.Psi = Psi - def _reduce_fast(self): - """ - Eliminate algebraic variables using floating-point linear solves. - - This method implements the same block elimination and equations as - :meth:`_reduce_symbolic`, but avoids symbolic matrix inversion: - - 1. Extract ``E_dd`` and the differential/algebraic blocks of ``G`` and - ``B`` as ``float`` arrays. - 2. Factor ``G_alg_a`` with MUMPS and solve for ``Phi`` and ``Psi``. - 3. Form the Schur-complement terms involving ``G_diff_a``. - 4. Solve the differential mass system ``E_dd`` for ``A`` and ``B_ss``. - - The algebraic MUMPS factorization is reused across right-hand sides. - A small diagonal regularization ladder is available for numerically - singular systems and is always announced with a warning. The state - mass solve uses numeric SymPy LU in this implementation to avoid - environment-specific native dense-linear-algebra failures. - - ``"fast"`` changes only construction cost and numerical precision; it - does not simplify the circuit, discard states, or alter output - equations. It is generally preferred for large transmission-line and - densely coupled RLC/K models. - """ - E, G, B = self.E, self.G, self.B - dr, ar = self.diff_rows, self.alg_rows - dc, ac = self.diff_cols, self.alg_cols - - if isinstance(E, np.ndarray): - E_dd = E[np.ix_(dr, dc)] - G_alg_d = G[np.ix_(ar, dc)] - G_alg_a = G[np.ix_(ar, ac)] - G_diff_d = G[np.ix_(dr, dc)] - G_diff_a = G[np.ix_(dr, ac)] - B_alg = B[np.ix_(ar, range(self.n_u))] - B_diff = B[np.ix_(dr, range(self.n_u))] - else: - E_dd = np.array(E[dr, dc].evalf(), dtype=float) - G_alg_d = np.array(G[ar, dc].evalf(), dtype=float) - G_alg_a = np.array(G[ar, ac].evalf(), dtype=float) - G_diff_d = np.array(G[dr, dc].evalf(), dtype=float) - G_diff_a = np.array(G[dr, ac].evalf(), dtype=float) - B_alg = np.array(B[ar, :].evalf(), dtype=float) - B_diff = np.array(B[dr, :].evalf(), dtype=float) - - try: - Phi = -self._solve_with_mumps(G_alg_a, G_alg_d, "algebraic subsystem") - Psi = self._solve_with_mumps(G_alg_a, B_alg, "algebraic subsystem") - except (ImportError, ValueError): - self._raise_algebraic_singular() - - try: - gphi = self._safe_matmul(G_diff_a, Phi) - gpsi = self._safe_matmul(G_diff_a, Psi) - A = -self._solve_with_mumps(E_dd, G_diff_d + gphi, "state/mass subsystem") - Bmat = self._solve_with_mumps(E_dd, B_diff - gpsi, "state/mass subsystem") - except (ImportError, ValueError): - self._raise_state_singular() - - self.A = A - self.B_ss = Bmat - self.Phi = Phi - self.Psi = Psi - def _build_output_rows( self, coefficient_rows: list[dict[str, float]], @@ -1023,8 +769,6 @@ class NetlistStateSpace(StateSpace): Names of R, L, C, voltage-source, or current-source dipoles exposed as current outputs. Positive current follows each netlist ``n1 -> n2`` declaration. - reduction_mode: - Circuit reduction backend passed to :class:`CircuitModel`. initial_value: Initial differential state. Defaults to zero for every state. @@ -1042,7 +786,6 @@ class NetlistStateSpace(StateSpace): ... "filter.net", ... output_voltages=["n1"], ... output_currents=["Rload"], - ... reduction_mode="fast", ... ) """ @@ -1051,7 +794,6 @@ def __init__( netlist: str | Path, output_voltages: list[str] | None = None, output_currents: list[str] | None = None, - reduction_mode: ReductionMode = "symbolic", initial_value: np.ndarray | None = None, ): if isinstance(netlist, Path): @@ -1066,7 +808,7 @@ def __init__( else: raise TypeError("netlist must be a string or pathlib.Path") - self.model = CircuitModel(elements, reduction_mode=reduction_mode) + self.model = CircuitModel(elements) for node_name in output_voltages or []: self.model.add_node_voltage_output(node_name) for dipole_name in output_currents or []: From aa711d7245ef6547882ed28db713aeec46de03eb Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Wed, 30 Sep 2026 15:20:48 +0200 Subject: [PATCH 13/18] Refactor and add tests for netlist parser functionality --- tests/test_netlist_parser.py | 26 +++++++------------------- 1 file changed, 7 insertions(+), 19 deletions(-) diff --git a/tests/test_netlist_parser.py b/tests/test_netlist_parser.py index 6304b91..4faa14c 100644 --- a/tests/test_netlist_parser.py +++ b/tests/test_netlist_parser.py @@ -96,19 +96,6 @@ def test_get_system_without_outputs_has_zero_rows(self): self.assertEqual(D.shape, (0, B.shape[1])) self.assertEqual(model.output_labels, []) - def test_fast_mode_falls_back_to_symbolic_without_mumps(self): - """Test class falling back to sympy if mumps not installed""" - elements = parse_netlist("V1 in 0 1\nR1 in out 10\nC1 out 0 1u") - with patch.object(pathsim_rf.netlist_to_statespace, "_MUMPS_AVAILABLE", False): - model = CircuitModel(elements, reduction_mode="fast") - - self.assertEqual(model.reduction_mode, "symbolic") - - def test_invalid_reduction_mode_raises(self): - """Test class raising an error for unknown reduction type""" - with self.assertRaises(ValueError): - CircuitModel(parse_netlist("R1 n1 0 1k"), reduction_mode="unknown") - def test_outputs_accumulate_for_all_supported_dipoles(self): """Test output currents can be selected as StateSpace output for all dipole types""" model = CircuitModel( @@ -121,7 +108,6 @@ def test_outputs_accumulate_for_all_supported_dipoles(self): I1 n2 0 1 """ ), - reduction_mode="symbolic", ) model.add_node_voltage_output("n2") @@ -138,12 +124,14 @@ def test_outputs_accumulate_for_all_supported_dipoles(self): np.testing.assert_allclose(C[-1], 0.0) np.testing.assert_allclose(D[-1], [0.0, 1.0]) + def test_model_without_inputs_raises(self): + """Test error in case of no input present in netlist (in the form of current or voltage source)""" + with self.assertRaisesRegex(ValueError, "at least one independent"): + CircuitModel(parse_netlist("R1 n1 0 1k\nC1 n1 0 1u")) + def test_unknown_node_and_dipole_raise(self): """Test error in case of unknown dipole or node selected as output""" - model = CircuitModel( - parse_netlist("V1 n1 0 1\nR1 n1 n2 1k\nC1 n2 0 1u"), - reduction_mode="symbolic", - ) + model = CircuitModel(parse_netlist("V1 n1 0 1\nR1 n1 n2 1k\nC1 n2 0 1u")) with self.assertRaisesRegex(ValueError, "Unknown node"): model.add_node_voltage_output("missing") with self.assertRaisesRegex(ValueError, "Unknown dipole"): @@ -152,4 +140,4 @@ def test_unknown_node_and_dipole_raise(self): # RUN TESTS LOCALLY ==================================================================== if __name__ == '__main__': - unittest.main(verbosity=2) \ No newline at end of file + unittest.main(verbosity=2) From f7a90497c4c02d69973fc984eca3d6c9c3193483 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Wed, 30 Sep 2026 15:25:01 +0200 Subject: [PATCH 14/18] Refactor test cases for NetlistStateSpace --- tests/test_netlist_statespace.py | 40 ++++---------------------------- 1 file changed, 5 insertions(+), 35 deletions(-) diff --git a/tests/test_netlist_statespace.py b/tests/test_netlist_statespace.py index 7a239f7..3d6463b 100644 --- a/tests/test_netlist_statespace.py +++ b/tests/test_netlist_statespace.py @@ -22,7 +22,7 @@ class TestNetlistStateSpace(unittest.TestCase): def test_path_string_matches_manual_circuit_model(self): """Test CircuitModel is correctly used within NetlistStateSpace""" netlist_path = "filter.net" - model = CircuitModel(parse_netlist_file(netlist_path), reduction_mode="symbolic") + model = CircuitModel(parse_netlist_file(netlist_path)) model.add_node_voltage_output("n1") model.add_dipole_current_output("Rload") A, B, C, D = model.get_system() @@ -30,9 +30,7 @@ def test_path_string_matches_manual_circuit_model(self): block = NetlistStateSpace( str(netlist_path), output_voltages=["n1"], - output_currents=["Rload"], - reduction_mode="symbolic", - ) + output_currents=["Rload"]) self.assertIsInstance(block, StateSpace) np.testing.assert_allclose(block.A, A) @@ -47,9 +45,7 @@ def test_path_object_preserves_feedback_source_input(self): """Test current sources feedback used to chain models""" block = NetlistStateSpace("filter_without_load.net", output_voltages=["n1"], - output_currents=["Lfilter"], - reduction_mode="symbolic", - ) + output_currents=["Lfilter"]) self.assertEqual(block.input_labels, ["V1", "Bload"]) self.assertEqual(block.output_labels, ["V(n1)", "I(Lfilter)"]) @@ -65,9 +61,7 @@ def test_inline_netlist_string(self): C1 out 0 1u """, output_voltages=["out"], - output_currents=["R1", "V1"], - reduction_mode="symbolic", - ) + output_currents=["R1", "V1"]) self.assertEqual(block.input_labels, ["V1"]) self.assertEqual(block.output_labels, ["V(out)", "I(R1)", "I(V1)"]) @@ -87,31 +81,7 @@ def test_invalid_output_names_raise(self): with self.assertRaisesRegex(ValueError, "Unknown dipole"): NetlistStateSpace(netlist_path, output_currents=["R404"]) -@unittest.skipUnless(importlib.util.find_spec("mumps"), "python-mumps is not installed") -class TestFastReduction(unittest.TestCase): - """Class validation the fast implementation for reduction mode using mumps - using python-mumps""" - def test_fast_mode_matches_symbolic_on_small_case(self): - """Test if fast reduction mode is equivalent to symbolic - skipped if mumps not installed""" - elements = parse_netlist( - """ - V1 n1 0 1 - R1 n1 n2 2k - C1 n2 0 2u - L1 n2 n3 3m - L2 n3 0 5m - R2 n3 0 4k - K1 L1 L2 0.1 - """ - ) - symbolic = CircuitModel(elements, reduction_mode="symbolic") - fast = CircuitModel(elements, reduction_mode="fast") - - np.testing.assert_allclose(fast.A, symbolic.A, rtol=1e-9, atol=1e-12) - np.testing.assert_allclose(fast.B_ss, symbolic.B_ss, rtol=1e-9, atol=1e-12) - # RUN TESTS LOCALLY ==================================================================== if __name__ == '__main__': - unittest.main(verbosity=2) \ No newline at end of file + unittest.main(verbosity=2) From d4095f123adfd59887287c3874e3305520cd64d7 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Wed, 30 Sep 2026 15:26:33 +0200 Subject: [PATCH 15/18] Updated implementation to only run with scipy, removing mumps and sympy --- src/pathsim_rf/netlist_to_statespace.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/pathsim_rf/netlist_to_statespace.py b/src/pathsim_rf/netlist_to_statespace.py index 5667b4a..0cb3743 100644 --- a/src/pathsim_rf/netlist_to_statespace.py +++ b/src/pathsim_rf/netlist_to_statespace.py @@ -351,6 +351,10 @@ def col_iV(name): # ---- inputs: one column per V source, then one column per I source ---- self.input_names = self.V_names + [e.name for e in self.I] self.n_u = len(self.input_names) + if self.n_u == 0: + raise ValueError( + "Circuit must contain at least one independent voltage or current source." + ) # ---- build E (dynamic) and G (algebraic) matrices, and B (input map) ---- n = self.n_total @@ -545,13 +549,13 @@ def _reduce(self): dr, ar = self.diff_rows, self.alg_rows dc, ac = self.diff_cols, self.alg_cols - E_dd = np.array(E[dr, dc], dtype=float) - G_alg_d = np.array(G[ar, dc], dtype=float) - G_alg_a = np.array(G[ar, ac], dtype=float) - G_diff_d = np.array(G[dr, dc], dtype=float) - G_diff_a = np.array(G[dr, ac], dtype=float) - B_alg = np.array(B[ar, :], dtype=float) - B_diff = np.array(B[dr, :], dtype=float) + E_dd = E[np.ix_(dr, dc)] + G_alg_d = G[np.ix_(ar, dc)] + G_alg_a = G[np.ix_(ar, ac)] + G_diff_d = G[np.ix_(dr, dc)] + G_diff_a = G[np.ix_(dr, ac)] + B_alg = B[ar, :] + B_diff = B[dr, :] try: Phi = -self._solve_with_scipy(G_alg_a, G_alg_d, "algebraic subsystem") From ce1cb7c05d32ce1589c4c78acf10397d93f00874 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Wed, 30 Sep 2026 15:29:11 +0200 Subject: [PATCH 16/18] Update example to remove mention of import mode --- docs/source/examples/netlist_import.ipynb | 2 -- 1 file changed, 2 deletions(-) diff --git a/docs/source/examples/netlist_import.ipynb b/docs/source/examples/netlist_import.ipynb index 4e79ce6..35fee01 100644 --- a/docs/source/examples/netlist_import.ipynb +++ b/docs/source/examples/netlist_import.ipynb @@ -138,14 +138,12 @@ " filter_netlist,\n", " output_voltages=[\"n1\"],\n", " output_currents=[\"Rload\"],\n", - " reduction_mode=\"fast\",\n", ")\n", "\n", "state_space_bis = NetlistStateSpace(\n", " filter_without_load_netlist,\n", " output_voltages=[\"n1\"],\n", " output_currents=[\"Lfilter\"],\n", - " reduction_mode=\"symbolic\",\n", ")\n", "\n", "# Dynamic\n", From 36d73b018dd95dd583b79cc777c2a3c909fafd5b Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Sun, 4 Oct 2026 11:13:57 +0200 Subject: [PATCH 17/18] Tests bugfix for NetlistStateSpace for path handling Added test for string input path handling and updated netlist path handling in existing tests. --- tests/test_netlist_statespace.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/tests/test_netlist_statespace.py b/tests/test_netlist_statespace.py index 3d6463b..bbca2fc 100644 --- a/tests/test_netlist_statespace.py +++ b/tests/test_netlist_statespace.py @@ -10,10 +10,12 @@ from pathlib import Path import numpy as np import importlib.util +TEST_DIR = Path(__file__).parent from pathsim_rf.netlist_to_statespace import CircuitModel, NetlistStateSpace, parse_netlist_file, parse_netlist from pathsim.blocks.lti import StateSpace + # TESTS ================================================================================ class TestNetlistStateSpace(unittest.TestCase): @@ -21,16 +23,13 @@ class TestNetlistStateSpace(unittest.TestCase): def test_path_string_matches_manual_circuit_model(self): """Test CircuitModel is correctly used within NetlistStateSpace""" - netlist_path = "filter.net" + netlist_path = TEST_DIR / "filter.net" model = CircuitModel(parse_netlist_file(netlist_path)) model.add_node_voltage_output("n1") model.add_dipole_current_output("Rload") A, B, C, D = model.get_system() - block = NetlistStateSpace( - str(netlist_path), - output_voltages=["n1"], - output_currents=["Rload"]) + block = NetlistStateSpace(netlist_path, output_voltages=["n1"], output_currents=["Rload"]) self.assertIsInstance(block, StateSpace) np.testing.assert_allclose(block.A, A) @@ -43,9 +42,9 @@ def test_path_string_matches_manual_circuit_model(self): def test_path_object_preserves_feedback_source_input(self): """Test current sources feedback used to chain models""" - block = NetlistStateSpace("filter_without_load.net", - output_voltages=["n1"], - output_currents=["Lfilter"]) + block = NetlistStateSpace(TEST_DIR / "filter_without_load.net", + output_voltages=["n1"], + output_currents=["Lfilter"]) self.assertEqual(block.input_labels, ["V1", "Bload"]) self.assertEqual(block.output_labels, ["V(n1)", "I(Lfilter)"]) @@ -71,11 +70,16 @@ def test_inline_netlist_string(self): def test_explicit_missing_path_raises(self): """Test missing .net file provided""" with self.assertRaises(FileNotFoundError): - NetlistStateSpace(Path("missing.net")) + NetlistStateSpace(TEST_DIR / "missing.net") + + def test_path_str_import(self): + """Test if str args allows file to be imported""" + netlist_path = TEST_DIR / "filter.net" + NetlistStateSpace(str(netlist_path)) def test_invalid_output_names_raise(self): """Test invalid output names provided for state space generated""" - netlist_path = "filter.net" + netlist_path = TEST_DIR / "filter.net" with self.assertRaisesRegex(ValueError, "Unknown node"): NetlistStateSpace(netlist_path, output_voltages=["missing"]) with self.assertRaisesRegex(ValueError, "Unknown dipole"): From 7f5e7a0060b5ac4fd12e0b5977c829b8cf9b5474 Mon Sep 17 00:00:00 2001 From: Pimss <44709590+Pimss@users.noreply.github.com> Date: Sun, 4 Oct 2026 11:19:23 +0200 Subject: [PATCH 18/18] typo fix --- tests/test_netlist_statespace.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_netlist_statespace.py b/tests/test_netlist_statespace.py index bbca2fc..abb6a09 100644 --- a/tests/test_netlist_statespace.py +++ b/tests/test_netlist_statespace.py @@ -73,7 +73,7 @@ def test_explicit_missing_path_raises(self): NetlistStateSpace(TEST_DIR / "missing.net") def test_path_str_import(self): - """Test if str args allows file to be imported""" + """Test str netlist_path arg allows file to be imported""" netlist_path = TEST_DIR / "filter.net" NetlistStateSpace(str(netlist_path))