import 'server-only'

import { solveAcPowerFlow, type AcPowerFlowResult } from '../engine/ac-power-flow'
import { runContingencyAnalysis, type ContingencyReport } from '../engine/contingency'
import { analyseCalibration, validateGridModel, type CalibrationReport, type ValidationReport } from '../engine/grid-validator'
import {
  buildNetwork,
  splitIntoIslands,
  type BuiltNetwork,
  type ModelAsset,
  type ModelEdge,
  type NetworkIsland,
} from '../engine/network-model'
import { confidenceFor } from '../engine/forecast'
import { solvePowerFlow, type PowerFlowResult } from '../engine/power-flow'
import { round } from '../engine/rng'
import type { Perturbation } from '../engine/telemetry'
import { loadAssetProfiles } from './grid-service'
import { loadGridGraph } from './graph-service'
import { getTwinState, listTwins } from './twin-service'

/**
 * Physics-grade grid analysis (5.0 §1, §2, §3, §13–§16, §20, §21).
 *
 * The server side of the AC solver. Everything here runs on the server and nothing is
 * exported to the browser: §34 is explicit that no physics computation runs in a tab, and
 * a Newton–Raphson factorisation of a few hundred buses would be a poor thing to hand a
 * phone regardless.
 *
 * The two flow methods sit side by side deliberately. `fast` is the demand allocation the
 * 4.2 live view runs on — milliseconds, good enough to pan a map with. `ac` is the solved
 * load flow — voltage magnitudes and angles, reactive flow, real losses, and a convergence
 * report. A caller picks; nothing here decides for it, and every result says which it is.
 */

export type FlowMode = 'fast' | 'ac'

const MIN = 60_000

// ─────────────────────────────────────────────────────────────────────────────
// Shared model assembly
// ─────────────────────────────────────────────────────────────────────────────

export interface GridModelInputs {
  twinCode: string
  ts: number
  assets: ModelAsset[]
  edges: ModelEdge[]
  /** Codes in the twin, so callers can scope without re-querying. */
  memberCodes: string[]
}

/**
 * Assemble the electrical model for a twin at an instant.
 *
 * The topology is the whole graph among the twin's members — not the subset a scene
 * happens to draw. A backbone line carries the demand of everything behind it whether or
 * not that is on screen, and solving over the drawn subset was the 4.2 defect this
 * inherits the fix for.
 */
export async function buildGridModel(options: {
  twinCode?: string
  ts?: number
  perturbation?: Perturbation
}): Promise<GridModelInputs | null> {
  const ts = options.ts ?? Date.now()

  const twins = await listTwins()
  const twin =
    twins.find((entry) => entry.code === options.twinCode) ??
    twins.find((entry) => entry.level === 'national') ??
    twins[0]
  if (!twin) return null

  const [state, graph, profiles] = await Promise.all([
    getTwinState(twin.code, { ts, perturbation: options.perturbation }),
    loadGridGraph(),
    loadAssetProfiles(),
  ])
  if (!state) return null

  const memberCodes = state.assets.map((asset) => asset.code)
  const members = new Set(memberCodes)

  // The per-unit base must be the *nameplate* voltage, not the measured one.
  //
  // `TwinAssetState.voltageKv` is a telemetry reading — it already carries the droop the
  // solver exists to compute. Using it as its own base makes every bus exactly 1.0 pu by
  // construction, and the load flow has nothing left to find. The nameplate comes from the
  // asset profile, which is where a rating belongs.
  const nominalKv = new Map(profiles.map((profile) => [profile.code, profile.voltageKv]))

  return {
    twinCode: twin.code,
    ts,
    memberCodes,
    assets: state.assets.map<ModelAsset>((asset) => ({
      code: asset.code,
      assetTypeKey: asset.assetTypeKey,
      voltageKv: nominalKv.get(asset.code) ?? asset.voltageKv,
      capacityMw: asset.capacityMw,
      loadPct: asset.loadPct,
      loadMw: asset.loadMw,
      powerFactor: asset.powerFactor,
      isOutage: asset.isOutage,
    })),
    edges: graph.edges
      .filter((edge) => members.has(edge.from) && members.has(edge.to))
      .map<ModelEdge>((edge) => ({
        from: edge.from,
        to: edge.to,
        kind: edge.kind,
        capacityMw: edge.capacityMw,
        lengthKm: edge.lengthKm,
      })),
  }
}

// ─────────────────────────────────────────────────────────────────────────────
// Flow (§1)
// ─────────────────────────────────────────────────────────────────────────────

export interface FlowSolution {
  mode: FlowMode
  twinCode: string
  ts: number
  /** Milliseconds the solve took. Stamped by the caller, which has a clock. */
  durationMs: number
  /** Present on `ac`. Null on `fast`, which is not a solved load flow. */
  ac: AcPowerFlowResult | null
  /** Present on `fast`. Null on `ac`. */
  fast: PowerFlowResult | null
  /**
   * One entry per electrically separate area. A fragmented model is solved area by area,
   * each against its own reference, because a load flow needs exactly one slack per
   * synchronous system and several areas have several.
   */
  islandsSolved: Array<{
    slackCode: string | null
    buses: number
    energised: boolean
    converged: boolean
    solverStatus: string
    iterations: number
    generationMw: number
    demandMw: number
    lossMw: number
    /** Below 1 when the island's sampled generation was scaled to its own demand. */
    dispatchScale: number
  }>
  /** Buses in areas with no generation at all — unsupplied, and reported as such. */
  unsuppliedBuses: number
  /** Where the impedance came from — 'typical' in this build, never claimed measured. */
  impedanceSource: BuiltNetwork['impedanceSource'] | null
  islanded: string[]
  slackReason: string | null
  excluded: BuiltNetwork['excluded']
  /**
   * True only when the chosen method produced a usable answer. On `ac` this is the
   * solver's own convergence flag; a non-converged case must not be shown as a result.
   */
  usable: boolean
  /**
   * Minutes ahead of the request instant this solve was run at (§17–§18).
   *
   * Zero solves the grid as it stands. Anything else solves the conditions the twin's own
   * forward walk projects, which makes the result a *prediction of a power flow* rather
   * than a power flow — and the two must never wear the same label.
   */
  offsetMin: number
  /** `simulated` at the present instant, `predicted` at any horizon. Never `observed`. */
  stateKind: 'simulated' | 'predicted'
  /**
   * Confidence in the *conditions* the solve ran on, 0–100, decaying with the horizon.
   *
   * Not a confidence in the solver: Newton–Raphson either converged or it did not, and
   * that is reported separately. This is how much to trust the load it was handed.
   */
  forecastConfidence: number
  note: string
}

export async function solveGridFlow(options: {
  twinCode?: string
  ts?: number
  mode?: FlowMode
  perturbation?: Perturbation
  outages?: string[]
  taps?: Record<string, number>
  slackCode?: string
  tolerance?: number
  /**
   * Solve the conditions this many minutes ahead (§17–§18).
   *
   * The load comes from the twin's own forward walk rather than a second forecast model,
   * so the flow at +60 min and the forecast panel at +60 min are the same projection seen
   * two ways. A separate model here would have let the two disagree.
   */
  offsetMin?: number
}): Promise<FlowSolution | null> {
  const started = Date.now()
  const offsetMin = Math.max(0, Math.round(options.offsetMin ?? 0))
  const requestedTs = options.ts ?? Date.now()
  const model = await buildGridModel({ ...options, ts: requestedTs + offsetMin * MIN })
  if (!model) return null

  const mode: FlowMode = options.mode ?? 'fast'
  const stateKind: 'simulated' | 'predicted' = offsetMin > 0 ? 'predicted' : 'simulated'
  // The same decay the forecast panels use, so one horizon has one confidence.
  const forecastConfidence = round(confidenceFor(offsetMin) * 100, 1)

  if (mode === 'fast') {
    const fast = solvePowerFlow({
      nodes: model.assets.map((asset) => ({
        code: asset.code,
        assetTypeKey: asset.assetTypeKey,
        capacityMw: asset.capacityMw,
        loadPct: asset.loadPct,
        loadMw: asset.loadMw,
        voltageKv: asset.voltageKv,
        powerFactor: asset.powerFactor,
        isOutage: asset.isOutage || (options.outages ?? []).includes(asset.code),
      })),
      edges: model.edges.map((edge) => ({
        from: edge.from,
        to: edge.to,
        kind: edge.kind,
        capacityMw: edge.capacityMw,
      })),
    })

    return {
      mode,
      twinCode: model.twinCode,
      ts: model.ts,
      durationMs: Date.now() - started,
      ac: null,
      fast,
      impedanceSource: null,
      islanded: [],
      slackReason: null,
      excluded: [],
      islandsSolved: [],
      unsuppliedBuses: 0,
      usable: true,
      offsetMin,
      stateKind,
      forecastConfidence,
      note:
        offsetMin > 0
          ? `Demand allocation over the topology, on conditions projected ${offsetMin} min ahead. A prediction of a flow, not a flow — and not a solved load flow either.`
          : 'Demand allocation over the topology. Fast, and not a solved load flow — there are no voltage angles in this answer.',
    }
  }

  const built = buildNetwork(model.assets, model.edges, {
    outages: options.outages,
    taps: options.taps,
    slackCode: options.slackCode,
  })

  // Solve each electrically separate area against its own reference and merge the
  // results. Forcing one slack onto twenty-three separate systems gives a matrix with no
  // solution, which is not the same thing as a grid with a problem.
  const islands = splitIntoIslands(built.model)
  const solved = islands.map((island) => ({
    island,
    result: island.energised
      ? solveAcPowerFlow(island.model, {
          tolerance: options.tolerance,
          maxIterations: 30,
          enforceQLimits: true,
        })
      : null,
  }))

  const ac = mergeIslandResults(solved, built.model.baseMva)
  const unsuppliedBuses = islands
    .filter((island) => !island.energised)
    .reduce((sum, island) => sum + island.codes.length, 0)

  return {
    mode,
    twinCode: model.twinCode,
    ts: model.ts,
    durationMs: Date.now() - started,
    ac,
    fast: null,
    impedanceSource: built.impedanceSource,
    islanded: built.islanded,
    slackReason: built.slackReason,
    excluded: built.excluded,
    islandsSolved: solved.map((entry) => ({
      slackCode: entry.island.slackCode,
      buses: entry.island.codes.length,
      energised: entry.island.energised,
      converged: entry.result?.converged ?? false,
      solverStatus: entry.result?.solverStatus ?? 'no_slack_bus',
      iterations: entry.result?.iterations ?? 0,
      generationMw: entry.island.generationMw,
      demandMw: entry.island.demandMw,
      lossMw: entry.result?.converged ? entry.result.totalLossMw : 0,
      dispatchScale: entry.island.dispatchScale,
    })),
    unsuppliedBuses,
    usable: ac.converged,
    offsetMin,
    stateKind,
    forecastConfidence,
    note: ac.converged
      ? `${ac.solveNote} Solved across ${solved.filter((e) => e.island.energised).length} energised areas. Impedance is typical for the voltage class, not measured on these assets.${
          offsetMin > 0
            ? ` Conditions are projected ${offsetMin} min ahead, at ${forecastConfidence} % confidence — this is a predicted flow, not a measured one.`
            : ''
        }`
      : ac.solveNote,
  }
}

/**
 * Merge per-island solutions into one result.
 *
 * `converged` on the merged object is true only when every energised island converged:
 * a partial answer must not read as a whole one. Unenergised islands contribute their
 * buses to nothing — they have no solution to contribute.
 */
function mergeIslandResults(
  solved: Array<{ island: NetworkIsland; result: AcPowerFlowResult | null }>,
  baseMva: number,
): AcPowerFlowResult {
  const energised = solved.filter((entry) => entry.island.energised && entry.result)
  const results = energised.map((entry) => entry.result!)

  const converged = results.length > 0 && results.every((result) => result.converged)
  const worstConvergence = results.reduce(
    (worst, result) => Math.max(worst, result.convergence),
    0,
  )
  const iterations = results.reduce((most, result) => Math.max(most, result.iterations), 0)
  const failed = results.filter((result) => !result.converged)

  return {
    method: 'newton_raphson',
    solverStatus: converged ? 'converged' : (failed[0]?.solverStatus ?? 'no_slack_bus'),
    converged,
    iterations,
    convergence: worstConvergence,
    tolerance: results[0]?.tolerance ?? 1e-6,
    baseMva,
    mismatchHistory: results[0]?.mismatchHistory ?? [],
    // Only from areas that actually solved: a bus from a diverged island has no solved
    // voltage to report.
    buses: results.filter((r) => r.converged).flatMap((r) => r.buses),
    branches: results.filter((r) => r.converged).flatMap((r) => r.branches),
    slackCode: energised[0]?.island.slackCode ?? null,
    totalGenerationMw: round(
      results.filter((r) => r.converged).reduce((sum, r) => sum + r.totalGenerationMw, 0),
      3,
    ),
    totalLoadMw: round(
      results.filter((r) => r.converged).reduce((sum, r) => sum + r.totalLoadMw, 0),
      3,
    ),
    totalLossMw: round(
      results.filter((r) => r.converged).reduce((sum, r) => sum + r.totalLossMw, 0),
      4,
    ),
    totalLossMvar: round(
      results.filter((r) => r.converged).reduce((sum, r) => sum + r.totalLossMvar, 0),
      4,
    ),
    voltageViolations: results.filter((r) => r.converged).flatMap((r) => r.voltageViolations),
    overloads: results.filter((r) => r.converged).flatMap((r) => r.overloads),
    solveNote: converged
      ? `Every one of ${results.length} energised areas converged, worst mismatch ${worstConvergence.toExponential(2)} pu.`
      : `${failed.length} of ${results.length} energised areas did not converge. Their buses are omitted; the numbers here describe only the areas that solved.`,
  }
}

// ─────────────────────────────────────────────────────────────────────────────
// Validation and calibration (§2, §3)
// ─────────────────────────────────────────────────────────────────────────────

export interface GridHealthReport {
  twinCode: string
  ts: number
  validation: ValidationReport
  calibration: CalibrationReport
  /** Never true when validation found an error: a broken model is not ready to solve. */
  readyToSolve: boolean
}

export async function inspectGridModel(options: {
  twinCode?: string
  ts?: number
  perturbation?: Perturbation
}): Promise<GridHealthReport | null> {
  const model = await buildGridModel(options)
  if (!model) return null

  const validation = validateGridModel(model.assets, model.edges)
  const calibration = analyseCalibration(model.assets, model.edges)

  return {
    twinCode: model.twinCode,
    ts: model.ts,
    validation,
    calibration,
    readyToSolve: validation.solvable,
  }
}

// ─────────────────────────────────────────────────────────────────────────────
// Contingency (§13, §14)
// ─────────────────────────────────────────────────────────────────────────────

export interface ContingencyScope {
  /** The electrical area the scan covered. */
  slackCode: string | null
  buses: number
  /** Share of the twin's assets inside that area. */
  coveragePct: number
  reason: string
}

/**
 * N-1 over a *solvable* area.
 *
 * A contingency scan needs a base case that solves. On a fragmented model the whole twin
 * has no single solution — there is no one reference for twenty-three separate systems —
 * so the scan is scoped to the largest area whose base case converges, and the scope is
 * reported. Scanning the whole model and refusing every case would be technically correct
 * and operationally useless.
 */
export async function runGridContingency(options: {
  twinCode?: string
  ts?: number
  perturbation?: Perturbation
  depth?: 1 | 2
  maxCases?: number
  candidates?: string[]
  /** Restrict to the area containing this asset. Otherwise the largest solvable one. */
  scopeAsset?: string
}): Promise<
  | (ContingencyReport & { twinCode: string; ts: number; durationMs: number; scope: ContingencyScope })
  | null
> {
  const started = Date.now()
  const model = await buildGridModel(options)
  if (!model) return null

  const built = buildNetwork(model.assets, model.edges)
  const islands = splitIntoIslands(built.model)

  // Which area to scan: the one holding the named asset if there is one, else the largest
  // whose base case actually converges.
  const solvable = islands
    .filter((island) => island.energised)
    .map((island) => ({
      island,
      converged: solveAcPowerFlow(island.model, { maxIterations: 25, enforceQLimits: true }).converged,
    }))

  const named = options.scopeAsset
    ? solvable.find((entry) => entry.island.codes.includes(options.scopeAsset!))
    : undefined
  const chosen = named ?? solvable.find((entry) => entry.converged) ?? solvable[0]

  if (!chosen) {
    return {
      ...runContingencyAnalysis([], [], {}),
      twinCode: model.twinCode,
      ts: model.ts,
      durationMs: Date.now() - started,
      scope: {
        slackCode: null,
        buses: 0,
        coveragePct: 0,
        reason: 'No area of this model has generation in service, so there is no base case to scan against.',
      },
    }
  }

  const members = new Set(chosen.island.codes)
  const report = runContingencyAnalysis(
    model.assets.filter((asset) => members.has(asset.code)),
    model.edges.filter((edge) => members.has(edge.from) && members.has(edge.to)),
    {
      depth: options.depth ?? 1,
      maxCases: options.maxCases,
      candidates: options.candidates,
      maxIterations: 20,
    },
  )

  return {
    ...report,
    twinCode: model.twinCode,
    ts: model.ts,
    durationMs: Date.now() - started,
    scope: {
      slackCode: chosen.island.slackCode,
      buses: chosen.island.codes.length,
      coveragePct: round((chosen.island.codes.length / Math.max(1, model.assets.length)) * 100, 1),
      reason: named
        ? `The electrical area containing ${options.scopeAsset}.`
        : chosen.converged
          ? 'The largest area whose base case solves. A scan needs a base case to compare against.'
          : 'No area solved its base case; scanning the largest energised one so the failure is visible rather than silent.',
    },
  }
}

// ─────────────────────────────────────────────────────────────────────────────
// Voltage stability (§15)
// ─────────────────────────────────────────────────────────────────────────────

export interface VoltageStability {
  twinCode: string
  ts: number
  converged: boolean
  solverStatus: string
  /** Share of the model's buses this assessment actually covers. */
  coveragePct: number
  busesAssessed: number
  busesUnsolved: number
  /** 0–100, higher is stronger. Derived from the weakest bus's margin to its limit. */
  stabilityIndex: number
  weakBuses: Array<{
    code: string
    vPu: number
    angleDeg: number
    /** Percentage points between this bus and its nearer statutory limit. */
    marginPct: number
    /** Reactive power the bus is drawing — the quantity that collapses a voltage. */
    reactiveDemandMvar: number
    band: 'sound' | 'watch' | 'weak' | 'critical'
  }>
  criticalReactiveMvar: number
  collapseRisk: 'low' | 'moderate' | 'elevated' | 'high'
  note: string
}

const V_LOW = 0.94
const V_HIGH = 1.06

export async function assessVoltageStability(options: {
  twinCode?: string
  ts?: number
  perturbation?: Perturbation
}): Promise<VoltageStability | null> {
  const solution = await solveGridFlow({ ...options, mode: 'ac' })
  if (!solution || !solution.ac) return null

  const ac = solution.ac

  // Areas that solved are assessed; areas that did not are reported as unassessed rather
  // than as stable or unstable. Discarding a hundred and thirty-nine solved buses because
  // two areas elsewhere failed would throw away the answer to keep the question tidy.
  const totalBuses = ac.buses.length + solution.unsuppliedBuses + unsolvedBusesIn(solution)
  const coveragePct = totalBuses > 0 ? round((ac.buses.length / totalBuses) * 100, 1) : 0

  if (ac.buses.length === 0) {
    // Nothing solved anywhere. A load flow that will not converge is the classic numerical
    // signature of a system past the nose of its PV curve — but it is equally the
    // signature of a model fault, and this says so rather than choosing.
    return {
      twinCode: solution.twinCode,
      ts: solution.ts,
      converged: false,
      solverStatus: ac.solverStatus,
      coveragePct: 0,
      busesAssessed: 0,
      busesUnsolved: totalBuses,
      stabilityIndex: 0,
      weakBuses: [],
      criticalReactiveMvar: 0,
      collapseRisk: 'high',
      note: 'No area of this model solved, so there is no voltage state to assess. On a real network that is what collapse looks like numerically; it is equally what a broken model looks like. Validate the model before concluding anything about the grid.',
    }
  }

  const weakBuses = ac.buses
    .map((bus) => {
      const marginPct = round(
        Math.min(Math.abs(bus.vPu - V_LOW), Math.abs(V_HIGH - bus.vPu)) * 100,
        3,
      )
      return {
        code: bus.code,
        vPu: bus.vPu,
        angleDeg: bus.angleDeg,
        marginPct,
        reactiveDemandMvar: round(Math.max(0, -bus.qMvar), 3),
        band: (marginPct <= 0
          ? 'critical'
          : marginPct < 1
            ? 'weak'
            : marginPct < 3
              ? 'watch'
              : 'sound') as 'sound' | 'watch' | 'weak' | 'critical',
      }
    })
    .sort((a, b) => a.marginPct - b.marginPct)

  const weakest = weakBuses[0]
  // Capped by the weakest bus rather than averaged across the fleet: a network is as
  // stable as its worst bus, and a mean would let a hundred healthy ones hide it.
  const stabilityIndex = round(Math.max(0, Math.min(100, (weakest?.marginPct ?? 0) * 16)), 1)

  const criticalReactiveMvar = round(
    weakBuses.filter((bus) => bus.band !== 'sound').reduce((sum, bus) => sum + bus.reactiveDemandMvar, 0),
    2,
  )

  return {
    twinCode: solution.twinCode,
    ts: solution.ts,
    converged: ac.converged,
    solverStatus: ac.solverStatus,
    coveragePct,
    busesAssessed: ac.buses.length,
    busesUnsolved: totalBuses - ac.buses.length,
    stabilityIndex,
    weakBuses: weakBuses.slice(0, 20),
    criticalReactiveMvar,
    collapseRisk:
      stabilityIndex < 15 ? 'high' : stabilityIndex < 35 ? 'elevated' : stabilityIndex < 60 ? 'moderate' : 'low',
    note:
      coveragePct >= 99.5
        ? `Margin measured to the nearer of the ${V_LOW}–${V_HIGH} pu limits. The index follows the weakest bus, not the average.`
        : `Covers ${coveragePct} % of the model — the areas that solved. The rest is unassessed, not stable.`,
  }
}

/** Buses inside energised areas whose own solve failed. */
function unsolvedBusesIn(solution: FlowSolution): number {
  return solution.islandsSolved
    .filter((island) => island.energised && !island.converged)
    .reduce((sum, island) => sum + island.buses, 0)
}

// ─────────────────────────────────────────────────────────────────────────────
// Thermal limits (§16)
// ─────────────────────────────────────────────────────────────────────────────

export interface ThermalEntry {
  code: string
  kind: 'asset' | 'branch'
  assetTypeKey: string
  /** Present only where a temperature is actually measured or modelled. */
  tempC: number | null
  ratedTempC: number | null
  headroomC: number | null
  currentLoadingPct: number
  thermalLoadingPct: number | null
  /** Minutes to the thermal limit at the present rate of rise. Null when not rising. */
  timeToLimitMin: number | null
  /** Stated whenever a figure could not be computed, rather than a zero. */
  unavailableReason: string | null
}

export interface ThermalReport {
  twinCode: string
  ts: number
  entries: ThermalEntry[]
  worst: ThermalEntry | null
  withoutData: number
  note: string
}

/**
 * Thermal loading and time to limit.
 *
 * Computed only where the data exists. §16 says "only when the required data is
 * available", and the honest expression of that is an explicit `unavailableReason` rather
 * than a zero that reads like a measurement.
 */
export async function assessThermalLimits(options: {
  twinCode?: string
  ts?: number
  perturbation?: Perturbation
}): Promise<ThermalReport | null> {
  const ts = options.ts ?? Date.now()

  const twins = await listTwins()
  const twin =
    twins.find((entry) => entry.code === options.twinCode) ??
    twins.find((entry) => entry.level === 'national') ??
    twins[0]
  if (!twin) return null

  // Two instants so a rate of rise can be measured rather than assumed.
  const [now, earlier] = await Promise.all([
    getTwinState(twin.code, { ts, perturbation: options.perturbation }),
    getTwinState(twin.code, { ts: ts - 15 * MIN, perturbation: options.perturbation }),
  ])
  if (!now) return null

  const before = new Map((earlier?.assets ?? []).map((asset) => [asset.code, asset]))

  const entries: ThermalEntry[] = now.assets.map((asset) => {
    const thermalKinds = ['transformer', 'substation', 'power_plant', 'battery_storage']
    const hasThermalModel = thermalKinds.includes(asset.assetTypeKey)

    if (!hasThermalModel) {
      return {
        code: asset.code,
        kind: 'asset',
        assetTypeKey: asset.assetTypeKey,
        tempC: null,
        ratedTempC: null,
        headroomC: null,
        currentLoadingPct: round(Math.abs(asset.loadPct), 2),
        thermalLoadingPct: null,
        timeToLimitMin: null,
        unavailableReason: 'No thermal model for this asset type; loading is reported without a temperature.',
      }
    }

    // The rated temperature the twin's thermal stress was measured against.
    const ratedTempC = round(asset.tempC - asset.thermalStress * 30 + 15, 1)
    const headroomC = round(ratedTempC - asset.tempC, 2)

    const previous = before.get(asset.code)
    const risePerMin = previous ? (asset.tempC - previous.tempC) / 15 : 0
    const timeToLimitMin =
      risePerMin > 0.01 && headroomC > 0 ? round(headroomC / risePerMin, 1) : null

    return {
      code: asset.code,
      kind: 'asset',
      assetTypeKey: asset.assetTypeKey,
      tempC: asset.tempC,
      ratedTempC,
      headroomC,
      currentLoadingPct: round(Math.abs(asset.loadPct), 2),
      thermalLoadingPct: round(asset.thermalStress * 100, 2),
      timeToLimitMin,
      unavailableReason:
        risePerMin <= 0.01 ? 'Temperature is not rising, so there is no time-to-limit to report.' : null,
    }
  })

  const withTemperature = entries.filter((entry) => entry.tempC !== null)
  const worst = withTemperature.sort((a, b) => (a.headroomC ?? 0) - (b.headroomC ?? 0))[0] ?? null

  return {
    twinCode: twin.code,
    ts,
    entries,
    worst,
    withoutData: entries.length - withTemperature.length,
    note: 'Rate of rise is measured across a fifteen-minute window. Assets with no thermal model report loading only, and say so.',
  }
}

// ─────────────────────────────────────────────────────────────────────────────
// Energy balance (§20)
// ─────────────────────────────────────────────────────────────────────────────

export interface EnergyBalance {
  twinCode: string
  ts: number
  method: 'newton_raphson' | 'demand_allocation'
  generationMw: number
  storageDischargeMw: number
  storageChargeMw: number
  loadMw: number
  lossMw: number
  /** Generation + discharge − load − losses − charge. Reported, never forced to zero. */
  balanceErrorMw: number
  balanceErrorPct: number
  /** Capped by the weakest contributing figure, not averaged. */
  dataQualityPct: number
  converged: boolean
  note: string
}

export async function computeEnergyBalance(options: {
  twinCode?: string
  ts?: number
  perturbation?: Perturbation
  mode?: FlowMode
}): Promise<EnergyBalance | null> {
  const mode = options.mode ?? 'ac'
  const [solution, model] = await Promise.all([
    solveGridFlow({ ...options, mode }),
    buildGridModel(options),
  ])
  if (!solution || !model) return null

  const storage = model.assets.filter((asset) => asset.assetTypeKey === 'battery_storage' && !asset.isOutage)
  const storageDischargeMw = round(
    storage.filter((asset) => asset.loadPct >= 0).reduce((sum, asset) => sum + Math.abs(asset.loadMw), 0),
    3,
  )
  const storageChargeMw = round(
    storage.filter((asset) => asset.loadPct < 0).reduce((sum, asset) => sum + Math.abs(asset.loadMw), 0),
    3,
  )

  const generationMw = solution.ac
    ? solution.ac.totalGenerationMw
    : (solution.fast?.totalGenerationMw ?? 0)
  const loadMw = solution.ac ? solution.ac.totalLoadMw : (solution.fast?.totalDemandMw ?? 0)
  const lossMw = solution.ac ? solution.ac.totalLossMw : (solution.fast?.totalLossMw ?? 0)

  const balanceErrorMw = round(generationMw + storageDischargeMw - loadMw - lossMw - storageChargeMw, 3)
  const denominator = Math.max(1, generationMw + storageDischargeMw)

  // Islanded members and excluded assets are both reasons the balance may not close, and
  // both reduce how much the figure can be trusted.
  const coverage = model.assets.length > 0
    ? ((model.assets.length - solution.islanded.length - solution.excluded.length) / model.assets.length) * 100
    : 0
  const convergenceQuality = solution.usable ? 100 : 25

  return {
    twinCode: solution.twinCode,
    ts: solution.ts,
    method: solution.ac ? 'newton_raphson' : 'demand_allocation',
    generationMw,
    storageDischargeMw,
    storageChargeMw,
    loadMw,
    lossMw,
    balanceErrorMw,
    balanceErrorPct: round((balanceErrorMw / denominator) * 100, 3),
    dataQualityPct: round(Math.min(coverage, convergenceQuality), 1),
    converged: solution.usable,
    note: solution.usable
      ? 'A balance that does not close is reported rather than absorbed into a slack term — the residual is the finding.'
      : 'The solve did not converge, so this balance is arithmetic over an unsolved state and should not be read as a measurement.',
  }
}

// ─────────────────────────────────────────────────────────────────────────────
// State estimation (§21)
// ─────────────────────────────────────────────────────────────────────────────

export type StateOrigin = 'measured' | 'estimated'

export interface EstimatedQuantity {
  assetCode: string
  quantity: 'voltage' | 'current' | 'active_power' | 'reactive_power' | 'temperature' | 'frequency'
  value: number
  unit: string
  origin: StateOrigin
  /** The instrument behind a measured value; null when the value was estimated. */
  sensorCode: string | null
  /** 0–100. An estimate is never given the confidence of a measurement. */
  confidence: number
}

export interface StateEstimate {
  twinCode: string
  ts: number
  quantities: EstimatedQuantity[]
  measuredCount: number
  estimatedCount: number
  /** Share of the state that rests on an instrument rather than on the model. */
  observabilityPct: number
  note: string
}

/**
 * Separate what is measured from what is inferred.
 *
 * §21 asks for the distinction to be explicit and never blurred, and that is exactly what
 * this returns: every quantity carries an `origin`, and an estimated one carries no
 * sensor code because there is no instrument behind it. A surface that renders these
 * identically is the failure mode the requirement exists to prevent.
 */
export async function estimateGridState(options: {
  twinCode?: string
  ts?: number
  perturbation?: Perturbation
}): Promise<StateEstimate | null> {
  const { prisma } = await import('../db')
  const ts = options.ts ?? Date.now()

  const twins = await listTwins()
  const twin =
    twins.find((entry) => entry.code === options.twinCode) ??
    twins.find((entry) => entry.level === 'national') ??
    twins[0]
  if (!twin) return null

  const state = await getTwinState(twin.code, { ts, perturbation: options.perturbation })
  if (!state) return null

  const codes = state.assets.map((asset) => asset.code)
  const sensors = await prisma.sensor.findMany({
    where: { asset: { code: { in: codes } } },
    select: {
      code: true,
      kind: true,
      status: true,
      accuracy: true,
      asset: { select: { code: true } },
      trust: { select: { score: true } },
    },
  })

  // Which quantities each asset actually has an instrument for.
  const instruments = new Map<string, Map<string, { code: string; trust: number }>>()
  for (const sensor of sensors) {
    const assetCode = sensor.asset?.code
    if (!assetCode) continue
    const forAsset = instruments.get(assetCode) ?? new Map()
    // A sensor that is not reporting is not an instrument for this purpose.
    if (sensor.status === 'healthy' || sensor.status === 'degraded') {
      forAsset.set(sensor.kind, { code: sensor.code, trust: sensor.trust?.score ?? 70 })
    }
    instruments.set(assetCode, forAsset)
  }

  type TwinAssetRow = (typeof state.assets)[number]
  const map: Array<{
    quantity: EstimatedQuantity['quantity']
    sensorKind: string
    unit: string
    pick: (asset: TwinAssetRow) => number
  }> = [
    { quantity: 'voltage', sensorKind: 'voltage', unit: 'kV', pick: (a) => a.voltageKv },
    { quantity: 'current', sensorKind: 'current', unit: 'A', pick: (a) => a.currentA },
    { quantity: 'active_power', sensorKind: 'power', unit: 'MW', pick: (a) => a.loadMw },
    { quantity: 'temperature', sensorKind: 'temperature', unit: '°C', pick: (a) => a.tempC },
  ]

  const quantities: EstimatedQuantity[] = []
  for (const asset of state.assets) {
    const forAsset = instruments.get(asset.code)
    for (const entry of map) {
      const sensor = forAsset?.get(entry.sensorKind)
      quantities.push({
        assetCode: asset.code,
        quantity: entry.quantity,
        value: round(entry.pick(asset), 4),
        unit: entry.unit,
        origin: sensor ? 'measured' : 'estimated',
        sensorCode: sensor?.code ?? null,
        // An estimate is capped well below a measurement: it is the model talking about
        // an asset nothing is watching.
        confidence: sensor ? round(Math.min(97, sensor.trust), 1) : 55,
      })
    }
  }

  const measuredCount = quantities.filter((entry) => entry.origin === 'measured').length

  return {
    twinCode: twin.code,
    ts,
    quantities,
    measuredCount,
    estimatedCount: quantities.length - measuredCount,
    observabilityPct: quantities.length > 0 ? round((measuredCount / quantities.length) * 100, 1) : 0,
    note: 'An estimated quantity has no sensor behind it and is never presented as measured. Confidence on an estimate is capped below any measurement.',
  }
}
