"""Offline, deterministic descriptive analysis of Lombardia reported cases.
Run: python analysis.py --input-dir . --output-dir rerun_outputs
No network, installation, model fitting, or random resampling occurs.
"""
from pathlib import Path
import argparse
import csv
import hashlib
import json
import platform
import sys
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

SCOPE = {
    "region": "Lombardia", "region_code": 3,
    "start": "2020-07-01", "end": "2020-12-31",
    "context_start": "2020-06-21", "context_end": "2021-01-10",
    "windows": [7, 21], "primary_window": 7,
    "commit": "fc82987b9ebda01fe8e5b59d5667e3e00721edb9",
    "max_download_bytes": 100000000,
    "missing_policy": "Full daily windows only; no imputation",
}
WEEKDAYS = ["Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"]

def write_csv(frame, path, index=False):
    frame.to_csv(path, index=index, float_format="%.12g", date_format="%Y-%m-%d",
                 lineterminator="\n")

def verify_inputs(root):
    manifest = pd.read_csv(root / "provenance/source_manifest.csv")
    assert manifest["commit"].eq(SCOPE["commit"]).all(), "Mixed source revisions"
    assert manifest["bytes"].sum() < SCOPE["max_download_bytes"], "Download cap exceeded"
    for item in manifest.itertuples(index=False):
        path = root / item.path
        payload = path.read_bytes()
        assert len(payload) == item.bytes, f"Byte-count mismatch: {item.path}"
        assert hashlib.sha256(payload).hexdigest() == item.sha256, f"Checksum mismatch: {item.path}"
    return manifest

def analyze(root, out, tables_only=False):
    root, out = Path(root).resolve(), Path(out).resolve()
    for protected in (root / "inputs", root / "provenance"):
        if out == protected or out.is_relative_to(protected):
            raise ValueError("Output directory must not overwrite frozen inputs or provenance")
    out.mkdir(parents=True, exist_ok=True)
    manifest = verify_inputs(root)
    pieces, schemas = [], set()
    for day in pd.date_range(SCOPE["context_start"], SCOPE["context_end"]):
        filename = f"dpc-covid19-ita-regioni-{day:%Y%m%d}.csv"
        frame = pd.read_csv(root / "inputs/dati-regioni" / filename)
        schemas.add(tuple(frame.columns))
        row = frame.loc[pd.to_numeric(frame["codice_regione"]) == SCOPE["region_code"]].copy()
        assert len(row) == 1, f"Region duplicate/missing: {filename}"
        assert row["denominazione_regione"].eq(SCOPE["region"]).all()
        row["date"] = pd.to_datetime(row["data"].str[:10])
        assert row["date"].eq(day).all(), f"Date mismatch: {filename}"
        row["source_file"] = filename
        pieces.append(row)
    full = pd.concat(pieces, ignore_index=True).sort_values("date").set_index("date")
    expected = pd.date_range(SCOPE["context_start"], SCOPE["context_end"], name="date")
    assert full.index.equals(expected) and not full.index.duplicated().any()
    y = pd.to_numeric(full["nuovi_positivi"])
    full["cumulative_difference"] = full["totale_casi"].diff()
    full["cumulative_discrepancy"] = y - full["cumulative_difference"]
    full["tests_daily"] = full["tamponi"].diff()
    full["tests_mean7"] = full["tests_daily"].rolling(7, center=True, min_periods=7).mean()
    full["tests_to_mean7"] = full["tests_daily"] / full["tests_mean7"].where(full["tests_mean7"] > 0)
    for window in SCOPE["windows"]:
        full[f"mean{window}"] = y.rolling(window, center=True, min_periods=window).mean()
        h = window // 2
        direct = np.array([sum(float(v) for v in y.iloc[i-h:i+h+1]) / window
                           if y.iloc[i-h:i+h+1].notna().all() else np.nan
                           for i in range(h, len(y)-h)])
        np.testing.assert_allclose(full[f"mean{window}"].iloc[h:-h], direct,
                                   rtol=1e-12, atol=1e-10, equal_nan=True)
    full["raw_minus_mean7"] = y - full["mean7"]
    full["raw_to_mean7"] = y / full["mean7"].where(full["mean7"] > 0)
    full["relative_deviation_pct"] = 100 * (full["raw_to_mean7"] - 1)
    full["raw_daily_change"] = y.diff()
    full["mean7_daily_change"] = full["mean7"].diff()
    full["weekday"] = full.index.day_name()
    target = full.loc[SCOPE["start"]:SCOPE["end"]].copy()
    assert len(target) == 184
    assert target[["mean7", "mean21"]].notna().all().all(), "Incomplete target smoothing windows"
    write_csv(full.reset_index(), out / "lombardia_with_context.csv")
    write_csv(target.reset_index(), out / "daily_calculations.csv")

    quality_items = {
        "target_expected_dates": 184,
        "target_observed_dates": len(target),
        "context_observed_dates": len(full),
        "duplicate_dates": int(full.index.duplicated().sum()),
        "target_missing_new_cases": int(target["nuovi_positivi"].isna().sum()),
        "target_negative_new_cases": int((target["nuovi_positivi"] < 0).sum()),
        "target_cumulative_discrepancies": int((target["cumulative_discrepancy"].abs() > 1e-9).sum()),
        "target_missing_cumulative_tests": int(target["tamponi"].isna().sum()),
        "target_missing_people_tested": int(target["casi_testati"].isna().sum()),
        "target_negative_test_differences": int((target["tests_daily"] < 0).sum()),
        "target_missing_mean7": int(target["mean7"].isna().sum()),
        "target_missing_mean21": int(target["mean21"].isna().sum()),
        "daily_csv_header_variants": len(schemas),
        "source_file_count": len(manifest),
        "downloaded_source_bytes": int(manifest["bytes"].sum()),
    }
    for col in ("note", "note_test", "note_casi"):
        quality_items[f"target_nonempty_{col}"] = int(target[col].fillna("").astype(str).str.strip().ne("").sum())
    notes = pd.read_csv(root / "inputs/note/dpc-covid19-ita-note.csv")
    note_dates = pd.to_datetime(notes["data"].str[:10])
    period_notes = notes.loc[note_dates.between(pd.Timestamp(SCOPE["context_start"]),
                                               pd.Timestamp(SCOPE["context_end"]))].copy()
    region_pattern = r"\b(?:lombard\w*|brescia|bergamo|milano|varese|como|cremona|mantova|monza|pavia|sondrio|lecco|lodi)\b"
    region_notes = period_notes.loc[period_notes["note"].str.contains(region_pattern, case=False, na=False)]
    quality_items["context_global_note_rows"] = len(period_notes)
    quality_items["context_regional_keyword_note_rows"] = len(region_notes)
    write_csv(period_notes, out / "official_notes_in_context.csv")
    write_csv(region_notes, out / "official_regional_keyword_notes.csv")
    write_csv(pd.DataFrame(quality_items.items(), columns=["check", "value"]), out / "quality_checks.csv")
    write_csv(target.loc[(target["cumulative_discrepancy"].abs() > 1e-9) |
                         (target["nuovi_positivi"] < 0) | (target["tests_daily"] < 0)].reset_index(),
              out / "flagged_rows.csv")

    week_start = target.index - pd.to_timedelta(target.index.dayofweek, unit="D")
    weekly_all = target.groupby(week_start).agg(
        observed_days=("nuovi_positivi", "count"),
        cases_total=("nuovi_positivi", "sum"),
        tests_total=("tests_daily", "sum"),
    )
    weekly_all.index.name = "week_start"
    weekly_all["week_end"] = weekly_all.index + pd.Timedelta(days=6)
    weekly_all["complete_week"] = weekly_all["observed_days"].eq(7)
    write_csv(weekly_all.reset_index(), out / "week_coverage.csv")
    weekly = weekly_all.loc[weekly_all["complete_week"]].copy()
    weekly["change_from_previous_week_pct"] = 100 * weekly["cases_total"].pct_change(fill_method=None)
    weekly["tests_change_from_previous_week_pct"] = 100 * weekly["tests_total"].pct_change(fill_method=None)
    for start, row in weekly.iterrows():
        vals = target.loc[start:start+pd.Timedelta(days=6), "nuovi_positivi"]
        assert len(vals) == 7 and row["cases_total"] == sum(int(v) for v in vals)
    assert weekly.index.to_series().diff().dropna().eq(pd.Timedelta(days=7)).all()
    write_csv(weekly.reset_index(), out / "weekly_summaries.csv")

    weekday_rows = []
    for day in WEEKDAYS:
        subset = target.loc[target["weekday"].eq(day)]
        ratios = subset["raw_to_mean7"].dropna()
        tests = subset["tests_to_mean7"].dropna()
        weekday_rows.append({
            "weekday": day, "n": len(ratios),
            "ratio_q25": ratios.quantile(.25), "ratio_median": ratios.median(),
            "ratio_q75": ratios.quantile(.75),
            "deviation_pct_q25": 100*(ratios.quantile(.25)-1),
            "deviation_pct_median": 100*(ratios.median()-1),
            "deviation_pct_q75": 100*(ratios.quantile(.75)-1),
            "tests_ratio_median": tests.median(), "tests_n": len(tests),
        })
    weekday = pd.DataFrame(weekday_rows)
    write_csv(weekday, out / "weekday_summaries.csv")

    peak_rows = []
    for column in ("nuovi_positivi", "mean7", "mean21"):
        values = target[column]
        peak = float(values.max())
        maxima = target.index[np.isclose(values, peak, rtol=0, atol=1e-10)]
        near = target.index[values >= .95*peak]
        peak_rows.append({
            "series": column, "maximum_cases_per_day": peak,
            "maximum_dates": ";".join(maxima.strftime("%Y-%m-%d")),
            "dates_at_least_95pct_maximum": ";".join(near.strftime("%Y-%m-%d")),
            "number_dates_at_least_95pct_maximum": len(near),
            "first_date_at_least_95pct_maximum": near.min().strftime("%Y-%m-%d"),
            "last_date_at_least_95pct_maximum": near.max().strftime("%Y-%m-%d"),
        })
    peaks = pd.DataFrame(peak_rows)
    write_csv(peaks, out / "smoothing_sensitivity.csv")

    summary = {
        "scope": SCOPE, "quality": quality_items,
        "total_reported_cases_in_window": int(target["nuovi_positivi"].sum()),
        "daily_cases_mean": float(target["nuovi_positivi"].mean()),
        "daily_cases_median": float(target["nuovi_positivi"].median()),
        "daily_cases_sd": float(target["nuovi_positivi"].std(ddof=1)),
        "daily_cases_min": int(target["nuovi_positivi"].min()),
        "daily_cases_max": int(target["nuovi_positivi"].max()),
        "median_absolute_raw_daily_change": float(target["raw_daily_change"].abs().median()),
        "median_absolute_mean7_daily_change": float(target["mean7_daily_change"].abs().median()),
        "median_absolute_deviation_pct": float(target["relative_deviation_pct"].abs().median()),
        "relative_deviation_pct_q25": float(target["relative_deviation_pct"].quantile(.25)),
        "relative_deviation_pct_q75": float(target["relative_deviation_pct"].quantile(.75)),
        "complete_week_count": len(weekly),
        "peak_week_start": weekly["cases_total"].idxmax().strftime("%Y-%m-%d"),
        "peak_week_cases": int(weekly["cases_total"].max()),
        "last_complete_week_start": weekly.index[-1].strftime("%Y-%m-%d"),
        "last_complete_week_cases": int(weekly["cases_total"].iloc[-1]),
        "last_week_vs_peak_pct": float(100*(weekly["cases_total"].iloc[-1]/weekly["cases_total"].max()-1)),
        "peaks": peak_rows,
        "validation": "All source checksums, date coverage, direct-window arithmetic and weekly sums passed.",
    }
    (out / "summary.json").write_text(json.dumps(summary, indent=2, ensure_ascii=False)+"\n", encoding="utf-8")
    versions = {"python": sys.version, "platform": platform.platform(),
                "numpy": np.__version__, "pandas": pd.__version__, "matplotlib": matplotlib.__version__}
    (out / "runtime_versions.json").write_text(json.dumps(versions, indent=2)+"\n", encoding="utf-8")
    if not tables_only:
        render_figure(target, weekday, out)
    return summary

def render_figure(target, weekday, out):
    with plt.rc_context({"font.size": 11, "axes.spines.top": False, "axes.spines.right": False,
                         "pdf.fonttype": 42, "savefig.facecolor": "white"}):
        fig = plt.figure(figsize=(12, 11), layout="constrained")
        grid = fig.add_gridspec(3, 1, height_ratios=[2.0, 1.05, 1.4])
        ax = fig.add_subplot(grid[0])
        dev = fig.add_subplot(grid[1], sharex=ax)
        dow = fig.add_subplot(grid[2])
        fig.suptitle("Lombardia: reported cases and reporting variation\n1 July–31 December 2020", fontsize=17)
        ax.plot(target.index, target["nuovi_positivi"], color="#777777", lw=.8,
                marker=".", markersize=2, label="Daily reported new cases")
        ax.plot(target.index, target["mean7"], color="#0072B2", lw=2.3,
                label="Centered 7-day mean (primary)")
        ax.plot(target.index, target["mean21"], color="#B45A00", lw=2.0, ls="--",
                label="Centered 21-day mean (sensitivity)")
        ax.set(title="A  The broad autumn wave persists after weekly smoothing",
               ylabel="Reported cases / day", ylim=(0, None))
        ax.legend(loc="upper left", frameon=False, fontsize=10)
        ax.tick_params(labelbottom=False)
        ax.grid(axis="y", alpha=.18)
        dev.axhline(0, color="#333333", lw=.8)
        dev.plot(target.index, target["relative_deviation_pct"], color="#555555",
                 lw=.7, marker=".", markersize=2.5)
        dev.set(title="B  Daily departures from the local 7-day mean",
                ylabel="Departure (%)")
        dev.xaxis.set_major_locator(mdates.MonthLocator())
        dev.xaxis.set_major_formatter(mdates.DateFormatter("%b"))
        dev.set_xlim(pd.Timestamp(SCOPE["start"]), pd.Timestamp(SCOPE["end"]))
        dev.grid(axis="y", alpha=.18)
        groups = [target.loc[target["weekday"].eq(day), "relative_deviation_pct"].dropna().values
                  for day in WEEKDAYS]
        dow.boxplot(groups, positions=np.arange(1, 8), widths=.5, patch_artist=True,
                    showfliers=True, whis=1.5,
                    boxprops={"facecolor":"#C7DDEC", "edgecolor":"#333333"},
                    medianprops={"color":"#111111", "linewidth":1.7},
                    whiskerprops={"color":"#444444"}, capprops={"color":"#444444"},
                    flierprops={"marker":"o", "markersize":3, "markeredgecolor":"#555555"})
        dow.axhline(0, color="#333333", lw=.8, ls="--")
        dow.set_xticks(np.arange(1, 8),
                       [f"{day[:3]}\nn={len(group)}" for day, group in zip(WEEKDAYS, groups)])
        dow.set(title="C  Weekday distributions of daily departures",
                ylabel="Departure (%)")
        dow.grid(axis="y", alpha=.18)
        fig.supxlabel("Departures = 100 × (daily cases / centered 7-day mean − 1).\n"
                      "Boxes: median and IQR; whiskers: ≤1.5 IQR; circles: outside whiskers. No confidence bands.\n"
                      "Source: Italian Civil Protection, pinned fc82987b9ebd; boundary data: 21 Jun–10 Jan.",
                      fontsize=9)
        fig.savefig(out / "comparison.png", dpi=220)
        fig.savefig(out / "comparison.pdf")
        plt.close(fig)

def main(argv=None):
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--input-dir", type=Path, default=Path(__file__).resolve().parent)
    parser.add_argument("--output-dir", type=Path, default=None)
    parser.add_argument("--tables-only", action="store_true")
    args = parser.parse_args(argv)
    out = args.output_dir or args.input_dir / "rerun_outputs"
    result = analyze(args.input_dir, out, tables_only=args.tables_only)
    print(json.dumps({"status":"verified", "total_cases":result["total_reported_cases_in_window"],
                      "complete_weeks":result["complete_week_count"],
                      "output_directory":str(out)}, indent=2))
    return result

if __name__ == "__main__":
    main()
