powerio_matrix/io/

gridfm.rs

//! GridFM output: a parsed case as the gridfm-datakit Parquet schema.
//!
//! [`gridfm-datakit`](https://github.com/gridfm) writes per-scenario Parquet
//! tables that [`gridfm-graphkit`](https://github.com/gridfm)'s
//! `HeteroGridDatasetDisk` trains a GNN on. This module emits the same four
//! tables — `bus_data`, `gen_data`, `branch_data`, `y_bus_data` — from one
//! parsed [`BalancedNetwork`], so graphkit can train on powerio output directly and the
//! scenario-batch path (issue #14) has its on-disk format.
//!
//! The reverse — `powerio::gridfm::read_gridfm_dataset` /
//! `read_gridfm_scenarios` / `read_gridfm_base_case`, with the pure
//! `read_gridfm_network` over in-memory batches — lives in the facade and
//! rebuilds a [`BalancedNetwork`] from such a dataset (lossy but
//! power flow complete), the ML→classical return leg
//! (issue #60). One reader plus the existing writers means gridfm → any classical
//! format. `y_bus_data` is ignored on read; branches carry raw `r/x/b`.
//!
//! # Snapshots and scenarios
//!
//! powerio has no power flow solver. One parsed case is one snapshot
//! (`scenario = 0`): voltages and generator dispatch are the case's stored
//! values. Branch flows `pf/qf/pt/qt` use an explicit [`crate::network::BranchSolution`]
//! when present; otherwise they are evaluated from those voltages and the
//! branch admittances (`branch_flows`). For an unsolved or flat-start case the
//! evaluated values are flows at the stored voltages, not a re-solved dispatch.
//!
//! A scenario batch ([`emit_gridfm_batch`] / [`to_gridfm_record_batches`])
//! row-stacks many snapshots into the four tables, keyed by the `scenario`
//! column. The snapshots share a base element set — the same bus/branch/gen
//! counts, system base, and bus/branch/generator row identities, so every dense
//! index means the same element across scenarios — enforced by the shape check
//! ([`Error::ScenarioShapeMismatch`]). Scenario ids are unique within a batch.
//! Within that, load, dispatch, voltages, branch status, bus type, and costs may
//! all differ per snapshot. This matches datakit, whose topology scenarios (N-K,
//! random component drop) toggle `BR_STATUS`/`GEN_STATUS` on a fixed element set,
//! and graphkit's `HeteroGridDatasetDisk`, which groups by `scenario` and
//! rebuilds the graph independently for each one. powerio doesn't generate the
//! perturbations; a caller (e.g. a scenario generator) supplies the snapshots.
//!
//! # Units
//!
//! - `Pd, Qd, Pg, Qg, p_mw, q_mvar` are MW/MVAr, passed through from the case
//!   (loads and generator setpoints are already MW/MVAr). The branch flows
//!   `pf, qf, pt, qt` are MW/MVAr too, computed in per-unit and scaled by
//!   `base_mva`.
//! - `Vm` per-unit, `Va` degrees; `r, x, b` and the `Y**` admittances per-unit.
//! - `GS, BS` are the MATPOWER shunt values (MW/MVAr at V = 1) divided by
//!   `base_mva`, matching datakit's normalization.
//! - Costs are the raw MATPOWER coefficients: `cp2 = c2`, `cp1 = c1`,
//!   `cp0 = c0`. A cost row gridfm can't represent (piecewise, missing,
//!   malformed, or cubic and higher) emits zeros — graphkit ignores the cost
//!   columns — and is counted in the manifest. The `_eur` suffixes are
//!   datakit's column names, not a unit powerio converts to.
//! - `bus`, `from_bus`, `to_bus` are dense `[0, n)` indices; `idx` is the
//!   0-based generator/branch row. An out of service branch keeps its physical
//!   `Y**` admittances but carries zero flows (its `br_status` is 0).

// Bus/branch indices and Y_bus nnz counts cast to `i64` for the Arrow columns;
// they are bounded far below `i64::MAX`, so the wrap clippy warns about can't
// happen.
#![allow(clippy::cast_possible_wrap)]

use std::path::{Path, PathBuf};
use std::sync::Arc;

use arrow::array::{Array, ArrayRef, Float64Array, Int64Array};
use arrow::datatypes::{Field, Schema};
use arrow::record_batch::RecordBatch;
use num_complex::Complex64;

use parquet::arrow::ArrowWriter;
use parquet::basic::Compression;
use parquet::file::properties::WriterProperties;
use serde::Serialize;

use crate::indexed::IndexedNetwork;
use crate::matrix::{
    BuildOptions, YbusFlags, branch_admittance, branch_flows, calc_admittance_matrix,
};
use crate::network::{BusId, BusType};
use crate::{BalancedNetwork, ElementCounts, Error, GenCost, Result, ScenarioMismatch};
use powerio_tx::{DEFAULT_BASE_FREQUENCY, GenCostPatch, MissingGenCostPolicy};

/// Options for the gridfm export — the batch-wide knobs. The scenario id is a
/// per-snapshot property (set via [`GridfmSnapshot::new`] / [`number_snapshots`],
/// or the explicit argument to the single-case [`emit_gridfm_dataset`] /
/// [`to_gridfm_record_batches_single`]), not an option here.
#[derive(Debug, Clone)]
pub struct GridfmOptions {
    /// Also write `y_bus_data.parquet`. graphkit reconstructs admittances from
    /// the branch table and ignores it, but datakit emits it, so the default is
    /// `true` for parity.
    pub include_y_bus: bool,
    /// Apply transformer tap ratios to the admittances. Default `true` (the
    /// physical admittances datakit stores).
    pub include_taps: bool,
    /// Apply phase shifts to the admittances. Default `true`.
    pub include_shifts: bool,
    /// How missing generator cost rows are handled before export.
    pub missing_gen_cost: MissingGenCostPolicy,
    /// Explicit generator cost patches applied to each snapshot before the missing
    /// cost policy.
    pub gen_cost_patches: Vec<GenCostPatch>,
}

impl Default for GridfmOptions {
    fn default() -> Self {
        Self {
            include_y_bus: true,
            include_taps: true,
            include_shifts: true,
            missing_gen_cost: MissingGenCostPolicy::Preserve,
            gen_cost_patches: Vec::new(),
        }
    }
}

impl GridfmOptions {
    /// The Y_bus build flags these options select (only taps and shifts matter
    /// for the gridfm admittances; the Bp/Bpp scheme doesn't apply). A zero
    /// impedance branch is always zeroed and counted, never a build error:
    /// the dataset export has its own reporting for it
    /// (`GridfmOutputs::dropped_zero_impedance`), matching every branch's
    /// row appearing in the columnar table regardless.
    fn build_options(&self) -> BuildOptions {
        BuildOptions {
            include_taps: self.include_taps,
            include_shifts: self.include_shifts,
            skip_zero_impedance: true,
            ..Default::default()
        }
    }

    fn cost_options_default(&self) -> bool {
        self.missing_gen_cost.is_preserve() && self.gen_cost_patches.is_empty()
    }
}

/// One snapshot in a gridfm scenario batch: a parsed [`BalancedNetwork`] and the scenario
/// id stamped into its rows.
///
/// powerio has no solver, so each snapshot is an operating point a caller (e.g. a
/// scenario generator) has already produced. Snapshots in one batch share a base
/// element set — the same bus/branch/gen counts, system base, and element row
/// identities — so every dense index means the same element across snapshots
/// and the tables stay schema-consistent. Scenario ids are unique within a
/// batch. The builders enforce these conditions and otherwise return
/// [`Error::ScenarioShapeMismatch`]. Within that, load, dispatch, voltages,
/// branch status, bus type, and costs may all vary per snapshot — this mirrors
/// gridfm-datakit, whose topology scenarios (N-K, random component drop) toggle
/// `BR_STATUS`/`GEN_STATUS` on a fixed element set, and gridfm-graphkit, which
/// rebuilds the graph independently for every scenario.
#[derive(Debug, Clone, Copy)]
pub struct GridfmSnapshot<'a> {
    /// The parsed case for this scenario.
    net: &'a BalancedNetwork,
    /// The scenario id stamped into the `scenario`/`load_scenario_idx` columns.
    scenario: i64,
}

impl<'a> GridfmSnapshot<'a> {
    /// A snapshot pairing a network with the scenario id stamped into its rows.
    /// For the common "k-th input is `base + k`" numbering, prefer
    /// [`number_snapshots`], which assigns ids with checked arithmetic.
    #[must_use]
    pub fn new(net: &'a BalancedNetwork, scenario: i64) -> Self {
        Self { net, scenario }
    }
}

/// The gridfm-datakit tables as Arrow record batches. The Parquet writer builds
/// from these; a deferred gridfm-schema Arrow C Data Interface export (issue #38)
/// would reuse them. (The raw network Arrow export that ships in powerio-capi is
/// a different, lighter schema.)
///
/// For a scenario batch the tables are row stacked: each table holds the rows of
/// every snapshot back-to-back, keyed by the `scenario` column (0-based dense bus
/// indices and generator/branch `idx` reset per scenario).
#[derive(Debug, Clone)]
#[non_exhaustive]
pub struct GridfmTables {
    pub bus: RecordBatch,
    pub generator: RecordBatch,
    pub branch: RecordBatch,
    /// `None` when [`GridfmOptions::include_y_bus`] is off — the table isn't
    /// built (graphkit reconstructs admittances from the branch table anyway).
    pub y_bus: Option<RecordBatch>,
}

/// What [`emit_gridfm_dataset`] wrote, plus the counts of columns it had to zero
/// (see the manifest) so a caller can surface them.
#[derive(Debug, Clone)]
pub struct GridfmOutputs {
    pub dir: PathBuf,
    pub files: Vec<PathBuf>,
    /// Branches whose impedance is too small to divide by, so their admittance
    /// and derived flow columns were zeroed.
    pub dropped_zero_impedance: usize,
    /// Generators whose cost row gridfm couldn't represent, whose `cp*` columns
    /// were zeroed.
    pub degenerate_cost_gens: usize,
    /// Generators with no cost after applying the export cost policy.
    pub missing_cost_gens: usize,
    /// Generators with a present cost row gridfm cannot represent.
    pub unsupported_cost_gens: usize,
    /// Generators whose missing cost was filled by the export policy.
    pub synthesized_gen_costs: usize,
    /// Generator costs replaced by explicit user patches.
    pub patched_gen_costs: usize,
    /// Projection losses and values GridFM required the writer to derive.
    pub diagnostics: Vec<powerio_core::Diagnostic>,
}

/// A complete GridFM directory ready for one [`powerio_core::Destination`].
///
/// Artifact names are relative to the requested output directory and begin
/// with `raw/`. The counts state every field the GridFM profile could not
/// carry without replacing it with a default.
#[derive(Debug)]
#[non_exhaustive]
pub struct GridfmDataset {
    pub artifacts: Vec<powerio_core::MemoryArtifact>,
    pub diagnostics: Vec<powerio_core::Diagnostic>,
    pub dropped_zero_impedance: usize,
    pub degenerate_cost_gens: usize,
    pub missing_cost_gens: usize,
    pub unsupported_cost_gens: usize,
    pub synthesized_gen_costs: usize,
    pub patched_gen_costs: usize,
}

#[derive(Serialize)]
struct GridfmMeta {
    case_name: String,
    base_mva: f64,
    /// The first snapshot's scenario id (the base for a batch).
    scenario: i64,
    /// Number of stacked scenarios (1 for a single case).
    n_scenarios: usize,
    schema: &'static str,
    /// Shared base element set (equal across all snapshots by the shape check).
    n_buses: usize,
    n_branches: usize,
    /// In-service branch count of the **first** snapshot; branch status may
    /// differ per scenario, so this describes scenario 0, not the whole batch.
    n_branches_in_service: usize,
    n_gens: usize,
    /// Reference (slack) bus of the **first** snapshot. Each snapshot resolves
    /// its own reference and carries it in the bus `REF` / gen `is_slack_gen`
    /// columns; this records scenario 0's.
    reference_bus: usize,
    /// Branches whose impedance is too small to divide by (admittance and
    /// derived flow columns zeroed), summed over every snapshot in the batch.
    dropped_zero_impedance: usize,
    /// Generators whose cost row gridfm can't represent (piecewise, missing,
    /// malformed, or cubic and higher; `cp*` columns zeroed), summed over every
    /// snapshot in the batch.
    degenerate_cost_gens: usize,
    /// Missing cost rows after applying the export cost policy.
    missing_cost_gens: usize,
    /// Present cost rows that gridfm cannot represent.
    unsupported_cost_gens: usize,
    /// Same as `degenerate_cost_gens`; named for the physical `cp*` columns.
    zeroed_cost_gens: usize,
    cost_policy: MissingGenCostPolicy,
    synthesized_gen_costs: usize,
    patched_gen_costs: usize,
    files: Vec<String>,
    powerio_version: String,
}

/// Build the four gridfm tables for one network, stamping `scenario` into the id
/// columns. Pure (no I/O). A thin wrapper over [`to_gridfm_record_batches`]
/// for one snapshot.
///
/// # Errors
/// As [`to_gridfm_record_batches`], minus the batch-only shape errors a single
/// snapshot cannot trigger.
pub fn to_gridfm_record_batches_single(
    net: &BalancedNetwork,
    scenario: i64,
    opts: &GridfmOptions,
) -> Result<GridfmTables> {
    let snap = GridfmSnapshot::new(net, scenario);
    to_gridfm_record_batches(std::slice::from_ref(&snap), opts)
}

/// Build the four gridfm tables for a batch of scenarios, row stacked and keyed
/// by the `scenario` column. Pure (no I/O). Each snapshot carries its own
/// scenario id; the `include_y_bus`/taps/shifts flags apply to every snapshot.
///
/// # Errors
/// [`Error::EmptyScenarioBatch`] for an empty batch,
/// [`Error::ScenarioShapeMismatch`] if the snapshots don't share one base element
/// set and system base or reuse a scenario id,
/// [`powerio_tx::Error::ReferenceBusCount`] unless each case has exactly one reference bus
/// (graphkit needs a slack), [`Error::NormalizedGridfmSnapshot`] for a normalized
/// input, [`Error::NonFiniteGridfmValue`] for a NaN/Inf physical quantity (or a
/// NaN limit — `±Inf` limits are the valid unbounded sentinel) that would reach
/// Parquet, [`powerio_tx::Error::NonFiniteSusceptance`] if a finite branch impedance still
/// yields a non-finite admittance, and [`powerio_tx::Error::UnknownBus`] if a generator or
/// branch references a bus the network doesn't define.
pub fn to_gridfm_record_batches(
    snapshots: &[GridfmSnapshot],
    opts: &GridfmOptions,
) -> Result<GridfmTables> {
    if opts.cost_options_default() {
        let views = snapshot_views(snapshots)?;
        return tables_from_views(&views, opts);
    }

    let (nets, _) = policy_adjusted_snapshots(snapshots, opts)?;
    let adjusted: Vec<_> = nets
        .iter()
        .zip(snapshots)
        .map(|(net, snap)| GridfmSnapshot::new(net, snap.scenario))
        .collect();
    let views = snapshot_views(&adjusted)?;
    tables_from_views(&views, opts)
}

/// The four tables from already-built, shape-checked snapshot views. The Y_bus
/// table is built only when [`GridfmOptions::include_y_bus`] is set — otherwise
/// it's `None` and the per-snapshot `calc_admittance_matrix` is skipped entirely.
fn tables_from_views(views: &[SnapshotView], opts: &GridfmOptions) -> Result<GridfmTables> {
    Ok(GridfmTables {
        bus: bus_batch(views)?,
        generator: gen_batch(views)?,
        branch: branch_batch(views, opts)?,
        y_bus: if opts.include_y_bus {
            Some(y_bus_batch(views, opts)?)
        } else {
            None
        },
    })
}

/// A resolved snapshot: its indexed view, scenario id, and reference bus.
struct SnapshotView<'a> {
    view: IndexedNetwork<'a>,
    scenario: i64,
    ref_bus: usize,
}

/// Build and shape-check the views for a scenario batch. Every snapshot must
/// resolve to exactly one reference bus and share the first snapshot's base
/// element set (bus / branch / generator counts and row identities) and system
/// base, so the row-stacked tables stay schema-consistent. Scenario ids must be
/// unique because they key every emitted row.
fn snapshot_views<'a>(snapshots: &'a [GridfmSnapshot<'a>]) -> Result<Vec<SnapshotView<'a>>> {
    let first = snapshots.first().ok_or(Error::EmptyScenarioBatch)?;
    let expected = shape_of(first.net);
    let mut scenario_indices = std::collections::HashMap::with_capacity(snapshots.len());

    let mut views = Vec::with_capacity(snapshots.len());
    for (k, snap) in snapshots.iter().enumerate() {
        validate_snapshot_inputs(snap.net, snap.scenario)?;
        if let Some(first_index) = scenario_indices.insert(snap.scenario, k) {
            return Err(Error::ScenarioShapeMismatch {
                index: k,
                reason: ScenarioMismatch::DuplicateScenarioId {
                    scenario: snap.scenario,
                    first_index,
                },
            });
        }
        let got = shape_of(snap.net);
        if got != expected {
            return Err(Error::ScenarioShapeMismatch {
                index: k,
                reason: ScenarioMismatch::Counts { expected, got },
            });
        }
        let ids_match = snap.net.buses().iter().map(|b| b.id).eq(first
            .net
            .buses()
            .iter()
            .map(|b| b.id));
        if !ids_match {
            return Err(Error::ScenarioShapeMismatch {
                index: k,
                reason: ScenarioMismatch::BusOrder,
            });
        }
        if !snap
            .net
            .branches()
            .iter()
            .zip(first.net.branches())
            .all(|(got, expected)| branch_identity_matches(got, expected))
        {
            return Err(Error::ScenarioShapeMismatch {
                index: k,
                reason: ScenarioMismatch::BranchOrder,
            });
        }
        if !snap
            .net
            .generators()
            .iter()
            .zip(first.net.generators())
            .all(|(got, expected)| generator_identity_matches(got, expected))
        {
            return Err(Error::ScenarioShapeMismatch {
                index: k,
                reason: ScenarioMismatch::GeneratorOrder,
            });
        }
        if snap.net.base_mva().to_bits() != first.net.base_mva().to_bits() {
            return Err(Error::ScenarioShapeMismatch {
                index: k,
                reason: ScenarioMismatch::BaseMva,
            });
        }
        let view = IndexedNetwork::new(snap.net);
        let ref_bus = view.reference_bus_index()?;
        views.push(SnapshotView {
            view,
            scenario: snap.scenario,
            ref_bus,
        });
    }
    Ok(views)
}

fn branch_identity_matches(got: &powerio_tx::Branch, expected: &powerio_tx::Branch) -> bool {
    got.from == expected.from
        && got.to == expected.to
        && got.uid.as_deref() == expected.uid.as_deref()
}

fn generator_identity_matches(
    got: &powerio_tx::Generator,
    expected: &powerio_tx::Generator,
) -> bool {
    got.bus == expected.bus && got.uid.as_deref() == expected.uid.as_deref()
}

fn validate_snapshot_inputs(net: &BalancedNetwork, scenario: i64) -> Result<()> {
    if net.is_normalized() {
        return Err(Error::NormalizedGridfmSnapshot { scenario });
    }
    net.check_base_mva()?;

    // Two tiers: physical quantities must be finite, while limit fields admit
    // `±Inf` as the unbounded sentinel (the PowerModels reader defaults absent
    // gen limits to `±Inf`, MATPOWER files carry literal `Inf` tokens, and
    // Parquet stores `±Inf` natively) — only `NaN` is corrupt there.
    for (row, b) in net.buses().iter().enumerate() {
        finite(scenario, "bus", row, "vm", b.vm)?;
        finite(scenario, "bus", row, "va", b.va)?;
        finite(scenario, "bus", row, "base_kv", b.base_kv)?;
        not_nan(scenario, "bus", row, "vmax", b.vmax)?;
        not_nan(scenario, "bus", row, "vmin", b.vmin)?;
    }
    for (row, l) in net.loads().iter().enumerate() {
        finite(scenario, "load", row, "p", l.p)?;
        finite(scenario, "load", row, "q", l.q)?;
    }
    for (row, s) in net.shunts().iter().enumerate() {
        finite(scenario, "shunt", row, "g", s.g)?;
        finite(scenario, "shunt", row, "b", s.b)?;
    }
    for (row, br) in net.branches().iter().enumerate() {
        finite(scenario, "branch", row, "r", br.r)?;
        finite(scenario, "branch", row, "x", br.x)?;
        finite(scenario, "branch", row, "b", br.calc_total_charging_b())?;
        finite(scenario, "branch", row, "tap", br.tap)?;
        finite(scenario, "branch", row, "shift", br.shift)?;
        not_nan(scenario, "branch", row, "angmin", br.angmin)?;
        not_nan(scenario, "branch", row, "angmax", br.angmax)?;
        not_nan(scenario, "branch", row, "rate_a", br.rate_a)?;
        if let Some(solution) = br.solution {
            finite(scenario, "branch", row, "pf", solution.pf)?;
            finite(scenario, "branch", row, "qf", solution.qf)?;
            finite(scenario, "branch", row, "pt", solution.pt)?;
            finite(scenario, "branch", row, "qt", solution.qt)?;
        }
    }
    for (row, g) in net.generators().iter().enumerate() {
        finite(scenario, "generator", row, "pg", g.pg)?;
        finite(scenario, "generator", row, "qg", g.qg)?;
        not_nan(scenario, "generator", row, "pmax", g.pmax)?;
        not_nan(scenario, "generator", row, "pmin", g.pmin)?;
        not_nan(scenario, "generator", row, "qmax", g.qmax)?;
        not_nan(scenario, "generator", row, "qmin", g.qmin)?;
        // Validate exactly what lands in the `cp0/cp1/cp2` columns, not the
        // raw coefficient list (non-representable rows export as zeros).
        let (cp0, cp1, cp2) = gridfm_cost(g.cost.as_ref());
        finite(scenario, "gencost", row, "cp0", cp0)?;
        finite(scenario, "gencost", row, "cp1", cp1)?;
        finite(scenario, "gencost", row, "cp2", cp2)?;
    }
    Ok(())
}

fn finite(
    scenario: i64,
    element: &'static str,
    row: usize,
    field: &'static str,
    value: f64,
) -> Result<()> {
    if value.is_finite() {
        Ok(())
    } else {
        Err(Error::NonFiniteGridfmValue {
            scenario,
            element,
            row,
            field,
            value,
        })
    }
}

/// Limit fields: `±Inf` is the valid unbounded sentinel, only `NaN` is corrupt.
fn not_nan(
    scenario: i64,
    element: &'static str,
    row: usize,
    field: &'static str,
    value: f64,
) -> Result<()> {
    if value.is_nan() {
        Err(Error::NonFiniteGridfmValue {
            scenario,
            element,
            row,
            field,
            value,
        })
    } else {
        Ok(())
    }
}

/// The base element shape a scenario batch shares.
fn shape_of(net: &BalancedNetwork) -> ElementCounts {
    ElementCounts {
        buses: net.buses().len(),
        branches: net.branches().len(),
        gens: net.generators().len(),
    }
}

#[derive(Debug, Clone, Copy, Default)]
struct CostPolicyBatchReport {
    synthesized: usize,
    patched: usize,
}

fn policy_adjusted_snapshots(
    snapshots: &[GridfmSnapshot],
    opts: &GridfmOptions,
) -> Result<(Vec<BalancedNetwork>, CostPolicyBatchReport)> {
    let mut report = CostPolicyBatchReport::default();
    let mut nets = Vec::with_capacity(snapshots.len());
    for snap in snapshots {
        let mut net = snap.net.clone();
        let r = net.apply_gen_cost_policy(&opts.gen_cost_patches, opts.missing_gen_cost)?;
        report.synthesized += r.synthesized;
        report.patched += r.patched;
        nets.push(net);
    }
    Ok((nets, report))
}

/// Number a list of networks into snapshots, stamping the k-th `base + k` — the
/// one place the "k-th input is scenario `base + k`" rule lives, so the CLI and
/// the Python binding can't drift. Checked: returns [`Error::ScenarioIdOverflow`]
/// rather than wrapping or panicking if a scenario id exceeds `i64`.
pub fn number_snapshots<'a>(
    nets: &[&'a BalancedNetwork],
    base: i64,
) -> Result<Vec<GridfmSnapshot<'a>>> {
    nets.iter()
        .enumerate()
        .map(|(k, &net)| {
            let scenario = i64::try_from(k)
                .ok()
                .and_then(|offset| base.checked_add(offset))
                .ok_or(Error::ScenarioIdOverflow { base, index: k })?;
            Ok(GridfmSnapshot::new(net, scenario))
        })
        .collect()
}

/// Write the gridfm-datakit Parquet dataset for one case under
/// `out_dir/<sanitized network name>/raw/` (see
/// [`sanitize_stem`](crate::sanitize_stem)), matching datakit's directory
/// layout. Stamps
/// `scenario` into the id columns. Writes `bus_data.parquet`, `gen_data.parquet`,
/// `branch_data.parquet`, optionally `y_bus_data.parquet`, and a
/// `gridfm_meta.json` manifest.
///
/// Expects a raw snapshot (powers in MW, angles in degrees); pass the parsed
/// `BalancedNetwork`, not a [`to_normalized`](powerio_tx::BalancedNetwork::to_normalized) per-unit
/// product, whose fields would be mislabeled.
///
/// # Errors
/// Propagates [`to_gridfm_record_batches_single`] and any filesystem/Parquet error.
pub fn emit_gridfm_dataset(
    net: &BalancedNetwork,
    scenario: i64,
    out_dir: impl AsRef<Path>,
    opts: &GridfmOptions,
) -> Result<GridfmOutputs> {
    let snap = GridfmSnapshot::new(net, scenario);
    emit_gridfm_batch(std::slice::from_ref(&snap), out_dir, opts)
}

/// Write a batch of scenarios as one gridfm-datakit dataset under
/// `out_dir/<sanitized network name>/raw/` (see
/// [`sanitize_stem`](crate::sanitize_stem)), row-stacking every snapshot's
/// tables and keying
/// them by the `scenario` column. The dataset name and the base element counts
/// come from the first snapshot (shared across the batch by the shape check); the
/// dropped/degenerate counts are summed over every snapshot, while `reference_bus`
/// / `n_branches_in_service` record the first snapshot only (they can differ per
/// scenario, so the manifest documents them as scenario 0's).
///
/// # Errors
/// Propagates [`to_gridfm_record_batches`] and any filesystem/Parquet error.
pub fn emit_gridfm_batch(
    snapshots: &[GridfmSnapshot],
    out_dir: impl AsRef<Path>,
    opts: &GridfmOptions,
) -> Result<GridfmOutputs> {
    let case_name = snapshots
        .first()
        .ok_or(Error::EmptyScenarioBatch)?
        .net
        .name();
    let case_root = out_dir.as_ref().join(crate::sanitize_stem(case_name));
    let dataset = build_gridfm_batch(snapshots, opts)?;
    let GridfmDataset {
        artifacts,
        diagnostics,
        dropped_zero_impedance,
        degenerate_cost_gens,
        missing_cost_gens,
        unsupported_cost_gens,
        synthesized_gen_costs,
        patched_gen_costs,
    } = dataset;
    let committed = powerio_core::Destination::path(&case_root).__commit_artifacts(
        true,
        powerio_core::Fidelity::Canonical,
        artifacts,
        Vec::new(),
    )?;
    let powerio_core::EmittedOutput::Path { root, artifacts } = committed.into_output() else {
        unreachable!("a path destination returns a path output")
    };

    Ok(GridfmOutputs {
        dir: root.join("raw"),
        files: artifacts,
        dropped_zero_impedance,
        degenerate_cost_gens,
        missing_cost_gens,
        unsupported_cost_gens,
        synthesized_gen_costs,
        patched_gen_costs,
        diagnostics,
    })
}

/// Build the GridFM artifacts for one balanced network without writing them.
///
/// The returned inventory can be committed through the same PowerIO
/// destination operation used by every other directory format.
pub fn build_gridfm_dataset(
    net: &BalancedNetwork,
    scenario: i64,
    opts: &GridfmOptions,
) -> Result<GridfmDataset> {
    let snapshot = GridfmSnapshot::new(net, scenario);
    build_gridfm_batch(std::slice::from_ref(&snapshot), opts)
}

/// Build one GridFM artifact inventory for a batch of numbered snapshots.
pub fn build_gridfm_batch(
    snapshots: &[GridfmSnapshot],
    opts: &GridfmOptions,
) -> Result<GridfmDataset> {
    if opts.cost_options_default() {
        return build_gridfm_batch_inner(snapshots, opts, CostPolicyBatchReport::default());
    }

    let (nets, cost_report) = policy_adjusted_snapshots(snapshots, opts)?;
    let adjusted: Vec<_> = nets
        .iter()
        .zip(snapshots)
        .map(|(net, snap)| GridfmSnapshot::new(net, snap.scenario))
        .collect();
    build_gridfm_batch_inner(&adjusted, opts, cost_report)
}

fn build_gridfm_batch_inner(
    snapshots: &[GridfmSnapshot],
    opts: &GridfmOptions,
    cost_report: CostPolicyBatchReport,
) -> Result<GridfmDataset> {
    let views = snapshot_views(snapshots)?;
    let tables = tables_from_views(&views, opts)?;

    // The shape check guarantees every snapshot shares the base element set, so
    // the name and structural counts come from the first.
    let net = views[0].view.network();
    let mut inventory: Vec<(&'static str, Vec<u8>)> = Vec::new();
    inventory.push(("bus_data.parquet", parquet_bytes(&tables.bus)?));
    inventory.push(("gen_data.parquet", parquet_bytes(&tables.generator)?));
    inventory.push(("branch_data.parquet", parquet_bytes(&tables.branch)?));
    if let Some(y_bus) = &tables.y_bus {
        inventory.push(("y_bus_data.parquet", parquet_bytes(y_bus)?));
    }

    // Branch status and costs may differ per scenario, so count the zeroed rows
    // across every snapshot — the totals describe the whole stacked dataset.
    let dropped_zero_impedance: usize = views
        .iter()
        .flat_map(|v| v.view.network().branches().iter())
        .filter(|br| br.r.hypot(br.x) < powerio_tx::dc::MIN_DIVISIBLE_MAGNITUDE)
        .count();
    let missing_cost_gens: usize = views
        .iter()
        .flat_map(|v| v.view.network().generators().iter())
        .filter(|g| g.cost.is_none())
        .count();
    let unsupported_cost_gens: usize = views
        .iter()
        .flat_map(|v| v.view.network().generators().iter())
        .filter(|g| g.cost.is_some() && !cost_representable(g.cost.as_ref()))
        .count();
    let degenerate_cost_gens = missing_cost_gens + unsupported_cost_gens;
    let diagnostics =
        projection_diagnostics(&views, opts, missing_cost_gens, unsupported_cost_gens);

    let meta = GridfmMeta {
        case_name: net.name().clone(),
        base_mva: net.base_mva(),
        scenario: views[0].scenario,
        n_scenarios: views.len(),
        schema: "gridfm-datakit",
        n_buses: net.buses().len(),
        n_branches: net.branches().len(),
        n_branches_in_service: net.branches().iter().filter(|b| b.in_service).count(),
        n_gens: net.generators().len(),
        reference_bus: views[0].ref_bus,
        dropped_zero_impedance,
        degenerate_cost_gens,
        missing_cost_gens,
        unsupported_cost_gens,
        zeroed_cost_gens: degenerate_cost_gens,
        cost_policy: opts.missing_gen_cost,
        synthesized_gen_costs: cost_report.synthesized,
        patched_gen_costs: cost_report.patched,
        // The manifest lists the table files it describes; it does not list
        // itself, matching the manifest shape consumers already read.
        files: inventory
            .iter()
            .map(|(name, _)| (*name).to_string())
            .collect(),
        powerio_version: env!("CARGO_PKG_VERSION").to_string(),
    };
    let json = serde_json::to_string_pretty(&meta).map_err(|e| Error::Parquet(e.to_string()))?;
    inventory.push(("gridfm_meta.json", json.into_bytes()));

    let artifacts = inventory
        .into_iter()
        .map(|(name, bytes)| {
            Ok(powerio_core::MemoryArtifact::new(
                powerio_core::ArtifactPath::new(format!("raw/{name}"))?,
                bytes,
            ))
        })
        .collect::<std::result::Result<Vec<_>, powerio_core::Error>>()?;
    Ok(GridfmDataset {
        artifacts,
        diagnostics,
        dropped_zero_impedance,
        degenerate_cost_gens,
        missing_cost_gens,
        unsupported_cost_gens,
        synthesized_gen_costs: cost_report.synthesized,
        patched_gen_costs: cost_report.patched,
    })
}

#[derive(Default)]
struct ProjectionCounts {
    renumbered_buses: usize,
    unsupported_bus_types: usize,
    bus_metadata: usize,
    collapsed_load_buses: usize,
    omitted_load_records: usize,
    load_metadata: usize,
    collapsed_shunt_buses: usize,
    omitted_shunt_records: usize,
    shunt_metadata: usize,
    branch_metadata: usize,
    branch_extra_ratings: usize,
    branch_charging: usize,
    unit_transformers: usize,
    derived_branch_flows: usize,
    zero_impedance_admittances: usize,
    generator_metadata: usize,
    generator_voltage: usize,
    generator_mbase: usize,
    reshaped_costs: usize,
    omitted_record_families: Vec<(&'static str, usize)>,
    network_metadata_snapshots: usize,
    base_frequency_snapshots: usize,
    taps_omitted_from_admittance: usize,
    shifts_omitted_from_admittance: usize,
}

#[allow(clippy::too_many_lines)] // one message mapping per GridFM projection count
fn projection_diagnostics(
    views: &[SnapshotView],
    opts: &GridfmOptions,
    missing_cost_gens: usize,
    unsupported_cost_gens: usize,
) -> Vec<powerio_core::Diagnostic> {
    let mut counts = ProjectionCounts::default();
    for snapshot in views {
        collect_projection_counts(snapshot, opts, &mut counts);
    }

    let mut diagnostics = Vec::new();
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_SUBSTITUTED,
        counts.renumbered_buses,
        "bus row(s) use source ids that GridFM replaces with dense zero-based bus indices",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_SUBSTITUTED,
        counts.unsupported_bus_types,
        "isolated bus row(s) have no GridFM bus-type column and read back as PQ",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
        counts.bus_metadata,
        "bus row(s) carry area, zone, name, uid, location, emergency voltage, or extra fields that GridFM bus_data does not store",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_COLLAPSED,
        counts.collapsed_load_buses,
        "bus row(s) combine several in-service load records into the Pd and Qd nodal totals",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
        counts.omitted_load_records,
        "load record(s) are out of service or zero-valued and therefore have no distinct GridFM row",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
        counts.load_metadata,
        "load record(s) carry a voltage model, uid, or extra fields that GridFM bus_data does not store",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_COLLAPSED,
        counts.collapsed_shunt_buses,
        "bus row(s) combine several in-service shunt records into the GS and BS nodal totals",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
        counts.omitted_shunt_records,
        "shunt record(s) are out of service or zero-valued and therefore have no distinct GridFM row",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
        counts.shunt_metadata,
        "shunt record(s) carry section, control, uid, or extra fields that GridFM bus_data does not store",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
        counts.branch_metadata,
        "branch row(s) carry name, uid, route, control, or extra fields that GridFM branch_data does not store",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
        counts.branch_extra_ratings,
        "branch row(s) carry rate_b, rate_c, named MVA ratings, or current ratings while GridFM stores only rate_a",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_SUBSTITUTED,
        counts.branch_charging,
        "branch row(s) have asymmetric or conductive terminal charging that GridFM replaces with one total susceptance",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_SUBSTITUTED,
        counts.unit_transformers,
        "unity-ratio, zero-shift transformer row(s) are indistinguishable from lines in GridFM branch_data",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_DEFAULTED,
        counts.derived_branch_flows,
        "branch row(s) have no stated solution, so GridFM pf, qf, pt, and qt are evaluated from the stored voltages",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_DEFAULTED,
        counts.zero_impedance_admittances,
        "zero-impedance branch row(s) use zero admittance values because finite Y entries cannot be derived",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
        counts.generator_metadata,
        "generator row(s) carry energy source, capability, regulation, active-power control, or uid fields that GridFM gen_data does not store",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_SUBSTITUTED,
        counts.generator_voltage,
        "generator voltage setpoint(s) differ from the connected bus Vm that GridFM uses on read",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_SUBSTITUTED,
        counts.generator_mbase,
        "generator mbase value(s) differ from the system base_mva that GridFM uses on read",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_DEFAULTED,
        missing_cost_gens,
        "generator row(s) have no cost, so the required GridFM cp0, cp1, and cp2 values are written as zero",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_SUBSTITUTED,
        unsupported_cost_gens,
        "generator cost row(s) are not polynomial curves of degree at most two and are replaced by zero GridFM coefficients",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_SUBSTITUTED,
        counts.reshaped_costs,
        "representable generator cost row(s) lose startup, shutdown, or their original coefficient count in GridFM's fixed quadratic columns",
    );
    for (family, count) in counts.omitted_record_families {
        diagnostics.push(powerio_core::Diagnostic::of(
            &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
            format!("{count} {family} record(s) have no GridFM table"),
        ));
    }
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
        counts.network_metadata_snapshots,
        "snapshot(s) carry geographic, case, detailed-connectivity, or solver metadata that GridFM does not store",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_FIELD_DROPPED,
        counts.base_frequency_snapshots,
        "snapshot(s) use a nondefault base frequency while GridFM has no base-frequency field",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_SUBSTITUTED,
        counts.taps_omitted_from_admittance,
        "branch row(s) retain tap values but the requested GridFM admittance tables ignore them",
    );
    push_count(
        &mut diagnostics,
        &crate::diagnostics::codes::EMIT_GRIDFM_VALUE_SUBSTITUTED,
        counts.shifts_omitted_from_admittance,
        "branch row(s) retain phase shifts but the requested GridFM admittance tables ignore them",
    );
    diagnostics
}

#[allow(clippy::too_many_lines)] // one count block per GridFM table and omitted record family
fn collect_projection_counts(
    snapshot: &SnapshotView,
    opts: &GridfmOptions,
    counts: &mut ProjectionCounts,
) {
    let net = snapshot.view.network();
    counts.base_frequency_snapshots +=
        usize::from((net.base_frequency() - DEFAULT_BASE_FREQUENCY).abs() > 1e-9);
    counts.network_metadata_snapshots += usize::from(
        net.geo().is_some()
            || net.case_metadata() != &powerio_tx::CaseMetadata::default()
            || net.detailed_connectivity().is_some()
            || net.solver().is_some(),
    );

    let omitted_families = [
        (
            "static VAR compensator",
            net.static_var_compensators().len(),
        ),
        ("switch", net.switches().len()),
        ("storage", net.storage().len()),
        ("HVDC", net.hvdc().len()),
        ("three-winding transformer", net.transformers_3w().len()),
        ("area", net.areas().len()),
    ];
    for (name, count) in omitted_families {
        if count > 0 {
            if let Some((_, total)) = counts
                .omitted_record_families
                .iter_mut()
                .find(|(family, _)| *family == name)
            {
                *total += count;
            } else {
                counts.omitted_record_families.push((name, count));
            }
        }
    }

    let mut load_counts = std::collections::HashMap::new();
    for load in net.loads() {
        if load.in_service && (load.p != 0.0 || load.q != 0.0) {
            *load_counts.entry(load.bus).or_insert(0usize) += 1;
        } else {
            counts.omitted_load_records += 1;
        }
        counts.load_metadata += usize::from(
            load.voltage_model.is_some()
                || (load.uid.is_some() && !net.uid_is_generated(load.uid.as_deref()))
                || !load.extras.is_empty(),
        );
    }
    counts.collapsed_load_buses += load_counts.values().filter(|&&count| count > 1).count();

    let mut shunt_counts = std::collections::HashMap::new();
    for shunt in net.shunts() {
        if shunt.in_service && (shunt.g != 0.0 || shunt.b != 0.0) {
            *shunt_counts.entry(shunt.bus).or_insert(0usize) += 1;
        } else {
            counts.omitted_shunt_records += 1;
        }
        counts.shunt_metadata += usize::from(
            shunt.section_count.is_some()
                || shunt.control.is_some()
                || (shunt.uid.is_some() && !net.uid_is_generated(shunt.uid.as_deref()))
                || !shunt.extras.is_empty(),
        );
    }
    counts.collapsed_shunt_buses += shunt_counts.values().filter(|&&count| count > 1).count();

    for (row, bus) in net.buses().iter().enumerate() {
        counts.renumbered_buses += usize::from(bus.id != BusId(row + 1));
        counts.unsupported_bus_types += usize::from(bus.kind == BusType::Isolated);
        counts.bus_metadata += usize::from(
            bus.area != 0
                || bus.zone != 0
                || bus.name.is_some()
                || (bus.uid.is_some() && !net.uid_is_generated(bus.uid.as_deref()))
                || bus.location.is_some()
                || bus.evhi.is_some()
                || bus.evlo.is_some()
                || !bus.extras.is_empty(),
        );
    }

    for branch in net.branches() {
        counts.branch_metadata += usize::from(
            branch.name.is_some()
                || branch.control.is_some()
                || (branch.uid.is_some() && !net.uid_is_generated(branch.uid.as_deref()))
                || branch.route.is_some()
                || !branch.extras.is_empty(),
        );
        counts.branch_extra_ratings += usize::from(
            branch.rate_b != 0.0
                || branch.rate_c != 0.0
                || !branch.rating_sets.is_empty()
                || branch.current_ratings.is_some(),
        );
        counts.branch_charging += usize::from(branch.has_non_matpower_charging());
        counts.unit_transformers +=
            usize::from(is_exact_unit_transformer(branch.tap, branch.shift));
        counts.derived_branch_flows += usize::from(branch.solution.is_none());
        counts.zero_impedance_admittances += usize::from(
            branch.in_service && branch.r.hypot(branch.x) < powerio_tx::dc::MIN_DIVISIBLE_MAGNITUDE,
        );
        counts.taps_omitted_from_admittance +=
            usize::from(!opts.include_taps && differs_from_exact_one(branch.calc_effective_tap()));
        counts.shifts_omitted_from_admittance +=
            usize::from(!opts.include_shifts && branch.shift != 0.0);
    }

    for generator in net.generators() {
        counts.generator_metadata += usize::from(
            generator.energy_source != powerio_tx::GeneratorEnergySource::Other
                || generator.has_caps()
                || !generator.voltage_regulation_on
                || generator.regulating_terminal.is_some()
                || generator.regulated_bus.is_some()
                || generator.active_power_control.is_some()
                || (generator.uid.is_some() && !net.uid_is_generated(generator.uid.as_deref())),
        );
        if let Some(bus) = snapshot
            .view
            .bus_index(generator.bus)
            .map(|index| &net.buses()[index])
        {
            counts.generator_voltage += usize::from((generator.vg - bus.vm).abs() > 1e-12);
        }
        counts.generator_mbase += usize::from((generator.mbase - net.base_mva()).abs() > 1e-9);
        counts.reshaped_costs += usize::from(generator.cost.as_ref().is_some_and(|cost| {
            cost_representable(Some(cost))
                && (cost.ncost != 3
                    || cost.coeffs.len() != 3
                    || cost.startup != 0.0
                    || cost.shutdown != 0.0)
        }));
    }
}

// Tap and shift sentinels classify an explicit unity transformer exactly.
#[allow(clippy::float_cmp)]
fn is_exact_unit_transformer(tap: f64, shift: f64) -> bool {
    tap == 1.0 && shift == 0.0
}

// The admittance changes for every effective tap except the exact identity ratio.
#[allow(clippy::float_cmp)]
fn differs_from_exact_one(value: f64) -> bool {
    value != 1.0
}

fn push_count(
    diagnostics: &mut Vec<powerio_core::Diagnostic>,
    code: &'static powerio_core::DiagnosticInfo,
    count: usize,
    message: &'static str,
) {
    if count > 0 {
        diagnostics.push(powerio_core::Diagnostic::of(
            code,
            format!("{count} {message}"),
        ));
    }
}

// --- table builders --------------------------------------------------------

fn bus_batch(snaps: &[SnapshotView]) -> Result<RecordBatch> {
    let total: usize = snaps.iter().map(|s| s.view.n()).sum();
    let mut scenario = Vec::with_capacity(total);
    let mut bus_idx = Vec::with_capacity(total);
    let (mut pd, mut qd) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut pg_col, mut qg_col) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut vm, mut va) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut pq, mut pv, mut refc) = (
        Vec::with_capacity(total),
        Vec::with_capacity(total),
        Vec::with_capacity(total),
    );
    let mut vn_kv = Vec::with_capacity(total);
    let (mut min_vm, mut max_vm) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut gs, mut bs) = (Vec::with_capacity(total), Vec::with_capacity(total));

    for s in snaps {
        let view = &s.view;
        let n = view.n();
        let base = view.base_mva();
        let buses = &view.network().buses();

        // Per-bus generation, summed over in-service generators (dense order).
        let mut pg = vec![0.0; n];
        let mut qg = vec![0.0; n];
        for (_, g) in view.in_service_gens() {
            if let Some(i) = view.bus_index(g.bus) {
                pg[i] += g.pg;
                qg[i] += g.qg;
            }
        }

        scenario.resize(scenario.len() + n, s.scenario);
        bus_idx.extend(0..n as i64);
        pd.extend_from_slice(view.pd());
        qd.extend_from_slice(view.qd());
        pg_col.extend(pg);
        qg_col.extend(qg);
        vm.extend(buses.iter().map(|b| b.vm));
        va.extend(buses.iter().map(|b| b.va));
        pq.extend(buses.iter().map(|b| i64::from(b.kind == BusType::Pq)));
        pv.extend(buses.iter().map(|b| i64::from(b.kind == BusType::Pv)));
        refc.extend(buses.iter().map(|b| i64::from(b.kind == BusType::Ref)));
        vn_kv.extend(buses.iter().map(|b| b.base_kv));
        min_vm.extend(buses.iter().map(|b| b.vmin));
        max_vm.extend(buses.iter().map(|b| b.vmax));
        gs.extend(view.gs().iter().map(|g| g / base));
        bs.extend(view.bs().iter().map(|b| b / base));
    }

    batch(with_scenario_pair(
        scenario,
        vec![
            ("bus", i64s(bus_idx)),
            ("Pd", f64s(pd)),
            ("Qd", f64s(qd)),
            ("Pg", f64s(pg_col)),
            ("Qg", f64s(qg_col)),
            ("Vm", f64s(vm)),
            ("Va", f64s(va)),
            ("PQ", i64s(pq)),
            ("PV", i64s(pv)),
            ("REF", i64s(refc)),
            ("vn_kv", f64s(vn_kv)),
            ("min_vm_pu", f64s(min_vm)),
            ("max_vm_pu", f64s(max_vm)),
            ("GS", f64s(gs)),
            ("BS", f64s(bs)),
        ],
    ))
}

fn gen_batch(snaps: &[SnapshotView]) -> Result<RecordBatch> {
    let total: usize = snaps.iter().map(|s| s.view.generators().len()).sum();
    let mut scenario = Vec::with_capacity(total);
    let mut idx = Vec::with_capacity(total);
    let mut bus = Vec::with_capacity(total);
    let (mut p_mw, mut q_mvar) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut min_p, mut max_p) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut min_q, mut max_q) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut cp0, mut cp1, mut cp2) = (
        Vec::with_capacity(total),
        Vec::with_capacity(total),
        Vec::with_capacity(total),
    );
    let mut in_service = Vec::with_capacity(total);
    let mut is_slack = Vec::with_capacity(total);

    for s in snaps {
        let view = &s.view;
        // One pass over the snapshot's generators: every column gets one push per
        // generator, in dense source order.
        for (row, g) in view.generators().iter().enumerate() {
            let i = view.bus_index(g.bus).ok_or(powerio_tx::Error::UnknownBus {
                bus_id: g.bus,
                element_index: row,
            })?;
            scenario.push(s.scenario);
            idx.push(row as i64);
            bus.push(i as i64);
            is_slack.push(i64::from(i == s.ref_bus));
            let (c0, c1, c2) = gridfm_cost(g.cost.as_ref());
            cp0.push(c0);
            cp1.push(c1);
            cp2.push(c2);
            p_mw.push(g.pg);
            q_mvar.push(g.qg);
            min_p.push(g.pmin);
            max_p.push(g.pmax);
            min_q.push(g.qmin);
            max_q.push(g.qmax);
            in_service.push(i64::from(g.in_service));
        }
    }

    batch(with_scenario_pair(
        scenario,
        vec![
            ("idx", i64s(idx)),
            ("bus", i64s(bus)),
            ("p_mw", f64s(p_mw)),
            ("q_mvar", f64s(q_mvar)),
            ("min_p_mw", f64s(min_p)),
            ("max_p_mw", f64s(max_p)),
            ("min_q_mvar", f64s(min_q)),
            ("max_q_mvar", f64s(max_q)),
            ("cp0_eur", f64s(cp0)),
            ("cp1_eur_per_mw", f64s(cp1)),
            ("cp2_eur_per_mw2", f64s(cp2)),
            ("in_service", i64s(in_service)),
            ("is_slack_gen", i64s(is_slack)),
        ],
    ))
}

#[allow(clippy::too_many_lines, clippy::many_single_char_names)]
fn branch_batch(snaps: &[SnapshotView], opts: &GridfmOptions) -> Result<RecordBatch> {
    let total: usize = snaps.iter().map(|s| s.view.branches().len()).sum();

    // Same flags the Y_bus builder derives, so the branch admittance columns and
    // y_bus_data come from one kernel. The taps/shifts flags are batch-wide (from
    // `opts`); the admittances themselves are recomputed per snapshot. Every
    // branch gets a row regardless of impedance, so a zero impedance branch
    // is zeroed here rather than erroring; the caller counts the drop.
    let flags = YbusFlags {
        unity_taps: !opts.include_taps,
        zero_shifts: !opts.include_shifts,
        skip_zero_impedance: true,
        ..Default::default()
    };

    let mut scenario = Vec::with_capacity(total);
    let mut idx = Vec::with_capacity(total);
    let (mut from_bus, mut to_bus) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut pf, mut qf, mut pt, mut qt) = (
        Vec::with_capacity(total),
        Vec::with_capacity(total),
        Vec::with_capacity(total),
        Vec::with_capacity(total),
    );
    let (mut yff_r, mut yff_i) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut yft_r, mut yft_i) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut ytf_r, mut ytf_i) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut ytt_r, mut ytt_i) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut r_col, mut x_col, mut b_col) = (
        Vec::with_capacity(total),
        Vec::with_capacity(total),
        Vec::with_capacity(total),
    );
    let (mut tap, mut shift) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let (mut ang_min, mut ang_max) = (Vec::with_capacity(total), Vec::with_capacity(total));
    let mut rate_a = Vec::with_capacity(total);
    let mut br_status = Vec::with_capacity(total);

    for s in snaps {
        let view = &s.view;
        let base = view.base_mva();
        let branches = view.branches();
        let buses = &view.network().buses();
        // Complex bus voltages `vm·e^{jθ}`, dense order, for the flow evaluation.
        let v: Vec<Complex64> = buses
            .iter()
            .map(|b| Complex64::from_polar(b.vm, b.va.to_radians()))
            .collect();

        scenario.resize(scenario.len() + branches.len(), s.scenario);
        idx.extend(0..branches.len() as i64);

        for (row, br) in branches.iter().enumerate() {
            let i = view
                .bus_index(br.from)
                .ok_or(powerio_tx::Error::UnknownBus {
                    bus_id: br.from,
                    element_index: row,
                })?;
            let j = view.bus_index(br.to).ok_or(powerio_tx::Error::UnknownBus {
                bus_id: br.to,
                element_index: row,
            })?;
            from_bus.push(i as i64);
            to_bus.push(j as i64);

            // Zero-impedance branch → `None` → zeroed admittance/flow columns (never NaN).
            let shift_rad = if flags.zero_shifts {
                0.0
            } else {
                view.to_radians(br.shift)
            };
            let block = branch_admittance(br, flags, shift_rad, row)?;
            let [y_ff, y_ft, y_tf, y_tt] = block.unwrap_or([Complex64::new(0.0, 0.0); 4]);
            yff_r.push(y_ff.re);
            yff_i.push(y_ff.im);
            yft_r.push(y_ft.re);
            yft_i.push(y_ft.im);
            ytf_r.push(y_tf.re);
            ytf_i.push(y_tf.im);
            ytt_r.push(y_tt.re);
            ytt_i.push(y_tt.im);

            let (sf, st) = if let Some(solution) = br.solution {
                (
                    Complex64::new(solution.pf / base, solution.qf / base),
                    Complex64::new(solution.pt / base, solution.qt / base),
                )
            } else if br.in_service && block.is_some() {
                branch_flows(&[y_ff, y_ft, y_tf, y_tt], v[i], v[j])
            } else {
                (Complex64::new(0.0, 0.0), Complex64::new(0.0, 0.0))
            };
            pf.push(sf.re * base);
            qf.push(sf.im * base);
            pt.push(st.re * base);
            qt.push(st.im * base);

            r_col.push(br.r);
            x_col.push(br.x);
            b_col.push(br.calc_total_charging_b());
            tap.push(br.calc_effective_tap());
            shift.push(br.shift);
            ang_min.push(br.angmin);
            ang_max.push(br.angmax);
            rate_a.push(br.rate_a);
            br_status.push(i64::from(br.in_service));
        }
    }

    batch(with_scenario_pair(
        scenario,
        vec![
            ("idx", i64s(idx)),
            ("from_bus", i64s(from_bus)),
            ("to_bus", i64s(to_bus)),
            ("pf", f64s(pf)),
            ("qf", f64s(qf)),
            ("pt", f64s(pt)),
            ("qt", f64s(qt)),
            ("r", f64s(r_col)),
            ("x", f64s(x_col)),
            ("b", f64s(b_col)),
            ("Yff_r", f64s(yff_r)),
            ("Yff_i", f64s(yff_i)),
            ("Yft_r", f64s(yft_r)),
            ("Yft_i", f64s(yft_i)),
            ("Ytf_r", f64s(ytf_r)),
            ("Ytf_i", f64s(ytf_i)),
            ("Ytt_r", f64s(ytt_r)),
            ("Ytt_i", f64s(ytt_i)),
            ("tap", f64s(tap)),
            ("shift", f64s(shift)),
            ("ang_min", f64s(ang_min)),
            ("ang_max", f64s(ang_max)),
            ("rate_a", f64s(rate_a)),
            ("br_status", i64s(br_status)),
        ],
    ))
}

fn y_bus_batch(snaps: &[SnapshotView], opts: &GridfmOptions) -> Result<RecordBatch> {
    // Upper bound on stacked nnz: each snapshot's Y_bus has at most 4 entries per
    // branch plus a diagonal. The exact count varies (lossless branches, dropped
    // zeros, per-scenario branch status), so this only sizes the allocation.
    let est: usize = snaps
        .iter()
        .map(|s| 4 * s.view.branches().len() + s.view.n())
        .sum();
    let mut scenario = Vec::with_capacity(est);
    let mut index1 = Vec::with_capacity(est);
    let mut index2 = Vec::with_capacity(est);
    let mut g_vals = Vec::with_capacity(est);
    let mut b_vals = Vec::with_capacity(est);

    for s in snaps {
        let parts = calc_admittance_matrix(&s.view, &opts.build_options())?;
        // G and B don't share a sparsity pattern: a lossless branch (r = 0) is a
        // pure reactance, so its G entries are structurally zero where B's aren't.
        // datakit keys y_bus rows on the complex value being nonzero, i.e. the
        // union of the G and B positions. Merge into a sorted (row, col) map so the
        // output is row-major like `np.nonzero`, then drop any all-zero position.
        let mut entries: std::collections::BTreeMap<(usize, usize), (f64, f64)> =
            std::collections::BTreeMap::new();
        for (row, g_row) in parts.g.outer_iterator().enumerate() {
            for (col, &gv) in g_row.iter() {
                entries.entry((row, col)).or_default().0 = gv;
            }
        }
        for (row, b_row) in parts.b.outer_iterator().enumerate() {
            for (col, &bv) in b_row.iter() {
                entries.entry((row, col)).or_default().1 = bv;
            }
        }

        for ((row, col), (gv, bv)) in entries {
            if gv == 0.0 && bv == 0.0 {
                continue;
            }
            scenario.push(s.scenario);
            index1.push(row as i64);
            index2.push(col as i64);
            g_vals.push(gv);
            b_vals.push(bv);
        }
    }

    batch(with_scenario_pair(
        scenario,
        vec![
            ("index1", i64s(index1)),
            ("index2", i64s(index2)),
            ("G", f64s(g_vals)),
            ("B", f64s(b_vals)),
        ],
    ))
}

// --- small helpers ---------------------------------------------------------

/// `(cp0, cp1, cp2)` = raw MATPOWER `(c0, c1, c2)` from a polynomial cost row.
/// Coefficients are highest-order first, so `ncost == 3` is `[c2, c1, c0]`.
/// Piecewise, missing, or malformed rows give zeros.
fn gridfm_cost(cost: Option<&GenCost>) -> (f64, f64, f64) {
    match cost {
        Some(c) if c.model == 2 && c.coeffs.len() >= c.ncost => match c.ncost {
            3 => (c.coeffs[2], c.coeffs[1], c.coeffs[0]),
            2 => (c.coeffs[1], c.coeffs[0], 0.0),
            1 => (c.coeffs[0], 0.0, 0.0),
            _ => (0.0, 0.0, 0.0),
        },
        _ => (0.0, 0.0, 0.0),
    }
}

/// Whether [`gridfm_cost`] can represent this cost row (used for the manifest's
/// degenerate-cost count).
fn cost_representable(cost: Option<&GenCost>) -> bool {
    matches!(cost, Some(c) if c.model == 2 && c.coeffs.len() >= c.ncost && (1..=3).contains(&c.ncost))
}

/// One table as complete Parquet bytes, for the atomic inventory commit.
fn parquet_bytes(batch: &RecordBatch) -> Result<Vec<u8>> {
    let mut bytes = Vec::new();
    let props = WriterProperties::builder()
        .set_compression(Compression::SNAPPY)
        .build();
    let mut writer = ArrowWriter::try_new(&mut bytes, batch.schema(), Some(props))
        .map_err(|e| Error::Parquet(e.to_string()))?;
    writer
        .write(batch)
        .map_err(|e| Error::Parquet(e.to_string()))?;
    writer.close().map_err(|e| Error::Parquet(e.to_string()))?;
    Ok(bytes)
}

/// Assemble a [`RecordBatch`] from named columns, in order; field types are read
/// off the arrays and all columns are non-null.
fn batch(columns: Vec<(&str, ArrayRef)>) -> Result<RecordBatch> {
    let fields: Vec<Field> = columns
        .iter()
        .map(|(name, arr)| Field::new(*name, arr.data_type().clone(), false))
        .collect();
    let arrays: Vec<ArrayRef> = columns.into_iter().map(|(_, arr)| arr).collect();
    RecordBatch::try_new(Arc::new(Schema::new(fields)), arrays)
        .map_err(|e| Error::Parquet(e.to_string()))
}

fn i64s(v: Vec<i64>) -> ArrayRef {
    Arc::new(Int64Array::from(v))
}

fn f64s(v: Vec<f64>) -> ArrayRef {
    Arc::new(Float64Array::from(v))
}

/// Prepend the `scenario` / `load_scenario_idx` id columns (which hold identical
/// values) to a table's other columns. The Int64 array is built once and the Arc
/// shared, so the duplicate column costs a pointer, not a second `Vec`.
fn with_scenario_pair(
    scenario: Vec<i64>,
    rest: Vec<(&'static str, ArrayRef)>,
) -> Vec<(&'static str, ArrayRef)> {
    let scenario = i64s(scenario);
    let mut cols = Vec::with_capacity(rest.len() + 2);
    cols.push(("scenario", scenario.clone()));
    cols.push(("load_scenario_idx", scenario));
    cols.extend(rest);
    cols
}

// === Reader: gridfm dataset to BalancedNetwork (issue #60) =========================
//
// The inverse of the writer above: select one `scenario` out of the GridFM tables
// and rebuild its native balanced representation. `y_bus_data` is derived and
// therefore is not needed to reconstruct the network.
#[cfg(test)]
mod tests {
    use super::*;
    use crate::network::{Branch, BranchCharging, BranchSolution, Bus, BusId, BusType, Generator};
    use arrow::array::{Float64Array, Int64Array};
    use parquet::arrow::arrow_reader::ParquetRecordBatchReaderBuilder;

    const BUS_COLS: &[&str] = &[
        "scenario",
        "load_scenario_idx",
        "bus",
        "Pd",
        "Qd",
        "Pg",
        "Qg",
        "Vm",
        "Va",
        "PQ",
        "PV",
        "REF",
        "vn_kv",
        "min_vm_pu",
        "max_vm_pu",
        "GS",
        "BS",
    ];
    const GEN_COLS: &[&str] = &[
        "scenario",
        "load_scenario_idx",
        "idx",
        "bus",
        "p_mw",
        "q_mvar",
        "min_p_mw",
        "max_p_mw",
        "min_q_mvar",
        "max_q_mvar",
        "cp0_eur",
        "cp1_eur_per_mw",
        "cp2_eur_per_mw2",
        "in_service",
        "is_slack_gen",
    ];
    const BRANCH_COLS: &[&str] = &[
        "scenario",
        "load_scenario_idx",
        "idx",
        "from_bus",
        "to_bus",
        "pf",
        "qf",
        "pt",
        "qt",
        "r",
        "x",
        "b",
        "Yff_r",
        "Yff_i",
        "Yft_r",
        "Yft_i",
        "Ytf_r",
        "Ytf_i",
        "Ytt_r",
        "Ytt_i",
        "tap",
        "shift",
        "ang_min",
        "ang_max",
        "rate_a",
        "br_status",
    ];
    const YBUS_COLS: &[&str] = &[
        "scenario",
        "load_scenario_idx",
        "index1",
        "index2",
        "G",
        "B",
    ];

    fn case14() -> BalancedNetwork {
        let path = concat!(env!("CARGO_MANIFEST_DIR"), "/../tests/data/case14.m");
        powerio_tx::parse(powerio_core::Source::open(path).unwrap())
            .unwrap()
            .into_value()
    }

    fn names(b: &RecordBatch) -> Vec<String> {
        b.schema()
            .fields()
            .iter()
            .map(|f| f.name().clone())
            .collect()
    }

    fn col_i64<'a>(b: &'a RecordBatch, name: &str) -> &'a Int64Array {
        b.column_by_name(name)
            .unwrap()
            .as_any()
            .downcast_ref()
            .unwrap()
    }

    fn col_f64<'a>(b: &'a RecordBatch, name: &str) -> &'a Float64Array {
        b.column_by_name(name)
            .unwrap()
            .as_any()
            .downcast_ref()
            .unwrap()
    }

    fn read(path: &Path) -> RecordBatch {
        let file = std::fs::File::open(path).unwrap();
        let mut reader = ParquetRecordBatchReaderBuilder::try_new(file)
            .unwrap()
            .build()
            .unwrap();
        // case14 fits in one row group / batch.
        reader.next().unwrap().unwrap()
    }

    #[test]
    fn schema_and_row_counts_match_case14() {
        let net = case14();
        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();

        assert_eq!(names(&tables.bus), BUS_COLS);
        assert_eq!(names(&tables.generator), GEN_COLS);
        assert_eq!(names(&tables.branch), BRANCH_COLS);
        assert_eq!(names(tables.y_bus.as_ref().unwrap()), YBUS_COLS);

        assert_eq!(tables.bus.num_rows(), net.buses().len()); // 14
        assert_eq!(tables.generator.num_rows(), net.generators().len()); // 5
        assert_eq!(tables.branch.num_rows(), net.branches().len()); // 20
    }

    #[test]
    fn branch_b_uses_terminal_charging_projection() {
        let mut net = BalancedNetwork::in_memory(
            "terminal-projection",
            100.0,
            vec![bus(1, BusType::Ref), bus(2, BusType::Pq)],
            Vec::new(),
        );
        let mut br = branch(1, 2, 0.01, 0.1);
        br.charging = Some(BranchCharging::new(0.01, 0.02, 0.03, 0.05));
        net.branches_mut().push(br);

        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();
        let b = col_f64(&tables.branch, "b").value(0);
        assert!((b - 0.07).abs() < 1e-12);

        let yff_i = col_f64(&tables.branch, "Yff_i").value(0);
        let ytt_i = col_f64(&tables.branch, "Ytt_i").value(0);
        let y_series_i = -0.1 / (0.01 * 0.01 + 0.1 * 0.1);
        let recovered_b = yff_i + ytt_i - 2.0 * y_series_i;
        assert!((recovered_b - b).abs() < 1e-12);
    }

    #[test]
    fn a_gridfm_write_never_replaces_an_existing_case_directory() {
        let net = case14();

        // A regular file at a produced name inside the case directory: the
        // write refuses and the file keeps its bytes. The dataset root itself
        // stays writable for other cases.
        let dir = tempfile::tempdir().unwrap();
        let case_dir = dir.path().join("case14");
        std::fs::create_dir_all(case_dir.join("raw")).unwrap();
        std::fs::write(case_dir.join("raw/bus_data.parquet"), b"precious").unwrap();
        let error =
            emit_gridfm_dataset(&net, 0, dir.path(), &GridfmOptions::default()).unwrap_err();
        assert!(matches!(error, crate::Error::Commit(_)), "{error:?}");
        assert_eq!(
            std::fs::read(case_dir.join("raw/bus_data.parquet")).unwrap(),
            b"precious"
        );

        // A symbolic link at the case directory name: the link survives and
        // the directory it designates keeps its contents.
        #[cfg(unix)]
        {
            let linked = tempfile::tempdir().unwrap();
            let designated = tempfile::tempdir().unwrap();
            std::fs::write(designated.path().join("keep.txt"), b"kept").unwrap();
            std::os::unix::fs::symlink(designated.path(), linked.path().join("case14")).unwrap();
            let error =
                emit_gridfm_dataset(&net, 0, linked.path(), &GridfmOptions::default()).unwrap_err();
            assert!(matches!(error, crate::Error::Commit(_)), "{error:?}");
            assert!(
                std::fs::symlink_metadata(linked.path().join("case14"))
                    .unwrap()
                    .file_type()
                    .is_symlink()
            );
            assert_eq!(
                std::fs::read(designated.path().join("keep.txt")).unwrap(),
                b"kept"
            );
            assert_eq!(
                std::fs::metadata(designated.path().join("keep.txt"))
                    .unwrap()
                    .len(),
                4
            );
        }

        // The same write into a fresh dataset root still produces the
        // complete inventory.
        let fresh = tempfile::tempdir().unwrap();
        let out = emit_gridfm_dataset(&net, 0, fresh.path(), &GridfmOptions::default()).unwrap();
        let raw = fresh.path().join("case14").join("raw");
        assert_eq!(out.dir, raw);
        for name in [
            "bus_data.parquet",
            "gen_data.parquet",
            "branch_data.parquet",
            "y_bus_data.parquet",
            "gridfm_meta.json",
        ] {
            assert!(raw.join(name).is_file(), "missing {name}");
        }
    }

    #[test]
    fn parquet_round_trips_through_reader() {
        let net = case14();
        let dir = tempfile::tempdir().unwrap();
        let out = emit_gridfm_dataset(&net, 0, dir.path(), &GridfmOptions::default()).unwrap();

        let raw = dir.path().join("case14").join("raw");
        assert_eq!(out.dir, raw);
        for f in ["bus_data", "gen_data", "branch_data", "y_bus_data"] {
            assert!(raw.join(format!("{f}.parquet")).is_file(), "missing {f}");
        }
        assert!(raw.join("gridfm_meta.json").is_file());

        let bus = read(&raw.join("bus_data.parquet"));
        assert_eq!(names(&bus), BUS_COLS);
        assert_eq!(bus.num_rows(), net.buses().len());
        assert_eq!(names(&read(&raw.join("gen_data.parquet"))), GEN_COLS);
        assert_eq!(names(&read(&raw.join("branch_data.parquet"))), BRANCH_COLS);
        assert_eq!(names(&read(&raw.join("y_bus_data.parquet"))), YBUS_COLS);
    }

    #[test]
    fn bus_table_values_are_consistent() {
        let net = case14();
        let view = IndexedNetwork::new(&net);
        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();
        let bus = &tables.bus;

        // Exactly one reference bus; PQ/PV/REF partition every bus.
        let (pq, pv, r) = (col_i64(bus, "PQ"), col_i64(bus, "PV"), col_i64(bus, "REF"));
        assert_eq!(r.values().iter().sum::<i64>(), 1);
        for i in 0..bus.num_rows() {
            assert_eq!(pq.value(i) + pv.value(i) + r.value(i), 1);
        }

        // GS/BS are the per-bus shunt aggregate divided by base_mva.
        let base = net.base_mva();
        let gs = col_f64(bus, "GS");
        for i in 0..bus.num_rows() {
            assert!((gs.value(i) - view.gs()[i] / base).abs() < 1e-12);
        }

        // `bus` column is the dense 0..n range.
        let bus_idx = col_i64(bus, "bus");
        for i in 0..bus.num_rows() {
            assert_eq!(bus_idx.value(i), i as i64);
        }
    }

    #[test]
    fn branch_admittance_columns_match_calculated_admittance_matrix() {
        // The branch table's Y** columns are the same kernel calc_admittance_matrix scatters,
        // so a known in-service branch's block must equal branch_admittance.
        let net = case14();
        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();
        let br = &tables.branch;

        let yff_r = col_f64(br, "Yff_r");
        let yff_i = col_f64(br, "Yff_i");
        for (row, branch) in net.branches().iter().enumerate() {
            // Raw fixture, so the shift is in degrees — convert as calc_admittance_matrix does.
            let shift_rad = branch.shift.to_radians();
            if let Some(block) =
                branch_admittance(branch, YbusFlags::default(), shift_rad, row).unwrap()
            {
                assert!((yff_r.value(row) - block[0].re).abs() < 1e-12);
                assert!((yff_i.value(row) - block[0].im).abs() < 1e-12);
            }
        }
    }

    #[test]
    fn is_slack_gen_marks_the_reference_bus() {
        let net = case14();
        let view = IndexedNetwork::new(&net);
        let ref_bus = view.reference_bus_index().unwrap();
        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();
        let g = &tables.generator;

        let bus = col_i64(g, "bus");
        let slack = col_i64(g, "is_slack_gen");
        for i in 0..g.num_rows() {
            assert_eq!(slack.value(i) == 1, bus.value(i) as usize == ref_bus);
        }
        assert!(slack.values().contains(&1), "no slack generator");
    }

    #[test]
    fn branch_flows_close_the_power_balance_on_a_solved_case() {
        // case14 ships a converged operating point, so the active flows must obey
        // KCL: total branch loss = generation - demand - shunt draw. This is the
        // value-level guard on branch_flows (a missing ×base, a sign flip, or a
        // wrong conj would break it). Every branch's real loss is also ≥ 0.
        let net = case14();
        let view = IndexedNetwork::new(&net);
        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();
        let br = &tables.branch;
        let (pf, pt, status) = (
            col_f64(br, "pf"),
            col_f64(br, "pt"),
            col_i64(br, "br_status"),
        );

        let mut loss = 0.0;
        for i in 0..br.num_rows() {
            if status.value(i) == 1 {
                let l = pf.value(i) + pt.value(i);
                assert!(l >= -1e-6, "branch {i} has negative real loss {l}");
                loss += l;
            }
        }
        assert!(loss > 1.0, "case14 has ~13 MW of real loss, got {loss}");

        let gen_p: f64 = net
            .generators()
            .iter()
            .filter(|g| g.in_service)
            .map(|g| g.pg)
            .sum();
        let load_p: f64 = net.loads().iter().map(|l| l.p).sum();
        // Real shunt draw at the stored voltages: gs() is MW at V = 1.
        let shunt_p: f64 = (0..view.n())
            .map(|i| view.gs()[i] * net.buses()[i].vm.powi(2))
            .sum();
        assert!(
            (loss - (gen_p - load_p - shunt_p)).abs() < 0.5,
            "power balance off: loss {loss} vs gen-load-shunt {}",
            gen_p - load_p - shunt_p
        );
    }

    #[test]
    fn subthreshold_impedance_branch_zeros_columns_and_is_counted() {
        let tiny = powerio_tx::dc::MIN_DIVISIBLE_MAGNITUDE / 2.0;
        let net = BalancedNetwork::in_memory(
            "zeroimp",
            100.0,
            vec![
                bus(1, BusType::Ref),
                bus(2, BusType::Pq),
                bus(3, BusType::Pq),
            ],
            vec![branch(1, 2, tiny, 0.0), branch(2, 3, 0.01, 0.1)],
        );
        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();
        let br = &tables.branch;
        for col in [
            "Yff_r", "Yff_i", "Yft_r", "Yft_i", "Ytf_r", "Ytf_i", "Ytt_r", "Ytt_i", "pf", "qf",
            "pt", "qt",
        ] {
            let v = col_f64(br, col).value(0);
            assert!(
                v == 0.0,
                "{col} should be 0 for the subthreshold-impedance branch, got {v}"
            );
        }

        let dir = tempfile::tempdir().unwrap();
        let out = emit_gridfm_dataset(&net, 0, dir.path(), &GridfmOptions::default()).unwrap();
        assert_eq!(out.dropped_zero_impedance, 1);
        let meta: serde_json::Value = serde_json::from_str(
            &std::fs::read_to_string(out.dir.join("gridfm_meta.json")).unwrap(),
        )
        .unwrap();
        assert_eq!(meta["dropped_zero_impedance"], 1);
    }

    #[test]
    fn gridfm_cost_maps_every_arm_to_raw_coefficients() {
        // Polynomial coeffs are highest-order first: [c2, c1, c0] -> (cp0, cp1, cp2).
        assert_eq!(
            gridfm_cost(Some(&gencost(2, 3, vec![2.0, 3.0, 4.0]))),
            (4.0, 3.0, 2.0)
        );
        assert_eq!(
            gridfm_cost(Some(&gencost(2, 2, vec![3.0, 4.0]))),
            (4.0, 3.0, 0.0)
        );
        assert_eq!(
            gridfm_cost(Some(&gencost(2, 1, vec![4.0]))),
            (4.0, 0.0, 0.0)
        );
        // Unrepresentable rows collapse to zeros and report as not representable.
        let piecewise = gencost(1, 2, vec![0.0, 0.0, 1.0, 1.0]);
        let malformed = gencost(2, 3, vec![1.0]); // fewer coeffs than ncost claims
        assert_eq!(gridfm_cost(Some(&piecewise)), (0.0, 0.0, 0.0));
        assert_eq!(gridfm_cost(Some(&malformed)), (0.0, 0.0, 0.0));
        assert_eq!(gridfm_cost(None), (0.0, 0.0, 0.0));
        assert!(!cost_representable(Some(&piecewise)));
        assert!(!cost_representable(Some(&malformed)));
        assert!(!cost_representable(None));
        assert!(cost_representable(Some(&gencost(
            2,
            3,
            vec![1.0, 2.0, 3.0]
        ))));
    }

    #[test]
    fn missing_reference_bus_errors() {
        // `to_gridfm_record_batches_single` requires exactly one reference bus.
        let net = BalancedNetwork::in_memory(
            "noref",
            100.0,
            vec![bus(1, BusType::Pq), bus(2, BusType::Pq)],
            vec![branch(1, 2, 0.01, 0.1)],
        );
        let err = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(
                err,
                Error::Transmission(powerio_tx::Error::ReferenceBusCount { .. })
            ),
            "got {err:?}"
        );
    }

    #[test]
    fn non_finite_bus_voltage_errors_before_parquet() {
        let mut net = case14();
        net.buses_mut()[0].vm = f64::NAN;
        let err = to_gridfm_record_batches_single(&net, 7, &GridfmOptions::default()).unwrap_err();
        match err {
            Error::NonFiniteGridfmValue {
                scenario,
                element,
                row,
                field,
                value,
            } => {
                assert_eq!(scenario, 7);
                assert_eq!(element, "bus");
                assert_eq!(row, 0);
                assert_eq!(field, "vm");
                assert!(value.is_nan());
            }
            other => panic!("expected NonFiniteGridfmValue, got {other:?}"),
        }
    }

    #[test]
    fn non_finite_tap_errors_even_without_y_bus_table() {
        let mut net = case14();
        net.branches_mut()[0].tap = f64::NAN;
        let opts = GridfmOptions {
            include_y_bus: false,
            ..Default::default()
        };
        let err = to_gridfm_record_batches_single(&net, 0, &opts).unwrap_err();
        assert!(
            matches!(
                err,
                Error::NonFiniteGridfmValue {
                    element: "branch",
                    row: 0,
                    field: "tap",
                    ..
                }
            ),
            "got {err:?}"
        );
    }

    #[test]
    fn non_finite_stated_branch_flow_errors_before_parquet() {
        let mut net = case14();
        net.branches_mut()[0].solution = Some(BranchSolution::new(f64::NAN, 0.0, 0.0, 0.0));
        let err = to_gridfm_record_batches_single(&net, 4, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(
                err,
                Error::NonFiniteGridfmValue {
                    scenario: 4,
                    element: "branch",
                    row: 0,
                    field: "pf",
                    ..
                }
            ),
            "got {err:?}"
        );
    }

    #[test]
    fn writer_reports_only_detected_gridfm_projection_changes() {
        let net = case14();
        let dataset = build_gridfm_dataset(&net, 0, &GridfmOptions::default()).unwrap();
        let lines = powerio_core::render_diagnostics(&dataset.diagnostics);
        assert!(
            lines.iter().any(|line| line.contains("area, zone")),
            "missing bus metadata diagnostic: {lines:?}"
        );
        assert!(
            lines.iter().any(|line| line.contains("no stated solution")),
            "missing derived branch flow diagnostic: {lines:?}"
        );
        assert!(
            lines.iter().all(|line| line.contains("EMIT.GRIDFM.")),
            "unexpected diagnostic family: {lines:?}"
        );
    }

    #[test]
    fn normalized_snapshot_is_rejected_in_release_builds() {
        let net = case14().to_normalized().unwrap();
        let err = to_gridfm_record_batches_single(&net, 3, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(err, Error::NormalizedGridfmSnapshot { scenario: 3 }),
            "got {err:?}"
        );
    }

    #[test]
    fn non_finite_representable_cost_errors() {
        let mut net = BalancedNetwork::in_memory(
            "badcost",
            100.0,
            vec![bus(1, BusType::Ref), bus(2, BusType::Pq)],
            vec![branch(1, 2, 0.01, 0.1)],
        );
        net.generators_mut()
            .push(gen_at(1, gencost(2, 3, vec![f64::NAN, 1.0, 0.0])));

        // coeffs are highest-order first, so the NaN lands in the cp2 column.
        let err = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(
                err,
                Error::NonFiniteGridfmValue {
                    element: "gencost",
                    row: 0,
                    field: "cp2",
                    ..
                }
            ),
            "got {err:?}"
        );
    }

    #[test]
    fn unbounded_limits_export_as_infinity() {
        // PowerModels defaults absent gen limits to ±Inf and MATPOWER carries
        // literal Inf tokens; the unbounded sentinel must reach the Parquet
        // columns, not abort the export (only NaN is corrupt in a limit).
        let mut net = case14();
        net.generators_mut()[0].qmax = f64::INFINITY;
        net.generators_mut()[0].qmin = f64::NEG_INFINITY;
        net.generators_mut()[1].pmax = f64::INFINITY;
        net.branches_mut()[0].angmin = f64::NEG_INFINITY;
        net.branches_mut()[0].angmax = f64::INFINITY;
        net.branches_mut()[1].rate_a = f64::INFINITY;
        net.buses_mut()[0].vmax = f64::INFINITY;
        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();
        let qmax = col_f64(&tables.generator, "max_q_mvar");
        assert!(qmax.value(0).is_infinite() && qmax.value(0) > 0.0);
    }

    #[test]
    fn nan_limit_still_errors() {
        let mut net = case14();
        net.generators_mut()[0].qmax = f64::NAN;
        let err = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(
                err,
                Error::NonFiniteGridfmValue {
                    element: "generator",
                    field: "qmax",
                    ..
                }
            ),
            "got {err:?}"
        );
    }

    /// case14 with every load and generator setpoint scaled — a perturbed
    /// operating point on the same topology, the scenario-batch invariant.
    fn scaled(net: &BalancedNetwork, factor: f64) -> BalancedNetwork {
        let mut s = net.clone();
        for l in s.loads_mut() {
            l.p *= factor;
            l.q *= factor;
        }
        for g in s.generators_mut() {
            g.pg *= factor;
            g.qg *= factor;
        }
        s
    }

    #[test]
    fn batch_stacks_scenarios_keyed_by_scenario_column() {
        let base = case14();
        let up = scaled(&base, 1.1);
        let down = scaled(&base, 0.9);
        let snaps = [
            GridfmSnapshot {
                net: &base,
                scenario: 0,
            },
            GridfmSnapshot {
                net: &up,
                scenario: 1,
            },
            GridfmSnapshot {
                net: &down,
                scenario: 2,
            },
        ];
        let tables = to_gridfm_record_batches(&snaps, &GridfmOptions::default()).unwrap();

        // Schema is unchanged; rows are 3× the single-snapshot counts.
        assert_eq!(names(&tables.bus), BUS_COLS);
        assert_eq!(names(&tables.branch), BRANCH_COLS);
        assert_eq!(tables.bus.num_rows(), 3 * base.buses().len());
        assert_eq!(tables.generator.num_rows(), 3 * base.generators().len());
        assert_eq!(tables.branch.num_rows(), 3 * base.branches().len());

        // The scenario column is blocked 0.., and the dense bus index resets to
        // 0..n within each scenario.
        let n = base.buses().len();
        let scen = col_i64(&tables.bus, "scenario");
        let lsi = col_i64(&tables.bus, "load_scenario_idx");
        let bus_idx = col_i64(&tables.bus, "bus");
        for k in 0..3 {
            for i in 0..n {
                let row = k * n + i;
                assert_eq!(scen.value(row), k as i64);
                assert_eq!(lsi.value(row), k as i64);
                assert_eq!(bus_idx.value(row), i as i64);
            }
        }

        // The first scenario's rows match the standalone single-case tables, so
        // batching is a pure row-stack over the established single-snapshot path.
        // Compare every column bit-exactly (not just one), so a per-column offset
        // or ordering regression in the row-stack can't slip through.
        let single = to_gridfm_record_batches_single(&base, 0, &GridfmOptions::default()).unwrap();
        let bit_exact = |b: &RecordBatch, s: &RecordBatch, col: &str, rows: usize| {
            let (bb, ss) = (col_f64(b, col), col_f64(s, col));
            for i in 0..rows {
                assert_eq!(
                    bb.value(i).to_bits(),
                    ss.value(i).to_bits(),
                    "scenario-0 {col}[{i}] differs from the single-case path"
                );
            }
        };
        for col in ["Pd", "Qd", "Pg", "Qg", "Vm", "Va", "GS", "BS"] {
            bit_exact(&tables.bus, &single.bus, col, n);
        }
        bit_exact(
            &tables.generator,
            &single.generator,
            "p_mw",
            base.generators().len(),
        );
        bit_exact(&tables.branch, &single.branch, "pf", base.branches().len());

        // The perturbed scenario's load really differs (guards against stamping
        // the same network three times).
        let pd_batch = col_f64(&tables.bus, "Pd");
        let pd_single = col_f64(&single.bus, "Pd");
        assert!((pd_batch.value(n) - 1.1 * pd_single.value(0)).abs() < 1e-9);
    }

    #[test]
    fn hostile_network_name_stays_inside_out_dir() {
        // The dataset directory comes from the network name, which is input
        // content; a separator-bearing name must not write outside out_dir.
        let mut net = case14();
        *net.name_mut() = "../escape".to_string();
        let dir = tempfile::tempdir().unwrap();
        let out = emit_gridfm_dataset(&net, 0, dir.path(), &GridfmOptions::default()).unwrap();
        assert!(
            out.dir.starts_with(dir.path()),
            "dataset dir {:?} escaped {:?}",
            out.dir,
            dir.path()
        );
        assert!(out.dir.exists());
    }

    #[test]
    fn batch_dataset_writes_stacked_parquet_with_scenario_count() {
        let base = case14();
        let up = scaled(&base, 1.25);
        let snaps = [
            GridfmSnapshot {
                net: &base,
                scenario: 0,
            },
            GridfmSnapshot {
                net: &up,
                scenario: 1,
            },
        ];
        let dir = tempfile::tempdir().unwrap();
        let out = emit_gridfm_batch(&snaps, dir.path(), &GridfmOptions::default()).unwrap();

        let bus = read(&out.dir.join("bus_data.parquet"));
        assert_eq!(bus.num_rows(), 2 * base.buses().len());
        let scen = col_i64(&bus, "scenario");
        assert_eq!(scen.value(0), 0);
        assert_eq!(scen.value(base.buses().len()), 1);

        let meta: serde_json::Value = serde_json::from_str(
            &std::fs::read_to_string(out.dir.join("gridfm_meta.json")).unwrap(),
        )
        .unwrap();
        assert_eq!(meta["n_scenarios"], 2);
        assert_eq!(meta["scenario"], 0);
    }

    #[test]
    fn empty_batch_errors() {
        let err = to_gridfm_record_batches(&[], &GridfmOptions::default()).unwrap_err();
        assert!(matches!(err, Error::EmptyScenarioBatch), "got {err:?}");
    }

    #[test]
    fn shape_mismatch_across_snapshots_errors() {
        let big = case14();
        let small = BalancedNetwork::in_memory(
            "small",
            100.0,
            vec![bus(1, BusType::Ref), bus(2, BusType::Pq)],
            vec![branch(1, 2, 0.01, 0.1)],
        );
        let snaps = [
            GridfmSnapshot {
                net: &big,
                scenario: 0,
            },
            GridfmSnapshot {
                net: &small,
                scenario: 1,
            },
        ];
        let err = to_gridfm_record_batches(&snaps, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(
                err,
                Error::ScenarioShapeMismatch {
                    index: 1,
                    reason: ScenarioMismatch::Counts { .. }
                }
            ),
            "got {err:?}"
        );
    }

    #[test]
    fn bus_order_mismatch_is_reported_distinctly() {
        // Same counts and the same bus-id set, but a different ordering: the dense
        // bus index would mean different buses across snapshots, so the batch is
        // rejected with the BusOrder reason (not the same-tuple Counts message).
        let base = case14();
        let mut reordered = base.clone();
        reordered.buses_mut().swap(0, 1);
        let snaps = [
            GridfmSnapshot {
                net: &base,
                scenario: 0,
            },
            GridfmSnapshot {
                net: &reordered,
                scenario: 1,
            },
        ];
        let err = to_gridfm_record_batches(&snaps, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(
                err,
                Error::ScenarioShapeMismatch {
                    index: 1,
                    reason: ScenarioMismatch::BusOrder
                }
            ),
            "got {err:?}"
        );
    }

    #[test]
    fn branch_order_mismatch_is_reported_distinctly() {
        let base = case14();
        let mut reordered = base.clone();
        reordered.branches_mut().swap(0, 1);
        let snaps = [
            GridfmSnapshot::new(&base, 0),
            GridfmSnapshot::new(&reordered, 1),
        ];
        let err = to_gridfm_record_batches(&snaps, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(
                err,
                Error::ScenarioShapeMismatch {
                    index: 1,
                    reason: ScenarioMismatch::BranchOrder
                }
            ),
            "got {err:?}"
        );
    }

    #[test]
    fn generator_order_mismatch_is_reported_distinctly() {
        let base = case14();
        let mut reordered = base.clone();
        reordered.generators_mut().swap(0, 1);
        let snaps = [
            GridfmSnapshot::new(&base, 0),
            GridfmSnapshot::new(&reordered, 1),
        ];
        let err = to_gridfm_record_batches(&snaps, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(
                err,
                Error::ScenarioShapeMismatch {
                    index: 1,
                    reason: ScenarioMismatch::GeneratorOrder
                }
            ),
            "got {err:?}"
        );
    }

    #[test]
    fn differing_system_base_is_rejected() {
        let base = case14();
        let mut changed = base.clone();
        *changed.base_mva_mut() = base.base_mva() * 2.0;
        let snaps = [
            GridfmSnapshot::new(&base, 0),
            GridfmSnapshot::new(&changed, 1),
        ];
        let err = to_gridfm_record_batches(&snaps, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(
                err,
                Error::ScenarioShapeMismatch {
                    index: 1,
                    reason: ScenarioMismatch::BaseMva
                }
            ),
            "got {err:?}"
        );
    }

    #[test]
    fn duplicate_scenario_id_is_rejected() {
        let base = case14();
        let snaps = [GridfmSnapshot::new(&base, 7), GridfmSnapshot::new(&base, 7)];
        let err = to_gridfm_record_batches(&snaps, &GridfmOptions::default()).unwrap_err();
        assert!(
            matches!(
                err,
                Error::ScenarioShapeMismatch {
                    index: 1,
                    reason: ScenarioMismatch::DuplicateScenarioId {
                        scenario: 7,
                        first_index: 0
                    }
                }
            ),
            "got {err:?}"
        );
    }

    #[test]
    fn manifest_counts_sum_over_the_batch() {
        // Two snapshots on the same element set, but only the second has a zeroed
        // branch impedance. The manifest's dropped count describes the whole
        // dataset (a total of 1), not just the first snapshot — branch status and
        // impedance may legitimately differ per scenario.
        let base = case14();
        let mut perturbed = base.clone();
        perturbed.branches_mut()[0].r = 0.0;
        perturbed.branches_mut()[0].x = 0.0;
        let snaps = [
            GridfmSnapshot {
                net: &base,
                scenario: 0,
            },
            GridfmSnapshot {
                net: &perturbed,
                scenario: 1,
            },
        ];
        let dir = tempfile::tempdir().unwrap();
        let out = emit_gridfm_batch(&snaps, dir.path(), &GridfmOptions::default()).unwrap();
        assert_eq!(out.dropped_zero_impedance, 1);
        let meta: serde_json::Value = serde_json::from_str(
            &std::fs::read_to_string(out.dir.join("gridfm_meta.json")).unwrap(),
        )
        .unwrap();
        assert_eq!(meta["dropped_zero_impedance"], 1);
    }

    #[test]
    fn y_bus_table_is_absent_when_disabled() {
        let net = case14();
        let opts = GridfmOptions {
            include_y_bus: false,
            ..Default::default()
        };
        let tables = to_gridfm_record_batches_single(&net, 0, &opts).unwrap();
        assert!(tables.y_bus.is_none(), "y_bus should not be built");

        let dir = tempfile::tempdir().unwrap();
        let out = emit_gridfm_dataset(&net, 0, dir.path(), &opts).unwrap();
        assert!(
            !out.dir.join("y_bus_data.parquet").exists(),
            "y_bus_data.parquet should not be written"
        );
    }

    #[test]
    fn number_snapshots_stamps_base_plus_k_and_checks_overflow() {
        // The shared builder both bindings use: the k-th network is scenario
        // `base + k`, in order.
        let net = case14();
        let nets = [&net, &net, &net];
        let snaps = number_snapshots(&nets, 5).unwrap();
        let released = number_snapshots(&nets, 5).unwrap();
        assert_eq!(
            snaps.iter().map(|snap| snap.scenario).collect::<Vec<_>>(),
            released
                .iter()
                .map(|snap| snap.scenario)
                .collect::<Vec<_>>()
        );
        assert_eq!(snaps.len(), 3);
        assert_eq!(snaps[0].scenario, 5);
        assert_eq!(snaps[1].scenario, 6);
        assert_eq!(snaps[2].scenario, 7);

        // Overflow is checked (not wrapped to a negative id, not a panic) and names
        // the offending index.
        let err = number_snapshots(&[&net, &net], i64::MAX).unwrap_err();
        assert!(
            matches!(err, Error::ScenarioIdOverflow { index: 1, .. }),
            "got {err:?}"
        );
    }

    #[test]
    fn out_of_service_generator_is_listed_but_excluded_from_bus_aggregate() {
        // Two paths react to `g.in_service`: gen_data emits an `in_service` column
        // for every generator (keeping its setpoint), while bus `Pg`/`Qg` aggregate
        // only in-service generation (`view.in_service_gens()`). Exercise the
        // `false` case on both.
        let mut net = BalancedNetwork::in_memory(
            "genoutage",
            100.0,
            vec![bus(1, BusType::Ref), bus(2, BusType::Pq)],
            vec![branch(1, 2, 0.01, 0.1)],
        );
        let mut g_on = gen_at(1, gencost(2, 3, vec![0.0, 1.0, 0.0]));
        g_on.pg = 50.0;
        let mut g_off = gen_at(2, gencost(2, 3, vec![0.0, 1.0, 0.0]));
        g_off.pg = 30.0;
        g_off.in_service = false;
        net.generators_mut().push(g_on);
        net.generators_mut().push(g_off);

        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();

        // gen_data lists both gens in source order, flags the out-of-service one,
        // and keeps its setpoint.
        let g = &tables.generator;
        assert_eq!(g.num_rows(), 2);
        let in_service = col_i64(g, "in_service");
        assert_eq!(in_service.value(0), 1, "in-service gen flagged 1");
        assert_eq!(in_service.value(1), 0, "out-of-service gen flagged 0");
        assert!(
            (col_f64(g, "p_mw").value(1) - 30.0).abs() < 1e-12,
            "gen_data keeps the out-of-service setpoint"
        );

        // bus Pg aggregates only in-service generation: bus 1 (dense 0) gets 50,
        // bus 2 (dense 1) excludes the out-of-service gen's 30.
        let pg = col_f64(&tables.bus, "Pg");
        assert!(
            (pg.value(0) - 50.0).abs() < 1e-12,
            "in-service gen folded into bus Pg"
        );
        assert!(
            pg.value(1) == 0.0,
            "out-of-service gen excluded from bus Pg, got {}",
            pg.value(1)
        );
    }

    #[test]
    fn out_of_service_branch_zeros_flows_but_keeps_admittance() {
        // An out-of-service branch keeps its physical Y** admittances but carries
        // zero flows and `br_status = 0` — the path datakit's topology scenarios
        // exercise. Use non-flat voltages so an *in-service* branch carries real
        // flow, which makes the zero on the tripped branch meaningful (not just an
        // artifact of a flat start).
        let mut net = BalancedNetwork::in_memory(
            "outage",
            100.0,
            vec![
                bus(1, BusType::Ref),
                bus(2, BusType::Pq),
                bus(3, BusType::Pq),
            ],
            vec![branch(1, 2, 0.01, 0.1), branch(2, 3, 0.02, 0.2)],
        );
        net.buses_mut()[1].va = -3.0;
        net.buses_mut()[2].va = -6.0;
        net.branches_mut()[0].in_service = false; // trip branch 0

        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();
        let br = &tables.branch;
        let status = col_i64(br, "br_status");
        assert_eq!(status.value(0), 0, "tripped branch reports br_status 0");
        assert_eq!(status.value(1), 1, "in-service branch reports br_status 1");

        for col in ["pf", "qf", "pt", "qt"] {
            let v = col_f64(br, col).value(0);
            assert!(
                v == 0.0,
                "{col} must be zero on the out-of-service branch, got {v}"
            );
        }
        // The in-service branch really carries flow at these voltages — guards
        // against a flat start false pass that would zero every branch anyway.
        assert!(
            col_f64(br, "pf").value(1).abs() > 1e-6,
            "in-service branch should carry nonzero flow"
        );
        // Admittances are retained for the tripped branch (unlike a zero-impedance
        // branch, which zeroes them).
        assert!(
            col_f64(br, "Yff_i").value(0).abs() > 0.0,
            "out-of-service branch keeps its physical Y** admittances"
        );
    }

    fn bus(id: usize, kind: BusType) -> Bus {
        Bus::new(BusId(id), kind, 1.0)
    }

    fn branch(from: usize, to: usize, r: f64, x: f64) -> Branch {
        Branch::new(BusId(from), BusId(to), r, x)
    }

    fn gencost(model: u8, ncost: usize, coeffs: Vec<f64>) -> GenCost {
        GenCost::with_ncost(model, 0.0, 0.0, ncost, coeffs)
    }

    fn gen_at(bus: usize, cost: GenCost) -> Generator {
        let mut generator = Generator::new(BusId(bus));
        generator.pmax = 100.0;
        generator.qmax = 50.0;
        generator.qmin = -50.0;
        generator.mbase = 100.0;
        generator.cost = Some(cost);
        generator
    }

    #[test]
    fn degenerate_cost_gen_zeros_columns_and_is_counted() {
        // Counterpart to the zero-impedance test: a piecewise (model 1) cost row
        // gets zeroed cp* columns and is counted; a polynomial row is kept.
        let mut net = BalancedNetwork::in_memory(
            "degen",
            100.0,
            vec![bus(1, BusType::Ref), bus(2, BusType::Pq)],
            vec![branch(1, 2, 0.01, 0.1)],
        );
        net.generators_mut()
            .push(gen_at(1, gencost(1, 2, vec![0.0, 0.0, 1.0, 1.0]))); // piecewise
        net.generators_mut()
            .push(gen_at(2, gencost(2, 3, vec![0.01, 5.0, 0.0]))); // polynomial

        let tables = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default()).unwrap();
        let g = &tables.generator;
        let (cp0, cp1, cp2) = (
            col_f64(g, "cp0_eur"),
            col_f64(g, "cp1_eur_per_mw"),
            col_f64(g, "cp2_eur_per_mw2"),
        );
        assert_eq!((cp0.value(0), cp1.value(0), cp2.value(0)), (0.0, 0.0, 0.0));
        assert_eq!((cp0.value(1), cp1.value(1), cp2.value(1)), (0.0, 5.0, 0.01));

        let dir = tempfile::tempdir().unwrap();
        let out = emit_gridfm_dataset(&net, 0, dir.path(), &GridfmOptions::default()).unwrap();
        assert_eq!(out.degenerate_cost_gens, 1);
        assert_eq!(out.missing_cost_gens, 0);
        assert_eq!(out.unsupported_cost_gens, 1);
        let meta: serde_json::Value = serde_json::from_str(
            &std::fs::read_to_string(out.dir.join("gridfm_meta.json")).unwrap(),
        )
        .unwrap();
        assert_eq!(meta["degenerate_cost_gens"], 1);
        assert_eq!(meta["missing_cost_gens"], 0);
        assert_eq!(meta["unsupported_cost_gens"], 1);
        assert_eq!(meta["zeroed_cost_gens"], 1);
    }

    #[test]
    fn missing_gridfm_costs_are_split_and_fill_policy_reduces_missing_count() {
        let mut net = BalancedNetwork::in_memory(
            "missing-cost",
            100.0,
            vec![bus(1, BusType::Ref), bus(2, BusType::Pq)],
            vec![branch(1, 2, 0.01, 0.1)],
        );
        let mut missing = gen_at(1, gencost(2, 3, vec![1.0, 2.0, 3.0]));
        missing.cost = None;
        net.generators_mut().push(missing);
        net.generators_mut()
            .push(gen_at(2, gencost(1, 2, vec![0.0, 0.0, 1.0, 1.0])));

        let dir = tempfile::tempdir().unwrap();
        let out = emit_gridfm_dataset(&net, 0, dir.path(), &GridfmOptions::default()).unwrap();
        assert_eq!(out.degenerate_cost_gens, 2);
        assert_eq!(out.missing_cost_gens, 1);
        assert_eq!(out.unsupported_cost_gens, 1);

        let opts = GridfmOptions {
            missing_gen_cost: MissingGenCostPolicy::zero(),
            ..Default::default()
        };
        let dir = tempfile::tempdir().unwrap();
        let out = emit_gridfm_dataset(&net, 0, dir.path(), &opts).unwrap();
        assert_eq!(out.synthesized_gen_costs, 1);
        assert_eq!(out.missing_cost_gens, 0);
        assert_eq!(out.unsupported_cost_gens, 1);
        assert_eq!(out.degenerate_cost_gens, 1);
        let meta: serde_json::Value = serde_json::from_str(
            &std::fs::read_to_string(out.dir.join("gridfm_meta.json")).unwrap(),
        )
        .unwrap();
        assert_eq!(meta["cost_policy"]["mode"], "fill");
        assert_eq!(meta["synthesized_gen_costs"], 1);
    }

    #[test]
    fn scenario_id_and_tap_toggle_take_effect() {
        let net = case14();

        // The scenario id (an explicit argument now) reaches both id columns.
        let bus = to_gridfm_record_batches_single(&net, 7, &GridfmOptions::default())
            .unwrap()
            .bus;
        assert_eq!(col_i64(&bus, "scenario").value(0), 7);
        assert_eq!(col_i64(&bus, "load_scenario_idx").value(0), 7);

        // Turning taps off changes a transformer's admittance columns.
        let on = to_gridfm_record_batches_single(&net, 0, &GridfmOptions::default())
            .unwrap()
            .branch;
        let off = to_gridfm_record_batches_single(
            &net,
            0,
            &GridfmOptions {
                include_taps: false,
                ..Default::default()
            },
        )
        .unwrap()
        .branch;
        let tap = col_f64(&on, "tap");
        let xfmr = (0..on.num_rows())
            .find(|&i| (tap.value(i) - 1.0).abs() > 1e-9)
            .expect("case14 has off-nominal transformers");
        // case14 transformers are lossless (r = 0), so compare the susceptance
        // (imaginary) part, which scales with 1/tap².
        assert!(
            (col_f64(&on, "Yff_i").value(xfmr) - col_f64(&off, "Yff_i").value(xfmr)).abs() > 1e-9,
            "taps off should change the transformer's Yff"
        );
    }

    // --- reader (issue #60) ------------------------------------------------
}