Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
5 changes: 4 additions & 1 deletion src/pipt/update_schemes/analysis/approx.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
37 changes: 37 additions & 0 deletions tests/assimilation/test_autoadaloc.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Loading