geointel/scripts/audit_detection_false_negative_evidence.py

from __future__ import annotations

import argparse
import hashlib
import json
import math
import statistics
from datetime import datetime, timezone
from pathlib import Path
from typing import Any


AREA_BUCKETS = (
    ("tiny_lt_25_m2", 0.0, 25.0),
    ("small_25_100_m2", 25.0, 100.0),
    ("medium_100_500_m2", 100.0, 500.0),
    ("large_gte_500_m2", 500.0, math.inf),
)


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Audit persisted GeoIntel false-negative QA evidence across model portfolios."
    )
    parser.add_argument(
        "--portfolio",
        action="append",
        required=True,
        help="Label and calibration_evidence_portfolio.json path as label=/path/file.json.",
    )
    parser.add_argument("--output-dir", required=True, type=Path)
    return parser.parse_args()


def load_json(path: Path) -> dict[str, Any]:
    if not path.is_file():
        raise SystemExit(f"JSON input is not readable: {path}")
    payload = json.loads(path.read_text(encoding="utf-8-sig"))
    if not isinstance(payload, dict):
        raise SystemExit(f"JSON input must be an object: {path}")
    return payload


def parse_portfolio_arg(raw: str) -> tuple[str, Path]:
    label, separator, path_raw = raw.partition("=")
    label = label.strip()
    if not separator or not label or not path_raw.strip():
        raise SystemExit("--portfolio must use label=/path/to/calibration_evidence_portfolio.json")
    return label, Path(path_raw.strip()).expanduser().resolve()


def resolve_evidence_path(portfolio_path: Path, sample: dict[str, Any]) -> Path:
    raw = str(sample.get("evidence_geojson_path") or "")
    configured = Path(raw).expanduser()
    candidates = [configured]
    if raw and not configured.is_absolute():
        candidates.append(portfolio_path.parent / configured)
    sample_slug = str(sample.get("sample_slug") or "").strip().lower()
    candidates.append(
        portfolio_path.parent
        / "samples"
        / sample_slug
        / "evidence"
        / "calibration_evidence.geojson"
    )
    for candidate in candidates:
        if candidate.is_file():
            return candidate.resolve()
    raise SystemExit(f"Evidence GeoJSON is not readable for {sample_slug}: {raw}")


def percentile(values: list[float], fraction: float) -> float | None:
    if not values:
        return None
    ordered = sorted(values)
    position = (len(ordered) - 1) * fraction
    lower = math.floor(position)
    upper = math.ceil(position)
    if lower == upper:
        return ordered[lower]
    return ordered[lower] + (ordered[upper] - ordered[lower]) * (position - lower)


def area_stats(values: list[float]) -> dict[str, float | int | None]:
    return {
        "count": len(values),
        "min": min(values) if values else None,
        "p10": percentile(values, 0.10),
        "median": statistics.median(values) if values else None,
        "p90": percentile(values, 0.90),
        "max": max(values) if values else None,
    }


def area_bucket(area_m2: float) -> str:
    for label, minimum, maximum in AREA_BUCKETS:
        if minimum <= area_m2 < maximum:
            return label
    raise AssertionError(f"No area bucket for {area_m2}")


def stable_reference_id(feature: dict[str, Any], geometry: Any) -> str:
    properties = feature.get("properties") or {}
    source_feature_id = properties.get("source_feature_id")
    if source_feature_id not in (None, ""):
        return f"source:{source_feature_id}"
    reference_feature_id = properties.get("reference_feature_id")
    if reference_feature_id not in (None, ""):
        return f"reference:{reference_feature_id}"
    normalized_wkb = getattr(geometry.normalize(), "wkb", geometry.wkb)
    return f"geometry:{hashlib.sha256(normalized_wkb).hexdigest()}"


def audit_feature_collection(payload: dict[str, Any], geod: Any, shape: Any) -> dict[str, Any]:
    if payload.get("type") != "FeatureCollection":
        raise SystemExit("Evidence GeoJSON must be a FeatureCollection")
    features = payload.get("features") or []
    false_negative_ids: set[str] = set()
    reference_ids: set[str] = set()
    false_negative_areas: list[float] = []
    matched_reference_areas: list[float] = []
    reference_records: dict[str, dict[str, Any]] = {}
    bucket_counts = {
        label: {"false_negative": 0, "matched_reference": 0, "total_reference": 0, "false_negative_rate": None}
        for label, _, _ in AREA_BUCKETS
    }

    for feature in features:
        if not isinstance(feature, dict):
            raise SystemExit("Evidence GeoJSON contains a non-object feature")
        role = str((feature.get("properties") or {}).get("qa_evidence_role") or "")
        if role not in {"false_negative", "match_reference"}:
            continue
        geometry_payload = feature.get("geometry")
        if not isinstance(geometry_payload, dict):
            raise SystemExit(f"Evidence feature {feature.get('id')} has no geometry")
        geometry = shape(geometry_payload)
        if geometry.is_empty or not geometry.is_valid:
            raise SystemExit(f"Evidence feature {feature.get('id')} has invalid geometry")
        if geometry.geom_type not in {"Polygon", "MultiPolygon"}:
            raise SystemExit(
                f"Evidence feature {feature.get('id')} must be Polygon or MultiPolygon"
            )
        area_m2 = abs(float(geod.geometry_area_perimeter(geometry)[0]))
        bucket = area_bucket(area_m2)
        reference_id = stable_reference_id(feature, geometry)
        reference_ids.add(reference_id)
        reference_records[reference_id] = {
            "area_m2": area_m2,
            "area_bucket": bucket,
            "feature": feature,
        }
        bucket_role = "false_negative" if role == "false_negative" else "matched_reference"
        bucket_counts[bucket][bucket_role] += 1
        bucket_counts[bucket]["total_reference"] += 1
        if role == "false_negative":
            false_negative_ids.add(reference_id)
            false_negative_areas.append(area_m2)
        else:
            matched_reference_areas.append(area_m2)

    for values in bucket_counts.values():
        total = values["total_reference"]
        values["false_negative_rate"] = values["false_negative"] / total if total else None

    total_reference = len(false_negative_areas) + len(matched_reference_areas)
    return {
        "false_negative_ids": false_negative_ids,
        "reference_ids": reference_ids,
        "reference_records": reference_records,
        "false_negative_count": len(false_negative_areas),
        "matched_reference_count": len(matched_reference_areas),
        "total_reference_count": total_reference,
        "false_negative_rate": len(false_negative_areas) / total_reference if total_reference else None,
        "false_negative_area_m2": area_stats(false_negative_areas),
        "matched_reference_area_m2": area_stats(matched_reference_areas),
        "area_buckets": bucket_counts,
    }


def build_recommendations(samples: list[dict[str, Any]]) -> list[str]:
    recommendations: list[str] = []
    persistent_total = sum(sample["persistent_false_negative_count"] for sample in samples)
    if persistent_total:
        recommendations.append(
            f"Prioritize targeted positive sampling for {persistent_total} reference buildings missed by every compared portfolio."
        )
    ranked = sorted(
        samples,
        key=lambda sample: max(
            (item.get("false_negative_rate") or 0.0 for item in sample["portfolios"]),
            default=0.0,
        ),
        reverse=True,
    )
    if ranked:
        recommendations.append(
            "Expand or oversample the weakest AOIs first: "
            + ", ".join(sample["sample_slug"] for sample in ranked[:3])
            + "."
        )
    small_bias_samples: list[str] = []
    for sample in samples:
        for portfolio in sample["portfolios"]:
            buckets = portfolio["area_buckets"]
            small_total = sum(
                buckets[key]["total_reference"]
                for key in ("tiny_lt_25_m2", "small_25_100_m2")
            )
            small_fn = sum(
                buckets[key]["false_negative"]
                for key in ("tiny_lt_25_m2", "small_25_100_m2")
            )
            larger_total = sum(
                buckets[key]["total_reference"]
                for key in ("medium_100_500_m2", "large_gte_500_m2")
            )
            larger_fn = sum(
                buckets[key]["false_negative"]
                for key in ("medium_100_500_m2", "large_gte_500_m2")
            )
            small_rate = small_fn / small_total if small_total else 0.0
            larger_rate = larger_fn / larger_total if larger_total else 0.0
            if small_total >= 5 and small_rate >= larger_rate + 0.15:
                small_bias_samples.append(sample["sample_slug"])
                break
    if small_bias_samples:
        recommendations.append(
            "Review small-building label retention and pixel-size gates for: "
            + ", ".join(sorted(set(small_bias_samples)))
            + "."
        )
    if not recommendations:
        recommendations.append(
            "No dominant false-negative pattern was found; preserve current data gates and collect more independent AOIs."
        )
    return recommendations


def run_audit(portfolio_args: list[str], output_dir: Path) -> tuple[Path, Path]:
    try:
        from pyproj import Geod
        from shapely.geometry import shape
    except ImportError as exc:
        raise SystemExit(
            "False-negative GIS audit requires the existing GeoIntel GIS extras (pyproj and shapely)."
        ) from exc

    parsed_portfolios = [parse_portfolio_arg(raw) for raw in portfolio_args]
    labels = [label for label, _ in parsed_portfolios]
    if len(labels) != len(set(labels)):
        raise SystemExit("Every --portfolio label must be unique")
    geod = Geod(ellps="WGS84")
    portfolio_samples: dict[str, dict[str, dict[str, Any]]] = {}
    portfolio_meta: list[dict[str, Any]] = []

    for label, portfolio_path in parsed_portfolios:
        portfolio = load_json(portfolio_path)
        samples = portfolio.get("samples") or []
        if not isinstance(samples, list) or not samples:
            raise SystemExit(f"Evidence portfolio has no samples: {portfolio_path}")
        sample_results: dict[str, dict[str, Any]] = {}
        for sample in samples:
            sample_slug = str(sample.get("sample_slug") or "").strip().lower()
            if not sample_slug:
                raise SystemExit(f"Evidence portfolio has a sample without sample_slug: {portfolio_path}")
            evidence_path = resolve_evidence_path(portfolio_path, sample)
            result = audit_feature_collection(load_json(evidence_path), geod, shape)
            result.update(
                {
                    "label": label,
                    "model_asset_id": portfolio.get("model_asset_id"),
                    "evidence_geojson_path": str(evidence_path),
                }
            )
            sample_results[sample_slug] = result
        portfolio_samples[label] = sample_results
        portfolio_meta.append(
            {
                "label": label,
                "path": str(portfolio_path),
                "model_asset_id": portfolio.get("model_asset_id"),
                "sample_count": len(sample_results),
            }
        )

    expected_slugs = set(next(iter(portfolio_samples.values())))
    for label, samples in portfolio_samples.items():
        if set(samples) != expected_slugs:
            raise SystemExit(
                f"Portfolio {label} has different AOIs; fixed-threshold comparisons require identical samples"
            )

    sample_reports: list[dict[str, Any]] = []
    persistent_evidence_features: list[dict[str, Any]] = []
    for sample_slug in sorted(expected_slugs):
        portfolio_rows = []
        false_negative_sets = []
        reference_sets = []
        for label, _ in parsed_portfolios:
            raw = portfolio_samples[label][sample_slug]
            false_negative_sets.append(raw["false_negative_ids"])
            reference_sets.append(raw["reference_ids"])
            portfolio_rows.append(
                {
                    key: value
                    for key, value in raw.items()
                    if key not in {"false_negative_ids", "reference_ids", "reference_records"}
                }
            )
        if any(reference_ids != reference_sets[0] for reference_ids in reference_sets[1:]):
            counts = ", ".join(
                f"{label}={len(reference_ids)}"
                for (label, _), reference_ids in zip(parsed_portfolios, reference_sets, strict=True)
            )
            raise SystemExit(
                f"Sample {sample_slug} has different reference populations across portfolios ({counts})"
            )
        persistent_ids = sorted(set.intersection(*false_negative_sets))
        reference_records = portfolio_samples[parsed_portfolios[0][0]][sample_slug]["reference_records"]
        persistent_areas = [float(reference_records[reference_id]["area_m2"]) for reference_id in persistent_ids]
        persistent_bucket_counts = {
            label: {"count": 0, "share": 0.0}
            for label, _, _ in AREA_BUCKETS
        }
        for reference_id in persistent_ids:
            record = reference_records[reference_id]
            persistent_bucket_counts[record["area_bucket"]]["count"] += 1
            source_feature = record["feature"]
            properties = dict(source_feature.get("properties") or {})
            properties.update(
                {
                    "qa_evidence_role": "persistent_false_negative",
                    "sample_slug": sample_slug,
                    "persistent_reference_id": reference_id,
                    "area_m2": record["area_m2"],
                    "area_bucket": record["area_bucket"],
                    "compared_portfolios": labels,
                }
            )
            persistent_evidence_features.append(
                {
                    "type": "Feature",
                    "id": f"{sample_slug}:{reference_id}",
                    "properties": properties,
                    "geometry": source_feature["geometry"],
                }
            )
        if persistent_ids:
            for values in persistent_bucket_counts.values():
                values["share"] = values["count"] / len(persistent_ids)
        sample_reports.append(
            {
                "sample_slug": sample_slug,
                "reference_population_count": len(reference_sets[0]),
                "persistent_false_negative_count": len(persistent_ids),
                "persistent_reference_ids": persistent_ids,
                "persistent_false_negative_area_m2": area_stats(persistent_areas),
                "persistent_area_buckets": persistent_bucket_counts,
                "portfolios": portfolio_rows,
            }
        )

    output_dir = output_dir.expanduser().resolve()
    output_dir.mkdir(parents=True, exist_ok=True)
    persistent_geojson_path = output_dir / "persistent_false_negatives.geojson"
    persistent_geojson_path.write_text(
        json.dumps(
            {"type": "FeatureCollection", "features": persistent_evidence_features},
            indent=2,
            sort_keys=True,
        ),
        encoding="utf-8",
    )
    report = {
        "generated_at": datetime.now(timezone.utc).isoformat(),
        "schema_version": 1,
        "input_crs": "EPSG:4326",
        "area_method": "WGS84 geodesic area via pyproj.Geod",
        "portfolio_count": len(portfolio_meta),
        "portfolios": portfolio_meta,
        "sample_count": len(sample_reports),
        "samples": sample_reports,
        "persistent_evidence_geojson_path": str(persistent_geojson_path),
        "recommendations": build_recommendations(sample_reports),
    }
    json_path = output_dir / "detection_false_negative_audit.json"
    json_path.write_text(json.dumps(report, indent=2, sort_keys=True), encoding="utf-8")

    lines = [
        "# Detection false-negative evidence audit",
        "",
        f"- Generated: {report['generated_at']}",
        f"- Portfolios: {report['portfolio_count']}",
        f"- AOIs: {report['sample_count']}",
        f"- Area method: {report['area_method']}",
        "",
        "## AOI comparison",
        "",
        "| AOI | Persistent misses | Persistent tiny/small | Persistent median m2 | "
        + " | ".join(f"{label} FN rate" for label, _ in parsed_portfolios)
        + " |",
        "|---|---:|---:|---:|" + "---:|" * len(parsed_portfolios),
    ]
    for sample in sample_reports:
        rates = [
            f"{(row['false_negative_rate'] or 0.0):.3f}"
            for row in sample["portfolios"]
        ]
        persistent_buckets = sample["persistent_area_buckets"]
        persistent_small_count = sum(
            persistent_buckets[key]["count"]
            for key in ("tiny_lt_25_m2", "small_25_100_m2")
        )
        persistent_median = sample["persistent_false_negative_area_m2"]["median"]
        persistent_median_text = f"{persistent_median:.1f}" if persistent_median is not None else "n/a"
        lines.append(
            f"| {sample['sample_slug']} | {sample['persistent_false_negative_count']} | "
            f"{persistent_small_count} | {persistent_median_text} | "
            + " | ".join(rates)
            + " |"
        )
    lines.extend(["", "## Recommended data actions", ""])
    lines.extend(f"- {item}" for item in report["recommendations"])
    markdown_path = output_dir / "detection_false_negative_audit.md"
    markdown_path.write_text("\n".join(lines) + "\n", encoding="utf-8")
    return json_path, markdown_path


def main() -> None:
    args = parse_args()
    json_path, markdown_path = run_audit(args.portfolio, args.output_dir)
    print(f"False-negative audit JSON: {json_path}")
    print(f"False-negative audit Markdown: {markdown_path}")


if __name__ == "__main__":
    main()