import asyncio
import ctypes
import json
import logging
import os
import re
import tempfile
from concurrent.futures import ProcessPoolExecutor

from fastapi import FastAPI, File, Form, HTTPException, Request, UploadFile
from fastapi.middleware.cors import CORSMiddleware
from fastapi.responses import JSONResponse, Response

from .batch import BatchError
from .columnar import to_columnar
from .gisio import read_vector
from . import dem as dem_store
from . import gisops
from .drawing_dxf import drawing_to_dxf
from urllib.parse import quote
from .jobs import get_job
from .epanet import EN, EpanetAPI, EpanetException, FLOW_UNITS, PRESS_UNITS, QUAL_TYPES

logger = logging.getLogger("spipe")

# Limits (overridable per deployment). The upload limit must match nginx's
# client_max_body_size; the period/cell caps bound CPU time and response size
# for inputs like DURATION 10000 hrs + HYDRAULIC TIMESTEP 1 sec.
MAX_UPLOAD_BYTES = int(float(os.environ.get("SPIPE_MAX_UPLOAD_MB", "20")) * 1024 * 1024)
# DEM rasters are much larger than INP/GIS files; streamed to disk, never held in memory.
MAX_DEM_BYTES = int(float(os.environ.get("SPIPE_MAX_DEM_MB", "300")) * 1024 * 1024)
ANALYZE_TIMEOUT_S = float(os.environ.get("SPIPE_ANALYZE_TIMEOUT_S", "120"))
MAX_PERIODS = int(os.environ.get("SPIPE_MAX_PERIODS", "10000"))
MAX_RESULT_CELLS = int(os.environ.get("SPIPE_MAX_RESULT_CELLS", "50000000"))
# Batch jobs (criticality, optimisation ...) solve the network many times.
JOB_TIMEOUT_S = float(os.environ.get("SPIPE_JOB_TIMEOUT_S", "900"))
WORKERS = int(os.environ.get("SPIPE_WORKERS", str(max(1, min(4, (os.cpu_count() or 2) - 1)))))

app = FastAPI()

app.add_middleware(
    CORSMiddleware,
    allow_origins=["http://localhost:5173", "http://127.0.0.1:5173"],
    allow_methods=["*"],
    allow_headers=["*"],
)

api = EpanetAPI()

# Analyses run in separate processes: the ctypes loop is CPU-bound and
# synchronous, so running it on the event loop froze every other request
# (health checks waited up to 16 s behind one large network).
_pool = None


def _get_pool() -> ProcessPoolExecutor:
    global _pool
    if _pool is None:
        _pool = ProcessPoolExecutor(max_workers=WORKERS)
    return _pool


@app.middleware("http")
async def reject_oversized_body(request: Request, call_next):
    # Starlette spools the whole multipart body before the handler runs, so
    # check the declared size first instead of after receiving it all.
    # scope["path"] is the routed path; request.url is rebuilt from the Host
    # header, so a malformed Host could make the path check miss (GHSA-86qp).
    length = request.headers.get("content-length")
    path = request.scope.get("path", "")
    if (path == "/api/analyze" or path.startswith(("/api/jobs/", "/api/gis/"))) and length and length.isdigit() and int(length) > MAX_UPLOAD_BYTES + 64 * 1024:
        return JSONResponse(status_code=413, content={"detail": f"File size exceeds {MAX_UPLOAD_BYTES // (1024 * 1024)} MiB limit"})
    return await call_next(request)


def run_job_bytes(name: str, contents: bytes, params: dict) -> dict:
    """Process-pool entry point for /api/jobs/{name}; same envelope as analyze_bytes."""
    try:
        result = get_job(name)(contents, params)
        return {"json": json.dumps(result, ensure_ascii=False, allow_nan=False, separators=(",", ":")).encode("utf-8")}
    except BatchError as e:
        return {"error": (400, str(e))}
    except EpanetException as e:
        return {"error": (422, str(e))}
    except ValueError as e:
        return {"error": (500, f"Result contains non-finite values: {e}")}


@app.post("/api/jobs/{name}")
async def run_job(name: str, file: UploadFile = File(...), params: str = Form("{}")):
    """Batch analysis over one INP. `params` is a JSON object (see backend/jobs.py)."""
    if not (file.filename or "").lower().endswith(".inp"):
        raise HTTPException(status_code=400, detail="File must be an .inp file")
    try:
        parsed = json.loads(params or "{}")
        if not isinstance(parsed, dict):
            raise ValueError
    except ValueError:
        raise HTTPException(status_code=400, detail="params must be a JSON object")
    contents = await file.read(MAX_UPLOAD_BYTES + 1)
    if len(contents) > MAX_UPLOAD_BYTES:
        raise HTTPException(status_code=413, detail=f"File size exceeds {MAX_UPLOAD_BYTES // (1024 * 1024)} MiB limit")
    loop = asyncio.get_running_loop()
    try:
        outcome = await asyncio.wait_for(
            loop.run_in_executor(_get_pool(), run_job_bytes, name, contents, parsed),
            timeout=JOB_TIMEOUT_S,
        )
    except asyncio.TimeoutError:
        raise HTTPException(status_code=504, detail=f"Analysis exceeded {JOB_TIMEOUT_S:.0f} s time limit")
    if "error" in outcome:
        status, detail = outcome["error"]
        raise HTTPException(status_code=status, detail=detail)
    return Response(content=outcome["json"], media_type="application/json")


def read_vector_bytes(contents: bytes, filename: str, layer: str | None, src_epsg: int | None, default_crs: str | None = None) -> dict:
    try:
        result = read_vector(contents, filename, layer, src_epsg, default_crs)
        return {"json": json.dumps(result, ensure_ascii=False, allow_nan=False, separators=(",", ":")).encode("utf-8")}
    except BatchError as e:
        return {"error": (400, str(e))}


@app.post("/api/gis/vector")
async def gis_vector(file: UploadFile = File(...), layer: str = Form(""), srcEpsg: str = Form(""), defaultCrs: str = Form("")):
    """Any GDAL vector format (zipped SHP, DXF, DGN, MIF, GML, KML, GPKG, E00 ...) -> GeoJSON (WGS84)."""
    src_epsg = None
    if srcEpsg.strip():
        if not srcEpsg.strip().isdigit():
            raise HTTPException(status_code=400, detail="srcEpsg must be an EPSG number")
        src_epsg = int(srcEpsg)
    contents = await file.read(MAX_UPLOAD_BYTES + 1)
    if len(contents) > MAX_UPLOAD_BYTES:
        raise HTTPException(status_code=413, detail=f"File size exceeds {MAX_UPLOAD_BYTES // (1024 * 1024)} MiB limit")
    loop = asyncio.get_running_loop()
    try:
        outcome = await asyncio.wait_for(
            loop.run_in_executor(_get_pool(), read_vector_bytes, contents, file.filename or "", layer or None, src_epsg, defaultCrs or None),
            timeout=ANALYZE_TIMEOUT_S,
        )
    except asyncio.TimeoutError:
        raise HTTPException(status_code=504, detail="GIS conversion timed out")
    if "error" in outcome:
        status, detail = outcome["error"]
        raise HTTPException(status_code=status, detail=detail)
    return Response(content=outcome["json"], media_type="application/json")


def _dem_call(fn, *args):
    """Run a DEM operation; BatchError -> 400 (pickle-safe tuple for the pool)."""
    try:
        return {"ok": fn(*args)}
    except dem_store.BatchError as e:
        return {"error": (400, str(e))}


async def _dem_run(fn, *args, timeout: float = 300):
    loop = asyncio.get_running_loop()
    try:
        outcome = await asyncio.wait_for(loop.run_in_executor(None, _dem_call, fn, *args), timeout=timeout)
    except asyncio.TimeoutError:
        raise HTTPException(status_code=504, detail="DEM processing timed out")
    if "error" in outcome:
        status, detail = outcome["error"]
        raise HTTPException(status_code=status, detail=detail)
    return outcome["ok"]


@app.post("/api/dem")
async def dem_upload(file: UploadFile = File(...)):
    """Store a DEM raster (GeoTIFF / ASC / IMG, or a zip holding one) and return its metadata."""
    fd, tmp = tempfile.mkstemp(prefix="spipe_dem_")
    size = 0
    try:
        with os.fdopen(fd, "wb") as f:
            while chunk := await file.read(1024 * 1024):
                size += len(chunk)
                if size > MAX_DEM_BYTES:
                    raise HTTPException(status_code=413, detail=f"DEM exceeds {MAX_DEM_BYTES // (1024 * 1024)} MiB limit")
                f.write(chunk)
        return await _dem_run(dem_store.import_dem, tmp, file.filename or "dem.tif")
    finally:
        try:
            os.remove(tmp)
        except OSError:
            pass


@app.get("/api/dem-base")
async def dem_base():
    """Server-provided DEMs (registered with `python -m backend.dem register ...`)."""
    return await _dem_run(dem_store.list_base)


@app.get("/api/dem/{dem_id}")
async def dem_get(dem_id: str):
    return await _dem_run(dem_store.dem_info, dem_id)


@app.get("/api/dem/{dem_id}/preview.png")
async def dem_preview(dem_id: str):
    png = await _dem_run(dem_store.preview_png, dem_id)
    return Response(content=png, media_type="image/png", headers={"Cache-Control": "private, max-age=86400"})


@app.get("/api/dem/{dem_id}/terrain/{z}/{x}/{y}.png")
async def dem_terrain(dem_id: str, z: int, x: int, y: int):
    """Terrarium-encoded terrain tiles for the 3D view."""
    png = await _dem_run(dem_store.terrain_tile, dem_id, z, x, y, timeout=60)
    return Response(content=png, media_type="image/png", headers={"Cache-Control": "private, max-age=86400"})


@app.post("/api/dem/{dem_id}/sample")
async def dem_sample(dem_id: str, request: Request):
    """Body: {"points": [[x, y], ...], "crs": "<EPSG or proj4 of the points>"}."""
    body = await request.json()
    points = body.get("points") or []
    if len(points) > 500_000:
        raise HTTPException(status_code=400, detail="too many points")
    return {"values": await _dem_run(dem_store.sample, dem_id, points, body.get("crs"))}


@app.post("/api/dem/{dem_id}/hydro")
async def dem_hydro(dem_id: str, request: Request):
    """Body: {"op": "streams", "minAreaKm2": 0.1} or {"op": "watershed", "points": [[x, y]], "crs": "...", "snap": m}."""
    body = await request.json()
    if body.get("op") == "streams":
        return await _dem_run(dem_store.streams, dem_id, float(body.get("minAreaKm2") or 0.1), body.get("bbox"), body.get("crs"), timeout=600)
    if body.get("op") == "watershed":
        return await _dem_run(dem_store.watershed, dem_id, body.get("points") or [], body.get("crs"), float(body.get("snap") or 50), body.get("bbox"), timeout=600)
    raise HTTPException(status_code=400, detail="op must be streams or watershed")


@app.post("/api/dem/from-points")
async def dem_from_points(request: Request):
    """Body: {"points": [[x, y, z], ...], "crs": "...", "cell": metres, "name": "..."} -> new DEM."""
    body = await request.json()
    points = body.get("points") or []
    if len(points) > 2_000_000:
        raise HTTPException(status_code=400, detail="too many points")
    return await _dem_run(dem_store.from_points, points, body.get("crs"), float(body.get("cell") or 0), str(body.get("name") or "점 표고 DEM"), body.get("targetCrs"))


@app.post("/api/gis/export")
async def gis_export(request: Request):
    """Body: {"featureCollection": {...WGS84}, "crs": "<target>", "format": "shp|gpkg|dxf|geojson|kml", "name": "..."}."""
    body = await request.json()
    data, filename, media = await _dem_run(gisops.write_vector, body.get("featureCollection"), body.get("crs"), str(body.get("format") or "geojson"), str(body.get("name") or "layer"))
    return Response(content=data, media_type=media, headers={"Content-Disposition": f"attachment; filename*=UTF-8''{quote(filename)}"})


@app.post("/api/gis/op")
async def gis_op(request: Request):
    """Body: {"op": "buffer|dissolve|intersection|clip|difference|voronoi|grid", "crs": "<planar>", ...}."""
    body = await request.json()
    return await _dem_run(gisops.vector_op, str(body.get("op") or ""), body)


@app.post("/api/gis/reproject")
async def gis_reproject(request: Request):
    """Body: {"featureCollection": {...}, "from": "<CRS it was read as>", "to": "<actual CRS>"}."""
    body = await request.json()
    return await _dem_run(gisops.reproject_fc, body.get("featureCollection"), body.get("from"), body.get("to"))


_UESC = re.compile(r"\\U\+([0-9A-Fa-f]{4})")
_CRLF = "\r\n"


def dxf_to_cp949(text: str) -> bytes:
    """R12 DXF from the browser (non-ASCII written as backslash-U+XXXX) -> Korean
    code page (ANSI_949), which Korean CAD reads; characters outside CP949 stay escaped."""
    def unesc(m: re.Match) -> str:
        ch = chr(int(m.group(1), 16))
        try:
            ch.encode("cp949")
            return ch
        except UnicodeEncodeError:
            return m.group(0)

    out = _UESC.sub(unesc, text)
    ver = _CRLF.join(["$ACADVER", "1", "AC1009", ""])
    out = out.replace(ver, ver + _CRLF.join(["9", "$DWGCODEPAGE", "3", "ANSI_949", ""]), 1)
    return out.encode("cp949", errors="replace")


@app.post("/api/drawing/dxf")
async def drawing_dxf(request: Request):
    """Body: {"drawing": {...items}, "name": ...} -> DXF R2018 (ezdxf, UTF-8 text);
    legacy {"dxf": "<R12 DXF text>"} -> CP949 R12."""
    body = await request.json()
    name = str(body.get("name") or "drawing") + ".dxf"
    if isinstance(body.get("drawing"), dict):
        data = await _dem_run(drawing_to_dxf, body["drawing"], timeout=120)
        return Response(content=data, media_type="application/dxf",
                        headers={"Content-Disposition": f"attachment; filename*=UTF-8''{quote(name)}"})
    text = str(body.get("dxf") or "")
    if not text.startswith("0" + _CRLF + "SECTION") or len(text) > 50_000_000:
        raise HTTPException(status_code=400, detail="DXF가 올바르지 않습니다.")
    name = str(body.get("name") or "drawing") + ".dxf"
    return Response(content=dxf_to_cp949(text), media_type="application/dxf",
                    headers={"Content-Disposition": f"attachment; filename*=UTF-8''{quote(name)}"})


@app.get("/api/health")
def health():
    return {"status": "ok"}


def analyze_bytes(contents: bytes, columnar: bool = False) -> dict:
    """Process-pool entry point. Returns {"json": bytes} or {"error": (status, detail)}
    because HTTPException does not survive pickling back to the parent."""
    try:
        with tempfile.TemporaryDirectory() as tmpdir:
            inp_path = os.path.join(tmpdir, "network.inp")
            rpt_path = os.path.join(tmpdir, "network.rpt")
            out_path = os.path.join(tmpdir, "network.out")
            with open(inp_path, "wb") as f:
                f.write(contents)
            result = run_epanet(inp_path, rpt_path, out_path)
        if columnar:
            result = to_columnar(result)
        # Serialise here, in the worker: encoding a 10 MB result on the event
        # loop stalled other requests by ~0.5 s per response.
        body = json.dumps(result, ensure_ascii=False, allow_nan=False, separators=(",", ":")).encode("utf-8")
        return {"json": body}
    except HTTPException as e:
        return {"error": (e.status_code, e.detail)}
    except ValueError as e:  # NaN/Inf in results
        return {"error": (500, f"Result contains non-finite values: {e}")}


@app.post("/api/analyze")
async def analyze(file: UploadFile = File(...), format: str = "object"):
    """Full EPS analysis. ``format=columnar`` returns the compact packed time
    series (backend/columnar.py); the default keeps the per-element objects."""
    if not (file.filename or "").lower().endswith(".inp"):
        raise HTTPException(status_code=400, detail="File must be an .inp file")

    contents = await file.read(MAX_UPLOAD_BYTES + 1)
    if len(contents) > MAX_UPLOAD_BYTES:
        raise HTTPException(status_code=413, detail=f"File size exceeds {MAX_UPLOAD_BYTES // (1024 * 1024)} MiB limit")

    loop = asyncio.get_running_loop()
    try:
        outcome = await asyncio.wait_for(
            loop.run_in_executor(_get_pool(), analyze_bytes, contents, format == "columnar"),
            timeout=ANALYZE_TIMEOUT_S,
        )
    except asyncio.TimeoutError:
        raise HTTPException(status_code=504, detail=f"Analysis exceeded {ANALYZE_TIMEOUT_S:.0f} s time limit")
    if "error" in outcome:
        status, detail = outcome["error"]
        raise HTTPException(status_code=status, detail=detail)
    return Response(content=outcome["json"], media_type="application/json")


def decode_id(raw_bytes: bytes) -> str:
    try:
        return raw_bytes.decode('utf-8')
    except UnicodeDecodeError:
        try:
            return raw_bytes.decode('cp949')
        except UnicodeDecodeError:
            return raw_bytes.decode('latin1', errors='replace')

def read_report_errors(rpt_path: str, limit: int = 10) -> list:
    """Return the 'Error NNN: ...' lines EPANET wrote to the report file."""
    if not os.path.exists(rpt_path):
        return []
    try:
        with open(rpt_path, "r", encoding="utf-8", errors="replace") as rf:
            lines = [ln.strip() for ln in rf if re.match(r"\s*Error \d+", ln)]
    except OSError:
        return []
    return lines[:limit]

def run_epanet(inp_path: str, rpt_path: str, out_path: str):
    ph = api.EN_Project()
    
    # 1. Version check
    version = ctypes.c_int()
    api.lib.EN_getversion(ctypes.byref(version))
    if version.value != 20305:
        raise HTTPException(status_code=500, detail=f"Unsupported EPANET version: {version.value}, expected 20305")
        
    ph_created = False
    ph_opened = False
    hydraulics_opened = False
        
    try:
        api.lib.EN_createproject(ctypes.byref(ph))
        ph_created = True
        
        # open -- a failed EN_open still leaves the input/report files open, so
        # EN_close must run in the finally block either way (EN_runproject does
        # the same); otherwise Windows can't delete the temp dir afterwards.
        ph_opened = True
        api.lib.EN_open(ph, inp_path.encode('utf-8'), rpt_path.encode('utf-8'), out_path.encode('utf-8'))
        
        # get units
        flow_unit = ctypes.c_int()
        api.lib.EN_getflowunits(ph, ctypes.byref(flow_unit))
        f_unit_str = FLOW_UNITS.get(flow_unit.value, "Unknown")
        
        press_unit = ctypes.c_double()
        api.lib.EN_getoption(ph, EN.PRESS_UNITS, ctypes.byref(press_unit))
        p_unit_str = PRESS_UNITS.get(int(press_unit.value), "Unknown")
        
        # get counts
        node_count = ctypes.c_int()
        api.lib.EN_getcount(ph, EN.NODECOUNT, ctypes.byref(node_count))
        
        link_count = ctypes.c_int()
        api.lib.EN_getcount(ph, EN.LINKCOUNT, ctypes.byref(link_count))

        # Reject resource-exhausting inputs before solving (a few-KB INP with
        # DURATION 10000 hrs / HYDRAULIC TIMESTEP 1 sec asks for 36M periods).
        duration, hstep = ctypes.c_long(), ctypes.c_long()
        api.lib.EN_gettimeparam(ph, EN.DURATION, ctypes.byref(duration))
        api.lib.EN_gettimeparam(ph, EN.HYDSTEP, ctypes.byref(hstep))
        report_step, report_start = ctypes.c_long(), ctypes.c_long()
        api.lib.EN_gettimeparam(ph, EN.REPORTSTEP, ctypes.byref(report_step))
        api.lib.EN_gettimeparam(ph, EN.REPORTSTART, ctypes.byref(report_start))
        planned_periods = duration.value // max(1, hstep.value) + 1
        planned_cells = planned_periods * (node_count.value + link_count.value)
        if planned_periods > MAX_PERIODS or planned_cells > MAX_RESULT_CELLS:
            raise HTTPException(
                status_code=413,
                detail=(
                    f"Simulation too large: {planned_periods} hydraulic periods x "
                    f"{node_count.value + link_count.value} elements "
                    f"(limits {MAX_PERIODS} periods, {MAX_RESULT_CELLS} result values). "
                    "Increase HYDRAULIC TIMESTEP or shorten DURATION."
                ),
            )
        
        # extract node metadata
        nodes = {}
        node_id_map = {}
        for i in range(1, node_count.value + 1):
            node_id = ctypes.create_string_buffer(256)
            api.lib.EN_getnodeid(ph, i, node_id)
            nid = decode_id(node_id.value)
            node_id_map[i] = nid
            
            node_type = ctypes.c_int()
            api.lib.EN_getnodetype(ph, i, ctypes.byref(node_type))
            ntype = node_type.value
            
            x, y = ctypes.c_double(), ctypes.c_double()
            coords = None
            try:
                api.lib.EN_getcoord(ph, i, ctypes.byref(x), ctypes.byref(y))
                coords = {"x": x.value, "y": y.value}
            except EpanetException:
                pass
                
            nodes[nid] = {
                "id": nid,
                "type": ntype,
                "coordinates": coords,
            }
            
        # extract link metadata
        links = {}
        link_id_map = {}
        for i in range(1, link_count.value + 1):
            link_id = ctypes.create_string_buffer(256)
            api.lib.EN_getlinkid(ph, i, link_id)
            lid = decode_id(link_id.value)
            link_id_map[i] = lid
            
            link_type = ctypes.c_int()
            api.lib.EN_getlinktype(ph, i, ctypes.byref(link_type))
            ltype = link_type.value
            
            n1, n2 = ctypes.c_int(), ctypes.c_int()
            api.lib.EN_getlinknodes(ph, i, ctypes.byref(n1), ctypes.byref(n2))
            
            n1_id = ctypes.create_string_buffer(256)
            api.lib.EN_getnodeid(ph, n1.value, n1_id)
            
            n2_id = ctypes.create_string_buffer(256)
            api.lib.EN_getnodeid(ph, n2.value, n2_id)
            
            v_count = ctypes.c_int()
            api.lib.EN_getvertexcount(ph, i, ctypes.byref(v_count))
            
            vertices = []
            for j in range(1, v_count.value + 1):
                vx, vy = ctypes.c_double(), ctypes.c_double()
                api.lib.EN_getvertex(ph, i, j, ctypes.byref(vx), ctypes.byref(vy))
                vertices.append({"x": vx.value, "y": vy.value})
            
            links[lid] = {
                "id": lid,
                "type": ltype,
                "fromNode": decode_id(n1_id.value),
                "toNode": decode_id(n2_id.value),
                "vertices": vertices,
            }

        # 3. Open hydraulics with EN_SAVE=1 (Save to hyd file so water quality can use it)
        api.lib.EN_openH(ph)
        hydraulics_opened = True
        api.lib.EN_initH(ph, EN.SAVE)

        time_steps = []
        node_series = {nid: {"demand": [], "head": [], "pressure": [], "quality": []} for nid in nodes}
        link_series = {lid: {"flow": [], "velocity": [], "headloss": [], "quality": [], "energy": [], "efficiency": []} for lid in links}

        t = ctypes.c_long()
        tstep = ctypes.c_long(1)

        engine_warnings = []
        while True:
            # Warning codes 1-6 (e.g. 1 unbalanced, 6 negative pressures) are
            # returned here, not raised; they only reach the report file when
            # [REPORT] MESSAGES is on, so collect them directly.
            code = api.lib.EN_runH(ph, ctypes.byref(t))
            if 0 < code <= 10:
                engine_warnings.append({"time": t.value, "code": code})
            time_steps.append(t.value)
            # Control/tank events add steps beyond DURATION/HYDSTEP; keep the cap hard.
            if len(time_steps) > MAX_PERIODS * 2:
                raise HTTPException(status_code=413, detail=f"Simulation produced more than {MAX_PERIODS * 2} hydraulic periods")

            for i in range(1, node_count.value + 1):
                nid = node_id_map[i]
                d, h, p = ctypes.c_double(), ctypes.c_double(), ctypes.c_double()
                api.lib.EN_getnodevalue(ph, i, EN.DEMAND, ctypes.byref(d))
                api.lib.EN_getnodevalue(ph, i, EN.HEAD, ctypes.byref(h))
                api.lib.EN_getnodevalue(ph, i, EN.PRESSURE, ctypes.byref(p))
                node_series[nid]["demand"].append(d.value)
                node_series[nid]["head"].append(h.value)
                node_series[nid]["pressure"].append(p.value)

            for i in range(1, link_count.value + 1):
                lid = link_id_map[i]
                f, v, hl = ctypes.c_double(), ctypes.c_double(), ctypes.c_double()
                api.lib.EN_getlinkvalue(ph, i, EN.FLOW, ctypes.byref(f))
                api.lib.EN_getlinkvalue(ph, i, EN.VELOCITY, ctypes.byref(v))
                api.lib.EN_getlinkvalue(ph, i, EN.HEADLOSS, ctypes.byref(hl))
                link_series[lid]["flow"].append(f.value)
                link_series[lid]["velocity"].append(v.value)
                link_series[lid]["headloss"].append(hl.value)

                if links[lid]["type"] == EN.PUMP:
                    e_val, ef_val = ctypes.c_double(), ctypes.c_double()
                    api.lib.EN_getlinkvalue(ph, i, EN.ENERGY, ctypes.byref(e_val))
                    api.lib.EN_getlinkvalue(ph, i, EN.PUMP_EFFIC, ctypes.byref(ef_val))
                    link_series[lid]["energy"].append(round(e_val.value, 2))
                    # EN_PUMP_EFFIC is a decimal ratio (0.75 for 75%), convert to percent
                    link_series[lid]["efficiency"].append(round(ef_val.value * 100, 1))
            api.lib.EN_nextH(ph, ctypes.byref(tstep))
            if tstep.value <= 0:
                break

        # Close hydraulics before running quality
        api.lib.EN_closeH(ph)
        hydraulics_opened = False

        # 4. Check and run Water Quality Simulation (Chlorine/Fluoride decay, Water Age, Chemical Tracer)
        qual_type = ctypes.c_int(0)
        trace_node = ctypes.c_int(0)
        chem_name = ctypes.create_string_buffer(256)
        chem_units = ctypes.create_string_buffer(256)
        quality_info = None
        quality_failed = False
        quality_error_message = None
        try:
            api.lib.EN_getqualinfo(ph, ctypes.byref(qual_type), chem_name, chem_units, ctypes.byref(trace_node))
            q_val = qual_type.value
            if q_val > 0:
                quality_info = {
                    "qualType": q_val,
                    "typeName": QUAL_TYPES.get(q_val, "CHEMICAL"),
                    "chemName": chem_name.value.decode("utf-8", errors="ignore").strip() or ("Age" if q_val == 2 else "Quality"),
                    "chemUnits": chem_units.value.decode("utf-8", errors="ignore").strip() or ("hours" if q_val == 2 else "mg/L"),
                    "traceNode": node_id_map.get(trace_node.value, "") if trace_node.value > 0 else "",
                    "status": "success",
                }
        except Exception as e:
            quality_failed = True
            quality_error_message = str(e)

        if quality_info and quality_info["qualType"] > 0:
            quality_opened = False
            try:
                api.lib.EN_openQ(ph)
                quality_opened = True
                api.lib.EN_initQ(ph, 0)
                tQ = ctypes.c_long()
                tstepQ = ctypes.c_long(1)

                while True:
                    api.lib.EN_runQ(ph, ctypes.byref(tQ))
                    for i in range(1, node_count.value + 1):
                        nid = node_id_map[i]
                        qv = ctypes.c_double()
                        api.lib.EN_getnodevalue(ph, i, EN.QUALITY, ctypes.byref(qv))
                        node_series[nid]["quality"].append(round(qv.value, 4))

                    for i in range(1, link_count.value + 1):
                        lid = link_id_map[i]
                        lqv = ctypes.c_double()
                        api.lib.EN_getlinkvalue(ph, i, EN.LINKQUAL, ctypes.byref(lqv))
                        link_series[lid]["quality"].append(round(lqv.value, 4))

                    api.lib.EN_nextQ(ph, ctypes.byref(tstepQ))
                    if tstepQ.value <= 0:
                        break
            except Exception as e:
                quality_failed = True
                quality_error_message = str(e)
                quality_info["status"] = "failed"
                quality_info["errorMessage"] = quality_error_message
            finally:
                if quality_opened:
                    try:
                        api.lib.EN_closeQ(ph)
                    except Exception:
                        pass

        # Assign initial values and full timeSeries
        num_steps = len(time_steps)
        has_valid_quality = (
            quality_info is not None
            and quality_info.get("status") == "success"
            and not quality_failed
        )

        for nid in nodes:
            if not has_valid_quality or len(node_series[nid]["quality"]) != num_steps:
                node_series[nid]["quality"] = []
            nodes[nid]["demand"] = node_series[nid]["demand"][0] if node_series[nid]["demand"] else 0.0
            nodes[nid]["head"] = node_series[nid]["head"][0] if node_series[nid]["head"] else 0.0
            nodes[nid]["pressure"] = node_series[nid]["pressure"][0] if node_series[nid]["pressure"] else 0.0
            nodes[nid]["quality"] = node_series[nid]["quality"][0] if has_valid_quality and node_series[nid]["quality"] else None
            nodes[nid]["timeSeries"] = node_series[nid]

        for lid in links:
            if not has_valid_quality or len(link_series[lid]["quality"]) != num_steps:
                link_series[lid]["quality"] = []
            links[lid]["flow"] = link_series[lid]["flow"][0] if link_series[lid]["flow"] else 0.0
            links[lid]["velocity"] = link_series[lid]["velocity"][0] if link_series[lid]["velocity"] else 0.0
            links[lid]["headloss"] = link_series[lid]["headloss"][0] if link_series[lid]["headloss"] else 0.0
            links[lid]["quality"] = link_series[lid]["quality"][0] if has_valid_quality and link_series[lid]["quality"] else None
            links[lid]["timeSeries"] = link_series[lid]

        # 5. Compute Pump Energy Metrics (kWh, peak kW, average efficiency, utilization %)
        pump_lids = [lid for lid, link in links.items() if link["type"] == EN.PUMP]
        pump_energy_map = {}
        total_duration_sec = (time_steps[-1] - time_steps[0]) if len(time_steps) > 1 else 0

        for lid in pump_lids:
            link = links[lid]
            energies = link_series[lid].get("energy", [])
            effics = link_series[lid].get("efficiency", [])
            active_effics = [ef for ef in effics if ef > 0]

            # dt-weighted kWh integration (excluding terminal reading which has no following duration)
            pump_kwh = 0.0
            pump_active_sec = 0.0
            if len(time_steps) > 1:
                for i in range(len(time_steps) - 1):
                    dt_hours = max(0.0, float(time_steps[i + 1] - time_steps[i])) / 3600.0
                    kw = energies[i] if i < len(energies) else 0.0
                    pump_kwh += kw * dt_hours
                    if kw > 0:
                        pump_active_sec += max(0.0, float(time_steps[i + 1] - time_steps[i]))
            elif energies:
                pump_kwh = float(energies[0])
                if energies[0] > 0:
                    pump_active_sec = 3600.0

            total_kwh = round(pump_kwh, 2)
            peak_kw = round(max(energies, default=0.0), 2)
            # Time-weighted (hydraulic steps are not uniform: control/tank events
            # insert short steps), i.e. kWh over hours actually running.
            avg_kw = round(pump_kwh / (pump_active_sec / 3600.0), 2) if pump_active_sec > 0 else 0.0
            avg_effic = round(sum(active_effics) / len(active_effics), 1) if active_effics else 0.0
            utilization = round((pump_active_sec / total_duration_sec) * 100.0, 1) if total_duration_sec > 0 else (100.0 if energies and energies[0] > 0 else 0.0)

            link["energy"] = energies[0] if energies else 0.0
            link["efficiency"] = effics[0] if effics else 0.0
            link["totalKwh"] = total_kwh
            link["peakKw"] = peak_kw
            link["avgKw"] = avg_kw
            link["utilization"] = utilization

            pump_energy_map[lid] = {
                "id": lid,
                "fromNode": link["fromNode"],
                "toNode": link["toNode"],
                "totalKwh": total_kwh,
                "peakKw": peak_kw,
                "avgKw": avg_kw,
                "avgEfficiency": avg_effic,
                "utilization": utilization,
            }

        # System peak kW is the maximum over all time steps of the simultaneous sum of pump kWs
        system_peak_kw = 0.0
        if pump_lids and time_steps:
            for i in range(len(time_steps)):
                step_sum = sum(link_series[lid]["energy"][i] for lid in pump_lids if i < len(link_series[lid]["energy"]))
                if step_sum > system_peak_kw:
                    system_peak_kw = step_sum

        energy_summary = {
            "pumps": pump_energy_map,
            "totalKwh": round(sum(p["totalKwh"] for p in pump_energy_map.values()), 2),
            "peakKw": round(system_peak_kw, 2),
        } if pump_energy_map else None

        if ph_opened:
            try:
                api.lib.EN_close(ph)
                ph_opened = False
            except Exception:
                pass
        report_text = ""
        if os.path.exists(rpt_path):
            try:
                with open(rpt_path, "r", encoding="utf-8", errors="ignore") as rf:
                    report_text = rf.read()
            except Exception:
                pass
        
        return {
            "engine": {
                "runtime": "native",
                "version": "2.3.5"
            },
            "timeSeconds": time_steps[0] if time_steps else 0,
            "timeSteps": time_steps,
            "units": {
                "flow": f_unit_str,
                "pressure": p_unit_str
            },
            "nodes": nodes,
            "links": links,
            "reportText": report_text,
            "qualityInfo": quality_info,
            "energySummary": energy_summary,
            "engineWarnings": engine_warnings,
            "reportStart": report_start.value,
            "reportStep": report_step.value,
        }
    except EpanetException as e:
        # The report file usually names the root cause (e.g. "Error 227: invalid
        # head curve") that the API code alone (e.g. 110) hides. Close first so
        # the report is flushed and readable.
        if hydraulics_opened:
            try:
                api.lib.EN_closeH(ph)
            except Exception:
                pass
            hydraulics_opened = False
        if ph_opened:
            try:
                api.lib.EN_close(ph)
            except Exception:
                pass
            ph_opened = False
        detail = str(e)
        causes = read_report_errors(rpt_path)
        if causes:
            detail += " | " + " / ".join(causes)
        raise HTTPException(status_code=500, detail=detail)
    finally:
        if hydraulics_opened:
            try:
                api.lib.EN_closeH(ph)
            except Exception:
                pass
        if ph_opened:
            try:
                api.lib.EN_close(ph)
            except Exception:
                pass
        if ph_created:
            try:
                api.lib.EN_deleteproject(ph)
            except Exception:
                pass