Metadata-Version: 2.4
Name: nowcast-bench
Version: 0.1.0
Summary: Evaluation framework for precipitation nowcasting: radar extrapolation, NWP, and blends
Author-email: Yuxiang Ren <xiaoheibuku@gmail.com>
License-Expression: MIT
Project-URL: Homepage, https://github.com/Yuxiang-Ren-HUB/nowcast-bench
Project-URL: Repository, https://github.com/Yuxiang-Ren-HUB/nowcast-bench
Project-URL: Issues, https://github.com/Yuxiang-Ren-HUB/nowcast-bench/issues
Keywords: nowcasting,precipitation,verification,radar,forecast-evaluation,meteorology
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Atmospheric Science
Classifier: Typing :: Typed
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Provides-Extra: dev
Requires-Dist: pytest>=7; extra == "dev"
Requires-Dist: ruff>=0.4; extra == "dev"
Dynamic: license-file

# nowcast-bench

Evaluation framework for precipitation nowcasting — radar extrapolation, NWP, and blends of
the two — built around one rule:

> **A skill score is not a result until you know how many samples it came from.**

## Why this exists

Blending radar extrapolation with NWP is supposed to beat either source alone. In practice it
is strongly case-dependent. On three real cases the ranking changed with the case, and with
lead time inside a case: ensemble blending beat every alternative at all six lead hours on a
fast-reorganising supercell, lost badly at +1h on a slow stratiform system before recovering
from +3h, and on a flash-flood case lost to both the simple blend and plain extrapolation at
the threshold that mattered.

A rule proposed after the first two cases — blending wins early when the storm is growing and
spatially concentrated — predicted that it should win from the start on the third. It did not.
That third case was a weaker test than the other two (verification only at 0/3/6h, and the NWP
field degraded to block-constant snapshots), so it is not a clean refutation of the rule. What
it is, is not enough evidence to encode the rule and act on it. So this package does not encode
it, or anything like it:

- There is **no `recommend_method(case_description)`** and there never will be. The only way to
  learn which method wins is `evaluate(methods, cases)` — an actual measured run.
- Every metric returns its **sample support** alongside its value. A CSI of 0.83 computed from
  11 hit pixels is not the same object as a CSI of 0.83 from 40 000, and the API refuses to let
  you conflate them by returning a bare float.
- Aggregating across cases requires an explicit choice of pooling (`"pool"` vs `"mean"`), because
  those two answer different questions and quietly disagree on skewed case sets.

## Status

Early. `nowcast_bench.metrics` (categorical scores + FSS), `nowcast_bench.protocol`
(cases, methods, `evaluate`/`aggregate`/`rank`), `nowcast_bench.blend` (fixed and
lead-time-dependent blends, plus weight fitting that refuses to be scored on its own fitting
cases), `nowcast_bench.uncertainty` (bootstrap intervals over *cases*, and rankings that say
when their own top two places are not distinguishable) and `nowcast_bench.report` (Markdown
rendering, no score without its support) are implemented. See `CONTRIBUTING.md` for the rules
any change has to follow.

## Install

```bash
pip install nowcast-bench
```

The only runtime dependency is NumPy. For development, clone the repository and install the test tools:

```bash
pip install -e ".[dev]"
```

## Usage

Categorical scores at a threshold. A point is an event when its value is `>= threshold`;
NaN points (and points where an optional boolean `mask` is `True`) are excluded from every
count rather than being treated as misses. `±inf` is *not* missing data — it is compared
against the threshold like any other value.

```python
import numpy as np
from nowcast_bench.metrics import categorical_scores, fss

obs = np.array([[0.0, 2.0, 3.0, 0.0],
                [0.0, 5.0, 4.0, 0.0],
                [0.0, 0.0, 0.0, 0.0],
                [0.0, 0.0, 0.0, np.nan]])   # bottom-right corner is missing
fct = np.array([[0.0, 0.0, 2.0, 1.0],
                [0.0, 6.0, 0.0, 3.0],
                [0.0, 0.0, 2.0, 0.0],
                [0.0, 0.0, 0.0, 0.0]])

scores = categorical_scores(obs, fct, threshold=1.0)
csi = scores["CSI"]
print(f"CSI = {csi.value:.3f}  (hits={csi.hits}, misses={csi.misses}, "
      f"false_alarms={csi.false_alarms}, of {csi.n_valid} valid points)")
# CSI = 0.286  (hits=2, misses=2, false_alarms=3, of 15 valid points)
```

There is no way to get the 0.286 without the counts next to it: every score is a frozen
dataclass carrying its contingency table.

An undefined score says so instead of returning a number:

```python
dry = categorical_scores(np.zeros((4, 4)), fct, threshold=1.0)
print(dry["POD"].value)        # Undefined(reason='no observed events: hits + misses == 0')
print(dry["POD"].is_defined)   # False
print(dry["POD"].false_alarms) # 5 -- the counts that made it undefined are still there
```

The Fractions Skill Score takes a square neighbourhood in grid points and reports how many
windows it was computed over, under which boundary convention, and whether it met missing data:

```python
f = fss(obs, fct, threshold=1.0, neighborhood=3)
print(f"FSS = {f.value:.3f} over {f.n_windows} windows of 3x3")
# FSS = 0.964 over 4 windows of 3x3
print(f.boundary, f.had_missing)                     # valid True

padded = fss(obs, fct, threshold=1.0, neighborhood=3, boundary="pad")
print(f"{padded.value:.3f} over {padded.n_windows} windows")
# 0.909 over 15 windows

print(f"{f.coverage_fraction:.1%} of the valid domain evaluated "
      f"({f.n_evaluated_points} of {f.n_valid_points} points)")
# 26.7% of the valid domain evaluated (4 of 15 points)
print(f"{padded.n_padded_windows} of {padded.n_windows} windows drew on padding")
# 11 of 15 windows drew on padding
```

`boundary` is optional and defaults to `"valid"`, which uses only the windows that fit entirely
inside the grid — `(rows - n + 1) * (cols - n + 1)` of them. `boundary="pad"` puts one window on
every grid point instead, padding the grid with *valid non-event* points, which is the
zero-filled uniform filter `pysteps.verification.spatialscores.fss` applies; pass it when you
need numbers comparable with pysteps or with published values.

The two are not two roundings of one number. On the 4x4 grid above they are 0.964 and 0.909, and
against pysteps on a 40x40 field the `"valid"` convention runs +0.0028 high at neighbourhood 3
and still +0.0003 at 15, while `"pad"` agrees to within 1e-12 absolute
(`tests/test_pysteps_crosscheck.py` asserts exactly that — 1e-12 is the tolerance in the
assertion, not a stand-in for "bit for bit" — and skips when pysteps is not installed, since it
is not a dependency). The default is `"valid"`, and that is a decision rather than a leftover:
treating the area outside the domain as "no rain" is a fabricated observation, and the edge
bias it introduces is invisible in the value. Whichever you pick, the choice is recorded on the
result and printed in every report, because an FSS number does not identify its own convention.

### What each convention costs

Neither convention is free, and neither price is in the number. `"valid"` never centres a
window closer than `neighborhood // 2` to an edge, so that frame of the domain is never
verified; `"pad"` covers the whole domain instead, by completing the windows on that same frame
with points nobody observed. They are the same frame, so the two costs are complements of one
another. Measured on a 60x60 domain:

| neighbourhood | `"valid"` evaluates | `"pad"` windows drawing on padding |
| --- | --- | --- |
| 3 | 93% of the domain | 7% |
| 9 | 75% | 25% |
| 15 | 59% | 41% |
| 31 | 25% | 75% |

At neighbourhood 31 three quarters of that domain takes no part in the FSS, and `n_windows`
says only `900`; that division is the one the result now does for you, which is why
`coverage_fraction` sits on it and why every rendering of it — `str(score)`, and every report
cell — prints the share beside the window count. On the same fields the two conventions'
*values* differ by at most 0.0018. **The difference between the conventions is not in the
number, it is in the area the number describes**, which is why
`coverage_fraction` (over the *valid* points, so masked and NaN points leave the denominator
rather than counting as unevaluated domain), `padded_fraction`, and the flags
`coverage_is_low` / `padding_is_heavy` sit on every `FractionsScore`. Both flags are taken
against one constant, `MAX_BOUNDARY_FRACTION` (10%), because they are one quantity; when either
fires, the result's `boundary_caveat` and the report's *FSS boundary cost* section say so with
the counts. A low coverage does not make the FSS wrong — it makes it incomparable with an FSS
measured over a domain of a different size, which is the part a value alone cannot tell you.

`had_missing` is the same idea one level down. A window's event fraction is its event count over
that window's *valid* points, not over `neighborhood ** 2`, so masked and NaN points are never
counted as non-events. With nothing missing the two denominators are the same number; when
something was missing, the flag says so and the report adds a line about it.

## Comparing methods across cases

`evaluate` runs every method on every case, `aggregate` combines the per-case results under a
strategy *you* have to name, and `rank` orders the methods while handing you the evidence for
and against the ordering. A "method" is just a name plus a function of the case's input
fields, so persistence, NWP and a blend are all the same kind of object to the protocol.

```python
import numpy as np
from nowcast_bench.protocol import Case, Method, aggregate, evaluate, rank

# Two cases of very different size, and two methods that each return one of the inputs.
small_obs = np.array([[5.0, 5.0, 5.0, 0.0],
                      [0.0, 0.0, 0.0, 0.0]])       # 3 observed event points out of 8
small_b = np.array([[0.0, 0.0, 0.0, 5.0],
                    [5.0, 5.0, 0.0, 0.0]])         # 3 forecast events, none of them right

big_obs = np.zeros((10, 10))
big_obs[:4] = 5.0                                  # 40 observed event points out of 100
big_a = np.zeros((10, 10))
big_a[:2] = 5.0                                    # 20 hits, 20 misses,  0 false alarms
big_b = np.zeros((10, 10))
big_b[:6] = 5.0                                    # 40 hits,  0 misses, 20 false alarms

cases = [Case("small", small_obs, {"a": small_obs, "b": small_b}),
         Case("big", big_obs, {"a": big_a, "b": big_b})]
methods = [Method("a", lambda case: case.inputs["a"]),
           Method("b", lambda case: case.inputs["b"])]

result = evaluate(methods, cases, threshold=1.0)

print("pool winner:", rank(result, "pool", "CSI").winner)
print("mean winner:", rank(result, "mean", "CSI").winner)
# pool winner: b
# mean winner: a
```

The two strategies name different winners, and that is why `aggregate(result, strategy)` has no
default. `"pool"` adds the contingency tables together (23/43 for `a` against 40/66 for `b`), so
the 100-point case outvotes the 8-point one; `"mean"` averages the per-case scores (0.75 for `a`
against 0.33 for `b`), giving each case one vote. Neither is the "real" answer — they answer
different questions, and picking one for you would be picking your conclusion.

A ranking never lets the headline stand on its own:

```python
ranking = rank(result, "mean", "CSI")
print(ranking.describe())
# CSI (higher is better), strategy='mean', n_cases=2
#   1. a: 0.7500  (2 cases used, 0 undefined, 0 failed)
#   2. b: 0.3333  (2 cases used, 0 undefined, 0 failed)
#   NOT unanimous: a loses on ['big']
#     case big: won by b
print(ranking.is_unanimous)        # False
print(dict(ranking.dissenting_cases))   # {'big': ('b',)}
```

Aggregated scores carry their case counts the way single scores carry their contingency table,
and a ranking carries the same counts per method, so two methods on neighbouring lines cannot
hide that one of them was scored on fewer cases:

```python
pod = aggregate(result, "mean").score("a", "POD")
print(pod.n_cases_used, pod.n_cases_undefined, pod.n_cases_failed)   # 2 0 0
print(ranking.support["b"].n_cases_used)                             # 2
```

A method that raises on a case is recorded as a failure, not skipped — `result.failures`,
`aggregate(...).has_failures` and `ranking.failures` all carry it, so a method cannot look good
by having quietly been scored on fewer cases than its rival.

There is **no** function that takes a description of a case and tells you which method to use.
That is not an omission to be filled in later: on the cases this package was built on, the
ranking changed from case to case and nothing knowable in advance predicted it.

## Blending, and where the weights come from

A blend is just another `Method`: a function of several input fields. `nowcast_bench.blend`
only builds them — `static_blend` for fixed weights, `lead_time_blend` for weights that hand
over from radar extrapolation to NWP as the lead time grows.

```python
import numpy as np
from nowcast_bench.blend import LeakageError, fit_static_weights, lead_time_blend
from nowcast_bench.protocol import Case, Method, evaluate, rank

def field(*points):
    """A 2x2 grid holding 2.0 at the given (row, col) points and 0.0 elsewhere."""
    grid = np.zeros((2, 2))
    for row, col in points:
        grid[row, col] = 2.0
    return grid

handover = lead_time_blend("radar", "nwp", half_life=30.0)
at_30 = Case("t+30", field((0, 0)), {"radar": field((0, 0)), "nwp": field((0, 1))},
             lead_time=30.0)
print("lead 30 min:", handover.fn(at_30).tolist())
# lead 30 min: [[1.0, 1.0], [0.0, 0.0]]
```

At one half-life the radar weight is `0.5 ** 1 = 0.5` and NWP takes the rest, so both fields
enter at half strength. A case without a `lead_time` makes this method raise: the lead time is
a physical parameter of the blend, not a case feature for choosing between methods.

Fixed weights have to come from somewhere. `fit_static_weights` sweeps them on cases you name,
and the result remembers which cases those were:

```python
fit_cases = [Case("f1", field((0, 0), (0, 1)), {"radar": field((0, 0)), "nwp": field((0, 1))}),
             Case("f2", field((1, 0), (1, 1)), {"radar": field((1, 0)), "nwp": field((1, 1))})]

fitted = fit_static_weights(fit_cases, ("radar", "nwp"),
                            threshold=1.0, strategy="mean", score="CSI", n_steps=5)

print("weights:", dict(fitted.weights))
print("fitted on:", sorted(fitted.fitted_on))
for point in fitted.trace:
    print(f"  w(radar)={point.weights['radar']:.2f}  CSI={point.score.value:.4f}")
# weights: {'radar': 0.5, 'nwp': 0.5}
# fitted on: ['f1', 'f2']
#   w(radar)=0.00  CSI=0.5000
#   w(radar)=0.25  CSI=0.5000
#   w(radar)=0.50  CSI=1.0000
#   w(radar)=0.75  CSI=0.5000
#   w(radar)=1.00  CSI=0.5000
```

The whole sweep comes back, not just its argmax, because a flat curve and a sharp peak mean
very different things: here the optimum is a single point, so the data really did choose the
weights. `fitted.training_score` is the CSI those weights reached *on the fitting cases*. It is
in-sample and systematically optimistic, and it is not this blend's performance.

For that you have to evaluate on other cases — and the package makes you:

```python
eval_cases = [
    Case("e1", field((0, 0), (0, 1)), {"radar": field((0, 0)), "nwp": field((0, 1))}),
    Case("e2", field((0, 0), (0, 1)), {"radar": field((0, 0), (0, 1)), "nwp": field((1, 1))}),
]

blend = fitted.method("blend", eval_cases)     # checked here, before anything is scored
methods = [Method("radar", lambda case: case.inputs["radar"]),
           Method("nwp", lambda case: case.inputs["nwp"]),
           blend]

ranking = rank(evaluate(methods, eval_cases, threshold=1.0), "mean", "CSI")
print(ranking.describe())
# CSI (higher is better), strategy='mean', n_cases=2
#   1. blend: 0.8333  (2 cases used, 0 undefined, 0 failed)
#   2. radar: 0.7500  (2 cases used, 0 undefined, 0 failed)
#   3. nwp: 0.2500  (2 cases used, 0 undefined, 0 failed)
#   NOT unanimous: blend loses on ['e2']
#     case e2: won by radar
```

The blend wins on average and still loses `e2`, where radar alone was already perfect and the
blend inherited NWP's false alarm — which is why the per-method case counts and the dissenting
case travel with the ordering.

Handing the fitting cases to `method()` instead is not a slightly worse experiment, it is a
different number entirely, so it is refused where it happens rather than recorded as a result:

```python
try:
    fitted.method("blend", fit_cases)
except LeakageError as exc:
    print("refused:", exc)
# refused: weights fitted on ['f1', 'f2'] would be evaluated on the same case(s) ['f1', 'f2'];
# evaluate a fitted blend only on cases it was not fitted on
```

The check is at construction time on purpose. `evaluate` records an exception raised by a
method as a `Failure`, so a leak caught only at run time would show up as a quiet "that method
failed" line inside an otherwise green evaluation instead of stopping you.

One known limitation, stated rather than corrected: **a linear blend damps extremes.** Two
strong echo cores a few grid points apart average into two weaker ones, so a blended field has
a flatter intensity distribution than either ingredient. This package does not silently apply
probability matching to hide that, because a correction stacked on the blend would make "did
blending help?" unanswerable. Apply one yourself if you want it, and evaluate it as its own
method.

## Scale decomposition, and why it is not the blend this package recommends

`static_blend` and `lead_time_blend` weight the whole field at once. `nowcast_bench.cascade`
splits both ingredients into scale bands in wavenumber space first, gives each band its own
half-life, and adds the bands back up — the reasoning being that small structures lose their
predictability before large ones do, so one weight for the whole spectrum is the wrong shape.

**This package does not hold the scale-decomposed blend to be better than the linear one, and
nothing here should be read as saying it is.** On the real cases this package was built on, a
scale-decomposed ensemble blend was clearly *behind* plain extrapolation at +1 h on a
slow-moving stratiform system — the ensemble mean smeared out a sharp, intense echo core that
extrapolation had kept — and only overtook it at +3 h. "It decomposes the field, so it must
score better" is exactly the assumption that result refutes. Which method wins on *your* cases
is an empirical question, and `evaluate` is the only thing here that answers it.

What the module does guarantee is that the decomposition itself is not quietly lossy: the
filters are a partition of unity in `log2` wavenumber, summing to exactly 1 at every
wavenumber, so the bands add back up to the field they came from.

```python
import numpy as np
from nowcast_bench import cascade_blend, decompose, intensity_ratio, recompose
from nowcast_bench.protocol import Case

cas_field = np.arange(64.0).reshape(8, 8) % 7
cas_parts = decompose(cas_field, 3)
print("levels:", cas_parts.n_levels, " shape:", cas_parts.shape)
print("perfect reconstruction:", bool(np.allclose(recompose(cas_parts), cas_field)))
# levels: 3  shape: (8, 8)
# perfect reconstruction: True
```

That invariant is not decoration. If the filters summed to 0.98 somewhere, then even a blend
that gave one ingredient *all* of the weight would hand back a quietly reshaped copy of it —
the forecast would have been filtered before any blending happened, with nothing in the output
to say so.

`half_lives` has no default and no "recommended" value: one half-life per band, largest scale
first, chosen by you. Small scales are usually held to decay faster, but that is a physical
assumption about your radar, your domain and your weather, and this package will not make it
on your behalf. The level count and the half-lives go into the method's name, so the
convention shows up in every table the blend appears in. Write them in any ordered container —
a list, a tuple or a `numpy` array; a `set` is refused, because the order is what decides which
half-life acts on which band.

```python
def cas_core(col):
    """An 8x8 field of zeros with one 2x2 block of 10.0 starting at that column."""
    grid = np.zeros((8, 8))
    grid[3:5, col:col + 2] = 10.0
    return grid

cas_fast, cas_slow = cas_core(1), cas_core(4)      # two cores, three columns apart
cas_case = Case("t+30", np.zeros((8, 8)),
                {"radar": cas_fast, "nwp": cas_slow}, lead_time=30.0)

cas_even = cascade_blend("radar", "nwp", n_levels=3, half_lives=(30.0, 30.0, 30.0))
cas_ratios = intensity_ratio(cas_even.fn(cas_case), cas_fast, quantiles=(1.0,))
print("method:", cas_even.name)
print(f"peak of blend / peak of radar: {cas_ratios.at(1.0).ratio:.2f}"
      f" over {cas_ratios.n_valid} points")
# method: cascade_blend(n_levels=3, half_lives=(30.0, 30.0, 30.0))
# peak of blend / peak of radar: 0.50 over 64 points
```

At one half-life on every band the weights are all `0.5 ** 1 = 0.5`, the decomposition cancels
out, and this *is* the linear blend — the cascade only differs from it when the half-lives do.
The ratio is `intensity_ratio`, which exists because "blending damps extremes" is otherwise
just a sentence: two 10.0 cores that do not overlap average into two 5.0 cores, so the peak of
the blend is half the peak of either ingredient. It is a diagnostic, not a score, and nothing
here corrects it — probability matching stacked on the blend would make "did blending help?"
unanswerable, as in the linear case.

Two costs come with the transform, and neither is visible in the output:

- **An FFT assumes the field is periodic, and precipitation is not.** The left edge of the
  domain is treated as the neighbour of the right edge, so a structure at one edge leaks into
  the opposite one. It is worst where nowcasting cares most: a strong echo core near a boundary
  puts a ghost of itself on the far side of the domain. The module does not window or taper the
  field to soften this, because a taper would alter the values over an unstated fraction of the
  domain and would be a second, unevaluated method stacked on the blend. Decompose a domain
  larger than the area you verify, and verify the interior.
- **An FFT cannot be taken over missing data, so missing data is refused rather than filled.**
  Zero, the field mean, an interpolation — every filling writes down an observation nobody
  made. A field with NaNs and a case carrying a mask are both turned away, for the same reason
  `"pad"` is not the default FSS boundary convention.

```python
cas_gappy = cas_fast.copy()
cas_gappy[0, 0] = np.nan
try:
    decompose(cas_gappy, 3)
except ValueError as exc:
    print("refused:", str(exc).split(":")[0])
# refused: field has 1 missing point(s) (NaN) of 64
```

What to do about those points is a modelling decision that belongs to you and has to be
evaluated as part of your method — not made silently inside a filter bank.

## Values the blend invented, and why nothing here clips them

The scale decomposition above was built first and measured afterwards, and measuring it turned
up a property nobody had designed for: **a per-scale blend produces values that neither of its
inputs holds, negative precipitation included.** Three constructions, all at `n_levels=3`,
`half_lives=(120.0, 30.0, 7.5)`, `lead_time=30.0`:

| construction | input range | equal half-lives: max | unequal: max / min |
| --- | --- | --- | --- |
| two offset strong centres | 0 - 40 | 40.000 | **47.262 / -5.426** |
| random field + a translation | 0.005 - 16.846 | 9.799 | **12.840 / -0.773** |
| large block + a nested core | 0 - 50 | 30.000 | **31.816 / -0.170** |

It is arithmetic rather than a bug. The output is `sum_i [w_i * fast_i + (1 - w_i) * slow_i]`.
With every `w_i` equal it collapses to `w * fast + (1 - w) * slow`, a pointwise convex
combination, which cannot leave the range of its two inputs — and that is precisely the case in
which the cascade degenerates into the linear blend. With the `w_i` different, which is the
only reason to decompose at all, it is not a convex combination and it overshoots on both
sides. The property shows up **only in the regime the method exists for**, which is why running
the method was what found it.

What makes it worth a section of its own is not the overshoot but the silence after it.
`metrics` compares every point against the threshold, and a negative rain rate is below every
positive threshold — so it is counted as *no event forecast*, exactly like a dry point. The
contingency table changes, the score changes, and nothing in the value or in its support says
the field held values that cannot exist.

Two diagnostics say it instead. `range_excursion` measures a field against a range **you**
name — this package assumes none, because "precipitation is not negative" is a statement about
precipitation and not about every field it will be handed — and `ingredient_excursion` measures
a blend against the pointwise range of its own ingredients, which is the property above,
directly:

```python
from nowcast_bench import ROUNDOFF_TOLERANCE, cascade_blend, ingredient_excursion, range_excursion

exc_uneven = cascade_blend("radar", "nwp", n_levels=3, half_lives=(120.0, 30.0, 7.5))
exc_field = exc_uneven.fn(cas_case)          # the two-core case from the section above
print(f"inputs: 0.000 to 10.000, blend: {exc_field.min():.3f} to {exc_field.max():.3f}")
floor = ROUNDOFF_TOLERANCE * 10.0             # one round-off floor at the inputs' scale (0-10)
print(range_excursion(exc_field, low=-floor).describe())
# inputs: 0.000 to 10.000, blend: -1.274 to 6.193
# 20 of 64 valid grid point(s) (31.25%) lie outside low=-1e-08, high=unchecked; 20 below, the
#   most extreme -1.27369 at (row 5, col 5), 1.27369 past the bound -1e-08

exc_even = cascade_blend("radar", "nwp", n_levels=3, half_lives=(30.0, 30.0, 30.0))
print("equal half-lives leave the ingredients' range:",
      not ingredient_excursion(exc_even.fn(cas_case), (cas_fast, cas_slow)).has_excursion)
# equal half-lives leave the ingredients' range: True
```

The lower bound sits one round-off floor below zero rather than at zero. `range_excursion`
compares exactly, so at `low=0.0` it also counts the points the two FFTs left within about
1e-16 of zero on the negative side, and how many of those there are depends on the platform's
floating-point arithmetic. The 20 points counted here are at least 0.086 below zero on every
platform; a count taken at `low=0.0` would not be reproducible.

`cascade_blend` takes an `excursion=` policy: `"allow"` (the default) hands the field back as
it is, and `"raise"` refuses it with the most extreme value and its position in the message.
The default is `"allow"` because all three constructions above went out of range — a default
that refused would leave the method unusable on the fields it exists for. What the default buys
is visibility, not silence.

**There is no `excursion="clip"`, and there will not be one.** Clipping is a correction applied
after the blend and invisible in every score computed afterwards: a gain could then come from
the blending or from the clipping, and "did blending help?" stops having an answer. It is the
same refusal `blend` makes for probability matching. Clip the returned array yourself if you
want it clipped — then the correction is part of your method, written where you can see it, and
evaluated with it.

Finally, `evaluate(methods, cases, threshold=..., valid_range=(0.0, None))` runs
`range_excursion` on every forecast field, hangs the result on the `(case, method)` it belongs
to, and `evaluation_report` prints a section naming each affected pair, how many points went
outside, how far, and how those points reached the contingency table. The excursion is counted
over exactly the points the table was counted over — a masked point, or one with no observation,
is excluded from both — so a point the score never saw is never reported as one that moved it.
An out-of-range point that *was* excluded is counted apart, as `n_excluded_outside`, and the
report names it in a sentence of its own: it moved no score, and a field whose impossible values
all sit under the mask is still not a clean field. Without `valid_range` nothing is checked and
nothing changes — including the report, which then does not claim a clean field it never looked
at.

## Reports

`nowcast_bench.report` renders those same result objects as Markdown and returns the text. It
writes no files: the core carries no I/O dependency, so where the report goes is your business.

Every score in a report goes through one internal formatter, so there is no rendering path that
prints a value without the support behind it, and an undefined score is spelled `undefined`
with the counts that made it undefined — never blanked out, dashed, or turned into `0.00`.

```python
from nowcast_bench import evaluation_report, fitted_weights_report, ranking_report

print(evaluation_report(result))   # `result` is the small/big evaluation from above
# # Evaluation report
#
# - Threshold: 1.0 (a grid point is an event when its value is >= threshold)
# - Scores: POD, FAR, CSI, BIAS, ETS
# - Neighborhood: not set, so no FSS was computed
# - Cases (2): small, big
# - Methods (2): a, b
# - Pairs: 4 of 4 scored, 0 failed
# ...
# | case | method | POD | FAR | CSI | BIAS | ETS |
# | --- | --- | --- | --- | --- | --- | --- |
# | small | a | 1.0000 [H3 M0 F0 C5] | 0.0000 [H3 M0 F0 C5] | 1.0000 [H3 M0 F0 C5] | ...
# | small | b | 0.0000 [H0 M3 F3 C2] | 1.0000 [H0 M3 F3 C2] | 0.0000 [H0 M3 F3 C2] | ...
# ...
# ## Failures
#
# None: every (case, method) pair ran to completion.
```

The threshold is in the header because a precipitation score report that does not state it is
unusable, and the failures get a section of their own because a method quietly missing from the
table is exactly how a comparison gets made to look good on fewer samples than its rival.

A ranking report carries the per-method case counts on every line and puts the per-case winners
in their own section, shouting when they contradict the aggregate:

```python
print(ranking_report(ranking, threshold=1.0))   # `ranking` is the blend ranking from above
# # Ranking report
#
# - Score: CSI (higher is better)
# - Strategy: 'mean'
# - Threshold: 1.0
# - Cases ranked: 2
# ...
# | rank | method | CSI | note |
# | --- | --- | --- | --- |
# | 1 | blend | 0.8333 [2 cases used, 0 undefined, 0 failed] |  |
# | 2 | radar | 0.7500 [2 cases used, 0 undefined, 0 failed] |  |
# | 3 | nwp | 0.2500 [2 cases used, 0 undefined, 0 failed] |  |
# ...
# **NOT unanimous.** The aggregate winner `blend` loses on 1 of 2 ranked case(s). ...
#
# | dissenting case | won there by |
# | --- | --- |
# | e2 | radar |
```

A method backed by a single case gets a `**(single case)**` mark on its line, and a `Ranking`
whose `support` is missing raises `ValueError` instead of degrading to a bare list of places:
an ordering with no case counts behind it is not something this package will print.

**Ties are reported, not broken.** Two methods that reach exactly the same score did not come
first and second, so the table numbers them with competition places — `1`, `1`, `3`, never a `2`
between two equal values — the header carries a bold `Tied at the top:` line, and the per-case
winners section refuses to print the usual `Aggregate winner ... **x**` sentence at all:

```python
tied = rank(evaluate([Method("radar", lambda c: c.inputs["radar"]),
                      Method("nwp", lambda c: c.inputs["nwp"])],
                     [Case("t1", field((0, 0), (0, 1)),
                           {"radar": field((0, 0)), "nwp": field((0, 1))})],
                     threshold=1.0), "pool", "CSI")
print(tied.is_tied_at_top, tied.winners)   # True ('nwp', 'radar')
print(tied.places)          # ((1, 'nwp', 0.5), (1, 'radar', 0.5))
```

`ordering` still puts equal values in name order, but only so that two runs give the same list;
it says nothing about which is better. Read `is_tied_at_top` and `winners` before `winner`,
which is documented as the first *line* of the ordering rather than the winner of anything.

Equality here is exact, not within a tolerance: the values being compared come out of one code
path over one set of cases, so bit-identical means "the same measurement", while a tolerance
would be a second convention — of a width nobody could justify — that merges methods which
really did differ. The price is that a gap of `1e-15` is ranked as a win; that gap is not
evidence of anything either, which is what the per-case winners and the per-method case counts
beside every line are for. When such a gap is too small to show in the table — two places whose
cells print the same four decimals — the ranking report adds a line under the table saying so,
because two identical-looking rows numbered `1` and `2` read as "the first one won" just as
loudly as a mis-numbered tie would.

The fitted-weights report keeps the in-sample number out of every table of evaluation scores,
and states how flat the sweep was — a flat curve means the data did not choose the weights:

```python
print(fitted_weights_report(fitted))
# # Fitted weights report
#
# - Fitted on (2 case(s)): f1, f2
# - Score swept: CSI (maximised)
# - Aggregation strategy: 'mean'
# - Threshold: 1.0
# - Fitted weights: radar=0.5000, nwp=0.5000
# ...
# ## Training score
#
# CSI of the fitted weights on the fitting cases: 1.0000 [2 cases used, 0 undefined, 0 failed]
#   -- **in-sample on the fitting cases**: the weights were chosen to maximise this very
#   number, so it is systematically optimistic and is **not** this blend's performance.
# ...
# ### Did the sweep determine the weights?
#
# - Best CSI: 1.0000 [2 cases used, 0 undefined, 0 failed]
# - Runner-up: 0.5000 [2 cases used, 0 undefined, 0 failed]
# - Margin between them: 0.5000 (a difference of two scores, not a score itself)
# - Candidates as good as the best (within the band stated below): 1 of 5 defined (5 swept)
# - That band is 0.0100, taken as 1% of the best score itself: a candidate counts as within
#   when its score differs from the best by at most that much. The band is a difference of two
#   scores, not a score itself.
# - Sharp optimum: no other candidate comes inside that band, so the fitting cases really did
#   pick this weighting out. Read that against the band above, which is 1% of a number this
#   sweep itself produced: when every candidate crowds into a narrow range the band narrows
#   with them, and this line can then call an optimum sharp that the margin above shows to be
#   one hair wide.
```

Had every candidate scored the same, that last line would instead read **"Flat sweep: the data
did not choose these weights"** — which is the honest reading of an argmax over a plateau.

The band is always stated, and the count line above it names no yardstick of its own, because
there are two: the band is the wider of 1% of the best score and 1% of the range the defined
candidates span, whichever that turns out to be on the sweep at hand. `best == 0.0` — the
normal state of a FAR sweep, which is minimised — is only the case where the range must win,
1% of nothing being nothing; a minimised score of 0.3 over widely spread candidates takes the
range as well. The report says which one it took, and the "sharp"/"flat" line is only as
trustworthy as that band: both yardsticks are 1% of a number the sweep itself produced, so a
sweep whose candidates all crowd together gets a band that crowds with them, and a difference
of a single false alarm can still be reported as a sharp optimum. The margin line is what to
read against it.

The common public names are re-exported at the top level (`from nowcast_bench import evaluate,
rank, ranking_report, ...`); `nowcast_bench.__all__` lists them, and each one is the very object
its submodule exports.

## How sure is that ranking?

Everything above returns point estimates, and a point estimate ranked 1 against a point
estimate ranked 2 reads like a finding whether or not the case set could tell them apart.
`nowcast_bench.uncertainty` bootstraps **cases** — never grid points, which are so strongly
spatially autocorrelated that resampling them treats one storm as hundreds of independent
observations and hands back a confidently narrow interval — and attaches the result to the
ranking.

`seed` is required and has no default, for the same reason `strategy` has none: an interval is
the output of a random procedure, and one that cannot say which draw produced it is neither
reproducible nor comparable. The seed is recorded on the result and printed in the report.

Here is what three cases actually buy you. Method `a` is perfect on two of them and hopeless on
the third; `b` is the other way round but weaker:

```python
import numpy as np
from nowcast_bench import Case, Method, evaluate, rank_with_intervals, ranking_report

obs_u = np.array([[5.0, 5.0], [0.0, 0.0]])        # two observed event points
perfect = obs_u                                    # CSI 1.0
half = np.array([[5.0, 0.0], [0.0, 0.0]])          # 1 hit, 1 miss  -> CSI 0.5
wrong = np.array([[0.0, 0.0], [5.0, 5.0]])         # 0 hits, 2 misses, 2 false alarms -> CSI 0.0

u_cases = [Case("u1", obs_u, {"a": perfect, "b": wrong}),
           Case("u2", obs_u, {"a": wrong,   "b": half}),
           Case("u3", obs_u, {"a": perfect, "b": wrong})]
u_methods = [Method("a", lambda case: case.inputs["a"]),
             Method("b", lambda case: case.inputs["b"])]

u_ranking = rank_with_intervals(evaluate(u_methods, u_cases, threshold=1.0),
                                "mean", "CSI", n_resamples=2000, seed=20240601)

print("a:", u_ranking.intervals["a"].low, u_ranking.intervals["a"].high)
print("b:", u_ranking.intervals["b"].low, u_ranking.intervals["b"].high)
print("overlap:", u_ranking.top_intervals_overlap)
print("informative:", u_ranking.intervals["a"].is_informative)
# a: 0.0 1.0
# b: 0.0 0.5
# overlap: True
# informative: False
```

Method `a` aggregates to a mean CSI of 0.6667 and its 95% interval is `[0.0000, 1.0000]` — the
entire range the score can take. That is not a badly chosen example, it is what three cases
support: resample them and you can draw `{u2, u2, u2}`, where `a` scores 0, or `{u1, u1, u3}`,
where it scores 1. Both are ordinary draws.

So the object says two separate things about it, and both are numbers rather than hedging:

```python
print(u_ranking.intervals["a"].caveat)
# 3 cases admit only 10 distinct resamples, so the endpoints are order statistics over at most
# 10 different values and cannot resolve a 2.50% tail; fewer than 5 cases is therefore not
# called informative here
```

`is_informative` goes `False` below `MIN_INFORMATIVE_CASES` (5) because a bootstrap of *n*
cases can only ever draw `C(2n-1, n)` different multisets — 3 of them at n=2, 10 at n=3, 35 at
n=4 — and a 95% interval needs ~40 to place a 2.5% tail anywhere but on the most extreme draw.
The interval is still computed and still returned: refusing to compute it would only push you
back to the bare 0.6667, which claims more, not less.

The ranking then answers the question the places cannot:

```python
print(u_ranking.describe().splitlines()[3])
# ORDER NOT RESOLVED at level 95%: ...
u_text = ranking_report(u_ranking, threshold=1.0)
print("names a winner:", "Aggregate winner" in u_text)
# names a winner: False
```

When the top two intervals overlap the report does **not** print its usual
`Aggregate winner ... **a**` sentence at all; it says the order was not resolved at that level,
and the table marks both cells `**NOT INFORMATIVE**` so they cannot be read as the same kind of
object as an interval from forty cases.

Overlap is a **conservative** verdict and only in one direction. Intervals that do *not* overlap
imply the difference is real at that level; intervals that *do* overlap imply nothing — the
difference may still be real, because the interval of a *difference* between two positively
correlated estimates is narrower than the gap between their two separate intervals. Nothing in
this package words overlap as "the two methods are the same", and nothing here computes a
p-value: no t-test, no permutation test, no scipy.

Resamples in which the score has no definition — POD over a draw of cases that all happen to
have no observed events — are counted in `n_resamples_undefined` and excluded from the
percentile rather than dropped in silence, and when *every* resample is undefined both endpoints
are `Undefined` with the reason attached, not `nan` and not `(0.0, 0.0)`.

## License

MIT
