Add a benchmark of the update schemes on MiniRes twin experiments

benchmarks/minires/RESULTS.md

@@ -0,0 +1,28 @@
+# MiniRes benchmark: results
+
+Generated by `run.py`; do not edit by hand. See `cases.py` for the cases, `methods.py` for the
+methods and their settings, and `scores.py` for what is scored. Each cell is the mean over
+the twin experiments (truths), ± its standard error. The rows are not ranked: which case
+and which score matters is for the reader to judge.
+
+## Case `coarse`: 25 x 25 cells
+
+20 truths, 50 members; history of 16 of 40 steps, observed at steps [4, 8, 12, 16]. PET at `b731be7`, 2026-10-27.
+
+| | misfit | RMSE forecast WWCT | logK RMSE | logK corr | logK spread | sor error | sor spread | vo error | vo spread | simulations |
+|---|---|---|---|---|---|---|---|---|---|---|
+| prior | 91.7 ± 4.2 | 0.221 ± 0.021 | 0.783 ± 0.019 | 0.0201 ± 0.027 | 0.8 ± 0.0028 | 0.0314 ± 0.0058 | 0.0383 ± 0.00077 | 0.676 ± 0.15 | 1.02 ± 0.023 |  |
+| ES-MDA (4 steps, α = 4) | 1.68 ± 0.083 | 0.0847 ± 0.013 | 0.829 ± 0.016 | 0.205 ± 0.037 | 0.526 ± 0.0041 | 0.0349 ± 0.0059 | 0.0221 ± 0.0007 | 0.602 ± 0.094 | 0.469 ± 0.012 | 250 |
+| ES-MDA (4 steps, α = 4), adaptive localization (universal threshold) | 29.6 ± 1.7 | 0.146 ± 0.02 | 0.779 ± 0.019 | 0.122 ± 0.024 | 0.773 ± 0.0031 | 0.0302 ± 0.0052 | 0.0339 ± 0.001 | 0.531 ± 0.08 | 0.894 ± 0.025 | 250 |
+| ES-MDA (4 steps, α = 4), adaptive localization (3σ threshold) | 17.7 ± 0.91 | 0.123 ± 0.018 | 0.777 ± 0.019 | 0.167 ± 0.024 | 0.752 ± 0.0029 | 0.0253 ± 0.0053 | 0.0329 ± 0.00071 | 0.43 ± 0.083 | 0.773 ± 0.022 | 250 |
+
+## Case `fine`: 50 x 50 cells
+
+8 truths, 50 members; history of 16 of 40 steps, observed at steps [4, 8, 12, 16]. PET at `b731be7`, 2026-10-27.
+
+| | misfit | RMSE forecast WWCT | logK RMSE | logK corr | logK spread | sor error | sor spread | vo error | vo spread | simulations |
+|---|---|---|---|---|---|---|---|---|---|---|
+| prior | 97 ± 7.6 | 0.243 ± 0.029 | 0.815 ± 0.022 | -0.0673 ± 0.03 | 0.793 ± 0.0038 | 0.0282 ± 0.008 | 0.0402 ± 0.0012 | 0.583 ± 0.12 | 1.03 ± 0.035 |  |
+| ES-MDA (4 steps, α = 4) | 1.9 ± 0.21 | 0.0867 ± 0.015 | 0.869 ± 0.035 | 0.179 ± 0.059 | 0.511 ± 0.0037 | 0.0326 ± 0.01 | 0.0253 ± 0.001 | 0.356 ± 0.14 | 0.47 ± 0.015 | 250 |
+| ES-MDA (4 steps, α = 4), adaptive localization (universal threshold) | 47.7 ± 5.4 | 0.2 ± 0.026 | 0.806 ± 0.021 | 0.0307 ± 0.042 | 0.777 ± 0.0026 | 0.0328 ± 0.0082 | 0.0364 ± 0.00088 | 0.476 ± 0.15 | 0.926 ± 0.039 | 250 |
+| ES-MDA (4 steps, α = 4), adaptive localization (3σ threshold) | 20.2 ± 3.5 | 0.144 ± 0.021 | 0.794 ± 0.019 | 0.151 ± 0.04 | 0.745 ± 0.0018 | 0.0358 ± 0.0099 | 0.0355 ± 0.0013 | 0.495 ± 0.12 | 0.766 ± 0.035 | 250 |

benchmarks/minires/cases.py

@@ -0,0 +1,132 @@
+"""The cases: one waterflood, at two resolutions, observed over its history and scored on its forecast.
+
+The reservoir is a 1000 x 1000 x 20 m box in metric units (m, day, bar, mD, cP), with
+one injector at the centre on rate, four producers towards the corners on bottom-hole
+pressure, 40 report steps of 60 days. The producers' rates, and when they water out, are
+therefore outcomes of the permeability around them.
+
+The unknowns are the log-permeability field and two scalars of the fluid, the residual oil
+saturation ``sor`` and the oil viscosity ``vo``. Each twin experiment draws its truth from
+the same prior as the ensemble, simulates it to the end, and perturbs the history part into
+observations; the forecast part is what the posterior is scored against, unseen.
+"""
+
+from dataclasses import dataclass
+
+import numpy as np
+
+from misc.sampling import random_stream
+from misc.structures.layout import StateLayout
+from pipt.misc_tools.extract_tools import extract_prior_info
+from simulator.minires import MiniRes
+
+LX = 1000.0  # extent [m]
+DZ = 20.0  # thickness [m]: MiniRes is areal, so the rates are per unit thickness
+
+DT = 60.0  # report step [days]
+NSTEPS = 40
+HISTORY = 16  # the observed steps, of which every fourth is assimilated
+OBS_STEPS = list(range(4, HISTORY + 1, 4))
+FORECAST_STEPS = list(range(HISTORY + 1, NSTEPS + 1))
+
+PRODUCERS = ["PROD1", "PROD2", "PROD3", "PROD4"]
+INJECTOR = "INJ"
+DATATYPES = [f"WWCT:{w}" for w in PRODUCERS] + [f"WOPR:{w}" for w in PRODUCERS] + [f"WBHP:{INJECTOR}"]
+
+# Observation error, as (relative, absolute): the larger applies.
+OBS_ERROR = {"WWCT": (0.05, 0.02), "WOPR": (0.10, 5.0), "WBHP": (0.01, 1.0)}
+
+# The prior. log K about log(200 mD), sd 0.8, spherical variogram of range RANGE metres.
+LOGK_MEAN, LOGK_SD, RANGE = np.log(200.0), 0.8, 300.0
+SCALARS = {
+    #       mean   sd     limits
+    "sor": (0.22, 0.04, [0.05, 0.40]),
+    "vo": (5.0, 1.0, [1.5, 10.0]),
+}
+
+
+@dataclass(frozen=True)
+class Case:
+    name: str
+    n: int  # cells per side
+    truths: int  # the number of twin experiments run by default
+
+    @property
+    def dx(self):
+        return LX / self.n
+
+    def model(self) -> dict:
+        """``[simulator.model]``, less what the state supplies (K, and the scalars' values)."""
+        wells = [{"name": INJECTOR, "xy": [LX / 2, LX / 2], "rate": 1000.0 / DZ, "rw": 0.15}]
+        for name, (fx, fy) in zip(PRODUCERS, [(0.15, 0.15), (0.85, 0.15), (0.15, 0.85), (0.85, 0.85)]):
+            wells.append({"name": name, "xy": [fx * LX, fy * LX], "bhp": 250.0, "rw": 0.15})
+        return {
+            "Lx": LX, "Ly": LX, "Nx": self.n, "Ny": self.n,
+            "por": 0.2, "ct": 1e-5, "cdarcy": 0.008527,
+            "fluid": {"vw": 1.0, "vo": SCALARS["vo"][0], "swc": 0.2, "sor": SCALARS["sor"][0]},
+            "wells": wells,
+        }
+
+    def simulator(self, reportpoint) -> dict:
+        """The ``[simulator]`` section, reporting at ``reportpoint``."""
+        return {
+            "dt": DT, "reporttype": "steps", "reportpoint": list(reportpoint),
+            "datatype": DATATYPES, "parallel": 1, "p0": 300.0, "s0": 0.2,
+            "model": self.model(),
+        }
+
+    def ensemble(self, ne: int, seed: int, mean_file: str) -> dict:
+        """The ``[ensemble]`` section; ``mean_file`` is where :meth:`write_prior_mean` put the field's mean."""
+        ens = {
+            "ne": ne, "seed": seed, "disable_tqdm": True, "save_prior": False,
+            "state": ["permx", *SCALARS],
+            "prior_permx": {
+                "vario": "sph", "mean": mean_file, "var": LOGK_SD**2,
+                "range": RANGE / self.dx, "aniso": 1.0, "angle": 0.0, "grid": [self.n, self.n, 1],
+            },
+        }
+        for name, (mean, sd, limits) in SCALARS.items():
+            ens[f"prior_{name}"] = {"mean": mean, "var": sd**2, "limits": limits}
+        return ens
+
+    def write_prior_mean(self, path: str) -> None:
+        np.savez(path, permx=np.full(self.n**2, LOGK_MEAN))
+
+    def run(self, state: dict) -> np.ndarray:
+        """Simulate one member to the end: ``(NSTEPS, len(DATATYPES))``, or ``nan`` if it failed."""
+        sim = MiniRes(self.simulator(range(1, NSTEPS + 1)))
+        records = sim.run_fwd_sim(state, 0)
+        if records is False:
+            return np.full((NSTEPS, len(DATATYPES)), np.nan)
+        return np.array([[row[d][0] for d in DATATYPES] for row in records])
+
+
+CASES = {c.name: c for c in [Case("coarse", 25, truths=20), Case("fine", 50, truths=8)]}
+
+
+@dataclass
+class Twin:
+    """One twin experiment: a truth, its full response, and the observations of its history."""
+
+    truth: dict  # the true state, {variable: values}
+    response: np.ndarray  # (NSTEPS, len(DATATYPES)), noiseless
+    obs: np.ndarray  # (len(OBS_STEPS), len(DATATYPES)), perturbed
+    std: np.ndarray  # the observation error, same shape
+
+
+def twin(case: Case, seed: int, mean_file: str) -> Twin:
+    """Draw a truth from the prior, simulate it, and perturb the history into observations."""
+    keys = case.ensemble(1, seed, mean_file)
+    enX, layout = StateLayout.from_prior_info(extract_prior_info(keys), 1, rng=random_stream(seed), save=False)
+    truth = {k: v[:, 0] for k, v in layout.to_dict(enX).items()}
+    response = case.run(truth)
+
+    exact = response[np.array(OBS_STEPS) - 1]
+    rel, ab = np.array([OBS_ERROR[d.split(":")[0]] for d in DATATYPES]).T
+    std = np.maximum(rel * np.abs(exact), ab)
+    obs = exact + std * np.random.default_rng(seed).standard_normal(exact.shape)
+    # A water cut outside [0, 1], or a negative rate, would be unphysical
+    quantity = [d.split(":")[0] for d in DATATYPES]
+    lower = np.where(np.isin(quantity, ["WWCT", "WOPR"]), 0.0, -np.inf)
+    upper = np.where(np.isin(quantity, ["WWCT"]), 1.0, np.inf)
+    return Twin(truth, response, obs.clip(lower, upper), std)

benchmarks/minires/methods.py

@@ -0,0 +1,36 @@
+"""The methods compared: each a label and its ``[dataassim]`` settings.
+
+The settings are PET's own, as a config would give them. Unless a row says it is tuned,
+they are defaults or the textbook choice, the same for every case.
+"""
+
+COMMON = {"energy": 99.0, "obsname": "steps"}
+
+METHODS = [
+    {
+        "label": "ES-MDA (4 steps, α = 4)",
+        "dataassim": {
+            "scheme": "esmda",
+            "analysis": "approx",
+            "mda": {"tot_assim_steps": 4, "inflation_param": [4, 4, 4, 4]},
+        },
+    },
+    {
+        "label": "ES-MDA (4 steps, α = 4), adaptive localization (universal threshold)",
+        "dataassim": {
+            "scheme": "esmda",
+            "analysis": "approx",
+            "mda": {"tot_assim_steps": 4, "inflation_param": [4, 4, 4, 4]},
+            "localization": {"name": "autoadaloc", "threshold": "universal", "type": "hard"},
+        },
+    },
+    {
+        "label": "ES-MDA (4 steps, α = 4), adaptive localization (3σ threshold)",
+        "dataassim": {
+            "scheme": "esmda",
+            "analysis": "approx",
+            "mda": {"tot_assim_steps": 4, "inflation_param": [4, 4, 4, 4]},
+            "localization": {"name": "autoadaloc", "threshold": "adaptive", "cutoff": 3.0, "type": "hard"},
+        },
+    },
+]