Benchmark: a row for ES, and the CRPS

benchmarks/minires/methods.py

@@ -9,6 +9,10 @@
 COMMON = {"energy": 99.0, "obsname": "steps"}
 
 METHODS = [
+    {
+        "label": "ES (one update)",
+        "dataassim": {"scheme": "es", "analysis": "approx"},
+    },
     {
         "label": "ES-MDA (4 steps, α = 4)",
         "dataassim": {

benchmarks/minires/run.py

@@ -107,10 +107,12 @@ def experiment(case_name: str, method: dict, i: int, ne: int) -> dict:
 COLUMNS = {
     "misfit": "misfit",
     "forecast_wwct_rmse": "RMSE forecast WWCT",
+    "forecast_wwct_crps": "CRPS forecast WWCT",
     "logk_rmse": "logK RMSE",
     "logk_corr": "logK corr",
     "logk_spread": "logK spread",
-    **{f"{n}_{s}": f"{n} {s}" for n in SCALARS for s in ("error", "spread")},
+    "logk_crps": "logK CRPS",
+    **{f"{n}_{s}": f"{n} {t}" for n in SCALARS for s, t in [("error",) * 2, ("spread",) * 2, ("crps", "CRPS")]},
     "runs": "simulations",
 }
 
@@ -134,10 +136,12 @@ def render_case(title: str, res: dict) -> list:
     ]
     any_method = next(iter(runs.values()))
     lines += ["", "| | " + " | ".join(COLUMNS.values()) + " |", "|---" * (len(COLUMNS) + 1) + "|"]
-    prior = [cell(any_method["prior"][k]) if k != "runs" else "" for k in COLUMNS]
+    prior = [cell(any_method["prior"][k]) if k in any_method["prior"] else "" for k in COLUMNS]
     lines += ["| prior | " + " | ".join(prior) + " |"]
     for label, r in runs.items():
-        post = [cell(r["posterior"][k]) if k != "runs" else f"{np.mean(r['runs']):.0f}" for k in COLUMNS]
+        # A score added since the official results were run is blank in them
+        post = [cell(r["posterior"][k]) if k in r["posterior"] else "" for k in COLUMNS]
+        post[list(COLUMNS).index("runs")] = f"{np.mean(r['runs']):.0f}"
         lines += [f"| {label} | " + " | ".join(post) + " |"]
     return lines
 
@@ -214,8 +218,10 @@ def main():
 
         if quick:
             print("\n".join(render_case(f"## Case `{case_name}`: this run", results)))
+            official = {"methods": {}}
             if old and old["meta"]["ne"] == args.ne and old["meta"]["truths"] >= truths:
                 official = first_truths(old, truths, results["methods"])
+            if official["methods"]:
                 print("\n".join(render_case(f"## Case `{case_name}`: official, on the same truths", official)))
             print()
             continue

benchmarks/minires/scores.py

@@ -12,6 +12,18 @@
 WWCT = [i for i, d in enumerate(DATATYPES) if d.startswith("WWCT")]
 
 
+def crps(ensemble: np.ndarray, truth) -> np.ndarray:
+    """The continuous ranked probability score of each row of ``ensemble`` (members along axis 1).
+
+    ``E|X - y| - E|X - X'| / 2``, the latter by sorting: the mean absolute difference of the
+    members is ``2/ne^2 sum_i (2i - ne - 1) x_(i)``. In the units of the variable; zero is perfect.
+    """
+    ne = ensemble.shape[1]
+    x = np.sort(ensemble, axis=1)
+    spread = 2 / ne**2 * (x @ (2 * np.arange(1, ne + 1) - ne - 1))
+    return np.mean(np.abs(ensemble - np.atleast_1d(truth)[:, None]), axis=1) - spread / 2
+
+
 def scores(states: dict, responses: np.ndarray, twin: Twin) -> dict:
     """Score one ensemble.
 
@@ -34,17 +46,21 @@ def scores(states: dict, responses: np.ndarray, twin: Twin) -> dict:
     forecast = np.nanmean(responses[:, np.array(FORECAST_STEPS) - 1][..., WWCT], axis=0)
     truth = twin.response[np.array(FORECAST_STEPS) - 1][:, WWCT]
     out["forecast_wwct_rmse"] = float(np.sqrt(np.mean((forecast - truth) ** 2)))
+    members = responses[:, np.array(FORECAST_STEPS) - 1][..., WWCT].reshape(len(responses), -1).T
+    out["forecast_wwct_crps"] = float(np.nanmean(crps(members, truth.ravel())))
 
     # The field: error and correlation of the ensemble mean, and the spread left
     logk, true_logk = states["permx"], twin.truth["permx"]
     mean = logk.mean(axis=1)
     out["logk_rmse"] = float(np.sqrt(np.mean((mean - true_logk) ** 2)))
     out["logk_corr"] = float(np.corrcoef(mean, true_logk)[0, 1])
     out["logk_spread"] = float(logk.std(axis=1, ddof=1).mean())
+    out["logk_crps"] = float(crps(logk, true_logk).mean())
 
     # The scalars: error of the ensemble mean, and the spread left
     for name in SCALARS:
         values = states[name][0]
         out[f"{name}_error"] = float(abs(values.mean() - twin.truth[name][0]))
         out[f"{name}_spread"] = float(values.std(ddof=1))
+        out[f"{name}_crps"] = float(crps(values[None], twin.truth[name][0])[0])
     return out