Making Reservoirs Continuous

My last post about ReservoirComputing.jl was the generic tour: what reservoir computing is, and the four models I'd contributed before. This one is about the actual project. This summer I was in the SciML fellowship, SciML being the Julia differential-equation ecosystem run out of MIT CSAIL, with Chris Rackauckas as one of my mentors alongside Francesco Martinuzzi, who maintains the package day to day. My problem was one specific thing: continuous-time reservoir computing. Twelve merged PRs between June and August. This post is the honest version of how it went: the design fights, the bug that made a perfect model look broken, the performance work, and the part where my maintainer caught my AI-generated code not being up to standard.

The problem in one sentence

Every reservoir in the package was a discrete-time map. Input column comes in, state updates in one step, next column. That's how the literature mostly works and it's fine, until you want your reservoir to live inside the SciML ecosystem, where everything else is an ODEProblem. The idea behind issue #397 was simple to state: let any SciML problem be a reservoir. Feed the input as a continuous signal, evolve the state with a real differential equation solver, sample the trajectory, train the readout the same as always.

Francesco (the maintainer) sketched the shape in March: a reservoir built around an AbstractSciMLProblem, an interpolated input \(\mathbf{u}(t)\) as the forcing term, and a sampling rule that turns the continuous trajectory back into discrete states for the readout. The first version targets ODEProblem specifically; SDEs, DDEs, and everything else slot in later because the prob field stays untyped.

What continuous actually buys you

Worth pausing on why this matters, because "just use a smaller dt" is a fair question. The discrete leaky ESN update is an Euler discretization of a continuous system with step \(\alpha\). The canonical continuous form is Lukoševičius 2012, §3.2.6, eq. (5):

$$\dot{\mathbf{x}}(t) = -\mathbf{x}(t) + \tanh\!\bigl(\mathbf{W}_{\text{in}}\,\mathbf{u}(t) + \mathbf{W}_r\,\mathbf{x}(t) + \mathbf{b}\bigr)$$

Notice what's missing: the leak rate. \(\alpha\) is not a parameter of this ODE. It only appears when you Euler-discretize with step \(\Delta t = \alpha\). That single fact rearranged my whole mental model mid-project (more on that in the ContinuousESN section). The practical upshot: with a continuous reservoir you get adaptive solvers for stiff systems, irregular sampling, events and callbacks, and interop with every ODE-based model in SciML for free. The tradeoff is that you inherit everything adaptive solvers do, including some things you didn't want.

Discrete-time stability is spectral radius. Continuous-time stability is operator norm. \(\rho(W_r) < 1\) is a discrete artifact; the continuous contraction condition is \(\|W_r\|_2\). Half the folklore doesn't transfer.

The plan, and the one question I refused to guess on

The fellowship is officially twelve weeks, and my instinct was to break it into one PR per deliverable rather than one per week: each PR a complete, testable unit, stacked and rebased as they merge. Core types in week 1, the ODE extension in weeks 2–4, the model plus validation in 5–7, performance in 8–9, liquid state machines in 10–11, delay equations as the stretch in 12. I posted the whole table on the issue before writing a line of code, with one question flagged above everything else: how should collectstates dispatch?

The catch was subtle. collectstates dispatches on the computer, the thing holding reservoir plus readout. But the thing that makes a model continuous is the reservoir field. Option (a): parameterize the computer on its reservoir type. Cleanest for users. Option (b): a separate CTESN type with its own collectstates, which is exactly what DeepESN already does, with in-repo precedent. I leaned (b) and said so, but the maintainer's original sketch clearly wanted (a), and I'd rather ask than guess wrong and rewrite.

The answer that came back was neither, and it was better than both. A two-level dispatch: the public collectstates(rc, data, ps, st) becomes a one-liner delegating to _collectstates(rc.reservoir, rc, ...), the existing eachcol loop becomes the generic fallback (byte-for-byte the old body), and continuous reservoirs get their own _collectstates method in the extension. No signature churn anywhere. It's not the cleanest solution, as Francesco put it himself, but it's the most flexible one, and the existing eight-odd models route through it without noticing.

PR1: types, and two failures worth learning from

The first PR landed the hierarchy. AbstractSciMLProblemReservoir <: AbstractLuxLayer is the supertype, covering any reservoir whose dynamics come from an AbstractSciMLProblem. Under it sits exactly one concrete type, SciMLProblemReservoir, with five fields lifted straight from DiffEqFlux's NeuralDE: (prob, sampler, tspan, args, kwargs). The prob field stays deliberately untyped, which resolves point 1 of the issue: ODEProblem now, SDE/DDE later without touching the core. Solver settings and tspan live on the struct instead of threading through as arguments, so the existing 4-arg collectstates(rc, data, ps, st) that both train! and the benchmarks library call keeps its signature; the time grid gets synthesized inside from size(data, 2) and the struct's tspan. And the sampling rule becomes a hierarchy, AbstractSampler with TerminalStateSampling as the only v1 method, so trajectory or mean samplers land later as pure dispatch, no API change. All of it is pure src/: the PR adds zero dependencies, because the DiffEq side waits for the extension.

The dispatch refactor itself is small enough to show. The public function becomes a one-liner:

collectstates(rc::AbstractReservoirComputer, data, ps, st) =
    _collectstates(rc.reservoir, rc, data, ps, st)

The old eachcol loop (one _partial_apply per input column, state modifiers threaded between steps) drops down unchanged as the generic _collectstates(res, rc, ...) fallback, and a method on AbstractSciMLProblemReservoir throws for now with an error pointing at the extension that doesn't exist yet. Every existing model routes through the fallback without noticing; DeepESN's own collectstates override keeps winning by most-specific dispatch and stays byte-for-byte as it was. I shipped a regression scaffold asserting the public two-level path and an inline replica of the pre-refactor body produce elementwise-identical state matrices, 18 checks, plus 14 more for construction, the Lux initialparameters/initialstates interface, and the missing-extension error path. Two things went wrong along the way, both educational.

First, my regression test for the dispatch refactor captured hashes of the state matrices on one machine and asserted against them. The hashes matched on exactly one Julia version, the 1.12 that produced them, and failed on 1.10, 1.11, and pre-release. Float summation order differs across versions, so hashing output is a portability trap. The fix was embarrassingly simple: run both the public dispatch and an inline replica of the pre-refactor body in a single Julia session and assert elementwise equality. Portable by construction. Every "golden value" test I write now has to answer: golden on which Julia?

Second, a typed outer constructor. I typed the sampler argument, and type inference recursed forever. The reason: a typed outer constructor is strictly more specific than the inner one @concrete generates, so the 5-arg recursive call inside it resolves to itself. DiffEqFlux's NeuralDE family takes the same args untyped for the same reason. That's now documented inline in the source, which is the silver lining of spending an afternoon on a stack overflow.

It merged in two days. I've since learned that's not the normal speed. Enjoy it when it happens.

PR2: the extension, where the real subtlety lived

The whole DifferentialEquations side lives in a package extension, RCODEReservoirExt, triggered by loading SciMLBase + DataInterpolations + any solver package. Julia's extension mechanism means the base package never touches the ~150 transitive dependencies that come with an ODE stack, so an ESN-only user never compiles a line of OrdinaryDiffEq. Prediction gets the same two-level treatment as state collection: public predict delegates to _predict on the reservoir type, with a teacher-forced method taking an input matrix and an autoregressive one taking an initial point plus a step count, both solving the reservoir ODE over per-window spans and hitting the readout at the window ends.

Inside _collectstates, the pipeline is: build the continuous input \(\mathbf{u}(t)\) from the discrete columns, remake the problem with the struct's tspan and the readout-agnostic reservoir parameters, solve with saveat on the synthesized sample grid, push the trajectory through the sampler, then apply state modifiers once on the resulting (state_dims, T) matrix. Because the extension owns the sample grid, three solve keywords (saveat, save_everystep, dense) are rejected at construction with an explicit error rather than silently overridden. Same defensive posture elsewhere: a user parameter named :input collides with the reserved key the extension injects for the forcing signal, so that errors loudly instead of shadowing; prob.p is validated to NamedTuple / nothing / NullParameters (anything else would explode later inside the solver, far from the cause); degenerate tspans are rejected up front. Guardrails at construction, not stack traces at midnight.

My first version used linear interpolation for the input signal. It passed all the tests. It was also wrong. With linear interpolation, the state saved at step \(k\) depends on data[:, k+1]: the input from the next window bleeds into the current state for any non-Euler solver, because the interpolant inside window \(k\) reaches toward the endpoint at \(k+1\). That contradicts the documented semantics: states[:, k] is the state after processing input column \(k\), full stop. The fix was zero-order hold, each column held constant across its window, which also returns a view so the input function is allocation-free in the ODE hot path. A one-line semantic bug that no shape check would ever catch. The Euler-equivalence test caught it because with Euler there is no "inside the window."

Other things review caught: a collect() on the initial state that silently destroyed SVectors and errored on scalar u0 (the loop only ever reads, so the copy was pointless); an output buffer allocated with the wrong eltype, which would quietly convert a readout's results; an @assert where the rest of the API raises ArgumentError. And a small human moment: after approval, I asked which version to bump (0.12.24, or 0.13.0 for a "major" feature) and got to make the call on the package's semver because the reviewer was busy. Felt like a bigger deal than it was. That's most things in open source.

ContinuousESN, the cold-start bug, and the best typo of my summer

PR3 was supposed to be the big validation push: a CTESN model, Mackey-Glass, Lorenz, the works. It collapsed into something much thinner, and the collapse was the interesting part.

I had assumed "continuous ESN" meant the Anantharaman 2021 CTESN paper, a parametric surrogate with RBF-interpolated readouts. When I actually dug up the Lukoševičius reference Francesco pointed me at, eq. (5) turned out to be the whole story: no \(\alpha\) in the ODE, bounded state by construction (\(\tanh\) keeps \(\mathbf{x} \in (-1,1)^N\) with no extra contraction argument), and the leak rate reappearing purely as the Euler step. So the "CTESN model" shrank from a research-paper implementation to a thin convenience constructor that pre-bakes eq. (5). The redesign dropped half the open questions in one move. Reading the actual paper beats speculating from the abstract; I keep having to relearn this.

What shipped: ContinuousESN, a 3-field model, (reservoir, states_modifiers, readout), that mirrors ESN exactly, with the ODE substance in a ContinuousESNCell whose default right-hand side is the leaky integrator, evaluated in place with mul! and fused broadcasts for zero allocations per call. The constructor takes tspan and the solver as positional arguments, ContinuousESN(3, 300, 3, (0.0, 100.0), Tsit5(); use_bias = true), and passes them through to the same SciMLProblemReservoir/extension machinery as any user-built problem. The load-bearing test is Euler equivalence again: solve the continuous cell with forward Euler at step \(\alpha\) and it must match the discrete ESNCell leaky update to machine precision, at \(\alpha \in \{1.0, 0.5\}\). That test is the bridge between the two formulations; if either implementation drifts, the other catches it.

Then the smoke test. Teacher-forced training fit beautifully at NRMSE 0.0004. Then the autoregressive rollout on Lorenz came back at NRMSE ~1.5 at every horizon. A model that fits that well cannot forecast that badly. I dug into the autoregressive path in the extension and found it:

current_state = res.prob.u0, which is zeros, every time. The terminal reservoir state after training never made it into the predict-time initial condition. Every rollout started cold, from \(\mathbf{0}\), regardless of how well the readout was trained. Four lines of manual fix (re-run collectstates, take the last column, remake the problem with it as u0) and the same model did this:

Horizon (Lyapunov times)	Steps	NRMSE (cold)	NRMSE (warmed)
1 \(t_\lambda\)	55	1.489	0.157
2 \(t_\lambda\)	110	1.496	0.118
3 \(t_\lambda\)	166	1.743	0.112
4 \(t_\lambda\)	221	1.696	0.210
6 \(t_\lambda\)	331	1.644	0.729
8 \(t_\lambda\)	442	1.471	0.981

(\(N = 300\), Lorenz-63 at \(\sigma = 10, \rho = 28, \beta = 8/3\), \(dt = 0.02\), \(\|W_r\|_2 = 0.9\), Tsit5 at \(10^{-6}\).) Three Lyapunov times of usable forecast from a 300-neuron continuous reservoir, before chaotic divergence eats it. And the failure mode before the fix is the instructive part: cold autoregressive forecasting looks like a broken model when it's really a missing warm start. Same weights, same readout, one swapped line: NRMSE 1.5 vs 0.11.

The review on this PR also produced the most useful feedback I got all summer. My docstrings explained routing and dispatch; my code was dense with comments; I'd been prototyping with autonomous coding agents. Francesco's response: docstrings are for users ("when I use + I am not interested in how the dispatch was made"), comments should mostly not exist, and, verbatim: "it's ok to use llms for coding, but make sure the code is up to standard." That last line stung exactly as much as it should have. I told him honestly that the agents were the problem and dropped back to writing everything myself for the rest of the PR. The agents are fine for exploration. They are not fine for the thing that gets merged, because "up to standard" is a taste-level judgment that neither SciML style guide prompts nor AGENTS.md files reliably produce. I now keep an agent setup with the SciML guide forced at project level, but I've stopped trusting it to write final code.

The warmup saga, or: why is this needed at all?

The cold-start fix above was a manual workaround. Making it a real API took another month, and the path there is a good example of a maintainer asking the right question instead of rubber-stamping my design.

My proposal: a warmup_data keyword on predict that internally teacher-forces and seeds \(\mathbf{u}_0\). I built an investigation-only harness (no package changes, red CI on purpose) and ran the full matrix at the ContinuousESN scale. The headline numbers: valid prediction time went from 0.22 \(t_\lambda\) cold to 4.87 \(t_\lambda\) with even a 10-step warmup. And the seed ablation, which I ran mostly for completeness: warm with the wrong state (random, or a shuffled terminal state) gives NRMSE around 13, worse than zeros. Warm start isn't just "some reasonable state," it's a dynamically consistent one. Chaos has no sense of humor about initial conditions.

Francesco's reply opened with: "my main question in this issue was why is the warmup needed in the first place. after all, ESN doesn't use it." Which is the correct question, and asking it dissolved the design. The discrete models don't need a warmup concept because after train! their state struct carries the terminal state forward automatically. Continuous didn't need a warmup feature. It was missing the thing discrete already does. Not a feature to add. A parity gap to close.

So the fix, when it landed as #499, had no new keywords at all. After collectstates or train, the continuous path writes the raw ODE terminal state into st.reservoir.carry, stored as the same one-element NamedTuple shape the discrete models use (which means resetcarry! works on it unchanged, no special-casing), and the next collectstates or autoregressive rollout seeds its initial condition from it instead of prob.u0. Discrete and continuous now have identical state-threading semantics; the "warmup" concept evaporates because nothing needs warming anymore. The whole API discussion ended in a PR with no API. My favorite kind.

The fallout was embarrassing in a useful way: the continuous tutorials in the docs had been silently plotting collapsed forecasts the whole time, because the tutorial code rebuilt a second model and copied only the readout, so every plot started cold. #506 fixed the tutorials, and now the Lorenz plots in the docs actually reconstruct the attractor.

PR4: performance, driven by numbers first

I profiled before optimizing, and the profile eviscerated my assumptions. The plan had said "in-place RHS with mul! + preallocated caches" as the headline perf lever. Reality: the RHS was already zero-allocations per call, and it was irrelevant. The real numbers on a single collectstates at \(N = 500\), \(T = 1000\):

• Adaptive Tsit5: ~13.5 s, with 330,537 RHS evaluations (21,262 accepted steps, 33,827 rejected).
• Fixed-step Euler at dt=1: ~42 ms, 1,001 RHS evaluations.
• Discrete ESN, same problem size: ~42 ms.
• Stage timers: solve is ~100% of wall time on every continuous path.

320× over Euler. 350× over the discrete model. That's not an implementation problem. It's adaptive error control doing what it's designed to do on a stiff-ish recurrent system. You cannot micro-optimize your way out of 330k RHS calls; you have to change what the solver is asked to do. Three levers, each its own PR, each with before/after measurements:

Integrator reuse (#479). The autoregressive path rebuilt the entire ODEProblem/integrator/interpolation object per step: ~76 KB per solve, ~3 KB per problem, ~0.7 KB per interpolation, times 250 steps. The refactor builds the integrator once and reinit!s it per window with erase_sol=true, reset_dt=true (each window must stay an independent solve, so no state may leak, including the adaptive stepper's cached dt), then advances with solve!(int) instead of allocating a fresh solution, and reads the terminal state straight off integrator.u. The per-window input signal moves into a tiny mutable ConstantInputWindow whose single field gets rebound in place, so the RHS's p.input(t) contract is untouched and the type-specialized solver caches never see a new parameter type. A bundled change swaps the cell-path constructions to ODEProblem{true, FullSpecialize}, which sidesteps AutoSpecialize's FunctionWrappersWrapper boxing; five distinct allocation types vanish from the profile. Result on the benchmark: memory 10.68 → 2.18 MiB (-80%), allocations 50,008 → 5,892 (-88%), wall-time σ cut 63%, GC eliminated, and the output bit-identical to master, max abs diff = 0.0, same hash. That last check was Francesco's first question ("have you checked consistency of results before and after?"), and reinit! with both flags genuinely reproduces the fresh-solve trajectory. Median wall barely moved, and I reported that honestly rather than claiming a speedup the data didn't show. The wins were variance, memory, and tail latency (18 KB/step instead of 82), which is what allocation removal actually buys.

Sparse Jacobians (#480). For implicit solvers like Rodas5, dense finite-difference Jacobians dominated: ~334k RHS calls and ~3k Jacobian evaluations at \(N = 100\). The fix is an opt-in use_jac_prototype = true flag that hands the solver the analytic sparsity pattern, which is exactly \(\mathrm{pattern}(W_r) \cup \mathrm{diag}\), because

$$\mathbf{J}(\mathbf{x}) = -\mathbf{I} + \operatorname{diag}\!\left(1 - \tanh^2(\mathbf{z})\right) \mathbf{W}_r$$

and the \(-\mathbf{I}\) means the diagonal must be in the pattern even when \(W_r\)'s diagonal is empty. Omit it and the colored finite-difference skips those columns and produces a wrong Jacobian. Result: 1.72× wall time on the Rodas5 workload, allocations down 99.4%. Explicit solvers are untouched; the Jacobian is dead information to Tsit5. One review exchange worth recording: post-toggle output differs from dense by ~2.3e-3, which looks alarming and isn't. Different FD coloring means different perturbation columns and step sizes, both trajectories are within tolerance of the truth, and Tsit5 (no Jacobian at all) is bit-identical between the two. Knowing why a numerical result is allowed to differ is the difference between a regression and a rounding story.

Extension hygiene. When the sparse path needed SparseArrays, the first draft just added it to the extension trigger list. Francesco pushed back on load weight, so it became a second extension, RCODEReservoirSparseArraysExt, activated only when SparseArrays is loaded, with a generic nothing fallback living in src/. OrdinaryDiffEq-only users see zero change in what triggers the base extension. Is this over-engineering for one type? Maybe. But dependency weight in Julia is a compile-time tax everyone pays forever, and once you've seen the trigger-list discipline it's hard to unsee.

LSM, out of order, with a naming fight

The roadmap said spiking reservoirs were weeks 10–11, after performance, before a DDE stretch. Mid-fellowship the priority shifted: the planned slot got redirected toward ESN + universal differential equations work, but the design pass for liquid state machines was already done and the continuous infrastructure from PRs 1–4 made it cheap to land anyway. So LSM arrived out of order, as its own issue plus PR, in August. The DDE stretch never happened. Plans are hypotheses.

The architecture: a 3-field LSM mirroring ContinuousESN, built on an LSMCell holding a LIFNeuron population (membrane time constant \(\tau_m\), rest/reset/threshold voltages, refractory period \(\tau_{\text{ref}}\)) with exponential synaptic currents so the recurrent coupling lives in the continuous state rather than as instantaneous jumps at spike times. That single choice matters more than it sounds: adaptive solvers can integrate a smooth synaptic current, but a step jump at a detected threshold breaks the error estimator's assumptions and interacts badly with the autoregressive path. Threshold crossings and refractory bookkeeping are enforced by a single VectorContinuousCallback, one root-find per spike per neuron with the reset applied atomically at the crossing, instead of polling in the RHS, where a spike could hide between steps. Around the cell, three pluggable pieces: input encoders (direct current injection, and Poisson rate coding with event times precomputed up front so the adaptive stepper can never re-draw the RNG mid-step; encoder state lives in st.encoder, never in globals), spike features read off the reservoir (binned spike counts, exponential filters, raw membrane voltage), and a Dale-compliant dale_sparse initializer that bakes the E/I sign structure and balance into the connectome rather than hoping training preserves it. And because it all rides the same reservoir path, collectstates returns the same (state_dims, T) matrix the ridge readout expects, so the LSM is shape-compatible with every downstream tool the ESN already has. All of those decisions were made in a design comment before any code, and all of them stuck.

The naming fight was real. LIFCell collides with LIFESN (Local Information Flow), which is a completely unrelated model already in the package. Two identical acronyms in one namespace. Francesco's first suggestion was "be extremely clear in the docs," then mid-review he saw how the code read together and flipped to LIFNeuron. And a bigger catch: my AbstractSpikeReadout family looked extensible (abstract types, multiple dispatch, the works) but the implementation assumed the built-in types through direct field access, so a custom subtype would construct fine and blow up in the solver. His call: either define real dispatch points or make it explicitly closed. We made it closed. Unknown subtypes now throw ArgumentError at construction, and AR support defaults to false for anything not explicitly verified. Fake extensibility is worse than none; the abstract type is a promise, and unkept promises belong in neither APIs nor docstrings.

Final naming: LIFNeuron, AbstractSpikeFeature, SpikeCountFeatures, ExponentialSpikeFilter, MembraneVoltageFeature. The only "readout" in the stack is the ridge-trained LinearReadout, same as every other model. Sixty-seven tests, merged in nine days.

The detour that wasn't on the plan

In July I also landed #474, which had nothing to do with continuous time. LinearSolve went from optional package extension (RCLinearSolveExt, deleted) to hard dependency with the ridge logic moved into src/train.jl; the default ridge path stops forming normal equations and instead solves the Tikhonov-augmented least-squares system with a multi-column right-hand side through LinearSolve's QRFactorization(), numerically better conditioned than \((RR^\top + \beta I)^{-1}\) and one factorization for all channels. The training entry point became non-mutating train, with the objective (RidgeRegression, renamed from StandardRidge on Francesco's mid-review suggestion because the old name read badly next to a solver keyword) cleanly separated from the solver algorithm, and train! kept as a depwarn wrapper that maps the old positional train_method over. The deprecations live in a dedicated deprecated.jl with matching tests, so the v1.0 cleanup is "delete one file" (Invia's deprecation layout, adopted wholesale). This PR technically wasn't mine to do; it came from community issues #473 and #367. But breaking-API work is easier to land while you're already the person bumping versions every week.

The tail of the fellowship was smaller polish. A StateSpaceSet package extension lets train and teacher-forced predict take trajectories from DynamicalSystems.jl directly, points as time, no manual Matrix round-trip, with predict handing back a StateSpaceSet too; the conversion reuses the existing (dims, T) matrix path underneath, and autoregressive predict stays unchanged since initialdata is already a point vector. Francesco's review: "very clean and quick solution" plus one inline doubt on the code, resolved same day. Alongside that, a README restructure into ESN variants / wrappers / continuous-time (the model list had grown past "one flat table works"), and a CONTRIBUTING.md capturing the conventions this codebase actually uses: conventional commits, the model/cell shape, the closed-until-proven-extensible rule. Writing down the unwritten rules was more work than expected and more satisfying than most of the code.

The scoreboard

PR	What	Merged
#446	Core types + two-level _collectstates dispatch	June
#450	RCODEReservoirExt: continuous collect/predict, ZOH input	June
#456	ContinuousESN + Lorenz tutorial	July
#474	LinearSolve default ridge, non-mutating train	July
#479	Integrator reuse across AR windows (-88% allocs, bit-identical)	July
#480	Opt-in sparse jac_prototype (1.72× on Rodas5)	August
#497	LSM: LIF neurons, encoders, spike features, Dale init	August
#499	Persist continuous ODE state in st.carry (the cold-start fix)	August
#506–#527	Docs: tutorial fixes, LSM tutorial, README, CONTRIBUTING	Aug – Sep
#529	StateSpaceSet extension	open

Against the original twelve-week table: PRs 1–4 landed roughly on schedule and in order. PR5 jumped the queue and PR6 never ran. The single best deliverable, the carry fix in #499, didn't exist in any version of the plan; it came out of a bug I found by accident and a design conversation where the right question killed a proposed API. If I scored the fellowship on plan adherence it'd be a B. On what actually landed, I'll take it.

What the fellowship actually taught me

Three things, in descending order of how long they took to learn:

Design conversations are the work. The two-level dispatch, the ZOH switch, the carry-instead-of-warmup reversal, the closed-LSM call: every one was settled in an issue comment before code, and every one would have been expensive to discover in review. The single most valuable habit I built was writing the design decision with its two alternatives and my lean, then asking to be overruled. It gives the maintainer a cheap decision to make instead of an expensive one to reverse.

Numbers or it didn't happen. Every perf claim in this post has a table behind it, every parity claim has a max abs diff, and the one optimization that didn't speed anything up got reported as not speeding anything up. The profile-first discipline killed my pet optimization (RHS tuning) in one afternoon and redirected the whole perf phase toward what actually mattered (NFE counts, integrator churn, Jacobian structure).

Taste is the bottleneck. The gap between "works and passes tests" and "belongs in the package" is entirely taste: docstrings written for the reader instead of the author, error messages that name the fix, naming that doesn't lie about scope, knowing when an abstraction is a promise you can keep. Style guides help at the margin. Agents help at the margin. The rest is iteration against someone whose taste is better than yours, which is the actual privilege of working under a good maintainer.

References

Lukoševičius, M. (2012). A practical guide to applying echo state networks. Neural Networks: Tricks of the Trade, 2nd ed. (§3.2.6, eq. (5)).
Maass, W., Natschläger, T., & Markram, H. (2002). Real-time computing without stable states: a framework for neural computation based on perturbations. Neural Computation 14(11).
Anantharaman, R. et al. (2021). Continuous-time echo state networks. CHAOS 31.
Invia (2022). Deprecating code in Julia. The deprecation layout used in #474.
The fellowship tracking issue, where every design decision in this post lives, with the arguments.