==> src/solver.c <==1/* kasuari's Solver, Row and Symbol, ported. Rows are small arrays of (symbol, coefficient) kept in
2 * insertion order; kasuari keeps them in hash maps, whose order only matters between equally
3 * good pivots, which leave the optimum (and so the layout) the same. */
4#include "tui/solver.h"
5
6#include <math.h>
7#include <stdio.h>
8#include <stdlib.h>
9#include <string.h>
10
11typedef enum { K_INVALID, K_EXTERNAL, K_SLACK, K_ERROR, K_DUMMY } kind;
12
13typedef struct {
14 uint32_t id;
15 uint8_t kind;
16} symbol;
17
18static const symbol INVALID = {0, K_INVALID};
19
20static bool sym_eq(symbol a, symbol b) { return a.id == b.id && a.kind == b.kind; }
21
22typedef struct {
23 symbol sym;
24 double coeff;
25} cell;
26
27typedef struct {
28 cell *cells;
29 size_t n, cap;
30 double constant;
31} row;
32
33typedef struct {
34 symbol sym; /* the basic symbol */
35 row *row;
36} basic;
37
38typedef struct {
39 uint32_t var;
40 symbol sym;
41 double value;
42 size_t count;
43} vardata;
44
45struct tui_solver {
46 vardata *vars;
47 size_t nvars, capvars;
48 basic *rows;
49 size_t nrows, caprows;
50 symbol *infeasible;
51 size_t ninf, capinf;
52 row objective;
53 row *artificial;
54 uint32_t id_tick;
55};
56
57static void *grow(void *p, size_t *cap, size_t need, size_t size) {
58 if (need <= *cap) return p;
59 size_t c = *cap ? *cap : 8;
60 while (c < need) c *= 2;
61 p = realloc(p, c * size);
62 if (!p) abort();
63 *cap = c;
64 return p;
65}
66
67/* ---- expressions ---- */
68
69static tui_expr expr_append(tui_expr a, tui_expr b, double sign) {
70 for (int i = 0; i < b.n; i++) {
71 if (a.n == TUI_EXPR_TERMS) {
72 fputs("tui-solver: too many terms in an expression\n", stderr);
73 abort();
74 }
75 a.terms[a.n++] = (tui_term){b.terms[i].var, sign * b.terms[i].coeff};
76 }
77 a.constant += sign * b.constant;
78 return a;
79}
80
81tui_expr tui_expr_sub(tui_expr a, tui_expr b) { return expr_append(a, b, -1.0); }
82tui_expr tui_expr_add(tui_expr a, tui_expr b) { return expr_append(a, b, 1.0); }
83tui_expr tui_expr_div(tui_expr a, double k) {
84 a.constant /= k;
85 for (int i = 0; i < a.n; i++) a.terms[i].coeff /= k;
86 return a;
87}
88tui_expr tui_expr_scale(tui_expr a, double k) {
89 a.constant *= k;
90 for (int i = 0; i < a.n; i++) a.terms[i].coeff *= k;
91 return a;
92}
93
94/* ---- rows ---- */
95
96static bool near_zero(double v) { return v < 0 ? -v < 1e-8 : v < 1e-8; }
97
98static row *row_new(double constant) {
99 row *r = calloc(1, sizeof *r);
100 if (!r) abort();
101 r->constant = constant;
102 return r;
103}
104
105static void row_free(row *r) {
106 if (!r) return;
107 free(r->cells);
108 free(r);
109}
110
111static row *row_clone(const row *r) {
112 row *c = row_new(r->constant);
113 c->cells = grow(NULL, &c->cap, r->n ? r->n : 1, sizeof *c->cells);
114 memcpy(c->cells, r->cells, r->n * sizeof *c->cells);
115 c->n = r->n;
116 return c;
117}
118
119static long row_find(const row *r, symbol s) {
120 for (size_t i = 0; i < r->n; i++)
121 if (sym_eq(r->cells[i].sym, s)) return (long)i;
122 return -1;
123}
124
125static void row_remove_at(row *r, size_t i) {
126 memmove(&r->cells[i], &r->cells[i + 1], (r->n - i - 1) * sizeof *r->cells);
127 r->n--;
128}
129
130static void row_insert_symbol(row *r, symbol s, double coeff) {
131 long i = row_find(r, s);
132 if (i < 0) {
133 if (!near_zero(coeff)) {
134 r->cells = grow(r->cells, &r->cap, r->n + 1, sizeof *r->cells);
135 r->cells[r->n++] = (cell){s, coeff};
136 }
137 return;
138 }
139 r->cells[i].coeff += coeff;
140 if (near_zero(r->cells[i].coeff)) row_remove_at(r, (size_t)i);
141}
142
143static bool row_insert_row(row *r, const row *other, double coeff) {
144 double diff = other->constant * coeff;
145 r->constant += diff;
146 for (size_t i = 0; i < other->n; i++) row_insert_symbol(r, other->cells[i].sym, other->cells[i].coeff * coeff);
147 return diff != 0.0;
148}
149
150static void row_remove(row *r, symbol s) {
151 long i = row_find(r, s);
152 if (i >= 0) row_remove_at(r, (size_t)i);
153}
154
155static void row_reverse_sign(row *r) {
156 r->constant = -r->constant;
157 for (size_t i = 0; i < r->n; i++) r->cells[i].coeff = -r->cells[i].coeff;
158}
159
160static void row_solve_for_symbol(row *r, symbol s) {
161 long i = row_find(r, s);
162 double c = r->cells[i].coeff;
163 row_remove_at(r, (size_t)i);
164 double coeff = -1.0 / c;
165 r->constant *= coeff;
166 for (size_t k = 0; k < r->n; k++) r->cells[k].coeff *= coeff;
167}
168
169static void row_solve_for_symbols(row *r, symbol lhs, symbol rhs) {
170 row_insert_symbol(r, lhs, -1.0);
171 row_solve_for_symbol(r, rhs);
172}
173
174static double row_coefficient(const row *r, symbol s) {
175 long i = row_find(r, s);
176 return i < 0 ? 0.0 : r->cells[i].coeff;
177}
178
179static bool row_substitute(row *r, symbol s, const row *with) {
180 long i = row_find(r, s);
181 if (i < 0) return false;
182 double coeff = r->cells[i].coeff;
183 row_remove_at(r, (size_t)i);
184 return row_insert_row(r, with, coeff);
185}
186
187/* ---- the solver ---- */
188
189tui_solver *tui_solver_new(void) {
190 tui_solver *s = calloc(1, sizeof *s);
191 if (!s) abort();
192 s->id_tick = 1;
193 return s;
194}
195
196void tui_solver_free(tui_solver *s) {
197 if (!s) return;
198 for (size_t i = 0; i < s->nrows; i++) row_free(s->rows[i].row);
199 free(s->rows);
200 free(s->vars);
201 free(s->infeasible);
202 free(s->objective.cells);
203 row_free(s->artificial);
204 free(s);
205}
206
207static long rows_find(const tui_solver *s, symbol sym) {
208 for (size_t i = 0; i < s->nrows; i++)
209 if (sym_eq(s->rows[i].sym, sym)) return (long)i;
210 return -1;
211}
212
213static row *rows_get(const tui_solver *s, symbol sym) {
214 long i = rows_find(s, sym);
215 return i < 0 ? NULL : s->rows[i].row;
216}
217
218static void rows_insert(tui_solver *s, symbol sym, row *r) {
219 long i = rows_find(s, sym);
220 if (i >= 0) {
221 row_free(s->rows[i].row);
222 s->rows[i].row = r;
223 return;
224 }
225 s->rows = grow(s->rows, &s->caprows, s->nrows + 1, sizeof *s->rows);
226 s->rows[s->nrows++] = (basic){sym, r};
227}
228
229static row *rows_take(tui_solver *s, symbol sym) {
230 long i = rows_find(s, sym);
231 if (i < 0) return NULL;
232 row *r = s->rows[i].row;
233 memmove(&s->rows[i], &s->rows[i + 1], (s->nrows - (size_t)i - 1) * sizeof *s->rows);
234 s->nrows--;
235 return r;
236}
237
238static void push_infeasible(tui_solver *s, symbol sym) {
239 s->infeasible = grow(s->infeasible, &s->capinf, s->ninf + 1, sizeof *s->infeasible);
240 s->infeasible[s->ninf++] = sym;
241}
242
243static symbol var_symbol(tui_solver *s, uint32_t var) {
244 for (size_t i = 0; i < s->nvars; i++)
245 if (s->vars[i].var == var) {
246 s->vars[i].count++;
247 return s->vars[i].sym;
248 }
249 s->vars = grow(s->vars, &s->capvars, s->nvars + 1, sizeof *s->vars);
250 symbol sym = {s->id_tick++, K_EXTERNAL};
251 s->vars[s->nvars++] = (vardata){var, sym, NAN, 1};
252 return sym;
253}
254
255typedef struct {
256 symbol marker, other;
257} tag;
258
259static row *create_row(tui_solver *s, const tui_expr *e, tui_relation op, double strength, tag *t) {
260 row *r = row_new(e->constant);
261 for (int i = 0; i < e->n; i++) {
262 if (near_zero(e->terms[i].coeff)) continue;
263 symbol sym = var_symbol(s, e->terms[i].var);
264 row *other = rows_get(s, sym);
265 if (other) row_insert_row(r, other, e->terms[i].coeff);
266 else row_insert_symbol(r, sym, e->terms[i].coeff);
267 }
268 if (op == TUI_LE || op == TUI_GE) {
269 double coeff = op == TUI_LE ? 1.0 : -1.0;
270 symbol slack = {s->id_tick++, K_SLACK};
271 row_insert_symbol(r, slack, coeff);
272 if (strength < TUI_STRENGTH_REQUIRED) {
273 symbol error = {s->id_tick++, K_ERROR};
274 row_insert_symbol(r, error, -coeff);
275 row_insert_symbol(&s->objective, error, strength);
276 *t = (tag){slack, error};
277 } else {
278 *t = (tag){slack, INVALID};
279 }
280 } else if (strength < TUI_STRENGTH_REQUIRED) {
281 symbol errplus = {s->id_tick++, K_ERROR};
282 symbol errminus = {s->id_tick++, K_ERROR};
283 row_insert_symbol(r, errplus, -1.0);
284 row_insert_symbol(r, errminus, 1.0);
285 row_insert_symbol(&s->objective, errplus, strength);
286 row_insert_symbol(&s->objective, errminus, strength);
287 *t = (tag){errplus, errminus};
288 } else {
289 symbol dummy = {s->id_tick++, K_DUMMY};
290 row_insert_symbol(r, dummy, 1.0);
291 *t = (tag){dummy, INVALID};
292 }
293 if (r->constant < 0.0) row_reverse_sign(r);
294 return r;
295}
296
297static symbol choose_subject(const row *r, const tag *t) {
298 for (size_t i = 0; i < r->n; i++)
299 if (r->cells[i].sym.kind == K_EXTERNAL) return r->cells[i].sym;
300 if ((t->marker.kind == K_SLACK || t->marker.kind == K_ERROR) && row_coefficient(r, t->marker) < 0.0) return t->marker;
301 if ((t->other.kind == K_SLACK || t->other.kind == K_ERROR) && row_coefficient(r, t->other) < 0.0) return t->other;
302 return INVALID;
303}
304
305static bool all_dummies(const row *r) {
306 for (size_t i = 0; i < r->n; i++)
307 if (r->cells[i].sym.kind != K_DUMMY) return false;
308 return true;
309}
310
311static void substitute(tui_solver *s, symbol sym, const row *r) {
312 for (size_t i = 0; i < s->nrows; i++) {
313 row_substitute(s->rows[i].row, sym, r);
314 if (s->rows[i].sym.kind != K_EXTERNAL && s->rows[i].row->constant < 0.0) push_infeasible(s, s->rows[i].sym);
315 }
316 row_substitute(&s->objective, sym, r);
317 if (s->artificial) row_substitute(s->artificial, sym, r);
318}
319
320static symbol entering_symbol(const row *objective) {
321 for (size_t i = 0; i < objective->n; i++)
322 if (objective->cells[i].sym.kind != K_DUMMY && objective->cells[i].coeff < 0.0) return objective->cells[i].sym;
323 return INVALID;
324}
325
326static bool leaving_row(tui_solver *s, symbol entering, symbol *leaving, row **r) {
327 double ratio = INFINITY;
328 long found = -1;
329 for (size_t i = 0; i < s->nrows; i++) {
330 if (s->rows[i].sym.kind == K_EXTERNAL) continue;
331 double temp = row_coefficient(s->rows[i].row, entering);
332 if (temp < 0.0) {
333 double tr = -s->rows[i].row->constant / temp;
334 if (tr < ratio) {
335 ratio = tr;
336 found = (long)i;
337 }
338 }
339 }
340 if (found < 0) return false;
341 *leaving = s->rows[found].sym;
342 *r = rows_take(s, *leaving);
343 return true;
344}
345
346static bool optimize(tui_solver *s, bool artificial) {
347 for (;;) {
348 symbol entering = entering_symbol(artificial ? s->artificial : &s->objective);
349 if (entering.kind == K_INVALID) return true;
350 symbol leaving;
351 row *r;
352 if (!leaving_row(s, entering, &leaving, &r)) return false; /* objective unbounded */
353 row_solve_for_symbols(r, leaving, entering);
354 substitute(s, entering, r);
355 rows_insert(s, entering, r);
356 }
357}
358
359static symbol any_pivotable_symbol(const row *r) {
360 for (size_t i = 0; i < r->n; i++)
361 if (r->cells[i].sym.kind == K_SLACK || r->cells[i].sym.kind == K_ERROR) return r->cells[i].sym;
362 return INVALID;
363}
364
365static int add_with_artificial_variable(tui_solver *s, const row *r) {
366 symbol art = {s->id_tick++, K_SLACK};
367 rows_insert(s, art, row_clone(r));
368 s->artificial = row_clone(r);
369 if (!optimize(s, true)) return -1;
370 bool success = near_zero(s->artificial->constant);
371 row_free(s->artificial);
372 s->artificial = NULL;
373 row *ar = rows_take(s, art);
374 if (ar) {
375 if (ar->n == 0) {
376 row_free(ar);
377 return success;
378 }
379 symbol entering = any_pivotable_symbol(ar);
380 if (entering.kind == K_INVALID) {
381 row_free(ar);
382 return 0;
383 }
384 row_solve_for_symbols(ar, art, entering);
385 substitute(s, entering, ar);
386 rows_insert(s, entering, ar);
387 }
388 for (size_t i = 0; i < s->nrows; i++) row_remove(s->rows[i].row, art);
389 row_remove(&s->objective, art);
390 return success;
391}
392
393bool tui_solver_add(tui_solver *s, tui_expr lhs, tui_relation op, double strength, tui_expr rhs) {
394 tui_expr e = tui_expr_sub(lhs, rhs);
395 tag t;
396 row *r = create_row(s, &e, op, strength, &t);
397 symbol subject = choose_subject(r, &t);
398 if (subject.kind == K_INVALID && all_dummies(r)) {
399 if (!near_zero(r->constant)) {
400 row_free(r);
401 return false;
402 }
403 subject = t.marker;
404 }
405 if (subject.kind == K_INVALID) {
406 int ok = add_with_artificial_variable(s, r);
407 row_free(r);
408 if (ok != 1) return false;
409 } else {
410 row_solve_for_symbol(r, subject);
411 substitute(s, subject, r);
412 rows_insert(s, subject, r);
413 }
414 return optimize(s, false);
415}
416
417double tui_solver_value(const tui_solver *s, uint32_t var) {
418 for (size_t i = 0; i < s->nvars; i++)
419 if (s->vars[i].var == var) {
420 row *r = rows_get(s, s->vars[i].sym);
421 return r ? r->constant : 0.0;
422 }
423 return 0.0;
424}