diff --git a/README.md b/README.md index 22a16cd..fc59fce 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,7 @@ Currently, Argon supports the following features: - Live reload of GUI upon changes in code editor - Parametric cells - Hierarchy -- General linear constraint solving (slow) +- Linear constraint solving: fast sparse elimination, with a general (dense) solver as fallback - Basic diagnostic reporting in the code editor - Basic detection of under/overconstrained systems diff --git a/bench/README.md b/bench/README.md index 327da0c..f4b449c 100644 --- a/bench/README.md +++ b/bench/README.md @@ -95,7 +95,7 @@ how an axis scales — without editing any source. Pass a comma-separated list: | `ARGON_BENCH_SHAPES` | shapes (recursion) | `500,1000,2000,4000,8000,16000,32000` | | `ARGON_BENCH_SHAPES_LOOP` | shapes (`for` loop) | `500,1000,2000,4000,8000,16000,32000` | | `ARGON_BENCH_INSTANCES` | instances | `500,…,64000` | -| `ARGON_BENCH_CONSTRAINTS` | coupled constraints | `32,64,128,256,512,1024` | +| `ARGON_BENCH_CONSTRAINTS` | coupled constraints | `32,64,128,256,512,1024,2048,4096,8192,16384` | | `ARGON_BENCH_HIER_SINGLE` | hierarchy (1 ref) | `4,8,16,32,48,64,96,128` | | `ARGON_BENCH_HIER_DOUBLE` | hierarchy (2 refs) | `4,8,16,32,48,64,96,128` | @@ -131,7 +131,7 @@ parameter; "peak" is peak heap allocated during compilation. | Shapes (recursion) | 32 000 rects | 1.52 s | 0.89 GiB | **~linear** (time `∝ n^1.2`, mem `∝ n^1.0`) | | Instances | 64 000 insts | 3.08 s | 1.26 GiB | **~linear** (time `∝ n^1.2`, mem `∝ n^1.0`) | | Hierarchy, 1 child ref | depth 128 | 0.005 s | 11 MiB | **linear** in depth | -| Coupled constraints | 1 024 rects | 22.0 s | 0.12 GiB | **super-cubic in time** (see below) | +| Coupled constraints | 16 384 rects | 1.76 s | 0.59 GiB | **~linear** (time `∝ n^1.04`, mem `∝ n^0.90`) | | Shapes (`for`-loop) | 32 000 rects | 1.06 s | 0.85 GiB | **~linear** (time `∝ n^1.2`, mem `∝ n^1.0`) | | Hierarchy, 2 child refs | depth 128 | 0.006 s | 11 MiB | **linear** in depth (was exponential before the shared-type fix) | @@ -145,16 +145,23 @@ parameter; "peak" is peak heap allocated during compilation. back-substitution without ever forming a matrix. This is the common case for real parametric cells and it scales comfortably to "thousands of rectangles". -- **Coupled constraints are the expensive axis.** When constraints form one - large connected component that *cannot* be back-substituted (here, a ring of - mutually-coupled edges), Argon falls back to its general linear solver, which - builds a dense matrix and takes an SVD. The per-doubling cost climbs from ~4× - at `n=64→128` to ~15× at `n=512→1024`, i.e. it steepens toward the `O(n^3)` - of dense factorization (and worse, because `solve()` is re-run as the system - is assembled). This is the "general linear constraint solving (slow)" caveat - in the top-level README, quantified: ~1 000 coupled editable variables take - ~22 s. Layouts whose constraints decompose into many small independent groups - (the typical case) avoid this entirely. +- **Coupled constraints now scale linearly too.** When constraints form one + large connected component that 1-variable back-substitution cannot crack + (here, a ring of mutually-coupled edges), the solver first runs a sparse + *elimination pre-pass*: it generalizes back-substitution to 2-variable + "definitional" constraints, expressing one variable in terms of another and + substituting it out of the few constraints that mention it. Because a + 2-variable constraint replaces one variable with one variable, this never + increases any constraint's size — the ring's leaf edges drop out and its + chain telescopes to a single closure, so the system is solved in `O(n)` and + the dense SVD never runs. The general dense solver is retained only as a + fallback for an *irreducible* coupled core (a block with no ≤2-variable + pivot). This axis used to be **super-cubic** — `~22 s` at `n=1024`, steepening + toward the `O(n^3)` of dense factorization — and is now `~n^1.04` in time + (`1.76 s` at `n=16384`, 16× larger, in less memory than the old `n=1024` + dense matrix used). The "general linear constraint solving (slow)" caveat in + the top-level README now bites only for genuinely dense coupled blocks, not + for the common sparse-but-coupled case. - **Hierarchy depth scales linearly.** A cell's static type (`CellTy`) records the structural type of every field, including instantiated sub-cells. That @@ -185,12 +192,14 @@ parameter; "peak" is peak heap allocated during compilation. constant — `shapes_loop` is even marginally faster, as the native `range` avoids the per-element recursion overhead of `emit_shapes`. -The takeaways for the paper: editable-object count, instance count, and -hierarchy depth all scale linearly; the one practically-relevant limit is the -dense general constraint solver on large *coupled* systems, which lines up with -the future-work item already listed in the project README (faster linear -constraint solving). The bullets above describe the build at the time of -measurement; because every axis is re-runnable (and size-configurable), the same -harness can be used to confirm improvements from compiler optimizations. +The takeaways for the paper: editable-object count, instance count, hierarchy +depth, and now coupled-constraint count all scale linearly. The dense general +solver is no longer a practical limit for sparse-but-coupled systems — the +elimination pre-pass reduces them first — and remains the fallback only for +irreducible dense constraint blocks with no ≤2-variable pivot, which lines up +with the "faster linear constraint solving" future-work item in the project +README. The bullets above describe the build at the time of measurement; because +every axis is re-runnable (and size-configurable), the same harness can be used +to confirm improvements from compiler optimizations. ![Argon scaling](argon_scaling.png) diff --git a/bench/argon_scaling.pdf b/bench/argon_scaling.pdf index 7012fd3..6b49e33 100644 Binary files a/bench/argon_scaling.pdf and b/bench/argon_scaling.pdf differ diff --git a/bench/argon_scaling.png b/bench/argon_scaling.png index a3d2761..7501131 100644 Binary files a/bench/argon_scaling.png and b/bench/argon_scaling.png differ diff --git a/bench/results/constraints.csv b/bench/results/constraints.csv index 4415917..d44fc70 100644 --- a/bench/results/constraints.csv +++ b/bench/results/constraints.csv @@ -1,7 +1,11 @@ size,time_s,peak_bytes,n_objects -32,0.003860995,2558554,33 -64,0.00719482,3802170,65 -128,0.029517606,6722670,129 -256,0.224260882,15369454,257 -512,1.435889679,42100206,513 -1024,21.963011181,133310430,1025 +32,0.003404379,2558933,33 +64,0.004213333,3802613,65 +128,0.007537572,6289973,129 +256,0.014580103,11264693,257 +512,0.02941705,21214133,513 +1024,0.063108165,41113013,1025 +2048,0.137085263,80910773,2049 +4096,0.306775942,160506293,4097 +8192,0.819698386,319697317,8193 +16384,1.762328533,638079381,16385 diff --git a/core/compiler/src/lib.rs b/core/compiler/src/lib.rs index d48da97..de0bf1e 100644 --- a/core/compiler/src/lib.rs +++ b/core/compiler/src/lib.rs @@ -371,8 +371,9 @@ mod tests { write_bench_csv("shapes_loop", &rows); } - /// Axis 2: number of mutually-coupled constraints solved by the general - /// (dense) linear-constraint solver. + /// Axis 2: number of mutually-coupled constraints. The coupled ring is reduced by + /// the solver's sparse elimination pre-pass (2-variable constraints telescope away); + /// the general dense SVD remains as a fallback for any irreducible coupled core. #[test] #[ignore = "scaling benchmark; run in release, serially: cargo test -p compiler --release -- --ignored --test-threads=1 bench_"] fn bench_constraints() { @@ -381,7 +382,10 @@ mod tests { assert!(o.static_errors().is_empty(), "{:?}", o.static_errors()); let ast = o.ast(); let mut rows = Vec::new(); - for &n in &bench_sizes("ARGON_BENCH_CONSTRAINTS", &[32, 64, 128, 256, 512, 1024]) { + for &n in &bench_sizes( + "ARGON_BENCH_CONSTRAINTS", + &[32, 64, 128, 256, 512, 1024, 2048, 4096, 8192, 16384], + ) { let (dt, mem, out) = measure(1, || { compile( &ast, diff --git a/core/compiler/src/solver.rs b/core/compiler/src/solver.rs index bab5cf8..6ad4b82 100644 --- a/core/compiler/src/solver.rs +++ b/core/compiler/src/solver.rs @@ -25,6 +25,13 @@ pub struct Solver { back_substitute_stack: Vec, inconsistent_constraints: IndexSet, invalid_rounding: IndexSet, + // Per-`solve()` scratch for the sparse elimination pre-pass (`eliminate_definitional`). + // `elim_worklist` holds constraints to (re)examine for a small pivot; `substitutions` + // records `var = expr` definitions for variables eliminated via a 2-variable + // constraint, resolved into numbers afterwards by `resolve_substitutions`. Both are + // cleared at the start of each elimination pass, so they hold no state between solves. + elim_worklist: VecDeque, + substitutions: Vec<(Var, LinearExpr)>, } fn round(x: f64) -> f64 { @@ -160,6 +167,17 @@ impl Solver { if self.unsolved_vars.is_empty() || self.constraints.is_empty() { return; } + // Sparsity-exploiting pre-pass: peel off variables that are uniquely defined + // by a constraint of size <= 2 (generalizing 1-variable back-substitution), + // shrinking the system before the dense SVD below. Variables eliminated via a + // 2-variable constraint are expressed in terms of another variable and recorded + // in `self.substitutions`; their numeric values are recovered by + // `resolve_substitutions` once the remaining (irreducible) core has been solved. + // For systems whose constraints are all <= 2 variables (e.g. the coupled ring in + // `bench_constraints`) this resolves everything in O(n) and the SVD never runs; + // for a genuinely dense block it is a no-op and behaviour is identical to before. + self.eliminate_definitional(); + for component in self.constraint_components() { self.solve_component(&component.vars, &component.constraints); } @@ -173,6 +191,164 @@ impl Solver { } self.constraints .retain(|_, constraint| !constraint.coeffs.is_empty()); + + self.resolve_substitutions(); + } + + /// Sparse elimination pre-pass. Repeatedly examines constraints with at most two + /// variables and uses each to eliminate one of its variables, substituting it out + /// of the (few) other constraints that mention it. Because a size-2 constraint + /// expresses a variable as (one variable + constant), substitution replaces one + /// term with one term and so never increases any constraint's variable count: the + /// system only shrinks, and the pass runs in O(nnz). Constraints with > 2 variables + /// are left untouched for the dense path in `solve_component`. + fn eliminate_definitional(&mut self) { + self.substitutions.clear(); + self.elim_worklist.clear(); + self.elim_worklist.extend(self.constraints.keys().copied()); + while let Some(id) = self.elim_worklist.pop_front() { + let Some(constraint) = self.constraints.get_mut(&id) else { + continue; + }; + constraint.simplify(&self.solved_vars); + coalesce_terms(constraint); + let len = constraint.coeffs.len(); + let constant = constraint.constant; + match len { + 0 => { + if relative_ne!(constant, 0., epsilon = EPSILON) { + self.inconsistent_constraints.insert(id); + } + self.remove_constraint(id); + } + 1 => self.eliminate_unary(id), + 2 => self.eliminate_binary(id), + _ => {} + } + } + } + + /// Removes a constraint from the system and unlinks it from `var_to_constraints`. + fn remove_constraint(&mut self, id: ConstraintId) { + if let Some(constraint) = self.constraints.swap_remove(&id) { + for (_, var) in &constraint.coeffs { + if let Some(set) = self.var_to_constraints.get_mut(var) { + set.swap_remove(&id); + } + } + } + } + + /// Re-queues every still-live constraint mentioning `var`; they may have just + /// shrunk to a size the pre-pass can act on. + fn requeue_neighbors(&mut self, var: Var) { + let neighbors: Vec = self + .var_to_constraints + .get(&var) + .into_iter() + .flatten() + .copied() + .collect(); + self.elim_worklist.extend(neighbors); + } + + /// Grounds the single variable of a 1-variable constraint (post-simplify). The + /// 1-variable analogue of `eliminate_binary`; mirrors `try_back_substitute`. + fn eliminate_unary(&mut self, id: ConstraintId) { + let (coeff, var) = self.constraints[&id].coeffs[0]; + let val = -self.constraints[&id].constant / coeff; + self.remove_constraint(id); + self.requeue_neighbors(var); + self.assign_var(var, val); + } + + /// Eliminates one variable of a 2-variable constraint (post-simplify) by expressing + /// it in terms of the other and substituting it out of every other constraint that + /// mentions it. Records the definition in `substitutions` for later resolution. + fn eliminate_binary(&mut self, id: ConstraintId) { + let (c0, v0) = self.constraints[&id].coeffs[0]; + let (c1, v1) = self.constraints[&id].coeffs[1]; + let constant = self.constraints[&id].constant; + // Pivot on the larger-magnitude coefficient (partial-pivoting analogue). + let ((a, v), (cw, w)) = if c0.abs() >= c1.abs() { + ((c0, v0), (c1, v1)) + } else { + ((c1, v1), (c0, v0)) + }; + if a.abs() <= EPSILON || v == w { + return; + } + // From `a*v + cw*w + constant = 0`: v = (-cw/a) * w + (-constant/a). + let v_expr = LinearExpr { + coeffs: vec![(-cw / a, w)], + constant: -constant / a, + }; + self.remove_constraint(id); + let neighbors: Vec = self + .var_to_constraints + .get(&v) + .into_iter() + .flatten() + .copied() + .collect(); + for nid in neighbors { + let Some(constraint) = self.constraints.get_mut(&nid) else { + continue; + }; + substitute_var(constraint, v, &v_expr); + self.var_to_constraints.entry(w).or_default().insert(nid); + self.elim_worklist.push_back(nid); + } + self.var_to_constraints.swap_remove(&v); + self.unsolved_vars.swap_remove(&v); + self.substitutions.push((v, v_expr)); + } + + /// Recovers numeric values for variables eliminated by `eliminate_binary`. Walks + /// `substitutions` in reverse (reverse-topological) order: by the time each entry is + /// reached, every variable its expression depends on is either solved by the core + /// SVD / back-substitution or resolved earlier in this walk. A variable whose + /// expression does not become ground belongs to an under-determined component; its + /// defining constraint is restored (a row-equivalent of the pivot row that was + /// removed) so the under-constrained diagnostics in `rowspace_vecs` are unchanged. + fn resolve_substitutions(&mut self) { + while let Some((var, mut expr)) = self.substitutions.pop() { + expr.simplify(&self.solved_vars); + if expr.coeffs.is_empty() { + self.assign_var(var, expr.constant); + } else { + self.unsolved_vars.insert(var); + let id = self.next_constraint; + self.next_constraint += 1; + let mut coeffs = Vec::with_capacity(expr.coeffs.len() + 1); + coeffs.push((1., var)); + for (c, v) in expr.coeffs { + coeffs.push((-c, v)); + } + let constraint = LinearExpr { + coeffs, + constant: -expr.constant, + }; + for (_, v) in &constraint.coeffs { + self.var_to_constraints.entry(*v).or_default().insert(id); + } + self.constraints.insert(id, constraint); + } + } + } + + /// Rounds `val` to the solver grid, flags off-grid values, and records the solution. + /// Shared by the elimination pre-pass; matches the rounding contract used by + /// `try_back_substitute` and `solve_component`. + fn assign_var(&mut self, var: Var, val: f64) { + if self.solved_vars.contains_key(&var) { + return; + } + let rounded = round(val); + if relative_ne!(val, rounded, epsilon = EPSILON) { + self.invalid_rounding.insert(var); + } + self.solve_var(var, rounded); } pub fn rowspace_vecs(&mut self) -> Vec> { @@ -345,6 +521,40 @@ impl Solver { } } +/// Replaces variable `v` in `expr` with `v_expr` (an expression equal to `v`), +/// coalescing any resulting duplicate terms. Used by the elimination pre-pass. +fn substitute_var(expr: &mut LinearExpr, v: Var, v_expr: &LinearExpr) { + let Some(pos) = expr.coeffs.iter().position(|(_, var)| *var == v) else { + return; + }; + let (cv, _) = expr.coeffs.remove(pos); + for &(c, var) in &v_expr.coeffs { + if let Some(term) = expr.coeffs.iter_mut().find(|(_, ev)| *ev == var) { + term.0 += cv * c; + } else { + expr.coeffs.push((cv * c, var)); + } + } + expr.constant += cv * v_expr.constant; +} + +/// Merges duplicate-variable terms and drops near-zero coefficients in place, so a +/// constraint's `coeffs.len()` faithfully reflects its number of distinct variables +/// (e.g. a chain closing onto itself collapses `x + x` to a single `2x` term, and a +/// cancelling substitution collapses to an empty/contradiction constraint). +fn coalesce_terms(expr: &mut LinearExpr) { + let mut merged: Vec<(f64, Var)> = Vec::with_capacity(expr.coeffs.len()); + for &(c, v) in &expr.coeffs { + if let Some(term) = merged.iter_mut().find(|(_, mv)| *mv == v) { + term.0 += c; + } else { + merged.push((c, v)); + } + } + merged.retain(|(c, _)| relative_ne!(*c, 0., epsilon = EPSILON)); + expr.coeffs = merged; +} + struct ConstraintComponent { vars: IndexSet, constraints: Vec, @@ -503,4 +713,159 @@ mod tests { assert!(!solver.unsolved_vars.contains(&y)); assert!(solver.unsolved_vars.contains(&z)); } + + fn c(coeffs: Vec<(f64, Var)>, constant: f64) -> LinearExpr { + LinearExpr { coeffs, constant } + } + + /// A consistent ring of 2-variable constraints with no 1-variable starting point + /// (the minimal `bench_constraints` shape). The pre-pass must break the cycle by + /// substitution, then telescope to a 1-variable closure that grounds the chain. + #[test] + fn three_cycle_determined() { + let mut s = Solver::new(); + let a = s.new_var(); + let b = s.new_var(); + let d = s.new_var(); + s.constrain_eq0(c(vec![(1., a), (-1., b)], -5.)); // a - b = 5 + s.constrain_eq0(c(vec![(1., b), (-1., d)], -5.)); // b - c = 5 + s.constrain_eq0(c(vec![(1., a), (1., d)], -100.)); // a + c = 100 + s.solve(); + assert_relative_eq!(s.value_of(a).unwrap(), 55., epsilon = EPSILON); + assert_relative_eq!(s.value_of(b).unwrap(), 50., epsilon = EPSILON); + assert_relative_eq!(s.value_of(d).unwrap(), 45., epsilon = EPSILON); + assert!(s.inconsistent_constraints().is_empty()); + assert!(s.invalid_rounding().is_empty()); + } + + /// A chain pinned at one end resolves transitively (reverse-topological order). + #[test] + fn chain_telescopes() { + let mut s = Solver::new(); + let a = s.new_var(); + let b = s.new_var(); + let d = s.new_var(); + let e = s.new_var(); + s.constrain_eq0(c(vec![(1., a), (-1., b)], -1.)); // a - b = 1 + s.constrain_eq0(c(vec![(1., b), (-1., d)], -1.)); // b - c = 1 + s.constrain_eq0(c(vec![(1., d), (-1., e)], -1.)); // c - d = 1 + s.constrain_eq0(c(vec![(1., a)], -10.)); // a = 10 + s.solve(); + assert_relative_eq!(s.value_of(a).unwrap(), 10., epsilon = EPSILON); + assert_relative_eq!(s.value_of(b).unwrap(), 9., epsilon = EPSILON); + assert_relative_eq!(s.value_of(d).unwrap(), 8., epsilon = EPSILON); + assert_relative_eq!(s.value_of(e).unwrap(), 7., epsilon = EPSILON); + assert!(s.inconsistent_constraints().is_empty()); + } + + /// An under-determined pair: neither variable is pinned, and the row space still + /// reports one constrained direction (the pre-pass must re-materialize its + /// constraint after failing to ground the eliminated variable). + #[test] + fn underdetermined_pair_unsolved() { + let mut s = Solver::new(); + let a = s.new_var(); + let b = s.new_var(); + s.constrain_eq0(c(vec![(1., a), (-1., b)], 0.)); // a - b = 0 + s.solve(); + assert!(s.value_of(a).is_none()); + assert!(s.value_of(b).is_none()); + assert!(s.unsolved_vars().contains(&a)); + assert!(s.unsolved_vars().contains(&b)); + assert_eq!(s.rowspace_vecs().len(), 1); + } + + /// Two contradictory 2-variable constraints: substitution drives one to `0 = -5`. + #[test] + fn inconsistent_pair() { + let mut s = Solver::new(); + let a = s.new_var(); + let b = s.new_var(); + s.constrain_eq0(c(vec![(1., a), (-1., b)], 0.)); // a - b = 0 + s.constrain_eq0(c(vec![(1., a), (-1., b)], -5.)); // a - b = 5 + s.solve(); + assert!(!s.inconsistent_constraints().is_empty()); + } + + /// An over-constrained cycle whose differences sum to a nonzero constant. + #[test] + fn inconsistent_cycle() { + let mut s = Solver::new(); + let a = s.new_var(); + let b = s.new_var(); + let d = s.new_var(); + s.constrain_eq0(c(vec![(1., a), (-1., b)], -5.)); // a - b = 5 + s.constrain_eq0(c(vec![(1., b), (-1., d)], -5.)); // b - c = 5 + s.constrain_eq0(c(vec![(1., d), (-1., a)], -5.)); // c - a = 5 (loop sum = 15) + s.solve(); + assert!(!s.inconsistent_constraints().is_empty()); + } + + /// A duplicate constraint inside a coupled core is dropped as redundant (not flagged + /// inconsistent), and the remaining system still solves. + #[test] + fn redundant_in_core() { + let mut s = Solver::new(); + let a = s.new_var(); + let b = s.new_var(); + s.constrain_eq0(c(vec![(1., a), (-1., b)], 0.)); // a - b = 0 + s.constrain_eq0(c(vec![(1., a), (-1., b)], 0.)); // a - b = 0 (duplicate) + s.constrain_eq0(c(vec![(1., a), (1., b)], -10.)); // a + b = 10 + s.solve(); + assert_relative_eq!(s.value_of(a).unwrap(), 5., epsilon = EPSILON); + assert_relative_eq!(s.value_of(b).unwrap(), 5., epsilon = EPSILON); + assert!(s.inconsistent_constraints().is_empty()); + } + + /// A value reached only through elimination + resolution that lands off the 0.1 + /// grid is flagged in `invalid_rounding`. + #[test] + fn off_grid_cycle() { + let mut s = Solver::new(); + let a = s.new_var(); + let b = s.new_var(); + let d = s.new_var(); + s.constrain_eq0(c(vec![(1., a), (-1., b)], 0.)); // a - b = 0 + s.constrain_eq0(c(vec![(1., b), (-1., d)], 0.)); // b - c = 0 + s.constrain_eq0(c(vec![(1., a), (1., b), (1., d)], -1.)); // a + b + c = 1 => 1/3 each + s.solve(); + assert!(!s.invalid_rounding().is_empty()); + } + + /// Variables eliminated in an earlier `solve()` (before the closing constraint + /// exists) are resolved once a later constraint pins the chain. + #[test] + fn incremental_resolution() { + let mut s = Solver::new(); + let a = s.new_var(); + let b = s.new_var(); + let d = s.new_var(); + s.constrain_eq0(c(vec![(1., a), (-1., b)], -1.)); // a - b = 1 + s.constrain_eq0(c(vec![(1., b), (-1., d)], -1.)); // b - c = 1 + s.solve(); // under-determined so far + assert!(s.value_of(a).is_none()); + s.constrain_eq0(c(vec![(1., a)], -10.)); // a = 10 + s.solve(); + assert_relative_eq!(s.value_of(a).unwrap(), 10., epsilon = EPSILON); + assert_relative_eq!(s.value_of(b).unwrap(), 9., epsilon = EPSILON); + assert_relative_eq!(s.value_of(d).unwrap(), 8., epsilon = EPSILON); + } + + /// A fully-coupled 3x3 block has no size-<=2 pivot: the pre-pass is a no-op and the + /// dense SVD path solves it, exactly as before. + #[test] + fn dense_block_falls_back() { + let mut s = Solver::new(); + let a = s.new_var(); + let b = s.new_var(); + let d = s.new_var(); + s.constrain_eq0(c(vec![(1., a), (1., b), (1., d)], -6.)); // a + b + c = 6 + s.constrain_eq0(c(vec![(1., a), (2., b), (3., d)], -14.)); // a + 2b + 3c = 14 + s.constrain_eq0(c(vec![(1., a), (3., b), (6., d)], -25.)); // a + 3b + 6c = 25 + s.solve(); + assert_relative_eq!(s.value_of(a).unwrap(), 1., epsilon = EPSILON); + assert_relative_eq!(s.value_of(b).unwrap(), 2., epsilon = EPSILON); + assert_relative_eq!(s.value_of(d).unwrap(), 3., epsilon = EPSILON); + assert!(s.inconsistent_constraints().is_empty()); + } }