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

pip install nowcast-bench

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

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.

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:

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:

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:

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.

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:

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:

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.

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:

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:

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:

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.

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.

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.

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:

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:

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.

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:

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:

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:

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:

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:

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:

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