pypfda

Particle filter data assimilation in pure Python, with first-class support for paleoclimate Observing System Simulation Experiments (OSSEs).

Status — work in progress. pypfda is the open-source companion to a paper currently in preparation (Fallah et al., 2026). The core particle filter, weight / ESS / resampling primitives, observation tempering, post-resample inflation, max-weight degeneracy cap and Gaspari–Cohn localization are in place and tested; higher-level diagnostics (genealogy tracking, Welch / Nyquist spectral tools) and the paleo forward-model subpackage are on the roadmap. Public APIs may evolve before v1.0. Pin a specific version in production code.

pypfda provides a clean, model-agnostic implementation of the sequential importance resampling (SIR) particle filter, plus the building blocks needed to deploy it on real Earth-system problems: spatial localization, ensemble inflation, degeneracy prevention, multi-year assimilation windows, and ensemble diagnostics (effective sample size, weight entropy, genealogy, spectral analysis).

The core filter knows nothing about climate models. You bring a forward model — a Lorenz-96 toy, a coupled GCM, anything callable — and pypfda runs the analysis cycle.

Architecture

pypfda is a model-agnostic engine. A dynamical core plugs in through one small interface — forecast · observe · get/set_state · inflate — and the same sequential-importance-resampling cycle (cost → weights → ESS check → systematic resample → inflate) drives it. The companion OSSE study runs this identical filter on two structurally independent ocean cores — CM2Mc-BLING (MOM5 z-level GCM) and CLIMBER-X (GOLDSTEIN frictional-geostrophic EMIC) — and recovers AMOC variability from sea-surface-temperature pseudo-observations in both, with the same diversity–memory trade-off emerging in each.

It works on Lorenz-96

A 40-variable Lorenz-96 twin experiment with 400 members, observations on every other variable, observation-error tempering (η = 4), post-resample Gaussian inflation, and a max-weight cap (0.3). The particle filter drives ensemble-mean RMSE down to roughly half of the free-ensemble RMSE and keeps the effective sample size well above the degeneracy threshold. Reproduce with

python examples/03_lorenz96_twin.py

This is the same combination of techniques the companion paper applies in an online coupled-climate-model OSSE; Lorenz-96 is the smallest chaotic benchmark on which those techniques can be exhibited end-to-end without a climate model in the loop.

Paper claims reproduced on Lorenz-96

Four of the paper's central methodological findings re-expressed on Lorenz-96. Each is a self-contained script that runs on a laptop in seconds to a few minutes and writes a single figure to docs/_static/. These are analogues, not reproductions — Lorenz-96 is not AMOC — but they show that the filter's qualitative behaviour is not an artefact of the coupled model.

1. Diversity–memory trade-off (§5, headline claim)

A post-resampling Gaussian inflation sweep, repeated across 5 independent rng seeds so that the U-curve is demonstrably not a single-seed artefact. Zero inflation collapses the ensemble to duplicates and the filter dies (RMSE ≈ 4.75 ± 0.09); too-much inflation washes out the accumulated information (RMSE rises back to ≈ 4.12 ± 0.21 at σ = 2). A clean U-curve with optimum near σ ≈ 0.35, RMSE ≈ 1.66 ± 0.41. The across-seed envelope (shaded blue) narrows at the endpoints (robust collapse and robust over-noise) and widens at the optimum (the filter is on the edge of its stable regime, where different seeds land in different minima of the same U). The right panel re-parameterises the seed-mean sweep as skill vs achieved diversity. Reproduce with python examples/05_l96_diversity_memory.py.

2. Cycle-length sensitivity — the two intrinsic bounds behind T10 / T11 / T12

The paper's cycle-length finding for AMOC is that two intrinsic timescales constrain the usable observation interval:

• the Lyapunov time (how fast forecast uncertainty grows without observations), and
• the Nyquist bound (half the period of the dominant spectral mode).

The tighter of the two determines where the filter breaks down. On AMOC, the coupled-system Lyapunov time is decades while Nyquist is ~6.7 years (from the 13.3-year AMOC mode), so Nyquist binds and a 5-year cycle (T11) is optimum.

This figure shows the same two bounds on Lorenz-96, correctly marked. On L96 the ordering is reversed: Lyapunov is ~0.42 and the Nyquist bound (from the spectral peak in figure 4 below) is ~0.83, so Lyapunov binds and the filter breaks down at ~0.4, well before any aliasing kicks in. We are not claiming L96 reproduces the AMOC U-curve — it doesn't, and it cannot, because L96's chaos horizon is short relative to its oscillation period. We are claiming both bounds exist here, they are measurable, they have the values the theory predicts, and which one matters is set by the system, not the algorithm. That is the point the paper makes, made visible on a benchmark that fits on a laptop. Reproduce with python examples/04_l96_cycle_sensitivity.py.

3. Diverse-IC ensemble (analogue of T13)

Each of 200 members is independently spun up on the attractor, so the initial ensemble is maximally diverse. Left: the first state variable over time. DA (blue) tracks the truth (black); the free ensemble mean (grey) collapses to climatology because averaging uncorrelated attractor trajectories returns the attractor mean. Right: ensemble-mean RMSE over time. DA pulls the scattered ensemble onto the truth trajectory without shared initial-state memory, matching the paper's T13 finding. Reproduce with python examples/06_l96_diverse_ics.py.

4. Welch spectrum and the Nyquist argument (§4, Figure 2 analogue)

A 2000-unit control integration of Lorenz-96 shows a genuine spectral peak at T ≈ 1.65 model time units with peak / median power ratio ≈ 134 — sharp, not broadband. The corresponding Nyquist bound (T/2 ≈ 0.83) is the same one overlaid on figure 2 above. The autocorrelation panel independently reports an e-folding time of 0.30, consistent with the Lyapunov-limited predictability horizon. Reproduce with python examples/07_l96_nyquist.py.

The walk-through with derivations and extra figures lives in the docs at Paper claims reproduced on Lorenz-96.

Coupled fast–slow benchmarks

Observe only the fast ring, reconstruct the slow large-scale index. The DA ensemble mean (orange, with its 5–95% band) tracks the truth (black) while the unassimilated free ensemble (blue, dashed) cannot — a clean, laptop-runnable analogue of inferring slow ocean overturning from fast surface fields. Reproduce with python examples/08_two_scale_l96_fast_slow.py.

The package ships two coupled fast–slow systems as first-class ForwardModels in pypfda.models.lorenz, each driven by the same CycleDriver that the companion paper points at its coupled GCM cores:

• TwoScaleLorenz96 — the two-level Lorenz-96 model (Lorenz 1996; Lorenz & Emanuel 1998): slow large-scale variables X_k coupled to a fast small-scale ring Y_{j,k}.
• CoupledLorenz63 — a two-timescale coupled Lorenz-63 (Peña & Kalnay 2004): a fast "atmosphere" two-way coupled to a slow "ocean".

In both, the twin OSSE observes only the fast variable and reconstructs only the slow one — the minimal, laptop-runnable analogue of reconstructing slow ocean overturning (AMOC) from fast surface temperature (SST). Each run is checked by pypfda.verify.scan_osse_result (no clone / stale-diagnostic artifact) and records the forecast diagnostic, so it cannot silently drift to the optimistic analysis convention. On two-scale Lorenz-96 the filter lifts the slow-index correlation from ≈ −0.1 (free) to ≈ +0.2 (Δr ≈ +0.3), at the honest cost of near-total genealogical collapse (effective ancestor size → 1) — the diversity–memory trade-off the inflation kernel exists to manage. The coupled Lorenz-63 case is the easier mirror image: the same operators cut slow-variable RMSE by ≈ 97 % (correlation +0.52 → +0.83) while the genealogy stays healthy (effective ancestor size ≈ 40), spanning both the degenerate and the well-conditioned regime without a single change to the filter. Example 10 pins down a clearly working regime as a reusable positive benchmark: with a chaotic slow subsystem whose unobserved variable genuinely oscillates and switches attractor lobes throughout the run (standard deviation ≈ 8, not a frozen fixed point), strong coupling, and a generously tempered likelihood that keeps the genealogy alive, the filter pulls the diverse free ensemble off the wrong attractor lobe and tracks the time-varying slow variable — flipping its correlation from ≈ −0.5 (free) to ≈ +0.7 (Δr ≈ +1.2; median r_DA ≈ +0.68, with r ≥ 0.5 on every representative seed and r ≥ 0.6 on ≈ 90 % of them) and cutting the anomaly RMSE (level-removed, so this is variability tracking, not a constant offset) from ≈ 8.3 to ≈ 6.1, while the effective sample size stays ≈ 0.85 N and the effective ancestor size stays ≈ 0.46 N and varies cycle to cycle — a genuine, non-degenerate constraint from the observations rather than a single particle memorising the truth.

Example 11 is the textbook positive control: a single canonical Lorenz-63 attractor (σ = 10, ρ = 28, β = 8/3), observing only x (with noise) and reconstructing the unobserved z. With just three state dimensions and a few hundred particles there is no curse of dimensionality, so the same SIR filter reconstructs the hidden dimension almost perfectly: the DA ensemble-mean tracks the truth's chaotic z oscillation while the diverse free ensemble collapses to the attractor centroid (correlation +0.15 → +0.93, Δr ≈ +0.77; anomaly RMSE cut by ≈ 52 %, 8.6 → 4.1), and it does so on every seed (7/7 with r ≥ 0.92). The genealogy stays healthy and exercised — effective sample size ≈ 0.66 N and effective ancestor size varying cycle to cycle (median ≈ 15, minimum ≈ 6, never collapsing to a single lineage) — so the gain is an observational constraint, not one particle memorising the truth.

python examples/08_two_scale_l96_fast_slow.py
python examples/09_coupled_l63_fast_slow.py
python examples/10_coupled_l63_clearly_works.py
python examples/11_classic_l63_partial_obs.py

Scope — what pypfda is and is not

pypfda ships the data-assimilation method, not a turnkey coupling to any specific Earth-system model. In concrete terms:

What is here, today. The particle filter itself (SIR), numerically stable weight and ESS computation, four resampling schemes, Gaspari–Cohn localization, observation-error tempering, a max-weight degeneracy cap, and post-resample inflation. Demonstrated end-to-end on Lorenz-63 and Lorenz-96.

What is not here, today, and will not arrive by magic.

• Integration glue for any specific climate model (CM2Mc-BLING, MITgcm, CESM, ICON, …). Nothing in this repository knows how to write an FV3 restart file, launch a SLURM job array, or cycle ensemble state between INPUT/ and RESTART/ directories.
• The HPC orchestration used by the companion paper — that is a separate, cluster-specific driver that calls pypfda functions but lives elsewhere.
• A coral δ¹⁸O proxy forward model, a PAGES 2k loader, or any paleo-specific I/O. These are planned for pypfda.paleo but are not yet implemented.
• The 100-member ensemble output underlying the paper's figures (that sits on the authors' cluster; see Paper data below).

Using pypfda with your own model is real engineering work. You are expected to write the driver that runs your forward model, reads its state into a NumPy array, calls pf.assimilate(...), writes the resampled state back, and handles your own HPC scheduling. For Lorenz-style ODEs that driver is a few dozen lines (see examples/). For a coupled climate model with distributed restart files and a queue system, expect weeks of integration effort per model. If you are attempting this and want guidance, open an issue or email the authors — we are happy to help, but there is no plug-and-play path.

Highlights

• Pure Python. No Fortran, no compilation. Works on Linux, macOS, and Windows.
• Model-agnostic. The filter is decoupled from any specific simulator; bring your own forward step.
• Numerically careful. Log-domain weight computation, well-tested resampling routines (systematic, stratified, residual, multinomial), numerically stable ESS.
• Diagnostics. ESS, weight entropy, genealogical diversity, rank histograms, spectral / Nyquist analysis for cycle-length design.
• Paleoclimate-ready. Optional paleo extra includes coral δ¹⁸O proxy system models and a PAGES 2k loader.
• Production tooling. Strict typing (mypy --strict), property-based tests, ≥ 80 % coverage, conventional commits, semantic versioning, reproducible builds via hatchling.

Installation

pip install pypfda                 # core only
pip install 'pypfda[io,plot]'      # + NetCDF and matplotlib helpers
pip install 'pypfda[paleo]'        # + coral PSM and PAGES 2k loader
pip install 'pypfda[all]'          # everything including dev + docs

pypfda requires Python ≥ 3.10.

60-second example

import numpy as np
from pypfda import ParticleFilter

rng = np.random.default_rng(0)
n_members, n_obs = 100, 5

# Linear Gaussian toy: x_{t+1} = 0.95 x_t + w,   y_t = H x_t + v
def forecast(state):
    return 0.95 * state + rng.normal(0, 0.5, state.shape)

H = rng.normal(size=(n_obs, 4))
truth = rng.normal(size=4)
ensemble = rng.normal(size=(n_members, 4))

pf = ParticleFilter(ess_threshold=0.5, resampling="systematic")

for _ in range(50):
    truth = 0.95 * truth + rng.normal(0, 0.5, 4)
    ensemble = np.array([forecast(m) for m in ensemble])
    obs = H @ truth + rng.normal(0, 0.1, n_obs)
    obs_pred = ensemble @ H.T
    ensemble, info = pf.assimilate(ensemble, obs_pred, obs, obs_err=0.1)
    print(f"ESS={info.ess:.1f}  resampled={info.resampled}")

See the quickstart and Lorenz-63 tutorial for full walk-throughs.

Documentation

• Theory — derivation of the SIR update, comparison of resampling schemes, the diversity/memory trade-off, choosing inflation and localization parameters.
• Tutorials — Lorenz-63 twin experiment and, for high-dimensional chaos, the Lorenz-96 demo used by the figure above.
• API reference — every public function and class, generated by Sphinx + autosummary.

Read the docs at https://bijanf.github.io/pypfda.

Citing pypfda

If you use pypfda in published work, please cite the software (via the CITATION.cff button on GitHub) and the methodological paper:

Fallah, B., Rostami, M., Huiskamp, W., Goosse, H., & Rahmstorf, S. (2026). Online Particle Filter Data Assimilation for AMOC Reconstruction: A Diversity–Memory Trade-off in Observing System Simulation Experiments. In preparation.

A BibTeX snippet is provided in CITATION.cff.

Paper data

The companion paper is an Observing System Simulation Experiment built on the coupled CM2Mc-BLING climate model with 100-member ensembles integrated for ~100 years each (several months of cluster wall time per experiment, of order one terabyte of netCDF output). That raw ensemble is not distributed with this repository and would be impractical to re-generate from scratch. It currently resides on the Potsdam Institute for Climate Impact Research (PIK) cluster; interested researchers are welcome to contact the authors for access or for processed diagnostics.

What is in this repository is the method: a model-agnostic implementation of the techniques the paper applies (SIR, observation tempering, post-resample inflation, max-weight cap, Gaspari–Cohn localization), validated on Lorenz-96 (see the figure above). That is what the paper's Code Availability statement points to.

Scripts that operate on small processed diagnostics — for example the Welch / Nyquist spectral analysis of the control AMOC time series — may be added under reproduce/ in a later release, driven by a small Zenodo deposit. Regeneration of the full figure set is not a goal of this package.

Related work

• DA_offline_PF — the Fortran offline particle filter from Dalaiden et al. that motivated this Python implementation.
• Dubinkina, S. et al. (2011), Testing a particle filter to reconstruct climate changes over the past centuries, IJBC 21, 3611.
• Goosse, H. et al. (2010), Reconstructing surface temperature changes over the past 600 years using climate model simulations with data assimilation, JGR 115.

Contributing

Bug reports, feature requests, and pull requests are welcome. Please read CONTRIBUTING.md and the Code of Conduct before opening an issue or PR.

License

pypfda is distributed under the MIT License.