oximo-baron 0.7.0

use oximo_solver::prepare::{LoweringContext, PreparedExpressions};
use oximo_solver::reconstruct::normalize_result;
use std::borrow::Cow;
use std::fmt::Write as FmtWrite;
use std::process::Stdio;
use std::sync::atomic::{AtomicU64, Ordering};
use std::time::Instant;
use std::{fs, io};

use oximo_core::{
    Constraint, ConstraintId, Domain, Model, ModelId, ModelMismatchError, Objective,
    ObjectiveSense, Sense, SocConstraint, SocConstraintId, VarId, Variable,
};
use oximo_expr::{ExprArena, ExprId, ExprNode, LinearTerms, UnaryOp};
use oximo_solver::{
    DualStatus, Iis, PrimalStatus, SolutionPoint, SolverError, SolverResult, TerminationStatus,
    VarBoundKind,
};
use rustc_hash::{FxHashMap, FxHashSet};

use crate::BaronOptions;
use crate::options::write_options;

static SOLVE_ID: AtomicU64 = AtomicU64::new(0);

const RES_NAME: &str = "res.lst";
const TIM_NAME: &str = "tim.lst";
const BAR_NAME: &str = "problem.bar";

/// Write `model` to a temporary BARON `.bar` file, execute the `baron`
/// executable, and return the parsed [`SolverResult`].
///
/// `exec` is an optional override for the BARON executable path. `None` uses
/// `"baron"` resolved from `PATH`. [`BaronOptions::baron_path`] takes precedence
/// over `exec`.
///
/// # Errors
///
/// Returns [`SolverError`] on constructs BARON's `.bar` format cannot represent
/// (`sin`/`cos`, semicontinuous/semi-integer variables), a missing BARON
/// executable, a BARON run that produced no times file, or I/O failures.
///
/// # Panics
///
/// Panics if variable indices overflow `u32`.
pub fn solve(
    model: &Model,
    opts: &BaronOptions,
    exec: Option<&str>,
) -> Result<SolverResult, SolverError> {
    if model.has_active_sos_constraints() {
        return Err(SolverError::UnsupportedSos);
    }
    if model.has_active_indicator_constraints() {
        return Err(SolverError::UnsupportedIndicator);
    }
    let sense = model.objective().as_ref().map_or(ObjectiveSense::Minimize, |o| o.sense);
    let (bar, var_order, con_order, soc_bounds) = build_bar(model, opts)?;
    let run = run_baron(&bar, opts, exec)?;
    Ok(parse_solution(
        &run.tim,
        &run.res,
        sense,
        run.elapsed,
        run.raw_log,
        &var_order,
        &con_order,
        &soc_bounds,
        model.id(),
    ))
}

/// Build `model`, run BARON with IIS search enabled, and when BARON reports the
/// model infeasible, parse the irreducible infeasible subsystem it prints in the
/// results file.
///
/// If the caller has not set `CompIIS`, it is enabled (`CompIIS: 1`). A user-chosen
/// `CompIIS` algorithm is left untouched.
///
/// # Errors
///
/// Returns a [`SolverError`] on the same conditions as [`solve`], or when the model is
/// not infeasible (there is nothing to diagnose).
///
/// # Panics
///
/// Panics if variable indices overflow `u32`.
pub fn compute_iis(
    model: &Model,
    opts: &BaronOptions,
    exec: Option<&str>,
) -> Result<Iis, SolverError> {
    let iis_opts = if opts.has_comp_iis() { opts.clone() } else { opts.clone().comp_iis(1) };
    let (bar, _var_order, _con_order, _soc_bounds) = build_bar(model, &iis_opts)?;
    let run = run_baron(&bar, &iis_opts, exec)?;

    let termination = map_status(&run.res, tim_model_status(&run.tim));
    if !termination.is_infeasible() {
        return Err(SolverError::Backend(format!(
            "cannot compute an IIS: model is not infeasible ({termination:?})"
        )));
    }

    let iis = parse_iis(&run.res);
    if iis.is_empty() {
        return Err(SolverError::Backend(
            "BARON reported the model infeasible but emitted no IIS block \
             (CompIIS may have failed or been unable to isolate a subsystem)"
                .to_string(),
        ));
    }
    Ok(iis)
}

/// The text files BARON leaves in its working directory after a run.
struct BaronRun {
    /// Contents of the times file (`tim.lst`): the single status line.
    tim: String,
    /// Contents of the results file (`res.lst`): the primal/dual blocks and, when
    /// `CompIIS` is on, the IIS listing.
    res: String,
    elapsed: std::time::Duration,
    /// Captured stdout/stderr, surfaced on failure (`None` when the run succeeded or
    /// output was streamed in verbose mode).
    raw_log: Option<String>,
}

/// Write `bar` to a fresh temp directory, execute BARON, read back its times and
/// results files, then remove the directory. Shared by [`solve`] and [`compute_iis`].
///
/// # Errors
///
/// Returns a [`SolverError`] if the temp directory or `.bar` file cannot be written,
/// the `baron` executable is missing, or BARON produced no times file (a syntax or
/// license error that never reached the solve).
fn run_baron(bar: &str, opts: &BaronOptions, exec: Option<&str>) -> Result<BaronRun, SolverError> {
    // - Temp directory. Combine a timestamp with a per-process atomic counter so
    //   concurrent invocations never share a directory.
    let ts = std::time::SystemTime::now()
        .duration_since(std::time::UNIX_EPOCH)
        .map_or(0, |d| d.as_millis());
    let id = SOLVE_ID.fetch_add(1, Ordering::Relaxed);
    let tmp_dir = std::env::temp_dir().join(format!("oximo_baron_{ts}_{id}"));
    fs::create_dir_all(&tmp_dir)
        .map_err(|e| SolverError::Backend(format!("cannot create temp dir: {e}")))?;

    let bar_path = tmp_dir.join(BAR_NAME);
    fs::write(&bar_path, bar)
        .map_err(|e| SolverError::Backend(format!("cannot write .bar file: {e}")))?;

    // - Execute BARON.
    let baron_exec =
        opts.baron_path.as_deref().and_then(std::path::Path::to_str).or(exec).unwrap_or("baron");
    let verbose = opts.universal.verbose.unwrap_or(false);

    let started = Instant::now();
    let mut cmd = std::process::Command::new(baron_exec);
    cmd.arg(BAR_NAME);
    cmd.current_dir(&tmp_dir);

    let launch_err = |e: io::Error| {
        let _ = fs::remove_dir_all(&tmp_dir);
        if e.kind() == io::ErrorKind::NotFound {
            SolverError::Backend(format!(
                "BARON executable '{baron_exec}' not found. \
                Install BARON and ensure it is on PATH, or set the 'baron_path' option."
            ))
        } else {
            SolverError::Backend(format!("failed to launch BARON: {e}"))
        }
    };

    // When verbose, stream BARON's output to the terminal. Otherwise capture it
    // so we can surface it on failure.
    let (exit_ok, raw_log) = if verbose {
        let status =
            cmd.stdout(Stdio::inherit()).stderr(Stdio::inherit()).status().map_err(launch_err)?;
        (status.success(), None)
    } else {
        let out = cmd.output().map_err(launch_err)?;
        let log = if out.status.success() {
            None
        } else {
            let mut s = String::from_utf8_lossy(&out.stdout).into_owned();
            if !out.stderr.is_empty() {
                s.push('\n');
                s.push_str(&String::from_utf8_lossy(&out.stderr));
            }
            Some(s)
        };
        (out.status.success(), log)
    };
    let elapsed = started.elapsed();

    // - Parse the times file (status) and results file (primal).
    let tim_path = tmp_dir.join(TIM_NAME);
    if !tim_path.exists() {
        // BARON never reached the solve (syntax/license error). Surface its log.
        let _ = fs::remove_dir_all(&tmp_dir);
        let detail = raw_log.unwrap_or_else(|| {
            if exit_ok {
                "BARON produced no times file and emitted no output.".to_string()
            } else {
                "BARON exited with a non-zero exit code and produced no times file.".to_string()
            }
        });
        return Err(SolverError::Backend(format!("BARON did not produce a solution.\n{detail}")));
    }

    let tim = fs::read_to_string(&tim_path)
        .map_err(|e| SolverError::Backend(format!("cannot read times file: {e}")))?;
    let res = fs::read_to_string(tmp_dir.join(RES_NAME)).unwrap_or_default();

    let _ = fs::remove_dir_all(&tmp_dir);
    Ok(BaronRun { tim, res, elapsed, raw_log })
}

/// Everything `solve` needs from the `.bar` writer: the file text, the
/// variable declaration order (see [`write_var_declarations`]) used to decode
/// BARON's numeric solution-pool blocks, the constraint emit-order (see
/// [`write_equations`]) used to fold range duals back onto their
/// `ConstraintId`, and the affine bound side of each explicit SOC constraint
/// used to rescale its squared-row price to the norm form.
type BarParts = (String, Vec<VarId>, Vec<ConstraintId>, Vec<LinearTerms<'static>>);

fn build_bar(model: &Model, opts: &BaronOptions) -> Result<BarParts, SolverError> {
    if model.has_active_sos_constraints() {
        return Err(SolverError::UnsupportedSos);
    }
    if model.has_active_indicator_constraints() {
        return Err(SolverError::UnsupportedIndicator);
    }
    let prepared = LoweringContext::new(model)?;
    let vars = prepared.variables();
    let model_constraints = prepared.constraints();
    let constraints = model_constraints.algebraic();
    let socs = prepared.constraints().second_order_cones();
    let objective = prepared.objective();

    let mut bar = String::with_capacity(4096);
    write_options(&mut bar, opts, RES_NAME, TIM_NAME);
    let var_order = write_var_declarations(&mut bar, vars)?;
    write_bounds(&mut bar, vars);
    let convex_equation_ids = opts
        .convex_equation_handles()
        .iter()
        .map(|handle| {
            if handle.model_id() == model.id() {
                Ok(handle.id())
            } else {
                let mismatch = ModelMismatchError::new(model.id(), handle.model_id());
                Err(SolverError::Backend(format!(
                    "BARON convex-equation assertion has a {mismatch}"
                )))
            }
        })
        .collect::<Result<Vec<_>, _>>()?;
    let con_order = write_equations(&mut bar, &prepared, constraints, socs, &convex_equation_ids)?;
    write_objective(&mut bar, &prepared, objective)?;
    write_starting_point(&mut bar, vars);
    let soc_bounds = socs
        .iter()
        .map(|s| prepared.linear(s.bound).expect("SOC bound is validated affine").into_owned())
        .collect();
    Ok((bar, var_order, con_order, soc_bounds))
}

// - .bar writers

/// Emit the `BINARY_VARIABLES` / `INTEGER_VARIABLES` / `POSITIVE_VARIABLES` /
/// `VARIABLES` declaration sections.
///
/// Returns the variable `VarId`s in the order they are declared (binary,
/// integer, positive, free). BARON numbers variables 1-based in this order, so
/// the result maps "Variable no. `k`" back to a `VarId` when parsing the
/// enumerated solution pool.
fn write_var_declarations(bar: &mut String, vars: &[Variable]) -> Result<Vec<VarId>, SolverError> {
    let (mut bin, mut int, mut pos, mut free) = (Vec::new(), Vec::new(), Vec::new(), Vec::new());
    for v in vars {
        match v.domain {
            Domain::Binary => bin.push(v),
            Domain::Integer => int.push(v),
            // A continuous variable with a zero lower bound is exactly BARON's
            // POSITIVE_VARIABLES; everything else (free, negative, finite lb) is
            // a general VARIABLES with explicit bounds emitted below.
            Domain::Real if v.lb == 0.0 => pos.push(v),
            Domain::Real => free.push(v),
            Domain::SemiContinuous { .. } | Domain::SemiInteger { .. } => {
                return Err(SolverError::Backend(format!(
                    "variable x{} has a semicontinuous/semi-integer domain, \
                    which BARON's .bar format cannot represent",
                    v.id.index()
                )));
            }
        }
    }
    write_var_section(bar, "BINARY_VARIABLES", &bin);
    write_var_section(bar, "INTEGER_VARIABLES", &int);
    write_var_section(bar, "POSITIVE_VARIABLES", &pos);
    write_var_section(bar, "VARIABLES", &free);
    writeln!(bar).unwrap();
    let order = bin.iter().chain(&int).chain(&pos).chain(&free).map(|v| v.id).collect();
    Ok(order)
}

fn write_var_section(bar: &mut String, header: &str, vars: &[&Variable]) {
    if vars.is_empty() {
        return;
    }
    write!(bar, "{header} ").unwrap();
    for (k, v) in vars.iter().enumerate() {
        if k > 0 {
            write!(bar, ", ").unwrap();
        }
        write!(bar, "x{}", v.id.index()).unwrap();
    }
    writeln!(bar, ";").unwrap();
}

fn write_bounds(bar: &mut String, vars: &[Variable]) {
    let mut lo = String::new();
    let mut hi = String::new();
    for v in vars {
        let i = v.id.index();
        if let Some(lb) = lower_bound_to_emit(v) {
            writeln!(lo, "x{i}: {};", fmt(lb)).unwrap();
        }
        if let Some(ub) = upper_bound_to_emit(v) {
            writeln!(hi, "x{i}: {};", fmt(ub)).unwrap();
        }
    }
    if !lo.is_empty() {
        writeln!(bar, "LOWER_BOUNDS{{").unwrap();
        bar.push_str(&lo);
        writeln!(bar, "}}").unwrap();
        writeln!(bar).unwrap();
    }
    if !hi.is_empty() {
        writeln!(bar, "UPPER_BOUNDS{{").unwrap();
        bar.push_str(&hi);
        writeln!(bar, "}}").unwrap();
        writeln!(bar).unwrap();
    }
}

/// The lower bound to write, or `None` when it equals the implied default for
/// the variable's declaration section (so we avoid redundant lines).
fn lower_bound_to_emit(v: &Variable) -> Option<f64> {
    if !v.lb.is_finite() {
        return None;
    }
    match v.domain {
        // Binary defaults to [0, 1]. Real with lb == 0 sits in POSITIVE_VARIABLES.
        Domain::Binary | Domain::Real if v.lb == 0.0 => None,
        _ => Some(v.lb),
    }
}

fn upper_bound_to_emit(v: &Variable) -> Option<f64> {
    if !v.ub.is_finite() {
        return None;
    }
    match v.domain {
        // Binary defaults to an upper bound of 1. Only emit when overridden.
        Domain::Binary if (v.ub - 1.0).abs() < f64::EPSILON => None,
        _ => Some(v.ub),
    }
}

type EquationRow = (&'static str, &'static str, f64);

fn equation_rows(c: &Constraint) -> [Option<EquationRow>; 2] {
    oximo_solver::prepare::row_sides(c).map(|side| {
        side.map(|(sense, rhs)| {
            let suffix =
                if c.is_range() { if sense == Sense::Ge { "_lo" } else { "_hi" } } else { "" };
            let op = match sense {
                Sense::Le => "<=",
                Sense::Ge => ">=",
                Sense::Eq => "==",
            };
            (suffix, op, rhs)
        })
    })
}

/// Write the `EQUATIONS` block, returning the emit-order map.
/// The `ConstraintId` behind each emitted equation, in BARON's 1-based
/// "Constraint no." order. BARON has no two-sided equation,
/// so a range emits two rows (`c{i}_lo`/`c{i}_hi`) that both map back
/// to the one `ConstraintId`. The dual parser folds them together so
/// positional indices stay correct.
fn write_equations(
    bar: &mut String,
    prepared: &PreparedExpressions,
    constraints: &[Constraint],
    socs: &[SocConstraint],
    convex_equation_ids: &[ConstraintId],
) -> Result<Vec<ConstraintId>, SolverError> {
    let arena = prepared.arena();
    let mut emit_map: Vec<ConstraintId> = Vec::with_capacity(constraints.len());
    let mut names: Vec<String> = Vec::with_capacity(constraints.len() + 2 * socs.len());
    let mut bodies = String::new();

    for (i, c) in constraints.iter().enumerate() {
        let rows = equation_rows(c);
        if rows.iter().all(Option::is_none) {
            // Free `[-inf, +inf]` row
            continue;
        }
        // BARON rejects a constraint whose expression evaluates to a constant,
        // so surface a clear error instead of emitting an invalid `.bar`.
        if !expr_has_var(arena, c.lhs) {
            return Err(SolverError::Backend(format!(
                "constraint '{}' has no variables (its left-hand side is constant); \
                BARON requires every constraint to contain at least one variable",
                c.name
            )));
        }
        let id = ConstraintId(u32::try_from(i).expect("constraint count overflow"));
        for (suffix, op, rhs) in rows.into_iter().flatten() {
            names.push(format!("c{i}{suffix}"));
            emit_map.push(id);
            write!(bodies, "c{i}{suffix}: ").unwrap();
            write_constraint_body(&mut bodies, prepared, c.lhs, op, rhs)?;
        }
    }
    for i in 0..socs.len() {
        names.push(format!("soc{i}"));
        names.push(format!("soc{i}_sign"));
    }
    let convex_names = convex_equation_names(constraints, convex_equation_ids)?;
    if names.is_empty() {
        return Ok(emit_map);
    }

    writeln!(bar, "EQUATIONS {};", names.join(", ")).unwrap();
    if !convex_names.is_empty() {
        writeln!(bar, "CONVEX_EQUATIONS {};", convex_names.join(", ")).unwrap();
    }
    bar.push_str(&bodies);
    write_soc_rows(bar, prepared, socs);
    writeln!(bar).unwrap();
    Ok(emit_map)
}

/// Resolve user-supplied algebraic IDs to the exact BARON row names emitted by
/// [`write_equations`]. A range becomes two BARON rows and is rejected: asserting
/// convexity for both sides is not implied by convexity of the combined range.
fn convex_equation_names(
    constraints: &[Constraint],
    requested: &[ConstraintId],
) -> Result<Vec<String>, SolverError> {
    let mut names = Vec::with_capacity(requested.len());
    for &id in requested {
        let Some(constraint) = constraints.get(id.index()) else {
            return Err(SolverError::Backend(format!(
                "BARON convex-equation assertion references unknown algebraic constraint {id:?}"
            )));
        };
        if !constraint.active {
            return Err(SolverError::Backend(format!(
                "BARON convex-equation assertion references inactive constraint '{}' ({id:?})",
                constraint.name
            )));
        }
        if constraint.is_range() {
            return Err(SolverError::Backend(format!(
                "BARON convex-equation assertion cannot mark two-sided range constraint '{}' \
                 ({id:?}); express it as separate single-sided constraints and mark only the \
                 convex sides",
                constraint.name
            )));
        }
        let Some((suffix, _, _)) = equation_rows(constraint).into_iter().flatten().next() else {
            return Err(SolverError::Backend(format!(
                "BARON convex-equation assertion references constraint '{}' ({id:?}), which \
                 emits no equation",
                constraint.name
            )));
        };
        names.push(format!("c{}{suffix}", id.index()));
    }
    Ok(names)
}

/// Emit each explicit SOC constraint `||terms||_2 <= bound` as the polynomial
/// row `(term_1)^2 + ... - (bound)^2 <= 0` plus the sign row `bound >= 0`
/// (squaring loses the sign of the bound side).
fn write_soc_rows(bar: &mut String, prepared: &PreparedExpressions, socs: &[SocConstraint]) {
    for (i, s) in socs.iter().enumerate() {
        write!(bar, "soc{i}: ").unwrap();
        for (k, &term) in s.terms.iter().enumerate() {
            if k > 0 {
                write!(bar, " + ").unwrap();
            }
            let t = prepared.linear(term).expect("SOC members are validated affine");
            write!(bar, "(").unwrap();
            write_linear(bar, &t, true);
            write!(bar, ")^2").unwrap();
        }
        let b = prepared.linear(s.bound).expect("SOC bound is validated affine");
        write!(bar, " - (").unwrap();
        write_linear(bar, &b, true);
        writeln!(bar, ")^2 <= 0;").unwrap();

        write!(bar, "soc{i}_sign: ").unwrap();
        write_linear(bar, &b, true);
        writeln!(bar, " >= 0;").unwrap();
    }
}

/// Emit one constraint body `<lhs> <op> <rhs>;`. Linear constraints fold the
/// constant into the RHS (canonical `lhs <op> rhs`). Nonlinear ones emit the
/// full LHS and keep the original RHS.
fn write_constraint_body(
    bar: &mut String,
    prepared: &PreparedExpressions,
    lhs: ExprId,
    op: &str,
    rhs: f64,
) -> Result<(), SolverError> {
    let arena = prepared.arena();
    if let Some(t) = prepared.linear(lhs) {
        let adjusted_rhs = rhs - t.constant;
        write_linear(bar, &t, false);
        writeln!(bar, " {op} {};", fmt(adjusted_rhs)).unwrap();
    } else {
        write_bar_expr(bar, arena, lhs)?;
        writeln!(bar, " {op} {};", fmt(rhs)).unwrap();
    }
    Ok(())
}

fn write_objective(
    bar: &mut String,
    prepared: &PreparedExpressions,
    objective: Option<&Objective>,
) -> Result<(), SolverError> {
    let arena = prepared.arena();
    write!(bar, "OBJ: ").unwrap();
    match objective {
        // BARON requires an objective. A feasibility problem minimizes a constant.
        None => writeln!(bar, "minimize 0;").unwrap(),
        Some(o) => {
            let kw = match o.sense {
                ObjectiveSense::Minimize => "minimize",
                ObjectiveSense::Maximize => "maximize",
            };
            write!(bar, "{kw} ").unwrap();
            if let Some(t) = prepared.linear(o.expr) {
                write_linear(bar, &t, true);
            } else {
                write_bar_expr(bar, arena, o.expr)?;
            }
            writeln!(bar, ";").unwrap();
        }
    }
    writeln!(bar).unwrap();
    Ok(())
}

fn write_starting_point(bar: &mut String, vars: &[Variable]) {
    if !vars.iter().any(|v| v.initial.is_some()) {
        return;
    }
    writeln!(bar, "STARTING_POINT{{").unwrap();
    for v in vars {
        if let Some(val) = v.initial {
            writeln!(bar, "x{}: {};", v.id.index(), fmt(val)).unwrap();
        }
    }
    writeln!(bar, "}}").unwrap();
    writeln!(bar).unwrap();
}

/// Append the linear expression `t` to `bar` as a standalone BARON expression.
fn write_linear(bar: &mut String, t: &LinearTerms<'_>, include_constant: bool) {
    let mut first = true;
    for &(v, coef) in t.coeffs.iter() {
        if coef == 0.0 {
            continue;
        }
        let idx = v.index();
        if first {
            write!(bar, "{}*x{idx}", fmt(coef)).unwrap();
            first = false;
        } else if coef < 0.0 {
            write!(bar, " - {}*x{idx}", fmt(-coef)).unwrap();
        } else {
            write!(bar, " + {}*x{idx}", fmt(coef)).unwrap();
        }
    }
    if include_constant && t.constant != 0.0 {
        if first {
            write!(bar, "{}", fmt(t.constant)).unwrap();
            first = false;
        } else if t.constant < 0.0 {
            write!(bar, " - {}", fmt(-t.constant)).unwrap();
        } else {
            write!(bar, " + {}", fmt(t.constant)).unwrap();
        }
    }
    if first {
        write!(bar, "0").unwrap();
    }
}

/// Recursive infix printer for a BARON-compatible expression.
#[expect(clippy::too_many_lines)]
fn write_bar_expr(bar: &mut String, arena: &ExprArena, id: ExprId) -> Result<(), SolverError> {
    match arena.get(id) {
        ExprNode::Const(c) => write!(bar, "{}", fmt(*c)).unwrap(),
        ExprNode::Var(v) => write!(bar, "x{}", v.index()).unwrap(),
        ExprNode::Param(p) => write!(bar, "{}", fmt(arena.param_value(*p))).unwrap(),
        ExprNode::Linear { coeffs, constant } => {
            let t = LinearTerms::borrowed(coeffs, *constant);
            write!(bar, "(").unwrap();
            write_linear(bar, &t, true);
            write!(bar, ")").unwrap();
        }
        ExprNode::Unary(UnaryOp::Neg, inner) => {
            write!(bar, "(-").unwrap();
            write_bar_expr(bar, arena, *inner)?;
            write!(bar, ")").unwrap();
        }
        ExprNode::Add(children) => {
            write!(bar, "(").unwrap();
            for (i, c) in children.iter().enumerate() {
                if i > 0 {
                    write!(bar, " + ").unwrap();
                }
                write_bar_expr(bar, arena, *c)?;
            }
            write!(bar, ")").unwrap();
        }
        ExprNode::Mul(children) => {
            write!(bar, "(").unwrap();
            for (i, c) in children.iter().enumerate() {
                if i > 0 {
                    write!(bar, " * ").unwrap();
                }
                write_bar_expr(bar, arena, *c)?;
            }
            write!(bar, ")").unwrap();
        }
        ExprNode::Pow(base, exp) => {
            // BARON natively supports `x^a` for a constant exponent `a` and
            // `b^x` for a constant base `b`. Only a variable-on-variable power
            // needs the `exp(y*log(x))` rewrite.
            let exp_is_const = matches!(arena.get(*exp), ExprNode::Const(_));
            let base_is_const = matches!(arena.get(*base), ExprNode::Const(_));
            if exp_is_const || base_is_const {
                write!(bar, "(").unwrap();
                write_bar_expr(bar, arena, *base)?;
                write!(bar, " ^ ").unwrap();
                write_bar_expr(bar, arena, *exp)?;
                write!(bar, ")").unwrap();
            } else {
                write!(bar, "exp((").unwrap();
                write_bar_expr(bar, arena, *exp)?;
                write!(bar, ") * log(").unwrap();
                write_bar_expr(bar, arena, *base)?;
                write!(bar, "))").unwrap();
            }
        }
        ExprNode::Div(num, den) => {
            write!(bar, "(").unwrap();
            write_bar_expr(bar, arena, *num)?;
            write!(bar, " / ").unwrap();
            write_bar_expr(bar, arena, *den)?;
            write!(bar, ")").unwrap();
        }
        ExprNode::Unary(UnaryOp::Exp, a) => {
            write!(bar, "exp(").unwrap();
            write_bar_expr(bar, arena, *a)?;
            write!(bar, ")").unwrap();
        }
        ExprNode::Unary(UnaryOp::Log, a) => {
            write!(bar, "log(").unwrap();
            write_bar_expr(bar, arena, *a)?;
            write!(bar, ")").unwrap();
        }
        ExprNode::Unary(UnaryOp::Sqrt, a) => {
            write!(bar, "((").unwrap();
            write_bar_expr(bar, arena, *a)?;
            write!(bar, ") ^ 0.5)").unwrap();
        }
        ExprNode::Unary(UnaryOp::Exp2, a) => {
            write!(bar, "(2 ^ (").unwrap();
            write_bar_expr(bar, arena, *a)?;
            write!(bar, "))").unwrap();
        }
        ExprNode::Unary(UnaryOp::Log10, a) => {
            write!(bar, "(0.4342944819032518 * log(").unwrap();
            write_bar_expr(bar, arena, *a)?;
            write!(bar, "))").unwrap();
        }
        ExprNode::Unary(UnaryOp::Abs, a) => {
            // BARON has no abs() intrinsic.
            // We reformulate: |x| = (x^2)^(1/2),
            // As suggested by the BARON user manual.
            write!(bar, "(((").unwrap();
            write_bar_expr(bar, arena, *a)?;
            write!(bar, ") ^ 2) ^ 0.5)").unwrap();
        }
        ExprNode::Unary(op, _) => {
            return Err(SolverError::UnsupportedNonlinearOperator {
                backend: "BARON",
                operator: op.name(),
            });
        }
        ExprNode::Atan2(_, _) => {
            return Err(SolverError::UnsupportedNonlinearOperator {
                backend: "BARON",
                operator: "atan2",
            });
        }
        ExprNode::Min(_) => {
            return Err(SolverError::UnsupportedNonlinearOperator {
                backend: "BARON",
                operator: "min",
            });
        }
        ExprNode::Max(_) => {
            return Err(SolverError::UnsupportedNonlinearOperator {
                backend: "BARON",
                operator: "max",
            });
        }
    }
    Ok(())
}

/// Whether `id` references at least one variable (vs. evaluating to a constant).
/// Used to enforce BARON's rule that every constraint must contain a non-constant
/// expression.
fn expr_has_var(arena: &ExprArena, id: ExprId) -> bool {
    match arena.get(id) {
        ExprNode::Var(_) => true,
        ExprNode::Const(_) | ExprNode::Param(_) => false,
        ExprNode::Linear { coeffs, .. } => coeffs.iter().any(|(_, c)| *c != 0.0),
        ExprNode::Unary(_, a) => expr_has_var(arena, *a),
        ExprNode::Pow(a, b) | ExprNode::Div(a, b) | ExprNode::Atan2(a, b) => {
            expr_has_var(arena, *a) || expr_has_var(arena, *b)
        }
        ExprNode::Add(children)
        | ExprNode::Mul(children)
        | ExprNode::Min(children)
        | ExprNode::Max(children) => children.iter().any(|c| expr_has_var(arena, *c)),
    }
}

/// Format an `f64` for use in a `.bar` file.
fn fmt(v: f64) -> String {
    if v == f64::INFINITY {
        return "1e51".into();
    }
    if v == f64::NEG_INFINITY {
        return "-1e51".into();
    }
    format!("{v}")
}

// - Result parsing

/// Parse BARON's times file (`tim.lst`) and results file (`res.lst`).
///
/// Read BARON's model status from the times file (`tim.lst`). The 9th whitespace
/// token of its single status line. Defaults to `5` (a BARON status oximo maps to an
/// `Other` termination) when the field is absent or unparseable.
fn tim_model_status(tim: &str) -> i64 {
    tim.split_whitespace().nth(8).and_then(|s| s.parse::<i64>().ok()).unwrap_or(5)
}

/// The times file is a single whitespace-separated line. Field positions follow
/// the BARON convention (0-indexed here):
/// `[5]` lower bound, `[6]` upper bound, `[8]` model status, `[10]`
/// branch-and-reduce iterations, `[11]` node where the optimum was found
/// (`-3` => no solution), last = wall time. The termination reason comes from
/// the `res.lst` banner (see [`map_status`]), not the numeric `[7]` solver status.
#[expect(clippy::too_many_arguments)]
fn parse_solution(
    tim: &str,
    res: &str,
    sense: ObjectiveSense,
    elapsed: std::time::Duration,
    raw_log: Option<String>,
    var_order: &[VarId],
    con_order: &[ConstraintId],
    soc_bounds: &[LinearTerms<'_>],
    model_id: ModelId,
) -> SolverResult {
    let tokens: Vec<&str> = tim.split_whitespace().collect();
    let int_at = |i: usize| tokens.get(i).and_then(|s| s.parse::<i64>().ok());
    let float_at = |i: usize| tokens.get(i).and_then(|s| parse_baron_float(s));

    let model_status = tim_model_status(tim);
    let lower = float_at(5);
    let upper = float_at(6);
    let node_count = int_at(10).and_then(|n| u64::try_from(n).ok());
    let iterations = node_count.unwrap_or(0);
    let nodeopt = int_at(11);

    let termination = map_status(res, model_status);
    // BARON has a usable point when the model status says optimal (`1`) or
    // intermediate-feasible (`4`). The solution node (`nodeopt`) only
    // controls whether a primal vector is available.
    let has_sol = matches!(model_status, 1 | 4);

    // For minimization the incumbent is the upper bound, for maximization the
    // lower bound (the other field is the dual/relaxation bound).
    let objective = match sense {
        ObjectiveSense::Minimize => upper,
        ObjectiveSense::Maximize => lower,
    };

    let mut solutions = if has_sol && nodeopt != Some(-3) {
        parse_results(res, var_order, sense)
    } else {
        Vec::new()
    };

    // BARON prints each solution's exact objective. Only when the status says a
    // solution exists but no primal block was parsed do we fall back to the times
    // file so `result_count` still reflects that a solution exists.
    if has_sol && solutions.is_empty() {
        solutions.push(SolutionPoint { model_id, primal: FxHashMap::default(), objective });
    }

    let (dual, reduced_costs, soc_prices) = if has_sol {
        parse_dual_solution(res, var_order, con_order)
    } else {
        (FxHashMap::default(), FxHashMap::default(), FxHashMap::default())
    };

    // Rescale each SOC squared-row price to the norm-form bound multiplier
    // using the bound's value at the best primal point.
    let mut soc_dual: FxHashMap<SocConstraintId, f64> = FxHashMap::default();
    if let Some(best) = solutions.first() {
        for (i, bound) in soc_bounds.iter().enumerate() {
            if let Some(price) = soc_prices.get(&i) {
                let Some(b_val) = bound.coeffs.iter().try_fold(bound.constant, |acc, &(v, c)| {
                    best.primal.get(&v).map(|value| acc + c * *value)
                }) else {
                    continue;
                };
                soc_dual.insert(
                    SocConstraintId(u32::try_from(i).expect("SOC count overflow")),
                    2.0 * b_val * price.abs(),
                );
            }
        }
    }

    let best_bound = match sense {
        ObjectiveSense::Minimize => lower,
        ObjectiveSense::Maximize => upper,
    };
    let raw_status = termination_banner(res)
        .map_or_else(|| format!("model_status={model_status}"), str::to_owned);
    let dual_available =
        has_sol && res.to_ascii_lowercase().contains("corresponding dual solution");

    normalize_result(
        SolverResult {
            model_id,
            solutions,
            dual,
            soc_dual,
            reduced_costs,
            termination,
            primal_status: PrimalStatus::NoSolution,
            dual_status: if dual_available {
                DualStatus::FeasiblePoint
            } else if has_sol {
                DualStatus::Unknown
            } else {
                DualStatus::NoSolution
            },
            best_bound,
            gap: None,
            solve_time: elapsed,
            iterations,
            node_count,
            raw_status: Some(raw_status.into()),
            raw_log,
            solver_name: Some(crate::NAME.into()),
            solver_version: baron_version(res),
        },
        var_order.len(),
    )
}

enum BaronBanner {
    Completion,
    NodeLimit,
    IterationLimit,
    TimeLimit,
    NumericError,
    Interrupted,
    Feasible,
    MemoryLimit,
    MissingBounds,
}

fn matched_banner(log: &str) -> Option<BaronBanner> {
    if log.contains("nodes in memory") {
        Some(BaronBanner::NodeLimit)
    } else if log.contains("bar iterations") {
        Some(BaronBanner::IterationLimit)
    } else if log.contains("time exceeded") {
        Some(BaronBanner::TimeLimit)
    } else if log.contains("numerically sensitive") {
        Some(BaronBanner::NumericError)
    } else if log.contains("search interrupted by user") || log.contains("access violation") {
        Some(BaronBanner::Interrupted)
    } else if log.contains("heuristic termination") {
        Some(BaronBanner::Feasible)
    } else if log.contains("insufficient memory") {
        Some(BaronBanner::MemoryLimit)
    } else if log.contains("appropriate variable bounds") {
        Some(BaronBanner::MissingBounds)
    } else if log.contains("completion") {
        Some(BaronBanner::Completion)
    } else {
        None
    }
}

/// Determine the [`TerminationStatus`] from BARON's results-file termination banner
/// (BARON manual, Section 5.2 "Termination messages, model and solver statuses").
///
/// BARON prints exactly one `*** ... ***` banner at the end of `res.lst`. We map
/// it to a termination reason. `Normal completion` carries no limit/error reason,
/// so it (and any unrecognized output) defers to the numeric model status via
/// [`model_status_termination`]. Whether a usable point exists is tracked
/// separately via [`PrimalStatus`], so a limit- or heuristic-terminated run still
/// keeps its incumbent.
fn map_status(res: &str, model_status: i64) -> TerminationStatus {
    let log = res.to_ascii_lowercase();

    match matched_banner(&log) {
        Some(BaronBanner::NodeLimit) => {
            // *** Max. allowable nodes in memory reached ***
            TerminationStatus::NodeLimit
        }
        Some(BaronBanner::IterationLimit) => {
            // *** Max. allowable BaR iterations reached ***
            TerminationStatus::IterationLimit
        }
        Some(BaronBanner::TimeLimit) => {
            // *** Max. allowable time exceeded ***
            TerminationStatus::TimeLimit
        }
        Some(BaronBanner::NumericError) => {
            // *** Problem is numerically sensitive ***
            TerminationStatus::NumericError
        }
        Some(BaronBanner::Interrupted) => {
            // *** Search interrupted by user *** / *** ... access violation ... ***
            TerminationStatus::Interrupted
        }
        Some(BaronBanner::Feasible) => {
            // *** Heuristic termination ***, feasible found, global optimality not
            // guaranteed (DeltaTerm).
            TerminationStatus::Feasible
        }
        Some(BaronBanner::MemoryLimit) => {
            // *** Insufficient Memory for Data structures ***
            TerminationStatus::MemoryLimit
        }
        Some(BaronBanner::MissingBounds) => {
            // *** User did not provide appropriate variable bounds ***, relaxation
            // bounds may be invalid, so neither globality nor infeasibility is
            // guaranteed. Feasibility, if any, is still reflected in the primal status.
            TerminationStatus::Other("baron_missing_bounds".into())
        }
        Some(BaronBanner::Completion) | None => {
            // *** Normal completion *** or no recognizable banner.
            // The model status is the authoritative outcome.
            model_status_termination(model_status)
        }
    }
}

fn termination_banner(res: &str) -> Option<&str> {
    res.lines().map(str::trim).find(|line| {
        line.starts_with("***") && matched_banner(&line.to_ascii_lowercase()).is_some()
    })
}

fn baron_version(res: &str) -> Option<Cow<'static, str>> {
    for line in res.lines() {
        let lower = line.to_ascii_lowercase();
        let Some(position) = lower.find("baron version") else {
            continue;
        };
        let version = line[position + "baron version".len()..]
            .trim()
            .trim_start_matches([':', '='])
            .split_whitespace()
            .next()?;
        let version = version.trim_matches(|c: char| !c.is_ascii_alphanumeric() && c != '.');
        if !version.is_empty() {
            return Some(version.to_owned().into());
        }
    }
    None
}

/// Outcome implied by BARON's numeric model status, used on normal completion or
/// when `res.lst` carries no recognizable termination banner. `1`=optimal,
/// `2`=infeasible, `3`=unbounded, `4`=intermediate feasible (a normal completion
/// that reports a feasible point is a solved global optimum), `5`=unknown.
fn model_status_termination(model_status: i64) -> TerminationStatus {
    match model_status {
        1 | 4 => TerminationStatus::Optimal,
        2 => TerminationStatus::Infeasible,
        3 => TerminationStatus::Unbounded,
        n => TerminationStatus::Other(format!("baron_model_status_{n}")),
    }
}

/// Parse the primal solutions from BARON's results file (`res.lst`), best first.
///
/// When BARON enumerates a pool (`NumSol > 1`) it prints the distinct points
/// after `*** Normal completion ***` as numeric blocks (see
/// [`parse_solution_pool`]). Each has its own objective, but in the order it
/// found them. We need to sort the pool by objective per the
/// optimization `sense`, so index `0` is the incumbent.
fn parse_results(res: &str, var_order: &[VarId], sense: ObjectiveSense) -> Vec<SolutionPoint> {
    let mut pool = parse_solution_pool(res, var_order);
    if pool.is_empty() {
        return parse_best_table(res).into_iter().collect();
    }
    pool.sort_by(|a, b| {
        let ord = match sense {
            ObjectiveSense::Maximize => b.objective.partial_cmp(&a.objective),
            ObjectiveSense::Minimize => a.objective.partial_cmp(&b.objective),
        };
        ord.unwrap_or(std::cmp::Ordering::Equal)
    });
    pool
}

/// Parse the enumerated solution pool BARON prints after
/// `*** Normal completion ***`:
///
/// ```text
///  >>> Objective value is:        <obj>
///  >>> Corresponding solution vector is:
///  >>> Variable no.       Value
///  >>>      1             <v1>
///  >>>      2             <v2>
/// ```
///
/// one block per distinct point, best first. Variable numbers are 1-based
/// positions in the `.bar` declaration order, so `var_order[k - 1]` is the
/// `VarId` for "Variable no. `k`". Only the post-completion region is scanned.
///
/// Returns empty when BARON printed no pool (e.g. a time-limited run).
fn parse_solution_pool(res: &str, var_order: &[VarId]) -> Vec<SolutionPoint> {
    let Some(pos) = res.find("Normal completion") else {
        return Vec::new();
    };
    let strip = |l: &str| l.trim_start().trim_start_matches(">>>").trim().to_string();

    let mut out: Vec<SolutionPoint> = Vec::new();
    let mut last_obj: Option<f64> = None;
    let mut lines = res[pos..].lines().peekable();

    while let Some(line) = lines.next() {
        let t = strip(line);
        if t.starts_with("Objective value is") {
            last_obj = t.rsplit(':').next().and_then(parse_baron_float);
        } else if t.starts_with("Corresponding dual solution")
            || t.starts_with("The best solution found")
        {
            break;
        } else if t.starts_with("Corresponding solution vector") {
            let mut primal: FxHashMap<VarId, f64> = FxHashMap::default();
            while let Some(peek) = lines.peek() {
                let p = strip(peek);
                if p.is_empty()
                    || p.starts_with("Objective value is")
                    || p.starts_with("Corresponding")
                    || p.starts_with("The best solution found")
                {
                    break;
                }
                lines.next();
                let parts: Vec<&str> = p.split_whitespace().collect();
                if let (Some(k), Some(val)) = (
                    parts.first().and_then(|s| s.parse::<usize>().ok()),
                    parts.get(1).and_then(|s| parse_baron_float(s)),
                ) && (1..=var_order.len()).contains(&k)
                {
                    primal.insert(var_order[k - 1], val);
                }
            }
            if !primal.is_empty() {
                out.push(SolutionPoint {
                    model_id: ModelId::UNASSIGNED,
                    primal,
                    objective: last_obj,
                });
            }
        }
    }
    out
}

/// `(dual, reduced_costs, soc_prices)` parsed from BARON's dual section.
/// `soc_prices` keys the raw price of each appended SOC quadratic row by cone
/// index (its sign row is skipped).
type DualParts = (FxHashMap<ConstraintId, f64>, FxHashMap<VarId, f64>, FxHashMap<usize, f64>);

/// Parse the dual solution BARON prints
///
/// Variable marginals are reduced costs, keyed by 1-based position in the
/// `.bar` declaration order (`var_order`, as in [`parse_solution_pool`]).
/// Constraint prices are duals, keyed by 1-based position in the `EQUATIONS`
/// order, which is the `ConstraintId` index plus one; positions past
/// `con_order` are the appended SOC `(quadratic row, sign row)` pairs in cone
/// order. Values are passed through with BARON's sign convention.
///
/// Returns empty maps when the section is absent (e.g. `WantDual` off, or
/// BARON reported "No dual information is available").
fn parse_dual_solution(res: &str, var_order: &[VarId], con_order: &[ConstraintId]) -> DualParts {
    let mut dual: FxHashMap<ConstraintId, f64> = FxHashMap::default();
    let mut reduced_costs: FxHashMap<VarId, f64> = FxHashMap::default();
    let mut soc_prices: FxHashMap<usize, f64> = FxHashMap::default();
    // BARON may print the dual vector more than once (e.g. for a solution
    // found during preprocessing and again after `*** Normal completion ***`);
    // the last block is the one for the final best point.
    let Some(pos) = res.rfind("Corresponding dual solution vector") else {
        return (dual, reduced_costs, soc_prices);
    };
    let strip = |l: &str| l.trim_start().trim_start_matches(">>>").trim().to_string();

    // Marginals first, then prices once the `Constraint no. Price` header is
    // seen. Both subsections are optional
    let mut in_prices = false;
    for line in res[pos..].lines().skip(1) {
        let t = strip(line);
        if t.is_empty() || t.starts_with("The best solution found") {
            break;
        }
        if t.contains("Price") {
            in_prices = true;
            continue;
        }
        if t.contains("Marginal") {
            continue;
        }
        let parts: Vec<&str> = t.split_whitespace().collect();
        let (Some(k), Some(val)) = (
            parts.first().and_then(|s| s.parse::<usize>().ok()),
            parts.get(1).and_then(|s| parse_baron_float(s)),
        ) else {
            break;
        };
        if k == 0 {
            continue;
        }
        if in_prices {
            if let Some(&id) = con_order.get(k - 1) {
                *dual.entry(id).or_insert(0.0) += val;
            } else {
                // Past the ConstraintId rows: the appended SOC pairs, two
                // rows per cone. Keep the quadratic row's price, skip the
                // sign row's.
                let rel = k - con_order.len() - 1;
                if rel.is_multiple_of(2) {
                    soc_prices.insert(rel / 2, val);
                }
            }
        } else if k <= var_order.len() {
            reduced_costs.insert(var_order[k - 1], val);
        }
    }
    (dual, reduced_costs, soc_prices)
}

/// Parse the single `"The best solution found is:"` table: rows of
/// `x{VarId} xlo xbest xup`, with the value in the third whitespace column and a
/// trailing `"... objective value of: X"` line. Returns `None` when the banner
/// is absent.
fn parse_best_table(res: &str) -> Option<SolutionPoint> {
    let pos = res.find("The best solution found")?;
    let mut primal: FxHashMap<VarId, f64> = FxHashMap::default();
    let mut objective = None;
    let mut started = false;
    for line in res[pos..].lines().skip(1) {
        if line.contains("objective value of") {
            objective = line.split_whitespace().last().and_then(parse_baron_float);
            break;
        }
        let parts: Vec<&str> = line.split_whitespace().collect();
        if let (Some(&p0), Some(&p2)) = (parts.first(), parts.get(2))
            && let (Some(idx), Some(val)) = (extract_index(p0), parse_baron_float(p2))
        {
            primal.insert(VarId(idx), val);
            started = true;
            continue;
        }
        if started {
            break;
        }
    }
    (!primal.is_empty() || objective.is_some()).then_some(SolutionPoint {
        model_id: ModelId::UNASSIGNED,
        primal,
        objective,
    })
}

/// Parse BARON's IIS listing from the results file (`res.lst`).
///
/// For an infeasible model solved with `CompIIS`, BARON prints:
///
/// ```text
///  IIS contains 1 row and 2 columns as follows:
///   c0   Upper
///   x0   Lower
///   x1   Lower
/// ```
///
/// Each entry is `<name>   <side>`. Names follow the writer's convention
/// ([`write_equations`]/[`write_var_section`]). `c{n}`/`c{n}_lo`/`c{n}_hi` are
/// constraint rows, `soc{n}`/`soc{n}_sign` are SOC rows, and `x{n}` are variable
/// columns whose `side` (`Lower`/`Upper`) selects the bound. Row entries carry a side
/// too, but it is informational. Returns an empty [`Iis`] when no IIS block is present.
///
/// The banner is matched on `"IS contains"` so both `"IIS contains ..."` (a proven
/// irreducible set) and `"(I)IS contains ..."` (BARON's label when irreducibility
/// cannot be guaranteed) are handled.
fn parse_iis(res: &str) -> Iis {
    let mut iis = Iis::default();
    let Some(pos) = res.find("IS contains") else {
        return iis;
    };
    let start = res[pos..].find("as follows:").map_or(res.len(), |a| pos + a + "as follows:".len());

    let mut seen_con: FxHashSet<u32> = FxHashSet::default();
    let mut seen_soc: FxHashSet<u32> = FxHashSet::default();
    let mut started = false;
    for line in res[start..].lines() {
        let t = line.trim();
        if t.is_empty() {
            if started {
                break;
            }
            continue;
        }
        if t.starts_with("***") {
            break;
        }
        let mut parts = t.split_whitespace();
        let Some(name) = parts.next() else { continue };
        let side = parts.next();
        // Classify by name prefix (check `x` and `soc` before the bare `c`).
        if name.starts_with('x') {
            if let Some(i) = extract_index(name) {
                iis.var_bounds.push((VarId(i), side_kind(side)));
                started = true;
            }
        } else if name.starts_with("soc") {
            if let Some(i) = extract_index(name) {
                if seen_soc.insert(i) {
                    iis.soc_constraints.push(SocConstraintId(i));
                }
                started = true;
            }
        } else if name.starts_with('c') {
            if let Some(i) = extract_index(name) {
                if seen_con.insert(i) {
                    iis.constraints.push(ConstraintId(i));
                }
                started = true;
            }
        } else if started {
            break;
        }
    }
    iis
}

/// Map BARON's IIS bound-side label (`Lower`/`Upper`) to a [`VarBoundKind`],
/// defaulting to the lower bound for anything unrecognized.
fn side_kind(side: Option<&str>) -> VarBoundKind {
    match side {
        Some(s) if s.eq_ignore_ascii_case("upper") => VarBoundKind::Upper,
        _ => VarBoundKind::Lower,
    }
}

/// Recover the variable index from a synthetic `x{i}` name by reading its digits.
fn extract_index(name: &str) -> Option<u32> {
    let digits: String =
        name.chars().skip_while(|c| !c.is_ascii_digit()).take_while(char::is_ascii_digit).collect();
    digits.parse().ok()
}

/// Parse a BARON-formatted float, tolerating its infinity sentinels and the
/// leading-decimal-point form it prints for values in `(-1, 1)` (e.g. `-.9999`).
fn parse_baron_float(s: &str) -> Option<f64> {
    match s.trim() {
        "" => None,
        "inf" | "Inf" | "+inf" | "+Inf" => Some(f64::INFINITY),
        "-inf" | "-Inf" => Some(f64::NEG_INFINITY),
        other => {
            if let Some(rest) = other.strip_prefix("-.") {
                format!("-0.{rest}").parse().ok()
            } else if let Some(rest) = other.strip_prefix('.') {
                format!("0.{rest}").parse().ok()
            } else {
                other.parse().ok()
            }
        }
    }
}

#[cfg(feature = "benchmark-support")]
#[doc(hidden)]
#[expect(clippy::cast_precision_loss)]
#[allow(clippy::wildcard_imports)]
pub mod benchmark_support {
    use oximo_core::constraint::Relate;
    use rayon::prelude::*;

    use super::*;

    /// Crossover candidate used only to size the preprocessing benchmark cases.
    pub const THRESHOLD: usize = 1_024;

    pub fn model(rows: usize, degree: usize) -> Model {
        let model = Model::new("baron_render_bench");
        let x = model.__var("x").lb(-5.0).ub(5.0).build();
        let y = model.__var("y").lb(-5.0).ub(5.0).build();
        let z = model.__var("z").lb(-5.0).ub(5.0).build();
        model.__minimize(x + y + z);
        for i in 0..rows {
            let lhs = match degree {
                1 => x + 2.0 * y - z,
                2 => x.powi(2) + y * z + x,
                _ => x * y * z + x.exp(),
            };
            model.__add_constraint_auto(lhs.le(i as f64 + 10.0));
        }
        model
    }

    pub fn render_equations(model: &Model, parallel: bool) -> Result<String, SolverError> {
        let arena = PreparedExpressions::new((*model.arena()).clone());
        let constraints = model.constraints().algebraic().to_vec();
        let socs = model.soc_constraints().clone();
        let mut out = String::new();
        let map = if parallel {
            render_parallel(&mut out, &arena, &constraints, &socs)?
        } else {
            write_equations(&mut out, &arena, &constraints, &socs, &[])?
        };
        std::hint::black_box(map);
        Ok(out)
    }

    fn render_parallel(
        out: &mut String,
        arena: &PreparedExpressions,
        constraints: &[Constraint],
        socs: &[SocConstraint],
    ) -> Result<Vec<ConstraintId>, SolverError> {
        let mut emit_map = Vec::with_capacity(constraints.len());
        if constraints.is_empty() && socs.is_empty() {
            return Ok(emit_map);
        }
        let mut names = Vec::with_capacity(constraints.len() + 2 * socs.len());
        for (i, c) in constraints.iter().enumerate() {
            let id = ConstraintId(u32::try_from(i).expect("constraint count overflow"));
            for (suffix, _, _) in equation_rows(c).into_iter().flatten() {
                names.push(format!("c{i}{suffix}"));
                emit_map.push(id);
            }
        }
        for i in 0..socs.len() {
            names.extend([format!("soc{i}"), format!("soc{i}_sign")]);
        }
        writeln!(out, "EQUATIONS {};", names.join(", ")).unwrap();
        let fragments: Vec<String> = constraints
            .par_iter()
            .enumerate()
            .map(|(i, c)| constraint_fragment(arena, i, c))
            .collect::<Result<_, _>>()?;
        for fragment in fragments {
            out.push_str(&fragment);
        }
        write_soc_rows(out, arena, socs);
        writeln!(out).unwrap();
        Ok(emit_map)
    }

    fn constraint_fragment(
        arena: &PreparedExpressions,
        index: usize,
        c: &Constraint,
    ) -> Result<String, SolverError> {
        if !expr_has_var(arena, c.lhs) {
            return Err(SolverError::Backend(format!(
                "constraint '{}' has no variables (its left-hand side is constant); BARON requires every constraint to contain at least one variable",
                c.name
            )));
        }
        let mut out = String::new();
        for (suffix, op, rhs) in equation_rows(c).into_iter().flatten() {
            write!(out, "c{index}{suffix}: ").unwrap();
            write_constraint_body(&mut out, arena, c.lhs, op, rhs)?;
        }
        Ok(out)
    }
}

#[cfg(test)]
mod tests {
    use oximo_core::prelude::*;

    use super::*;

    fn render(model: &Model) -> String {
        build_bar(model, &BaronOptions::default()).expect("build_bar").0
    }

    #[test]
    fn lp_emits_minimize_and_positive_vars() {
        let m = Model::new("lp");
        variable!(m, 0.0 <= x <= 10.0);
        variable!(m, 0.0 <= y <= 10.0);
        constraint!(m, c, x + y <= 5.0);
        objective!(m, Min, x + 2.0 * y);
        let bar = render(&m);
        assert!(bar.contains("POSITIVE_VARIABLES x0, x1;"), "{bar}");
        assert!(bar.contains("OBJ: minimize"), "{bar}");
        assert!(bar.contains("EQUATIONS c0;"), "{bar}");
        assert!(bar.contains("<= 5"), "{bar}");
        assert!(bar.contains("UPPER_BOUNDS{"), "{bar}");
    }

    #[test]
    fn convex_equation_hint_names_selected_rows_before_their_definitions() {
        let m = Model::new("convex_rows");
        variable!(m, -1.0 <= x <= 1.0);
        variable!(m, -1.0 <= y <= 1.0);
        let disk = constraint!(m, disk, x.powi(2) + y.powi(2) <= 1.0);
        constraint!(m, cap, x + y <= 1.0);
        objective!(m, Min, x + y);

        let options = BaronOptions::default().convex_equations([disk, disk]);
        let bar = build_bar(&m, &options).expect("build convex equation hint").0;
        assert!(bar.contains("EQUATIONS c0, c1;"), "{bar}");
        assert!(bar.contains("CONVEX_EQUATIONS c0;"), "{bar}");
        let declaration = bar.find("CONVEX_EQUATIONS c0;").unwrap();
        let definition = bar.find("c0:").unwrap();
        assert!(declaration < definition, "declaration must precede equation body:\n{bar}");
    }

    #[test]
    fn convex_equation_hint_rejects_unknown_id() {
        let m = Model::new("unknown_convex_row");
        variable!(m, 0.0 <= x <= 1.0);
        constraint!(m, cap, x <= 1.0);
        objective!(m, Min, x);

        let error =
            convex_equation_names(m.constraints().algebraic(), &[ConstraintId(99)]).unwrap_err();
        assert!(error.to_string().contains("unknown algebraic constraint"), "{error}");
    }

    #[test]
    fn convex_equation_hint_rejects_foreign_handle_with_colliding_id() {
        let target = Model::new("target");
        variable!(target, x);
        constraint!(target, target_row, x <= 1.0);
        objective!(target, Min, x);

        let foreign = Model::new("foreign");
        variable!(foreign, y);
        let foreign_row = constraint!(foreign, foreign_row, y <= 2.0);

        let error =
            build_bar(&target, &BaronOptions::default().convex_equation(foreign_row)).unwrap_err();
        assert!(error.to_string().contains("model mismatch"), "{error}");
        assert_eq!(foreign_row.id(), ConstraintId(0));
    }

    #[test]
    fn convex_equation_hint_rejects_free_and_two_sided_rows() {
        let m = Model::new("invalid_convex_rows");
        variable!(m, -1.0 <= x <= 1.0);
        let free = m.__add_constraint_interval("free", x, f64::NEG_INFINITY, f64::INFINITY);
        let range = m.__add_constraint_interval("range", x, -0.5, 0.5);
        objective!(m, Min, x);

        let free = m.constraint_handle_from_id(free).unwrap();
        let range = m.constraint_handle_from_id(range).unwrap();

        let free_error = build_bar(&m, &BaronOptions::default().convex_equation(free)).unwrap_err();
        assert!(free_error.to_string().contains("emits no equation"), "{free_error}");

        let range_error =
            build_bar(&m, &BaronOptions::default().convex_equation(range)).unwrap_err();
        assert!(range_error.to_string().contains("two-sided range"), "{range_error}");
    }

    #[test]
    fn convex_equation_hint_rejects_inactive_row() {
        let m = Model::new("inactive_convex_row");
        variable!(m, 0.0 <= x <= 1.0);
        constraint!(m, cap, x <= 1.0);
        let mut constraints = m.constraints().algebraic().to_vec();
        constraints[0].active = false;

        let error = convex_equation_names(&constraints, &[ConstraintId(0)]).unwrap_err();
        assert!(error.to_string().contains("inactive constraint"), "{error}");
    }

    #[test]
    fn free_variable_emits_lower_and_upper_bounds() {
        let m = Model::new("free");
        variable!(m, -5.0 <= x <= 5.0);
        objective!(m, Min, x * x);
        let bar = render(&m);
        assert!(bar.contains("VARIABLES x0;"), "{bar}");
        assert!(bar.contains("LOWER_BOUNDS{"), "{bar}");
        assert!(bar.contains("x0: -5;"), "{bar}");
        assert!(bar.contains("x0: 5;"), "{bar}");
    }

    #[test]
    fn nlp_emits_exp_and_log() {
        let m = Model::new("nlp");
        variable!(m, 0.1 <= x <= 10.0);
        objective!(m, Min, (1.0 + x).log() + x.exp());
        let bar = render(&m);
        assert!(bar.contains("log("), "{bar}");
        assert!(bar.contains("exp("), "{bar}");
    }

    #[test]
    fn minlp_partitions_binary_integer_and_continuous() {
        let m = Model::new("minlp");
        variable!(m, b, Bin);
        variable!(m, 0.0 <= n <= 5.0, Int);
        variable!(m, 0.0 <= y <= 10.0);
        constraint!(m, budget, b + n + y <= 8.0);
        objective!(m, Max, (1.0 + y).log() + 2.0 * b + n);
        let bar = render(&m);
        assert!(bar.contains("BINARY_VARIABLES x0;"), "{bar}");
        assert!(bar.contains("INTEGER_VARIABLES x1;"), "{bar}");
        assert!(bar.contains("POSITIVE_VARIABLES x2;"), "{bar}");
        assert!(bar.contains("OBJ: maximize"), "{bar}");
    }

    #[test]
    fn explicit_soc_emits_squared_rows_and_sign_row() {
        let m = Model::new("socp");
        variable!(m, -10.0 <= x <= 10.0);
        variable!(m, -10.0 <= y <= 10.0);
        variable!(m, t >= 0.0);
        m.add_soc_constraint("cone", [x, y], t);
        constraint!(m, c, x + y >= 1.0);
        objective!(m, Min, t);
        assert_eq!(m.kind(), ModelKind::SOCP);

        let bar = render(&m);
        assert!(bar.contains("EQUATIONS c0, soc0, soc0_sign;"), "declares SOC rows last:\n{bar}");
        assert!(bar.contains("soc0: "), "emits SOC row:\n{bar}");
        assert!(bar.contains(")^2"), "squares members:\n{bar}");
        assert!(bar.contains(")^2 <= 0;"), "row against 0:\n{bar}");
        assert!(bar.contains("soc0_sign: "), "emits sign row:\n{bar}");
        assert!(bar.contains(">= 0;"), "sign row nonneg:\n{bar}");
    }

    #[test]
    fn abs_reformulated_as_square_root() {
        let m = Model::new("absbar");
        variable!(m, -10.0 <= x <= 10.0);
        objective!(m, Min, x.abs());
        let bar = render(&m);
        // BARON has no abs(), reformulate |x| = (x^2)^(1/2).
        assert!(bar.contains(") ^ 2) ^ 0.5)"), "expected abs rewrite:\n{bar}");
        assert!(!bar.contains("abs("), "must not emit a literal abs():\n{bar}");
    }

    #[test]
    fn integer_power_uses_caret() {
        let m = Model::new("pow");
        variable!(m, -10.0 <= x <= 10.0);
        objective!(m, Min, x.powi(3));
        let bar = render(&m);
        assert!(bar.contains(" ^ 3)"), "expected caret power:\n{bar}");
    }

    #[test]
    fn constant_base_uses_native_caret() {
        let m = Model::new("cbpow");
        variable!(m, 0.0 <= x <= 5.0);
        let two = Expr::constant(x.arena, 2.0);
        objective!(m, Min, two.pow(x));
        let bar = render(&m);
        assert!(bar.contains("2 ^ x0"), "expected native b^x:\n{bar}");
        assert!(!bar.contains("exp("), "constant base must not rewrite to exp/log:\n{bar}");
    }

    #[test]
    fn variable_exponent_rewrites_to_exp_log() {
        let m = Model::new("vpow");
        variable!(m, 0.1 <= x <= 10.0);
        variable!(m, 0.1 <= y <= 10.0);
        objective!(m, Min, x.pow(y));
        let bar = render(&m);
        assert!(bar.contains("exp("), "{bar}");
        assert!(bar.contains("log("), "{bar}");
        assert!(!bar.contains('^'), "must not emit caret for variable exponent:\n{bar}");
    }

    #[test]
    fn added_functions_use_supported_powers_and_documented_log10_mapping() {
        let m = Model::new("added_functions");
        variable!(m, 0.1 <= x <= 10.0);
        objective!(m, Min, x.sqrt() + x.exp2() + x.log10());
        let bar = render(&m);
        assert!(bar.contains(") ^ 0.5)"), "sqrt mapping:\n{bar}");
        assert!(bar.contains("2 ^ ("), "exp2 mapping:\n{bar}");
        assert!(bar.contains("0.4342944819032518 * log("), "log10 mapping:\n{bar}");
    }

    #[test]
    fn quadratic_constraint_keeps_rhs() {
        let m = Model::new("qcp");
        variable!(m, 0.0 <= x <= 5.0);
        variable!(m, 0.0 <= y <= 5.0);
        constraint!(m, xy, x * y <= 4.0);
        objective!(m, Min, x + y);
        let bar = render(&m);
        assert!(bar.contains("x0") && bar.contains("x1"), "{bar}");
        assert!(bar.contains("<= 4;"), "{bar}");
    }

    #[test]
    fn feasibility_problem_minimizes_zero() {
        let m = Model::new("feas");
        variable!(m, 0.0 <= x <= 1.0);
        constraint!(m, c, x <= 1.0);
        objective!(m, Feasibility);
        let bar = render(&m);
        assert!(bar.contains("OBJ: minimize 0;"), "{bar}");
    }

    #[test]
    fn sin_is_rejected() {
        let m = Model::new("trig");
        variable!(m, -1.0 <= x <= 1.0);
        objective!(m, Min, x.sin());
        let err = build_bar(&m, &BaronOptions::default()).unwrap_err();
        assert!(matches!(
            err,
            SolverError::UnsupportedNonlinearOperator { backend: "BARON", operator: "sin" }
        ));
    }

    #[test]
    fn structural_extrema_are_rejected_without_automatic_reformulation() {
        let m = Model::new("extrema");
        variable!(m, -1.0 <= x <= 1.0);
        objective!(m, Min, x.min(-x));
        assert!(matches!(
            build_bar(&m, &BaronOptions::default()),
            Err(SolverError::UnsupportedNonlinearOperator { backend: "BARON", operator: "min" })
        ));
    }

    #[test]
    fn semicontinuous_is_rejected() {
        let m = Model::new("semi");
        variable!(m, x <= 10.0, SemiCont(1.0));
        objective!(m, Min, x);
        let err = build_bar(&m, &BaronOptions::default()).unwrap_err();
        assert!(matches!(err, SolverError::Backend(_)));
    }

    #[test]
    fn constant_constraint_is_rejected() {
        let m = Model::new("constc");
        variable!(m, 0.0 <= x <= 5.0);
        // x - x folds to a constant left-hand side, which BARON would reject.
        constraint!(m, trivial, x - x <= 1.0);
        objective!(m, Min, x);
        let err = build_bar(&m, &BaronOptions::default()).unwrap_err();
        match err {
            SolverError::Backend(msg) => assert!(msg.contains("no variables"), "{msg}"),
            other => panic!("expected Backend error, got {other:?}"),
        }
    }

    #[test]
    fn binary_fixed_to_one_emits_lower_bound() {
        let m = Model::new("fix1");
        variable!(m, b, Bin);
        m.fix(b, 1.0).unwrap();
        objective!(m, Min, b);
        let bar = render(&m);
        assert!(bar.contains("BINARY_VARIABLES x0;"), "{bar}");
        // lb=1 differs from the binary default 0, so it must be emitted.
        assert!(bar.contains("LOWER_BOUNDS{"), "{bar}");
        assert!(bar.contains("x0: 1;"), "fixed-to-1 binary must pin lb:\n{bar}");
    }

    #[test]
    fn binary_fixed_to_zero_emits_upper_bound() {
        let m = Model::new("fix0");
        variable!(m, b, Bin);
        m.fix(b, 0.0).unwrap();
        objective!(m, Min, b);
        let bar = render(&m);
        // ub=0 differs from the binary default 1, so it must be emitted.
        assert!(bar.contains("UPPER_BOUNDS{"), "{bar}");
        assert!(bar.contains("x0: 0;"), "fixed-to-0 binary must pin ub:\n{bar}");
    }

    #[test]
    fn starting_point_emitted_when_initial_set() {
        let m = Model::new("start");
        variable!(m, 0.0 <= x <= 10.0);
        m.set_initial(x, 3.5).unwrap();
        objective!(m, Min, x * x);
        let bar = render(&m);
        assert!(bar.contains("STARTING_POINT{"), "{bar}");
        assert!(bar.contains("x0: 3.5;"), "{bar}");
    }

    #[test]
    fn termination_from_banner() {
        // `map_status(res, model_status)`. The results-file banner
        // (manual, Section 5.2) drives the termination reason,
        // the model status is the fallback.
        let nc = "*** Normal completion ***";
        assert_eq!(map_status(nc, 1), TerminationStatus::Optimal);
        assert_eq!(map_status(nc, 2), TerminationStatus::Infeasible);
        assert_eq!(map_status(nc, 3), TerminationStatus::Unbounded);
        // Limit / interrupt banners are authoritative regardless of model status,
        // and a feasible incumbent is kept via the primal status.
        assert_eq!(
            map_status("*** Max. allowable nodes in memory reached ***", 4),
            TerminationStatus::NodeLimit
        );
        assert_eq!(
            map_status("*** Max. allowable BaR iterations reached ***", 4),
            TerminationStatus::IterationLimit
        );
        assert_eq!(
            map_status("*** Max. allowable time exceeded ***", 4),
            TerminationStatus::TimeLimit
        );
        assert_eq!(
            map_status("*** Problem is numerically sensitive ***", 4),
            TerminationStatus::NumericError
        );
        assert_eq!(
            map_status("*** Search interrupted by user ***", 4),
            TerminationStatus::Interrupted
        );
        assert_eq!(map_status("*** Heuristic termination ***", 4), TerminationStatus::Feasible);
        assert_eq!(
            map_status("*** Insufficient Memory for Data structures ***", 4),
            TerminationStatus::MemoryLimit
        );
        assert_eq!(
            map_status("*** User did not provide appropriate variable bounds ***", 4),
            TerminationStatus::Other("baron_missing_bounds".into())
        );
        // No recognizable banner falls back to the model status.
        assert_eq!(map_status("", 1), TerminationStatus::Optimal);
    }

    #[test]
    fn mapped_banners_preserve_raw_status() {
        let cases = [
            ("*** Normal completion ***", 1, TerminationStatus::Optimal),
            ("*** Max. allowable nodes in memory reached ***", 4, TerminationStatus::NodeLimit),
            ("*** Max. allowable BaR iterations reached ***", 4, TerminationStatus::IterationLimit),
            ("*** Max. allowable time exceeded ***", 4, TerminationStatus::TimeLimit),
            ("*** Problem is numerically sensitive ***", 4, TerminationStatus::NumericError),
            ("*** Search interrupted by user ***", 4, TerminationStatus::Interrupted),
            ("*** BARON stopped after an access violation ***", 4, TerminationStatus::Interrupted),
            ("*** Heuristic termination ***", 4, TerminationStatus::Feasible),
            ("*** Insufficient Memory for Data structures ***", 4, TerminationStatus::MemoryLimit),
            (
                "*** User did not provide appropriate variable bounds ***",
                4,
                TerminationStatus::Other("baron_missing_bounds".into()),
            ),
        ];

        for (banner, model_status, expected) in cases {
            assert_eq!(map_status(banner, model_status), expected);
            assert_eq!(termination_banner(banner), Some(banner));

            let tim = format!("m 0 0 0 0 0 0 1 {model_status} 0 0 0 0 0 0.0");
            let result = parse_solution(
                &tim,
                banner,
                ObjectiveSense::Minimize,
                std::time::Duration::ZERO,
                None,
                &[],
                &[],
                &[],
                ModelId::UNASSIGNED,
            );
            assert_eq!(result.termination, expected);
            assert_eq!(result.raw_status.as_deref(), Some(banner));
        }
    }

    #[test]
    fn parse_tim_picks_objective_by_sense() {
        // name ncon nvar a b lower upper solver model c iters nodeopt ... wall
        let tim = "m 1 2 0 0 1.5 9.5 1 1 0 42 7 0 0 0.42";
        let r = parse_solution(
            tim,
            "",
            ObjectiveSense::Minimize,
            std::time::Duration::ZERO,
            None,
            &[],
            &[],
            &[],
            ModelId::UNASSIGNED,
        );
        assert_eq!(r.termination, TerminationStatus::Optimal);
        assert_eq!(r.objective(), Some(9.5)); // upper bound for minimize
        assert_eq!(r.iterations, 42); // branch-and-reduce iterations from tim[10]
        assert_eq!(r.node_count, Some(42));
        assert!(r.raw_status.as_deref().unwrap().contains("model_status=1"));
        let r = parse_solution(
            tim,
            "",
            ObjectiveSense::Maximize,
            std::time::Duration::ZERO,
            None,
            &[],
            &[],
            &[],
            ModelId::UNASSIGNED,
        );
        assert_eq!(r.objective(), Some(1.5)); // lower bound for maximize
    }

    #[test]
    fn parse_res_extracts_primal() {
        let res = "\
junk line
The best solution found is:

  variable     xlo      xbest    xup
  x0           0.0      1.25     10
  x1           0.0      3.50     10

The above solution has an objective value of:    0.0
";
        let blocks = parse_results(res, &[VarId(0), VarId(1)], ObjectiveSense::Minimize);
        assert_eq!(blocks.len(), 1);
        assert_eq!(blocks[0].primal.get(&VarId(0)), Some(&1.25));
        assert_eq!(blocks[0].primal.get(&VarId(1)), Some(&3.5));
    }

    #[test]
    fn parse_res_extracts_multiple_solutions() {
        let res = "\
                         *** Normal completion ***

 >>> Objective value is:           0.0000000000000000000000
 >>> Corresponding solution vector is:
 >>> Variable no.              Value
 >>>       1             1.0000000010231691049967
 >>>       2             1.0000000010231691049967

 >>> Objective value is:           0.0000000000000000000000
 >>> Corresponding solution vector is:
 >>> Variable no.              Value
 >>>       1            -.99999999996972754878755
 >>>       2            0.99999999997690125486116

 >>> Corresponding dual solution vector is:
 >>> Variable no.              Marginal
 >>>       1             0.0000000000000000000000

The best solution found is:

  variable     xlo      xbest                          xup
  x0           -2       1.00000000102316910499667e+00  2
  x1           -2       1.00000000102316910499667e+00  2

The above solution has an objective value of:  0.0000000000000000000000
";
        let blocks = parse_results(res, &[VarId(0), VarId(1)], ObjectiveSense::Minimize);
        assert_eq!(blocks.len(), 2, "expected two solution blocks: {blocks:?}");
        assert!((blocks[0].primal[&VarId(0)] - 1.0).abs() < 1e-6);
        assert!((blocks[0].primal[&VarId(1)] - 1.0).abs() < 1e-6);
        assert!((blocks[1].primal[&VarId(0)] + 1.0).abs() < 1e-6);
        assert!((blocks[1].primal[&VarId(1)] - 1.0).abs() < 1e-6);
        assert_eq!(blocks[1].objective, Some(0.0));
    }

    #[test]
    fn parse_dual_extracts_marginals_and_prices() {
        let res = "\
                         *** Normal completion ***

 >>> Objective value is:           5.0000000000000000000000
 >>> Corresponding solution vector is:
 >>> Variable no.              Value
 >>>       1             1.0000000000000000000000
 >>>       2             4.0000000000000000000000

 >>> Corresponding dual solution vector is:
 >>> Variable no.              Marginal
 >>>       1             0.0000000000000000000000
 >>>       2            -.50000000000000000000000
 >>> Constraint no.             Price
 >>>       1             2.0000000000000000000000
 >>>       2            -.25000000000000000000000

The best solution found is:

  variable     xlo      xbest    xup
  x0           0.0      1.0      10
  x1           0.0      4.0      10

The above solution has an objective value of:    5.0
";
        let (dual, rc, _soc) =
            parse_dual_solution(res, &[VarId(0), VarId(1)], &[ConstraintId(0), ConstraintId(1)]);
        assert_eq!(rc.get(&VarId(0)), Some(&0.0));
        assert_eq!(rc.get(&VarId(1)), Some(&-0.5));
        assert_eq!(dual.get(&ConstraintId(0)), Some(&2.0));
        assert_eq!(dual.get(&ConstraintId(1)), Some(&-0.25));
        assert_eq!(dual.len(), 2);
        assert_eq!(rc.len(), 2);
    }

    #[test]
    fn parse_dual_takes_last_block_when_printed_twice() {
        let res = "\
 Solving bounding LP
 >>> Preprocessing found feasible solution
 >>> Objective value is:           6.0
 >>> Corresponding solution vector is:
 >>> Variable no.              Value
 >>>       1             6.0000000000000000000000

 >>> Corresponding dual solution vector is:
 >>> Variable no.              Marginal
 >>>       1             9.0000000000000000000000
 >>> Constraint no.            Price
 >>>       1             9.0000000000000000000000


                         *** Normal completion ***

 >>> Objective value is:           5.0
 >>> Corresponding solution vector is:
 >>> Variable no.              Value
 >>>       1             5.0000000000000000000000

 >>> Corresponding dual solution vector is:
 >>> Variable no.              Marginal
 >>>       1             0.0000000000000000000000
 >>> Constraint no.            Price
 >>>       1             1.0000000000000000000000

The best solution found is:
";
        let (dual, rc, _soc) = parse_dual_solution(res, &[VarId(0)], &[ConstraintId(0)]);
        assert_eq!(rc.get(&VarId(0)), Some(&0.0));
        assert_eq!(dual.get(&ConstraintId(0)), Some(&1.0));
    }

    #[test]
    fn parse_dual_folds_range_rows_onto_one_constraint() {
        let res = "\
 >>> Corresponding dual solution vector is:
 >>> Variable no.              Marginal
 >>>       1             0.0000000000000000000000
 >>> Constraint no.             Price
 >>>       1             2.0000000000000000000000
 >>>       2            -.50000000000000000000000
";
        let (dual, _rc, _soc) =
            parse_dual_solution(res, &[VarId(0)], &[ConstraintId(0), ConstraintId(0)]);
        assert_eq!(dual.len(), 1);
        assert_eq!(dual.get(&ConstraintId(0)), Some(&1.5)); // 2.0 + (-0.5)
    }

    #[test]
    fn parse_dual_maps_soc_rows_past_constraint_order() {
        let res = "\
 >>> Corresponding dual solution vector is:
 >>> Variable no.              Marginal
 >>>       1             0.0000000000000000000000
 >>> Constraint no.             Price
 >>>       1             2.0000000000000000000000
 >>>       2            -.75000000000000000000000
 >>>       3             0.0000000000000000000000
";
        let (dual, _rc, soc) = parse_dual_solution(res, &[VarId(0)], &[ConstraintId(0)]);
        assert_eq!(dual.get(&ConstraintId(0)), Some(&2.0));
        assert_eq!(dual.len(), 1);
        assert_eq!(soc.get(&0), Some(&-0.75));
        assert_eq!(soc.len(), 1);
    }

    #[test]
    fn free_row_from_infinite_bound_emits_no_equation() {
        let capacity = [f64::INFINITY, 10.0];
        let m = Model::new("freerow");
        variable!(m, 0.0 <= x[i in 0..2] <= 10.0);
        constraint!(m, cap[i in 0..2], x[i] <= capacity[i]);
        objective!(m, Max, x[0] + x[1]);

        let (bar, _vars, con_order, _soc) =
            build_bar(&m, &BaronOptions::default()).expect("build_bar");

        assert!(bar.contains("EQUATIONS c1;"), "{bar}");
        assert!(!bar.contains("c0"), "free row must emit nothing: {bar}");
        assert_eq!(con_order, vec![ConstraintId(1)], "{con_order:?}");
    }

    #[cfg(feature = "benchmark-support")]
    #[test]
    fn benchmark_render_uses_equation_rows_for_ranges() {
        let m = Model::new("benchmark-range");
        let x = m.__var("x").lb(-5.0).ub(5.0).build();
        m.__add_constraint_interval("free", x, f64::NEG_INFINITY, f64::INFINITY);
        m.__add_constraint_interval("range", x, 1.0, 3.0);

        let parallel = super::benchmark_support::render_equations(&m, true).unwrap();
        let serial = super::benchmark_support::render_equations(&m, false).unwrap();

        assert_eq!(parallel, serial, "benchmark rendering must preserve BARON output");
        assert!(parallel.contains("EQUATIONS c1_lo, c1_hi;"), "{parallel}");
        assert!(parallel.contains("c1_lo:") && parallel.contains("c1_hi:"), "{parallel}");
        assert!(!parallel.contains("c0"), "unconstrained row must emit nothing: {parallel}");
    }

    #[test]
    fn parse_dual_absent_section_leaves_maps_empty() {
        let res = "\
                         *** Normal completion ***

 >>> Objective value is:           5.0
 >>> Corresponding solution vector is:
 >>> Variable no.              Value
 >>>       1             1.0

No dual information is available.

The best solution found is:
";
        let (dual, rc, soc) = parse_dual_solution(res, &[VarId(0)], &[ConstraintId(0)]);
        assert!(dual.is_empty());
        assert!(rc.is_empty());
        assert!(soc.is_empty());
    }

    #[test]
    fn parse_solution_populates_duals_alongside_pool() {
        let tim = "m 1 2 0 0 0 0 1 1 0 42 7 0 0 0.42";
        let res = "\
                         *** Normal completion ***

 >>> Objective value is:           0.0000000000000000000000
 >>> Corresponding solution vector is:
 >>> Variable no.              Value
 >>>       1             1.0000000000000000000000
 >>>       2             1.0000000000000000000000

 >>> Objective value is:           0.0000000000000000000000
 >>> Corresponding solution vector is:
 >>> Variable no.              Value
 >>>       1            -.99999999996972754878755
 >>>       2            0.99999999997690125486116

 >>> Corresponding dual solution vector is:
 >>> Variable no.              Marginal
 >>>       1             0.0000000000000000000000
 >>>       2             1.5000000000000000000000
 >>> Constraint no.             Price
 >>>       1            -3.0000000000000000000000

The best solution found is:

  variable     xlo      xbest    xup
  x0           -2       1.0      2
  x1           -2       1.0      2

The above solution has an objective value of:  0.0
";
        let r = parse_solution(
            tim,
            res,
            ObjectiveSense::Minimize,
            std::time::Duration::ZERO,
            None,
            &[VarId(0), VarId(1)],
            &[ConstraintId(0)],
            &[],
            ModelId::UNASSIGNED,
        );
        assert_eq!(r.result_count(), 2);
        assert_eq!(r.reduced_costs.get(&VarId(1)), Some(&1.5));
        assert_eq!(r.dual.get(&ConstraintId(0)), Some(&-3.0));
        assert_eq!(r.dual.len(), 1);
    }

    #[test]
    fn parse_iis_extracts_rows_and_columns() {
        // Captured from a real BARON run (CompIIS: 1) on an infeasible model.
        let res = "\
 >>> Problem solved during preprocessing
 Lower bound is   infinity
 IIS contains 1 row and 2 columns as follows:
  c0   Upper
  x0   Lower
  x1   Lower


                         *** Normal completion ***
";
        let iis = parse_iis(res);
        assert_eq!(iis.constraints, vec![ConstraintId(0)]);
        assert!(iis.soc_constraints.is_empty());
        assert_eq!(
            iis.var_bounds,
            vec![(VarId(0), VarBoundKind::Lower), (VarId(1), VarBoundKind::Lower)]
        );
    }

    #[test]
    fn parse_iis_maps_ranges_and_socs_and_dedups() {
        let res = "\
 IIS contains 3 rows and 1 column as follows:
  c2_lo   Lower
  c2_hi   Upper
  soc0   Upper
  soc0_sign   Lower
  x5   Upper
";
        let iis = parse_iis(res);
        assert_eq!(iis.constraints, vec![ConstraintId(2)]);
        assert_eq!(iis.soc_constraints, vec![SocConstraintId(0)]);
        assert_eq!(iis.var_bounds, vec![(VarId(5), VarBoundKind::Upper)]);
    }

    #[test]
    fn parse_iis_handles_non_guaranteed_banner() {
        let res = "\
 >>> Problem solved during preprocessing
 Lower bound is   infinity
 Irreducibility of the IS provided cannot be guaranteed
 (I)IS contains 2 rows and 0 columns as follows:
  c0   Lower
  c1   Upper


                         *** Normal completion ***
";
        let iis = parse_iis(res);
        assert_eq!(iis.constraints, vec![ConstraintId(0), ConstraintId(1)]);
        assert!(iis.var_bounds.is_empty());
    }

    #[test]
    fn parse_iis_absent_block_is_empty() {
        let res = "*** Normal completion ***\nThe best solution found is:\n";
        assert!(parse_iis(res).is_empty());
    }

    #[test]
    fn no_solution_node_minus_three_leaves_primal_empty() {
        // model=1 (optimal) but nodeopt = -3 => no solution vector. We surface a
        // result without a point, because it carries no primal data,
        // so it must NOT count as a usable solution.
        let tim = "m 1 1 0 0 0 0 1 1 0 0 -3 0 0 0.01";
        let res = "The best solution found is:\n\n\n  x0  0  9.9\n";
        let r = parse_solution(
            tim,
            res,
            ObjectiveSense::Minimize,
            std::time::Duration::ZERO,
            None,
            &[VarId(0)],
            &[],
            &[],
            ModelId::UNASSIGNED,
        );
        assert_eq!(r.result_count(), 0);
        assert_eq!(r.primal_status, PrimalStatus::NoSolution);
        assert!(r.primal().is_none());
        assert!(r.value(VarId(0)).is_none(), "nodeopt -3 must not yield a variable value");
    }
}