#include "thermal_transient_algebra.h"

#include <algorithm>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <string>
#include <utility>

#include "fesa/assembly/thermal_assembler.h"
#include "thermal_transient_bounds.h"
#include "math/dense_blas_internal.h"
#include "math/finite_values.h"

namespace fesa::thermal_transient_internal {
namespace {
using bounds::Abs;
using bounds::Add;
using bounds::Div;
using bounds::Gamma;
using bounds::I;
using bounds::Magnitude;
using bounds::Mul;
using bounds::Roundoff;
using bounds::Sub;
using bounds::Tiny;
using bounds::Zero;

/// @brief Reports certification failure without publishing a candidate.
Status Failure(const std::string& message) {
  return Status::Failure(FailureCategory::kSolver, {{Severity::kError,
                                                     "thermal-increment-failed",
                                                     {},
                                                     "HEAT TRANSFER",
                                                     "",
                                                     message}});
}
/// @brief Checks every scalar while preserving signed temperatures and powers.
bool Finite(const Vector& vector) {
  return math_internal::IsFinite(vector);
}
/// @brief Checks aligned full operators and their coefficient-error storage.
bool SamePattern(const SparseMatrix& a, const SparseMatrix& b) {
  return a.Rows() == b.Rows() && a.Columns() == b.Columns() &&
         a.RowOffsets() == b.RowOffsets() &&
         a.ColumnIndices() == b.ColumnIndices();
}
/// @brief Encloses a componentwise maximum for infinity norms.
I Max(I a, I b) {
  return {std::max(a.lower, b.lower), std::max(a.upper, b.upper)};
}
/// @brief Encloses the exact infinity norm of stored scalar entries.
I Norm(const Vector& v) {
  I result;
  for (std::size_t i = 0; i < v.Size(); ++i)
    result = Max(result, Magnitude(v[i]));
  return result;
}
/// @brief Bounds construction error against the independently enclosed value.
double Error(I exact, double rounded) {
  return Abs(Sub(exact, I{rounded})).upper;
}
/// @brief Checks optional entrywise coefficient bounds without inventing
/// defaults.
bool Nonnegative(const SparseMatrix* error, const SparseMatrix& matrix) {
  if (!error) return true;
  if (!error->Validate().IsOk() || !SamePattern(*error, matrix)) return false;
  for (const double x : error->Values())
    if (x < 0) return false;
  return true;
}
/// @brief Validates the private generic algebra boundary before indexing.
Status Validate(const ThermalSystem& s, const ThermalPartition& p,
                const Vector& old, double delta) {
  const auto m = s.k.Rows();
  if (!s.k.Validate().IsOk() || !s.c.Validate().IsOk() ||
      !SamePattern(s.k, s.c) || s.k.Columns() != m || s.q.Size() != m ||
      old.Size() != m || !Finite(s.q) || !Finite(old) ||
      !Finite(p.prescribed) || !std::isfinite(delta) || delta <= 0 ||
      p.prescribed.Size() != p.constrained.size() ||
      p.free.size() + p.constrained.size() != m ||
      !Nonnegative(s.k_error, s.k) || !Nonnegative(s.c_error, s.c))
    return Failure(
        "Invalid finite thermal operators, partition, or actual increment.");
  if (s.q_error) {
    if (s.q_error->Size() != m || !Finite(*s.q_error))
      return Failure("Invalid power error bounds.");
    for (std::size_t i = 0; i < m; ++i)
      if ((*s.q_error)[i] < 0) return Failure("Negative power error bound.");
  }
  std::vector<bool> used(m, false);
  for (const auto* rows : {&p.free, &p.constrained}) {
    if (!std::is_sorted(rows->begin(), rows->end()))
      return Failure("Partition order is invalid.");
    for (const auto row : *rows) {
      if (row >= m || used[row])
        return Failure("Partition is not disjoint and exhaustive.");
      used[row] = true;
    }
  }
  return Status::Ok();
}

/// @brief Keeps rounded solve data and separate exact-algebra enclosures.
struct Equations {
  SparseMatrix a;
  Vector b;
  std::vector<I> exact_a;
  std::vector<I> exact_b;
};
/// @brief Forms the reduced operator and raw-old-state backward Euler RHS.
Result<Equations> Build(const ThermalSystem& s, const ThermalPartition& p,
                        const Vector& old, double delta) {
  const auto n = p.free.size(), m = old.Size();
  std::vector<std::size_t> free(m, n), constrained(m, p.constrained.size());
  for (std::size_t i = 0; i < n; ++i) free[p.free[i]] = i;
  for (std::size_t i = 0; i < p.constrained.size(); ++i)
    constrained[p.constrained[i]] = i;
  SparsePattern pattern{{0}, {}};
  std::vector<CooContribution> coo;
  const auto entries = s.c.Values().size();
  Vector capacity_diagonal{entries}, conductivity_diagonal{entries};
  Vector capacity_input{entries}, conductivity_input{entries};
  Vector b{n};
  std::vector<I> exact_a, exact_b(n);
  for (std::size_t f = 0; f < n; ++f) {
    const auto row = p.free[f];
    b[f] = s.q[row];
    exact_b[f] = I{s.q[row]};
    for (auto j = s.c.RowOffsets()[row]; j < s.c.RowOffsets()[row + 1]; ++j) {
      const auto col = s.c.ColumnIndices()[j];
      const double c = s.c.Values()[j], k = s.k.Values()[j];
      const I cd = Div(I{c}, I{delta});
      capacity_diagonal[j] = c / delta;
      if (free[col] != n) {
        capacity_input[j] = old[col];
        coo.push_back({f, free[col], c / delta + k, 0, coo.size()});
        pattern.column_indices.push_back(free[col]);
        exact_a.push_back(Add(cd, I{k}));
        exact_b[f] = Add(exact_b[f], Mul(cd, I{old[col]}));
      } else {
        const double g = p.prescribed[constrained[col]];
        capacity_input[j] = old[col] - g;
        conductivity_diagonal[j] = k;
        conductivity_input[j] = g;
        // Preserve raw old constrained capacity; never project frame zero.
        exact_b[f] = Add(exact_b[f],
                         Sub(Mul(cd, Sub(I{old[col]}, I{g})), Mul(I{k}, I{g})));
      }
    }
    pattern.row_offsets.push_back(coo.size());
  }
  const auto capacity_action = dense_blas_internal::MultiplyDiagonal(
      capacity_diagonal, capacity_input);
  const auto conductivity_action = dense_blas_internal::MultiplyDiagonal(
      conductivity_diagonal, conductivity_input);
  for (std::size_t f = 0; f < n; ++f) {
    const auto row = p.free[f];
    for (auto j = s.c.RowOffsets()[row]; j < s.c.RowOffsets()[row + 1]; ++j) {
      const auto col = s.c.ColumnIndices()[j];
      // Preserve the original per-entry subtraction before the row sum.
      // Separately summing C*delta and K*g loses small stable RHS terms.
      b[f] += free[col] != n ? capacity_action[j]
                             : capacity_action[j] - conductivity_action[j];
    }
  }
  auto a = SparseMatrix::FromCoo(n, n, std::move(coo), pattern);
  if (!a.HasValue()) return Result<Equations>::Failure(a.GetStatus());
  if (!Finite(b))
    return Result<Equations>::Failure(Failure("Nonfinite backward Euler RHS."));
  return Result<Equations>::Success({std::move(a.Value()), std::move(b),
                                     std::move(exact_a), std::move(exact_b)});
}

/// @brief Certifies a residual using the reduced denominator and exact-zero
/// rule.
bool Certificate(double numerator, I denominator, double limit, double& beta) {
  if (!std::isfinite(numerator) || !denominator.IsFinite()) return false;
  if (denominator.lower <= 0) {
    beta = 0;
    return Zero(denominator) && numerator == 0;
  }
  beta = Div(I{numerator}, I{denominator.lower}).upper;
  return std::isfinite(beta) && beta <= limit;
}

/// @brief Checks reduced/full residuals, physical heat balance, and BE
/// dissipation.
/// @note Coefficient errors enter forward evidence only, never the solve
/// tolerance.
Status Check(const ThermalSystem& s, const ThermalPartition& p,
             const Vector& old, double delta, const Equations& eq,
             ThermalStepCandidate& candidate) {
  const auto m = old.Size(), n = p.free.size();
  const double md = static_cast<double>(m), nd = static_cast<double>(n);
  auto& ev = candidate.evidence;
  const I beta_limit =
      Add(Mul(I{nd}, Gamma(nd + 1)), Add(Gamma(2 * nd + 4), Gamma(4 * nd + 8)));
  ev.beta_limit = beta_limit.lower;
  I a_norm, t_norm, b_norm = Norm(eq.b);
  Vector red{n}, red_error{n};
  std::vector<I> construction(n);
  Vector free_temperature{n};
  for (std::size_t f = 0; f < n; ++f)
    free_temperature[f] = candidate.temperature[p.free[f]];
  const auto reduced_action = eq.a.Multiply(free_temperature);
  for (std::size_t f = 0; f < n; ++f) {
    I row_norm, scale, ea, identity;
    for (auto j = eq.a.RowOffsets()[f]; j < eq.a.RowOffsets()[f + 1]; ++j) {
      const double a = eq.a.Values()[j],
                   t = candidate.temperature[p.free[eq.a.ColumnIndices()[j]]];
      row_norm = Add(row_norm, Magnitude(a));
      scale = Add(scale, Mul(Magnitude(a), Magnitude(t)));
      const I error{Error(eq.exact_a[j], a)};
      ea = Add(ea, error);
      identity = Add(identity, Mul(Sub(eq.exact_a[j], I{a}), I{t}));
    }
    a_norm = Max(a_norm, row_norm);
    t_norm = Max(t_norm, Magnitude(candidate.temperature[p.free[f]]));
    red[f] = reduced_action[f] - eq.b[f];
    scale = Add(scale, Magnitude(eq.b[f]));
    const double nu =
        static_cast<double>(eq.a.RowOffsets()[f + 1] - eq.a.RowOffsets()[f]);
    red_error[f] = Roundoff(2 * nu + 2, scale).upper;
    ev.free_residual_cert = std::max(
        ev.free_residual_cert, Add(Magnitude(red[f]), I{red_error[f]}).upper);
    ev.operator_construction_error =
        std::max(ev.operator_construction_error, ea.upper);
    const double eb = Error(eq.exact_b[f], eq.b[f]);
    ev.rhs_construction_error = std::max(ev.rhs_construction_error, eb);
    construction[f] = Sub(identity, Sub(eq.exact_b[f], I{eq.b[f]}));
  }
  const I denominator = Add(Mul(a_norm, t_norm), b_norm);
  if (n && !Certificate(ev.free_residual_cert, denominator, ev.beta_limit,
                        ev.beta_cert))
    return Failure(
        "Reduced residual is not certified on the reduced-system scale.");

  Vector change{m}, full_error{m};
  for (std::size_t i = 0; i < m; ++i)
    change[i] = candidate.temperature[i] - old[i];
  if (!Finite(change)) return Failure("Nonfinite temperature change.");
  const auto capacity_change = s.c.Multiply(change);
  const auto conductivity_new = s.k.Multiply(candidate.temperature);
  I c_change, k_new, q_abs, c_temperatures, underflow_sum;
  std::vector<I> k_columns(m);
  double energy = 0;
  for (std::size_t i = 0; i < m; ++i) {
    const double cd = capacity_change[i], kt = conductivity_new[i];
    double ksum = 0;
    I cs, ks, ts, cabs, nullscale, nullerror, coefficient;
    for (auto j = s.c.RowOffsets()[i]; j < s.c.RowOffsets()[i + 1]; ++j) {
      const auto col = s.c.ColumnIndices()[j];
      const double c = s.c.Values()[j], k = s.k.Values()[j],
                   t = candidate.temperature[col];
      ksum += k;
      cs = Add(cs, Mul(Magnitude(c), Magnitude(change[col])));
      ks = Add(ks, Mul(Magnitude(k), Magnitude(t)));
      ts = Add(ts, Mul(Magnitude(c), Add(Magnitude(t), Magnitude(old[col]))));
      cabs = Add(cabs, Magnitude(c));
      nullscale = Add(nullscale, Magnitude(k));
      k_columns[col] = Add(k_columns[col], I{k});
      if (s.k_error) {
        nullerror = Add(nullerror, I{s.k_error->Values()[j]});
        coefficient =
            Add(coefficient, Mul(I{s.k_error->Values()[j]}, Magnitude(t)));
      }
      if (s.c_error)
        coefficient =
            Add(coefficient,
                Div(Mul(I{s.c_error->Values()[j]}, Magnitude(change[col])),
                    I{delta}));
    }
    energy += cd;
    candidate.residual[i] = (cd / delta + kt) - s.q[i];
    const double nu =
        static_cast<double>(s.c.RowOffsets()[i + 1] - s.c.RowOffsets()[i]);
    const I scale = Add(Div(cs, I{delta}), Add(ks, Magnitude(s.q[i])));
    // Propagate underflow in the capacity dot product through division, and
    // underflow in the old/new subtraction through each capacity coefficient.
    I under = Add(Div(Tiny(2 * nu + 2, cs), I{delta}),
                  Tiny(2 * nu + 4, Add(ks, Magnitude(s.q[i]))));
    under = Add(under, Div(Mul(cabs, Tiny(1, ts)), I{delta}));
    full_error[i] = Add(Add(Mul(Gamma(4 * nu + 6), scale),
                            Div(Mul(Gamma(1), ts), I{delta})),
                        under)
                        .upper;
    if (s.q_error) coefficient = Add(coefficient, I{(*s.q_error)[i]});
    candidate.certified_residual[i] =
        Add(Add(Magnitude(candidate.residual[i]), I{full_error[i]}),
            coefficient)
            .upper;
    const I nullbound = Add(nullerror, Roundoff(nu, nullscale));
    if (!nullbound.IsFinite() || std::abs(ksum) > nullbound.upper)
      return Failure("Conductivity does not certify its constant null mode.");
    c_change = Add(c_change, cs);
    k_new = Add(k_new, ks);
    q_abs = Add(q_abs, Magnitude(s.q[i]));
    c_temperatures = Add(c_temperatures, ts);
    underflow_sum = Add(underflow_sum, under);
  }
  if (!Finite(candidate.residual) || !Finite(candidate.certified_residual) ||
      !std::isfinite(energy))
    return Failure("Nonfinite full residual or residual enclosure.");
  for (std::size_t f = 0; f < n; ++f) {
    const auto row = p.free[f];
    const I direct =
        Add(I{candidate.residual[row]}, I{-full_error[row], full_error[row]});
    const I transferred =
        Add(Add(I{red[f]}, I{-red_error[f], red_error[f]}), construction[f]);
    const I intersection{std::max(direct.lower, transferred.lower),
                         std::min(direct.upper, transferred.upper)};
    if (!direct.IsFinite() || !transferred.IsFinite() ||
        !intersection.IsFinite())
      return Failure(
          "Independent full residual enclosures have empty intersection.");
    ev.full_free_residual_cert =
        std::max(ev.full_free_residual_cert, Abs(intersection).upper);
    const I bound =
        Add(Abs(construction[f]), Add(I{red_error[f]}, I{full_error[row]}));
    ev.construction_identity_bound =
        std::max(ev.construction_identity_bound, bound.upper);
    if (!bound.IsFinite() ||
        Abs(Sub(I{candidate.residual[row]}, I{red[f]})).lower > bound.upper)
      return Failure(
          "Full/reduced equation construction identity is not enclosed.");
  }
  if (n && !Certificate(ev.full_free_residual_cert, denominator, ev.beta_limit,
                        ev.full_beta_cert))
    return Failure(
        "Full free residual is not certified on the reduced-system scale.");

  double qsum = 0, rsum = 0;
  I rabs;
  for (std::size_t i = 0; i < m; ++i) qsum += s.q[i];
  for (std::size_t i = 0; i < p.constrained.size(); ++i) {
    candidate.reaction[i] = candidate.residual[p.constrained[i]];
    rsum += candidate.reaction[i];
    rabs = Add(rabs, Magnitude(candidate.reaction[i]));
  }
  I column_norm;
  for (const auto x : k_columns) column_norm = Add(column_norm, Abs(x));
  const I se = Add(c_change, Mul(I{delta}, Add(q_abs, rabs)));
  I be = Mul(Mul(Mul(I{delta}, I{nd}), beta_limit), denominator);
  be = Add(be, Mul(Mul(I{delta}, column_norm), Norm(candidate.temperature)));
  be = Add(be, Roundoff(4 * md + 12, se));
  be = Add(be, Mul(Mul(I{delta}, Gamma(2 * md + 6)),
                   Add(Div(c_change, I{delta}), Add(k_new, q_abs))));
  be = Add(be, Mul(Gamma(1), c_temperatures));
  be = Add(be, Mul(I{delta}, underflow_sum));
  ev.energy_change = energy;
  ev.energy_defect = std::abs(energy - delta * (qsum + rsum));
  ev.energy_bound = be.upper;
  if (!be.IsFinite() || !std::isfinite(ev.energy_defect) ||
      ev.energy_defect > be.upper)
    return Failure(
        "Physical heat balance failed independently of the free residual.");

  if (p.constrained.empty() && Zero(q_abs)) {
    ev.dissipation_checked = true;
    double cn = 0, co = 0, cd = 0, kn = 0;
    const auto capacity_new = s.c.Multiply(candidate.temperature);
    const auto capacity_old = s.c.Multiply(old);
    I sc, subscale, t_one, d_square, o_square, skew_square, quadratic_underflow;
    for (std::size_t i = 0; i < m; ++i) {
      const double t = candidate.temperature[i];
      const double row_cn = capacity_new[i], row_co = capacity_old[i],
                   row_cd = capacity_change[i], row_kn = conductivity_new[i];
      I scale_cn, scale_co, scale_cd, scale_kn;
      t_one = Add(t_one, Magnitude(t));
      d_square = Add(d_square, Mul(Magnitude(change[i]), Magnitude(change[i])));
      o_square = Add(o_square, Mul(Magnitude(old[i]), Magnitude(old[i])));
      for (auto j = s.c.RowOffsets()[i]; j < s.c.RowOffsets()[i + 1]; ++j) {
        const auto col = s.c.ColumnIndices()[j];
        const double c = s.c.Values()[j], k = s.k.Values()[j];
        scale_cn = Add(
            scale_cn, Mul(Magnitude(c), Magnitude(candidate.temperature[col])));
        scale_co = Add(scale_co, Mul(Magnitude(c), Magnitude(old[col])));
        scale_cd = Add(scale_cd, Mul(Magnitude(c), Magnitude(change[col])));
        scale_kn = Add(
            scale_kn, Mul(Magnitude(k), Magnitude(candidate.temperature[col])));
        const I nc = Mul(Mul(Magnitude(t), Magnitude(c)),
                         Magnitude(candidate.temperature[col]));
        const I oc =
            Mul(Mul(Magnitude(old[i]), Magnitude(c)), Magnitude(old[col]));
        const I dc = Mul(Mul(Magnitude(change[i]), Magnitude(c)),
                         Magnitude(change[col]));
        sc = Add(sc, Mul(I{.5}, Add(nc, Add(oc, dc))));
        sc = Add(sc, Mul(I{delta}, Mul(Mul(Magnitude(t), Magnitude(k)),
                                       Magnitude(candidate.temperature[col]))));
        subscale =
            Add(subscale,
                Mul(Mul(Add(Magnitude(t), Magnitude(old[i])), Magnitude(c)),
                    Magnitude(change[col])));
        const auto begin = s.c.ColumnIndices().begin() +
                           static_cast<std::ptrdiff_t>(s.c.RowOffsets()[col]);
        const auto end = s.c.ColumnIndices().begin() +
                         static_cast<std::ptrdiff_t>(s.c.RowOffsets()[col + 1]);
        const auto at = std::lower_bound(begin, end, i);
        const double ct = (at != end && *at == i)
                              ? s.c.Values()[static_cast<std::size_t>(
                                    at - s.c.ColumnIndices().begin())]
                              : 0;
        const I skew = c == ct ? I{0} : Abs(Sub(I{c}, I{ct}));
        skew_square = Add(skew_square, Mul(skew, skew));
      }
      cn += t * row_cn;
      co += old[i] * row_co;
      cd += change[i] * row_cd;
      kn += t * row_kn;
      const double operations =
          2 *
          static_cast<double>(s.c.RowOffsets()[i + 1] - s.c.RowOffsets()[i]);
      I c_under =
          Add(Mul(Magnitude(t), Tiny(operations, scale_cn)),
              Add(Mul(Magnitude(old[i]), Tiny(operations, scale_co)),
                  Mul(Magnitude(change[i]), Tiny(operations, scale_cd))));
      c_under = Add(c_under, Tiny(6, Add(scale_cn, Add(scale_co, scale_cd))));
      const I k_under =
          Add(Mul(Magnitude(t), Tiny(operations, scale_kn)), Tiny(2, scale_kn));
      quadratic_underflow =
          Add(quadratic_underflow,
              Add(Mul(I{.5}, c_under), Mul(I{delta}, k_under)));
    }
    I bc = Mul(Mul(Mul(I{delta}, t_one), beta_limit), denominator);
    bc = Add(bc, Roundoff(6 * md + 20, sc));
    // Frobenius bounds the requested spectral skew norm without dense storage.
    if (!Zero(skew_square))
      bc = Add(bc, Mul(I{.5}, Mul(thermal_internal::Sqrt(d_square),
                                  Mul(thermal_internal::Sqrt(skew_square),
                                      thermal_internal::Sqrt(o_square)))));
    bc = Add(bc, Mul(Gamma(1), subscale));
    bc = Add(bc, quadratic_underflow);
    ev.dissipation_defect = std::abs(.5 * ((cn - co) + cd) + delta * kn);
    ev.dissipation_bound = bc.upper;
    if (!bc.IsFinite() || !std::isfinite(ev.dissipation_defect) ||
        ev.dissipation_defect > bc.upper || !std::isfinite(cn) ||
        !std::isfinite(co) || Sub(I{cn}, I{co}).lower > Mul(I{2}, bc).upper)
      return Failure(
          "Insulated backward Euler capacity-norm dissipation failed.");
  }
  return Status::Ok();
}
}  // namespace

Result<ThermalStepCandidate> SolveThermalIncrement(
    const ThermalSystem& s, const ThermalPartition& p, const Vector& old,
    double delta, LinearSolver& solver, ThermalFactorization& cache) try {
  auto status = Validate(s, p, old, delta);
  if (!status.IsOk()) return Result<ThermalStepCandidate>::Failure(status);
  auto equations = Build(s, p, old, delta);
  if (!equations.HasValue())
    return Result<ThermalStepCandidate>::Failure(equations.GetStatus());
  const auto n = p.free.size(), m = old.Size();
  ThermalStepCandidate candidate{
      Vector{m}, Vector{m}, Vector{p.constrained.size()}, Vector{m}, {}};
  for (std::size_t i = 0; i < p.constrained.size(); ++i)
    candidate.temperature[p.constrained[i]] = p.prescribed[i];
  if (n) {
    if (!cache.matrix || cache.k != &s.k || cache.c != &s.c ||
        cache.solver != &solver || cache.free != p.free ||
        cache.delta != delta) {
      cache = ThermalFactorization{};
      cache.matrix = equations.Value().a;
      status = solver.Factorize(*cache.matrix);
      if (!status.IsOk()) {
        cache = ThermalFactorization{};
        return Result<ThermalStepCandidate>::Failure(status);
      }
      cache.k = &s.k;
      cache.c = &s.c;
      cache.solver = &solver;
      cache.free = p.free;
      cache.delta = delta;
    }
    Vector free{n};
    status = solver.Solve(equations.Value().b, free);
    if (!status.IsOk()) return Result<ThermalStepCandidate>::Failure(status);
    if (free.Size() != n || !Finite(free))
      return Result<ThermalStepCandidate>::Failure(
          Failure("Solver returned invalid temperatures."));
    for (std::size_t i = 0; i < n; ++i)
      candidate.temperature[p.free[i]] = free[i];
  }
  status = Check(s, p, old, delta, equations.Value(), candidate);
  if (!status.IsOk()) return Result<ThermalStepCandidate>::Failure(status);
  return Result<ThermalStepCandidate>::Success(std::move(candidate));
} catch (const std::invalid_argument&) {
  return Result<ThermalStepCandidate>::Failure(
      Failure("Nonfinite backward Euler RHS."));
} catch (const std::length_error&) {
  return Result<ThermalStepCandidate>::Failure(
      Failure("Thermal sparse action exceeds the backend integer range."));
} catch (const std::runtime_error&) {
  return Result<ThermalStepCandidate>::Failure(
      Failure("Thermal sparse action failed."));
}

Status RunThermalHistory(const ThermalAssembly& assembly,
                         const DofManager& dofs, ThermalAnalysisState& state,
                         LinearSolver& solver, const ThermalRecovery& recovery,
                         std::vector<ThermalStepEvidence>* evidence) {
  const BackwardEulerHeatIntegrator integrator;
  return RunThermalHistory(assembly, dofs, state, solver, integrator, recovery,
                            evidence);
}

Status RunThermalHistory(const ThermalAssembly& assembly,
                         const DofManager& dofs, ThermalAnalysisState& state,
                         LinearSolver& solver, const HeatStepIntegrator& integrator,
                         const ThermalRecovery& recovery,
                         std::vector<ThermalStepEvidence>* evidence) {
  if (!dofs.ValidateInvariants().IsOk() ||
      dofs.Profile() != DofProfile::kThermal ||
      dofs.FullDofCount() != state.CommittedState().temperature.Size() ||
      assembly.ElementOrder() != state.ElementOrder() || !recovery)
    return Failure("History does not match the thermal operators and mapping.");
  ThermalSystem system{assembly.Conductivity(),   assembly.Capacity(),
                       assembly.HeatInput(),      &assembly.ConductivityError(),
                       &assembly.CapacityError(), &assembly.HeatInputError()};
  ThermalPartition partition{dofs.FreeDofs(), dofs.ConstrainedDofs(),
                             dofs.PrescribedValues()};
  ThermalFactorization cache;
  const std::uint64_t count = state.TimeGrid().Count();
  for (std::uint64_t i =
           static_cast<std::uint64_t>(state.AcceptedFrames().size());
       i <= count; ++i) {
    const auto frame = state.TimeGrid().Frame(static_cast<std::uint32_t>(i));
    if (!frame.HasValue()) return frame.GetStatus();
    const Status step_status = RunThermalStep(
        integrator, system, partition, state.CommittedState().temperature,
        frame.Value().delta, solver, cache, recovery,
        [&](ThermalStepCandidate&& candidate,
            std::vector<ThermalElementFlux>&& flux) -> Status {
          auto completed = state.MakeCompletedFrame(
              static_cast<std::uint32_t>(i), std::move(candidate.temperature),
              std::move(candidate.reaction), std::move(flux));
          if (!completed.HasValue()) return completed.GetStatus();
          const auto committed = state.CommitAcceptedFrame(completed.Value());
          if (!committed.IsOk()) return committed;
          if (evidence) evidence->push_back(candidate.evidence);
          return Status::Ok();
        });
    if (!step_status.IsOk()) return step_status;
  }
  return Status::Ok();
}
}  // namespace fesa::thermal_transient_internal