diff --git a/CHANGELOG.md b/CHANGELOG.md index b3bc8876..b6d3d815 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ and versions follow [Semantic Versioning](https://semver.org/spec/v2.0.0.html). - **`analysis = "subspace2"`**, the ensemble-transform IES of Raanes, Stordal & Evensen (2019), on ES-MDA, LM-EnRML and GN-EnRML. It solves for the `ne × ne` transform `W` directly, starting from `W = I`, and uses the analytic data covariance through `scale_data` rather than the ensemble representation `E Eᵀ` that `subspace` uses — so it takes no SVD and reads neither `energy` nor `iteration.energy`. It is exactly `margis` with the marginalised error scale `Ratio` fixed at 1, i.e. with the data uncertainty taken as known; `tests/assimilation/test_subspace2.py` pins that identity bit-for-bit. Ported from upstream 6f313d7, with two departures: the transform is initialised at iteration 0 rather than 1 (schemes here count from 0, so the reference version never initialises and dies on the first call), and the whitened observations are recomputed every call instead of cached on the first. ES-MDA redraws its observations and their scale at every assimilation step, so the cache would drive later steps with the first step's observations whitened by the first step's factor — the run still completes and the misfit still falls, which is why it would have gone unnoticed. On the characterisation case it is the difference between a posterior misfit of 175.7 and 201.1; the cached, lower number is the wrong one. The 39 existing golden arrays are unchanged. ### Fixed +- **`autoadaloc` with `projection = "ensemble"` runs, and computes the update it describes.** The `approx` analysis built the projected anomalies as `X2 @ D_anom`, an `(nr, ne) @ (nd, ne)` product that raises `ValueError` unless `nd == ne`, and that would count the innovations twice even when it ran. They are now `V_r Σ_r X2`, the transform the non-localized branch already uses, so an all-ones taper gives the non-localized step exactly (#170). - **A `LOCALIZATION` block that names no mode runs again.** Localization was selected by which keyword appeared in the block -- `autoadaloc`, `localanalysis`, `dist_loc`, a pickled mask file, or none of them for the parallel update -- and the rewrite replaced that with a required `name` key. Every config written before the rewrite therefore stopped at startup with `Localization config has no 'name'`, naming three strings its author had never seen, and `pet migrate` did not cover the block. The mode is now inferred from the keyword that used to select it, and an explicit `name` still wins. - **A crashed simulation costs the optimizer its trial point, not the whole run.** `EnsembleOptimizationBase.function` raised `RuntimeError` whenever `calc_prediction` reported failure, so one control vector the simulator could not run ended the optimization and discarded every iteration before it. It now reports `inf`, which no backtracking step can improve on, so the point is rejected and the run continues from the last good iterate — the behaviour `765d637` removed. A crashed *single-point* evaluation additionally leaves `stateF` at its last good value: the gradient is `enF - repeat(stateF, nr)`, so writing `inf` there would have poisoned every later gradient with `inf`/`NaN` instead of rejecting one point. The separate abort when *every* member of a forecast fails is unchanged — there is nothing left to compute a gradient from. - **A pickled localization file can be used again.** Those 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. They were loaded and returned unconverted, so the first code to ask for `.taper` raised `AttributeError: 'dict' object has no attribute 'taper'` and no pickled mask file worked at all. They are now converted to `LocalizationEntry`, including older files that wrote `range` as the radius alone. diff --git a/src/pipt/update_schemes/analysis/approx.py b/src/pipt/update_schemes/analysis/approx.py index de5f1324..c864b1de 100644 --- a/src/pipt/update_schemes/analysis/approx.py +++ b/src/pipt/update_schemes/analysis/approx.py @@ -102,7 +102,10 @@ def update(self, enX, enY, enE, **kwargs): return AnalysisResult(step=Cxy_loc @ X2) # shape: (nx, ne) elif y_proj == 'ensemble': - Y_anom_proj = X2 @ D_anom # shape: (ne, ne) + # W = V_r S_r (I + S_r^2)^-1 U_r^T D: the innovations enter + # once, through X2, so an all-ones taper gives the + # non-localized step X_anom @ X3. + Y_anom_proj = (VrT.T * Sr[None, :]) @ X2 # shape: (ne, ne) T_loc = localization( # shape: (nx, ne) X = scx[:, None]*X_anom, # shape: (nx, ne) Y = Y_anom_proj diff --git a/tests/assimilation/test_autoadaloc.py b/tests/assimilation/test_autoadaloc.py index 58cde6ca..e09cfbcd 100644 --- a/tests/assimilation/test_autoadaloc.py +++ b/tests/assimilation/test_autoadaloc.py @@ -172,3 +172,40 @@ def __init__(self, localization): np.testing.assert_allclose(step_loc, step_loc_expected) np.testing.assert_allclose(step_no_loc, step_expected_no_loc) assert not np.array_equal(step_loc, step_no_loc) + + +def test_approx_update_autoadaloc_ensemble_projection(): + # Regression for #170: the 'ensemble' projection multiplied X2 by D_anom a + # second time, which fails unless nd == ne. Use nd != ne != nr, and an + # all-ones taper so the step must equal the non-localized one. + rng = np.random.default_rng(170) + nx, nd, ne = 6, 15, 10 + enX = rng.normal(size=(nx, ne)) + enY = rng.normal(size=(nd, ne)) + enE = enY.mean(axis=1)[:, None] + rng.normal(0, 0.1, size=enY.shape) + + class FakeScheme: + lam = 1.0 + trunc_energy = 0.9 + cov_data = 0.1*np.ones(nd) + keys_da = {"emp_cov": False} + + def __init__(self, localization): + self.localization = localization + + loc_info = { + "name": "autoadaloc", + "field": [nx], + "actnum": None, + "threshold": "fixed", + "cutoff": 0.0, + "type": "hard", + "projection": "ensemble", + } + step_loc = approx_update(FakeScheme(AutoAdaptiveLocalization(loc_info))).update(enX, enY, enE).step + step_no_loc = approx_update( + FakeScheme(type('localization', (object,), {'name': None})()) + ).update(enX, enY, enE).step + + assert step_loc.shape == (nx, ne) + np.testing.assert_allclose(step_loc, step_no_loc)