Conversation
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>
Prior realisations, perturbed observations (AssimilationEnsemble, ES-MDA, EnKF, multilevel), outlier and crash replacement, the auto-adaptive localization's shuffle and popt's control perturbations all drew from NumPy's global functions, so a run could only be reproduced by seeding the whole process, and any other code drawing in between changed the result. The base ensemble now owns `self.rng`, built from `keys_en['seed']`: a private RandomState when a seed is given, otherwise `GlobalRandomStream`, a picklable stand-in for the global functions (the numpy.random module itself cannot be pickled, and the emergency dump pickles the ensemble). Every draw site takes it; the localization factory passes it to builders. The geostat sampler is replicated draw for draw in `misc.sampling.gen_real` with the stream as an argument (tests compare it to geostat under identical seeds for vector, diagonal, full and scalar variances, with and without limits). geostat stays a dependency for its covariance builder. Verification: ruff clean; tests/test_sampling.py and tests/assimilation/test_seed_option.py (seeded run reproduces under different global states, leaves the global state untouched; unseeded run still governed by np.random.seed); full suite 430 passed with all thirteen characterisation goldens unchanged. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…diction `restart_sim_results.pkl` is a hand-placed file: a saved forecast copied to that name so a restarted run can skip the forecast it had already finished. Nothing in PET writes it, yet the forecast consumed it on any run that found it in the working directory, so a forgotten file silently replaced a fresh forecast. It is now used only when the ensemble's `restart` flag is on, and once used it is moved into the results folder as `sim_results.pkl` (the working directory when saving is disabled, as before). `calc_prediction(..., save_prediction=name)`, which only popt calls, read `self.ensemble.keys_da` -- an attribute the base ensemble never had -- so the option raised AttributeError, and it wrote into a folder it never created. The folder now comes from the ensemble's own options (`savefolder` or `save_folder`, default `Predictions`) and is created first. Verification: ruff clean; tests/assimilation/test_restart_forecast_file.py and test_save_prediction.py fail 6 of 7 on the parent commit in a detached worktree (the passing one is the nosave restart, unchanged by design) and pass here; full suite 437 passed with the characterisation goldens unchanged. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
`minimize`, the constructor tail (restore-or-evaluate, banner, first log row, first result), `_accept_step`'s bookkeeping, the callback, nine `if self.saveit:` sites, four `_log_iteration` bodies and three identical projected-gradient checks were copied across EnOpt, LineSearch, TrustRegion and SmcOpt. `update_step` returned a bool and had to mutate six attributes listed only in a comment; forgetting one gave a run that iterated and logged normally while never converging. The base now defines `StepReport(accepted, message)` and `_commit_step(x, f, jac=, hess=)`; `update_step` commits through the latter and returns the former. `run_optimization` evaluates the starting point (`_start`), and after every accepted step runs the callback, `_record_results`, the log row built from `log_columns()`, and the function/state/gradient checks. The four optimizers keep their step, their state updates and their columns. `EnOpt`/`SmcOpt` constructors take `(x0, fun, ...)` like the other two and every `minimize`; SmcOpt's `autorun` is gone. Steihaug's unconditional prints are DEBUG log records. Also fixed: `LineSearch(recompute_jac=n)` cleared the gradient on retry and then took `-None` as the next direction. Verification: 21 deterministic cases (LS GD/BFGS/Newton-CG, TR iterative/CG-Steihaug with exact and BFGS Hessians, EnOpt GD/Adam/AdaMax/ Steihaug plus hessian+nesterov+resample, seeded ensemble LS and SmcOpt, each with/without the unit-cube transform) give bit-identical x, fun, nit, nfev, njev, nhev and message before and after; the recompute_jac reproduction raises TypeError on the parent commit and converges here; ruff clean; full suite 443 passed. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…eads Two restart mechanisms coexisted and neither worked as a restart. The scheme's RestartMixin never saw the config: every scheme passed only zero tolerances to its base, so `restart`/`restartsave` were always False, `save_restart()` was unreachable, every `if self.restart is False:` guard was dead-true, and ES-MDA's restart branch referenced an undefined `loop_ind`. The ensemble read `restart` from the merged config and loaded a pickle of its own `__dict__` (`emergency_dump`, also written by `restartsave`), which restored the state but left the scheme re-initialised from scratch: iteration counter, damping, misfit history all reset. Now `restart_options(keys_da)` carries the three keys from `[dataassim]` into the base; construction always initialises; `run_assimilation()` overlays the checkpoint. The payload is the loop's bookkeeping plus `RESTART_ATTRIBUTES` declared per scheme (perturbed observations, ES-MDA's un-inflated draw, the subspace E and W, the misfit the acceptance test uses, lam/gamma) plus the ensemble's `restart_state()` (state, prior, forecast, scaling, SVD of the scaled prior, iteration, and its random stream). The ensemble no longer loads anything; `save()` remains the crash dump. `PETStateArray` gained `__reduce__`/`__setstate__`: an ndarray subclass loses its attributes on unpickle, so `indices` was absent on any state array read back from either file. Verification: tests/assimilation/test_restart_resume.py interrupts ES-MDA, LM-EnRML (approx) and GN-EnRML (subspace) after their second accepted iteration, resumes from the checkpoint in a process seeded differently, and gets bit-identical x and data misfit to an uninterrupted run; checkpoint written after the prior forecast; `restart = yes` no longer demands an emergency_dump; ruff clean; full suite passed with the characterisation goldens unchanged. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
The forecast was one 180-line method at five nesting levels: member input construction, three execution backends, adjoint splitting, output coercion for two simulator return shapes, scaling and saving, all inline. It now reads as the sequence it is -- `_simulator_input`, `_run_members`, `_replace_failed_simulations`, `_collect_adjoints`, `_collect_sim_data` -- with the same operations in the same order. A backend, a return shape or a forecast mode is now a change to one step rather than to the loop. Verification: ruff clean; the thirteen characterisation goldens are unchanged (serial backend); tests/assimilation/test_forecast_backends.py pins the process-pool backend against the serial one for the first time; full suite passed. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Four reader functions and their two string-to-array helpers, 470 of the file's 709 lines, had no callers in src, tests or the tutorials; DataReader, the reader the ensemble uses, called none of them. The module docstring described the dead functions as the module's main API and is rewritten around DataReader. Verification: ruff clean (the pandas import went with them); reader tests and full suite passed. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
The order of every data array -- observations, variance, predictions, adjoints -- was implied by PETDataFrame.to_matrix(): label-major, then data type, all-NaN cells dropped. Observations and predictions were flattened separately, so a single unobserved cell left the observation vector one row shorter than the prediction matrix with nothing raised (probe: 14 vs 15 rows on the tiny case). DataLayout is that order computed once from the observed frame: rows of (label, datatype, start, stop). The pipt ensemble builds it after scaling and exposes obs_vector; ESMDA, EnKF, the iterative schemes and multilevel read that instead of flattening data_df. to_frame() gives the frame view back from a vector or an (nd, ne) matrix with the cell shapes the merged prediction frames use, so the legacy flatten of a view reproduces the matrix. Predictions still travel as frames; they move in the next step. Verification: ruff clean; tests/test_data_layout.py (walk and skipping, vector == to_matrix with and without an empty cell, matrix <-> view round trips), reader test that obs_vector == data_df.to_matrix() on the tiny case; full suite passed with the goldens unchanged. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
pred_data is now a PredictedData: the (nd, ne) matrix in the layout's row order, plus the layout. The pipt forecast fills it from what each member's simulation returned -- records placed by the simulator's true_order, or a frame per member -- scaled as the observations were, instead of merging the members into an object-dtype frame, filtering it by label, and flattening it on every analysis attempt. The schemes read pred_data.matrix; the frame is a view (to_frame) for QA/QC and inspection. Observations and predictions share one row order by construction: a member that lacks an observed cell or returns the wrong length is reported, where the frame path silently produced a shorter observation vector. Moved onto the matrix: the multilevel model-error correction (per-row means), outlier detection (get_outlier_index takes arrays), and the `scale` option, which had never worked -- it iterated the characters of the column names. The seismic compression path still runs on the frame, built from sim_data as before, and is wrapped into the container afterwards; it moves once it has a test. sim_data, the full forecast, stays a frame and stays what is saved. Verification: ruff clean; tests/test_predicted_data.py (records and frames fill identically, missing type and wrong size are reported, scaling equals the frame's, view and member selection); tests/assimilation/ test_prediction_fill.py (matrix == legacy flatten on the tiny case with and without scaling; a blanked observation no longer misaligns); outlier, multilevel, QA/QC, pipeline, restart-file, backend and scheme-base tests adapted or unchanged; full suite passed with the thirteen characterisation goldens untouched. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
The full forecast was merged into an object-dtype frame on every forecast, whether or not anything read it; on the analysis path nothing does since pred_data is filled directly. The members' raw outputs are kept (member_outputs) and the frame view (sim_data) is built on first use -- saving, QA/QC, popt's objective functions -- and cached until the next forecast. Outlier replacement reorders the raw outputs instead of mapping a function over every cell of two frames; a forecast loaded from a file, which exists only as a frame, is still mapped. Adjoints reach the analyses as an (nd, nx, ne) array in layout order, filled from the members' adjoint frames and scaled with the data, instead of a merged frame flattened on every analysis attempt. The frame path stacked every row the simulator reported rather than the observed ones. Values are identical to the legacy stack; that stack was non-contiguous, so its member mean summed in a different order, and adjoint-based updates move at the 1e-13 level (Van der Pol ES-MDA: max |dx| 1.9e-13). No golden covers adjoint runs. Verification: ruff clean; tests/assimilation/test_adjoints.py (array == legacy stack, None without adjoints), test_remove_outliers on raw outputs and the array (plus the frame fallback), failed-member, save_prediction, prediction-fill, restart and popt ensemble tests; full suite passed with the characterisation goldens untouched. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…_cov The data covariance was assembled by walking the variance frame cell by cell and then dropping every NaN entry, so a NaN variance for an observed cell left the covariance one entry shorter than the observation vector with nothing raised -- the same silent misalignment the predictions had. EnKF, the multilevel ensemble and the observation perturbation each called the walk themselves. The ensemble now builds obs_variance once, in layout order: an (nd,) vector, or the (nd, ne) error ensemble when the variance is empirical. A NaN for an observed cell is reported with the cell. The three call sites and outlier detection read it; construct_data_cov and the outlier helper that rebuilt the same array are gone. Verification: ruff clean; reader test that obs_variance equals the frame flatten on the tiny case and that a NaN variance is reported; outlier, multilevel, step-and-score and prediction-fill tests; full suite passed with the characterisation goldens untouched. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
The state ensemble was an ndarray subclass carrying the {variable: (start,
stop)} row map. __array_finalize__ copied that map onto every slice and
view, so x[3:5].to_dict() reported the full layout on a five-row array;
unpickling dropped it (patched with reduce hooks in fff9583); twenty
operator overrides existed only so a type checker inferred the subclass.
The state is now a plain (nx, ne) array. The row map lives once, as the
ensemble's idX dictionary -- which popt already maintained itself -- and
misc.structures.StateLayout wraps it with the conversions the boundary
needs: to_dict(enX), member_dicts(enX) for the simulator, clip(enX,
limits), and the constructors from_dict and from_prior_info, both returning
(matrix, layout). BaseEnsemble.state_layout derives the object from idX, so
there is one source of truth for pipt and popt alike. The analyses were
already doing their arithmetic on arrays; the schemes, the scheme base's
saving and QA/QC hand-off, and the multilevel scheme go through the layout.
Verification: ruff clean; tests/test_structures.py covers the layout
(rows, dict views, member dicts round trip, from_dict with ne, the three
clip forms, prior generation with consecutive indices); failed-member,
save_prediction, multilevel, QA/QC, scheme-base, step-and-score and popt
ensemble tests; full suite passed with the thirteen characterisation
goldens untouched. The state tutorial cell and dev guide describe the
layout.
Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
The scheme-writing tutorial's example read pred_data.to_matrix(), enX.indices and construct_data_cov, none of which exist after the data-structure work; it now reads pred_data.matrix, idX, obs_vector and obs_variance. The frame tutorial's table and paragraph describe how the data reach the analyses: one DataLayout, matrices built in its order, the frame as the view. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
The state's row map and the data vector's row order are the same kind of thing -- the one fixed ordering everything else is built in -- so the two layouts live in one module, with the prior-limits helper StateLayout uses. No behaviour change; imports updated. Verification: ruff clean; structures, layout and predicted-data tests; full suite passed. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Compressed data types (the `compress` section) were the last frame on the analysis path: post_process_forecast rebuilt the prediction frame from sim_data, deep-copied the compressed columns, compressed every cell member by member, wrote the coefficients back into frame cells, and was wrapped into the container afterwards. Compression only happened when post_process_forecast was also enabled, so a config with `compress` alone left raw predictions against compressed observations. Compression is now a per-row transform inside PredictedData.from_members: a member's raw vintage goes through the same SparseRepresentation the reader built for the observed vintage -- the n-th compressed layout row is vintage n, both walk the frame the same way -- and enters the matrix as its leading coefficients. It follows from `compress` alone; post_process_forecast only enables the sim2seis scaling. Reconstructions are computed only when saveforecast writes them. CompressionMixin and compress_manager are gone. `use_ensemble` is refused with the reason. For it the reader kept the observed vintage raw while giving it the compressed-length variance, and schemes perturb observations at construction, so the option failed on the shape mismatch before any forecast; it is listed under Known issues next to screendata, which needs the same change. Verification: ruff clean; tests/assimilation/test_compression.py on a synthetic masked-grid case shaped like the AVO one (two vintages, db2 level 2, universal hard thresholding): observed rows equal the leading coefficient counts with est_noise^2 as variance, each member's filled row equals its raw vintage compressed by hand, an uncompressed type is untouched, ES-MDA runs and reduces the misfit, use_ensemble is refused, reconstructions are written only with saveforecast; full suite 471 passed with the characterisation goldens untouched. On the real AVO case (293x60x60 grid, two vintages): the reader reproduces the previous run's saved compressed observations and variances bit for bit (7376 and 7122 coefficients), and four real member vintages filled through the new path in 0.3 s equal their direct compression, aligned with the 14498-row observation vector. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Every scheme set max_iter and then maxiter = max_iter - 1, and nothing read max_iter afterwards. The pair dates from the legacy loop, which counted the prior forecast as iteration 0, so a config's MAX_ITER of 5 meant four updates; the base loop counts updates, and the subtraction kept existing configs running unchanged. The subtraction stays where it means something (the iterative schemes, from the config key); ES-MDA and EnKF take one update per assimilation step, ES one, without the +1 -1 detour. The config key max_iter and the number of updates are unchanged. Verification: ruff clean; full suite passed with the characterisation goldens untouched, which pin the number of updates for every scheme. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
The legacy loop numbered the prior forecast iteration 0, so a config's max_iter of 5 meant four updates, and the iterative schemes subtracted one to keep that meaning under the base loop, which counts updates. The setting now means what it says: max_iter updates. Existing configs get one more update than before; lowering max_iter by one restores the old run. The run table, the convergence message and the assimilation_result_i files already numbered updates from 1 with the prior as 0. Verification: ruff clean; tests/assimilation/test_max_iter.py (LM-EnRML and GN-EnRML take exactly max_iter updates when no tolerance stops them, for 1 and 3; ES-MDA takes one per assimilation step); the characterisation config's max_iter lowered from 3 to 2 and the linear-model test's from 5 to 4, with every golden and pinned number bit-identical, which pins that only the meaning changed; full suite passed. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Config sections were plain dicts that every consumer read with its own fallbacks (savefolder/save_folder in three places, truedata/data, datavar/var, importstaticvar/importstate, restart_file/restartfile), converted yes/no strings itself (eight is_enabled sites) or read them by truthiness (scale_data = "no" enabled scaling), turned legacy row blocks into dicts in place (extract_maxiter rewrote keys['iteration'], organize_sparse_representation rewrote the caller's compress block), and had written back into: the ensemble stored datatype, truedataindex and assimindex in the caller's dictionary. Validation was assert statements run only by the legacy text reader and `pet validate`; a TOML user with a missing key got a KeyError inside the run. input_output.config is the boundary. normalize() gives the three sections in the one form PET reads -- canonical names, boolean flags, dict blocks, field conversions -- as copies; every reader returns it, and both ensemble constructors pass their sections through it, so a script-built dict gets the same treatment as a file. Consumers read one name. validate() reports problems by section and key; `pet validate` prints them all plus keys nothing reads, and AssimilationEnsemble raises ConfigError listing the fatal ones. The legacy reader returns three sections like the others. is_enabled and list_to_dict remain as names for the boundary's helpers. Verification: ruff clean; tests/test_config_boundary.py (aliases, flags, row blocks, caller untouched, idempotent, canonical spelling wins, validation messages and the fatal subset, unknown keys, all three readers agree), tests/assimilation/test_config_boundary_ensemble.py (the ensemble's copy carries the canonical keys; nothing is written back); CLI, parser, reader, migrate and pipeline tests; full suite passed with the characterisation goldens untouched. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Two pages the documentation lacked. `docs/configuration.md` lists every key the code reads from the dataassim, ensemble, optim and simulator sections -- meaning, default, which schemes use it, the sub-blocks (iteration, mda, localization, compress, prior_<name>, controls), the legacy spellings the boundary accepts, and what `pet validate` checks. `docs/architecture.md` describes the three layers, how a run is put together, the scheme, analysis and optimizer contracts, the data layouts on the analysis path, forecast, restart, random numbers and logging, and a table of where a new scheme, analysis, localization, optimizer, simulator or config key goes. Both are in the site navigation and linked from the README and the developer guide; the site builds. `popt` exported nothing; it now exports the four optimizers, the base and StepReport, and the two ensembles. Every public function and class has a docstring (116 were missing, most of them in the generalized ensemble and its marginals). `add_synthetic_noise`, which nothing reads, is no longer a known key, so `pet validate` points it out. Verification: ruff clean; mkdocs build passes with both pages rendered; full suite passed. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…ommits Upstream moved little since the merge base 8cd86cc -- 15 commits over 12 files -- while this line restructured most of pipt and popt. Eleven paths conflicted: six textually and five as modify/delete, because upstream edited files this branch had moved or removed. Every conflict is resolved in favour of this branch. The five modify/delete paths are kept deleted; their content lives on at the new locations: update_methods_ns/{margIS,subspace}_update.py -> update_schemes/analysis/ popt/loop/{ensemble_base,optimize}.py -> popt/{ensembles,optimization_methods}/ popt/update_schemes/linesearch.py -> popt/optimization_methods/ One path needed care that git could not give it. subspace2_update.py is new upstream and absent here, so it merged clean -- into the update_methods_ns package this branch had deleted, importing a mixin contract no scheme here satisfies. It is removed in this commit and ported onto AnalysisBase separately. The only content this merge brings in is upstream's two new README sections, Simulation wrappers and Visualization. The four other upstream changes that are not already covered here -- the UTF-8 log file, the HPC extraction guard, the EPF convergence criterion and subspace2 -- are ported as their own commits rather than resolved into this one. Suite unchanged at 482 passed; ruff clean. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The timestamp format embeds U+2502 and the tables are drawn with box characters. The file handler opened in the OS default encoding, so on a cp1252 Windows every record raised UnicodeEncodeError inside the handler: the log file stayed empty and the run carried on without saying why. Ports upstream ab0293c onto the per-file logger this branch introduced in 2a1b976, which constructs the same FileHandler and so takes the argument unchanged. The stream handler is left alone; upstream did not touch it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
A member whose results could not be extracted took the whole forecast batch down with it: extract_data raised out of run_on_HPC and nothing after it ran. Ports the intent of upstream c629c0f, not its text. Upstream appends to en_pred inside the try and again in the except, so a raise from store_ensemble_sim_information -- which it also moved inside the try -- appends twice for one member. en_pred is positional, so every later member then reads one slot too early: for a four-member ensemble with member 1 failing to save, upstream returns five entries and members 2 and 3 silently take each other's predictions. Here extraction and saving are guarded separately, there is exactly one append per member on every path, and a failure to store diagnostic information is logged without discarding the prediction that was extracted successfully. The message goes to the logger only; upstream also prints, and this module has no other print. tests/assimilation/test_hpc_extraction.py pins all three paths; the second test fails against upstream's version with `assert 5 == 4`. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
… size The loop stopped once no control had moved more than `conv_crit` relative to its previous value. That asks the wrong question twice over: it reported success whenever the inner optimizer stalled, however badly the constraints were still violated, and it refused to finish while one control kept jittering at a feasible point. It now stops once `mean(epf['penalty']) / epf['r']` falls below `conv_crit` -- the penalty with the factor divided back out, so the test reads the violation itself. The objective is handed the epf dict by reference (`_wrap_callable` sets `kwargs['epf']`) and is responsible for writing `penalty` into it; an objective that does not now raises KeyError rather than quietly converging on the step size. The `1e-5` default is kept, so a config written for the old criterion still loads, but the number now carries the units of the objective instead of being dimensionless. CHANGELOG records that under Breaking changes and docs/configuration.md states the contract. Ports upstream 5358e07 and 4b8d878, and the optimize.py hunk of 97fa87a (the empty-penalty guard, and `np.mean(p)/r` rather than `np.mean(p/r)`). `np.asarray` is used before `.size` so a scalar penalty reports the empty case instead of AttributeError. The two other upstream commits against that file, 77ecbed (outer iteration count) and 1d039bf (`ftol` rather than the never-assigned `obj_func_tol`), are already fixed here; the end-to-end test pins the first by asserting exactly `max_epf_iter` outer passes. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Localization was selected by which keyword appeared in the block -- `autoadaloc`, `localanalysis`, `dist_loc` as a key or a bare value, a pickled mask file, or none of them for the parallel update. The rewrite replaced that with a required `name`, so every config written before it stopped at startup with "Localization config has no 'name'", naming three strings its author had never seen. `pet migrate` does not cover the block, so there was no way forward but to hand-edit. `infer_name` works the mode out from the keyword that used to select it, and `normalize_parsed_info` fills it in. An explicit `name` still wins. Two of the five modes cannot run here. `localanalysis` and the parallel update both update each parameter against its own subset of the data, which needs the per-subset observation machinery -- `_ext_obs`, `current_state`, `pert_preddata` -- that the scheme rewrite replaced with one DataLayout built at setup. They are refused as the config is read, with a message that says why and names `autoadaloc` and `distance_loc` as the alternatives. `localanalysis` is no longer registered, so it is not advertised by `available_localizations()` either; it did not previously warn and return None as the changelog claimed, but raised TypeError on construction. `analysis_tools.parallel_upd` is left where it is: the dormant GIES schemes still call it, and they are not being touched. Also reads the auto-adaptive cutoff from wherever the config put it. The value is how many noise standard deviations a correlation must clear -- `nstd` in the old code -- and it was carried as the value of the `autoadaloc` keyword itself. Only `cutoff` was read, so `autoadaloc = 2` ran at the default of 0.3 with no error, a different taper and a different posterior. The default stays 0.3 rather than reverting to the old 1, so only a block giving `autoadaloc` as a valueless flag tapers differently than it used to; the changelog records that. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`_resolve_mask` returned a mask over the whole grid while `_zero_mask` reduced to the active cells, so with an `actnum` the two disagreed: a localized parameter contributed one row per grid cell and an unlocalized one a row per active cell. On a 1x10x10 field with 60 of 100 cells active and two parameters, the operator came out (160, 3) where the state has 120 rows. Nothing checked, so the mismatch surfaced far from here. The kernel is still placed on the full grid -- it has to be, the positions are grid coordinates -- and reduced at the end, on the same C-order flattening the caller applies. An all-active `actnum` now gives exactly what passing none gives, which is the property the second test pins. No test covered `actnum` with distance localization at all; two do now. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Pickled mask files hold plain dicts keyed by (data_type, time, parameter):
`taper_func`, `position`, `range` as [radius, z_range], `anisotropi` as
[ratio, rotation], and `file` for the `import` taper. `_parse_config` loaded such a
file and returned its values untouched, but everything downstream expects
`LocalizationEntry`, so the first attribute lookup raised
AttributeError: 'dict' object has no attribute 'taper'
and no pickled mask file could be used at all.
`_entry_from_legacy` converts one entry, mapping the fields the way `_parse_rows`
maps the equivalent row: range[0] is the radius and range[1] the z-range,
anisotropi[0] the ratio and anisotropi[1] the rotation in degrees. A `taper_func` of
None is a skeleton entry and stays empty; `import` keeps its file and no geometry.
Older files that wrote `range` as the radius alone are read as covering every layer.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The one piece of upstream functionality this branch lacked. `subspace2` solves for
the ne x ne transform W directly, starting from W = I, minimising
J(W) = 0.5 (ne-1) ||W - I||_F^2 + 0.5 ||D - g(xbar + Xp W)||^2_{Cd^-1}
and uses the analytic data covariance through `scale_data` rather than the ensemble
representation E E.T that `subspace` uses, so it takes no SVD and consults neither
`energy` nor `iteration.energy`. Registered on ES-MDA, LM-EnRML and GN-EnRML, which
is where upstream registered it; the sequential schemes cannot apply a transform one
datum at a time, and `subspace` already fails there.
No scheme-side plumbing was needed: `propose_state` has reconstructed
mean(prior_enX) + prior_anomalies * sqrt(ne-1) @ W since ad595a8 added it for margis,
and `IterativeEnRML.RESTART_ATTRIBUTES` already carries W and current_W.
Two departures from upstream 6f313d7, both deliberate:
- The transform is initialised at iteration 0, not 1. Schemes here count from 0, so
the reference version never initialises and dies with AttributeError on
`current_W`. The same correction was already needed for margis.
- The whitened observations are recomputed every call rather than cached on the
first. ES-MDA redraws enObs and scale_data at every assimilation step from
alpha[iteration] * cov_data, so a cache built at iteration 0 drives every later
step with the first step's observations whitened by the first step's Cholesky
factor, while sY uses the current one -- the two in different units. Nothing
reports it: the run completes and the misfit still falls. On the characterisation
case the cache gives a posterior misfit of 175.7 against 201.1 without it, and
the lower number is the wrong one. For LM- and GN-EnRML, where enE and scale_data
are fixed across iterations, both forms give bit-identical results (1437.0 and
830.4), so this costs one solve and changes nothing there. `margis` still caches
the same way; it is registered only on GN-EnRML, so it is inert there, but the
trap is the same one and worth removing separately.
With no reference output to check against, the anchor is an identity: subspace2 is
exactly margIS with `Ratio` forced to 1, i.e. the data error scale taken as known
rather than marginalised over an inverse-chi2 prior. tests/assimilation/test_subspace2.py
asserts that bit-for-bit, for diagonal and full-matrix scale_data, at lambda 0 and 3,
and away from the first iteration where W is no longer I.
The three cases are added to the characterisation suite. All 39 pre-existing golden
arrays are bit-identical; the file grows to 48.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
GenOpt went with a693852 during the popt restructuring, which rewrote EnOpt, LineSearch and TrustRegion into optimization_methods/ and dropped GenOpt with nothing in its place; cma.py followed later as code with no remaining caller. Neither removal was recorded, and src/popt/README.md went on listing GenOpt as implemented. Meanwhile GeneralizedEnsemble.mutation_gradient and .mutation_hessian survived with no consumer anywhere in src/ -- they were written for this. GenOpt draws from the generalized ensemble's marginals and moves the sampling distribution along with the controls: the controls from `jac` with backtracking, then `theta` from `jac_mut` at `alpha_theta`, then the correlation from `corr_adapt` -- either a CMA instance, which needs the ensemble, or any callable, whose result is descended along at `alpha_corr`. Written against OptimizerBase rather than transcribed, so it inherits the choreography every other optimizer here gets: `_commit_step`, the callback, the result record, the log row and convergence. Two behavioural fixes over the version that was removed: - `alpha_corr` is read from `alpha_corr`. It read `options['alpha_theta']` under both names, so setting `alpha_corr` did nothing and the correlation always moved at the theta step size. - It no longer runs itself from its own constructor. `run_loop()` as the last statement of `__init__` is what the `autorun` removal took out of SmcOpt; call `run_optimization()` or `GenOpt.minimize(...)`. `mutation_gradient` gains `return_ensembles=True`, returning the Gaussian samples and their objective values alongside the gradient, so the CMA path adapts from the ensemble the gradient was built from instead of drawing and simulating a second one. cma.py is upstream's file with trailing whitespace stripped for the lint gate; its `hansen2006` citation, orphaned in docs/references.md since the removal, resolves again. README also gains TrustRegion, which was missing from the same list. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
`EnsembleOptimizationBase.function` raised RuntimeError whenever `calc_prediction` reported failure. The optimizer proposes control vectors and some of them are ones the simulator cannot run, so a single such trial point ended the optimization and threw away every iteration before it. It reports inf instead. No backtracking comparison can improve on inf, so the point is rejected and the next trial is taken from the last good iterate -- the behaviour 765d637 removed when it cleaned up the popt ensembles. A crashed single-point evaluation additionally leaves `stateF` alone. `gradient` computes `enF - repeat(stateF, nr)`, so writing inf there poisons every later gradient with inf/NaN rather than rejecting the one point; upstream parked the value in `enF` instead, which corrupts the ensemble objective for the next gradient in the same way. Leaving the last good value in place is what the caller actually needs, and the inf still reaches the optimizer as the return value. The separate abort when every member of a forecast fails is left as it is: there is nothing to compute a gradient from, and tests/test_logging_and_paths.py pins it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
rolfjl
approved these changes
Sep 21, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
4 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
This is large — 252 files, +59.6k / −24.0k — and it merges cleanly: every commit currently on main is already in this branch's history, including all fifteen since 8cd86cc. Suite is 559 passing, ruff check src tests clean.
It is two bodies of work. The first ~199 commits are the POPT restructuring done on this fork in July (loop/ → ensembles/ + optimization_methods/, every optimizer on a shared OptimizerBase). The remainder does the same for PIPT: schemes hold an analysis rather than inheriting one, so eighteen hand-written classes became five algorithms × four flavours through a registry.
What is new here that main does not have
A characterisation suite pinning 16 (scheme, analysis) pairs to stored numerics at rtol 1e-9, plus ~550 tests where main has 6.
A typed configuration boundary (input_output/config.py), a pet CLI with migrate/convert, and a 950-line CHANGELOG.
ensemble.protocols.ForwardSimulator — the simulator contract written down.
What I took from main, and what I didn't
subspace2 (#154) is ported to the new analysis package. Two deliberate departures: the transform initialises at iteration 0 (schemes here count from 0, so the original never initialises), and the whitened observations are recomputed rather than cached — ES-MDA redraws them every step with alpha[iteration] * cov_data, so the cache drives later steps with the first step's observations. On our test case that is the difference between a posterior misfit of 175.7 and 201.1, and the cached, lower number is the wrong one.
Our margis fixes three things in #154's version: the hardcoded range(70) loop (an IndexError below 70 data rows, and otherwise one data type per row), an iteration == 1 off-by-one, and a duplicated scale() helper.
I did not take the truncSVD → eps-floored SVD change to subspace. It makes the documented energy key inert and retains every direction down to 1e-8·σ_max. We fixed a different bug on that line — whitening Y, without which the weights depend on the units of the data. tests/assimilation/test_subspace_scale_invariance.py is the evidence. Happy to add the eps-floor as an option.
c629c0f is ported by intent, not by text. It appends to en_pred inside the try and in the except, so a raise from store_ensemble_sim_information appends twice for one member and every later member reads one slot early. A four-member ensemble returns five entries. tests/assimilation/test_hpc_extraction.py fails against the current version with assert 5 == 4.
Removals
cov_regularization.py → the pipt.localization package; data_tools.py → PETDataFrame methods; basic_tools.py → byte-identical duplicates of input_output/get_ecl_key_val.py; loop/assimilation.py → scheme.run_assimilation(). esmda_geo is gone — it took the wrong constructor arguments and read an attribute never set, so it could not be built. GenOpt and CMA were removed during the July restructuring and are restored here.
Known issues, stated plainly
Local analysis and the parallel update are refused at config time rather than half-working — both need per-subset observation machinery this rewrite replaced, and local analysis previously left the posterior equal to the prior while reporting a misfit. es/enkf with subspace raises on the sequential path. The GIES schemes cannot be constructed (pre-existing). Notebooks are updated but not re-executed — they need OPM flow. ES-MDA restart does not carry the subspace transform.
Backwards compatibility for existing config files is covered: a LOCALIZATION block that names no mode still parses, autoadaloc = is read again, and daalg → scheme has a pet migrate path.
I'm happy to split this along the July/September boundary into two PRs if that is easier to review.
🤖 Generated with Claude Code