pipekit-cycle

Summary¶

Decision. pipekit-cycle 0.0.2 stays as it is: a separate package beside pipekit-array, depending only on pipekit, with no JAX and no domain code. Its protocols, cycles and helpers are the contract algorithm libraries target. Two additions close the gaps the running example hits, and one more is already proposed in Report 13.

• A. Stateful pipelines as cycle steps. pipekit.Sequential threads a carry state when it holds a stateful step, but it is not a StatefulOperator, so Cycle and its siblings reject it with a TypeError. Accept any Sequential that is stateful as a step.

• B. Observation inputs on the cycle’s clock. obs_source(i) receives a call-local index that restarts at 0 on every call, and R comes only from the state. Let the index follow state.cycle_count (B1), and let a source set R per step (B2, the overview’s obs_err_source).

• Report 13 D1. Ensemble histories. EnsembleCycle and EnsembleDACycle keep no trajectory; Report 13 adds save_history to both, so cycle outputs can be written and indexed.

End state. A script composes a step from stateless and stateful operators, cycles it, and restarts it from the last state without re-reading old observations. Algorithm libraries supply forward models, observation operators and analysis steps by structural typing. Every cycle can keep its history, which Report 13’s state profile indexes.

Ownership rules¶

pipekit-cycle owns the shape of a run over time; the carrier, the physics and the algorithms belong to whoever supplies the steps (Figure 1).

pipekit-cycle owns the cycles, the DA cycles, the protocols, the generic observation and forward helpers, the carry states and the adjoint specs. It imports only pipekit, and never an array library, JAX or a domain package.

Neighbours follow the stack rules. Two matter here: algorithm libraries (filterax, vardax, plumax) satisfy the protocols by structural typing without importing pipekit-cycle, and scripts own the carrier, the initial state, the observations and every file a run writes (Report 13).

What ships¶

The eleven further protocols (the italic rows of Table 1) are the shared vocabulary for the algorithm libraries: a tangent-linear model and its adjoint (M′, M*), a low-rank error subspace B ≈ LLᵀ, a reduced-order model with encode, decode and step, covariance localisers, ensemble inflators, adaptive step schedulers, and ensemble Kalman inversion as an IterativeProcess. None is used by a cycle in pipekit-cycle; they exist so that, say, filterax’s inflator and vardax’s error subspace agree on one signature.

The running example: a methane plume advected along a transect¶

The problem. A methane enhancement $x \in \mathbb{R}^N$, in ppb above background, sits on a periodic transect of $N = 40$ cells. Wind carries it downstream one hour at a time. A column sensor sees every fourth cell, so there are $p = 10$ observations per hour, for $T = 10$ hours. The plume starts as a Gaussian of 50 ppb peak at cell 8.

The model. Let $M \in \mathbb{R}^{N \times N}$ be first-order upwind advection with Courant number 0.4, $M = 0.6 I + 0.4 S$ with $S$ the one-cell downstream shift; $Q = 0.01 I$ the model-error covariance; $H \in \{0, 1\}^{p \times N}$ the selection of every fourth cell; and $R = 0.1 I$ the observation-error covariance. The truth and the observations follow (1):

The filter starts from $m_0 = 0$ with $P_0 = 100 I$: it knows nothing about the plume. The data are synthetic, so every stage has an exact answer to check against.

Stage 1 · Forecast only (pipekit-cycle)¶

TL;DR. Push the plume downstream for six hours and keep every second state.

Problem. With no observations the forecast is $x_{t} = M x_{t-1}$, repeated $n$ times; a history with stride $s$ keeps the states after steps $1, 1 + s, 1 + 2s, \ldots$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# Shapes: N = 40 cells
import pipekit_cycle as pc

class Advect:                                  # satisfies pc.ForwardModel    user
    dt = 1.0
    state_signature = None
    # x_t = M x_{t−1}
    def step(self, x, dt):
        return M @ x                           # (N,) → (N,)

fc = pc.Cycle(Advect(), n_steps=6, save_history=True, history_stride=2)  # today
x6, state = fc(x_true_0, pc.DAState())         # (N,) → (N,)                 today
len(fc.history)                                # → 3: after steps 1, 3, 5    today
[int(np.argmax(c)) for c, _ in fc.history]     # → [9, 10, 10]: the peak moves
state.t                                        # → 0.0: Cycle keeps no clock

A plain Cycle wraps the ForwardModel in an internal adapter and threads the state through untouched, so DAState.t stays at 0; only the DA cycles advance it.

Stage 2 · A Kalman filter as a DACycle (pipekit-cycle)¶

TL;DR. Assimilate the ten hours of observations with an exact Kalman filter, and check that the cycle changes nothing about the arithmetic.

Problem. The Kalman filter carries a mean $m_t$ and a covariance $P_t$:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Shapes: N = 40, p = 10, T = 10; the carrier is (m, P): (N,), (N, N)
class KFForward:                               # pc.ForwardModel              user
    dt = 1.0
    state_signature = None
    # m_b = M m;  P_b = M P Mᵀ + Q
    def step(self, carrier, dt):
        m, P = carrier
        return M @ m, M @ P @ M.T + Q          # ((N,), (N, N)) → same

class KalmanAnalysis:                          # pc.AnalysisStep              user
    # K = P_b Hᵀ (H P_b Hᵀ + R)⁻¹
    # m = m_b + K (y − H m_b);  P = (I − K H) P_b
    def __call__(self, forecast, obs, *, obs_op, obs_err_cov):
        m, P = forecast
        Hm = obs_op.linearize(m).H                                   # (p, N)
        K = P @ Hm.T @ np.linalg.inv(Hm @ P @ Hm.T + obs_err_cov)  # (N, p)
        return m + K @ (obs - obs_op(m)), (np.eye(N) - K @ Hm) @ P

da = pc.DACycle(KFForward(), pc.LinearObs(H), KalmanAnalysis(),
                obs_source=pk.Lambda(lambda i: y[i]), n_steps=T)    # today
(m, P), state = da((m0, P0), pc.DAState(obs_err_cov=R))           # today
# → max |m − m_hand| = 0.0, max |P − P_hand| = 0.0
# → RMSE against the truth 0.125 ppb (prior RMSE 17.1 ppb)
state.t, state.cycle_count, len(da.history)   # → 10.0, 10, 10       today

Each cycle steps the forward model, calls obs_source(i) for the step’s observations, and hands both to the analysis with obs_err_cov read from the state. The result is bit-for-bit the hand-written filter. Two conventions follow from the code and are worth stating once:

• The covariance travels in the carrier. AnalysisStep.__call__ receives the forecast, the observations, the observation operator and R, but never the state. Anything the analysis updates, such as $P$ or an ensemble, must be part of the carrier; here it is a (m, P) tuple.

• ObservationOperator requires linearize. The protocol check fails for a class with only __call__. LinearObs and IdentityObs return themselves; CallableObs takes an optional linearize_fn and raises NotImplementedError without it.

Stage 3 · The same problem as a stochastic EnKF (pipekit-cycle)¶

TL;DR. Replace the exact covariance with an ensemble and check that the ensemble mean approaches the Kalman mean as the ensemble grows.

Problem. With $N_e$ members $x^{(i)}$, the forecast covariance is the sample covariance of the anomalies, and every member is updated against its own perturbed observations, as in the overview ((11)). For a linear-Gaussian model the EnKF mean converges to the Kalman mean $m_t$, and its error falls as $N_e^{-1/2}$ once the ensemble is large enough that the gain is well estimated.

1
2
3
4
5
6
7
8
9
10
# Shapes: N_e members, each (N,); the carrier is a list of N_e arrays
class EnKF:                                    # pc.AnalysisStep              user
    # K = P_b Hᵀ (H P_b Hᵀ + R)⁻¹;  x⁽ⁱ⁾ ← x⁽ⁱ⁾ + K (y + ε⁽ⁱ⁾ − H x⁽ⁱ⁾)
    def __call__(self, forecast, obs, *, obs_op, obs_err_cov): ...   # list → list

eda = pc.EnsembleDACycle(Member(), pc.LinearObs(H), EnKF(),
                         obs_source=obs_source, n_steps=T, n_members=N_e)  # today
members, state = eda(members_0, pc.DAState(obs_err_cov=R))   # N_e × (N,)     today
# RMS |mean(members) − m| , median over 10 seeds; KF posterior std 0.193 ppb
# → N_e =  50: 0.367   N_e = 200: 0.063   N_e = 800: 0.027

From 200 to 800 members the gap shrinks by a factor of 2.3, close to the factor of 2 that $N_e^{-1/2}$ predicts: the log-log slope is −0.61, and a seeded bootstrap over the ten seeds puts its 99 % interval at [−0.96, −0.38], which contains −0.5. At 50 members it is six times larger than at 200: with no localisation, 50 members estimate a 40 × 40 covariance poorly, and the gain carries that error. Inflation and localisation are the analysis step’s business; pipekit-cycle names their signatures (EnsembleInflator, CovarianceLocalizer) but does not apply them.

Stage 4 · Swap the forward model (pipekit-cycle)¶

TL;DR. Replace the physical model with an emulator and leave the cycle unchanged.

Problem. The cycle only needs something that satisfies ForwardModel. An emulator $\hat{M}$ trained to reproduce $M$ should therefore drop in with no other change, and an exact emulator must give the same analysis.

1
2
3
4
5
6
7
8
# Shapes: carrier (m, P)
emulator = pc.CallableForward(lambda carrier, dt: step_fn(carrier), dt=1.0)  # today
neural = pc.NeuralForward(model_op=pk.Lambda(step_fn), dt=1.0)               # today
for fwd in (emulator, neural):
    da = pc.DACycle(fwd, pc.LinearObs(H), KalmanAnalysis(),
                    obs_source=obs_source, n_steps=T)                       # today
    (m_swap, _), _ = da((m0, P0), pc.DAState(obs_err_cov=R))
# → max |m_swap − m| = 0.0 for both

NeuralForward(model_op, dt, state_signature=None) is a thin, working wrapper, not a stub: step calls model_op(state) and ignores dt, so any pipekit Operator (a registry-loaded model, a Sequential) becomes a forward model.

Stage 5 · A Gauss–Newton retrieval as a Recurrence (pipekit-cycle)¶

TL;DR. Retrieve a methane column enhancement from five-band transmittances by iterating until the update is negligible.

Problem. Beer–Lambert absorption gives transmittances $f_b(\alpha) = \exp(-k_b \alpha)$ in bands $b = 1, \ldots, 5$ with absorption coefficients $k_b$, for a column enhancement $\alpha$. With measured $y_b$ and Jacobian $J_b = -k_b \exp(-k_b \alpha)$, Gauss–Newton iterates (3),

and stops when $|\alpha_{j+1} - \alpha_j| < 10^{-8}$.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# Shapes: α scalar; y, f(α), J: (5,)
f = pc.CallableObs(lambda a: np.exp(-k * a),
                   linearize_fn=lambda a: -k * np.exp(-k * a))       # today

class GaussNewton(pk.StatefulOperator):                               # user
    # α ← α + Jᵀ (y − f(α)) / (Jᵀ J)
    def _apply(self, a, s):
        J, r = f.linearize(a), y_b - f(a)                     # (5,), (5,)
        step = float(J @ r / (J @ J))
        return a + step, pc.IterationState(iter_count=s.iter_count + 1,
                                           residual_norm=abs(step))

gn = pc.Recurrence(GaussNewton(),
                   condition_op=pk.Lambda(lambda s: s.residual_norm < 1e-8),
                   max_iters=50)                                      # today
a_hat, s = gn(0.0, pc.IterationState())                               # today
# → α̂ = 120.014 (truth 120), gn.iters = 6, gn.converged = True
# → Jᵀ (y − f(α̂)) = 9.9e-20: the normal equations hold at α̂

Recurrence(step_op, condition_op, max_iters=1000) calls condition_op on the state after each step and stops on a truthy result; it records iters and converged on the operator. The condition must be an Operator, so a bare function is wrapped in pk.Lambda.

Stage 6 · Restart the filter, and compose a step (the two gaps)¶

TL;DR. Run the ten hours as two five-hour calls, as an operational cycle restarted from its last state would, and build a cycle step from a Sequential.

Problem. Restarting a filter from its last carrier and state must give the same analysis as one uninterrupted run, provided the second call reads observations $y_5, \ldots, y_9$. A step that mixes a stateless operator and a stateful one is a Sequential, which pipekit core threads state through.

1
2
3
4
5
6
7
8
9
10
11
12
# Shapes: as stage 2
da5 = pc.DACycle(KFForward(), pc.LinearObs(H), KalmanAnalysis(),
                 obs_source=obs_source, n_steps=5)                   # today
c, s = da5((m0, P0), pc.DAState(obs_err_cov=R))                     # hours 0–4
(m_restart, _), s = da5(c, s)                                       # hours 5–9?
# indices passed to obs_source → [0, 1, 2, 3, 4, 0, 1, 2, 3, 4]
# → max |m_restart − m| = 15.4 ppb, although s.cycle_count = 10

step = pk.Sequential([pk.Lambda(lambda x: M @ x), Count()])         # today
step(x, pc.DAState())          # → threads the state: cycle_count 1  today
pc.Cycle(step, n_steps=3)
# → TypeError: step_op must be a StatefulOperator or a ForwardModel, got Sequential.

The restart silently re-assimilates the first five hours of observations and misplaces the plume by up to 15 ppb, with no warning: the state’s cycle_count is 10, but obs_source sees the call-local index. BFNCycle’s docstring already tells users to “supply a stateful or clock-aware obs_source” for a global clock. Workstream B makes the clock part of the cycle. The Sequential case is workstream A.

The scoping report’s names, checked¶

Table 2 checks the stack names the scoping report’s examples use.

Current-state audit¶

Two of the scoping report’s pieces need action; the rest ship and work (Table 3), and Figure 2 counts the page’s calls by marker. The audit read references/geostack/pipekit at ef3c17f on 13 October 2026.

Smaller findings:

• The README is a release behind. Its table says “What’s shipped (v0.0.1)” and omits BFNCycle, the adjoint specs and the eleven further protocols. The CHANGELOG for 0.0.2 (11 July 2026) lists them.

• The latent-space design is not built. docs/design/latent.md (draft 0.1.0) proposes a latent.py with Encoder, Decoder, LatentForwardModel, LiftedObservationOperator, EncodedForwardModel, LatentDAState and a LatentDACycle in three flavours (strong-latent, prior-only, hybrid). Only the LatentMap protocol ships.

• SmootherCycle cannot see the window’s start state. It steps the model window times and passes the forecasts after each step, so a strong-constraint 4D-Var analysis, whose control is the start state, must carry it inside the carrier.

• WindowState is unused, and Recurrence never touches IterationState: the step owns the counters, as in stage 5.

• Sequential threading works the other way round. A cycle inside a Sequential is fine (Inc() | Cycle(Inc(), 3) runs); only a Sequential inside a cycle is rejected.

Workstream A · Stateful pipelines as cycle steps¶

Decision. A Sequential that holds at least one stateful step is a valid step_op for every cycle (Table 4).

1. Accept it. In _coerce_step, return a Sequential unchanged when step_op._is_stateful is true, so Cycle, EnsembleCycle, WindowedCycle and Recurrence take it.

2. Name the rule in the error. A stateless Sequential gets a message that says which wrapper to use.

3. Make the private check public if the cycles are not the only callers: an is_stateful property on Sequential in pipekit core.

1
2
3
4
5
6
7
# Shapes: x (N,)
step = pk.Sequential([pk.Lambda(lambda x: M @ x), Count()])   # stateful   today
# Today:
pc.Cycle(step, n_steps=3)                     # → TypeError                today
# After:
x3, s = pc.Cycle(step, n_steps=3)(x, pc.DAState())            # proposed (A)
# → x3 = M³ x, s.cycle_count = 3

Workstream B · Observation inputs on the cycle’s clock¶

Decision. The DA cycles take their observation index, and optionally R, from the cycle’s own clock in the state (Table 5).

B1 · An index that survives restarts¶

The cycle already advances cycle_count in the state. With obs_clock="state" it passes that count to obs_source instead of the call-local i. Nothing else changes, and the default stays "call" so existing scripts keep their meaning.

B2 · A per-step observation error¶

Today a script can return (y_t, R_t) from obs_source and unpack it in its own analysis step, because the analysis is user code. That works, but it hides R from the cycle and from every analysis step that expects obs_err_cov. obs_err_source puts $R_t$ where the protocol already expects it.

1
2
3
4
5
6
7
8
9
10
11
# Shapes: as stage 2
# Today:
da5 = pc.DACycle(KFForward(), pc.LinearObs(H), KalmanAnalysis(),
                 obs_source=obs_source, n_steps=5)                    # today
# → second call reads y₀ … y₄ again; max |m_restart − m| = 15.4 ppb

# After:
da5 = pc.DACycle(KFForward(), pc.LinearObs(H), KalmanAnalysis(),
                 obs_source=obs_source, obs_err_source=r_source,
                 obs_clock="state", n_steps=5)                  # proposed (B1, B2)
# → second call reads y₅ … y₉; m_restart equals stage 2's m

Testing gates¶

Each workstream merges only when its gates pass; the existing gates are checks the code passes today and must keep passing (Table 6).

Sequencing and risks¶

The two workstreams are independent and small, so they ship together with Report 13’s ensemble histories; the documentation follows (Figure 3).

• A default that changes meaning. Making obs_clock="state" the default would change every existing script that restarts a cycle on purpose with per-call observations. Fallback: keep "call" as the default and warn when a DA cycle starts with cycle_count > 0 and the default clock.

• Stateful detection leaks a private attribute. Workstream A reads Sequential._is_stateful. Fallback: add the public is_stateful property in pipekit core in the same release, or duck-type on _apply_stateful.

• Analysis steps that need the state. Some analyses want the cycle count or the time (adaptive inflation, time-dependent localisation), and today the analysis never sees the state. Fallback: keep the carrier convention and pass time through obs_source’s output; revisit only if two adapters need it (open question 2).

• Unverified adapters. filterax, vardax and plumax are not pinned in references/, so whether their adapters satisfy the protocols is unchecked. Fallback: pin one of them and add an adapter test that runs stage 2 with its analysis step.

Open questions¶

1. Should a plain Cycle advance the clock? DACycle advances DAState.t by dt; Cycle leaves it at 0, so Report 13 computes valid times from the history index. Advancing t when the step is a ForwardModel and the state has t would make the two agree.

2. Should AnalysisStep receive the state? Passing state= as a keyword would let analyses read the clock, but it changes the protocol every adapter satisfies.

3. A cycle checkpoint. The scoping report’s Checkpoint(cycle_op) does not exist. A checkpoint is the triple (config, CarryState.to_dict(), carrier); is it a pipekit-cycle operator, a pipekit-experiment artifact, or a script’s job?

4. The window’s start state in SmootherCycle. Should the analysis receive the state the window started from, so a 4D-Var analysis need not carry it inside the carrier?

5. CompositeForward and dt. Summing the components’ steps is one reading of composition (operator splitting with sequential sub-steps). Should a composite that advances all components over one shared dt exist as well?

6. The latent design. Build latent.py as docs/design/latent.md drafts it, or retire the draft in favour of the shipped LatentMap and ReducedOrderModel protocols?

Upcoming work¶

Table 7 lists the follow-ups.

References¶

Sources¶

• The scoping report, “Report 10 — pipekit-cycle: time-stepping, DA, and observation operators” (supplied by the author).

• pipekit-cycle 0.0.2 in references/geostack/pipekit/packages/pipekit-cycle: its source, README, CHANGELOG, docs/design/latent.md and tests.

• The running example’s script and outputs, described in Section ??.