Codecov Report

❌ Patch coverage is 12.75720% with 212 lines in your changes missing coverage. Please review.

📢 Thoughts on this report? Let us know!

Re: the Codecov patch-coverage warning — traced the cause: pyproject.toml's default addopts = "-m 'not slow'" combined with the conftest.py hook that auto-marks everything under tests/evals/ as slow means CI's pytest --cov=... (no -m override in tests_codecov.yml) never runs tests/evals/test_trim_and_linearize_system.py or tests/evals/test_linearize_system_python_control.py — hence utils/linearization.py shows ~9% instead of the ~88-100% these tests actually give it when run explicitly. This isn't new to this PR: all 97 tests under tests/evals/ (including the pre-existing test_steadystate_transient_system.py) are excluded from CI's default run the same way.

Both new test files run in well under a second, so they aren't "slow" in the sense the marker implies — they landed in evals/ just to match test_steadystate_transient_system.py's precedent. Happy to move them into tests/pathsim/ so they run by default in CI (fixing the coverage number for real and giving actual regression protection on this new code), but wanted to check first since it touches the existing evals/ vs regular test-tree convention rather than just this PR. How would you like to handle it?

Thanks for this, and sorry for the slow response — this touches the part of the
core I'm most protective about, so I wanted to go through it properly rather
than give you a quick reply.

Let me start with what I want to take right away, then explain where I want to
take the rest, because the second part is a roadmap decision on my side rather
than anything wrong with your work.

Merging now

Three of the fixes in here are independent of the linearization design and I'd
like them in as a separate, small PR:

1. blocks/_block.py — the DynamicOperator branch in Block.linearize().
   Real TypeError, and the regression test in test_function.py is exactly
   how I'd have written it.
2. utils/register.py — update_from_array on a 0-d array. Confirmed: this
   raises TypeError: len() of unsized object on master. (The second half, the
   1-element array in __setitem__, is a DeprecationWarning under NumPy 2.3
   rather than an error — still worth hardening, just not the same severity as
   the PR description suggests.)
3. tests/pathsim/blocks/test_table.py — the duplicate TestLUT1D class
   shadowing 8 tests is a genuine catch, pathsim.optim.value.Value really is
   gone from the codebase, and 2.5 is the correct expected value: LUT1D
   interpolates linearly through interp1d, so x=1.5 sits between (1,1) and
   (2,4).

If you open that as its own PR I'll merge it.

Where I want to take linearize_system()

Here's the thing I need to be upfront about: pathsim and fastsim are converging
on a shared architecture, and linearization is one of the places where fastsim
has already made the call. Its Operator (src/blocks/operator.rs) derives the
Jacobian from the SSA graph as a pure query — the doc comment there says it
explicitly: "the Jacobian is recomputed per Newton step rather than
linearized-and-reused". And Simulation.linearize() is deliberately
NotImplementedError on that side.

So the target shape for pathsim is: Jacobian is a query, never a mode. That's a
decision that sits above your PR, and it's why I don't want to merge this path
even though the assembly itself is correct — I checked it against analytic
results for gain cascades, closed loops, pure feedthrough and a 2x2 MIMO
StateSpace, and it reproduces the matrices exactly.

Two things in the current implementation are what convinced me the query layer
has to come first, rather than being a later refactor:

The simulation stays linearized. linearize_block() calls block.linearize(t)
for its side effect in order to read op_dyn.Jx/.Ju, and nothing reverts it. So
after sim.linearize_system(...), a subsequent sim.run() silently integrates
the linear surrogate:

nonlinear run,               y(end) = 1.99955
run AFTER linearize_system,  y(end) = 2.47252

The tests here don't catch it because they happen to call Sim.linearize()
before and Sim.delinearize() after.

Worth knowing: the pure API you need already exists one layer down —
Operator.jac(x) and DynamicOperator.jac_x(x,u,t) / .jac_u(x,u,t) in
optim/operator.py, which blocks/dae.py:253 already uses in the hot solve
path. I rewrote linearize_block() against those and got bit-identical matrices
with no mode switch at all.

Blanket numerical fallback. When a block has no operator, the code
finite-differences it. 22 of 36 block modules have no operator, and that set
includes Comparator, Relay, Switch, RNG and Counter — for those a
numerical "Jacobian" is meaningless rather than approximate. Linearizability
should be an explicit per-block property, and blocks that don't have a valid
linear model should raise rather than return a number. That's also how fastsim
separates traceable operators from genuinely opaque ones.

Two smaller things I found while testing, in case they're useful wherever this
lands:

• _state_labels() and _port_ref_labels() count independently, so the label
  sets collide across categories. On a two-integrator chain I get
  state_labels = ['Integrator_0', 'Integrator_1'] and
  output_labels = ['Integrator_0'] where the tapped block is actually the
  second integrator. Since the labels are meant as python-control identifiers,
  that one matters.
• assemble_linear_system resolves break points with inputs._map(port) instead
  of PortReference._get_input_indices(), which is what resizes the register.
  Marking an input port that isn't wired yet (a disturbance input, say) raises
  IndexError: index 1 is out of bounds for axis 0 with size 1.

The plan

Three levels, to_statespace() throughout:

• Operator — jac() / jac_x() / jac_u(), pure. Already there.
• Block — Block.to_statespace(t), defaulting to the operator queries, with
  exact overrides where an operator is absent (Integrator is A=0, B=I, C=I, D=0, no differencing needed) and a hard failure where no valid linear model
  exists. Block.linearize(t) stays as the mode switch but starts returning the
  matrices it computes anyway.
• Simulation — Simulation.to_statespace(inputs, outputs, t) assembling via
  the interconnection form rather than DAG substitution, which handles surviving
  algebraic loops instead of rejecting them.

Returning a StateSpace block at every level also removes the need for the
as_block flag and the label attributes attached from outside, and
Subsystem.to_statespace() falls out for free.

I'm writing this up properly in #195. If you'd like to take the block level —
level 2 above — it's the most self-contained piece and I'd be glad to have it.
No obligation either way; the fixes at the top are worth merging on their own.

Your open questions

evals/ vs tests/: you're right about the cause, and yes — move them into
tests/pathsim/. They run in under two seconds, they aren't slow in the sense
the marker means, and evals/ being excluded from the coverage run isn't a
convention worth preserving for them.

trim(): I'd like to keep the capability but not as a separate solver stack.
steadystate() works by switching all blocks to the SteadyState solver and
reusing _solve() — no new numerics. Trim is that same solve with the unknown
vector extended by the free parameters and the residual extended by the output
targets, so it belongs there rather than beside it. That also settles the
python-control question: it stays pathsim's own, because it's an extension of
machinery pathsim already has.

Splitting the independent bugfixes out as #240 per your review above — thanks for the thorough look. Closing this one since Block.to_statespace()/Simulation.to_statespace() in #239 already cover the linearization design. Following up on trim() in #195.