Add comprehensive test suite for simulate_cascade and all components

# Powersheds Test Documentation

**Total: 117 tests across 7 test files** (all passing)

## Table of Contents

1. [Test Infrastructure](#test-infrastructure)

2. [test_hpf.py — HPF Interpolation (74 tests)](#test_hpfpy--hpf-interpolation-74-tests)

3. [test_simulate_timestep.py — Reservoir Constraints (10 tests)](#test_simulate_timesteppy--reservoir-constraints-10-tests)

4. [test_cascade.py — Integration Tests (13 tests)](#test_cascadepy--integration-tests-13-tests)

5. [test_river_confluence.py — Flow Routing (6 tests)](#test_river_confluencepy--flow-routing-6-tests)

6. [test_interpolation.py — Storage-Elevation (10 tests)](#test_interpolationpy--storage-elevation-10-tests)

7. [test_regression.py — Regression Tests (2 tests)](#test_regressionpy--regression-tests-2-tests)

8. [test_unit_conversion.py — Unit Conversion (2 tests)](#test_unit_conversionpy--unit-conversion-2-tests)

9. [Running the Tests](#running-the-tests)

---

## Test Infrastructure

### Shared Modules

**`tests/helpers.py`** — Shared dataclasses, factory functions, and helper utilities used across all test files.

| Export | Description |

|--------|-------------|

| `ReservoirData` | Dataclass mirroring the Rust `ReservoirData` struct for PyO3 bridge. |

| `RiverData` | Dataclass for river objects (lag-based flow routing). |

| `ConfluenceData` | Dataclass for confluence objects (instantaneous flow merging). |

| `CascadeData` | Wrapper holding dicts of reservoirs, rivers, and confluences. |

| `make_reservoir(n_hours, **overrides)` | Factory function creating a `ReservoirData` with sensible defaults. Returns a reservoir with capacity=500, initial_pool_elevation=200, head~50m, and a simple 2x3 HPF table. |

| `make_cascade(reservoirs, rivers, confluences)` | Wraps dicts into a `CascadeData`. |

| `run_single_reservoir(n_hours, **overrides)` | Shortcut: builds a one-reservoir cascade, runs `simulate_cascade`, and returns that reservoir's results. |

| `expected_storage(elevation)` | Analytically computes storage from the default linear storage-elevation curve. |

| `expected_elevation(storage)` | Analytically computes elevation from the default linear storage-elevation curve. |

| `DEFAULT_SET_STORAGE`, `DEFAULT_SET_ELEVATION` | `[0, 250, 500]` -> `[180, 195, 210]` — uniform slope of 0.06 m/Mm3. |

| `DEFAULT_HPF_H`, `DEFAULT_HPF_P`, `DEFAULT_HPF_Q` | 2 heads (40, 60) x 3 powers (0, 50, 100). At H=40: Q=1.2P cumecs. At H=60: Q=0.9P cumecs. |

**`tests/conftest.py`** — Pytest fixtures and configuration.

| Export | Scope | Description |

|--------|-------|-------------|

| `cumberland_cascade` | session | Loads the full Cumberland 8-reservoir cascade from `examples/Cumberland/` (YAML config, CSV time series, CSV storage-elevation tables, parquet HPF tables). Returns a `CascadeData` ready for `simulate_cascade()`. |

| `cumberland_results` | session | Runs `simulate_cascade(cumberland_cascade)` once and caches the result dict for all tests in the session. |

### Synthetic Test Data

The default synthetic reservoir used by `make_reservoir()` has these properties:

- **Storage-elevation curve**: Linear, 3 breakpoints. `elev = 180 + storage * 0.06`. Allows exact analytical verification.

- **HPF table**: 2x3 regular grid. At H=50m (midpoint), P=50MW: Q = 52.5 cumecs = 0.189 Mm3/hr.

- **Initial state**: pool_elevation=200m -> storage=333.33 Mm3, head=50m (tailwater=150m).

- **Constraints**: capacity=500, min_power_pool=190, max_release=2.0, min_release=0.0.

---

## `test_hpf.py` — HPF Interpolation (74 tests)

Tests the reversibility and accuracy of HPF bilinear interpolation: `compute_release` (P -> Q forward) and `compute_power` (Q -> P reverse). Uses real HPF table data from all 8 Cumberland reservoirs.

### Fixtures

| Fixture | Scope | Description |

|---------|-------|-------------|

| `hpf_tables` | module | Loads all 8 Cumberland HPF parquet files into a dict with raw and NaN-cleaned variants. |

| `unique_grids` | module | Extracts sorted unique H and P values for each reservoir with min/max bounds. |

### Helper Functions

| Function | Description |

|----------|-------------|

| `roundtrip_error(power, head, hpf_h, hpf_p, hpf_q)` | Performs P -> Q -> P' roundtrip. Returns `(recovered_power, absolute_error)`. Filters NaN, zero-release, and errors > 2 MW (clamping artifacts). |

| `get_grid_cell_corners(unique_h, unique_p, h_idx, p_idx)` | Returns (H, P) corner values bounding a grid cell. |

### `TestGridCorners` (9 tests)

Tests roundtrip at **exact grid points** (no interpolation needed).

| Test | Reservoirs | Description | Tolerance |

|------|-----------|-------------|-----------|

| `test_corner_roundtrip_wolfcreek` | WolfCreek | First 5x5 corners (skipping P=0). | < 2.0 MW |

| `test_corner_roundtrip_all_reservoirs` | All 8 (parametrized) | 4x4 sampled corners at 0%, 25%, 50%, 100% positions. | < 2.0 MW |

### `TestGridEdges` (16 tests)

Tests roundtrip along **grid edges** (one dimension interpolated).

| Test | Reservoirs | Description | Tolerance |

|------|-----------|-------------|-----------|

| `test_h_edge_interpolation` | All 8 | Exact P, interpolated H midpoints. Up to 5x5 tests. | < 1.5 MW |

| `test_p_edge_interpolation` | All 8 | Exact H, interpolated P midpoints. Up to 5x5 tests. | < 1.5 MW |

### `TestGridInterior` (8 tests)

Tests roundtrip at **interior points** (full bilinear interpolation).

| Test | Reservoirs | Description | Tolerance |

|------|-----------|-------------|-----------|

| `test_interior_interpolation` | All 8 | Both H and P interpolated at cell midpoints. Up to 5x5 cells. | < 5.0 MW |

### `TestClampedRegions` (16 tests)

Tests behavior **outside grid bounds** (clamping logic).

| Test | Reservoirs | Description | Tolerance |

|------|-----------|-------------|-----------|

| `test_below_h_min` | All 8 | H = grid_min - 5m, mid-range P. | < 5.0 MW |

| `test_above_p_max` | All 8 | P = grid_max + 10MW, mid-range H. | recovered P <= p_max + 0.1 |

### `TestStatisticalRobustness` (9 tests)

Statistical tests with large random samples.

| Test | Reservoirs | Description | Tolerance |

|------|-----------|-------------|-----------|

| `test_random_roundtrip_distribution` | All 8 | 1000 random (H, P) pairs per reservoir. Reports percentile statistics. | 95th percentile < 1.5 MW |

| `test_error_by_region` | WolfCreek | 10x10 spatial error grid, 50 samples/cell. Identifies worst region. | Diagnostic (no assertion) |

### `TestPerformance` (10 tests)

Benchmarks for HPF function speed.

| Test | Reservoirs | Description | Tolerance |

|------|-----------|-------------|-----------|

| `test_forward_function_speed` | WolfCreek | 10,000 `compute_release` calls. | < 15,000 us/call |

| `test_reverse_function_speed` | WolfCreek | 10,000 `compute_power` calls. | < 15,000 us/call |

| `test_speed_by_table_size` | All 8 | 1,000 forward calls per reservoir. | Diagnostic |

### `TestEdgeCases` (4 tests)

Boundary inputs and degenerate cases.

| Test | Description | Expected |

|------|-------------|----------|

| `test_zero_power` | `compute_release(0.0, 50.0, ...)` | Returns `0.0` |

| `test_zero_release` | `compute_power(0.0, 50.0, ...)` | Returns `0.0` |

| `test_negative_head` | `compute_power(0.01, -5.0, ...)` | Returns `0.0` |

| `test_empty_table` | Both functions with empty `[]` vectors | Returns `NaN` |

### `TestDiagnostics` (2 tests)

Physical monotonicity constraints.

| Test | Description | Expected |

|------|-------------|----------|

| `test_forward_values_progression` | Q at fixed H, increasing P | Q monotonically non-decreasing |

| `test_reverse_values_progression` | P at fixed H, increasing Q | P monotonically non-decreasing (0.001 MW tolerance) |

---

## `test_simulate_timestep.py` — Reservoir Constraints (10 tests)

Tests the per-timestep reservoir simulation logic indirectly through `simulate_cascade` with single-reservoir, short-duration configurations. Covers every constraint path in the `simulate_timestep` function.

### `TestUnconstrainedOperation` (2 tests)

| Test | Description | Key Assertions |

|------|-------------|----------------|

| `test_basic_unconstrained_operation` | Default setup: head=50m, P=50MW. Target release (~0.189 Mm3/hr) is within all bounds. | `release[0] == 0.189`, `actual_power[0] == 50.0` (exact), `spill[0] == 0.0` |

| `test_actual_power_equals_target_when_unconstrained` | 4-hour run. Verifies the Rust code skips power recomputation when release is unconstrained. | `actual_power[t] == 50.0` (exact float equality, not approximate) |

### `TestBelowPowerPool` (1 test)

| Test | Description | Key Assertions |

|------|-------------|----------------|

| `test_below_power_pool_zero_release` | initial_pool_elevation=185 (below min_power_pool=190). | `release[0] == 0.0`, `actual_power[0] == 0.0`, `storage[0] == init_storage + inflow` |

### `TestReleaseConstraints` (2 tests)

| Test | Description | Key Assertions |

|------|-------------|----------------|

| `test_max_release_binding` | target_power=100 with max_release=0.1 (far below target release of ~0.378). | `release[0] == 0.1`, `actual_power[0] < 100` (recomputed) |

| `test_min_release_binding` | target_power=5 with min_release=0.1 (above target release of ~0.019). | `release[0] == 0.1`, `actual_power[0] > 5` (more water = more power) |

### `TestInsufficientWater` (1 test)

| Test | Description | Key Assertions |

|------|-------------|----------------|

| `test_insufficient_water` | initial_pool_elevation=180.01 (storage=0.167), min_power_pool=179 (so pool check doesn't fire), zero inflow. Target release (0.216) > available water (0.167). | `release[0] == init_storage`, `storage[0] == 0.0` |

### `TestSpill` (2 tests)

| Test | Description | Key Assertions |

|------|-------------|----------------|

| `test_spill_from_capacity_overflow` | capacity=350, inflow=50, target_power=0 (release=0). potential_storage=383.33 > 350. | `storage[0] == 350.0`, `spill[0] == 33.33` |

| `test_no_spill_normal_operation` | Default setup, storage well below capacity. | `spill[0] == 0.0` |

### `TestMassBalance` (2 tests)

| Test | Description | Key Assertions |

|------|-------------|----------------|

| `test_mass_balance_single_timestep` | Default 1-hour run. | `init_storage + inflow - storage - release - spill == 0` (abs < 1e-9) |

| `test_mass_balance_multi_timestep` | 24-hour run. Checks every timestep. | Mass balance holds at all 24 timesteps (abs < 1e-9) |

---

## `test_cascade.py` — Integration Tests (13 tests)

Tests the full `simulate_cascade()` function: multi-object routing, simulation ordering, output structure, and physical invariants.

### `TestOutputStructure` (2 tests)

| Test | Description |

|------|-------------|

| `test_single_reservoir_output_structure` | 1 reservoir, 24h. Verifies result has all 10 expected fields, each of length 24. |

| `test_output_structure_rivers_confluences` | Cascade with reservoir + river + confluence. Verifies river/confluence results have `inflow` and `outflow` fields. |

### `TestErrorHandling` (1 test)

| Test | Description |

|------|-------------|

| `test_empty_reservoirs_raises` | Empty reservoirs dict raises an exception. |

### `TestTwoReservoirChain` (1 test)

| Test | Description |

|------|-------------|

| `test_two_reservoir_chain` | ResA -> River (lag=2) -> ResB, 12h. Verifies ResB sees zero upstream flow for t<2 (legacy), then lagged ResA releases for t>=2. |

### `TestSimulationOrder` (1 test)

| Test | Description |

|------|-------------|

| `test_simulation_order_respected` | ResA (order=1) -> ResB (order=2). Verifies ResB total_inflow includes ResA release+spill at each timestep (because A runs first). |

### `TestDownstreamRouting` (2 tests)

| Test | Description |

|------|-------------|

| `test_downstream_routing_reservoir_to_river` | Verifies river inflow[t] == upstream reservoir release[t] + spill[t] for all t. |

| `test_confluence_merges_two_branches` | Y-topology: ResA->RivA + ResB->RivB -> Confluence -> ResC. Verifies confluence outflow == sum of river outflows. |

### `TestMassBalance` (1 test)

| Test | Description |

|------|-------------|

| `test_mass_balance_all_timesteps` | 24-hour single reservoir. Conservation law holds at every timestep (tolerance 1e-9). |

### `TestPhysicalInvariants` (4 tests)

| Test | Description |

|------|-------------|

| `test_storage_never_negative` | 168h with high power, low inflow. `storage[t] >= 0` for all t. |

| `test_release_within_bounds` | 24h run. `min_release <= release[t] <= max_release` for all non-zero releases. |

| `test_spill_only_at_capacity` | 48h with high inflow. `spill[t] > 0` only when `storage[t] == capacity`. |

| `test_pool_elevation_within_table_range` | 24h run. `180 <= pool_elevation[t] <= 210` for all t. |

### `TestCumberland` (1 test)

| Test | Description |

|------|-------------|

| `test_no_nan_in_cumberland` | Full 8-reservoir Cumberland cascade. No NaN values in any output field for any object. |

---

## `test_river_confluence.py` — Flow Routing (6 tests)

Tests river lag-based flow routing and confluence instantaneous flow summation.

### `TestRiverLagRouting` (4 tests)

| Test | Topology | Description |

|------|----------|-------------|

| `test_river_lag_basic` | ResA -> River (lag=3, legacy=[0,0,0]) -> ResB | `outflow[t<3] == 0` (legacy), `outflow[t>=3] == inflow[t-3]`. |

| `test_river_legacy_flows_used` | ResA -> River (lag=3, legacy=[1,2,3]) -> ResB | `outflow[0:3] == [1.0, 2.0, 3.0]`. |

| `test_river_lag_zero_passthrough` | ResA -> River (lag=0, legacy=[]) -> ResB | `outflow[t] == inflow[t]` for all t. |

| `test_river_legacy_length_mismatch_raises` | River with lag=3, legacy=[1,2] | Raises exception matching "legacy_flows must equal lag". |

### `TestConfluenceRouting` (2 tests)

| Test | Topology | Description |

|------|----------|-------------|

| `test_confluence_sums_two_upstreams` | ResA->RivA + ResB->RivB -> Conf -> ResC | `conf.outflow[t] == RivA.outflow[t] + RivB.outflow[t]` |

| `test_confluence_single_upstream` | ResA -> RivA -> Conf -> ResB | `conf.outflow[t] == RivA.outflow[t]` |

---

## `test_interpolation.py` — Storage-Elevation (10 tests)

Tests `interpolate_elevation()` and `interpolate_storage()` indirectly through `simulate_cascade`. Uses `target_power=0` and `min_release=0` so release=0, making `storage[t] = storage[t-1] + inflow[t]` and `pool_elevation[t] = interpolate_elevation(storage[t])` directly verifiable.

### `TestElevationAtBreakpoints` (3 tests)

| Test | Description | Key Assertion |

|------|-------------|---------------|

| `test_elevation_at_exact_breakpoint` | initial_pool_elevation=195 (storage=250), zero inflow. | `pool_elevation[0] == 195.0` |

| `test_elevation_below_minimum_clamped` | initial_pool_elevation=180 (storage=0), zero inflow. | `pool_elevation[0] == 180.0` |

| `test_elevation_above_maximum_clamped` | initial_pool_elevation=210 (storage=500), high inflow. | `pool_elevation[0] == 210.0`, spill > 0 |

### `TestInterpolationMidpoints` (1 test)

| Test | Description | Key Assertion |

|------|-------------|---------------|

| `test_elevation_interpolation_midpoint` | Start at storage=250, add 125 Mm3. Expected: elev = 195 + 125*15/250 = 202.5. | `pool_elevation[1] == 202.5` |

### `TestRoundtrip` (1 test)

| Test | Description | Key Assertion |

|------|-------------|---------------|

| `test_initial_pool_elevation_roundtrip` | Tests 6 elevations: 182, 190, 195, 200, 205, 208. Each goes through `interpolate_storage` (at t=0 init) then `interpolate_elevation`. | All recover to within 0.01m of original. |

### `TestStorageTracking` (2 tests)

| Test | Description | Key Assertion |

|------|-------------|---------------|

| `test_storage_tracks_inflow_with_zero_release` | Start at storage=250, constant inflow=1.0, 10 hours. | `storage[t] == 250 + (t+1)*1.0` |

| `test_pool_elevation_follows_linear_curve` | Same setup. Verify elevation matches analytical curve. | `pool_elevation[t] == expected_elevation(storage[t])` |

### `TestCumberlandInterpolation` (1 test)

| Test | Description | Key Assertion |

|------|-------------|---------------|

| `test_cumberland_wolfcreek_initial_elevation` | Real WolfCreek data, initial_pool_elevation=220. | `pool_elevation[0] == 220.0` (within 0.5m) |

### `TestMonotonicity` (2 tests)

| Test | Description | Key Assertion |

|------|-------------|---------------|

| `test_monotonic_elevation_with_filling` | Low initial storage, constant inflow, zero release. | `pool_elevation[t] >= pool_elevation[t-1]` for all t |

| `test_monotonic_elevation_with_draining` | High initial storage, zero inflow, moderate power target. | `pool_elevation[t] <= pool_elevation[t-1]` for all t |

---

## `test_regression.py` — Regression Tests (2 tests)

Golden baseline tests against known-good Cumberland cascade outputs. Baseline stored in `tests/fixtures/cumberland_baseline.json`.

| Test | Description | Tolerance |

|------|-------------|-----------|

| `test_cumberland_baseline` | Runs full Cumberland cascade. Compares `final_storage`, `sum_actual_power`, `sum_spill`, and `mean_pool_elevation` for all 8 reservoirs against baseline. | rel=1e-6 (or abs=1e-6 for spill) |

| `test_smoke_test_from_docs` | Reproduces the CLAUDE.md smoke test exactly. | No NaN in `actual_power` |

To update the baseline after intentional Rust changes, regenerate `tests/fixtures/cumberland_baseline.json` by running the Cumberland cascade and saving per-reservoir summary statistics.

---

## `test_unit_conversion.py` — Unit Conversion (2 tests)

Verifies the cumecs (m3/s) to Mm3/hr conversion: `release_Mm3hr = Q_cumecs * 3600 / 1,000,000`.

| Test | Description | Key Assertion |

|------|-------------|---------------|

| `test_cumecs_to_mm3hr_at_known_grid_point` | At H=40, P=50: Q=60 cumecs -> 0.216 Mm3/hr. | `compute_release(50, 40, ...) == 0.216` |

| `test_known_conversion_1000_cumecs` | Synthetic table: Q=1000 cumecs -> 3.6 Mm3/hr. | `compute_release(100, 50, ...) == 3.6` |

---

## Running the Tests

```bash

# Run all tests

uv run python -m pytest tests/ -v

# Run only fast tests (exclude performance benchmarks)

uv run python -m pytest tests/ -v -m "not slow"

# Run a specific test file

uv run python -m pytest tests/test_cascade.py -v

# Run a specific test class or method

uv run python -m pytest tests/test_simulate_timestep.py::TestMassBalance -v

```

### Test Count by File

| File | Tests |

|------|-------|

| `test_hpf.py` | 74 |

| `test_cascade.py` | 13 |

| `test_simulate_timestep.py` | 10 |

| `test_interpolation.py` | 10 |

| `test_river_confluence.py` | 6 |

| `test_regression.py` | 2 |

| `test_unit_conversion.py` | 2 |

| **Total** | **117** |

# powersheds test suite

"""

Pytest configuration for powersheds tests.

Pytest configuration and shared fixtures for powersheds tests.

"""

from pathlib import Path

import pandas as pd

import pytest

import yaml

import powersheds

from helpers import CascadeData, ConfluenceData, ReservoirData, RiverData

# =============================================================================

# Custom Markers

# =============================================================================

def pytest_configure(config):

    """Configure custom markers."""

    config.addinivalue_line(

        "markers", "slow: marks tests as slow (deselect with '-m \"not slow\"')"

    )

# =============================================================================

# Cumberland Fixtures

# =============================================================================

CUMBERLAND_DIR = Path(__file__).parent.parent / "examples" / "Cumberland"

@pytest.fixture(scope="session")

def cumberland_cascade():

    """Load the full Cumberland 8-reservoir cascade from example data."""

    config_path = CUMBERLAND_DIR / "cascade_config.yaml"

    with open(config_path) as f:

        config_dict = yaml.safe_load(f)

    reservoir_dict = {}

    river_dict = {}

    confluence_dict = {}

    for name, specs in config_dict.items():

        if specs["object_type"] == "reservoir":

            ts = pd.read_csv(CUMBERLAND_DIR / "time_series" / f"{name}.csv")

            se = pd.read_csv(

                CUMBERLAND_DIR / "storage_HWelevation_tables" / f"{name}.csv"

            )

            hpf = pd.read_parquet(

                CUMBERLAND_DIR / "head_power_flow_tables" / name

            )

            reservoir_dict[name] = ReservoirData(

                **specs,

                catchment_inflow=ts["catchment_inflow"].tolist(),

                target_power=ts["target_power"].tolist(),

                set_storage=se["storage_Mm3"].tolist(),

                set_elevation=se["elevation_m"].tolist(),

                hpf_h=hpf["H_m"].astype(float).tolist(),

                hpf_p=hpf["P_MW"].astype(float).tolist(),

                hpf_q=hpf["Q_cumecs"].astype(float).tolist(),

            )

        elif specs["object_type"] == "river":

            legacy = pd.read_csv(

                CUMBERLAND_DIR / "time_series" / f"{name}.csv"

            )

            river_dict[name] = RiverData(

                **specs,

                legacy_flows=legacy["legacy_flow"].tolist(),

            )

        elif specs["object_type"] == "confluence":

            confluence_dict[name] = ConfluenceData(**specs)

    return CascadeData(

        reservoirs=reservoir_dict,

        rivers=river_dict,

        confluences=confluence_dict,

    )

@pytest.fixture(scope="session")

def cumberland_results(cumberland_cascade):

    """Run the Cumberland cascade and cache results for the session."""

    return powersheds.simulate_cascade(cumberland_cascade)

{

  "WolfCreek": {

    "final_storage": 5086.765250757235,

    "sum_actual_power": 14509.02,

    "sum_spill": 0.0,

    "mean_pool_elevation": 220.33568136836752

  },

  "DaleHollow": {

    "final_storage": 2047.3999999999971,

    "sum_actual_power": 0.0,

    "sum_spill": 0.0,

    "mean_pool_elevation": 201.13975687529637

  },

  "CordellHull": {

    "final_storage": 540.3762189139461,

    "sum_actual_power": 4693.02,

    "sum_spill": 121.46417359346776,

    "mean_pool_elevation": 157.1643329717268

  },

  "CenterHill": {

    "final_storage": 2719.2148359068806,

    "sum_actual_power": 3635.94,

    "sum_spill": 0.0,

    "mean_pool_elevation": 210.18284075212327

  },

  "OldHickory": {

    "final_storage": 857.0,

    "sum_actual_power": 0.0,

    "sum_spill": 419.99814714439424,

    "mean_pool_elevation": 138.7

  },

  "JPercyPriest": {

    "final_storage": 620.7600000000002,

    "sum_actual_power": 0.0,

    "sum_spill": 0.0,

    "mean_pool_elevation": 150.9903683783322

  },

  "Cheatham": {

    "final_storage": 594.2905274140278,

    "sum_actual_power": 2270.82,

    "sum_spill": 0.0,

    "mean_pool_elevation": 123.3456746352141

  },

  "Barkley": {

    "final_storage": 1423.766365644019,

    "sum_actual_power": 5460.54,

    "sum_spill": 0.0,

    "mean_pool_elevation": 110.13269757075129

  }

}

 No newline at end of file