#!/usr/bin/env python3
from __future__ import annotations

import argparse
import json
from pathlib import Path
from typing import Any


def load_json(path: Path) -> Any:
    return json.loads(path.read_text(encoding="utf-8"))


def write_json(path: Path, obj: Any) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(json.dumps(obj, indent=2), encoding="utf-8")


def write_text(path: Path, text: str) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(text, encoding="utf-8")


def interval_hhmm(interval_min: int) -> str:
    if interval_min <= 0:
        raise ValueError("--interval-min must be > 0")
    hours = interval_min // 60
    minutes = interval_min % 60
    return f"{hours}:{minutes:02d}"


def build_raingage_text(
    *,
    gage_id: str,
    rain_format: str,
    interval_min: int,
    scf: float,
    series_name: str,
) -> str:
    hhmm = interval_hhmm(interval_min)
    lines = [
        "[RAINGAGES]",
        ";;Name             Format     Interval   SCF      Source",
        f"{gage_id:<18} {rain_format:<10} {hhmm:<10} {scf:<8g} TIMESERIES {series_name}",
    ]
    return "\n".join(lines) + "\n"


def parse_series_by_station(climate_obj: Any) -> dict[str, str]:
    stations = climate_obj.get("stations")
    if not isinstance(stations, list):
        return {}
    out: dict[str, str] = {}
    for idx, item in enumerate(stations, start=1):
        if not isinstance(item, dict):
            continue
        station_id = str(item.get("station_id") or "").strip()
        series_name = str(item.get("series_name") or "").strip()
        if not station_id or not series_name:
            continue
        if station_id in out:
            raise ValueError(f"Rainfall JSON has duplicate station_id '{station_id}' in stations[{idx}]")
        out[station_id] = series_name
    return out


def main() -> None:
    ap = argparse.ArgumentParser(description="Build a deterministic SWMM [RAINGAGES] helper snippet.")
    ap.add_argument("--gage-id", default="RG1")
    ap.add_argument("--series-name", default=None)
    ap.add_argument(
        "--station-id",
        default=None,
        help="Optional station ID when --rainfall-json contains multiple stations.",
    )
    ap.add_argument("--rainfall-json", type=Path, default=None, help="Optional JSON from format_rainfall.py")
    ap.add_argument("--rain-format", default="INTENSITY", choices=["INTENSITY", "VOLUME", "CUMULATIVE"])
    ap.add_argument("--interval-min", type=int, default=5)
    ap.add_argument("--scf", type=float, default=1.0, help="Snow catch deficiency factor")
    ap.add_argument("--out-text", type=Path, required=True)
    ap.add_argument("--out-json", type=Path, required=True)
    args = ap.parse_args()

    series_name = args.series_name
    if args.rainfall_json is not None:
        climate = load_json(args.rainfall_json)
        series_by_station = parse_series_by_station(climate)

        json_series = ""
        if args.station_id is not None:
            requested_station = str(args.station_id).strip()
            if not requested_station:
                raise ValueError("--station-id cannot be blank when provided")
            if not series_by_station:
                raise ValueError(
                    "--station-id was provided, but rainfall JSON does not contain a stations[] mapping"
                )
            if requested_station not in series_by_station:
                known = ", ".join(sorted(series_by_station.keys()))
                raise ValueError(
                    f"Station '{requested_station}' not found in rainfall JSON stations[]; known stations: {known}"
                )
            json_series = series_by_station[requested_station]
        elif len(series_by_station) == 1:
            json_series = next(iter(series_by_station.values()))
        else:
            json_series = str(climate.get("series_name") or "").strip()
            if not json_series and len(series_by_station) > 1:
                known = ", ".join(sorted(series_by_station.keys()))
                raise ValueError(
                    "Rainfall JSON contains multiple stations. Provide --station-id or --series-name. "
                    f"Known stations: {known}"
                )

        if series_name is None:
            series_name = json_series
        elif json_series and json_series != series_name:
            raise ValueError(
                f"--series-name ({series_name}) does not match rainfall JSON series_name ({json_series})"
            )

    if not series_name:
        raise ValueError("A series name is required via --series-name or --rainfall-json")

    snippet = build_raingage_text(
        gage_id=args.gage_id,
        rain_format=args.rain_format,
        interval_min=args.interval_min,
        scf=args.scf,
        series_name=series_name,
    )
    write_text(args.out_text, snippet)

    payload = {
        "ok": True,
        "skill": "swmm-climate",
        "gage": {
            "id": args.gage_id,
            "rain_format": args.rain_format,
            "interval_min": args.interval_min,
            "scf": args.scf,
            "source": {
                "kind": "TIMESERIES",
                "series_name": series_name,
            },
        },
        "source_rainfall_json": str(args.rainfall_json) if args.rainfall_json is not None else None,
        "source_station_id": args.station_id,
        "outputs": {
            "text": str(args.out_text),
        },
    }
    write_json(args.out_json, payload)

    print(
        json.dumps(
            {
                "ok": True,
                "out_text": str(args.out_text),
                "out_json": str(args.out_json),
                "gage_id": args.gage_id,
                "series_name": series_name,
            },
            indent=2,
        )
    )


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Generate synthetic design-storm hyetographs (Chicago / alternating-block).

Output contract matches ``format_rainfall.py`` ``--out-json`` + ``--out-timeseries`` shape
so that ``build_swmm_inp.py --rainfall-json`` consumes the result unchanged.

Stdlib-only; zero ``agentic_swmm`` imports (portability constraint identical to the
other two scripts in this directory).

Determinism guarantee: no ``datetime.now()``, no ``random``, no network I/O.
Same args always produce byte-identical output files.
"""
from __future__ import annotations

import argparse
import csv
import json
import math
import sys
from pathlib import Path
from typing import Any

# ---------------------------------------------------------------------------
# IDF intensity helpers
# ---------------------------------------------------------------------------

def _idf_cn_form(t_min: float, *, A1: float, C: float, lgP: float, b: float, n: float) -> float:
    """Chicago/Chinese CN formula.

    q = 167 · A1 · (1 + C · lgP) / (t + b)^n   [L/s/ha]

    Convert to mm/hr:  mm/hr = q × 0.36
    (because L/s/ha × 3600 s/hr × 1 mm/1 L·m⁻² × 1/10 000 ha/m² … works out to ×0.36)
    """
    q_lsha = 167.0 * A1 * (1.0 + C * lgP) / ((t_min + b) ** n)
    return q_lsha * 0.36  # mm/hr


def _idf_generic_form(t_min: float, *, a: float, b: float, c: float) -> float:
    """Generic IDF formula.

    i = a / (t + b)^c   [mm/hr]
    """
    return a / ((t_min + b) ** c)


# ---------------------------------------------------------------------------
# Chicago / Keifer-Chu hyetograph
# ---------------------------------------------------------------------------

def chicago_hyetograph(
    *,
    coefficients: dict[str, float],
    form: str,
    return_period_yr: float,
    duration_min: float,
    dt_min: float,
    r: float,
) -> list[float]:
    """Return a list of depth increments (mm per timestep) for the Chicago hyetograph.

    The hyetograph is discretised at ``dt_min`` resolution; ``len(result) == n_steps``
    where ``n_steps = round(duration_min / dt_min)``.

    The peak block lands at index ``floor(r * n_steps)`` or ``floor(r * n_steps) + 1``,
    which is within one dt of ``r * duration_min`` (standard Chicago convention).

    **Mass conservation (Keifer-Chu defining property)**: the sum of all
    depth increments equals ``IDF_depth(duration_min)`` — the IDF cumulative
    depth for the design duration itself. The limbs carry r-weighted shares:
    the rising limb totals ``r * IDF_depth(T)`` and the falling limb
    ``(1 - r) * IDF_depth(T)``, because any window of duration tau centred
    on the peak (r*tau before, (1-r)*tau after) must accumulate exactly
    ``IDF_depth(tau)``. Limb cumulatives are therefore
    ``C_pre(s) = r * IDF_depth(s / r)`` and
    ``C_post(s) = (1 - r) * IDF_depth(s / (1 - r))``.

    Parameters
    ----------
    coefficients:
        For ``form="CN"``:  keys ``A1``, ``C``, ``b``, ``n``
        (``lgP`` is derived from ``return_period_yr``).
        For ``form="generic"``: keys ``a``, ``b``, ``c``.
    form:
        ``"CN"`` or ``"generic"``.
    return_period_yr:
        Return period in years (used only for CN form).
    duration_min:
        Total storm duration in minutes.
    dt_min:
        Timestep in minutes.
    r:
        Peak-position ratio (0 < r < 1, default 0.4).

    Returns
    -------
    list[float]
        Depth increments in mm per timestep, length == ``round(duration_min / dt_min)``.
    """
    n_steps = round(duration_min / dt_min)
    if n_steps < 1:
        raise ValueError(f"duration_min={duration_min} / dt_min={dt_min} must yield >= 1 step")

    form_upper = form.upper()
    lgP = math.log10(return_period_yr) if return_period_yr > 0 else 0.0

    def idf_depth(t_branch: float) -> float:
        """Cumulative IDF depth (mm) for branch duration t_branch (minutes).

        D(t) = i(t) * t / 60   where i(t) is average intensity over duration t [mm/hr].
        """
        if t_branch <= 0.0:
            return 0.0
        if form_upper == "CN":
            i_t = _idf_cn_form(
                t_branch,
                A1=coefficients["A1"],
                C=coefficients["C"],
                lgP=lgP,
                b=coefficients["b"],
                n=coefficients["n"],
            )
        elif form_upper == "GENERIC":
            i_t = _idf_generic_form(
                t_branch,
                a=coefficients["a"],
                b=coefficients["b"],
                c=coefficients["c"],
            )
        else:
            raise ValueError(f"Unknown IDF form '{form}'. Use 'CN' or 'generic'.")
        return i_t * t_branch / 60.0  # mm

    # Keifer-Chu limb cumulatives. A window of duration tau centred on the
    # peak (r*tau before, (1-r)*tau after) must accumulate exactly D(tau),
    # so the limb cumulative at distance s from the peak is the r-weighted
    # share of the window it closes: C_pre(s) = r * D(s/r) and
    # C_post(s) = (1-r) * D(s/(1-r)).
    def limb_pre(s: float) -> float:
        return r * idf_depth(s / r)

    def limb_post(s: float) -> float:
        return (1.0 - r) * idf_depth(s / (1.0 - r))

    # Blocks are differences of the storm cumulative F at bin edges:
    #   F(t) = C_pre(t_pre) - C_pre(t_pre - t)   for t <= t_pre
    #   F(t) = C_pre(t_pre) + C_post(t - t_pre)  for t >  t_pre
    # Mass is exact at any discretisation (the differences telescope to
    # F(T) = C_pre(rT) + C_post((1-r)T) = D(T)), and the block containing
    # the peak carries both its rising and falling slivers — no zero-depth
    # hole when r * duration lands exactly on a bin boundary.
    # Twin implementation: agentic_swmm/agent/swmm_runtime/design_storm.py
    # ``_chicago_from_idf`` (this script's portability constraint forbids
    # sharing code; tests/test_chicago_idf_parity.py locks the two equal).
    t_pre = r * duration_min

    def cumulative(t: float) -> float:
        if t <= t_pre:
            return limb_pre(t_pre) - limb_pre(t_pre - t)
        return limb_pre(t_pre) + limb_post(t - t_pre)

    depths: list[float] = []
    for k in range(n_steps):
        lo = k * dt_min
        # The last block ends at exactly duration_min so the storm total is
        # exactly D(duration_min) even when duration is not a multiple of dt.
        hi = duration_min if k == n_steps - 1 else (k + 1) * dt_min
        depths.append(max(0.0, cumulative(hi) - cumulative(lo)))

    return depths


# ---------------------------------------------------------------------------
# Alternating-block hyetograph
# ---------------------------------------------------------------------------

def alternating_block_hyetograph(
    *,
    idf_table: list[dict[str, float]],
    duration_min: float,
    dt_min: float,
) -> list[float]:
    """Return a list of depth increments (mm per timestep) using the alternating-block method.

    Parameters
    ----------
    idf_table:
        Sorted list of dicts with keys ``duration_min`` and ``intensity_mm_per_hr``.
        Must cover durations from dt_min up to duration_min at dt_min resolution
        (or be interpolable up to that resolution).
    duration_min:
        Total storm duration in minutes.
    dt_min:
        Timestep in minutes.

    Returns
    -------
    list[float]
        Depth increments in mm per timestep, with the largest block at the center.
    """
    n_steps = round(duration_min / dt_min)
    if n_steps < 1:
        raise ValueError(f"duration_min={duration_min} / dt_min={dt_min} must yield >= 1 step")

    # Sort IDF table by duration for interpolation
    sorted_table = sorted(idf_table, key=lambda r: r["duration_min"])
    durations_min = [r["duration_min"] for r in sorted_table]
    intensities = [r["intensity_mm_per_hr"] for r in sorted_table]

    def lookup_intensity(t: float) -> float:
        """Linearly interpolate (or extrapolate at ends) intensity for duration t."""
        if t <= durations_min[0]:
            return intensities[0]
        if t >= durations_min[-1]:
            return intensities[-1]
        for k in range(len(durations_min) - 1):
            if durations_min[k] <= t <= durations_min[k + 1]:
                frac = (t - durations_min[k]) / (durations_min[k + 1] - durations_min[k])
                return intensities[k] + frac * (intensities[k + 1] - intensities[k])
        return intensities[-1]

    # Compute incremental depths: delta_d[k] = D(k*dt) - D((k-1)*dt)
    # where D(t) = i(t) * t / 60  (mm, total depth over t minutes)
    incremental: list[float] = []
    prev_depth = 0.0
    for k in range(1, n_steps + 1):
        t = k * dt_min
        total_depth = lookup_intensity(t) * t / 60.0  # mm
        delta = total_depth - prev_depth
        incremental.append(max(0.0, delta))
        prev_depth = total_depth

    # Sort increments descending to assign to positions
    sorted_increments = sorted(incremental, reverse=True)

    # Alternating-block: place largest at center, alternate left/right
    depths: list[float] = [0.0] * n_steps
    center = n_steps // 2
    left = center - 1
    right = center + 1
    # Place largest at center
    depths[center] = sorted_increments[0]
    toggle = True  # True = place right next, False = place left
    for val in sorted_increments[1:]:
        if toggle and right < n_steps:
            depths[right] = val
            right += 1
            toggle = False
        elif not toggle and left >= 0:
            depths[left] = val
            left -= 1
            toggle = True
        elif right < n_steps:
            depths[right] = val
            right += 1
        elif left >= 0:
            depths[left] = val
            left -= 1

    return depths


# ---------------------------------------------------------------------------
# Output formatting helpers  (match format_rainfall.py contract)
# ---------------------------------------------------------------------------

def _format_number(value: float) -> str:
    """Format a float with up to 6 decimal places, stripping trailing zeros."""
    return f"{value:.6f}".rstrip("0").rstrip(".") or "0"


def build_timeseries_lines(series_name: str, depths_mm: list[float], dt_min: int) -> list[str]:
    """Render SWMM [TIMESERIES] body lines in the same style as format_rainfall.py.

    SWMM expects calendar-date rows:
        <name>  mm/dd/yyyy  HH:MM  value

    We anchor at a synthetic epoch: 01/01/2000.
    The value is depth per timestep in mm (VOLUME format in the gage section).
    The builder reads them as INTENSITY by default, but this script emits
    mm/hr so that the generated gage can use INTENSITY / the same unit as
    format_rainfall.py.

    Wait — format_rainfall.py emits mm/hr (intensity) values.  We must match.
    Convert depth_mm (per dt) back to mm/hr intensity for the TIMESERIES body.
    """
    lines = [";;Name             Date         Time       Value"]
    base_year = 2000
    base_month = 1
    base_day = 1
    total_minutes_offset = 0
    for depth_mm in depths_mm:
        # Convert depth (mm per dt_min) → intensity (mm/hr)
        intensity_mm_hr = (depth_mm / dt_min) * 60.0

        # Compute calendar time from offset
        total_minutes = total_minutes_offset
        days_offset = total_minutes // 1440
        remaining_minutes = total_minutes % 1440
        hh = remaining_minutes // 60
        mm_time = remaining_minutes % 60

        # Simple date arithmetic: roll day count from 01/01/2000
        # Using Gregorian day count for correctness
        _date = _minutes_to_date(base_year, base_month, base_day, total_minutes)
        date_str = f"{_date[1]:02d}/{_date[2]:02d}/{_date[0]}"
        time_str = f"{_date[3]:02d}:{_date[4]:02d}"
        lines.append(
            f"{series_name:<18} {date_str} {time_str} {_format_number(intensity_mm_hr)}"
        )
        total_minutes_offset += dt_min
    return lines


def _minutes_to_date(
    base_year: int, base_month: int, base_day: int, total_minutes: int
) -> tuple[int, int, int, int, int]:
    """Return (year, month, day, hour, minute) adding total_minutes to the base date."""
    days_add = total_minutes // 1440
    remaining = total_minutes % 1440
    hour = remaining // 60
    minute = remaining % 60

    # Day arithmetic using a simple Julian-day approach
    def is_leap(y: int) -> bool:
        return (y % 4 == 0 and y % 100 != 0) or (y % 400 == 0)

    def days_in_month(y: int, m: int) -> int:
        dom = [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
        if m == 2 and is_leap(y):
            return 29
        return dom[m]

    year = base_year
    month = base_month
    day = base_day + days_add

    # Normalise day overflow
    while True:
        dim = days_in_month(year, month)
        if day <= dim:
            break
        day -= dim
        month += 1
        if month > 12:
            month = 1
            year += 1

    return year, month, day, hour, minute


def build_out_json(
    *,
    ok: bool,
    out_json: str,
    out_timeseries: str,
    series_name: str,
    rows: int,
    interval_minutes: int,
    method: str,
    return_period_yr: float,
    coefficients: dict[str, float],
    form: str,
    duration_min: float,
    dt_min: float,
    r: float | None,
) -> dict[str, Any]:
    """Construct the metadata JSON matching format_rainfall.py's stdout contract.

    Mandatory keys (superset of format_rainfall.py's stdout):
      ok, out_json, out_timeseries, series_name, series_names, rows,
      stations, interval_minutes
    Plus new design-storm-specific keys:
      method, return_period_yr, coefficients, form, duration_min, dt_min, r
    """
    return {
        "ok": ok,
        "skill": "swmm-climate",
        "method": method,
        "form": form,
        "return_period_yr": return_period_yr,
        "coefficients": dict(coefficients),
        "duration_min": duration_min,
        "dt_min": dt_min,
        "r": r,
        "series_name": series_name,
        "series_names": [series_name],
        "rows": rows,
        "stations": 1,
        "interval_minutes": interval_minutes,
        "range": {
            "start": "2000-01-01T00:00",
            "end": None,
            "interval_minutes": interval_minutes,
        },
        "outputs": {
            "timeseries_text": out_timeseries,
        },
        "out_json": out_json,
        "out_timeseries": out_timeseries,
    }


# ---------------------------------------------------------------------------
# IDF table loading helpers
# ---------------------------------------------------------------------------

def load_idf_table_csv(path: Path) -> list[dict[str, float]]:
    """Load IDF table from CSV file with columns duration_min,intensity_mm_per_hr."""
    rows: list[dict[str, float]] = []
    with path.open("r", encoding="utf-8", newline="") as f:
        reader = csv.DictReader(f)
        for idx, row in enumerate(reader, start=2):
            try:
                dur = float(row["duration_min"])
                intensity = float(row["intensity_mm_per_hr"])
            except (KeyError, ValueError) as exc:
                raise ValueError(
                    f"IDF CSV row {idx}: expected columns 'duration_min' and 'intensity_mm_per_hr'. "
                    f"Error: {exc}"
                ) from exc
            rows.append({"duration_min": dur, "intensity_mm_per_hr": intensity})
    if not rows:
        raise ValueError(f"IDF table CSV is empty: {path}")
    return rows


def load_idf_table_json(raw: str) -> list[dict[str, float]]:
    """Load IDF table from inline JSON string (list of {duration_min, intensity_mm_per_hr})."""
    try:
        data = json.loads(raw)
    except json.JSONDecodeError as exc:
        raise ValueError(f"IDF table JSON parse error: {exc}") from exc
    if not isinstance(data, list):
        raise ValueError("IDF table JSON must be a list of objects")
    rows: list[dict[str, float]] = []
    for idx, item in enumerate(data):
        if not isinstance(item, dict):
            raise ValueError(f"IDF table JSON entry {idx} must be an object")
        try:
            dur = float(item["duration_min"])
            intensity = float(item["intensity_mm_per_hr"])
        except (KeyError, ValueError) as exc:
            raise ValueError(
                f"IDF table JSON entry {idx}: expected 'duration_min' and 'intensity_mm_per_hr'. "
                f"Error: {exc}"
            ) from exc
        rows.append({"duration_min": dur, "intensity_mm_per_hr": intensity})
    if not rows:
        raise ValueError("IDF table JSON list is empty")
    return rows


def write_json(path: Path, obj: Any) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(json.dumps(obj, indent=2), encoding="utf-8")


def write_text(path: Path, text: str) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(text, encoding="utf-8")


# ---------------------------------------------------------------------------
# CLI
# ---------------------------------------------------------------------------

def _build_cn_coefficients(args: argparse.Namespace) -> dict[str, float]:
    required = ["A1", "C", "b", "n"]
    missing = [k for k in required if getattr(args, k.lower(), None) is None]
    if missing:
        raise ValueError(f"CN form requires --A1, --C, --b, --n. Missing: {missing}")
    return {
        "A1": float(args.a1),
        "C": float(args.c),
        "b": float(args.b),
        "n": float(args.n),
    }


def _build_generic_coefficients(args: argparse.Namespace) -> dict[str, float]:
    required_attrs = [("a_coeff", "--a"), ("b", "--b"), ("c", "--c")]
    missing = []
    for attr, flag in required_attrs:
        if getattr(args, attr, None) is None:
            missing.append(flag)
    if missing:
        raise ValueError(f"generic form requires --a, --b, --c. Missing: {missing}")
    return {
        "a": float(args.a_coeff),
        "b": float(args.b),
        "c": float(args.c),
    }


def main(argv: list[str] | None = None) -> None:
    ap = argparse.ArgumentParser(
        description=(
            "Generate a synthetic design-storm hyetograph (Chicago or alternating-block) "
            "and write SWMM-compatible outputs matching format_rainfall.py's contract."
        )
    )

    # --- Method selection ---
    ap.add_argument(
        "--method",
        choices=["chicago", "alternating_block"],
        required=True,
        help="Hyetograph method: 'chicago' (Keifer-Chu) or 'alternating_block'.",
    )

    # --- IDF form (chicago only) ---
    ap.add_argument(
        "--form",
        choices=["CN", "generic"],
        default="generic",
        help=(
            "IDF formula form for --method chicago. "
            "'CN': q=167·A1·(1+C·lgP)/(t+b)^n [L/s/ha]; "
            "'generic': i=a/(t+b)^c [mm/hr]."
        ),
    )

    # --- Coefficient arguments ---
    ap.add_argument("--a1", type=float, default=None, help="CN form: coefficient A1.")
    ap.add_argument("--C", dest="c", type=float, default=None,
                    help="CN form: coefficient C. (lower-case dest avoids conflict)")
    ap.add_argument("--a-coeff", dest="a_coeff", type=float, default=None,
                    help="generic form: coefficient a.")
    ap.add_argument("--b", type=float, default=None, help="Both forms: coefficient b (time offset, min).")
    ap.add_argument("--n", type=float, default=None, help="CN form: exponent n.")
    # For generic form 'c' exponent - but --C is already used for CN 'C', so use a different name
    ap.add_argument("--c-exp", dest="c_exp", type=float, default=None,
                    help="generic form: exponent c in i=a/(t+b)^c.")

    # --- Alternating-block IDF table ---
    ap.add_argument(
        "--idf-csv",
        type=Path,
        default=None,
        help="CSV file with columns 'duration_min,intensity_mm_per_hr' for alternating-block method.",
    )
    ap.add_argument(
        "--idf-json",
        type=str,
        default=None,
        help=(
            "Inline JSON string (list of {duration_min, intensity_mm_per_hr}) "
            "for alternating-block method. Alternative to --idf-csv."
        ),
    )

    # --- Storm parameters ---
    ap.add_argument(
        "--return-period",
        type=float,
        default=2.0,
        help="Return period in years (default: 2).",
    )
    ap.add_argument(
        "--duration",
        type=float,
        required=True,
        help="Storm duration in minutes.",
    )
    ap.add_argument(
        "--dt",
        type=float,
        default=5.0,
        help="Timestep in minutes (default: 5).",
    )
    ap.add_argument(
        "--r",
        type=float,
        default=0.4,
        help="Peak-position ratio for Chicago method (default: 0.4). Ignored for alternating-block.",
    )

    # --- Output ---
    ap.add_argument("--out-json", type=Path, required=True, help="Output metadata JSON path.")
    ap.add_argument(
        "--out-timeseries", type=Path, required=True,
        help="Output text path for SWMM [TIMESERIES] body."
    )
    ap.add_argument(
        "--series-name",
        default=None,
        help=(
            "Series name token (default: TS_DESIGN_P<P>Y_<duration>MIN). "
            "Must not contain whitespace."
        ),
    )

    args = ap.parse_args(argv)

    # --- Validate duration / dt ---
    if args.duration <= 0:
        raise ValueError("--duration must be > 0")
    if args.dt <= 0:
        raise ValueError("--dt must be > 0")
    if args.return_period <= 0:
        raise ValueError("--return-period must be > 0")
    if not (0.0 < args.r < 1.0):
        raise ValueError("--r must be strictly between 0 and 1")

    n_steps = round(args.duration / args.dt)
    if n_steps < 1:
        raise ValueError(f"duration ({args.duration} min) / dt ({args.dt} min) < 1 step")

    # --- Derive series name ---
    P_str = f"{int(args.return_period)}" if args.return_period == int(args.return_period) else f"{args.return_period}"
    D_str = f"{int(args.duration)}" if args.duration == int(args.duration) else f"{args.duration}"
    default_series_name = f"TS_DESIGN_P{P_str}Y_{D_str}MIN"
    series_name = args.series_name if args.series_name is not None else default_series_name

    # --- Compute depths ---
    if args.method == "chicago":
        form = args.form.upper()
        if form == "CN":
            coefficients = _build_cn_coefficients(args)
        else:
            # Generic form: c exponent from --c-exp
            if args.c_exp is None:
                raise ValueError("generic form requires --c-exp (exponent c in i=a/(t+b)^c)")
            raw_coeff = {"a": args.a_coeff, "b": args.b, "c": args.c_exp}
            missing_g = [k for k, v in raw_coeff.items() if v is None]
            if missing_g:
                flag_map = {"a": "--a-coeff", "b": "--b", "c": "--c-exp"}
                raise ValueError(
                    f"generic form requires --a-coeff, --b, --c-exp. Missing: {[flag_map[k] for k in missing_g]}"
                )
            coefficients = {"a": float(args.a_coeff), "b": float(args.b), "c": float(args.c_exp)}

        depths = chicago_hyetograph(
            coefficients=coefficients,
            form=args.form,
            return_period_yr=args.return_period,
            duration_min=args.duration,
            dt_min=args.dt,
            r=args.r,
        )
    else:
        # Alternating block
        if args.idf_csv is not None and args.idf_json is not None:
            raise ValueError("Provide either --idf-csv or --idf-json, not both.")
        if args.idf_csv is not None:
            idf_table = load_idf_table_csv(args.idf_csv)
        elif args.idf_json is not None:
            idf_table = load_idf_table_json(args.idf_json)
        else:
            raise ValueError("--method alternating_block requires --idf-csv or --idf-json.")

        depths = alternating_block_hyetograph(
            idf_table=idf_table,
            duration_min=args.duration,
            dt_min=args.dt,
        )
        coefficients = {}
        form = "table"

    # --- Build timeseries text ---
    dt_int = int(round(args.dt))
    ts_lines = build_timeseries_lines(series_name, depths, dt_int)
    timeseries_text = "\n".join(ts_lines) + "\n"

    # --- Build out-json ---
    out_payload = build_out_json(
        ok=True,
        out_json=str(args.out_json),
        out_timeseries=str(args.out_timeseries),
        series_name=series_name,
        rows=n_steps,
        interval_minutes=dt_int,
        method=args.method,
        return_period_yr=args.return_period,
        coefficients=coefficients,
        form=args.form if args.method == "chicago" else "table",
        duration_min=args.duration,
        dt_min=args.dt,
        r=args.r if args.method == "chicago" else None,
    )

    # --- Write outputs ---
    write_text(args.out_timeseries, timeseries_text)
    write_json(args.out_json, out_payload)

    # stdout: compact summary matching format_rainfall.py's stdout shape
    print(
        json.dumps(
            {
                "ok": True,
                "out_json": str(args.out_json),
                "out_timeseries": str(args.out_timeseries),
                "series_name": series_name,
                "series_names": [series_name],
                "rows": n_steps,
                "stations": 1,
                "interval_minutes": dt_int,
            },
            indent=2,
        )
    )


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
from __future__ import annotations

import argparse
import csv
import glob
import hashlib
import json
import re
from collections import defaultdict
from dataclasses import dataclass
from datetime import datetime
from pathlib import Path
from typing import Any


@dataclass(frozen=True)
class RainRecord:
    station_id: str
    timestamp: datetime
    rainfall_mm_per_hr: float
    source_file: Path
    source_row: int


ALLOWED_UNITS_CANONICAL = ("mm_per_hr", "in_per_hr")
UNIT_ALIASES = {
    "mm_per_hr": "mm_per_hr",
    "mm/hr": "mm_per_hr",
    "mm/h": "mm_per_hr",
    "mmhr": "mm_per_hr",
    "in_per_hr": "in_per_hr",
    "in/hr": "in_per_hr",
    "in/h": "in_per_hr",
    "inhr": "in_per_hr",
}
UNIT_POLICY_CHOICES = ("strict", "convert_to_mm_per_hr")
TIMESTAMP_POLICY_CHOICES = ("strict", "sort")


def sha256_file(path: Path) -> str:
    digest = hashlib.sha256()
    with path.open("rb") as f:
        for chunk in iter(lambda: f.read(1024 * 1024), b""):
            digest.update(chunk)
    return digest.hexdigest()


def parse_timestamp(value: str, explicit_format: str) -> datetime:
    formats = [
        explicit_format,
        "%Y-%m-%d %H:%M:%S",
        "%Y-%m-%dT%H:%M",
        "%Y-%m-%dT%H:%M:%S",
        "%Y/%m/%d %H:%M",
        "%Y/%m/%d %H:%M:%S",
    ]
    for fmt in formats:
        try:
            return datetime.strptime(value, fmt)
        except ValueError:
            pass
    raise ValueError(f"Unsupported timestamp format: '{value}'")


def parse_window_timestamp(value: str | None, *, timestamp_format: str) -> datetime | None:
    if value is None:
        return None
    token = value.strip()
    if not token:
        return None
    return parse_timestamp(token, timestamp_format)


def normalize_units(value: str) -> str:
    token = value.strip().lower().replace(" ", "")
    if token not in UNIT_ALIASES:
        accepted = ", ".join(sorted(ALLOWED_UNITS_CANONICAL))
        raise ValueError(f"Unsupported --value-units '{value}'. Accepted canonical units: {accepted}")
    return UNIT_ALIASES[token]


def convert_to_mm_per_hr(value: float, *, units: str, policy: str) -> float:
    if policy not in UNIT_POLICY_CHOICES:
        raise ValueError(f"Unsupported unit policy: {policy}")
    if units not in ALLOWED_UNITS_CANONICAL:
        raise ValueError(f"Unsupported canonical units: {units}")

    if policy == "strict":
        if units != "mm_per_hr":
            raise ValueError(
                f"Unit policy 'strict' requires mm_per_hr input, got {units}. "
                "Use --unit-policy convert_to_mm_per_hr to convert supported non-SI units."
            )
        return value

    if units == "mm_per_hr":
        return value
    if units == "in_per_hr":
        return value * 25.4
    raise ValueError(f"Cannot convert units: {units}")


def format_location(rec: RainRecord) -> str:
    return f"{rec.source_file}:{rec.source_row}"


def sanitize_series_token(value: str) -> str:
    token = re.sub(r"[^A-Za-z0-9_]+", "_", value.strip())
    token = token.strip("_")
    return token or "STATION"


def derive_station_id_for_file(
    *,
    input_csv: Path,
    input_count: int,
    user_default_station_id: str | None,
) -> str:
    if input_count == 1:
        if user_default_station_id is not None and not user_default_station_id.strip():
            raise ValueError("--default-station-id cannot be blank")
        return user_default_station_id.strip() if user_default_station_id is not None else "STATION1"

    if user_default_station_id is not None:
        raise ValueError(
            "--default-station-id can only be used with a single --input when --station-column is omitted"
        )
    return input_csv.stem


def resolve_input_paths(*, explicit_inputs: list[Path], input_globs: list[str]) -> list[Path]:
    resolved: list[Path] = []
    for item in explicit_inputs:
        resolved.append(item)

    for pattern in input_globs:
        matches = sorted(glob.glob(pattern))
        if not matches:
            raise ValueError(f"--input-glob pattern matched no files: {pattern}")
        for matched in matches:
            resolved.append(Path(matched))

    if not resolved:
        raise ValueError("At least one input CSV is required via --input (and optionally --input-glob)")

    deduped: list[Path] = []
    seen_real: set[str] = set()
    for path in resolved:
        real = str(path.resolve())
        if real in seen_real:
            continue
        seen_real.add(real)
        if not path.exists():
            raise ValueError(f"Input file not found: {path}")
        if not path.is_file():
            raise ValueError(f"Input path is not a file: {path}")
        deduped.append(path)
    return deduped


DAT_UNIT_ALIASES = {
    "mm_per_hr": ("mm_per_hr", None),
    "mm/hr": ("mm_per_hr", None),
    "mm/h": ("mm_per_hr", None),
    "in_per_hr": ("in_per_hr", None),
    "in/hr": ("in_per_hr", None),
    "in/h": ("in_per_hr", None),
    "mm_per_day": ("mm_per_hr", 24.0),
    "mm/day": ("mm_per_hr", 24.0),
    "in_per_day": ("in_per_hr", 24.0),
    "in/day": ("in_per_hr", 24.0),
}


def parse_dat_value_units(raw: str) -> tuple[str, float | None]:
    token = raw.strip().lower().replace(" ", "")
    if token not in DAT_UNIT_ALIASES:
        raise ValueError(
            f"Unsupported --dat-value-units '{raw}'. Accepted: "
            f"{sorted(DAT_UNIT_ALIASES.keys())}"
        )
    return DAT_UNIT_ALIASES[token]


def read_records_from_dat(
    input_dat: Path,
    *,
    dat_value_units_raw: str,
    unit_policy: str,
) -> tuple[list[RainRecord], int]:
    canonical_units, interval_hours = parse_dat_value_units(dat_value_units_raw)
    records: list[RainRecord] = []
    interval_minutes_from_dat: int | None = (
        int(round(interval_hours * 60)) if interval_hours is not None else None
    )
    lines = input_dat.read_text(encoding="utf-8", errors="ignore").splitlines()
    for idx, raw_line in enumerate(lines, start=1):
        line = raw_line.strip()
        if not line or line.startswith(";"):
            continue
        parts = line.split()
        if len(parts) < 7:
            raise ValueError(
                f"SWMM .dat row at {input_dat}:{idx} expected 7 tokens "
                f"'<series> YYYY M D HH MM value', got {len(parts)}: {line}"
            )
        series_token, year, month, day, hour, minute, value = parts[:7]
        try:
            ts = datetime(int(year), int(month), int(day), int(hour), int(minute))
        except ValueError as exc:
            raise ValueError(f"Invalid date/time at {input_dat}:{idx}: {line}") from exc
        try:
            raw_value = float(value)
        except ValueError as exc:
            raise ValueError(f"Invalid value at {input_dat}:{idx}: {line}") from exc
        if raw_value < 0:
            raise ValueError(f"Rainfall value must be >= 0 at {input_dat}:{idx}")

        if interval_hours is not None:
            intensity = raw_value / interval_hours
        else:
            intensity = raw_value
        mm_per_hr = convert_to_mm_per_hr(intensity, units=canonical_units, policy=unit_policy)

        records.append(
            RainRecord(
                station_id=str(series_token),
                timestamp=ts,
                rainfall_mm_per_hr=mm_per_hr,
                source_file=input_dat,
                source_row=idx,
            )
        )
    if not records:
        raise ValueError(f"No data rows found in SWMM .dat: {input_dat}")
    if interval_minutes_from_dat is None:
        interval_minutes_from_dat = estimate_interval_minutes(records)
    return records, interval_minutes_from_dat or 0


def read_records_from_file(
    input_csv: Path,
    *,
    input_count: int,
    timestamp_column: str,
    value_column: str,
    station_column: str | None,
    default_station_id: str | None,
    input_units: str,
    unit_policy: str,
    timestamp_format: str,
) -> list[RainRecord]:
    with input_csv.open("r", encoding="utf-8", newline="") as f:
        reader = csv.DictReader(f)
        fieldnames = reader.fieldnames
        if fieldnames is None:
            raise ValueError(f"CSV is missing header row: {input_csv}")
        rows = list(reader)

    if not rows:
        raise ValueError(f"CSV has no rows: {input_csv}")

    if timestamp_column not in fieldnames:
        raise ValueError(f"Missing required column '{timestamp_column}' in {input_csv}")
    if value_column not in fieldnames:
        raise ValueError(f"Missing required column '{value_column}' in {input_csv}")
    if station_column is not None and station_column not in fieldnames:
        raise ValueError(f"Missing required station column '{station_column}' in {input_csv}")

    implicit_station_id = None
    if station_column is None:
        implicit_station_id = derive_station_id_for_file(
            input_csv=input_csv,
            input_count=input_count,
            user_default_station_id=default_station_id,
        )

    records: list[RainRecord] = []
    for idx, row in enumerate(rows, start=2):
        raw_ts = (row.get(timestamp_column) or "").strip()
        raw_value = (row.get(value_column) or "").strip()
        if not raw_ts:
            raise ValueError(f"Blank timestamp at {input_csv}:{idx}")
        if not raw_value:
            raise ValueError(f"Blank rainfall value at {input_csv}:{idx}")

        ts = parse_timestamp(raw_ts, timestamp_format)
        try:
            value = float(raw_value)
        except ValueError as exc:
            raise ValueError(f"Invalid rainfall value '{raw_value}' at {input_csv}:{idx}") from exc
        if value < 0:
            raise ValueError(f"Rainfall intensity must be >= 0 at {input_csv}:{idx}")

        if station_column is not None:
            station_raw = (row.get(station_column) or "").strip()
            if not station_raw:
                raise ValueError(f"Blank station id in column '{station_column}' at {input_csv}:{idx}")
            station_id = station_raw
        else:
            station_id = str(implicit_station_id)

        converted_value = convert_to_mm_per_hr(value, units=input_units, policy=unit_policy)
        records.append(
            RainRecord(
                station_id=station_id,
                timestamp=ts,
                rainfall_mm_per_hr=converted_value,
                source_file=input_csv,
                source_row=idx,
            )
        )

    return records


def format_number(value: float) -> str:
    return f"{value:.6f}".rstrip("0").rstrip(".") or "0"


def assign_series_names_by_station(
    *,
    station_ids: list[str],
    base_series_name: str,
    series_name_template: str | None,
) -> dict[str, str]:
    if not station_ids:
        raise ValueError("No stations available for series assignment")

    mapping: dict[str, str] = {}
    if len(station_ids) == 1 and series_name_template is None:
        station = station_ids[0]
        mapping[station] = base_series_name
        return mapping

    for station_id in station_ids:
        station_safe = sanitize_series_token(station_id)
        if series_name_template is None:
            series_name = f"{base_series_name}_{station_safe}"
        else:
            try:
                series_name = series_name_template.format(station=station_id, station_safe=station_safe)
            except KeyError as exc:
                raise ValueError(
                    f"--series-name-template contains unsupported placeholder '{exc.args[0]}'. "
                    "Allowed placeholders: {station}, {station_safe}"
                ) from exc
        series_name = series_name.strip()
        if not series_name:
            raise ValueError(f"Derived blank series name for station '{station_id}'")
        mapping[station_id] = series_name

    reverse: dict[str, str] = {}
    for station_id, series_name in mapping.items():
        if series_name in reverse:
            other_station = reverse[series_name]
            raise ValueError(
                f"Derived duplicate series name '{series_name}' for stations '{other_station}' and '{station_id}'. "
                "Adjust --series-name or --series-name-template."
            )
        reverse[series_name] = station_id
    return mapping


def render_timeseries_lines(
    *,
    station_order: list[str],
    series_by_station: dict[str, str],
    records_by_station: dict[str, list[RainRecord]],
) -> list[str]:
    lines = [
        ";;Name             Date         Time       Value",
    ]
    for station_id in station_order:
        series_name = series_by_station[station_id]
        for rec in records_by_station[station_id]:
            lines.append(
                # SWMM expects calendar dates in mm/dd/yyyy format in [TIMESERIES].
                f"{series_name:<18} {rec.timestamp.strftime('%m/%d/%Y')} {rec.timestamp.strftime('%H:%M')} {format_number(rec.rainfall_mm_per_hr)}"
            )
    return lines


def estimate_interval_minutes(records: list[RainRecord]) -> int | None:
    if len(records) < 2:
        return None
    deltas = [
        int((records[i + 1].timestamp - records[i].timestamp).total_seconds() / 60)
        for i in range(len(records) - 1)
    ]
    if any(d <= 0 for d in deltas):
        raise ValueError("Timestamps must be strictly increasing")
    if len(set(deltas)) == 1:
        return deltas[0]
    return None


def filter_records_by_window(
    records: list[RainRecord],
    *,
    window_start: datetime | None,
    window_end: datetime | None,
) -> list[RainRecord]:
    if window_start is None and window_end is None:
        return list(records)

    out: list[RainRecord] = []
    for rec in records:
        if window_start is not None and rec.timestamp < window_start:
            continue
        if window_end is not None and rec.timestamp > window_end:
            continue
        out.append(rec)
    return out


def validate_temporal_consistency(
    records: list[RainRecord],
    *,
    timestamp_policy: str,
) -> dict[str, bool]:
    if timestamp_policy not in TIMESTAMP_POLICY_CHOICES:
        raise ValueError(f"Unsupported --timestamp-policy: {timestamp_policy}")

    previous_by_station: dict[str, RainRecord] = {}
    seen_timestamps_by_station: dict[str, dict[datetime, RainRecord]] = defaultdict(dict)
    input_sorted_by_station: dict[str, bool] = {}

    for rec in records:
        station_id = rec.station_id
        input_sorted_by_station.setdefault(station_id, True)

        seen_for_station = seen_timestamps_by_station[station_id]
        prior_same_ts = seen_for_station.get(rec.timestamp)
        if prior_same_ts is not None:
            raise ValueError(
                f"Duplicate timestamp for station '{station_id}' at {format_location(rec)} and "
                f"{format_location(prior_same_ts)} ({rec.timestamp.isoformat(timespec='minutes')})"
            )
        seen_for_station[rec.timestamp] = rec

        prev = previous_by_station.get(station_id)
        if prev is not None and rec.timestamp <= prev.timestamp:
            input_sorted_by_station[station_id] = False
            if timestamp_policy == "strict":
                raise ValueError(
                    f"Non-monotonic timestamp for station '{station_id}': "
                    f"{format_location(prev)} has {prev.timestamp.isoformat(timespec='minutes')} and "
                    f"{format_location(rec)} has {rec.timestamp.isoformat(timespec='minutes')}. "
                    "Timestamps must be strictly increasing per station."
                )
        previous_by_station[station_id] = rec

    return input_sorted_by_station


def write_json(path: Path, obj: Any) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(json.dumps(obj, indent=2), encoding="utf-8")


def write_text(path: Path, text: str) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(text, encoding="utf-8")


def main() -> None:
    ap = argparse.ArgumentParser(
        description="Format rainfall CSV to SWMM-friendly timeseries text + JSON metadata (deterministic)."
    )
    ap.add_argument(
        "--input",
        type=Path,
        action="append",
        default=[],
        help="Input rainfall CSV. Repeat for batch mode.",
    )
    ap.add_argument(
        "--input-glob",
        action="append",
        default=[],
        help="Optional glob pattern for additional input files (e.g. 'data/rain_*.csv'). Repeat as needed.",
    )
    ap.add_argument(
        "--input-dat",
        type=Path,
        action="append",
        default=[],
        help=(
            "SWMM-style rainfall .dat file with rows '<series> YYYY M D HH MM value'. "
            "Use --dat-value-units to declare the per-row value units."
        ),
    )
    ap.add_argument(
        "--dat-value-units",
        default="mm_per_hr",
        help=(
            "Units of the value column in --input-dat rows. Accepts mm_per_hr, in_per_hr, "
            "or interval-volume aliases mm_per_day (24h volume) and in_per_day."
        ),
    )
    ap.add_argument("--out-json", type=Path, required=True, help="Output metadata JSON path.")
    ap.add_argument("--out-timeseries", type=Path, required=True, help="Output text path for SWMM [TIMESERIES] body.")
    ap.add_argument("--series-name", default="TS_RAIN")
    ap.add_argument(
        "--series-name-template",
        default=None,
        help="Optional template for per-station series names. Placeholders: {station}, {station_safe}.",
    )
    ap.add_argument("--timestamp-column", default="timestamp")
    ap.add_argument("--value-column", default="rainfall_mm_per_hr")
    ap.add_argument(
        "--station-column",
        default=None,
        help="Optional station/gage ID column. If omitted, single-file mode uses one station and multi-file mode derives station ids from file stems.",
    )
    ap.add_argument(
        "--default-station-id",
        default=None,
        help="Optional station id when --station-column is omitted and a single input file is provided.",
    )
    ap.add_argument("--timestamp-format", default="%Y-%m-%d %H:%M")
    ap.add_argument(
        "--window-start",
        default=None,
        help="Optional inclusive event window start timestamp.",
    )
    ap.add_argument(
        "--window-end",
        default=None,
        help="Optional inclusive event window end timestamp.",
    )
    ap.add_argument(
        "--value-units",
        default="mm_per_hr",
        help="Input rainfall units. Accepted canonical aliases: mm_per_hr/mm/hr, in_per_hr/in/hr.",
    )
    ap.add_argument(
        "--unit-policy",
        choices=list(UNIT_POLICY_CHOICES),
        default="strict",
        help="strict: only mm_per_hr accepted. convert_to_mm_per_hr: converts supported units to mm_per_hr.",
    )
    ap.add_argument(
        "--timestamp-policy",
        choices=list(TIMESTAMP_POLICY_CHOICES),
        default="strict",
        help="strict: reject non-monotonic timestamps per station. sort: allow and sort per station.",
    )
    args = ap.parse_args()

    use_dat_mode = bool(args.input_dat)
    if use_dat_mode and (args.input or args.input_glob):
        raise ValueError("--input-dat cannot be combined with --input or --input-glob")

    if use_dat_mode:
        input_paths = list(args.input_dat)
        for path in input_paths:
            if not path.exists():
                raise ValueError(f"Input file not found: {path}")
            if not path.is_file():
                raise ValueError(f"Input path is not a file: {path}")
        normalized_input_units = "mm_per_hr"
        dat_window_format = "%Y-%m-%d"
    else:
        input_paths = resolve_input_paths(explicit_inputs=args.input, input_globs=args.input_glob)
        normalized_input_units = normalize_units(args.value_units)
        dat_window_format = args.timestamp_format

    window_start = parse_window_timestamp(args.window_start, timestamp_format=dat_window_format)
    window_end = parse_window_timestamp(args.window_end, timestamp_format=dat_window_format)
    if window_start is not None and window_end is not None and window_start > window_end:
        raise ValueError("--window-start must be <= --window-end")

    all_records: list[RainRecord] = []
    dat_interval_minutes: int | None = None
    if use_dat_mode:
        if args.default_station_id is not None and len(input_paths) > 1:
            raise ValueError(
                "--default-station-id can only be used with a single --input-dat when "
                "--station-column is not used (multi .dat batch mode derives station ids "
                "from the series column)"
            )
        for input_dat in input_paths:
            file_records, interval_min = read_records_from_dat(
                input_dat=input_dat,
                dat_value_units_raw=args.dat_value_units,
                unit_policy="convert_to_mm_per_hr",
            )
            if len(input_paths) == 1 and args.default_station_id is not None:
                override = args.default_station_id.strip()
                if not override:
                    raise ValueError("--default-station-id cannot be blank")
                file_records = [
                    RainRecord(
                        station_id=override,
                        timestamp=rec.timestamp,
                        rainfall_mm_per_hr=rec.rainfall_mm_per_hr,
                        source_file=rec.source_file,
                        source_row=rec.source_row,
                    )
                    for rec in file_records
                ]
            all_records.extend(file_records)
            if interval_min:
                dat_interval_minutes = (
                    interval_min if dat_interval_minutes is None else min(dat_interval_minutes, interval_min)
                )
    else:
        for input_csv in input_paths:
            file_records = read_records_from_file(
                input_csv=input_csv,
                input_count=len(input_paths),
                timestamp_column=args.timestamp_column,
                value_column=args.value_column,
                station_column=args.station_column,
                default_station_id=args.default_station_id,
                input_units=normalized_input_units,
                unit_policy=args.unit_policy,
                timestamp_format=args.timestamp_format,
            )
            all_records.extend(file_records)

    rows_before_window = len(all_records)
    records = filter_records_by_window(
        all_records,
        window_start=window_start,
        window_end=window_end,
    )
    if not records:
        raise ValueError("No records remain after applying optional event window")

    input_sorted_by_station = validate_temporal_consistency(records, timestamp_policy=args.timestamp_policy)

    records_by_station: dict[str, list[RainRecord]] = defaultdict(list)
    for rec in records:
        records_by_station[rec.station_id].append(rec)

    for station_records in records_by_station.values():
        station_records.sort(key=lambda r: r.timestamp)

    station_ids = sorted(records_by_station.keys())
    series_by_station = assign_series_names_by_station(
        station_ids=station_ids,
        base_series_name=args.series_name,
        series_name_template=args.series_name_template,
    )

    station_payloads: list[dict[str, Any]] = []
    all_intervals: list[int | None] = []
    for station_id in station_ids:
        station_records = records_by_station[station_id]
        interval_minutes = estimate_interval_minutes(station_records)
        all_intervals.append(interval_minutes)
        source_files = sorted({str(rec.source_file) for rec in station_records})
        station_payloads.append(
            {
                "station_id": station_id,
                "series_name": series_by_station[station_id],
                "counts": {
                    "rows": len(station_records),
                    "input_sorted": input_sorted_by_station.get(station_id, True),
                },
                "range": {
                    "start": station_records[0].timestamp.isoformat(timespec="minutes"),
                    "end": station_records[-1].timestamp.isoformat(timespec="minutes"),
                    "interval_minutes": interval_minutes,
                },
                "source_files": source_files,
            }
        )

    all_sorted_records = sorted(records, key=lambda r: r.timestamp)
    interval_minutes_global: int | None = None
    if len(station_payloads) == 1:
        interval_minutes_global = all_intervals[0]

    timeseries_lines = render_timeseries_lines(
        station_order=station_ids,
        series_by_station=series_by_station,
        records_by_station=records_by_station,
    )
    timeseries_text = "\n".join(timeseries_lines) + "\n"

    input_sources = [{"path": str(path), "sha256": sha256_file(path)} for path in input_paths]
    series_name_legacy = station_payloads[0]["series_name"] if len(station_payloads) == 1 else None

    payload = {
        "ok": True,
        "skill": "swmm-climate",
        "input_csv": str(input_paths[0]) if len(input_paths) == 1 else None,
        "input_sha256": input_sources[0]["sha256"] if len(input_sources) == 1 else None,
        "inputs": input_sources,
        "schema": {
            "timestamp_column": args.timestamp_column,
            "value_column": args.value_column,
            "station_column": args.station_column,
            "value_units": "mm_per_hr",
            "input_value_units": normalized_input_units,
            "unit_policy": args.unit_policy,
        },
        "window": {
            "start": window_start.isoformat(timespec="minutes") if window_start is not None else None,
            "end": window_end.isoformat(timespec="minutes") if window_end is not None else None,
        },
        "series_name": series_name_legacy,
        "series_names": [station["series_name"] for station in station_payloads],
        "stations": station_payloads,
        "counts": {
            "rows": len(records),
            "rows_before_window": rows_before_window,
            "rows_after_window": len(records),
            "stations": len(station_payloads),
            "input_sorted": all(input_sorted_by_station.get(station_id, True) for station_id in station_ids),
            "input_sorted_by_station": {station_id: input_sorted_by_station.get(station_id, True) for station_id in station_ids},
        },
        "range": {
            "start": all_sorted_records[0].timestamp.isoformat(timespec="minutes"),
            "end": all_sorted_records[-1].timestamp.isoformat(timespec="minutes"),
            "interval_minutes": interval_minutes_global,
        },
        "outputs": {
            "timeseries_text": str(args.out_timeseries),
        },
    }

    write_text(args.out_timeseries, timeseries_text)
    write_json(args.out_json, payload)

    print(
        json.dumps(
            {
                "ok": True,
                "out_json": str(args.out_json),
                "out_timeseries": str(args.out_timeseries),
                "series_name": series_name_legacy,
                "series_names": payload["series_names"],
                "rows": len(records),
                "stations": len(station_payloads),
                "interval_minutes": interval_minutes_global,
            },
            indent=2,
        )
    )


if __name__ == "__main__":
    main()