mirror of
https://github.com/storytold/gridcraft.git
synced 2026-10-29 21:26:48 -06:00
refactor(engine): move cmd/solver.rs items to cmd/solver/
This commit is contained in:
aiongg
+15
-1848
File diff suppressed because it is too large
Load Diff
@@ -0,0 +1,255 @@
|
||||
use super::lp::{Outcome, Row, Sense};
|
||||
use super::model::{EvoOptions, Goal};
|
||||
use super::nlp::{project, shape, violation};
|
||||
use super::problem::Evaluator;
|
||||
use super::*;
|
||||
|
||||
// ---------------------------------------------------------------- evolutionary
|
||||
|
||||
/// A small deterministic generator (SplitMix64), so a seed reproduces a run on every platform.
|
||||
struct Rng(u64);
|
||||
|
||||
impl Rng {
|
||||
fn new(seed: u64) -> Rng {
|
||||
Rng(seed ^ 0x2545_F491_4F6C_DD1D)
|
||||
}
|
||||
fn next(&mut self) -> u64 {
|
||||
self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
|
||||
let mut z = self.0;
|
||||
z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
|
||||
z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
|
||||
z ^ (z >> 31)
|
||||
}
|
||||
/// Uniform in [0, 1).
|
||||
fn uniform(&mut self) -> f64 {
|
||||
(self.next() >> 11) as f64 / (1u64 << 53) as f64
|
||||
}
|
||||
fn between(&mut self, a: f64, b: f64) -> f64 {
|
||||
a + (b - a) * self.uniform()
|
||||
}
|
||||
fn below(&mut self, n: usize) -> usize {
|
||||
if n == 0 { 0 } else { usize::try_from(self.next() % n as u64).unwrap_or(0) }
|
||||
}
|
||||
/// Standard normal (Box–Muller).
|
||||
fn normal(&mut self) -> f64 {
|
||||
let u = self.uniform().max(f64::MIN_POSITIVE);
|
||||
(-2.0 * u.ln()).sqrt() * (std::f64::consts::TAU * self.uniform()).cos()
|
||||
}
|
||||
}
|
||||
|
||||
/// A candidate: its point, the objective to minimise and its constraint violation.
|
||||
#[derive(Clone, Debug)]
|
||||
struct Member {
|
||||
x: Vec<f64>,
|
||||
phi: f64,
|
||||
viol: f64,
|
||||
}
|
||||
|
||||
impl Member {
|
||||
/// Feasible beats infeasible; then the lower objective, or the smaller violation.
|
||||
fn better(&self, o: &Member, tol: f64) -> bool {
|
||||
match (self.viol <= tol, o.viol <= tol) {
|
||||
(true, false) => true,
|
||||
(false, true) => false,
|
||||
(true, true) => self.phi < o.phi,
|
||||
(false, false) => self.viol < o.viol,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Why the Evolutionary search stopped.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub(super) enum EvoEnd {
|
||||
/// Almost all of the population is within Convergence of the best member.
|
||||
Converged,
|
||||
/// The best member didn't improve for the set number of generations.
|
||||
Stalled,
|
||||
/// The evaluation budget ran out.
|
||||
Budget,
|
||||
}
|
||||
|
||||
/// The search box, integer flags and scoring shared by the genetic and pattern phases.
|
||||
struct EvoSearch<'a> {
|
||||
ev: &'a mut Evaluator,
|
||||
rows: Vec<Row>,
|
||||
senses: Vec<Sense>,
|
||||
goal: Goal,
|
||||
lo: Vec<f64>,
|
||||
hi: Vec<f64>,
|
||||
integer: &'a [bool],
|
||||
tol: f64,
|
||||
}
|
||||
|
||||
impl EvoSearch<'_> {
|
||||
/// Clamps into the box and rounds integer variables.
|
||||
fn repair(&self, x: &mut [f64]) {
|
||||
for (j, v) in x.iter_mut().enumerate() {
|
||||
let (l, h) = (self.lo.get(j).copied().unwrap_or(0.0), self.hi.get(j).copied().unwrap_or(0.0));
|
||||
*v = if v.is_finite() { v.clamp(l, h) } else { l };
|
||||
if self.integer.get(j).copied().unwrap_or(false) {
|
||||
*v = v.round().clamp(l.ceil(), h.floor());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Scores `x`; `None` once the budget is spent.
|
||||
fn score(&mut self, x: Vec<f64>) -> Option<Member> {
|
||||
let p = self.ev.eval(&x)?;
|
||||
let (phi, cons) = shape(&p, &self.rows, self.goal);
|
||||
let viol = violation(&cons, &self.senses);
|
||||
let phi = if phi.is_finite() { phi } else { f64::INFINITY };
|
||||
let viol = if viol.is_finite() && cons.iter().all(|c| c.is_finite()) { viol } else { f64::INFINITY };
|
||||
Some(Member { x, phi, viol })
|
||||
}
|
||||
}
|
||||
|
||||
/// The Evolutionary method within `lo..hi` from `x0`.
|
||||
pub(super) fn evolutionary(
|
||||
ev: &mut Evaluator,
|
||||
x0: &[f64],
|
||||
lo: &[f64],
|
||||
hi: &[f64],
|
||||
integer: &[bool],
|
||||
goal: Goal,
|
||||
o: &EvoOptions,
|
||||
) -> (Outcome, EvoEnd) {
|
||||
let rows = ev.rows.clone();
|
||||
let senses: Vec<Sense> = rows.iter().map(|r| r.sense).collect();
|
||||
let scale = rows.iter().map(|r| r.rhs_const.abs()).fold(1.0, f64::max);
|
||||
let mut start = x0.to_vec();
|
||||
project(&mut start, lo, hi);
|
||||
// Without a bound, search within 100 of the start.
|
||||
let blo: Vec<f64> = lo.iter().zip(&start).map(|(l, s)| if l.is_finite() { *l } else { s - 100.0 }).collect();
|
||||
let bhi: Vec<f64> = hi.iter().zip(&start).map(|(h, s)| if h.is_finite() { *h } else { s + 100.0 }).collect();
|
||||
let no_integer = blo.iter().zip(&bhi).zip(integer).any(|((l, h), int)| *int && l.ceil() > h.floor());
|
||||
if blo.iter().zip(&bhi).any(|(l, h)| l > h) || no_integer {
|
||||
return (Outcome::Infeasible, EvoEnd::Converged);
|
||||
}
|
||||
let span: Vec<f64> = blo.iter().zip(&bhi).map(|(l, h)| (h - l).max(1e-9)).collect();
|
||||
let mut se = EvoSearch { ev, rows, senses, goal, lo: blo, hi: bhi, integer, tol: PRECISION * scale };
|
||||
let mut rng = Rng::new(o.seed);
|
||||
let size = o.population.clamp(MIN_POPULATION, MAX_POPULATION);
|
||||
let tol = se.tol;
|
||||
let mut end = EvoEnd::Budget;
|
||||
// The first generation: the start and random points in the box.
|
||||
let mut pop: Vec<Member> = Vec::with_capacity(size);
|
||||
se.repair(&mut start);
|
||||
let mut next_x = Some(start);
|
||||
while pop.len() < size {
|
||||
let x = next_x.take().unwrap_or_else(|| {
|
||||
let mut x: Vec<f64> = se.lo.iter().zip(&se.hi).map(|(l, h)| rng.between(*l, *h)).collect();
|
||||
se.repair(&mut x);
|
||||
x
|
||||
});
|
||||
match se.score(x) {
|
||||
Some(m) => pop.push(m),
|
||||
None => break,
|
||||
}
|
||||
}
|
||||
let mut best: Option<Member> = None;
|
||||
for m in &pop {
|
||||
if best.as_ref().is_none_or(|b| m.better(b, tol)) {
|
||||
best = Some(m.clone());
|
||||
}
|
||||
}
|
||||
let mut last_improved = 0u64;
|
||||
let mut generation = 0u64;
|
||||
'search: while pop.len() == size {
|
||||
let Some(b) = best.clone() else { break };
|
||||
// Converged: 99% of the population is feasible and within Convergence of the best.
|
||||
let close = pop.iter().filter(|m| m.viol <= tol && (m.phi - b.phi).abs() <= o.convergence * b.phi.abs().max(1.0)).count();
|
||||
if b.viol <= tol && generation >= 5 && close * 100 >= size * 99 {
|
||||
end = EvoEnd::Converged;
|
||||
break;
|
||||
}
|
||||
if generation.saturating_sub(last_improved) >= o.stall {
|
||||
end = EvoEnd::Stalled;
|
||||
break;
|
||||
}
|
||||
// The next generation: the best survives, the rest are children of two tournaments.
|
||||
let mut next = Vec::with_capacity(size);
|
||||
next.push(b.clone());
|
||||
while next.len() < size {
|
||||
let mut parent = || {
|
||||
let (i, k) = (rng.below(size), rng.below(size));
|
||||
match (pop.get(i), pop.get(k)) {
|
||||
(Some(p), Some(q)) => Some(if q.better(p, tol) { q.x.clone() } else { p.x.clone() }),
|
||||
_ => None,
|
||||
}
|
||||
};
|
||||
let (Some(pa), Some(pb)) = (parent(), parent()) else { break 'search };
|
||||
let mut x: Vec<f64> = pa
|
||||
.iter()
|
||||
.zip(&pb)
|
||||
.zip(&span)
|
||||
.map(|((a, c), w)| {
|
||||
let mut v = a + rng.between(-0.25, 1.25) * (c - a);
|
||||
if rng.uniform() < o.mutation {
|
||||
v += 0.1 * w * rng.normal();
|
||||
}
|
||||
v
|
||||
})
|
||||
.collect();
|
||||
se.repair(&mut x);
|
||||
match se.score(x) {
|
||||
Some(m) => next.push(m),
|
||||
None => break 'search,
|
||||
}
|
||||
}
|
||||
generation += 1;
|
||||
for m in &next {
|
||||
let Some(cur) = &best else { continue };
|
||||
if m.better(cur, tol) {
|
||||
let improved = cur.viol > tol || cur.phi - m.phi > 1e-9 * cur.phi.abs().max(1.0);
|
||||
best = Some(m.clone());
|
||||
if improved {
|
||||
last_improved = generation;
|
||||
}
|
||||
}
|
||||
}
|
||||
pop = next;
|
||||
}
|
||||
// Polish: a compass search on the continuous variables around the best member.
|
||||
if let Some(mut b) = best.clone() {
|
||||
let mut step: Vec<f64> = span.iter().map(|w| 0.01 * w).collect();
|
||||
let n = b.x.len();
|
||||
'polish: for _ in 0..10_000 {
|
||||
if step.iter().zip(&span).zip(integer).all(|((s, w), int)| *int || *s <= 1e-9 * w.max(1.0)) {
|
||||
break;
|
||||
}
|
||||
let mut moved = false;
|
||||
for j in 0..n {
|
||||
if integer.get(j).copied().unwrap_or(false) {
|
||||
continue;
|
||||
}
|
||||
for dir in [1.0, -1.0] {
|
||||
let mut x = b.x.clone();
|
||||
if let (Some(v), Some(s)) = (x.get_mut(j), step.get(j)) {
|
||||
*v += dir * s;
|
||||
}
|
||||
se.repair(&mut x);
|
||||
if x == b.x {
|
||||
continue;
|
||||
}
|
||||
let Some(m) = se.score(x) else { break 'polish };
|
||||
if m.better(&b, tol) {
|
||||
b = m;
|
||||
moved = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
if !moved {
|
||||
for s in &mut step {
|
||||
*s *= 0.5;
|
||||
}
|
||||
}
|
||||
}
|
||||
best = Some(b);
|
||||
}
|
||||
match best {
|
||||
Some(b) if b.viol <= tol => (Outcome::Solved(b.x, b.phi, end != EvoEnd::Budget), end),
|
||||
_ if end == EvoEnd::Budget => (Outcome::Limit, end),
|
||||
_ => (Outcome::Infeasible, end),
|
||||
}
|
||||
}
|
||||
@@ -0,0 +1,447 @@
|
||||
use super::problem::Evaluator;
|
||||
use super::*;
|
||||
|
||||
// ---------------------------------------------------------------- evaluation
|
||||
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub(super) enum Sense {
|
||||
Le,
|
||||
Ge,
|
||||
Eq,
|
||||
}
|
||||
|
||||
/// A scalar constraint `lhs (sense) rhs`.
|
||||
#[derive(Clone, Debug)]
|
||||
pub(super) struct Row {
|
||||
pub(super) lhs: CellRef,
|
||||
pub(super) sense: Sense,
|
||||
pub(super) rhs: Option<CellRef>,
|
||||
pub(super) rhs_const: f64,
|
||||
}
|
||||
|
||||
/// What a relaxation (or the whole solve) ended with.
|
||||
#[derive(Clone, Debug)]
|
||||
pub(super) enum Outcome {
|
||||
/// x, objective to minimise, and whether the solver converged (vs. ran out of budget).
|
||||
Solved(Vec<f64>, f64, bool),
|
||||
Infeasible,
|
||||
Unbounded,
|
||||
Limit,
|
||||
}
|
||||
|
||||
/// Depth-first branch and bound over a relaxation solver; returns the best integer point.
|
||||
pub(super) fn branch_and_bound(
|
||||
lo: &[f64],
|
||||
hi: &[f64],
|
||||
integer: &[bool],
|
||||
max_nodes: usize,
|
||||
gap: f64,
|
||||
relax: &mut dyn FnMut(&[f64], &[f64]) -> Outcome,
|
||||
) -> (Outcome, usize) {
|
||||
let mut stack: Vec<(Vec<f64>, Vec<f64>)> = vec![(lo.to_vec(), hi.to_vec())];
|
||||
let mut best: Option<(Vec<f64>, f64)> = None;
|
||||
let mut nodes = 0usize;
|
||||
let mut limited = false;
|
||||
let mut root = true;
|
||||
while let Some((l, h)) = stack.pop() {
|
||||
if nodes >= max_nodes {
|
||||
limited = true;
|
||||
break;
|
||||
}
|
||||
nodes += 1;
|
||||
let out = relax(&l, &h);
|
||||
let was_root = std::mem::replace(&mut root, false);
|
||||
let (x, obj, converged) = match out {
|
||||
Outcome::Solved(x, obj, c) => (x, obj, c),
|
||||
Outcome::Unbounded if was_root => return (Outcome::Unbounded, nodes),
|
||||
Outcome::Limit => {
|
||||
limited = true;
|
||||
break;
|
||||
}
|
||||
_ => continue,
|
||||
};
|
||||
if !converged {
|
||||
limited = true;
|
||||
}
|
||||
if let Some((_, b)) = &best
|
||||
&& obj >= b - gap * b.abs().max(1e-12) - 1e-9 * (1.0 + b.abs())
|
||||
{
|
||||
continue;
|
||||
}
|
||||
// The most fractional integer variable.
|
||||
let frac = x
|
||||
.iter()
|
||||
.zip(integer)
|
||||
.enumerate()
|
||||
.filter(|(_, (_, int))| **int)
|
||||
.map(|(j, (v, _))| (j, (v - v.round()).abs()))
|
||||
.filter(|(_, f)| *f > 1e-6)
|
||||
.max_by(|a, b| a.1.total_cmp(&b.1));
|
||||
match frac {
|
||||
None => {
|
||||
let snapped: Vec<f64> = x.iter().zip(integer).map(|(v, int)| if *int { v.round() } else { *v }).collect();
|
||||
best = Some((snapped, obj));
|
||||
}
|
||||
Some((j, _)) => {
|
||||
let v = x.get(j).copied().unwrap_or(0.0);
|
||||
let (mut down_h, mut up_l) = (h.clone(), l.clone());
|
||||
if let Some(e) = down_h.get_mut(j) {
|
||||
*e = e.min(v.floor());
|
||||
}
|
||||
if let Some(e) = up_l.get_mut(j) {
|
||||
*e = e.max(v.ceil());
|
||||
}
|
||||
// Explore the nearer side first (pushed last).
|
||||
if v - v.floor() < 0.5 {
|
||||
stack.push((up_l, h));
|
||||
stack.push((l, down_h));
|
||||
} else {
|
||||
stack.push((l, down_h));
|
||||
stack.push((up_l, h));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
match best {
|
||||
Some((x, obj)) => (Outcome::Solved(x, obj, !limited), nodes),
|
||||
None if limited => (Outcome::Limit, nodes),
|
||||
None => (Outcome::Infeasible, nodes),
|
||||
}
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------- simplex
|
||||
|
||||
/// Pivots the tableau on (r, j), updating the objective row.
|
||||
fn pivot(t: &mut [Vec<f64>], obj: &mut [f64], r: usize, j: usize) {
|
||||
let Some(prow) = t.get_mut(r) else { return };
|
||||
let pv = prow.get(j).copied().unwrap_or(0.0);
|
||||
if pv == 0.0 {
|
||||
return;
|
||||
}
|
||||
for v in prow.iter_mut() {
|
||||
*v /= pv;
|
||||
}
|
||||
let prow = prow.clone();
|
||||
for (i, row) in t.iter_mut().enumerate() {
|
||||
if i == r {
|
||||
continue;
|
||||
}
|
||||
let f = row.get(j).copied().unwrap_or(0.0);
|
||||
if f != 0.0 {
|
||||
for (v, p) in row.iter_mut().zip(&prow) {
|
||||
*v -= f * p;
|
||||
}
|
||||
}
|
||||
}
|
||||
let f = obj.get(j).copied().unwrap_or(0.0);
|
||||
if f != 0.0 {
|
||||
for (v, p) in obj.iter_mut().zip(&prow) {
|
||||
*v -= f * p;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
enum Step {
|
||||
Optimal,
|
||||
Unbounded,
|
||||
Limit,
|
||||
}
|
||||
|
||||
/// Minimises with the reduced costs in `obj` (last entry: −objective). Columns not `allowed`
|
||||
/// never enter.
|
||||
fn simplex(t: &mut [Vec<f64>], basis: &mut [usize], obj: &mut [f64], allowed: &[bool], budget: &mut u64) -> Step {
|
||||
let n = obj.len().saturating_sub(1);
|
||||
let cost = (t.len() as u64 + 1).saturating_mul(n as u64 + 1).max(1);
|
||||
let mut stalled = 0u32;
|
||||
loop {
|
||||
if *budget < cost {
|
||||
return Step::Limit;
|
||||
}
|
||||
*budget -= cost;
|
||||
let bland = stalled > 50;
|
||||
let mut enter = None;
|
||||
let mut best = -EPS;
|
||||
for (j, d) in obj.iter().take(n).enumerate() {
|
||||
if !allowed.get(j).copied().unwrap_or(false) || *d >= -EPS {
|
||||
continue;
|
||||
}
|
||||
if bland {
|
||||
enter = Some(j);
|
||||
break;
|
||||
}
|
||||
if *d < best {
|
||||
best = *d;
|
||||
enter = Some(j);
|
||||
}
|
||||
}
|
||||
let Some(j) = enter else { return Step::Optimal };
|
||||
let mut leave: Option<(usize, f64)> = None;
|
||||
for (i, row) in t.iter().enumerate() {
|
||||
let a = row.get(j).copied().unwrap_or(0.0);
|
||||
if a <= EPS {
|
||||
continue;
|
||||
}
|
||||
let ratio = row.get(n).copied().unwrap_or(0.0) / a;
|
||||
let better = match leave {
|
||||
None => true,
|
||||
Some((li, lr)) => ratio < lr - 1e-12 || (ratio <= lr + 1e-12 && basis.get(i) < basis.get(li)),
|
||||
};
|
||||
if better {
|
||||
leave = Some((i, ratio));
|
||||
}
|
||||
}
|
||||
let Some((r, ratio)) = leave else { return Step::Unbounded };
|
||||
stalled = if ratio <= EPS { stalled + 1 } else { 0 };
|
||||
pivot(t, obj, r, j);
|
||||
if let Some(b) = basis.get_mut(r) {
|
||||
*b = j;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Minimises `c·x` subject to `rows` (coefficients, sense, rhs) and `lo ≤ x ≤ hi`.
|
||||
pub(super) fn solve_lp(c: &[f64], rows: &[(Vec<f64>, Sense, f64)], lo: &[f64], hi: &[f64], budget: &mut u64) -> Outcome {
|
||||
let n = c.len();
|
||||
if lo.iter().zip(hi).any(|(l, h)| l > h) {
|
||||
return Outcome::Infeasible;
|
||||
}
|
||||
// x_j = offset + Σ sign·y_k with y ≥ 0.
|
||||
struct Map {
|
||||
var: usize,
|
||||
sign: f64,
|
||||
}
|
||||
let mut cols: Vec<Map> = Vec::new();
|
||||
let mut offset = vec![0.0; n];
|
||||
let mut extra: Vec<(Vec<(usize, f64)>, Sense, f64)> = Vec::new();
|
||||
for j in 0..n {
|
||||
let (l, h) = (lo.get(j).copied().unwrap_or(f64::NEG_INFINITY), hi.get(j).copied().unwrap_or(f64::INFINITY));
|
||||
if l.is_finite() {
|
||||
if let Some(o) = offset.get_mut(j) {
|
||||
*o = l;
|
||||
}
|
||||
cols.push(Map { var: j, sign: 1.0 });
|
||||
if h.is_finite() {
|
||||
extra.push((vec![(cols.len() - 1, 1.0)], Sense::Le, h - l));
|
||||
}
|
||||
} else if h.is_finite() {
|
||||
if let Some(o) = offset.get_mut(j) {
|
||||
*o = h;
|
||||
}
|
||||
cols.push(Map { var: j, sign: -1.0 });
|
||||
} else {
|
||||
cols.push(Map { var: j, sign: 1.0 });
|
||||
cols.push(Map { var: j, sign: -1.0 });
|
||||
}
|
||||
}
|
||||
let ny = cols.len();
|
||||
// Rows over y: (coefficients, sense, rhs).
|
||||
let mut yrows: Vec<(Vec<f64>, Sense, f64)> = Vec::with_capacity(rows.len() + extra.len());
|
||||
for (a, sense, b) in rows {
|
||||
let mut coef = vec![0.0; ny];
|
||||
for (k, m) in cols.iter().enumerate() {
|
||||
if let Some(v) = coef.get_mut(k) {
|
||||
*v = a.get(m.var).copied().unwrap_or(0.0) * m.sign;
|
||||
}
|
||||
}
|
||||
let shift: f64 = a.iter().zip(&offset).map(|(x, o)| x * o).sum();
|
||||
yrows.push((coef, *sense, b - shift));
|
||||
}
|
||||
for (entries, sense, b) in extra {
|
||||
let mut coef = vec![0.0; ny];
|
||||
for (k, v) in entries {
|
||||
if let Some(e) = coef.get_mut(k) {
|
||||
*e = v;
|
||||
}
|
||||
}
|
||||
yrows.push((coef, sense, b));
|
||||
}
|
||||
let m = yrows.len();
|
||||
let n_slack = yrows.iter().filter(|r| r.1 != Sense::Eq).count();
|
||||
// Rows with b < 0 are negated so every right side is ≥ 0.
|
||||
let normalized: Vec<(Vec<f64>, Sense, f64)> = yrows
|
||||
.into_iter()
|
||||
.map(|(a, s, b)| {
|
||||
if b < 0.0 {
|
||||
let s = match s {
|
||||
Sense::Le => Sense::Ge,
|
||||
Sense::Ge => Sense::Le,
|
||||
Sense::Eq => Sense::Eq,
|
||||
};
|
||||
(a.iter().map(|v| -v).collect(), s, -b)
|
||||
} else {
|
||||
(a, s, b)
|
||||
}
|
||||
})
|
||||
.collect();
|
||||
let n_art = normalized.iter().filter(|r| r.1 != Sense::Le).count();
|
||||
let width = ny + n_slack + n_art;
|
||||
let mut t: Vec<Vec<f64>> = Vec::with_capacity(m);
|
||||
let mut basis = vec![0usize; m];
|
||||
let (mut si, mut ai) = (ny, ny + n_slack);
|
||||
for (i, (a, s, b)) in normalized.iter().enumerate() {
|
||||
let mut row = vec![0.0; width + 1];
|
||||
for (v, x) in row.iter_mut().zip(a) {
|
||||
*v = *x;
|
||||
}
|
||||
match s {
|
||||
Sense::Le => {
|
||||
if let Some(v) = row.get_mut(si) {
|
||||
*v = 1.0;
|
||||
}
|
||||
if let Some(bi) = basis.get_mut(i) {
|
||||
*bi = si;
|
||||
}
|
||||
si += 1;
|
||||
}
|
||||
Sense::Ge => {
|
||||
if let Some(v) = row.get_mut(si) {
|
||||
*v = -1.0;
|
||||
}
|
||||
si += 1;
|
||||
if let Some(v) = row.get_mut(ai) {
|
||||
*v = 1.0;
|
||||
}
|
||||
if let Some(bi) = basis.get_mut(i) {
|
||||
*bi = ai;
|
||||
}
|
||||
ai += 1;
|
||||
}
|
||||
Sense::Eq => {
|
||||
if let Some(v) = row.get_mut(ai) {
|
||||
*v = 1.0;
|
||||
}
|
||||
if let Some(bi) = basis.get_mut(i) {
|
||||
*bi = ai;
|
||||
}
|
||||
ai += 1;
|
||||
}
|
||||
}
|
||||
if let Some(v) = row.get_mut(width) {
|
||||
*v = *b;
|
||||
}
|
||||
t.push(row);
|
||||
}
|
||||
let is_art = |j: usize| j >= ny + n_slack && j < width;
|
||||
// Phase 1: minimise the sum of the artificials.
|
||||
let mut obj = vec![0.0; width + 1];
|
||||
for j in (ny + n_slack)..width {
|
||||
if let Some(v) = obj.get_mut(j) {
|
||||
*v = 1.0;
|
||||
}
|
||||
}
|
||||
for (i, row) in t.iter().enumerate() {
|
||||
if basis.get(i).is_some_and(|b| is_art(*b)) {
|
||||
for (v, x) in obj.iter_mut().zip(row) {
|
||||
*v -= x;
|
||||
}
|
||||
}
|
||||
}
|
||||
let all: Vec<bool> = vec![true; width];
|
||||
if n_art > 0 {
|
||||
match simplex(&mut t, &mut basis, &mut obj, &all, budget) {
|
||||
Step::Limit => return Outcome::Limit,
|
||||
Step::Unbounded | Step::Optimal => {}
|
||||
}
|
||||
let scale = normalized.iter().map(|r| r.2.abs()).fold(1.0f64, f64::max);
|
||||
if -obj.get(width).copied().unwrap_or(0.0) > 1e-7 * scale {
|
||||
return Outcome::Infeasible;
|
||||
}
|
||||
// Drive zero-valued artificials out of the basis where possible.
|
||||
for i in 0..m {
|
||||
if !basis.get(i).is_some_and(|b| is_art(*b)) {
|
||||
continue;
|
||||
}
|
||||
let j = t.get(i).and_then(|row| (0..ny + n_slack).find(|&j| row.get(j).is_some_and(|v| v.abs() > EPS)));
|
||||
if let Some(j) = j {
|
||||
pivot(&mut t, &mut obj, i, j);
|
||||
if let Some(b) = basis.get_mut(i) {
|
||||
*b = j;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Phase 2: the real costs; artificials may not re-enter.
|
||||
let mut cost = vec![0.0; width + 1];
|
||||
for (k, mp) in cols.iter().enumerate() {
|
||||
if let Some(v) = cost.get_mut(k) {
|
||||
*v = c.get(mp.var).copied().unwrap_or(0.0) * mp.sign;
|
||||
}
|
||||
}
|
||||
let mut obj = cost.clone();
|
||||
for (i, row) in t.iter().enumerate() {
|
||||
let cb = basis.get(i).and_then(|b| cost.get(*b)).copied().unwrap_or(0.0);
|
||||
if cb != 0.0 {
|
||||
for (v, x) in obj.iter_mut().zip(row) {
|
||||
*v -= cb * x;
|
||||
}
|
||||
}
|
||||
}
|
||||
let allowed: Vec<bool> = (0..width).map(|j| !is_art(j)).collect();
|
||||
match simplex(&mut t, &mut basis, &mut obj, &allowed, budget) {
|
||||
Step::Limit => return Outcome::Limit,
|
||||
Step::Unbounded => return Outcome::Unbounded,
|
||||
Step::Optimal => {}
|
||||
}
|
||||
let mut y = vec![0.0; width];
|
||||
for (i, row) in t.iter().enumerate() {
|
||||
if let (Some(b), Some(v)) = (basis.get(i), row.get(width))
|
||||
&& let Some(e) = y.get_mut(*b)
|
||||
{
|
||||
*e = *v;
|
||||
}
|
||||
}
|
||||
let mut x = offset;
|
||||
for (k, mp) in cols.iter().enumerate() {
|
||||
if let Some(e) = x.get_mut(mp.var) {
|
||||
*e += mp.sign * y.get(k).copied().unwrap_or(0.0);
|
||||
}
|
||||
}
|
||||
let val: f64 = c.iter().zip(&x).map(|(a, b)| a * b).sum();
|
||||
Outcome::Solved(x, val, true)
|
||||
}
|
||||
|
||||
/// Linearises the model at `x0` (n + 1 evaluations) and checks it at one more point. Returns the
|
||||
/// objective gradient and constant, and per row the gradient and constant of `lhs − rhs`.
|
||||
#[allow(clippy::type_complexity)]
|
||||
pub(super) fn linearize(ev: &mut Evaluator, x0: &[f64]) -> std::result::Result<((Vec<f64>, f64), Vec<(Vec<f64>, f64)>), &'static str> {
|
||||
let n = x0.len();
|
||||
let p0 = ev.eval(x0).ok_or("budget")?;
|
||||
if !p0.f.is_finite() || p0.g.iter().any(|v| !v.is_finite()) {
|
||||
return Err("error");
|
||||
}
|
||||
let mut cf = vec![0.0; n];
|
||||
let mut ga = vec![vec![0.0; n]; p0.g.len()];
|
||||
for j in 0..n {
|
||||
let mut x = x0.to_vec();
|
||||
if let Some(v) = x.get_mut(j) {
|
||||
*v += 1.0;
|
||||
}
|
||||
let p = ev.eval(&x).ok_or("budget")?;
|
||||
if let Some(v) = cf.get_mut(j) {
|
||||
*v = p.f - p0.f;
|
||||
}
|
||||
for (row, (gi, g0)) in ga.iter_mut().zip(p.g.iter().zip(&p0.g)) {
|
||||
if let Some(v) = row.get_mut(j) {
|
||||
*v = gi - g0;
|
||||
}
|
||||
}
|
||||
}
|
||||
let dot = |a: &[f64], x: &[f64]| a.iter().zip(x).map(|(p, q)| p * q).sum::<f64>();
|
||||
let f_const = p0.f - dot(&cf, x0);
|
||||
let rows: Vec<(Vec<f64>, f64)> = ga.into_iter().zip(&p0.g).map(|(a, g0)| (a.clone(), g0 - dot(&a, x0))).collect();
|
||||
// Linearity check at a spread-out point.
|
||||
let xt: Vec<f64> =
|
||||
x0.iter().enumerate().map(|(j, v)| v + if j % 2 == 0 { 1.7 + 0.3 * (j % 5) as f64 } else { -0.6 - 0.2 * (j % 3) as f64 }).collect();
|
||||
let pt = ev.eval(&xt).ok_or("budget")?;
|
||||
let close = |actual: f64, pred: f64, terms: f64| actual.is_finite() && (actual - pred).abs() <= 1e-6 * (1.0 + actual.abs() + terms);
|
||||
let terms = |a: &[f64]| a.iter().zip(&xt).map(|(p, q)| (p * q).abs()).sum::<f64>();
|
||||
if !close(pt.f, f_const + dot(&cf, &xt), terms(&cf)) {
|
||||
return Err("nonlinear");
|
||||
}
|
||||
for ((a, k), g) in rows.iter().zip(&pt.g) {
|
||||
if !close(*g, k + dot(a, &xt), terms(a)) {
|
||||
return Err("nonlinear");
|
||||
}
|
||||
}
|
||||
Ok(((cf, f_const), rows))
|
||||
}
|
||||
@@ -0,0 +1,459 @@
|
||||
use super::*;
|
||||
|
||||
// ---------------------------------------------------------------- model
|
||||
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub(super) enum Rel {
|
||||
Le,
|
||||
Eq,
|
||||
Ge,
|
||||
Int,
|
||||
Bin,
|
||||
}
|
||||
|
||||
impl Rel {
|
||||
fn parse(s: &str) -> Option<Rel> {
|
||||
match s.trim().to_ascii_lowercase().as_str() {
|
||||
"<=" | "≤" | "=<" | "le" => Some(Rel::Le),
|
||||
">=" | "≥" | "=>" | "ge" => Some(Rel::Ge),
|
||||
"=" | "==" | "eq" => Some(Rel::Eq),
|
||||
"int" | "integer" => Some(Rel::Int),
|
||||
"bin" | "binary" => Some(Rel::Bin),
|
||||
_ => None,
|
||||
}
|
||||
}
|
||||
pub(super) fn text(self) -> &'static str {
|
||||
match self {
|
||||
Rel::Le => "<=",
|
||||
Rel::Eq => "=",
|
||||
Rel::Ge => ">=",
|
||||
Rel::Int => "int",
|
||||
Rel::Bin => "bin",
|
||||
}
|
||||
}
|
||||
/// The relation code of the saved model (`solver_relN`), as in Solver's VBA `SolverAdd`.
|
||||
fn code(self) -> u8 {
|
||||
match self {
|
||||
Rel::Le => 1,
|
||||
Rel::Eq => 2,
|
||||
Rel::Ge => 3,
|
||||
Rel::Int => 4,
|
||||
Rel::Bin => 5,
|
||||
}
|
||||
}
|
||||
fn from_code(c: &str) -> Option<Rel> {
|
||||
match c.trim() {
|
||||
"1" => Some(Rel::Le),
|
||||
"2" => Some(Rel::Eq),
|
||||
"3" => Some(Rel::Ge),
|
||||
"4" => Some(Rel::Int),
|
||||
"5" => Some(Rel::Bin),
|
||||
_ => None,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub(super) enum Rhs {
|
||||
None,
|
||||
Const(f64),
|
||||
Cells(Vec<CellRef>),
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub(super) struct Constraint {
|
||||
pub(super) lhs: Vec<CellRef>,
|
||||
lhs_text: String,
|
||||
pub(super) rel: Rel,
|
||||
pub(super) rhs: Rhs,
|
||||
rhs_text: String,
|
||||
}
|
||||
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub(super) enum Goal {
|
||||
Max,
|
||||
Min,
|
||||
Value(f64),
|
||||
}
|
||||
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub(super) enum Method {
|
||||
Simplex,
|
||||
Grg,
|
||||
Evolutionary,
|
||||
}
|
||||
|
||||
/// The Evolutionary method's options (Solver Options › Evolutionary). The search keeps no
|
||||
/// clock, so the time without improvement is counted in generations.
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub(super) struct EvoOptions {
|
||||
pub(super) population: usize,
|
||||
pub(super) mutation: f64,
|
||||
pub(super) seed: u64,
|
||||
pub(super) convergence: f64,
|
||||
pub(super) stall: u64,
|
||||
pub(super) require_bounds: bool,
|
||||
}
|
||||
|
||||
impl Default for EvoOptions {
|
||||
fn default() -> EvoOptions {
|
||||
EvoOptions { population: 100, mutation: 0.075, seed: 0, convergence: 1e-4, stall: 50, require_bounds: true }
|
||||
}
|
||||
}
|
||||
|
||||
fn evo_options(p: &Json) -> EvoOptions {
|
||||
let d = EvoOptions::default();
|
||||
let num = |k: &str| p.get(k).and_then(Json::as_f64).or_else(|| str_param(p, k).and_then(num_text)).filter(|v| v.is_finite());
|
||||
EvoOptions {
|
||||
population: num("populationSize").map_or(d.population, |v| v.clamp(MIN_POPULATION as f64, MAX_POPULATION as f64) as usize),
|
||||
mutation: num("mutationRate").map_or(d.mutation, |v| v.clamp(0.0, 1.0)),
|
||||
seed: num("randomSeed").map_or(d.seed, |v| v.clamp(0.0, f64::from(u32::MAX)) as u64),
|
||||
convergence: num("convergence").map_or(d.convergence, |v| v.clamp(0.0, 1.0)),
|
||||
stall: num("generationsWithoutImprovement").map_or(d.stall, |v| v.clamp(1.0, 100_000.0) as u64),
|
||||
require_bounds: bool_param(p, "requireBounds").unwrap_or(d.require_bounds),
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub(super) struct Model {
|
||||
pub(super) sheet: usize,
|
||||
pub(super) objective: CellRef,
|
||||
pub(super) goal: Goal,
|
||||
pub(super) vars: Vec<CellRef>,
|
||||
vars_text: String,
|
||||
pub(super) cons: Vec<Constraint>,
|
||||
pub(super) non_negative: bool,
|
||||
pub(super) method: Method,
|
||||
pub(super) evo: EvoOptions,
|
||||
}
|
||||
|
||||
/// Cells of `"B2:B4,D1"` on one sheet (row-major per area, in order), at most `cap`.
|
||||
fn cells_of(wb: &Workbook, t: &str, default_sheet: usize, cap: usize) -> Result<(usize, Vec<CellRef>)> {
|
||||
let mut sheet = None;
|
||||
let mut out = Vec::new();
|
||||
for part in t.split(',').map(str::trim).filter(|x| !x.is_empty()) {
|
||||
let (sh, body) = split_sheet(part);
|
||||
let si = match sh {
|
||||
Some(n) => wb.sheet_index(&n).ok_or_else(|| bad(CMD, format!("no sheet named `{n}`")))?,
|
||||
None => default_sheet,
|
||||
};
|
||||
if *sheet.get_or_insert(si) != si {
|
||||
return Err(bad(CMD, "all cells of a reference must be on one sheet"));
|
||||
}
|
||||
let r = RangeRef::parse(body.trim()).ok_or_else(|| bad(CMD, format!("not a reference: `{part}`")))?;
|
||||
if r.count() > cap as u64 || out.len() + r.count() as usize > cap {
|
||||
return Err(bad(CMD, format!("too many cells (at most {cap})")));
|
||||
}
|
||||
out.extend(r.iter());
|
||||
}
|
||||
if out.is_empty() {
|
||||
return Err(bad(CMD, format!("not a reference: `{t}`")));
|
||||
}
|
||||
Ok((sheet.unwrap_or(default_sheet), out))
|
||||
}
|
||||
|
||||
fn num_text(t: &str) -> Option<f64> {
|
||||
gridcraft_core::parse::parse_number_text(t.trim()).filter(|v| v.is_finite())
|
||||
}
|
||||
|
||||
/// Splits `"B5 <= 10"` / `"B2 int"` into (lhs, relation, rhs).
|
||||
fn split_constraint(t: &str) -> Option<(String, Rel, String)> {
|
||||
for op in ["<=", ">=", "=<", "=>", "≤", "≥", "=="] {
|
||||
if let Some((a, b)) = t.split_once(op) {
|
||||
return Some((a.trim().into(), Rel::parse(op)?, b.trim().into()));
|
||||
}
|
||||
}
|
||||
if let Some((a, b)) = t.split_once('=') {
|
||||
return Some((a.trim().into(), Rel::Eq, b.trim().into()));
|
||||
}
|
||||
let lower = t.trim().to_ascii_lowercase();
|
||||
for (suffix, rel) in [("integer", Rel::Int), ("binary", Rel::Bin), ("int", Rel::Int), ("bin", Rel::Bin)] {
|
||||
if let Some(a) = lower.strip_suffix(suffix) {
|
||||
let a = a.trim();
|
||||
return Some((t.trim().get(..a.len())?.trim().to_string(), rel, String::new()));
|
||||
}
|
||||
}
|
||||
None
|
||||
}
|
||||
|
||||
fn parse_constraint(wb: &Workbook, sheet: usize, c: &Json) -> Result<Constraint> {
|
||||
let (lhs_text, rel, rhs_json) = match c {
|
||||
Json::String(t) => {
|
||||
let (l, r, rhs) = split_constraint(t).ok_or_else(|| bad(CMD, format!("not a constraint: `{t}`")))?;
|
||||
(l, r, Json::String(rhs))
|
||||
}
|
||||
Json::Object(_) => {
|
||||
let lhs = str_param(c, "lhs").or_else(|| str_param(c, "cell")).or_else(|| str_param(c, "range")).unwrap_or("").to_string();
|
||||
let op = str_param(c, "op").or_else(|| str_param(c, "relation")).unwrap_or("");
|
||||
let rel = Rel::parse(op).ok_or_else(|| bad(CMD, format!("unknown constraint relation `{op}`")))?;
|
||||
(lhs, rel, c.get("rhs").or_else(|| c.get("value")).cloned().unwrap_or(Json::Null))
|
||||
}
|
||||
_ => return Err(bad(CMD, "a constraint is a string like \"B5 <= 10\" or {lhs, op, rhs}")),
|
||||
};
|
||||
let (lsi, lhs) = cells_of(wb, &lhs_text, sheet, MAX_CONSTRAINTS)?;
|
||||
if lsi != sheet {
|
||||
return Err(EngineError::Other("Constraint cells must be on the objective cell's sheet.".into()));
|
||||
}
|
||||
let (rhs, rhs_text) = match rel {
|
||||
Rel::Int | Rel::Bin => (Rhs::None, String::new()),
|
||||
_ => match &rhs_json {
|
||||
Json::Number(n) => {
|
||||
let v = n.as_f64().filter(|v| v.is_finite()).ok_or_else(|| bad(CMD, "the constraint's right side must be a number"))?;
|
||||
(Rhs::Const(v), v.to_string())
|
||||
}
|
||||
Json::String(t) => {
|
||||
let t = t.trim().trim_start_matches('=').trim();
|
||||
if let Some(v) = num_text(t) {
|
||||
(Rhs::Const(v), t.to_string())
|
||||
} else {
|
||||
let (rsi, cells) = cells_of(wb, t, sheet, MAX_CONSTRAINTS)?;
|
||||
if rsi != sheet {
|
||||
return Err(EngineError::Other("Constraint cells must be on the objective cell's sheet.".into()));
|
||||
}
|
||||
if cells.len() != 1 && cells.len() != lhs.len() {
|
||||
return Err(EngineError::Other("The constraint's two sides must be the same size (or the right side one cell).".into()));
|
||||
}
|
||||
(Rhs::Cells(cells), t.to_string())
|
||||
}
|
||||
}
|
||||
_ => return Err(bad(CMD, "missing the constraint's right side")),
|
||||
},
|
||||
};
|
||||
Ok(Constraint { lhs, lhs_text, rel, rhs, rhs_text })
|
||||
}
|
||||
|
||||
pub(super) fn parse_model(s: &Session, p: &Json) -> Result<Model> {
|
||||
let d = s.doc()?;
|
||||
let wb = &d.wb;
|
||||
let obj_text = str_param(p, "objective").or_else(|| str_param(p, "setObjective")).or_else(|| str_param(p, "set")).unwrap_or("");
|
||||
let (sheet, obj) = cells_of(wb, obj_text, wb.active_sheet(), 1).map_err(|_| EngineError::Other("The objective must be a single cell.".into()))?;
|
||||
let objective = *obj.first().ok_or_else(|| EngineError::Other("The objective must be a single cell.".into()))?;
|
||||
let sh = wb.sheet(sheet).ok_or(EngineError::NoDocument)?;
|
||||
if sh.cell(objective).is_none_or(|c| c.formula.is_none()) {
|
||||
return Err(EngineError::Other("The objective cell must contain a formula.".into()));
|
||||
}
|
||||
let goal = match str_param(p, "goal").or_else(|| str_param(p, "to")).unwrap_or("max").to_ascii_lowercase().as_str() {
|
||||
"max" | "maximize" | "maximum" => Goal::Max,
|
||||
"min" | "minimize" | "minimum" => Goal::Min,
|
||||
"value" | "valueof" | "value of" => {
|
||||
let v = f64_param(p, "value")
|
||||
.or_else(|| str_param(p, "value").and_then(num_text))
|
||||
.filter(|v| v.is_finite())
|
||||
.ok_or_else(|| bad(CMD, "goal `value` needs a number in `value`"))?;
|
||||
Goal::Value(v)
|
||||
}
|
||||
g => return Err(bad(CMD, format!("unknown goal `{g}` (max, min or value)"))),
|
||||
};
|
||||
let vars_text = str_param(p, "variables").or_else(|| str_param(p, "changing")).or_else(|| str_param(p, "byChanging")).unwrap_or("").to_string();
|
||||
let (vsi, mut vars) = cells_of(wb, &vars_text, sheet, MAX_VARS).map_err(|e| match e {
|
||||
EngineError::BadParams { msg, .. } if msg.starts_with("too many") => {
|
||||
EngineError::Other(format!("Too many variable cells (at most {MAX_VARS})."))
|
||||
}
|
||||
_ => EngineError::Other("The variable cells must be a reference like B2:B4.".into()),
|
||||
})?;
|
||||
if vsi != sheet {
|
||||
return Err(EngineError::Other("The variable cells must be on the objective cell's sheet.".into()));
|
||||
}
|
||||
let mut seen = std::collections::HashSet::new();
|
||||
vars.retain(|c| seen.insert(*c));
|
||||
if vars.iter().any(|c| sh.cell(*c).is_some_and(|c| c.formula.is_some())) {
|
||||
return Err(EngineError::Other("Variable cells must hold constants, not formulas.".into()));
|
||||
}
|
||||
let mut cons = Vec::new();
|
||||
let mut scalar = 0usize;
|
||||
if let Some(list) = p.get("constraints") {
|
||||
let list = list.as_array().ok_or_else(|| bad(CMD, "`constraints` must be an array"))?;
|
||||
if list.len() > MAX_CONSTRAINTS {
|
||||
return Err(EngineError::Other(format!("Too many constraints (at most {MAX_CONSTRAINTS}).")));
|
||||
}
|
||||
for c in list {
|
||||
let c = parse_constraint(wb, sheet, c)?;
|
||||
scalar = scalar.saturating_add(c.lhs.len());
|
||||
if scalar > MAX_CONSTRAINTS {
|
||||
return Err(EngineError::Other(format!("Too many constraints (at most {MAX_CONSTRAINTS} cells).")));
|
||||
}
|
||||
if matches!(c.rel, Rel::Int | Rel::Bin) && c.lhs.iter().any(|x| !vars.contains(x)) {
|
||||
return Err(EngineError::Other("Integer and binary constraints must be on variable cells.".into()));
|
||||
}
|
||||
cons.push(c);
|
||||
}
|
||||
}
|
||||
let non_negative = bool_param(p, "nonNegative").unwrap_or(true);
|
||||
let method = match str_param(p, "method").unwrap_or("grg").to_ascii_lowercase().replace([' ', '_', '-'], "").as_str() {
|
||||
"grg" | "grgnonlinear" | "nonlinear" => Method::Grg,
|
||||
"simplex" | "simplexlp" | "lp" | "linear" => Method::Simplex,
|
||||
"evolutionary" | "evo" | "genetic" => Method::Evolutionary,
|
||||
m => return Err(bad(CMD, format!("unknown method `{m}` (grg, simplex or evolutionary)"))),
|
||||
};
|
||||
Ok(Model { sheet, objective, goal, vars, vars_text, cons, non_negative, method, evo: evo_options(p) })
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------- saved model
|
||||
|
||||
fn ref_text(sheet_name: &str, cells: &str) -> String {
|
||||
let q = crate::pivot::quote_sheet(sheet_name);
|
||||
cells
|
||||
.split(',')
|
||||
.map(str::trim)
|
||||
.filter(|x| !x.is_empty())
|
||||
.filter_map(|part| {
|
||||
let r = RangeRef::parse(split_sheet(part).1.trim())?;
|
||||
let abs = |c: CellRef| format!("${}${}", gridcraft_core::col_to_letters(c.col), c.row as u64 + 1);
|
||||
Some(if r.is_single() { format!("{q}!{}", abs(r.start)) } else { format!("{q}!{}:{}", abs(r.start), abs(r.end)) })
|
||||
})
|
||||
.collect::<Vec<_>>()
|
||||
.join(",")
|
||||
}
|
||||
|
||||
/// The model as hidden names scoped to its sheet (replacing any saved before).
|
||||
pub(super) fn save_model(wb: &mut Workbook, m: &Model) {
|
||||
let Some(sheet_name) = wb.sheet(m.sheet).map(|s| s.name.clone()) else { return };
|
||||
wb.names.retain(|n| !(n.scope == Some(m.sheet) && n.name.to_ascii_lowercase().starts_with("solver_")));
|
||||
let mut put = |name: String, formula: String| {
|
||||
wb.names.push(DefinedName { name, scope: Some(m.sheet), formula, comment: String::new(), hidden: true });
|
||||
};
|
||||
put("solver_opt".into(), ref_text(&sheet_name, &m.objective.a1()));
|
||||
let (typ, val) = match m.goal {
|
||||
Goal::Max => (1, 0.0),
|
||||
Goal::Min => (2, 0.0),
|
||||
Goal::Value(v) => (3, v),
|
||||
};
|
||||
put("solver_typ".into(), typ.to_string());
|
||||
put("solver_val".into(), val.to_string());
|
||||
put("solver_adj".into(), ref_text(&sheet_name, &m.vars_text));
|
||||
for (i, c) in m.cons.iter().enumerate() {
|
||||
let k = i + 1;
|
||||
put(format!("solver_lhs{k}"), ref_text(&sheet_name, &c.lhs_text));
|
||||
put(format!("solver_rel{k}"), c.rel.code().to_string());
|
||||
let rhs = match &c.rhs {
|
||||
Rhs::None => {
|
||||
if c.rel == Rel::Bin {
|
||||
"\"binary\"".to_string()
|
||||
} else {
|
||||
"\"integer\"".to_string()
|
||||
}
|
||||
}
|
||||
Rhs::Const(v) => v.to_string(),
|
||||
Rhs::Cells(_) => ref_text(&sheet_name, &c.rhs_text),
|
||||
};
|
||||
put(format!("solver_rhs{k}"), rhs);
|
||||
}
|
||||
put("solver_num".into(), m.cons.len().to_string());
|
||||
put("solver_neg".into(), if m.non_negative { "1" } else { "2" }.into());
|
||||
let eng = match m.method {
|
||||
Method::Grg => "1",
|
||||
Method::Simplex => "2",
|
||||
Method::Evolutionary => "3",
|
||||
};
|
||||
put("solver_eng".into(), eng.into());
|
||||
if m.method == Method::Evolutionary {
|
||||
put("solver_ssz".into(), m.evo.population.to_string());
|
||||
put("solver_mrt".into(), m.evo.mutation.to_string());
|
||||
put("solver_rsd".into(), m.evo.seed.to_string());
|
||||
put("solver_rbv".into(), if m.evo.require_bounds { "1" } else { "2" }.into());
|
||||
}
|
||||
put("solver_ver".into(), "3".into());
|
||||
}
|
||||
|
||||
/// The saved model of a sheet as `data.solver` params, if any.
|
||||
pub fn saved_model(wb: &Workbook, sheet: usize) -> Option<Json> {
|
||||
let get = |n: &str| {
|
||||
wb.names
|
||||
.iter()
|
||||
.find(|d| d.name.eq_ignore_ascii_case(n) && (d.scope == Some(sheet) || d.scope.is_none()))
|
||||
.map(|d| d.formula.trim().trim_start_matches('=').to_string())
|
||||
};
|
||||
let local = |t: String| -> String { t.split(',').map(|p| split_sheet(p.trim()).1.replace('$', "")).collect::<Vec<_>>().join(",") };
|
||||
let objective = local(get("solver_opt")?);
|
||||
if objective.is_empty() {
|
||||
return None;
|
||||
}
|
||||
let goal = match get("solver_typ").as_deref().map(str::trim) {
|
||||
Some("2") => "min",
|
||||
Some("3") => "value",
|
||||
_ => "max",
|
||||
};
|
||||
let value = get("solver_val").and_then(|v| num_text(&v)).unwrap_or(0.0);
|
||||
let count = get("solver_num").and_then(|v| num_text(&v)).unwrap_or(0.0).clamp(0.0, MAX_CONSTRAINTS as f64) as usize;
|
||||
let mut constraints = Vec::new();
|
||||
for k in 1..=count {
|
||||
let (Some(lhs), Some(rel)) = (get(&format!("solver_lhs{k}")), get(&format!("solver_rel{k}")).and_then(|r| Rel::from_code(&r))) else {
|
||||
continue;
|
||||
};
|
||||
let lhs = local(lhs);
|
||||
constraints.push(match rel {
|
||||
Rel::Int | Rel::Bin => json!(format!("{lhs} {}", rel.text())),
|
||||
_ => {
|
||||
let rhs = get(&format!("solver_rhs{k}")).unwrap_or_default();
|
||||
let rhs = if num_text(&rhs).is_some() { rhs } else { local(rhs) };
|
||||
json!(format!("{lhs} {} {rhs}", rel.text()))
|
||||
}
|
||||
});
|
||||
}
|
||||
let method = match get("solver_eng").as_deref().map(str::trim) {
|
||||
Some("2") => "simplex",
|
||||
Some("3") => "evolutionary",
|
||||
_ => "grg",
|
||||
};
|
||||
let mut model = json!({
|
||||
"objective": objective,
|
||||
"goal": goal,
|
||||
"value": value,
|
||||
"variables": get("solver_adj").map(local).unwrap_or_default(),
|
||||
"constraints": constraints,
|
||||
"nonNegative": get("solver_neg").as_deref().map(str::trim) != Some("2"),
|
||||
"method": method,
|
||||
});
|
||||
for (name, key) in [("solver_ssz", "populationSize"), ("solver_mrt", "mutationRate"), ("solver_rsd", "randomSeed")] {
|
||||
if let Some(v) = get(name).and_then(|v| num_text(&v)) {
|
||||
model[key] = json!(v);
|
||||
}
|
||||
}
|
||||
if let Some(v) = get("solver_rbv") {
|
||||
model["requireBounds"] = json!(v.trim() != "2");
|
||||
}
|
||||
Some(model)
|
||||
}
|
||||
|
||||
/// The model as `data.solver` params, in the form [`saved_model`] returns (the objective
|
||||
/// qualified by its sheet).
|
||||
pub(super) fn model_params(wb: &Workbook, m: &Model) -> Json {
|
||||
let sheet = wb.sheet(m.sheet).map(|s| crate::pivot::quote_sheet(&s.name)).unwrap_or_default();
|
||||
let (goal, value) = match m.goal {
|
||||
Goal::Max => ("max", 0.0),
|
||||
Goal::Min => ("min", 0.0),
|
||||
Goal::Value(v) => ("value", v),
|
||||
};
|
||||
let constraints: Vec<Json> = m
|
||||
.cons
|
||||
.iter()
|
||||
.map(|c| match c.rel {
|
||||
Rel::Int | Rel::Bin => json!(format!("{} {}", c.lhs_text, c.rel.text())),
|
||||
_ => json!(format!("{} {} {}", c.lhs_text, c.rel.text(), c.rhs_text)),
|
||||
})
|
||||
.collect();
|
||||
let method = match m.method {
|
||||
Method::Simplex => "simplex",
|
||||
Method::Grg => "grg",
|
||||
Method::Evolutionary => "evolutionary",
|
||||
};
|
||||
let mut out = json!({
|
||||
"objective": format!("{sheet}!{}", m.objective.a1()),
|
||||
"goal": goal,
|
||||
"value": value,
|
||||
"variables": m.vars_text,
|
||||
"constraints": constraints,
|
||||
"nonNegative": m.non_negative,
|
||||
"method": method,
|
||||
});
|
||||
if m.method == Method::Evolutionary {
|
||||
let e = &m.evo;
|
||||
out["populationSize"] = json!(e.population);
|
||||
out["mutationRate"] = json!(e.mutation);
|
||||
out["randomSeed"] = json!(e.seed);
|
||||
out["convergence"] = json!(e.convergence);
|
||||
out["generationsWithoutImprovement"] = json!(e.stall);
|
||||
out["requireBounds"] = json!(e.require_bounds);
|
||||
}
|
||||
out
|
||||
}
|
||||
@@ -0,0 +1,197 @@
|
||||
use super::lp::{Outcome, Row, Sense};
|
||||
use super::model::Goal;
|
||||
use super::problem::{Evaluator, Point};
|
||||
use super::*;
|
||||
|
||||
// ---------------------------------------------------------------- nonlinear
|
||||
|
||||
/// Augmented-Lagrangian value of a point: objective `phi` plus the constraint terms.
|
||||
fn aug(phi: f64, cons: &[f64], senses: &[Sense], lam: &[f64], mu: f64) -> f64 {
|
||||
let mut v = phi;
|
||||
for ((c, s), l) in cons.iter().zip(senses).zip(lam) {
|
||||
v += match s {
|
||||
Sense::Eq => l * c + 0.5 * mu * c * c,
|
||||
_ => {
|
||||
let t = (l + mu * c).max(0.0);
|
||||
(t * t - l * l) / (2.0 * mu)
|
||||
}
|
||||
};
|
||||
}
|
||||
v
|
||||
}
|
||||
|
||||
/// Constraint values in the `≤ 0` / `= 0` form, and the objective to minimise.
|
||||
pub(super) fn shape(p: &Point, rows: &[Row], goal: Goal) -> (f64, Vec<f64>) {
|
||||
let phi = match goal {
|
||||
Goal::Max => -p.f,
|
||||
Goal::Min => p.f,
|
||||
Goal::Value(v) => (p.f - v) * (p.f - v),
|
||||
};
|
||||
let cons = p.g.iter().zip(rows).map(|(g, r)| if r.sense == Sense::Ge { -g } else { *g }).collect();
|
||||
(phi, cons)
|
||||
}
|
||||
|
||||
pub(super) fn violation(cons: &[f64], senses: &[Sense]) -> f64 {
|
||||
cons.iter().zip(senses).map(|(c, s)| if *s == Sense::Eq { c.abs() } else { c.max(0.0) }).fold(0.0, f64::max)
|
||||
}
|
||||
|
||||
pub(super) fn project(x: &mut [f64], lo: &[f64], hi: &[f64]) {
|
||||
for ((v, l), h) in x.iter_mut().zip(lo).zip(hi) {
|
||||
if *v < *l {
|
||||
*v = *l;
|
||||
}
|
||||
if *v > *h {
|
||||
*v = *h;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Local minimisation of the model within bounds `lo..hi` from `x0`.
|
||||
pub(super) fn nlp(ev: &mut Evaluator, x0: &[f64], lo: &[f64], hi: &[f64], goal: Goal) -> Outcome {
|
||||
let n = x0.len();
|
||||
let rows = ev.rows.clone();
|
||||
let senses: Vec<Sense> = rows.iter().map(|r| r.sense).collect();
|
||||
let mut x = x0.to_vec();
|
||||
project(&mut x, lo, hi);
|
||||
let mut lam = vec![0.0; rows.len()];
|
||||
let mut mu = 10.0;
|
||||
let Some(p) = ev.eval(&x) else { return Outcome::Limit };
|
||||
let (mut phi, mut cons) = shape(&p, &rows, goal);
|
||||
if !phi.is_finite() {
|
||||
return Outcome::Infeasible;
|
||||
}
|
||||
let phi0 = phi;
|
||||
let mut prev_viol = violation(&cons, &senses);
|
||||
let mut converged = false;
|
||||
let mut out_of_budget = false;
|
||||
let mut carried: Option<f64> = None;
|
||||
for _outer in 0..40 {
|
||||
// Inner: projected gradient on the augmented Lagrangian.
|
||||
let lval = |phi: f64, cons: &[f64], lam: &[f64], mu: f64| aug(phi, cons, &senses, lam, mu);
|
||||
let grad = |ev: &mut Evaluator, x: &[f64], lam: &[f64], mu: f64| -> Option<Vec<f64>> {
|
||||
let mut gr = vec![0.0; n];
|
||||
for j in 0..n {
|
||||
let xj = x.get(j).copied().unwrap_or(0.0);
|
||||
let h = 1e-6 * xj.abs().max(1.0);
|
||||
let (l, u) = (lo.get(j).copied().unwrap_or(f64::NEG_INFINITY), hi.get(j).copied().unwrap_or(f64::INFINITY));
|
||||
let (a, b) = (if xj - h >= l { xj - h } else { xj }, if xj + h <= u { xj + h } else { xj });
|
||||
if b <= a {
|
||||
continue;
|
||||
}
|
||||
let mut xa = x.to_vec();
|
||||
let mut xb = x.to_vec();
|
||||
if let Some(v) = xa.get_mut(j) {
|
||||
*v = a;
|
||||
}
|
||||
if let Some(v) = xb.get_mut(j) {
|
||||
*v = b;
|
||||
}
|
||||
let (pa, ca) = shape(&ev.eval(&xa)?, &rows, goal);
|
||||
let (pb, cb) = shape(&ev.eval(&xb)?, &rows, goal);
|
||||
let d = (lval(pb, &cb, lam, mu) - lval(pa, &ca, lam, mu)) / (b - a);
|
||||
if let Some(v) = gr.get_mut(j) {
|
||||
*v = if d.is_finite() { d } else { 0.0 };
|
||||
}
|
||||
}
|
||||
Some(gr)
|
||||
};
|
||||
let mut l_cur = lval(phi, &cons, &lam, mu);
|
||||
let Some(mut g) = grad(ev, &x, &lam, mu) else {
|
||||
out_of_budget = true;
|
||||
break;
|
||||
};
|
||||
// The step length carries over between outer iterations once known.
|
||||
let gnorm = g.iter().map(|v| v.abs()).fold(0.0, f64::max);
|
||||
let mut alpha = carried.unwrap_or(if gnorm > 0.0 { 1.0 / gnorm } else { 1.0 });
|
||||
let mut inner_done = false;
|
||||
let mut quiet = 0;
|
||||
for _inner in 0..500 {
|
||||
let mut accepted = None;
|
||||
let mut a = alpha;
|
||||
for _ in 0..40 {
|
||||
let mut xn: Vec<f64> = x.iter().zip(&g).map(|(v, d)| v - a * d).collect();
|
||||
project(&mut xn, lo, hi);
|
||||
let Some(pn) = ev.eval(&xn) else { break };
|
||||
let (phn, cn) = shape(&pn, &rows, goal);
|
||||
let ln = lval(phn, &cn, &lam, mu);
|
||||
let decrease: f64 = g.iter().zip(x.iter().zip(&xn)).map(|(d, (u, v))| d * (u - v)).sum();
|
||||
if ln.is_finite() && ln <= l_cur - 1e-4 * decrease {
|
||||
accepted = Some((xn, phn, cn, ln));
|
||||
break;
|
||||
}
|
||||
a *= 0.25;
|
||||
}
|
||||
if ev.exhausted() {
|
||||
out_of_budget = true;
|
||||
}
|
||||
let Some((xn, phn, cn, ln)) = accepted else {
|
||||
inner_done = true;
|
||||
break;
|
||||
};
|
||||
// Variables or objective running off towards infinity: unbounded.
|
||||
if phn.abs() > 1e15 * (1.0 + phi0.abs()) || xn.iter().zip(x0).any(|(v, s)| v.abs() > 1e15 * (1.0 + s.abs())) {
|
||||
return Outcome::Unbounded;
|
||||
}
|
||||
let s: Vec<f64> = xn.iter().zip(&x).map(|(a, b)| a - b).collect();
|
||||
let step = s.iter().map(|v| v.abs()).fold(0.0, f64::max);
|
||||
let xmax = xn.iter().map(|v| v.abs()).fold(1.0, f64::max);
|
||||
let small = step <= 1e-10 * xmax || (l_cur - ln).abs() <= 1e-13 * (1.0 + l_cur.abs());
|
||||
quiet = if small { quiet + 1 } else { 0 };
|
||||
x = xn;
|
||||
phi = phn;
|
||||
cons = cn;
|
||||
l_cur = ln;
|
||||
if quiet >= 3 {
|
||||
inner_done = true;
|
||||
break;
|
||||
}
|
||||
let Some(gn) = grad(ev, &x, &lam, mu) else {
|
||||
out_of_budget = true;
|
||||
break;
|
||||
};
|
||||
// Barzilai–Borwein step.
|
||||
let yv: Vec<f64> = gn.iter().zip(&g).map(|(a, b)| a - b).collect();
|
||||
let sy: f64 = s.iter().zip(&yv).map(|(a, b)| a * b).sum();
|
||||
let ss: f64 = s.iter().map(|v| v * v).sum();
|
||||
alpha = if sy > 0.0 { (ss / sy).clamp(1e-12, 1e12) } else { (a * 4.0).min(1e12) };
|
||||
g = gn;
|
||||
if g.iter().all(|v| *v == 0.0) {
|
||||
inner_done = true;
|
||||
break;
|
||||
}
|
||||
}
|
||||
carried = Some(alpha);
|
||||
if out_of_budget {
|
||||
break;
|
||||
}
|
||||
let viol = violation(&cons, &senses);
|
||||
let scale = rows.iter().map(|r| r.rhs_const.abs()).fold(1.0, f64::max);
|
||||
if viol <= PRECISION * scale && inner_done {
|
||||
converged = true;
|
||||
if rows.is_empty() {
|
||||
break;
|
||||
}
|
||||
}
|
||||
// Multiplier and penalty updates.
|
||||
for ((l, c), s) in lam.iter_mut().zip(&cons).zip(&senses) {
|
||||
*l = match s {
|
||||
Sense::Eq => *l + mu * c,
|
||||
_ => (*l + mu * c).max(0.0),
|
||||
};
|
||||
}
|
||||
if viol > 0.25 * prev_viol {
|
||||
mu = (mu * 10.0).min(1e10);
|
||||
}
|
||||
if converged && viol <= PRECISION * scale * 1e-2 {
|
||||
break;
|
||||
}
|
||||
prev_viol = viol;
|
||||
converged = false;
|
||||
}
|
||||
let viol = violation(&cons, &senses);
|
||||
let scale = rows.iter().map(|r| r.rhs_const.abs()).fold(1.0, f64::max);
|
||||
if viol > PRECISION * scale * 10.0 {
|
||||
return if out_of_budget { Outcome::Limit } else { Outcome::Infeasible };
|
||||
}
|
||||
Outcome::Solved(x, phi, (converged || viol <= PRECISION * scale) && !out_of_budget)
|
||||
}
|
||||
@@ -0,0 +1,144 @@
|
||||
use super::lp::{Row, Sense};
|
||||
use super::model::{Model, Rel, Rhs};
|
||||
use super::*;
|
||||
|
||||
/// The values the solvers see at a point: the objective and `lhs − rhs` per row.
|
||||
#[derive(Clone, Debug)]
|
||||
pub(super) struct Point {
|
||||
pub(super) f: f64,
|
||||
pub(super) g: Vec<f64>,
|
||||
}
|
||||
|
||||
pub(super) struct Evaluator {
|
||||
pub(super) trial: Trial,
|
||||
pub(super) sheet: usize,
|
||||
pub(super) vars: Vec<CellRef>,
|
||||
pub(super) objective: CellRef,
|
||||
pub(super) rows: Vec<Row>,
|
||||
pub(super) evals: u64,
|
||||
pub(super) budget: u64,
|
||||
/// The Limits report: the objective is this variable itself.
|
||||
pub(super) target: Option<usize>,
|
||||
}
|
||||
|
||||
pub(super) fn as_num(v: &Value) -> f64 {
|
||||
match v {
|
||||
Value::Number(n) => *n,
|
||||
Value::Empty => 0.0,
|
||||
Value::Bool(b) => f64::from(u8::from(*b)),
|
||||
_ => f64::NAN,
|
||||
}
|
||||
}
|
||||
|
||||
impl Evaluator {
|
||||
pub(super) fn exhausted(&self) -> bool {
|
||||
self.evals >= self.budget
|
||||
}
|
||||
|
||||
/// Evaluates the model at `x`; `None` once the budget is spent.
|
||||
pub(super) fn eval(&mut self, x: &[f64]) -> Option<Point> {
|
||||
if self.exhausted() {
|
||||
return None;
|
||||
}
|
||||
self.evals += 1;
|
||||
let sets: Vec<(usize, CellRef, Value)> = self.vars.iter().zip(x).map(|(c, v)| (self.sheet, *c, Value::number(*v))).collect();
|
||||
let mut reads = Vec::with_capacity(1 + 2 * self.rows.len());
|
||||
reads.push((self.sheet, self.objective));
|
||||
for r in &self.rows {
|
||||
reads.push((self.sheet, r.lhs));
|
||||
if let Some(c) = r.rhs {
|
||||
reads.push((self.sheet, c));
|
||||
}
|
||||
}
|
||||
let vals = self.trial.run(&sets, &reads);
|
||||
let mut it = vals.iter().map(as_num);
|
||||
let f = it.next().unwrap_or(f64::NAN);
|
||||
let f = match self.target {
|
||||
Some(j) => x.get(j).copied().unwrap_or(f64::NAN),
|
||||
None => f,
|
||||
};
|
||||
let mut g = Vec::with_capacity(self.rows.len());
|
||||
for r in &self.rows {
|
||||
let l = it.next().unwrap_or(f64::NAN);
|
||||
let rv = if r.rhs.is_some() { it.next().unwrap_or(f64::NAN) } else { r.rhs_const };
|
||||
g.push(l - rv);
|
||||
}
|
||||
Some(Point { f, g })
|
||||
}
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------- problem
|
||||
|
||||
pub(super) struct Problem {
|
||||
pub(super) lo: Vec<f64>,
|
||||
pub(super) hi: Vec<f64>,
|
||||
pub(super) integer: Vec<bool>,
|
||||
pub(super) rows: Vec<Row>,
|
||||
}
|
||||
|
||||
/// Splits the constraints into variable bounds and general rows.
|
||||
pub(super) fn build_problem(m: &Model) -> Problem {
|
||||
let n = m.vars.len();
|
||||
let mut lo = vec![f64::NEG_INFINITY; n];
|
||||
let mut hi = vec![f64::INFINITY; n];
|
||||
let mut has_lower = vec![false; n];
|
||||
let mut integer = vec![false; n];
|
||||
let mut rows = Vec::new();
|
||||
let var_index = |c: &CellRef| m.vars.iter().position(|v| v == c);
|
||||
for c in &m.cons {
|
||||
for (k, cell) in c.lhs.iter().enumerate() {
|
||||
let j = var_index(cell);
|
||||
match c.rel {
|
||||
Rel::Int | Rel::Bin => {
|
||||
if let Some(j) = j {
|
||||
if let Some(b) = integer.get_mut(j) {
|
||||
*b = true;
|
||||
}
|
||||
if c.rel == Rel::Bin
|
||||
&& let (Some(l), Some(h), Some(hl)) = (lo.get_mut(j), hi.get_mut(j), has_lower.get_mut(j))
|
||||
{
|
||||
*l = l.max(0.0);
|
||||
*h = h.min(1.0);
|
||||
*hl = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
rel => {
|
||||
let sense = match rel {
|
||||
Rel::Le => Sense::Le,
|
||||
Rel::Ge => Sense::Ge,
|
||||
_ => Sense::Eq,
|
||||
};
|
||||
let (rhs, rhs_const) = match &c.rhs {
|
||||
Rhs::Const(v) => (None, *v),
|
||||
Rhs::Cells(cells) => (cells.get(k).or_else(|| cells.first()).copied(), 0.0),
|
||||
Rhs::None => (None, 0.0),
|
||||
};
|
||||
match (j, rhs) {
|
||||
// A constant bound on a variable cell is a box bound.
|
||||
(Some(j), None) => {
|
||||
if let (Some(l), Some(h), Some(hl)) = (lo.get_mut(j), hi.get_mut(j), has_lower.get_mut(j)) {
|
||||
if sense != Sense::Ge {
|
||||
*h = h.min(rhs_const);
|
||||
}
|
||||
if sense != Sense::Le {
|
||||
*l = l.max(rhs_const);
|
||||
*hl = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
_ => rows.push(Row { lhs: *cell, sense, rhs, rhs_const }),
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if m.non_negative {
|
||||
for (l, hl) in lo.iter_mut().zip(&has_lower) {
|
||||
if !hl {
|
||||
*l = l.max(0.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
Problem { lo, hi, integer, rows }
|
||||
}
|
||||
@@ -0,0 +1,373 @@
|
||||
use super::lp::Outcome;
|
||||
use super::model::{Constraint, Goal, Method, Model, Rel, Rhs, parse_model, saved_model};
|
||||
use super::problem::{Evaluator, Problem, as_num, build_problem};
|
||||
use super::*;
|
||||
|
||||
// ---------------------------------------------------------------- reports
|
||||
|
||||
#[derive(Clone, Copy, Debug, PartialEq)]
|
||||
pub(super) enum ReportKind {
|
||||
Answer,
|
||||
Limits,
|
||||
}
|
||||
|
||||
impl ReportKind {
|
||||
fn title(self) -> &'static str {
|
||||
match self {
|
||||
ReportKind::Answer => "Answer",
|
||||
ReportKind::Limits => "Limits",
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// The `reports` list (`None` when absent).
|
||||
pub(super) fn report_kinds(cmd: &'static str, p: &Json) -> Result<Option<Vec<ReportKind>>> {
|
||||
let Some(list) = p.get("reports") else { return Ok(None) };
|
||||
let list = list.as_array().ok_or_else(|| bad(cmd, "`reports` must be an array of answer and limits"))?;
|
||||
let mut kinds = Vec::new();
|
||||
for k in list {
|
||||
let kind = match k.as_str().map(|t| t.trim().to_ascii_lowercase()).as_deref() {
|
||||
Some("answer") => ReportKind::Answer,
|
||||
Some("limits") => ReportKind::Limits,
|
||||
_ => return Err(bad(cmd, format!("unknown report {k} (answer or limits)"))),
|
||||
};
|
||||
if !kinds.contains(&kind) {
|
||||
kinds.push(kind);
|
||||
}
|
||||
}
|
||||
Ok(Some(kinds))
|
||||
}
|
||||
|
||||
/// A report's cells: (row, column, value, bold).
|
||||
type Lines = Vec<(u32, u32, Value, bool)>;
|
||||
|
||||
fn abs_ref(c: CellRef) -> String {
|
||||
format!("${}${}", gridcraft_core::col_to_letters(c.col), u64::from(c.row) + 1)
|
||||
}
|
||||
|
||||
/// A cell's name in the reports: the nearest text to its left and the nearest text above it.
|
||||
fn cell_name(sh: &Sheet, c: CellRef) -> String {
|
||||
let text = |c: CellRef| match sh.value(c) {
|
||||
Value::Text(t) if !t.trim().is_empty() => Some(t.to_string()),
|
||||
_ => None,
|
||||
};
|
||||
let left = (0..c.col).rev().find_map(|col| text(CellRef::new(c.row, col)));
|
||||
let top = (0..c.row).rev().find_map(|row| text(CellRef::new(row, c.col)));
|
||||
[left, top].into_iter().flatten().collect::<Vec<_>>().join(" ")
|
||||
}
|
||||
|
||||
fn method_name(m: Method) -> &'static str {
|
||||
match m {
|
||||
Method::Simplex => "Simplex LP",
|
||||
Method::Grg => "GRG Nonlinear",
|
||||
Method::Evolutionary => "Evolutionary",
|
||||
}
|
||||
}
|
||||
|
||||
/// Header labels in bold from column 1, then the rows below them; returns the next free row.
|
||||
fn table(out: &mut Lines, row: u32, headers: &[&str], rows: Vec<Vec<Value>>) -> u32 {
|
||||
for (k, h) in headers.iter().enumerate() {
|
||||
out.push((row, 1 + k as u32, Value::text(*h), true));
|
||||
}
|
||||
let mut r = row;
|
||||
for cells in rows {
|
||||
r = r.saturating_add(1);
|
||||
for (k, v) in cells.into_iter().enumerate() {
|
||||
out.push((r, 1 + k as u32, v, false));
|
||||
}
|
||||
}
|
||||
r.saturating_add(2)
|
||||
}
|
||||
|
||||
/// Everything a report describes: the model, its starting and final values, and the result.
|
||||
pub(super) struct ReportInput<'a> {
|
||||
pub(super) model: &'a Model,
|
||||
pub(super) prob: &'a Problem,
|
||||
pub(super) original: Vec<f64>,
|
||||
pub(super) solution: Vec<f64>,
|
||||
pub(super) message: String,
|
||||
pub(super) budget: u64,
|
||||
}
|
||||
|
||||
fn title(out: &mut Lines, kind: ReportKind, sheet: &str) -> u32 {
|
||||
out.push((0, 0, Value::text(format!("GridCraft {} Report", kind.title())), true));
|
||||
out.push((1, 0, Value::text(format!("Worksheet: {sheet}")), false));
|
||||
3
|
||||
}
|
||||
|
||||
fn format_constraint_row(sh: &Sheet, c: &Constraint, cell: CellRef, v: f64, rhs: &[f64], k: usize) -> Vec<Value> {
|
||||
let (formula, status, slack) = match (c.rel, &c.rhs) {
|
||||
(Rel::Int, _) => (format!("{}=Integer", abs_ref(cell)), "Binding", 0.0),
|
||||
(Rel::Bin, _) => (format!("{}=Binary", abs_ref(cell)), "Binding", 0.0),
|
||||
(rel, rhs_cells) => {
|
||||
let b = rhs.get(k).or_else(|| rhs.first()).copied().unwrap_or(f64::NAN);
|
||||
let rhs_text = match rhs_cells {
|
||||
Rhs::Cells(cells) => cells.get(k).or_else(|| cells.first()).map(|c| abs_ref(*c)).unwrap_or_default(),
|
||||
_ => gridcraft_core::number_to_text(b),
|
||||
};
|
||||
let slack = (v - b).abs();
|
||||
let binding = rel == Rel::Eq || slack <= PRECISION * (1.0 + b.abs());
|
||||
let slack = if binding { 0.0 } else { slack };
|
||||
(format!("{}{}{rhs_text}", abs_ref(cell), rel.text()), if binding { "Binding" } else { "Not Binding" }, slack)
|
||||
}
|
||||
};
|
||||
let num = |v: f64| if v.is_finite() { Value::number(v) } else { Value::Error(gridcraft_core::CellError::Value) };
|
||||
vec![Value::text(abs_ref(cell)), Value::text(cell_name(sh, cell)), num(v), Value::text(formula), Value::text(status), num(slack)]
|
||||
}
|
||||
|
||||
fn answer_report(d: &DocState, sh: &Sheet, inp: &ReportInput) -> Lines {
|
||||
let m = inp.model;
|
||||
let mut trial = Trial::new(d);
|
||||
let read = |trial: &mut Trial, x: &[f64], cells: &[CellRef]| -> Vec<f64> {
|
||||
let sets: Vec<(usize, CellRef, Value)> = m.vars.iter().zip(x).map(|(c, v)| (m.sheet, *c, Value::number(*v))).collect();
|
||||
let reads: Vec<(usize, CellRef)> = cells.iter().map(|c| (m.sheet, *c)).collect();
|
||||
trial.run(&sets, &reads).iter().map(as_num).collect()
|
||||
};
|
||||
let f0 = read(&mut trial, &inp.original, &[m.objective]).first().copied().unwrap_or(f64::NAN);
|
||||
// Every cell the constraints read, at the solution.
|
||||
let mut cells = vec![m.objective];
|
||||
for c in &m.cons {
|
||||
cells.extend(&c.lhs);
|
||||
if let Rhs::Cells(r) = &c.rhs {
|
||||
cells.extend(r);
|
||||
}
|
||||
}
|
||||
let vals = read(&mut trial, &inp.solution, &cells);
|
||||
let mut at = vals.iter().copied();
|
||||
let f1 = at.next().unwrap_or(f64::NAN);
|
||||
let num = |v: f64| if v.is_finite() { Value::number(v) } else { Value::Error(gridcraft_core::CellError::Value) };
|
||||
let mut out = Lines::new();
|
||||
let mut r = title(&mut out, ReportKind::Answer, &sh.name);
|
||||
out.push((r, 0, Value::text(format!("Result: {}", inp.message)), false));
|
||||
out.push((r + 1, 0, Value::text(format!("Engine: {}", method_name(m.method))), false));
|
||||
r += 3;
|
||||
let goal = match m.goal {
|
||||
Goal::Max => "Max".to_string(),
|
||||
Goal::Min => "Min".to_string(),
|
||||
Goal::Value(v) => format!("Value Of {v}"),
|
||||
};
|
||||
out.push((r, 0, Value::text(format!("Objective Cell ({goal})")), true));
|
||||
let obj = vec![vec![Value::text(abs_ref(m.objective)), Value::text(cell_name(sh, m.objective)), num(f0), num(f1)]];
|
||||
r = table(&mut out, r + 1, &["Cell", "Name", "Original Value", "Final Value"], obj);
|
||||
out.push((r, 0, Value::text("Variable Cells"), true));
|
||||
let kind = |j: usize| {
|
||||
let bin = m.cons.iter().any(|c| c.rel == Rel::Bin && m.vars.get(j).is_some_and(|v| c.lhs.contains(v)));
|
||||
match (bin, inp.prob.integer.get(j).copied().unwrap_or(false)) {
|
||||
(true, _) => "Binary",
|
||||
(false, true) => "Integer",
|
||||
_ => "Contin",
|
||||
}
|
||||
};
|
||||
let vars = m
|
||||
.vars
|
||||
.iter()
|
||||
.enumerate()
|
||||
.map(|(j, c)| {
|
||||
let (a, b) = (inp.original.get(j).copied().unwrap_or(f64::NAN), inp.solution.get(j).copied().unwrap_or(f64::NAN));
|
||||
vec![Value::text(abs_ref(*c)), Value::text(cell_name(sh, *c)), num(a), num(b), Value::text(kind(j))]
|
||||
})
|
||||
.collect();
|
||||
r = table(&mut out, r + 1, &["Cell", "Name", "Original Value", "Final Value", "Integer"], vars);
|
||||
out.push((r, 0, Value::text("Constraints"), true));
|
||||
let mut rows = Vec::new();
|
||||
for c in &m.cons {
|
||||
let lhs: Vec<f64> = c.lhs.iter().map(|_| at.next().unwrap_or(f64::NAN)).collect();
|
||||
let rhs: Vec<f64> = match &c.rhs {
|
||||
Rhs::Cells(cells) => cells.iter().map(|_| at.next().unwrap_or(f64::NAN)).collect(),
|
||||
Rhs::Const(v) => vec![*v],
|
||||
Rhs::None => vec![],
|
||||
};
|
||||
for (k, (cell, v)) in c.lhs.iter().zip(&lhs).enumerate() {
|
||||
rows.push(format_constraint_row(sh, c, *cell, *v, &rhs, k));
|
||||
}
|
||||
}
|
||||
table(&mut out, r + 1, &["Cell", "Name", "Cell Value", "Formula", "Status", "Slack"], rows);
|
||||
out
|
||||
}
|
||||
|
||||
fn limits_report(d: &DocState, sh: &Sheet, inp: &ReportInput) -> Result<Lines> {
|
||||
let m = inp.model;
|
||||
if inp.prob.integer.iter().any(|b| *b) {
|
||||
return Err(EngineError::Other("The Limits report doesn't apply to models with integer or binary constraints.".into()));
|
||||
}
|
||||
if m.method == Method::Evolutionary {
|
||||
return Err(EngineError::Other("The Limits report is available for Simplex LP and GRG Nonlinear.".into()));
|
||||
}
|
||||
let n = m.vars.len();
|
||||
let x = &inp.solution;
|
||||
// Each variable is minimised and maximised with the others held at the solution.
|
||||
let per = (inp.budget / (2 * n as u64).max(1)).max(n as u64 + 3);
|
||||
let mut ev = Evaluator {
|
||||
trial: Trial::new(d),
|
||||
sheet: m.sheet,
|
||||
vars: m.vars.clone(),
|
||||
objective: m.objective,
|
||||
rows: inp.prob.rows.clone(),
|
||||
evals: 0,
|
||||
budget: 0,
|
||||
target: None,
|
||||
};
|
||||
let objective_at = |ev: &mut Evaluator, x: &[f64]| {
|
||||
ev.target = None;
|
||||
ev.budget = ev.evals.saturating_add(1);
|
||||
ev.eval(x).map(|p| p.f).unwrap_or(f64::NAN)
|
||||
};
|
||||
let num = |v: f64| if v.is_finite() { Value::number(v) } else { Value::Error(gridcraft_core::CellError::Value) };
|
||||
let f = objective_at(&mut ev, x);
|
||||
let mut rows = Vec::new();
|
||||
for (j, c) in m.vars.iter().enumerate() {
|
||||
let (mut lo, mut hi) = (x.clone(), x.clone());
|
||||
if let (Some(l), Some(h), Some(pl), Some(ph)) = (lo.get_mut(j), hi.get_mut(j), inp.prob.lo.get(j), inp.prob.hi.get(j)) {
|
||||
*l = *pl;
|
||||
*h = *ph;
|
||||
}
|
||||
let mut ends = Vec::new();
|
||||
for goal in [Goal::Min, Goal::Max] {
|
||||
ev.target = Some(j);
|
||||
ev.evals = 0;
|
||||
ev.budget = per;
|
||||
let solve = optimise(&mut ev, m, goal, x, &lo, &hi, &inp.prob.integer, 0.0);
|
||||
ends.push(match solve.outcome {
|
||||
Outcome::Solved(p, _, _) => {
|
||||
let v = p.get(j).copied().unwrap_or(f64::NAN);
|
||||
(num(v), num(objective_at(&mut ev, &p)))
|
||||
}
|
||||
Outcome::Unbounded => (Value::text("Unbounded"), Value::Empty),
|
||||
_ => (Value::Error(gridcraft_core::CellError::NA), Value::Empty),
|
||||
});
|
||||
}
|
||||
let mut row = vec![Value::text(abs_ref(*c)), Value::text(cell_name(sh, *c)), num(x.get(j).copied().unwrap_or(f64::NAN))];
|
||||
for (v, o) in ends {
|
||||
row.push(v);
|
||||
row.push(o);
|
||||
}
|
||||
rows.push(row);
|
||||
}
|
||||
let mut out = Lines::new();
|
||||
let r = title(&mut out, ReportKind::Limits, &sh.name);
|
||||
let obj = vec![vec![Value::text(abs_ref(m.objective)), Value::text(cell_name(sh, m.objective)), num(f)]];
|
||||
let r = table(&mut out, r, &["Objective Cell", "Name", "Value"], obj);
|
||||
table(&mut out, r, &["Variable Cell", "Name", "Value", "Lower Limit", "Objective Result", "Upper Limit", "Objective Result"], rows);
|
||||
Ok(out)
|
||||
}
|
||||
|
||||
/// Builds every report before the solution or saved model can change.
|
||||
pub(super) fn prepare_reports(d: &DocState, kinds: &[ReportKind], inp: &ReportInput) -> Result<Vec<(String, Lines)>> {
|
||||
if kinds.is_empty() {
|
||||
return Ok(Vec::new());
|
||||
}
|
||||
if d.wb.protected_structure {
|
||||
return Err(EngineError::Other("The workbook is protected and cannot be changed.".into()));
|
||||
}
|
||||
let sh = d.wb.sheet(inp.model.sheet).ok_or(EngineError::NoDocument)?;
|
||||
let mut built = Vec::new();
|
||||
for k in kinds {
|
||||
let lines = match k {
|
||||
ReportKind::Answer => answer_report(d, sh, inp),
|
||||
ReportKind::Limits => limits_report(d, sh, inp)?,
|
||||
};
|
||||
let name = (1..=10_000u32)
|
||||
.map(|i| format!("{} Report {i}", k.title()))
|
||||
.find(|n| d.wb.sheet_index(n).is_none() && !built.iter().any(|(b, _): &(String, Lines)| b == n))
|
||||
.ok_or_else(|| EngineError::Other("Too many report sheets.".into()))?;
|
||||
built.push((name, lines));
|
||||
}
|
||||
Ok(built)
|
||||
}
|
||||
|
||||
/// Inserts prepared reports before the model's sheet, which stays active.
|
||||
pub(super) fn write_reports(cx: &mut Ctx, at: usize, built: Vec<(String, Lines)>) -> Vec<String> {
|
||||
let names: Vec<String> = built.iter().map(|(n, _)| n.clone()).collect();
|
||||
let count = built.len();
|
||||
if count == 0 {
|
||||
return names;
|
||||
}
|
||||
let bold = cx.wb.styles.derive(gridcraft_model::StyleId::DEFAULT, |st| st.font.bold = true);
|
||||
for (i, (name, lines)) in built.into_iter().enumerate() {
|
||||
let mut sheet = Sheet::new(name);
|
||||
for (r, c, v, b) in lines {
|
||||
let style = if b { bold } else { gridcraft_model::StyleId::DEFAULT };
|
||||
sheet.set_cell(CellRef::new(r, c), Cell { value: v, formula: None, style });
|
||||
}
|
||||
for (c, w) in [(0u32, 14.0f32), (1, 110.0), (2, 140.0), (3, 120.0), (4, 120.0), (5, 120.0), (6, 120.0), (7, 120.0)] {
|
||||
sheet.cols.insert(c, gridcraft_model::LineInfo { size: Some(w), ..Default::default() });
|
||||
}
|
||||
cx.wb.sheets.insert(at + i, std::sync::Arc::new(sheet));
|
||||
}
|
||||
for nm in cx.wb.names.iter_mut() {
|
||||
if let Some(sc) = nm.scope.as_mut()
|
||||
&& *sc >= at
|
||||
{
|
||||
*sc += count;
|
||||
}
|
||||
}
|
||||
if cx.wb.active_sheet() >= at {
|
||||
cx.wb.set_active_sheet(cx.wb.active_sheet().saturating_add(count));
|
||||
}
|
||||
cx.structural = true;
|
||||
names
|
||||
}
|
||||
|
||||
pub(super) fn insert_report_selections(s: &mut Session, at: usize, count: usize) {
|
||||
if let Ok(d) = s.doc_mut() {
|
||||
for i in 0..count {
|
||||
d.sheet_selections.insert((at + i).min(d.sheet_selections.len()), crate::Selection::default());
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// `data.solverReport`: reports for a solution of the sheet's saved model.
|
||||
pub(super) fn report(s: &mut Session, p: &Json) -> Result<Json> {
|
||||
let kinds = report_kinds(REPORT_CMD, p)?.ok_or_else(|| bad(REPORT_CMD, "missing `reports` (answer, limits)"))?;
|
||||
if kinds.is_empty() {
|
||||
return Err(bad(REPORT_CMD, "`reports` is empty (answer, limits)"));
|
||||
}
|
||||
let d = s.doc()?;
|
||||
// The given model as is, else the saved model with any params given here over it.
|
||||
let params = match p.get("model") {
|
||||
Some(m @ Json::Object(_)) => m.clone(),
|
||||
Some(Json::Null) | None => {
|
||||
let mut params = saved_model(&d.wb, d.wb.active_sheet()).unwrap_or_else(|| json!({}));
|
||||
if let (Some(base), Some(over)) = (params.as_object_mut(), p.as_object()) {
|
||||
for (k, v) in over {
|
||||
base.insert(k.clone(), v.clone());
|
||||
}
|
||||
}
|
||||
params
|
||||
}
|
||||
Some(_) => return Err(bad(REPORT_CMD, "`model` must be data.solver params")),
|
||||
};
|
||||
let model = parse_model(s, ¶ms)?;
|
||||
let prob = build_problem(&model);
|
||||
let d = s.doc()?;
|
||||
let sh = d.wb.sheet(model.sheet).ok_or(EngineError::NoDocument)?;
|
||||
let current: Vec<f64> = model.vars.iter().map(|c| as_num(&sh.value(*c))).collect();
|
||||
// `{cell: number}` maps; cells they leave out keep their current values.
|
||||
let point = |key: &str| -> Result<Vec<f64>> {
|
||||
let mut x = current.clone();
|
||||
let Some(map) = p.get(key) else { return Ok(x) };
|
||||
let map = map.as_object().ok_or_else(|| bad(REPORT_CMD, format!("`{key}` must map cells to numbers")))?;
|
||||
for (cell, v) in map {
|
||||
let c = CellRef::parse(&cell.replace('$', "")).ok_or_else(|| bad(REPORT_CMD, format!("not a cell: `{cell}`")))?;
|
||||
let j = model.vars.iter().position(|v| *v == c).ok_or_else(|| bad(REPORT_CMD, format!("`{cell}` is not a variable cell")))?;
|
||||
let v = v.as_f64().filter(|v| v.is_finite()).ok_or_else(|| bad(REPORT_CMD, format!("`{key}.{cell}` must be a number")))?;
|
||||
if let Some(e) = x.get_mut(j) {
|
||||
*e = v;
|
||||
}
|
||||
}
|
||||
Ok(x)
|
||||
};
|
||||
let solution = point("values")?;
|
||||
let original = if p.get("original").is_some() { point("original")? } else { solution.clone() };
|
||||
let budget = params
|
||||
.get("maxEvaluations")
|
||||
.and_then(Json::as_f64)
|
||||
.filter(|v| v.is_finite() && *v >= 1.0)
|
||||
.map_or(DEFAULT_EVALS, |v| (v as u64).min(MAX_EVALS));
|
||||
let message = str_param(p, "message").unwrap_or("").chars().take(500).collect();
|
||||
let inp = ReportInput { model: &model, prob: &prob, original, solution, message, budget };
|
||||
let built = prepare_reports(s.doc()?, &kinds, &inp)?;
|
||||
let sheets = edit(s, |cx| Ok(write_reports(cx, model.sheet, built)))?;
|
||||
insert_report_selections(s, model.sheet, sheets.len());
|
||||
Ok(json!({"sheets": sheets}))
|
||||
}
|
||||
Reference in New Issue
Block a user