qmatlib.H: remove a pole that was never really there

[remislp]

The bug

qmatlib.H evaluates

H(s) = QAA + QAF · X⁻¹ (I − exp(−X·tres)) · QFA,      X = sI − QFF

and computed X⁻¹ directly. That breaks down whenever s reaches an eigenvalue of QFF — and the asymptotic root search that calls H walks s across exactly that range, because the roots of the determinant equation are interlaced with those eigenvalues.

It was never really a pole

X⁻¹ (I − exp(−X·t)) is the integral of exp(−Xu) over [0, t]. The offending mode contributes (1 − e^{−μt})/μ, whose limit as μ → 0 is simply t. The singularity is removable; only the factored form has a pole.

integral_expQt evaluates the integral spectrally, taking that limit explicitly with expm1 for accuracy near zero, and H uses it.

How it was found

In the reproduction of Colquhoun, Hatton & Hawkes (2003), 2 of the first 12 fits of one scenario raised LinAlgError("Singular matrix") out of HJC_shut_time_pdf at 10 µM. The rates looked entirely ordinary — α₂ = 2119 against a true 2000.

The pattern turned out to be that both fits had a rate constant sitting at its upper limit of 1e6, which puts an eigenvalue of QAA exactly at −1e6 for the shut-time problem to walk into. All twelve now work.

The worse failure is the silent one

Whether the old form raised or merely returned nonsense depended on whether the eigenvalue was exactly representable:

• round-numbered rates — a limit of exactly 1e6 — land on the pole and raise;
• otherwise the matrix is simply catastrophically ill-conditioned and the answer is quietly wrong.

That second case is why this is a fix rather than a guard: a try/except would have papered over the loud failure and left the silent one.

dW is deliberately untouched

dW contains a similar-looking inv(sI − QFF), and its pole is genuine, not removable — W really does have poles at the eigenvalues, which is why the bisection brackets roots between them by design. Changing it would be wrong.

Tests

19 in tests/test_integral_expQt.py. The integral is checked

• against inv(M)(exp(Mt) − I) on 8 random matrices — agrees to 1e-8;
• against a quadrature, so the identity is verified rather than asserted;
• for continuity as an eigenvalue passes through zero (μ = 1e-3 … 0);
• for giving exactly t at the zero mode.

and H is checked unchanged away from the pole, and finite and continuous through it.

One test failed honestly while I wrote it: I had asserted the old form raises at any computed eigenvalue, and for CH82 it doesn't — the eigenvalue isn't exactly representable, so inv merely becomes ill-conditioned. The test now asserts the property that actually matters (finiteness and continuity) plus a measured condition number, and the docstrings were corrected to match.

Nothing else moves

• 635 passed, 1 skipped — the full suite including the slow simulator tests.
• The asymptotic areas still agree with the 2003 Fortran to 5×10⁻⁷ on all seven time constants and areas, which is that program's print precision.

🤖 Generated with Claude Code