?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Compilation: fix up tools includes.
[gnumeric.git] / src / tools / gnm-solver.c
blob6e3e9c3ed6b667839a5d9dc59d3929a47664db60
1 #include <gnumeric-config.h>
2 #include "gnumeric.h"
3 #include "value.h"
4 #include "cell.h"
5 #include "expr.h"
6 #include "expr-deriv.h"
7 #include "sheet.h"
8 #include "workbook.h"
9 #include "rangefunc.h"
10 #include "ranges.h"
11 #include "gutils.h"
12 #include "mathfunc.h"
13 #include "gnm-solver.h"
14 #include "workbook-view.h"
15 #include "workbook-control.h"
16 #include "application.h"
17 #include "gnm-marshalers.h"
18 #include <tools/dao.h>
19 #include "gui-util.h"
20 #include "gnm-i18n.h"
21 #include <gsf/gsf-impl-utils.h>
22 #include <gsf/gsf-output-stdio.h>
23 #include <glib/gi18n-lib.h>
24 #include <glib/gstdio.h>
25 #include <unistd.h>
26 #include <signal.h>
27 #include <string.h>
28 #ifdef HAVE_SYS_WAIT_H
29 #include <sys/wait.h>
30 #endif
31
32 #ifdef G_OS_WIN32
33 #include <windows.h>
34 #ifndef WIFEXITED
35 #define WIFEXITED(x) ((x) != STILL_ACTIVE)
36 #endif
37 #ifndef WEXITSTATUS
38 #define WEXITSTATUS(x) (x)
39 #endif
40 #endif
41
42 /* ------------------------------------------------------------------------- */
43
44 gboolean
45 gnm_solver_debug (void)
46 {
47 static int debug = -1;
48 if (debug == -1)
49 debug = gnm_debug_flag ("solver");
50 return debug;
51 }
52
53 static char *
54 gnm_solver_cell_name (GnmCell const *cell, Sheet *origin)
55 {
56 GnmConventionsOut out;
57 GnmCellRef cr;
58 GnmParsePos pp;
59
60 gnm_cellref_init (&cr, cell->base.sheet,
61 cell->pos.col, cell->pos.row,
62 TRUE);
63 out.accum = g_string_new (NULL);
64 out.pp = parse_pos_init_sheet (&pp, origin);
65 out.convs = origin->convs;
66 cellref_as_string (&out, &cr, cell->base.sheet == origin);
67 return g_string_free (out.accum, FALSE);
68 }
69
70 /* ------------------------------------------------------------------------- */
71
72 GType
73 gnm_solver_status_get_type (void)
74 {
75 static GType etype = 0;
76 if (etype == 0) {
77 static GEnumValue const values[] = {
78 { GNM_SOLVER_STATUS_READY,
79 "GNM_SOLVER_STATUS_READY",
80 "ready"
81 },
82 { GNM_SOLVER_STATUS_PREPARING,
83 "GNM_SOLVER_STATUS_PREPARING",
84 "preparing"
85 },
86 { GNM_SOLVER_STATUS_PREPARED,
87 "GNM_SOLVER_STATUS_PREPARED",
88 "prepared"
89 },
90 { GNM_SOLVER_STATUS_RUNNING,
91 "GNM_SOLVER_STATUS_RUNNING",
92 "running"
93 },
94 { GNM_SOLVER_STATUS_DONE,
95 "GNM_SOLVER_STATUS_DONE",
96 "done"
97 },
98 { GNM_SOLVER_STATUS_ERROR,
99 "GNM_SOLVER_STATUS_ERROR",
100 "error"
101 },
102 { GNM_SOLVER_STATUS_CANCELLED,
103 "GNM_SOLVER_STATUS_CANCELLED",
104 "cancelled"
105 },
106 { 0, NULL, NULL }
107 };
108 etype = g_enum_register_static ("GnmSolverStatus", values);
109 }
110 return etype;
111 }
112
113 GType
114 gnm_solver_problem_type_get_type (void)
115 {
116 static GType etype = 0;
117 if (etype == 0) {
118 static GEnumValue const values[] = {
119 { GNM_SOLVER_MINIMIZE,
120 "GNM_SOLVER_MINIMIZE",
121 "minimize"
122 },
123 { GNM_SOLVER_MAXIMIZE,
124 "GNM_SOLVER_MAXIMIZE",
125 "maximize"
126 },
127 { 0, NULL, NULL }
128 };
129 etype = g_enum_register_static ("GnmSolverProblemType", values);
130 }
131 return etype;
132 }
133
134 /* ------------------------------------------------------------------------- */
135
136 static GObjectClass *gnm_solver_parent_class;
137
138 GnmSolverConstraint *
139 gnm_solver_constraint_new (Sheet *sheet)
140 {
141 GnmSolverConstraint *res = g_new0 (GnmSolverConstraint, 1);
142 dependent_managed_init (&res->lhs, sheet);
143 dependent_managed_init (&res->rhs, sheet);
144 return res;
145 }
146
147 void
148 gnm_solver_constraint_free (GnmSolverConstraint *c)
149 {
150 gnm_solver_constraint_set_lhs (c, NULL);
151 gnm_solver_constraint_set_rhs (c, NULL);
152 g_free (c);
153 }
154
155 GnmSolverConstraint *
156 gnm_solver_constraint_dup (GnmSolverConstraint *c, Sheet *sheet)
157 {
158 GnmSolverConstraint *res = gnm_solver_constraint_new (sheet);
159 res->type = c->type;
160 dependent_managed_set_expr (&res->lhs, c->lhs.texpr);
161 dependent_managed_set_expr (&res->rhs, c->rhs.texpr);
162 return res;
163 }
164
165 static GnmSolverConstraint *
166 gnm_solver_constraint_dup1 (GnmSolverConstraint *c)
167 {
168 return gnm_solver_constraint_dup (c, c->lhs.sheet);
169 }
170
171 GType
172 gnm_solver_constraint_get_type (void)
173 {
174 static GType t = 0;
175
176 if (t == 0)
177 t = g_boxed_type_register_static ("GnmSolverConstraint",
178 (GBoxedCopyFunc)gnm_solver_constraint_dup1,
179 (GBoxedFreeFunc)gnm_solver_constraint_free);
180 return t;
181 }
182
183 gboolean
184 gnm_solver_constraint_equal (GnmSolverConstraint const *a,
185 GnmSolverConstraint const *b)
186 {
187 return (a->type == b->type &&
188 gnm_expr_top_equal (a->lhs.texpr, b->lhs.texpr) &&
189 (!gnm_solver_constraint_has_rhs (a) ||
190 gnm_expr_top_equal (a->rhs.texpr, b->rhs.texpr)));
191 }
192
193 gboolean
194 gnm_solver_constraint_has_rhs (GnmSolverConstraint const *c)
195 {
196 g_return_val_if_fail (c != NULL, FALSE);
197
198 switch (c->type) {
199 case GNM_SOLVER_LE:
200 case GNM_SOLVER_GE:
201 case GNM_SOLVER_EQ:
202 return TRUE;
203 case GNM_SOLVER_INTEGER:
204 case GNM_SOLVER_BOOLEAN:
205 default:
206 return FALSE;
207 }
208 }
209
210 gboolean
211 gnm_solver_constraint_valid (GnmSolverConstraint const *c,
212 GnmSolverParameters const *sp)
213 {
214 GnmValue const *lhs;
215
216 g_return_val_if_fail (c != NULL, FALSE);
217
218 lhs = gnm_solver_constraint_get_lhs (c);
219 if (lhs == NULL || !VALUE_IS_CELLRANGE (lhs))
220 return FALSE;
221
222 if (gnm_solver_constraint_has_rhs (c)) {
223 GnmValue const *rhs = gnm_solver_constraint_get_rhs (c);
224 if (rhs == NULL)
225 return FALSE;
226 if (VALUE_IS_CELLRANGE (rhs)) {
227 GnmSheetRange srl, srr;
228
229 gnm_sheet_range_from_value (&srl, lhs);
230 gnm_sheet_range_from_value (&srr, rhs);
231
232 if (range_width (&srl.range) != range_width (&srr.range) ||
233 range_height (&srl.range) != range_height (&srr.range))
234 return FALSE;
235 } else if (VALUE_IS_FLOAT (rhs)) {
236 /* Nothing */
237 } else
238 return FALSE;
239 }
240
241 switch (c->type) {
242 case GNM_SOLVER_INTEGER:
243 case GNM_SOLVER_BOOLEAN: {
244 GnmValue const *vinput = gnm_solver_param_get_input (sp);
245 GnmSheetRange sr_input, sr_c;
246
247 if (!vinput)
248 break; /* No need to blame constraint. */
249
250 gnm_sheet_range_from_value (&sr_input, vinput);
251 gnm_sheet_range_from_value (&sr_c, lhs);
252
253 if (eval_sheet (sr_input.sheet, sp->sheet) !=
254 eval_sheet (sr_c.sheet, sp->sheet) ||
255 !range_contained (&sr_c.range, &sr_input.range))
256 return FALSE;
257 break;
258 }
259
260 default:
261 break;
262 }
263
264 return TRUE;
265 }
266
267 /**
268 * gnm_solver_constraint_get_lhs:
269 * @c: GnmSolverConstraint
270 *
271 * Returns: (transfer none) (nullable): Get left-hand side of constraint @c.
272 */
273 GnmValue const *
274 gnm_solver_constraint_get_lhs (GnmSolverConstraint const *c)
275 {
276 GnmExprTop const *texpr = c->lhs.texpr;
277 return texpr ? gnm_expr_top_get_constant (texpr) : NULL;
278 }
279
280 /**
281 * gnm_solver_constraint_set_lhs:
282 * @c: GnmSolverConstraint
283 * @v: (transfer full) (nullable): new left-hand side
284 */
285 void
286 gnm_solver_constraint_set_lhs (GnmSolverConstraint *c, GnmValue *v)
287 {
288 GnmExprTop const *texpr = v ? gnm_expr_top_new_constant (v) : NULL;
289 dependent_managed_set_expr (&c->lhs, texpr);
290 if (texpr) gnm_expr_top_unref (texpr);
291 }
292
293 /**
294 * gnm_solver_constraint_get_rhs:
295 * @c: GnmSolverConstraint
296 *
297 * Returns: (transfer none) (nullable): Get right-hand side of constraint @c.
298 */
299 GnmValue const *
300 gnm_solver_constraint_get_rhs (GnmSolverConstraint const *c)
301 {
302 GnmExprTop const *texpr = c->rhs.texpr;
303 return texpr ? gnm_expr_top_get_constant (texpr) : NULL;
304 }
305
306 /**
307 * gnm_solver_constraint_set_rhs:
308 * @c: GnmSolverConstraint
309 * @v: (transfer full) (nullable): new right-hand side
310 */
311 void
312 gnm_solver_constraint_set_rhs (GnmSolverConstraint *c, GnmValue *v)
313 {
314 GnmExprTop const *texpr = v ? gnm_expr_top_new_constant (v) : NULL;
315 dependent_managed_set_expr (&c->rhs, texpr);
316 if (texpr) gnm_expr_top_unref (texpr);
317 }
318
319 /**
320 * gnm_solver_constraint_get_part:
321 * @c: GnmSolverConstraint
322 * @sp: GnmSolverParameters
323 * @i: part index
324 * @lhs: (optional) (out): #GnmCell of left-hand side
325 * @cl: (optional) (out): constant value of left-hand side
326 * @rhs: (optional) (out): #GnmCell of right-hand side
327 * @cr: (optional) (out): constant value of left-hand side
328 *
329 * This splits @c into parts and returns information about the @i'th part.
330 * There will be multiple parts when the left-hand side is a cell range.
331 *
332 * Returns: %TRUE if the @i'th part exists.
333 */
334 gboolean
335 gnm_solver_constraint_get_part (GnmSolverConstraint const *c,
336 GnmSolverParameters const *sp, int i,
337 GnmCell **lhs, gnm_float *cl,
338 GnmCell **rhs, gnm_float *cr)
339 {
340 GnmSheetRange sr;
341 int h, w, dx, dy;
342 GnmValue const *vl, *vr;
343
344 if (cl) *cl = 0;
345 if (cr) *cr = 0;
346 if (lhs) *lhs = NULL;
347 if (rhs) *rhs = NULL;
348
349 if (!gnm_solver_constraint_valid (c, sp))
350 return FALSE;
351
352 vl = gnm_solver_constraint_get_lhs (c);
353 vr = gnm_solver_constraint_get_rhs (c);
354
355 gnm_sheet_range_from_value (&sr, vl);
356 w = range_width (&sr.range);
357 h = range_height (&sr.range);
358
359 dy = i / w;
360 dx = i % w;
361 if (dy >= h)
362 return FALSE;
363
364 if (lhs)
365 *lhs = sheet_cell_get (eval_sheet (sr.sheet, sp->sheet),
366 sr.range.start.col + dx,
367 sr.range.start.row + dy);
368
369 if (!gnm_solver_constraint_has_rhs (c)) {
370 /* Nothing */
371 } else if (VALUE_IS_FLOAT (vr)) {
372 if (cr)
373 *cr = value_get_as_float (vr);
374 } else {
375 gnm_sheet_range_from_value (&sr, vr);
376 if (rhs)
377 *rhs = sheet_cell_get (eval_sheet (sr.sheet, sp->sheet),
378 sr.range.start.col + dx,
379 sr.range.start.row + dy);
380 }
381
382 return TRUE;
383 }
384
385 void
386 gnm_solver_constraint_set_old (GnmSolverConstraint *c,
387 GnmSolverConstraintType type,
388 int lhs_col, int lhs_row,
389 int rhs_col, int rhs_row,
390 int cols, int rows)
391 {
392 GnmRange r;
393
394 c->type = type;
395
396 range_init (&r,
397 lhs_col, lhs_row,
398 lhs_col + (cols - 1), lhs_row + (rows - 1));
399 gnm_solver_constraint_set_lhs
400 (c, value_new_cellrange_r (NULL, &r));
401
402 if (gnm_solver_constraint_has_rhs (c)) {
403 range_init (&r,
404 rhs_col, rhs_row,
405 rhs_col + (cols - 1), rhs_row + (rows - 1));
406 gnm_solver_constraint_set_rhs
407 (c, value_new_cellrange_r (NULL, &r));
408 } else
409 gnm_solver_constraint_set_rhs (c, NULL);
410 }
411
412 void
413 gnm_solver_constraint_side_as_str (GnmSolverConstraint const *c,
414 Sheet const *sheet,
415 GString *buf, gboolean lhs)
416 {
417 GnmExprTop const *texpr;
418
419 texpr = lhs ? c->lhs.texpr : c->rhs.texpr;
420 if (texpr) {
421 GnmConventionsOut out;
422 GnmParsePos pp;
423
424 out.accum = buf;
425 out.pp = parse_pos_init_sheet (&pp, sheet);
426 out.convs = sheet->convs;
427 gnm_expr_top_as_gstring (texpr, &out);
428 } else
429 g_string_append (buf,
430 value_error_name (GNM_ERROR_REF,
431 sheet->convs->output.translated));
432 }
433
434 char *
435 gnm_solver_constraint_as_str (GnmSolverConstraint const *c, Sheet *sheet)
436 {
437 const char * const type_str[] = {
438 "\xe2\x89\xa4" /* "<=" */,
439 "\xe2\x89\xa5" /* ">=" */,
440 "=",
441 N_("Int"),
442 N_("Bool")
443 };
444 const char *type = type_str[c->type];
445 gboolean translate = (c->type >= GNM_SOLVER_INTEGER);
446 GString *buf = g_string_new (NULL);
447
448 gnm_solver_constraint_side_as_str (c, sheet, buf, TRUE);
449 g_string_append_c (buf, ' ');
450 g_string_append (buf, translate ? _(type) : type);
451 if (gnm_solver_constraint_has_rhs (c)) {
452 g_string_append_c (buf, ' ');
453 gnm_solver_constraint_side_as_str (c, sheet, buf, FALSE);
454 }
455
456 return g_string_free (buf, FALSE);
457 }
458
459 char *
460 gnm_solver_constraint_part_as_str (GnmSolverConstraint const *c, int i,
461 GnmSolverParameters *sp)
462 {
463 const char * const type_str[] = {
464 "\xe2\x89\xa4" /* "<=" */,
465 "\xe2\x89\xa5" /* ">=" */,
466 "=",
467 N_("Int"),
468 N_("Bool")
469 };
470 const char *type = type_str[c->type];
471 gboolean translate = (c->type >= GNM_SOLVER_INTEGER);
472 GString *buf;
473 gnm_float cl, cr;
474 GnmCell *lhs, *rhs;
475
476 if (!gnm_solver_constraint_get_part (c, sp, i, &lhs, &cl, &rhs, &cr))
477 return NULL;
478
479 buf = g_string_new (NULL);
480
481 g_string_append (buf, cell_name (lhs));
482 g_string_append_c (buf, ' ');
483 g_string_append (buf, translate ? _(type) : type);
484 if (gnm_solver_constraint_has_rhs (c)) {
485 g_string_append_c (buf, ' ');
486 g_string_append (buf, cell_name (rhs));
487 }
488
489 return g_string_free (buf, FALSE);
490 }
491
492 /* ------------------------------------------------------------------------- */
493
494 enum {
495 SOLP_PROP_0,
496 SOLP_PROP_SHEET,
497 SOLP_PROP_PROBLEM_TYPE
498 };
499
500 static GObjectClass *gnm_solver_param_parent_class;
501
502 GnmSolverParameters *
503 gnm_solver_param_new (Sheet *sheet)
504 {
505 return g_object_new (GNM_SOLVER_PARAMETERS_TYPE,
506 "sheet", sheet,
507 NULL);
508 }
509
510 /**
511 * gnm_solver_param_dup:
512 * @src: #GnmSolverParameters
513 * @new_sheet: #Sheet
514 *
515 * Returns: (transfer full): duplicate @src, but for @new_sheet.
516 */
517 GnmSolverParameters *
518 gnm_solver_param_dup (GnmSolverParameters *src, Sheet *new_sheet)
519 {
520 GnmSolverParameters *dst = gnm_solver_param_new (new_sheet);
521 GSList *l;
522
523 dst->problem_type = src->problem_type;
524 dependent_managed_set_expr (&dst->target, src->target.texpr);
525 dependent_managed_set_expr (&dst->input, src->input.texpr);
526
527 g_free (dst->options.scenario_name);
528 dst->options = src->options;
529 dst->options.algorithm = NULL;
530 dst->options.scenario_name = g_strdup (src->options.scenario_name);
531 gnm_solver_param_set_algorithm (dst, src->options.algorithm);
532
533 /* Copy the constraints */
534 for (l = src->constraints; l; l = l->next) {
535 GnmSolverConstraint *old = l->data;
536 GnmSolverConstraint *new =
537 gnm_solver_constraint_dup (old, new_sheet);
538
539 dst->constraints = g_slist_prepend (dst->constraints, new);
540 }
541 dst->constraints = g_slist_reverse (dst->constraints);
542
543 return dst;
544 }
545
546 gboolean
547 gnm_solver_param_equal (GnmSolverParameters const *a,
548 GnmSolverParameters const *b)
549 {
550 GSList *la, *lb;
551
552 if (a->sheet != b->sheet ||
553 a->problem_type != b->problem_type ||
554 !gnm_expr_top_equal (a->target.texpr, b->target.texpr) ||
555 !gnm_expr_top_equal (a->input.texpr, b->input.texpr) ||
556 a->options.max_time_sec != b->options.max_time_sec ||
557 a->options.max_iter != b->options.max_iter ||
558 a->options.algorithm != b->options.algorithm ||
559 a->options.model_type != b->options.model_type ||
560 a->options.assume_non_negative != b->options.assume_non_negative ||
561 a->options.assume_discrete != b->options.assume_discrete ||
562 a->options.automatic_scaling != b->options.automatic_scaling ||
563 a->options.program_report != b->options.program_report ||
564 a->options.sensitivity_report != b->options.sensitivity_report ||
565 a->options.add_scenario != b->options.add_scenario ||
566 strcmp (a->options.scenario_name, b->options.scenario_name) ||
567 a->options.gradient_order != b->options.gradient_order)
568 return FALSE;
569
570 for (la = a->constraints, lb = b->constraints;
571 la && lb;
572 la = la->next, lb = lb->next) {
573 GnmSolverConstraint *ca = la->data;
574 GnmSolverConstraint *cb = lb->data;
575 if (!gnm_solver_constraint_equal (ca, cb))
576 return FALSE;
577 }
578 return la == lb;
579 }
580
581 /**
582 * gnm_solver_param_get_input:
583 * @sp: #GnmSolverParameters
584 *
585 * Returns: (transfer none) (nullable): the input cell area.
586 */
587 GnmValue const *
588 gnm_solver_param_get_input (GnmSolverParameters const *sp)
589 {
590 return sp->input.texpr
591 ? gnm_expr_top_get_constant (sp->input.texpr)
592 : NULL;
593 }
594
595 /**
596 * gnm_solver_param_set_input:
597 * @sp: #GnmSolverParameters
598 * @v: (transfer full) (nullable): new input area
599 */
600 void
601 gnm_solver_param_set_input (GnmSolverParameters *sp, GnmValue *v)
602 {
603 GnmExprTop const *texpr = v ? gnm_expr_top_new_constant (v) : NULL;
604 dependent_managed_set_expr (&sp->input, texpr);
605 if (texpr) gnm_expr_top_unref (texpr);
606 }
607
608 static GnmValue *
609 cb_grab_cells (GnmCellIter const *iter, gpointer user)
610 {
611 GPtrArray *input_cells = user;
612 GnmCell *cell;
613
614 if (NULL == (cell = iter->cell))
615 cell = sheet_cell_create (iter->pp.sheet,
616 iter->pp.eval.col, iter->pp.eval.row);
617 g_ptr_array_add (input_cells, cell);
618 return NULL;
619 }
620
621 /**
622 * gnm_solver_param_get_input_cells:
623 * @sp: #GnmSolverParameters
624 *
625 * Returns: (element-type GnmCell) (transfer container):
626 */
627 GPtrArray *
628 gnm_solver_param_get_input_cells (GnmSolverParameters const *sp)
629 {
630 GnmValue const *vr = gnm_solver_param_get_input (sp);
631 GPtrArray *input_cells = g_ptr_array_new ();
632
633 if (vr) {
634 GnmEvalPos ep;
635 eval_pos_init_sheet (&ep, sp->sheet);
636 workbook_foreach_cell_in_range (&ep, vr, CELL_ITER_ALL,
637 cb_grab_cells,
638 input_cells);
639 }
640
641 return input_cells;
642 }
643
644 void
645 gnm_solver_param_set_target (GnmSolverParameters *sp, GnmCellRef const *cr)
646 {
647 if (cr) {
648 GnmExprTop const *texpr;
649 GnmCellRef cr2 = *cr;
650 /* Make reference absolute to avoid tracking problems on row/col
651 insert. */
652 cr2.row_relative = FALSE;
653 cr2.col_relative = FALSE;
654
655 texpr = gnm_expr_top_new (gnm_expr_new_cellref (&cr2));
656 dependent_managed_set_expr (&sp->target, texpr);
657 gnm_expr_top_unref (texpr);
658 } else
659 dependent_managed_set_expr (&sp->target, NULL);
660 }
661
662 const GnmCellRef *
663 gnm_solver_param_get_target (GnmSolverParameters const *sp)
664 {
665 return sp->target.texpr
666 ? gnm_expr_top_get_cellref (sp->target.texpr)
667 : NULL;
668 }
669
670 GnmCell *
671 gnm_solver_param_get_target_cell (GnmSolverParameters const *sp)
672 {
673 const GnmCellRef *cr = gnm_solver_param_get_target (sp);
674 if (!cr)
675 return NULL;
676
677 return sheet_cell_get (eval_sheet (cr->sheet, sp->sheet),
678 cr->col, cr->row);
679 }
680
681 void
682 gnm_solver_param_set_algorithm (GnmSolverParameters *sp,
683 GnmSolverFactory *algo)
684 {
685 sp->options.algorithm = algo;
686 }
687
688 gboolean
689 gnm_solver_param_valid (GnmSolverParameters const *sp, GError **err)
690 {
691 GSList *l;
692 int i;
693 GnmCell *target_cell;
694 GPtrArray *input_cells;
695 unsigned ui;
696
697 target_cell = gnm_solver_param_get_target_cell (sp);
698 if (!target_cell) {
699 g_set_error (err,
700 go_error_invalid (),
701 0,
702 _("Invalid solver target"));
703 return FALSE;
704 }
705 gnm_cell_eval (target_cell);
706
707 if (!gnm_cell_has_expr (target_cell) ||
708 target_cell->value == NULL ||
709 !VALUE_IS_FLOAT (target_cell->value)) {
710 char *tcname = gnm_solver_cell_name (target_cell, sp->sheet);
711 g_set_error (err,
712 go_error_invalid (),
713 0,
714 _("Target cell, %s, must contain a formula that evaluates to a number"),
715 tcname);
716 g_free (tcname);
717 return FALSE;
718 }
719
720 if (!gnm_solver_param_get_input (sp)) {
721 g_set_error (err,
722 go_error_invalid (),
723 0,
724 _("Invalid solver input range"));
725 return FALSE;
726 }
727 input_cells = gnm_solver_param_get_input_cells (sp);
728 for (ui = 0; ui < input_cells->len; ui++) {
729 GnmCell *cell = g_ptr_array_index (input_cells, ui);
730 if (gnm_cell_has_expr (cell)) {
731 char *cname = gnm_solver_cell_name (cell, sp->sheet);
732 g_set_error (err,
733 go_error_invalid (),
734 0,
735 _("Input cell %s contains a formula"),
736 cname);
737 g_free (cname);
738 g_ptr_array_free (input_cells, TRUE);
739 return FALSE;
740 }
741 }
742 g_ptr_array_free (input_cells, TRUE);
743
744 for (i = 1, l = sp->constraints; l; i++, l = l->next) {
745 GnmSolverConstraint *c = l->data;
746 if (!gnm_solver_constraint_valid (c, sp)) {
747 g_set_error (err,
748 go_error_invalid (),
749 0,
750 _("Solver constraint #%d is invalid"),
751 i);
752 return FALSE;
753 }
754 }
755
756 return TRUE;
757 }
758
759 static GObject *
760 gnm_solver_param_constructor (GType type,
761 guint n_construct_properties,
762 GObjectConstructParam *construct_params)
763 {
764 GObject *obj;
765 GnmSolverParameters *sp;
766
767 obj = gnm_solver_param_parent_class->constructor
768 (type, n_construct_properties, construct_params);
769 sp = GNM_SOLVER_PARAMETERS (obj);
770
771 dependent_managed_init (&sp->target, sp->sheet);
772 dependent_managed_init (&sp->input, sp->sheet);
773
774 sp->options.model_type = GNM_SOLVER_LP;
775 sp->options.max_iter = 1000;
776 sp->options.max_time_sec = 60;
777 sp->options.assume_non_negative = TRUE;
778 sp->options.scenario_name = g_strdup ("Optimal");
779 sp->options.gradient_order = 10;
780
781 return obj;
782 }
783
784 static void
785 gnm_solver_param_finalize (GObject *obj)
786 {
787 GnmSolverParameters *sp = GNM_SOLVER_PARAMETERS (obj);
788
789 dependent_managed_set_expr (&sp->target, NULL);
790 dependent_managed_set_expr (&sp->input, NULL);
791 g_slist_free_full (sp->constraints,
792 (GFreeFunc)gnm_solver_constraint_free);
793 g_free (sp->options.scenario_name);
794
795 gnm_solver_param_parent_class->finalize (obj);
796 }
797
798 static void
799 gnm_solver_param_get_property (GObject *object, guint property_id,
800 GValue *value, GParamSpec *pspec)
801 {
802 GnmSolverParameters *sp = (GnmSolverParameters *)object;
803
804 switch (property_id) {
805 case SOLP_PROP_SHEET:
806 g_value_set_object (value, sp->sheet);
807 break;
808
809 case SOLP_PROP_PROBLEM_TYPE:
810 g_value_set_enum (value, sp->problem_type);
811 break;
812
813 default:
814 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
815 break;
816 }
817 }
818
819 static void
820 gnm_solver_param_set_property (GObject *object, guint property_id,
821 GValue const *value, GParamSpec *pspec)
822 {
823 GnmSolverParameters *sp = (GnmSolverParameters *)object;
824
825 switch (property_id) {
826 case SOLP_PROP_SHEET:
827 /* We hold no ref. */
828 sp->sheet = g_value_get_object (value);
829 break;
830
831 case SOLP_PROP_PROBLEM_TYPE:
832 sp->problem_type = g_value_get_enum (value);
833 break;
834
835 default:
836 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
837 break;
838 }
839 }
840
841 static void
842 gnm_solver_param_class_init (GObjectClass *object_class)
843 {
844 gnm_solver_param_parent_class = g_type_class_peek_parent (object_class);
845
846 object_class->constructor = gnm_solver_param_constructor;
847 object_class->finalize = gnm_solver_param_finalize;
848 object_class->set_property = gnm_solver_param_set_property;
849 object_class->get_property = gnm_solver_param_get_property;
850
851 g_object_class_install_property (object_class, SOLP_PROP_SHEET,
852 g_param_spec_object ("sheet",
853 P_("Sheet"),
854 P_("Sheet"),
855 GNM_SHEET_TYPE,
856 GSF_PARAM_STATIC |
857 G_PARAM_CONSTRUCT_ONLY |
858 G_PARAM_READWRITE));
859
860 g_object_class_install_property (object_class, SOLP_PROP_PROBLEM_TYPE,
861 g_param_spec_enum ("problem-type",
862 P_("Problem Type"),
863 P_("Problem Type"),
864 GNM_SOLVER_PROBLEM_TYPE_TYPE,
865 GNM_SOLVER_MAXIMIZE,
866 GSF_PARAM_STATIC |
867 G_PARAM_READWRITE));
868 }
869
870 GSF_CLASS (GnmSolverParameters, gnm_solver_param,
871 gnm_solver_param_class_init, NULL, G_TYPE_OBJECT)
872
873 /* ------------------------------------------------------------------------- */
874
875 enum {
876 SOL_SIG_PREPARE,
877 SOL_SIG_START,
878 SOL_SIG_STOP,
879 SOL_SIG_LAST
880 };
881
882 static guint solver_signals[SOL_SIG_LAST] = { 0 };
883
884 enum {
885 SOL_PROP_0,
886 SOL_PROP_STATUS,
887 SOL_PROP_REASON,
888 SOL_PROP_PARAMS,
889 SOL_PROP_RESULT,
890 SOL_PROP_SENSITIVITY,
891 SOL_PROP_STARTTIME,
892 SOL_PROP_ENDTIME,
893 SOL_PROP_FLIP_SIGN
894 };
895
896 static GObjectClass *gnm_solver_parent_class;
897
898 static void gnm_solver_update_derived (GnmSolver *sol);
899
900 static void
901 gnm_solver_dispose (GObject *obj)
902 {
903 GnmSolver *sol = GNM_SOLVER (obj);
904
905 if (sol->status == GNM_SOLVER_STATUS_RUNNING) {
906 gboolean ok = gnm_solver_stop (sol, NULL);
907 if (ok) {
908 g_warning ("Failed to stop solver -- now what?");
909 }
910 }
911
912 g_free (sol->reason);
913 sol->reason = NULL;
914
915 if (sol->result) {
916 g_object_unref (sol->result);
917 sol->result = NULL;
918 }
919
920 if (sol->sensitivity) {
921 g_object_unref (sol->sensitivity);
922 sol->sensitivity = NULL;
923 }
924
925 if (sol->params) {
926 g_object_unref (sol->params);
927 sol->params = NULL;
928 gnm_solver_update_derived (sol);
929 }
930
931 if (sol->gradient) {
932 sol->gradient_status = 0;
933 g_ptr_array_unref (sol->gradient);
934 sol->gradient = NULL;
935 }
936
937 if (sol->hessian) {
938 sol->hessian_status = 0;
939 g_ptr_array_unref (sol->hessian);
940 sol->hessian = NULL;
941 }
942
943 gnm_solver_parent_class->dispose (obj);
944 }
945
946 static void
947 gnm_solver_get_property (GObject *object, guint property_id,
948 GValue *value, GParamSpec *pspec)
949 {
950 GnmSolver *sol = (GnmSolver *)object;
951
952 switch (property_id) {
953 case SOL_PROP_STATUS:
954 g_value_set_enum (value, sol->status);
955 break;
956
957 case SOL_PROP_REASON:
958 g_value_set_string (value, sol->reason);
959 break;
960
961 case SOL_PROP_PARAMS:
962 g_value_set_object (value, sol->params);
963 break;
964
965 case SOL_PROP_RESULT:
966 g_value_set_object (value, sol->result);
967 break;
968
969 case SOL_PROP_SENSITIVITY:
970 g_value_set_object (value, sol->sensitivity);
971 break;
972
973 case SOL_PROP_STARTTIME:
974 g_value_set_double (value, sol->starttime);
975 break;
976
977 case SOL_PROP_ENDTIME:
978 g_value_set_double (value, sol->endtime);
979 break;
980
981 case SOL_PROP_FLIP_SIGN:
982 g_value_set_boolean (value, sol->flip_sign);
983 break;
984
985 default:
986 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
987 break;
988 }
989 }
990
991 static void
992 gnm_solver_set_property (GObject *object, guint property_id,
993 GValue const *value, GParamSpec *pspec)
994 {
995 GnmSolver *sol = (GnmSolver *)object;
996
997 switch (property_id) {
998 case SOL_PROP_STATUS:
999 gnm_solver_set_status (sol, g_value_get_enum (value));
1000 break;
1001
1002 case SOL_PROP_REASON:
1003 gnm_solver_set_reason (sol, g_value_get_string (value));
1004 break;
1005
1006 case SOL_PROP_PARAMS: {
1007 GnmSolverParameters *p = g_value_dup_object (value);
1008 if (sol->params) g_object_unref (sol->params);
1009 sol->params = p;
1010 gnm_solver_update_derived (sol);
1011 break;
1012 }
1013
1014 case SOL_PROP_RESULT: {
1015 GnmSolverResult *r = g_value_dup_object (value);
1016 if (sol->result) g_object_unref (sol->result);
1017 sol->result = r;
1018 break;
1019 }
1020
1021 case SOL_PROP_SENSITIVITY: {
1022 GnmSolverSensitivity *s = g_value_dup_object (value);
1023 if (sol->sensitivity) g_object_unref (sol->sensitivity);
1024 sol->sensitivity = s;
1025 break;
1026 }
1027
1028 case SOL_PROP_STARTTIME:
1029 sol->starttime = g_value_get_double (value);
1030 break;
1031
1032 case SOL_PROP_ENDTIME:
1033 sol->endtime = g_value_get_double (value);
1034 break;
1035
1036 case SOL_PROP_FLIP_SIGN:
1037 sol->flip_sign = g_value_get_boolean (value);
1038 break;
1039
1040 default:
1041 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
1042 break;
1043 }
1044 }
1045
1046 /**
1047 * gnm_solver_prepare: (virtual prepare)
1048 * @sol: solver
1049 * @wbc: control for user interaction
1050 * @err: location to store error
1051 *
1052 * Prepare for solving. Preparation need not do anything, but may include
1053 * such tasks as checking that the model is valid for the solver and
1054 * locating necessary external programs.
1055 *
1056 * Returns: %TRUE ok success, %FALSE on error.
1057 */
1058 gboolean
1059 gnm_solver_prepare (GnmSolver *sol, WorkbookControl *wbc, GError **err)
1060 {
1061 gboolean res;
1062
1063 g_return_val_if_fail (GNM_IS_SOLVER (sol), FALSE);
1064 g_return_val_if_fail (sol->status == GNM_SOLVER_STATUS_READY, FALSE);
1065
1066 g_signal_emit (sol, solver_signals[SOL_SIG_PREPARE], 0, wbc, err, &res);
1067 return res;
1068 }
1069
1070 /**
1071 * gnm_solver_start: (virtual start)
1072 * @sol: solver
1073 * @wbc: control for user interaction
1074 * @err: location to store error
1075 *
1076 * Start the solving process. If needed, the solver will be prepared first.
1077 *
1078 * Returns: %TRUE ok success, %FALSE on error.
1079 */
1080 gboolean
1081 gnm_solver_start (GnmSolver *sol, WorkbookControl *wbc, GError **err)
1082 {
1083 gboolean res;
1084
1085 g_return_val_if_fail (sol->status == GNM_SOLVER_STATUS_READY ||
1086 sol->status == GNM_SOLVER_STATUS_PREPARED,
1087 FALSE);
1088
1089 if (sol->status == GNM_SOLVER_STATUS_READY) {
1090 res = gnm_solver_prepare (sol, wbc, err);
1091 if (!res)
1092 return FALSE;
1093 }
1094
1095 g_return_val_if_fail (sol->status == GNM_SOLVER_STATUS_PREPARED, FALSE);
1096
1097 g_signal_emit (sol, solver_signals[SOL_SIG_START], 0, wbc, err, &res);
1098 return res;
1099 }
1100
1101 /**
1102 * gnm_solver_stop: (virtual stop)
1103 * @sol: solver
1104 * @err: location to store error
1105 *
1106 * Terminate the currently-running solver.
1107 *
1108 * Returns: %TRUE ok success, %FALSE on error.
1109 */
1110 gboolean
1111 gnm_solver_stop (GnmSolver *sol, GError **err)
1112 {
1113 gboolean res;
1114
1115 g_return_val_if_fail (GNM_IS_SOLVER (sol), FALSE);
1116
1117 g_signal_emit (sol, solver_signals[SOL_SIG_STOP], 0, err, &res);
1118 return res;
1119 }
1120
1121 static double
1122 current_time (void)
1123 {
1124 GTimeVal now;
1125 g_get_current_time (&now);
1126 return now.tv_sec + (now.tv_usec / 1e6);
1127 }
1128
1129
1130 double
1131 gnm_solver_elapsed (GnmSolver *solver)
1132 {
1133 double endtime;
1134
1135 g_return_val_if_fail (GNM_IS_SOLVER (solver), 0);
1136
1137 if (solver->starttime < 0)
1138 return 0;
1139
1140 endtime = (solver->endtime < 0)
1141 ? current_time ()
1142 : solver->endtime;
1143
1144 return endtime - solver->starttime;
1145 }
1146
1147 gboolean
1148 gnm_solver_check_timeout (GnmSolver *solver)
1149 {
1150 GnmSolverParameters *sp;
1151
1152 g_return_val_if_fail (GNM_IS_SOLVER (solver), FALSE);
1153
1154 sp = solver->params;
1155
1156 if (solver->status != GNM_SOLVER_STATUS_RUNNING)
1157 return FALSE;
1158
1159 if (gnm_solver_elapsed (solver) <= sp->options.max_time_sec)
1160 return FALSE;
1161
1162 gnm_solver_stop (solver, NULL);
1163 gnm_solver_set_reason (solver, _("Timeout"));
1164
1165 return TRUE;
1166 }
1167
1168 void
1169 gnm_solver_store_result (GnmSolver *sol)
1170 {
1171 gnm_float const *solution;
1172 unsigned ui, n = sol->input_cells->len;
1173
1174 g_return_if_fail (GNM_IS_SOLVER (sol));
1175 g_return_if_fail (sol->result != NULL);
1176 g_return_if_fail (sol->result->solution);
1177
1178 solution = gnm_solver_has_solution (sol)
1179 ? sol->result->solution
1180 : NULL;
1181
1182 for (ui = 0; ui < n; ui++) {
1183 GnmCell *cell = g_ptr_array_index (sol->input_cells, ui);
1184 GnmValue *v = solution ? value_new_float (solution[ui]) : value_new_error_NA (NULL);
1185 gnm_cell_set_value (cell, v);
1186 cell_queue_recalc (cell);
1187 }
1188 }
1189
1190 gboolean
1191 gnm_solver_finished (GnmSolver *sol)
1192 {
1193 g_return_val_if_fail (GNM_IS_SOLVER (sol), TRUE);
1194
1195 switch (sol->status) {
1196
1197 case GNM_SOLVER_STATUS_READY:
1198 case GNM_SOLVER_STATUS_PREPARING:
1199 case GNM_SOLVER_STATUS_PREPARED:
1200 case GNM_SOLVER_STATUS_RUNNING:
1201 return FALSE;
1202 case GNM_SOLVER_STATUS_DONE:
1203 default:
1204 case GNM_SOLVER_STATUS_ERROR:
1205 case GNM_SOLVER_STATUS_CANCELLED:
1206 return TRUE;
1207 }
1208 }
1209
1210 void
1211 gnm_solver_set_status (GnmSolver *solver, GnmSolverStatus status)
1212 {
1213 GnmSolverStatus old_status;
1214
1215 g_return_if_fail (GNM_IS_SOLVER (solver));
1216
1217 if (status == solver->status)
1218 return;
1219
1220 gnm_solver_set_reason (solver, NULL);
1221
1222 old_status = solver->status;
1223 solver->status = status;
1224 g_object_notify (G_OBJECT (solver), "status");
1225
1226 if (status == GNM_SOLVER_STATUS_RUNNING)
1227 g_object_set (G_OBJECT (solver),
1228 "starttime", current_time (),
1229 "endtime", (double)-1,
1230 NULL);
1231 else if (old_status == GNM_SOLVER_STATUS_RUNNING)
1232 g_object_set (G_OBJECT (solver),
1233 "endtime", current_time (),
1234 NULL);
1235 }
1236
1237 void
1238 gnm_solver_set_reason (GnmSolver *solver, const char *reason)
1239 {
1240 g_return_if_fail (GNM_IS_SOLVER (solver));
1241
1242 if (g_strcmp0 (reason, solver->reason) == 0)
1243 return;
1244
1245 g_free (solver->reason);
1246 solver->reason = g_strdup (reason);
1247
1248 g_object_notify (G_OBJECT (solver), "reason");
1249 }
1250
1251
1252 gboolean
1253 gnm_solver_has_solution (GnmSolver *solver)
1254 {
1255 if (solver->result == NULL)
1256 return FALSE;
1257
1258 switch (solver->result->quality) {
1259 case GNM_SOLVER_RESULT_NONE:
1260 case GNM_SOLVER_RESULT_INFEASIBLE:
1261 case GNM_SOLVER_RESULT_UNBOUNDED:
1262 default:
1263 return FALSE;
1264 case GNM_SOLVER_RESULT_FEASIBLE:
1265 case GNM_SOLVER_RESULT_OPTIMAL:
1266 return TRUE;
1267 }
1268 }
1269
1270 gboolean
1271 gnm_solver_check_constraints (GnmSolver *solver)
1272 {
1273 GSList *l;
1274 GnmSolverParameters *sp = solver->params;
1275 GnmCell *target_cell;
1276
1277 if (sp->options.assume_non_negative ||
1278 sp->options.assume_discrete) {
1279 unsigned ui;
1280 gboolean bad;
1281
1282 for (ui = 0; ui < solver->input_cells->len; ui++) {
1283 GnmCell *cell = g_ptr_array_index (solver->input_cells, ui);
1284 gnm_float val;
1285
1286 gnm_cell_eval (cell);
1287 val = value_get_as_float (cell->value);
1288 if (sp->options.assume_non_negative && val < 0)
1289 break;
1290 if (sp->options.assume_discrete &&
1291 val != gnm_floor (val))
1292 break;
1293 }
1294 bad = (ui < solver->input_cells->len);
1295
1296 if (bad)
1297 return FALSE;
1298 }
1299
1300 for (l = sp->constraints; l; l = l->next) {
1301 GnmSolverConstraint *c = l->data;
1302 int i;
1303 gnm_float cl, cr;
1304 GnmCell *lhs, *rhs;
1305
1306 for (i = 0;
1307 gnm_solver_constraint_get_part (c, sp, i,
1308 &lhs, &cl,
1309 &rhs, &cr);
1310 i++) {
1311 if (lhs) {
1312 gnm_cell_eval (lhs);
1313 cl = value_get_as_float (lhs->value);
1314 }
1315 if (rhs) {
1316 gnm_cell_eval (rhs);
1317 cr = value_get_as_float (rhs->value);
1318 }
1319
1320 switch (c->type) {
1321 case GNM_SOLVER_INTEGER:
1322 if (cl == gnm_floor (cl))
1323 continue;
1324 return FALSE;
1325 case GNM_SOLVER_BOOLEAN:
1326 if (cl == 0 || cl == 1)
1327 continue;
1328 return FALSE;
1329 case GNM_SOLVER_LE:
1330 if (cl <= cr)
1331 continue;
1332 return FALSE;
1333 case GNM_SOLVER_GE:
1334 if (cl >= cr)
1335 continue;
1336 return FALSE;
1337 case GNM_SOLVER_EQ:
1338 if (cl == cr)
1339 continue;
1340 return FALSE;
1341 default:
1342 g_assert_not_reached ();
1343 return FALSE;
1344 }
1345 }
1346 }
1347
1348 target_cell = gnm_solver_param_get_target_cell (sp);
1349 gnm_cell_eval (target_cell);
1350 if (!target_cell || !VALUE_IS_NUMBER (target_cell->value))
1351 return FALSE;
1352
1353 return TRUE;
1354 }
1355
1356 gboolean
1357 gnm_solver_saveas (GnmSolver *solver, WorkbookControl *wbc,
1358 GOFileSaver *fs,
1359 const char *templ, char **filename,
1360 GError **err)
1361 {
1362 int fd;
1363 GsfOutput *output;
1364 FILE *file;
1365 GOIOContext *io_context;
1366 gboolean ok;
1367 WorkbookView *wbv = wb_control_view (wbc);
1368
1369 fd = g_file_open_tmp (templ, filename, err);
1370 if (fd == -1) {
1371 g_set_error (err, G_FILE_ERROR, 0,
1372 _("Failed to create file for linear program"));
1373 return FALSE;
1374 }
1375
1376 file = fdopen (fd, "wb");
1377 if (!file) {
1378 /* This shouldn't really happen. */
1379 close (fd);
1380 g_set_error (err, G_FILE_ERROR, 0,
1381 _("Failed to create linear program file"));
1382 return FALSE;
1383 }
1384
1385 /* Give the saver a way to talk to the solver. */
1386 g_object_set_data_full (G_OBJECT (fs),
1387 "solver", g_object_ref (solver),
1388 (GDestroyNotify)g_object_unref);
1389
1390 output = gsf_output_stdio_new_FILE (*filename, file, TRUE);
1391 io_context = go_io_context_new (GO_CMD_CONTEXT (wbc));
1392 workbook_view_save_to_output (wbv, fs, output, io_context);
1393 ok = !go_io_error_occurred (io_context);
1394 g_object_unref (io_context);
1395 g_object_unref (output);
1396
1397 g_object_set_data (G_OBJECT (fs), "solver", NULL);
1398
1399 if (!ok) {
1400 g_set_error (err, G_FILE_ERROR, 0,
1401 _("Failed to save linear program"));
1402 return FALSE;
1403 }
1404
1405 return TRUE;
1406 }
1407
1408 int
1409 gnm_solver_cell_index (GnmSolver *solver, GnmCell const *cell)
1410 {
1411 gpointer idx;
1412
1413 if (g_hash_table_lookup_extended (solver->index_from_cell,
1414 (gpointer)cell,
1415 NULL, &idx))
1416 return GPOINTER_TO_INT (idx);
1417 else
1418 return -1;
1419 }
1420
1421 static int
1422 cell_in_cr (GnmSolver *sol, GnmCell const *cell, gboolean follow)
1423 {
1424 int idx;
1425
1426 if (!cell)
1427 return -1;
1428
1429 idx = gnm_solver_cell_index (sol, cell);
1430 if (idx < 0 && follow) {
1431 /* If the expression is just =X42 then look at X42 instead.
1432 This is because the mps loader uses such a level of
1433 indirection. Note: we follow only one such step. */
1434 GnmCellRef const *cr = gnm_expr_top_get_cellref (cell->base.texpr);
1435 GnmCellRef cr2;
1436 GnmCell const *new_cell;
1437 GnmEvalPos ep;
1438
1439 if (!cr)
1440 return -1;
1441
1442 eval_pos_init_cell (&ep, cell);
1443 gnm_cellref_make_abs (&cr2, cr, &ep);
1444 new_cell = sheet_cell_get (eval_sheet (cr2.sheet, cell->base.sheet),
1445 cr2.col, cr2.row);
1446 return cell_in_cr (sol, new_cell, FALSE);
1447 }
1448
1449 return idx;
1450 }
1451
1452 static gboolean
1453 cell_is_constant (GnmCell *cell, gnm_float *pc)
1454 {
1455 if (!cell)
1456 return TRUE;
1457
1458 if (cell->base.texpr)
1459 return FALSE;
1460
1461 gnm_cell_eval (cell);
1462 *pc = value_get_as_float (cell->value);
1463 return TRUE;
1464 }
1465
1466 #define SET_LOWER(l_) \
1467 do { \
1468 sol->min[idx] = MAX (sol->min[idx], (l_)); \
1469 } while (0)
1470
1471 #define SET_UPPER(l_) \
1472 do { \
1473 sol->max[idx] = MIN (sol->max[idx], (l_)); \
1474 } while (0)
1475
1476
1477
1478 static void
1479 gnm_solver_update_derived (GnmSolver *sol)
1480 {
1481 GnmSolverParameters *params = sol->params;
1482
1483 if (sol->input_cells) {
1484 g_ptr_array_free (sol->input_cells, TRUE);
1485 sol->input_cells = NULL;
1486 }
1487
1488 if (sol->index_from_cell) {
1489 g_hash_table_destroy (sol->index_from_cell);
1490 sol->index_from_cell = NULL;
1491 }
1492 sol->target = NULL;
1493
1494 g_free (sol->min);
1495 sol->min = NULL;
1496
1497 g_free (sol->max);
1498 sol->max = NULL;
1499
1500 g_free (sol->discrete);
1501 sol->discrete = NULL;
1502
1503 if (params) {
1504 unsigned ui, n;
1505 GSList *l;
1506
1507 sol->target = gnm_solver_param_get_target_cell (params);
1508
1509 sol->input_cells = gnm_solver_param_get_input_cells (params);
1510
1511 n = sol->input_cells->len;
1512 sol->index_from_cell = g_hash_table_new (g_direct_hash, g_direct_equal);
1513 for (ui = 0; ui < n; ui++) {
1514 GnmCell *cell = g_ptr_array_index (sol->input_cells, ui);
1515 g_hash_table_insert (sol->index_from_cell, cell, GUINT_TO_POINTER (ui));
1516 }
1517
1518 sol->min = g_new (gnm_float, n);
1519 sol->max = g_new (gnm_float, n);
1520 sol->discrete = g_new (guint8, n);
1521 for (ui = 0; ui < n; ui++) {
1522 sol->min[ui] = params->options.assume_non_negative ? 0 : gnm_ninf;
1523 sol->max[ui] = gnm_pinf;
1524 sol->discrete[ui] = params->options.assume_discrete;
1525 }
1526
1527 for (l = params->constraints; l; l = l->next) {
1528 GnmSolverConstraint *c = l->data;
1529 int i;
1530 gnm_float cl, cr;
1531 GnmCell *lhs, *rhs;
1532
1533 for (i = 0;
1534 gnm_solver_constraint_get_part (c, params, i,
1535 &lhs, &cl,
1536 &rhs, &cr);
1537 i++) {
1538 int idx = cell_in_cr (sol, lhs, TRUE);
1539
1540 if (idx < 0)
1541 continue;
1542 if (!cell_is_constant (rhs, &cr))
1543 continue;
1544
1545 switch (c->type) {
1546 case GNM_SOLVER_INTEGER:
1547 sol->discrete[idx] = TRUE;
1548 break;
1549 case GNM_SOLVER_BOOLEAN:
1550 sol->discrete[idx] = TRUE;
1551 SET_LOWER (0.0);
1552 SET_UPPER (1.0);
1553 break;
1554 case GNM_SOLVER_LE:
1555 SET_UPPER (cr);
1556 break;
1557 case GNM_SOLVER_GE:
1558 SET_LOWER (cr);
1559 break;
1560 case GNM_SOLVER_EQ:
1561 SET_LOWER (cr);
1562 SET_UPPER (cr);
1563 break;
1564
1565 default:
1566 g_assert_not_reached ();
1567 break;
1568 }
1569 }
1570 }
1571
1572 /*
1573 * If parameters are discrete, narrow the range by eliminating
1574 * the fractional part of the limits.
1575 */
1576 for (ui = 0; ui < n; ui++) {
1577 if (sol->discrete[ui]) {
1578 sol->min[ui] = gnm_ceil (sol->min[ui]);
1579 sol->max[ui] = gnm_floor (sol->max[ui]);
1580 }
1581 }
1582 }
1583 }
1584
1585 #undef SET_LOWER
1586 #undef SET_UPPER
1587
1588
1589 #define ADD_HEADER(txt_) do { \
1590 dao_set_bold (dao, 0, R, 0, R); \
1591 dao_set_cell (dao, 0, R, (txt_)); \
1592 R++; \
1593 } while (0)
1594
1595 #define AT_LIMIT(s_,l_) \
1596 (gnm_finite (l_) ? gnm_abs ((s_) - (l_)) <= (gnm_abs ((s_)) + gnm_abs ((l_))) / 1e10 : (s_) == (l_))
1597
1598 #define MARK_BAD(col_) \
1599 do { \
1600 int c = (col_); \
1601 dao_set_colors (dao, c, R, c, R, \
1602 gnm_color_new_rgb8 (255, 0, 0), \
1603 NULL); \
1604 } while (0)
1605
1606
1607 static void
1608 add_value_or_special (data_analysis_output_t *dao, int col, int row,
1609 gnm_float x)
1610 {
1611 if (gnm_finite (x))
1612 dao_set_cell_float (dao, col, row, x);
1613 else {
1614 dao_set_cell (dao, col, row, "-");
1615 dao_set_align (dao, col, row, col, row,
1616 GNM_HALIGN_CENTER, GNM_VALIGN_TOP);
1617 }
1618 }
1619
1620 static void
1621 print_vector (const char *name, const gnm_float *v, int n)
1622 {
1623 int i;
1624
1625 if (name)
1626 g_printerr ("%s:\n", name);
1627 for (i = 0; i < n; i++)
1628 g_printerr ("%15.8" GNM_FORMAT_f " ", v[i]);
1629 g_printerr ("\n");
1630 }
1631
1632 static void
1633 gnm_solver_create_program_report (GnmSolver *solver, const char *name)
1634 {
1635 GnmSolverParameters *params = solver->params;
1636 int R = 0;
1637 data_analysis_output_t *dao;
1638 GSList *l;
1639
1640 dao = dao_init_new_sheet (NULL);
1641 dao->sheet = params->sheet;
1642 dao_prepare_output (NULL, dao, name);
1643
1644 /* ---------------------------------------- */
1645
1646 {
1647 char *tmp;
1648
1649 ADD_HEADER (_("Target"));
1650 dao_set_cell (dao, 1, R, _("Cell"));
1651 dao_set_cell (dao, 2, R, _("Value"));
1652 dao_set_cell (dao, 3, R, _("Type"));
1653 dao_set_cell (dao, 4, R, _("Status"));
1654 R++;
1655
1656 tmp = gnm_solver_cell_name
1657 (gnm_solver_param_get_target_cell (params),
1658 params->sheet);
1659 dao_set_cell (dao, 1, R, tmp);
1660 g_free (tmp);
1661
1662 dao_set_cell_float (dao, 2, R, solver->result->value);
1663
1664 switch (params->problem_type) {
1665 case GNM_SOLVER_MINIMIZE:
1666 dao_set_cell (dao, 3, R, _("Minimize"));
1667 break;
1668 case GNM_SOLVER_MAXIMIZE:
1669 dao_set_cell (dao, 3, R, _("Maximize"));
1670 break;
1671 }
1672
1673 switch (solver->result->quality) {
1674 default:
1675 break;
1676 case GNM_SOLVER_RESULT_FEASIBLE:
1677 dao_set_cell (dao, 4, R, _("Feasible"));
1678 break;
1679 case GNM_SOLVER_RESULT_OPTIMAL:
1680 dao_set_cell (dao, 4, R, _("Optimal"));
1681 break;
1682 }
1683
1684 R++;
1685 R++;
1686 }
1687
1688 /* ---------------------------------------- */
1689
1690 if (solver->input_cells->len > 0) {
1691 unsigned ui;
1692
1693 ADD_HEADER (_("Variables"));
1694
1695 dao_set_cell (dao, 1, R, _("Cell"));
1696 dao_set_cell (dao, 2, R, _("Value"));
1697 dao_set_cell (dao, 3, R, _("Lower"));
1698 dao_set_cell (dao, 4, R, _("Upper"));
1699 dao_set_cell (dao, 5, R, _("Slack"));
1700 R++;
1701
1702 for (ui = 0; ui < solver->input_cells->len; ui++) {
1703 GnmCell *cell = g_ptr_array_index (solver->input_cells, ui);
1704 gnm_float L = solver->min[ui];
1705 gnm_float H = solver->max[ui];
1706 gnm_float s = solver->result->solution[ui];
1707 gnm_float slack = MIN (s - L, H - s);
1708
1709 char *cname = gnm_solver_cell_name (cell, params->sheet);
1710 dao_set_cell (dao, 1, R, cname);
1711 g_free (cname);
1712 dao_set_cell_value (dao, 2, R, value_new_float (s));
1713 add_value_or_special (dao, 3, R, L);
1714 add_value_or_special (dao, 4, R, H);
1715
1716 add_value_or_special (dao, 5, R, slack);
1717 if (slack < 0)
1718 MARK_BAD (5);
1719
1720 if (AT_LIMIT (s, L) || AT_LIMIT (s, H))
1721 dao_set_cell (dao, 6, R, _("At limit"));
1722
1723 if (s < L || s > H) {
1724 dao_set_cell (dao, 7, R, _("Outside bounds"));
1725 MARK_BAD (7);
1726 }
1727
1728 R++;
1729 }
1730
1731 R++;
1732 }
1733
1734 /* ---------------------------------------- */
1735
1736 ADD_HEADER (_("Constraints"));
1737
1738 if (params->constraints) {
1739 dao_set_cell (dao, 1, R, _("Condition"));
1740 dao_set_cell (dao, 2, R, _("Value"));
1741 dao_set_cell (dao, 3, R, _("Limit"));
1742 dao_set_cell (dao, 4, R, _("Slack"));
1743 } else {
1744 dao_set_cell (dao, 1, R, _("No constraints"));
1745 }
1746 R++;
1747
1748 for (l = params->constraints; l; l = l->next) {
1749 GnmSolverConstraint *c = l->data;
1750 int i;
1751 gnm_float cl, cr;
1752 GnmCell *lhs, *rhs;
1753
1754 for (i = 0;
1755 gnm_solver_constraint_get_part (c, params, i,
1756 &lhs, &cl,
1757 &rhs, &cr);
1758 i++) {
1759 gnm_float slack = 0;
1760 char *ctxt = gnm_solver_constraint_part_as_str (c, i, params);
1761 dao_set_cell (dao, 1, R, ctxt);
1762 g_free (ctxt);
1763
1764 if (lhs) {
1765 gnm_cell_eval (lhs);
1766 cl = value_get_as_float (lhs->value);
1767 }
1768 if (rhs) {
1769 gnm_cell_eval (rhs);
1770 cr = value_get_as_float (rhs->value);
1771 }
1772
1773 switch (c->type) {
1774 case GNM_SOLVER_INTEGER: {
1775 gnm_float c = gnm_fake_round (cl);
1776 slack = 0 - gnm_abs (c - cl);
1777 break;
1778 }
1779 case GNM_SOLVER_BOOLEAN: {
1780 gnm_float c = (cl > 0.5 ? 1 : 0);
1781 slack = 0 - gnm_abs (c - cl);
1782 break;
1783 }
1784 case GNM_SOLVER_LE:
1785 slack = cr - cl;
1786 break;
1787 case GNM_SOLVER_GE:
1788 slack = cl - cr;
1789 break;
1790 case GNM_SOLVER_EQ:
1791 slack = 0 - gnm_abs (cl - cr);
1792 break;
1793 default:
1794 g_assert_not_reached ();
1795 }
1796
1797 add_value_or_special (dao, 2, R, cl);
1798 if (rhs)
1799 add_value_or_special (dao, 3, R, cr);
1800
1801 add_value_or_special (dao, 4, R, slack);
1802 if (slack < 0)
1803 MARK_BAD (4);
1804
1805 R++;
1806 }
1807 }
1808
1809 /* ---------------------------------------- */
1810
1811 dao_autofit_columns (dao);
1812 dao_redraw_respan (dao);
1813
1814 dao_free (dao);
1815 }
1816
1817 static void
1818 gnm_solver_create_sensitivity_report (GnmSolver *solver, const char *name)
1819 {
1820 GnmSolverParameters *params = solver->params;
1821 GnmSolverSensitivity *sols = solver->sensitivity;
1822 int R = 0;
1823 data_analysis_output_t *dao;
1824 GSList *l;
1825
1826 if (!sols)
1827 return;
1828
1829 dao = dao_init_new_sheet (NULL);
1830 dao->sheet = params->sheet;
1831 dao_prepare_output (NULL, dao, name);
1832
1833 /* ---------------------------------------- */
1834
1835 if (solver->input_cells->len > 0) {
1836 unsigned ui;
1837
1838 ADD_HEADER (_("Variables"));
1839
1840 dao_set_cell (dao, 1, R, _("Cell"));
1841 dao_set_cell (dao, 2, R, _("Final\nValue"));
1842 dao_set_cell (dao, 3, R, _("Reduced\nCost"));
1843 dao_set_cell (dao, 4, R, _("Lower\nLimit"));
1844 dao_set_cell (dao, 5, R, _("Upper\nLimit"));
1845 dao_set_align (dao, 1, R, 5, R, GNM_HALIGN_CENTER, GNM_VALIGN_BOTTOM);
1846 dao_autofit_these_rows (dao, R, R);
1847 R++;
1848
1849 for (ui = 0; ui < solver->input_cells->len; ui++) {
1850 GnmCell *cell = g_ptr_array_index (solver->input_cells, ui);
1851 gnm_float L = sols->vars[ui].low;
1852 gnm_float H = sols->vars[ui].high;
1853 gnm_float red = sols->vars[ui].reduced_cost;
1854 gnm_float s = solver->result->solution[ui];
1855
1856 char *cname = gnm_solver_cell_name (cell, params->sheet);
1857 dao_set_cell (dao, 1, R, cname);
1858 g_free (cname);
1859 dao_set_cell_value (dao, 2, R, value_new_float (s));
1860 add_value_or_special (dao, 3, R, red);
1861 add_value_or_special (dao, 4, R, L);
1862 add_value_or_special (dao, 5, R, H);
1863
1864 R++;
1865 }
1866
1867 R++;
1868 }
1869
1870 /* ---------------------------------------- */
1871
1872 ADD_HEADER (_("Constraints"));
1873
1874 if (params->constraints) {
1875 dao_set_cell (dao, 1, R, _("Constraint"));
1876 dao_set_cell (dao, 2, R, _("Shadow\nPrice"));
1877 dao_set_cell (dao, 3, R, _("Constraint\nLHS"));
1878 dao_set_cell (dao, 4, R, _("Constraint\nRHS"));
1879 dao_set_cell (dao, 5, R, _("Lower\nLimit"));
1880 dao_set_cell (dao, 6, R, _("Upper\nLimit"));
1881 dao_set_align (dao, 1, R, 6, R, GNM_HALIGN_CENTER, GNM_VALIGN_BOTTOM);
1882 dao_autofit_these_rows (dao, R, R);
1883 } else {
1884 dao_set_cell (dao, 1, R, _("No constraints"));
1885 }
1886 R++;
1887
1888 for (l = params->constraints; l; l = l->next) {
1889 GnmSolverConstraint *c = l->data;
1890 int i, cidx = 0;
1891 gnm_float cl, cr;
1892 GnmCell *lhs, *rhs;
1893
1894 for (i = 0;
1895 gnm_solver_constraint_get_part (c, params, i,
1896 &lhs, &cl,
1897 &rhs, &cr);
1898 i++, cidx++) {
1899 char *ctxt;
1900
1901 switch (c->type) {
1902 case GNM_SOLVER_INTEGER:
1903 case GNM_SOLVER_BOOLEAN:
1904 continue;
1905 default:
1906 ; // Nothing
1907 }
1908
1909 ctxt = gnm_solver_constraint_part_as_str (c, i, params);
1910 dao_set_cell (dao, 1, R, ctxt);
1911 g_free (ctxt);
1912
1913 if (lhs) {
1914 gnm_cell_eval (lhs);
1915 cl = value_get_as_float (lhs->value);
1916 }
1917 if (rhs) {
1918 gnm_cell_eval (rhs);
1919 cr = value_get_as_float (rhs->value);
1920 }
1921
1922 add_value_or_special (dao, 2, R, sols->constraints[cidx].shadow_price);
1923 add_value_or_special (dao, 3, R, cl);
1924 add_value_or_special (dao, 4, R, cr);
1925 add_value_or_special (dao, 5, R, sols->constraints[cidx].low);
1926 add_value_or_special (dao, 6, R, sols->constraints[cidx].high);
1927
1928 R++;
1929 }
1930 }
1931
1932 /* ---------------------------------------- */
1933
1934 dao_autofit_columns (dao);
1935 dao_redraw_respan (dao);
1936
1937 dao_free (dao);
1938 }
1939
1940 void
1941 gnm_solver_create_report (GnmSolver *solver, const char *base)
1942 {
1943 GnmSolverParameters *params = solver->params;
1944
1945 if (params->options.program_report) {
1946 char *name = g_strdup_printf (base, _("Program"));
1947 gnm_solver_create_program_report (solver, name);
1948 g_free (name);
1949 }
1950
1951 if (params->options.sensitivity_report) {
1952 char *name = g_strdup_printf (base, _("Sensitivity"));
1953 gnm_solver_create_sensitivity_report (solver, name);
1954 g_free (name);
1955 }
1956 }
1957
1958 #undef AT_LIMIT
1959 #undef ADD_HEADER
1960 #undef MARK_BAD
1961
1962 /**
1963 * gnm_solver_get_target_value:
1964 * @solver: solver
1965 *
1966 * Returns: the current value of the target cell, possibly with the sign
1967 * flipped.
1968 */
1969 gnm_float
1970 gnm_solver_get_target_value (GnmSolver *solver)
1971 {
1972 GnmValue const *v;
1973
1974 gnm_cell_eval (solver->target);
1975 v = solver->target->value;
1976
1977 if (VALUE_IS_NUMBER (v) || VALUE_IS_EMPTY (v)) {
1978 gnm_float y = value_get_as_float (v);
1979 return solver->flip_sign ? 0 - y : y;
1980 } else
1981 return gnm_nan;
1982 }
1983
1984 void
1985 gnm_solver_set_var (GnmSolver *sol, int i, gnm_float x)
1986 {
1987 GnmCell *cell = g_ptr_array_index (sol->input_cells, i);
1988
1989 if (cell->value &&
1990 VALUE_IS_FLOAT (cell->value) &&
1991 value_get_as_float (cell->value) == x)
1992 return;
1993
1994 gnm_cell_set_value (cell, value_new_float (x));
1995 cell_queue_recalc (cell);
1996 }
1997
1998 void
1999 gnm_solver_set_vars (GnmSolver *sol, gnm_float const *xs)
2000 {
2001 const int n = sol->input_cells->len;
2002 int i;
2003
2004 for (i = 0; i < n; i++)
2005 gnm_solver_set_var (sol, i, xs[i]);
2006 }
2007
2008 /**
2009 * gnm_solver_save_vars:
2010 * @sol: #GnmSolver
2011 *
2012 * Returns: (transfer full) (element-type GnmValue):
2013 */
2014 GPtrArray *
2015 gnm_solver_save_vars (GnmSolver *sol)
2016 {
2017 GPtrArray *vals = g_ptr_array_new ();
2018 unsigned ui;
2019
2020 for (ui = 0; ui < sol->input_cells->len; ui++) {
2021 GnmCell *cell = g_ptr_array_index (sol->input_cells, ui);
2022 g_ptr_array_add (vals, value_dup (cell->value));
2023 }
2024
2025 return vals;
2026 }
2027
2028 /**
2029 * gnm_solver_restore_vars:
2030 * @sol: #GnmSolver
2031 * @vals: (transfer full) (element-type GnmValue): values to restore
2032 */
2033 void
2034 gnm_solver_restore_vars (GnmSolver *sol, GPtrArray *vals)
2035 {
2036 unsigned ui;
2037
2038 for (ui = 0; ui < sol->input_cells->len; ui++) {
2039 GnmCell *cell = g_ptr_array_index (sol->input_cells, ui);
2040 gnm_cell_set_value (cell, g_ptr_array_index (vals, ui));
2041 cell_queue_recalc (cell);
2042 }
2043
2044 g_ptr_array_free (vals, TRUE);
2045 }
2046
2047 /**
2048 * gnm_solver_has_analytic_gradient:
2049 * @sol: the solver
2050 *
2051 * Returns: %TRUE if the gradient can be computed analytically.
2052 */
2053 gboolean
2054 gnm_solver_has_analytic_gradient (GnmSolver *sol)
2055 {
2056 const int n = sol->input_cells->len;
2057
2058 if (sol->gradient_status == 0) {
2059 int i;
2060
2061 sol->gradient_status++;
2062
2063 sol->gradient = g_ptr_array_new_with_free_func ((GDestroyNotify)gnm_expr_top_unref);
2064 for (i = 0; i < n; i++) {
2065 GnmCell *cell = g_ptr_array_index (sol->input_cells, i);
2066 GnmExprTop const *te =
2067 gnm_expr_cell_deriv (sol->target, cell);
2068 if (te)
2069 g_ptr_array_add (sol->gradient, (gpointer)te);
2070 else {
2071 if (gnm_solver_debug ())
2072 g_printerr ("Unable to compute analytic gradient\n");
2073 sol->gradient_status++;
2074 break;
2075 }
2076 }
2077 }
2078
2079 return sol->gradient_status == 1;
2080 }
2081
2082 static gnm_float *
2083 gnm_solver_compute_gradient_analytically (GnmSolver *sol, gnm_float const *xs)
2084 {
2085 const int n = sol->input_cells->len;
2086 int i;
2087 gnm_float *g = g_new (gnm_float, n);
2088 GnmEvalPos ep;
2089
2090 eval_pos_init_cell (&ep, sol->target);
2091 for (i = 0; i < n; i++) {
2092 GnmExprTop const *te = g_ptr_array_index (sol->gradient, i);
2093 GnmValue *v = gnm_expr_top_eval
2094 (te, &ep, GNM_EXPR_EVAL_SCALAR_NON_EMPTY);
2095 g[i] = VALUE_IS_NUMBER (v) ? value_get_as_float (v) : gnm_nan;
2096 if (sol->flip_sign)
2097 g[i] = 0 - g[i];
2098 value_release (v);
2099 }
2100
2101 if (gnm_solver_debug ())
2102 print_vector ("Analytic gradient", g, n);
2103
2104 return g;
2105 }
2106
2107 /**
2108 * gnm_solver_compute_gradient: (skip)
2109 * @sol: Solver
2110 * @xs: Point to compute gradient at
2111 *
2112 * Returns: xxx(transfer full): A vector containing the gradient. This
2113 * function takes the flip-sign property into account. This will use
2114 * analytic gradient, if possible, and a numerical approximation otherwise.
2115 */
2116 gnm_float *
2117 gnm_solver_compute_gradient (GnmSolver *sol, gnm_float const *xs)
2118 {
2119 gnm_float *g;
2120 gnm_float y0;
2121 const int n = sol->input_cells->len;
2122 int i;
2123 const int order = sol->params->options.gradient_order;
2124
2125 gnm_solver_set_vars (sol, xs);
2126 y0 = gnm_solver_get_target_value (sol);
2127
2128 if (gnm_solver_has_analytic_gradient (sol))
2129 return gnm_solver_compute_gradient_analytically (sol, xs);
2130
2131 g = g_new (gnm_float, n);
2132 for (i = 0; i < n; i++) {
2133 gnm_float x0 = xs[i];
2134 gnm_float dx, dy;
2135 int j;
2136
2137 /*
2138 * This computes the least-squares fit of an affine function
2139 * based on 2*order equidistant points symmetrically around
2140 * the value we compute the derivative for.
2141 *
2142 * We use an even number of ULPs for the step size. This
2143 * ensures that the x values are computed without rounding
2144 * error except, potentially, a single step that crosses an
2145 * integer power of 2.
2146 */
2147 dx = 16 * (go_add_epsilon (x0) - x0);
2148 dy = 0;
2149 for (j = -order; j <= order; j++) {
2150 gnm_float y;
2151
2152 if (j == 0)
2153 continue;
2154
2155 gnm_solver_set_var (sol, i, x0 + j * dx);
2156 y = gnm_solver_get_target_value (sol);
2157 dy += j * (y - y0);
2158 }
2159 dy /= 2 * (order * (2 * order * order + 3 * order + 1) / 6);
2160 g[i] = dy / dx;
2161
2162 gnm_solver_set_var (sol, i, x0);
2163 }
2164
2165 if (gnm_solver_debug ())
2166 print_vector ("Numerical gradient", g, n);
2167
2168 return g;
2169 }
2170
2171 /**
2172 * gnm_solver_has_analytic_hessian:
2173 * @sol: the solver
2174 *
2175 * Returns: %TRUE if the Hessian can be computed analytically.
2176 */
2177 gboolean
2178 gnm_solver_has_analytic_hessian (GnmSolver *sol)
2179 {
2180 const int n = sol->input_cells->len;
2181 int i, j;
2182 GnmEvalPos ep;
2183 GnmExprDeriv *info;
2184
2185 if (!gnm_solver_has_analytic_gradient (sol))
2186 sol->hessian_status = sol->gradient_status;
2187
2188 if (sol->hessian_status)
2189 return sol->hessian_status == 1;
2190
2191 sol->hessian_status++;
2192 sol->hessian = g_ptr_array_new_with_free_func ((GDestroyNotify)gnm_expr_top_unref);
2193
2194 eval_pos_init_cell (&ep, sol->target);
2195 info = gnm_expr_deriv_info_new ();
2196 for (i = 0; i < n && sol->hessian_status == 1; i++) {
2197 GnmExprTop const *gi = g_ptr_array_index (sol->gradient, i);
2198 for (j = i; j < n; j++) {
2199 GnmCell *cell;
2200 GnmExprTop const *te;
2201 GnmEvalPos var;
2202
2203 cell = g_ptr_array_index (sol->input_cells, j);
2204 eval_pos_init_cell (&var, cell);
2205 gnm_expr_deriv_info_set_var (info, &var);
2206 te = gnm_expr_top_deriv (gi, &ep, info);
2207
2208 if (te)
2209 g_ptr_array_add (sol->hessian, (gpointer)te);
2210 else {
2211 if (gnm_solver_debug ())
2212 g_printerr ("Unable to compute analytic hessian\n");
2213 sol->hessian_status++;
2214 break;
2215 }
2216 }
2217 }
2218
2219 gnm_expr_deriv_info_free (info);
2220
2221 return sol->hessian_status == 1;
2222 }
2223
2224 /**
2225 * gnm_solver_compute_hessian: (skip)
2226 * @sol: Solver
2227 * @xs: Point to compute Hessian at
2228 *
2229 * Returns: (transfer full): A vector containing the Hessian. This
2230 * function takes the flip-sign property into account. The result vector
2231 * will be n+(n-1)+...+2+1 elements long containing the triangular
2232 * Hessian. Use symmetry to obtain the full Hessian.
2233 */
2234 GnmMatrix *
2235 gnm_solver_compute_hessian (GnmSolver *sol, gnm_float const *xs)
2236 {
2237 int i, j, k;
2238 GnmMatrix *H;
2239 GnmEvalPos ep;
2240 int const n = sol->input_cells->len;
2241
2242 if (!gnm_solver_has_analytic_hessian (sol))
2243 return NULL;
2244
2245 gnm_solver_set_vars (sol, xs);
2246
2247 H = gnm_matrix_new (n, n);
2248 eval_pos_init_cell (&ep, sol->target);
2249 for (i = k = 0; i < n; i++) {
2250 for (j = i; j < n; j++, k++) {
2251 GnmExprTop const *te =
2252 g_ptr_array_index (sol->hessian, k);
2253 GnmValue *v = gnm_expr_top_eval
2254 (te, &ep, GNM_EXPR_EVAL_SCALAR_NON_EMPTY);
2255 gnm_float x = VALUE_IS_NUMBER (v)
2256 ? value_get_as_float (v)
2257 : gnm_nan;
2258 if (sol->flip_sign)
2259 x = 0 - x;
2260 value_release (v);
2261
2262 H->data[i][j] = x;
2263 H->data[j][i] = x;
2264 }
2265 }
2266
2267 return H;
2268 }
2269
2270
2271 static gnm_float
2272 try_step (GnmSolver *sol, gnm_float const *x0, gnm_float const *dir, gnm_float step)
2273 {
2274 int const n = sol->input_cells->len;
2275 gnm_float *x = g_new (gnm_float, n);
2276 int i;
2277 gnm_float y;
2278
2279 for (i = 0; i < n; i++)
2280 x[i] = x0[i] + step * dir[i];
2281 gnm_solver_set_vars (sol, x);
2282 y = gnm_solver_get_target_value (sol);
2283 g_free (x);
2284 return y;
2285 }
2286
2287 /**
2288 * gnm_solver_line_search:
2289 * @sol: Solver
2290 * @x0: Starting point
2291 * @dir: direction
2292 * @try_reverse: whether to try reverse direction at all
2293 * @step: initial step size
2294 * @max_step: largest allowed step
2295 * @eps: tolerance for optimal step
2296 * @py: (out): location to store resulting objective function value
2297 *
2298 * Returns: optimal step size.
2299 */
2300 gnm_float
2301 gnm_solver_line_search (GnmSolver *sol,
2302 gnm_float const *x0, gnm_float const *dir,
2303 gboolean try_reverse,
2304 gnm_float step, gnm_float max_step, gnm_float eps,
2305 gnm_float *py)
2306 {
2307 /*
2308 * 0: Initial
2309 * 1: Found improvement, but not far side of it
2310 * 2: Have (s0,s1,s2) with s1 lowest
2311 */
2312 int phase = 0;
2313 gnm_float s, s0, s1, s2;
2314 gnm_float y, y0, y1, y2;
2315 gnm_float const phi = (gnm_sqrt (5) + 1) / 2;
2316 gboolean rbig;
2317 gboolean debug = gnm_solver_debug ();
2318
2319 g_return_val_if_fail (eps >= 0, gnm_nan);
2320 g_return_val_if_fail (step > 0, gnm_nan);
2321 g_return_val_if_fail (max_step >= step, gnm_nan);
2322
2323 if (debug) {
2324 g_printerr ("LS: step=%" GNM_FORMAT_g ", max=%" GNM_FORMAT_g ", eps=%" GNM_FORMAT_g "\n", step, max_step, eps);
2325 print_vector (NULL, dir, sol->input_cells->len);
2326 }
2327
2328 gnm_solver_set_vars (sol, x0);
2329 y0 = gnm_solver_get_target_value (sol);
2330 s0 = 0;
2331
2332 s = step;
2333 while (phase == 0) {
2334 gboolean flat = TRUE;
2335
2336 y = try_step (sol, x0, dir, s);
2337 if (0 && debug)
2338 g_printerr ("LS0: s:%.6" GNM_FORMAT_g " y:%.6" GNM_FORMAT_g "\n",
2339 s, y);
2340 if (y < y0 && gnm_solver_check_constraints (sol)) {
2341 y1 = y;
2342 s1 = s;
2343 phase = 1;
2344 break;
2345 } else if (y != y0)
2346 flat = FALSE;
2347
2348 if (try_reverse) {
2349 y = try_step (sol, x0, dir, -s);
2350 if (0 && debug)
2351 g_printerr ("LS0: s:%.6" GNM_FORMAT_g " y:%.6" GNM_FORMAT_g "\n",
2352 -s, y);
2353 if (y < y0 && gnm_solver_check_constraints (sol)) {
2354 y1 = y;
2355 s1 = -s;
2356 phase = 1;
2357 break;
2358 } else if (y != y0)
2359 flat = FALSE;
2360 }
2361
2362 s /= 32;
2363
2364 if (s <= 0 || flat)
2365 return gnm_nan;
2366 }
2367
2368 while (phase == 1) {
2369 s = s1 * (phi + 1);
2370
2371 if (gnm_abs (s) >= max_step)
2372 goto bail;
2373
2374 y = try_step (sol, x0, dir, s);
2375 if (!gnm_finite (y) || !gnm_solver_check_constraints (sol))
2376 goto bail;
2377
2378 if (y < y1) {
2379 y1 = y;
2380 s1 = s;
2381 continue;
2382 }
2383
2384 y2 = y;
2385 s2 = s;
2386 phase = 2;
2387 }
2388
2389 /*
2390 * Phase 2: we have three steps, s0/s1/s2, in order (descending or ascending) such
2391 * that
2392 * 1. y1<=y0 (equality unlikely)
2393 * 2. y1<=y2 (equality unlikely)
2394 * 3a. (s2-s1)=phi*(s1-s0) or
2395 * 3b. (s2-s1)*phi=(s1-s0)
2396 */
2397 rbig = TRUE; /* Initially 3a holds. */
2398 while (phase == 2) {
2399 if (0 && debug) {
2400 gnm_float s01 = s1 - s0;
2401 gnm_float s12 = s2 - s1;
2402 g_printerr ("LS2: s0:%.6" GNM_FORMAT_g " s01=%.6" GNM_FORMAT_g " s12=%.6" GNM_FORMAT_g
2403 " r=%" GNM_FORMAT_g
2404 " y:%.10" GNM_FORMAT_g "/%.10" GNM_FORMAT_g "/%.10" GNM_FORMAT_g "\n",
2405 s0, s01, s12, s12 / s01, y0, y1, y2);
2406 }
2407
2408 s = rbig ? s1 + (s1 - s0) * (phi - 1) : s1 - (s2 - s1) * (phi - 1);
2409 if (s <= s0 || s >= s2 || gnm_abs (s - s1) <= eps)
2410 break;
2411
2412 y = try_step (sol, x0, dir, s);
2413 if (!gnm_finite (y) || !gnm_solver_check_constraints (sol))
2414 goto bail;
2415
2416 if (y < y1) {
2417 if (rbig) {
2418 y0 = y1;
2419 s0 = s1;
2420 } else {
2421 y2 = y1;
2422 s2 = s1;
2423 }
2424 y1 = y;
2425 s1 = s;
2426 } else {
2427 if (rbig) {
2428 y2 = y;
2429 s2 = s;
2430 } else {
2431 y0 = y;
2432 s0 = s;
2433 }
2434 rbig = !rbig;
2435
2436 if (y0 == y1 && y1 == y2)
2437 break;
2438 }
2439 }
2440
2441 bail:
2442 if (debug)
2443 g_printerr ("LS: step %.6" GNM_FORMAT_g "\n", s1);
2444
2445 *py = y1;
2446 return s1;
2447 }
2448
2449 /**
2450 * gnm_solver_pick_lp_coords:
2451 * @sol: Solver
2452 * @px1: (out): first coordinate value
2453 * @px2: (out): second coordinate value
2454 *
2455 * Pick two good values for each coordinate. We prefer 0 and 1
2456 * when they are valid.
2457 */
2458 void
2459 gnm_solver_pick_lp_coords (GnmSolver *sol,
2460 gnm_float **px1, gnm_float **px2)
2461 {
2462 const unsigned n = sol->input_cells->len;
2463 gnm_float *x1 = *px1 = g_new (gnm_float, n);
2464 gnm_float *x2 = *px2 = g_new (gnm_float, n);
2465 unsigned ui;
2466
2467 for (ui = 0; ui < n; ui++) {
2468 const gnm_float L = sol->min[ui], H = sol->max[ui];
2469
2470 if (L == H) {
2471 x1[ui] = x2[ui] = L;
2472 } else if (sol->discrete[ui] && H - L == 1) {
2473 x1[ui] = L;
2474 x2[ui] = H;
2475 } else {
2476 if (L <= 0 && H >= 0)
2477 x1[ui] = 0;
2478 else if (gnm_finite (L))
2479 x1[ui] = L;
2480 else
2481 x1[ui] = H;
2482
2483 if (x1[ui] + 1 <= H)
2484 x2[ui] = x1[ui] + 1;
2485 else if (x1[ui] - 1 >= H)
2486 x2[ui] = x1[ui] - 1;
2487 else if (x1[ui] != H)
2488 x2[ui] = (x1[ui] + H) / 2;
2489 else
2490 x2[ui] = (x1[ui] + L) / 2;
2491 }
2492 }
2493 }
2494
2495 /**
2496 * gnm_solver_get_lp_coeffs: (skip)
2497 * @sol: Solver
2498 * @ycell: Cell for which to compute coefficients
2499 * @x1: first coordinate value
2500 * @x2: second coordinate value
2501 * @err: error location
2502 *
2503 * Returns: xxx(transfer full) (nullable): coordinates, or %NULL in case of error.
2504 * Note: this function is not affected by the flip-sign property, even
2505 * if @ycell happens to coindice with the solver target cell.
2506 */
2507 gnm_float *
2508 gnm_solver_get_lp_coeffs (GnmSolver *sol, GnmCell *ycell,
2509 gnm_float const *x1, gnm_float const *x2,
2510 GError **err)
2511 {
2512 const unsigned n = sol->input_cells->len;
2513 unsigned ui;
2514 gnm_float *res = g_new (gnm_float, n);
2515 gnm_float y0;
2516
2517 gnm_solver_set_vars (sol, x1);
2518 gnm_cell_eval (ycell);
2519 y0 = VALUE_IS_NUMBER (ycell->value) ? value_get_as_float (ycell->value) : gnm_nan;
2520 if (!gnm_finite (y0))
2521 goto fail_calc;
2522
2523 for (ui = 0; ui < n; ui++) {
2524 gnm_float dx = x2[ui] - x1[ui], dy, y1;
2525
2526 if (dx <= 0) {
2527 res[ui] = 0;
2528 continue;
2529 }
2530
2531 gnm_solver_set_var (sol, ui, x2[ui]);
2532 gnm_cell_eval (ycell);
2533 y1 = VALUE_IS_NUMBER (ycell->value) ? value_get_as_float (ycell->value) : gnm_nan;
2534
2535 dy = y1 - y0;
2536 res[ui] = dy / dx;
2537 if (!gnm_finite (res[ui]))
2538 goto fail_calc;
2539
2540 if (!sol->discrete[ui] || dx != 1) {
2541 gnm_float x01, y01, e, emax;
2542
2543 x01 = (x1[ui] + x2[ui]) / 2;
2544 if (sol->discrete[ui]) x01 = gnm_floor (x01);
2545 gnm_solver_set_var (sol, ui, x01);
2546 gnm_cell_eval (ycell);
2547 y01 = VALUE_IS_NUMBER (ycell->value) ? value_get_as_float (ycell->value) : gnm_nan;
2548 if (!gnm_finite (y01))
2549 goto fail_calc;
2550
2551 emax = dy == 0 ? 1e-10 : gnm_abs (dy) / 1e-10;
2552 e = dy - 2 * (y01 - y0);
2553 if (gnm_abs (e) > emax)
2554 goto fail_linear;
2555 }
2556
2557 gnm_solver_set_var (sol, ui, x1[ui]);
2558 }
2559
2560 return res;
2561
2562 fail_calc:
2563 g_set_error (err,
2564 go_error_invalid (),
2565 0,
2566 _("Target cell did not evaluate to a number."));
2567 g_free (res);
2568 return NULL;
2569
2570 fail_linear:
2571 g_set_error (err,
2572 go_error_invalid (),
2573 0,
2574 _("Target cell does not appear to depend linearly on input cells."));
2575 g_free (res);
2576 return NULL;
2577 }
2578
2579
2580 static void
2581 gnm_solver_class_init (GObjectClass *object_class)
2582 {
2583 gnm_solver_parent_class = g_type_class_peek_parent (object_class);
2584
2585 object_class->dispose = gnm_solver_dispose;
2586 object_class->set_property = gnm_solver_set_property;
2587 object_class->get_property = gnm_solver_get_property;
2588
2589 g_object_class_install_property (object_class, SOL_PROP_STATUS,
2590 g_param_spec_enum ("status",
2591 P_("status"),
2592 P_("The solver's current status"),
2593 GNM_SOLVER_STATUS_TYPE,
2594 GNM_SOLVER_STATUS_READY,
2595 GSF_PARAM_STATIC |
2596 G_PARAM_READWRITE));
2597
2598 g_object_class_install_property (object_class, SOL_PROP_REASON,
2599 g_param_spec_string ("reason",
2600 P_("reason"),
2601 P_("The reason behind the solver's status"),
2602 NULL,
2603 GSF_PARAM_STATIC |
2604 G_PARAM_READWRITE));
2605
2606 g_object_class_install_property (object_class, SOL_PROP_PARAMS,
2607 g_param_spec_object ("params",
2608 P_("Parameters"),
2609 P_("Solver parameters"),
2610 GNM_SOLVER_PARAMETERS_TYPE,
2611 GSF_PARAM_STATIC |
2612 G_PARAM_CONSTRUCT_ONLY |
2613 G_PARAM_READWRITE));
2614
2615 g_object_class_install_property (object_class, SOL_PROP_RESULT,
2616 g_param_spec_object ("result",
2617 P_("Result"),
2618 P_("Current best feasible result"),
2619 GNM_SOLVER_RESULT_TYPE,
2620 GSF_PARAM_STATIC |
2621 G_PARAM_READWRITE));
2622
2623 g_object_class_install_property (object_class, SOL_PROP_SENSITIVITY,
2624 g_param_spec_object ("sensitivity",
2625 P_("Sensitivity"),
2626 P_("Sensitivity results"),
2627 GNM_SOLVER_SENSITIVITY_TYPE,
2628 GSF_PARAM_STATIC |
2629 G_PARAM_READWRITE));
2630
2631 g_object_class_install_property (object_class, SOL_PROP_STARTTIME,
2632 g_param_spec_double ("starttime",
2633 P_("Start Time"),
2634 P_("Time the solver was started"),
2635 -1, 1e10, -1,
2636 GSF_PARAM_STATIC |
2637 G_PARAM_READWRITE));
2638
2639 g_object_class_install_property (object_class, SOL_PROP_ENDTIME,
2640 g_param_spec_double ("endtime",
2641 P_("End Time"),
2642 P_("Time the solver finished"),
2643 -1, 1e10, -1,
2644 GSF_PARAM_STATIC |
2645 G_PARAM_READWRITE));
2646
2647 g_object_class_install_property (object_class, SOL_PROP_FLIP_SIGN,
2648 g_param_spec_boolean ("flip-sign",
2649 P_("Flip Sign"),
2650 P_("Flip sign of target value"),
2651 FALSE,
2652 GSF_PARAM_STATIC | G_PARAM_READWRITE));
2653
2654 solver_signals[SOL_SIG_PREPARE] =
2655 g_signal_new ("prepare",
2656 G_OBJECT_CLASS_TYPE (object_class),
2657 G_SIGNAL_RUN_LAST,
2658 G_STRUCT_OFFSET (GnmSolverClass, prepare),
2659 NULL, NULL,
2660 gnm__BOOLEAN__OBJECT_POINTER,
2661 G_TYPE_BOOLEAN, 2,
2662 G_TYPE_OBJECT,
2663 G_TYPE_POINTER);
2664
2665 solver_signals[SOL_SIG_START] =
2666 g_signal_new ("start",
2667 G_OBJECT_CLASS_TYPE (object_class),
2668 G_SIGNAL_RUN_LAST,
2669 G_STRUCT_OFFSET (GnmSolverClass, start),
2670 NULL, NULL,
2671 gnm__BOOLEAN__OBJECT_POINTER,
2672 G_TYPE_BOOLEAN, 2,
2673 G_TYPE_OBJECT,
2674 G_TYPE_POINTER);
2675
2676 solver_signals[SOL_SIG_STOP] =
2677 g_signal_new ("stop",
2678 G_OBJECT_CLASS_TYPE (object_class),
2679 G_SIGNAL_RUN_LAST,
2680 G_STRUCT_OFFSET (GnmSolverClass, stop),
2681 NULL, NULL,
2682 gnm__BOOLEAN__POINTER,
2683 G_TYPE_BOOLEAN, 1,
2684 G_TYPE_POINTER);
2685 }
2686
2687 GSF_CLASS (GnmSolver, gnm_solver,
2688 &gnm_solver_class_init, NULL, G_TYPE_OBJECT)
2689
2690 /* ------------------------------------------------------------------------- */
2691
2692 static GObjectClass *gnm_solver_result_parent_class;
2693
2694 static void
2695 gnm_solver_result_finalize (GObject *obj)
2696 {
2697 GnmSolverResult *r = GNM_SOLVER_RESULT (obj);
2698 g_free (r->solution);
2699 gnm_solver_result_parent_class->finalize (obj);
2700 }
2701
2702 static void
2703 gnm_solver_result_class_init (GObjectClass *object_class)
2704 {
2705 gnm_solver_result_parent_class =
2706 g_type_class_peek_parent (object_class);
2707
2708 object_class->finalize = gnm_solver_result_finalize;
2709 }
2710
2711 GSF_CLASS (GnmSolverResult, gnm_solver_result,
2712 &gnm_solver_result_class_init, NULL, G_TYPE_OBJECT)
2713
2714 /* ------------------------------------------------------------------------- */
2715
2716 static GObjectClass *gnm_solver_sensitivity_parent_class;
2717
2718 enum {
2719 SOLS_PROP_0,
2720 SOLS_PROP_SOLVER
2721 };
2722
2723 static void
2724 gnm_solver_sensitivity_constructed (GObject *obj)
2725 {
2726 GnmSolverSensitivity *sols = GNM_SOLVER_SENSITIVITY (obj);
2727 GnmSolver *sol = sols->solver;
2728 GnmSolverParameters *sp = sol->params;
2729 const int n = sol->input_cells->len;
2730 int i, cn;
2731 GSList *l;
2732
2733 /* Chain to parent first */
2734 gnm_solver_sensitivity_parent_class->constructed (obj);
2735
2736 sols->vars = g_new (struct GnmSolverSensitivityVars_, n);
2737 for (i = 0; i < n; i++) {
2738 sols->vars[i].low = gnm_nan;
2739 sols->vars[i].high = gnm_nan;
2740 sols->vars[i].reduced_cost = gnm_nan;
2741 }
2742
2743 cn = 0;
2744 for (l = sp->constraints; l; l = l->next) {
2745 GnmSolverConstraint *c = l->data;
2746 int i;
2747 gnm_float cl, cr;
2748 GnmCell *lhs, *rhs;
2749
2750 for (i = 0;
2751 gnm_solver_constraint_get_part (c, sp, i,
2752 &lhs, &cl,
2753 &rhs, &cr);
2754 i++) {
2755 cn++;
2756 }
2757 }
2758 sols->constraints = g_new (struct GnmSolverSensitivityConstraints_, cn);
2759 for (i = 0; i < cn; i++) {
2760 sols->constraints[i].low = gnm_nan;
2761 sols->constraints[i].high = gnm_nan;
2762 sols->constraints[i].shadow_price = gnm_nan;
2763 }
2764 }
2765
2766 static void
2767 gnm_solver_sensitivity_finalize (GObject *obj)
2768 {
2769 GnmSolverSensitivity *r = GNM_SOLVER_SENSITIVITY (obj);
2770 g_free (r->vars);
2771 g_free (r->constraints);
2772 gnm_solver_sensitivity_parent_class->finalize (obj);
2773 }
2774
2775 static void
2776 gnm_solver_sensitivity_get_property (GObject *object, guint property_id,
2777 GValue *value, GParamSpec *pspec)
2778 {
2779 GnmSolverSensitivity *sols = (GnmSolverSensitivity *)object;
2780
2781 switch (property_id) {
2782 case SOLS_PROP_SOLVER:
2783 g_value_set_object (value, sols->solver);
2784 break;
2785
2786 default:
2787 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
2788 break;
2789 }
2790 }
2791
2792 static void
2793 gnm_solver_sensitivity_set_property (GObject *object, guint property_id,
2794 GValue const *value, GParamSpec *pspec)
2795 {
2796 GnmSolverSensitivity *sols = (GnmSolverSensitivity *)object;
2797
2798 switch (property_id) {
2799 case SOLS_PROP_SOLVER:
2800 /* We hold no ref. */
2801 sols->solver = g_value_get_object (value);
2802 break;
2803
2804 default:
2805 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
2806 break;
2807 }
2808 }
2809
2810 static void
2811 gnm_solver_sensitivity_class_init (GObjectClass *object_class)
2812 {
2813 gnm_solver_sensitivity_parent_class =
2814 g_type_class_peek_parent (object_class);
2815
2816 object_class->finalize = gnm_solver_sensitivity_finalize;
2817 object_class->constructed = gnm_solver_sensitivity_constructed;
2818 object_class->set_property = gnm_solver_sensitivity_set_property;
2819 object_class->get_property = gnm_solver_sensitivity_get_property;
2820
2821 g_object_class_install_property
2822 (object_class, SOLS_PROP_SOLVER,
2823 g_param_spec_object ("solver",
2824 P_("Solver"),
2825 P_("Solver"),
2826 GNM_SOLVER_TYPE,
2827 GSF_PARAM_STATIC |
2828 G_PARAM_CONSTRUCT_ONLY |
2829 G_PARAM_READWRITE));
2830 }
2831
2832 GSF_CLASS (GnmSolverSensitivity, gnm_solver_sensitivity,
2833 gnm_solver_sensitivity_class_init, NULL, G_TYPE_OBJECT)
2834
2835 GnmSolverSensitivity *
2836 gnm_solver_sensitivity_new (GnmSolver *sol)
2837 {
2838 return g_object_new (GNM_SOLVER_SENSITIVITY_TYPE, "solver", sol, NULL);
2839 }
2840
2841 /* ------------------------------------------------------------------------- */
2842
2843 static GObjectClass *gnm_sub_solver_parent_class;
2844
2845 enum {
2846 SUB_SOL_SIG_CHILD_EXIT,
2847 SUB_SOL_SIG_LAST
2848 };
2849
2850 static guint sub_solver_signals[SUB_SOL_SIG_LAST] = { 0 };
2851
2852 void
2853 gnm_sub_solver_clear (GnmSubSolver *subsol)
2854 {
2855 int i;
2856
2857 if (subsol->child_watch) {
2858 g_source_remove (subsol->child_watch);
2859 subsol->child_watch = 0;
2860 }
2861
2862 if (subsol->child_pid) {
2863 #ifdef G_OS_WIN32
2864 TerminateProcess (subsol->child_pid, 127);
2865 #else
2866 kill (subsol->child_pid, SIGKILL);
2867 #endif
2868 g_spawn_close_pid (subsol->child_pid);
2869 subsol->child_pid = (GPid)0;
2870 }
2871
2872 for (i = 0; i <= 2; i++) {
2873 if (subsol->channel_watches[i]) {
2874 g_source_remove (subsol->channel_watches[i]);
2875 subsol->channel_watches[i] = 0;
2876 }
2877 if (subsol->channels[i]) {
2878 g_io_channel_unref (subsol->channels[i]);
2879 subsol->channels[i] = NULL;
2880 }
2881 if (subsol->fd[i] != -1) {
2882 close (subsol->fd[i]);
2883 subsol->fd[i] = -1;
2884 }
2885 }
2886
2887 if (subsol->program_filename) {
2888 g_unlink (subsol->program_filename);
2889 g_free (subsol->program_filename);
2890 subsol->program_filename = NULL;
2891 }
2892
2893 if (subsol->cell_from_name)
2894 g_hash_table_remove_all (subsol->cell_from_name);
2895
2896 if (subsol->name_from_cell)
2897 g_hash_table_remove_all (subsol->name_from_cell);
2898
2899 if (subsol->constraint_from_name)
2900 g_hash_table_remove_all (subsol->constraint_from_name);
2901 }
2902
2903 static void
2904 gnm_sub_solver_dispose (GObject *obj)
2905 {
2906 GnmSubSolver *subsol = GNM_SUB_SOLVER (obj);
2907
2908 gnm_sub_solver_clear (subsol);
2909
2910 gnm_sub_solver_parent_class->dispose (obj);
2911 }
2912
2913 static void
2914 gnm_sub_solver_finalize (GObject *obj)
2915 {
2916 GnmSubSolver *subsol = GNM_SUB_SOLVER (obj);
2917
2918 /*
2919 * The weird finalization in gnm_lpsolve_final makes it important that
2920 * we leave the object in a state that gnm_sub_solver_clear is happy
2921 * with.
2922 */
2923
2924 g_hash_table_destroy (subsol->cell_from_name);
2925 subsol->cell_from_name = NULL;
2926
2927 g_hash_table_destroy (subsol->name_from_cell);
2928 subsol->name_from_cell = NULL;
2929
2930 g_hash_table_destroy (subsol->constraint_from_name);
2931 subsol->constraint_from_name = NULL;
2932
2933 gnm_sub_solver_parent_class->finalize (obj);
2934 }
2935
2936 static void
2937 gnm_sub_solver_init (GnmSubSolver *subsol)
2938 {
2939 int i;
2940
2941 for (i = 0; i <= 2; i++)
2942 subsol->fd[i] = -1;
2943
2944 subsol->cell_from_name =
2945 g_hash_table_new_full (g_str_hash, g_str_equal,
2946 g_free, NULL);
2947 subsol->name_from_cell =
2948 g_hash_table_new (g_direct_hash, g_direct_equal);
2949
2950 subsol->constraint_from_name =
2951 g_hash_table_new_full (g_str_hash, g_str_equal,
2952 g_free, NULL);
2953 }
2954
2955 static void
2956 cb_child_exit (G_GNUC_UNUSED GPid pid, gint status, GnmSubSolver *subsol)
2957 {
2958 gboolean normal = WIFEXITED (status);
2959 int code;
2960
2961 subsol->child_watch = 0;
2962
2963 if (normal) {
2964 code = WEXITSTATUS (status);
2965 if (gnm_solver_debug ())
2966 g_printerr ("Solver process exited with code %d\n",
2967 code);
2968 #ifndef G_OS_WIN32
2969 } else if (WIFSIGNALED (status)) {
2970 code = WTERMSIG (status);
2971 if (gnm_solver_debug ())
2972 g_printerr ("Solver process received signal %d\n",
2973 code);
2974 #endif
2975 } else {
2976 code = -1;
2977 g_printerr ("Solver process exited with status 0x%x\n",
2978 status);
2979 }
2980
2981 g_signal_emit (subsol, sub_solver_signals[SUB_SOL_SIG_CHILD_EXIT], 0,
2982 normal, code);
2983
2984 if (subsol->child_pid) {
2985 g_spawn_close_pid (subsol->child_pid);
2986 subsol->child_pid = (GPid)0;
2987 }
2988 }
2989
2990 /**
2991 * gnm_sub_solver_spawn: (skip)
2992 */
2993 gboolean
2994 gnm_sub_solver_spawn (GnmSubSolver *subsol,
2995 char **argv,
2996 GSpawnChildSetupFunc child_setup, gpointer setup_data,
2997 GIOFunc io_stdout, gpointer stdout_data,
2998 GIOFunc io_stderr, gpointer stderr_data,
2999 GError **err)
3000 {
3001 GnmSolver *sol = GNM_SOLVER (subsol);
3002 gboolean ok;
3003 GSpawnFlags spflags = G_SPAWN_DO_NOT_REAP_CHILD;
3004 int fd;
3005
3006 g_return_val_if_fail (subsol->child_watch == 0, FALSE);
3007 g_return_val_if_fail (sol->status == GNM_SOLVER_STATUS_PREPARED, FALSE);
3008
3009 if (!g_path_is_absolute (argv[0]))
3010 spflags |= G_SPAWN_SEARCH_PATH;
3011
3012 if (io_stdout == NULL && !gnm_solver_debug ())
3013 spflags |= G_SPAWN_STDOUT_TO_DEV_NULL;
3014
3015 if (gnm_solver_debug ()) {
3016 GString *msg = g_string_new ("Spawning");
3017 int i;
3018 for (i = 0; argv[i]; i++) {
3019 g_string_append_c (msg, ' ');
3020 g_string_append (msg, argv[i]);
3021 }
3022 g_printerr ("%s\n", msg->str);
3023 g_string_free (msg, TRUE);
3024 }
3025
3026 #ifdef G_OS_WIN32
3027 /* Hope for the best... */
3028 child_setup = NULL;
3029 setup_data = NULL;
3030 #endif
3031
3032 ok = g_spawn_async_with_pipes
3033 (g_get_home_dir (), /* PWD */
3034 argv,
3035 NULL, /* environment */
3036 spflags,
3037 child_setup, setup_data,
3038 &subsol->child_pid,
3039 NULL, /* stdin */
3040 io_stdout ? &subsol->fd[1] : NULL, /* stdout */
3041 io_stdout ? &subsol->fd[2] : NULL, /* stderr */
3042 err);
3043 if (!ok)
3044 goto fail;
3045
3046 subsol->child_watch =
3047 g_child_watch_add (subsol->child_pid,
3048 (GChildWatchFunc)cb_child_exit, subsol);
3049
3050 subsol->io_funcs[1] = io_stdout;
3051 subsol->io_funcs_data[1] = stdout_data;
3052 subsol->io_funcs[2] = io_stderr;
3053 subsol->io_funcs_data[2] = stderr_data;
3054
3055 for (fd = 1; fd <= 2; fd++) {
3056 GIOFlags ioflags;
3057
3058 if (subsol->io_funcs[fd] == NULL)
3059 continue;
3060
3061 /*
3062 * Despite the name these are documented to work on Win32.
3063 * Let us hope that is actually true.
3064 */
3065 subsol->channels[fd] = g_io_channel_unix_new (subsol->fd[fd]);
3066 ioflags = g_io_channel_get_flags (subsol->channels[fd]);
3067 g_io_channel_set_flags (subsol->channels[fd],
3068 ioflags | G_IO_FLAG_NONBLOCK,
3069 NULL);
3070 subsol->channel_watches[fd] =
3071 g_io_add_watch (subsol->channels[fd],
3072 G_IO_IN,
3073 subsol->io_funcs[fd],
3074 subsol->io_funcs_data[fd]);
3075 }
3076
3077 gnm_solver_set_status (sol, GNM_SOLVER_STATUS_RUNNING);
3078 return TRUE;
3079
3080 fail:
3081 gnm_sub_solver_clear (subsol);
3082 gnm_solver_set_status (sol, GNM_SOLVER_STATUS_ERROR);
3083 return FALSE;
3084 }
3085
3086 const char *
3087 gnm_sub_solver_name_cell (GnmSubSolver *subsol, GnmCell const *cell,
3088 const char *name)
3089 {
3090 char *name_copy = g_strdup (name);
3091
3092 g_hash_table_insert (subsol->cell_from_name,
3093 name_copy,
3094 (gpointer)cell);
3095 g_hash_table_insert (subsol->name_from_cell,
3096 (gpointer)cell,
3097 name_copy);
3098
3099 return name_copy;
3100 }
3101
3102 GnmCell *
3103 gnm_sub_solver_find_cell (GnmSubSolver *subsol, const char *name)
3104 {
3105 return g_hash_table_lookup (subsol->cell_from_name, name);
3106 }
3107
3108 const char *
3109 gnm_sub_solver_get_cell_name (GnmSubSolver *subsol,
3110 GnmCell const *cell)
3111 {
3112 return g_hash_table_lookup (subsol->name_from_cell, (gpointer)cell);
3113 }
3114
3115 const char *
3116 gnm_sub_solver_name_constraint (GnmSubSolver *subsol,
3117 int cidx,
3118 const char *name)
3119 {
3120 char *name_copy = g_strdup (name);
3121
3122 g_hash_table_insert (subsol->constraint_from_name,
3123 name_copy,
3124 GINT_TO_POINTER (cidx));
3125
3126 return name_copy;
3127 }
3128
3129 int
3130 gnm_sub_solver_find_constraint (GnmSubSolver *subsol, const char *name)
3131 {
3132 gpointer idx;
3133
3134 if (g_hash_table_lookup_extended (subsol->constraint_from_name,
3135 (gpointer)name,
3136 NULL, &idx))
3137 return GPOINTER_TO_INT (idx);
3138 else
3139 return -1;
3140 }
3141
3142
3143 char *
3144 gnm_sub_solver_locate_binary (const char *binary, const char *solver,
3145 const char *url,
3146 WBCGtk *wbcg)
3147 {
3148 GtkWindow *parent;
3149 GtkWidget *dialog;
3150 char *path = NULL;
3151 int res;
3152 GtkFileChooser *fsel;
3153 char *title;
3154
3155 parent = wbcg ? wbcg_toplevel (wbcg) : NULL;
3156 dialog = gtk_message_dialog_new_with_markup
3157 (parent,
3158 GTK_DIALOG_DESTROY_WITH_PARENT,
3159 GTK_MESSAGE_QUESTION,
3160 GTK_BUTTONS_YES_NO,
3161 _("Gnumeric is unable to locate the program <i>%s</i> needed "
3162 "for the <i>%s</i> solver. For more information see %s.\n\n"
3163 "Would you like to locate it yourself?"),
3164 binary, solver, url);
3165 title = g_strdup_printf (_("Unable to locate %s"), binary);
3166 g_object_set (G_OBJECT (dialog),
3167 "title", title,
3168 NULL);
3169 g_free (title);
3170
3171 res = go_gtk_dialog_run (GTK_DIALOG (dialog), parent);
3172 switch (res) {
3173 case GTK_RESPONSE_NO:
3174 case GTK_RESPONSE_DELETE_EVENT:
3175 default:
3176 return NULL;
3177 case GTK_RESPONSE_YES:
3178 break;
3179 }
3180
3181 title = g_strdup_printf (_("Locate the %s program"), binary);
3182 fsel = GTK_FILE_CHOOSER
3183 (g_object_new (GTK_TYPE_FILE_CHOOSER_DIALOG,
3184 "action", GTK_FILE_CHOOSER_ACTION_OPEN,
3185 "local-only", TRUE,
3186 "title", title,
3187 NULL));
3188 g_free (title);
3189 go_gtk_dialog_add_button (GTK_DIALOG (fsel), GNM_STOCK_CANCEL,
3190 "gtk-cancel", GTK_RESPONSE_CANCEL);
3191 go_gtk_dialog_add_button (GTK_DIALOG (fsel), GNM_STOCK_OK,
3192 "system-run", GTK_RESPONSE_OK);
3193 g_object_ref (fsel);
3194 if (go_gtk_file_sel_dialog (parent, GTK_WIDGET (fsel))) {
3195 path = gtk_file_chooser_get_filename (fsel);
3196 if (!g_file_test (path, G_FILE_TEST_IS_EXECUTABLE)) {
3197 g_free (path);
3198 path = NULL;
3199 }
3200 }
3201
3202 gtk_widget_destroy (GTK_WIDGET (fsel));
3203 g_object_unref (fsel);
3204
3205 return path;
3206 }
3207
3208
3209 void
3210 gnm_sub_solver_flush (GnmSubSolver *subsol)
3211 {
3212 int fd;
3213
3214 for (fd = 1; fd <= 2; fd++) {
3215 if (subsol->io_funcs[fd] == NULL)
3216 continue;
3217
3218 subsol->io_funcs[fd] (subsol->channels[fd],
3219 G_IO_IN,
3220 subsol->io_funcs_data[fd]);
3221 }
3222 }
3223
3224 static void
3225 gnm_sub_solver_class_init (GObjectClass *object_class)
3226 {
3227 gnm_sub_solver_parent_class = g_type_class_peek_parent (object_class);
3228
3229 object_class->dispose = gnm_sub_solver_dispose;
3230 object_class->finalize = gnm_sub_solver_finalize;
3231
3232 sub_solver_signals[SUB_SOL_SIG_CHILD_EXIT] =
3233 g_signal_new ("child-exit",
3234 G_OBJECT_CLASS_TYPE (object_class),
3235 G_SIGNAL_RUN_LAST,
3236 G_STRUCT_OFFSET (GnmSubSolverClass, child_exit),
3237 NULL, NULL,
3238 gnm__VOID__BOOLEAN_INT,
3239 G_TYPE_NONE, 2,
3240 G_TYPE_BOOLEAN, G_TYPE_INT);
3241 }
3242
3243 GSF_CLASS (GnmSubSolver, gnm_sub_solver,
3244 gnm_sub_solver_class_init, gnm_sub_solver_init, GNM_SOLVER_TYPE)
3245
3246 /* ------------------------------------------------------------------------- */
3247
3248 enum {
3249 SOL_ITER_SIG_ITERATE,
3250 SOL_ITER_SIG_LAST
3251 };
3252
3253 static guint solver_iterator_signals[SOL_ITER_SIG_LAST] = { 0 };
3254
3255 static void
3256 gnm_solver_iterator_class_init (GObjectClass *object_class)
3257 {
3258 solver_iterator_signals[SOL_ITER_SIG_ITERATE] =
3259 g_signal_new ("iterate",
3260 G_OBJECT_CLASS_TYPE (object_class),
3261 G_SIGNAL_RUN_LAST,
3262 G_STRUCT_OFFSET (GnmSolverIteratorClass, iterate),
3263 NULL, NULL,
3264 gnm__BOOLEAN__VOID,
3265 G_TYPE_BOOLEAN, 0);
3266 }
3267
3268 gboolean
3269 gnm_solver_iterator_iterate (GnmSolverIterator *iter)
3270 {
3271 gboolean progress = FALSE;
3272 g_signal_emit (iter, solver_iterator_signals[SOL_ITER_SIG_ITERATE], 0, &progress);
3273 return progress;
3274 }
3275
3276 /**
3277 * gnm_solver_iterator_new_func: (skip)
3278 */
3279 GnmSolverIterator *
3280 gnm_solver_iterator_new_func (GCallback iterate, gpointer user)
3281 {
3282 GnmSolverIterator *iter;
3283
3284 iter = g_object_new (GNM_SOLVER_ITERATOR_TYPE, NULL);
3285 g_signal_connect (iter, "iterate", G_CALLBACK (iterate), user);
3286 return iter;
3287 }
3288
3289 static gboolean
3290 cb_polish_iter (GnmSolverIterator *iter, GnmIterSolver *isol)
3291 {
3292 GnmSolver *sol = GNM_SOLVER (isol);
3293 const int n = sol->input_cells->len;
3294 gnm_float *dir;
3295 gboolean progress = FALSE;
3296 int c;
3297
3298 dir = g_new0 (gnm_float, n);
3299 for (c = 0; c < n; c++) {
3300 gnm_float s, y, s0, sm, xc = isol->xk[c];
3301
3302 if (xc == 0) {
3303 s0 = 0.5;
3304 sm = 1;
3305 } else {
3306 int e;
3307 (void)gnm_frexp (xc, &e);
3308 s0 = gnm_ldexp (1, e - 10);
3309 if (s0 == 0) s0 = GNM_MIN;
3310 sm = gnm_abs (xc);
3311 }
3312
3313 dir[c] = 1;
3314 s = gnm_solver_line_search (sol, isol->xk, dir, TRUE,
3315 s0, sm, 0.0, &y);
3316 dir[c] = 0;
3317
3318 if (gnm_finite (s) && s != 0) {
3319 isol->xk[c] += s;
3320 isol->yk = y;
3321 progress = TRUE;
3322 }
3323 }
3324 g_free (dir);
3325
3326 if (progress)
3327 gnm_iter_solver_set_solution (isol);
3328
3329 return progress;
3330 }
3331
3332 /**
3333 * gnm_solver_iterator_new_polish:
3334 * @isol: the solver to operate on
3335 *
3336 * Returns: (transfer full): an iterator object that can be used to polish
3337 * a solution by simple axis-parallel movement.
3338 */
3339 GnmSolverIterator *
3340 gnm_solver_iterator_new_polish (GnmIterSolver *isol)
3341 {
3342 return gnm_solver_iterator_new_func (G_CALLBACK (cb_polish_iter), isol);
3343 }
3344
3345
3346 static gboolean
3347 cb_gradient_iter (GnmSolverIterator *iter, GnmIterSolver *isol)
3348 {
3349 GnmSolver *sol = GNM_SOLVER (isol);
3350 const int n = sol->input_cells->len;
3351 gboolean progress = FALSE;
3352 gnm_float s, y;
3353 gnm_float *g;
3354 int i;
3355
3356 /* Search in opposite direction of gradient. */
3357 g = gnm_solver_compute_gradient (sol, isol->xk);
3358 for (i = 0; i < n; i++)
3359 g[i] = -g[i];
3360
3361 s = gnm_solver_line_search (sol, isol->xk, g, FALSE,
3362 1, gnm_pinf, 0.0, &y);
3363 if (s > 0) {
3364 for (i = 0; i < n; i++)
3365 isol->xk[i] += s * g[i];
3366 isol->yk = y;
3367 progress = TRUE;
3368 }
3369
3370 g_free (g);
3371
3372 if (progress)
3373 gnm_iter_solver_set_solution (isol);
3374
3375 return progress;
3376 }
3377
3378 /**
3379 * gnm_solver_iterator_new_gradient:
3380 * @isol: the solver to operate on
3381 *
3382 * Returns: (transfer full): an iterator object that can be used to perform
3383 * a gradient descent step.
3384 */
3385 GnmSolverIterator *
3386 gnm_solver_iterator_new_gradient (GnmIterSolver *isol)
3387 {
3388 return gnm_solver_iterator_new_func (G_CALLBACK (cb_gradient_iter), isol);
3389 }
3390
3391 GSF_CLASS (GnmSolverIterator, gnm_solver_iterator,
3392 gnm_solver_iterator_class_init, NULL, G_TYPE_OBJECT)
3393
3394 /* ------------------------------------------------------------------------- */
3395
3396 enum {
3397 SOLIC_PROP_0,
3398 SOLIC_PROP_CYCLES
3399 };
3400
3401 static GObjectClass *gnm_solver_iterator_compound_parent_class;
3402
3403 /**
3404 * gnm_solver_iterator_compound_add:
3405 * @ic: Compound iterator
3406 * @iter: (transfer full): sub-iterator
3407 * @count: repeat count
3408 *
3409 * Add an iterator to a compound iterator with a given repeat count. As a
3410 * special case, a repeat count of zero means to try the iterator once
3411 * in a cycle, but only if no other sub-iterator has shown any progress so far.
3412 */
3413 void
3414 gnm_solver_iterator_compound_add (GnmSolverIteratorCompound *ic,
3415 GnmSolverIterator *iter,
3416 unsigned count)
3417 {
3418 g_ptr_array_add (ic->iterators, iter);
3419 ic->counts = g_renew (unsigned, ic->counts, ic->iterators->len);
3420 ic->counts[ic->iterators->len - 1] = count;
3421 }
3422
3423 static gboolean
3424 gnm_solver_iterator_compound_iterate (GnmSolverIterator *iter)
3425 {
3426 GnmSolverIteratorCompound *ic = (GnmSolverIteratorCompound *)iter;
3427 gboolean progress;
3428
3429 while (TRUE) {
3430 if (ic->cycle >= ic->cycles)
3431 return FALSE;
3432
3433 if (ic->next >= ic->iterators->len) {
3434 /* We've been through all iterators. */
3435 if (!ic->cycle_progress)
3436 return FALSE;
3437 ic->cycle_progress = FALSE;
3438 ic->next = 0;
3439 ic->next_counter = 0;
3440 ic->cycle++;
3441 continue;
3442 }
3443
3444 if (ic->next_counter < ic->counts[ic->next])
3445 break;
3446
3447 /* Special case: when count==0, use only if no progress. */
3448 if (!ic->cycle_progress && ic->next_counter == 0)
3449 break;
3450
3451 ic->next++;
3452 ic->next_counter = 0;
3453 }
3454
3455 progress = gnm_solver_iterator_iterate (g_ptr_array_index (ic->iterators, ic->next));
3456 if (progress) {
3457 ic->cycle_progress = TRUE;
3458 ic->next_counter++;
3459 } else {
3460 /* No progress, so don't retry. */
3461 ic->next++;
3462 ic->next_counter = 0;
3463 }
3464
3465 /* Report progress as long as we have stuff to try. */
3466 return TRUE;
3467 }
3468
3469 static void
3470 gnm_solver_iterator_compound_init (GnmSolverIteratorCompound *ic)
3471 {
3472 ic->iterators = g_ptr_array_new ();
3473 ic->cycles = G_MAXUINT;
3474 }
3475
3476 static void
3477 gnm_solver_iterator_compound_finalize (GObject *obj)
3478 {
3479 GnmSolverIteratorCompound *ic = (GnmSolverIteratorCompound *)obj;
3480 g_ptr_array_foreach (ic->iterators, (GFunc)g_object_unref, NULL);
3481 g_ptr_array_free (ic->iterators, TRUE);
3482 g_free (ic->counts);
3483 gnm_solver_iterator_compound_parent_class->finalize (obj);
3484 }
3485
3486 static void
3487 gnm_solver_iterator_compound_get_property (GObject *object, guint property_id,
3488 GValue *value, GParamSpec *pspec)
3489 {
3490 GnmSolverIteratorCompound *it = (GnmSolverIteratorCompound *)object;
3491
3492 switch (property_id) {
3493 case SOLIC_PROP_CYCLES:
3494 g_value_set_uint (value, it->cycles);
3495 break;
3496
3497 default:
3498 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
3499 break;
3500 }
3501 }
3502
3503 static void
3504 gnm_solver_iterator_compound_set_property (GObject *object, guint property_id,
3505 GValue const *value, GParamSpec *pspec)
3506 {
3507 GnmSolverIteratorCompound *it = (GnmSolverIteratorCompound *)object;
3508
3509 switch (property_id) {
3510 case SOLIC_PROP_CYCLES:
3511 it->cycles = g_value_get_uint (value);
3512 break;
3513
3514 default:
3515 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
3516 break;
3517 }
3518 }
3519
3520 static void
3521 gnm_solver_iterator_compound_class_init (GObjectClass *object_class)
3522 {
3523 GnmSolverIteratorClass *iclass = (GnmSolverIteratorClass *)object_class;
3524
3525 gnm_solver_iterator_compound_parent_class = g_type_class_peek_parent (object_class);
3526
3527 object_class->finalize = gnm_solver_iterator_compound_finalize;
3528 object_class->set_property = gnm_solver_iterator_compound_set_property;
3529 object_class->get_property = gnm_solver_iterator_compound_get_property;
3530 iclass->iterate = gnm_solver_iterator_compound_iterate;
3531
3532 g_object_class_install_property (object_class, SOLIC_PROP_CYCLES,
3533 g_param_spec_uint ("cycles",
3534 P_("Cycles"),
3535 P_("Maximum number of cycles"),
3536 0, G_MAXUINT, G_MAXUINT,
3537 GSF_PARAM_STATIC | G_PARAM_READWRITE));
3538 }
3539
3540 GSF_CLASS (GnmSolverIteratorCompound, gnm_solver_iterator_compound,
3541 gnm_solver_iterator_compound_class_init, gnm_solver_iterator_compound_init, GNM_SOLVER_ITERATOR_TYPE)
3542
3543 /* ------------------------------------------------------------------------- */
3544
3545 static GObjectClass *gnm_iter_solver_parent_class;
3546
3547 void
3548 gnm_iter_solver_set_iterator (GnmIterSolver *isol, GnmSolverIterator *iterator)
3549 {
3550 GnmSolverIterator *old_iterator;
3551
3552 g_return_if_fail (GNM_IS_ITER_SOLVER (isol));
3553
3554 old_iterator = isol->iterator;
3555 isol->iterator = iterator ? g_object_ref (iterator) : NULL;
3556 if (old_iterator)
3557 g_object_unref (old_iterator);
3558 }
3559
3560 gboolean
3561 gnm_iter_solver_get_initial_solution (GnmIterSolver *isol, GError **err)
3562 {
3563 GnmSolver *sol = GNM_SOLVER (isol);
3564 const int n = sol->input_cells->len;
3565 int i;
3566
3567 if (gnm_solver_check_constraints (sol))
3568 goto got_it;
3569
3570 /* More? */
3571
3572 g_set_error (err,
3573 go_error_invalid (),
3574 0,
3575 _("The initial values do not satisfy the constraints."));
3576 return FALSE;
3577
3578 got_it:
3579 for (i = 0; i < n; i++) {
3580 GnmCell *cell = g_ptr_array_index (sol->input_cells, i);
3581 isol->xk[i] = value_get_as_float (cell->value);
3582 }
3583 isol->yk = gnm_solver_get_target_value (sol);
3584
3585 gnm_iter_solver_set_solution (isol);
3586
3587 return TRUE;
3588 }
3589
3590 void
3591 gnm_iter_solver_set_solution (GnmIterSolver *isol)
3592 {
3593 GnmSolver *sol = GNM_SOLVER (isol);
3594 GnmSolverResult *result = g_object_new (GNM_SOLVER_RESULT_TYPE, NULL);
3595 const int n = sol->input_cells->len;
3596
3597 result->quality = GNM_SOLVER_RESULT_FEASIBLE;
3598 result->value = sol->flip_sign ? 0 - isol->yk : isol->yk;
3599 result->solution = g_memdup (isol->xk, n * sizeof (gnm_float));
3600 g_object_set (sol, "result", result, NULL);
3601 g_object_unref (result);
3602
3603 if (!gnm_solver_check_constraints (sol)) {
3604 g_printerr ("Infeasible solution set\n");
3605 }
3606 }
3607
3608 static void
3609 gnm_iter_solver_clear (GnmIterSolver *isol)
3610 {
3611 if (isol->idle_tag) {
3612 g_source_remove (isol->idle_tag);
3613 isol->idle_tag = 0;
3614 }
3615 }
3616
3617 static void
3618 gnm_iter_solver_dispose (GObject *obj)
3619 {
3620 GnmIterSolver *isol = GNM_ITER_SOLVER (obj);
3621 gnm_iter_solver_clear (isol);
3622 gnm_iter_solver_parent_class->dispose (obj);
3623 }
3624
3625 static void
3626 gnm_iter_solver_finalize (GObject *obj)
3627 {
3628 GnmIterSolver *isol = GNM_ITER_SOLVER (obj);
3629 g_free (isol->xk);
3630 gnm_iter_solver_parent_class->finalize (obj);
3631 }
3632
3633 static void
3634 gnm_iter_solver_constructed (GObject *obj)
3635 {
3636 GnmIterSolver *isol = GNM_ITER_SOLVER (obj);
3637 GnmSolver *sol = GNM_SOLVER (obj);
3638
3639 /* Chain to parent first */
3640 gnm_iter_solver_parent_class->constructed (obj);
3641
3642 isol->xk = g_new0 (gnm_float, sol->input_cells->len);
3643 }
3644
3645 static void
3646 gnm_iter_solver_init (GnmIterSolver *isol)
3647 {
3648 }
3649
3650 static gint
3651 gnm_iter_solver_idle (gpointer data)
3652 {
3653 GnmIterSolver *isol = data;
3654 GnmSolver *sol = &isol->parent;
3655 GnmSolverParameters *params = sol->params;
3656 gboolean progress;
3657
3658 progress = isol->iterator && gnm_solver_iterator_iterate (isol->iterator);
3659 isol->iterations++;
3660
3661 if (!gnm_solver_finished (sol)) {
3662 if (!progress) {
3663 gnm_solver_set_status (sol, GNM_SOLVER_STATUS_DONE);
3664 } else if (isol->iterations >= params->options.max_iter) {
3665 gnm_solver_stop (sol, NULL);
3666 gnm_solver_set_reason (sol, _("Iteration limit exceeded"));
3667 }
3668 }
3669
3670 if (gnm_solver_finished (sol)) {
3671 isol->idle_tag = 0;
3672
3673 gnm_app_recalc ();
3674
3675 return FALSE;
3676 } else {
3677 /* Call again. */
3678 return TRUE;
3679 }
3680 }
3681
3682 static gboolean
3683 gnm_iter_solver_start (GnmSolver *solver, WorkbookControl *wbc, GError **err)
3684 {
3685 GnmIterSolver *isol = GNM_ITER_SOLVER (solver);
3686
3687 g_return_val_if_fail (isol->idle_tag == 0, FALSE);
3688
3689 isol->idle_tag = g_idle_add (gnm_iter_solver_idle, solver);
3690 gnm_solver_set_status (solver, GNM_SOLVER_STATUS_RUNNING);
3691
3692 return TRUE;
3693 }
3694
3695 static gboolean
3696 gnm_iter_solver_stop (GnmSolver *solver, GError **err)
3697 {
3698 GnmIterSolver *isol = GNM_ITER_SOLVER (solver);
3699 GnmSolver *sol = &isol->parent;
3700
3701 gnm_iter_solver_clear (isol);
3702
3703 g_clear_object (&isol->iterator);
3704
3705 gnm_solver_set_status (sol, GNM_SOLVER_STATUS_CANCELLED);
3706
3707 return TRUE;
3708 }
3709
3710 static void
3711 gnm_iter_solver_class_init (GObjectClass *object_class)
3712 {
3713 GnmSolverClass *sclass = (GnmSolverClass *)object_class;
3714
3715 gnm_iter_solver_parent_class = g_type_class_peek_parent (object_class);
3716
3717 object_class->dispose = gnm_iter_solver_dispose;
3718 object_class->finalize = gnm_iter_solver_finalize;
3719 object_class->constructed = gnm_iter_solver_constructed;
3720 sclass->start = gnm_iter_solver_start;
3721 sclass->stop = gnm_iter_solver_stop;
3722 }
3723
3724 GSF_CLASS (GnmIterSolver, gnm_iter_solver,
3725 gnm_iter_solver_class_init, gnm_iter_solver_init, GNM_SOLVER_TYPE)
3726
3727 /* ------------------------------------------------------------------------- */
3728
3729 static GObjectClass *gnm_solver_factory_parent_class;
3730
3731 static void
3732 gnm_solver_factory_finalize (GObject *obj)
3733 {
3734 GnmSolverFactory *factory = GNM_SOLVER_FACTORY (obj);
3735
3736 if (factory->notify)
3737 factory->notify (factory->data);
3738
3739 g_free (factory->id);
3740 g_free (factory->name);
3741
3742 gnm_solver_factory_parent_class->finalize (obj);
3743 }
3744
3745 static void
3746 gnm_solver_factory_class_init (GObjectClass *object_class)
3747 {
3748 gnm_solver_factory_parent_class =
3749 g_type_class_peek_parent (object_class);
3750
3751 object_class->finalize = gnm_solver_factory_finalize;
3752 }
3753
3754 GSF_CLASS (GnmSolverFactory, gnm_solver_factory,
3755 gnm_solver_factory_class_init, NULL, G_TYPE_OBJECT)
3756
3757
3758 static GSList *solvers;
3759
3760 /**
3761 * gnm_solver_db_get:
3762 *
3763 * Returns: (transfer none) (element-type GnmSolverFactory): list of
3764 * registered solver factories.
3765 */
3766 GSList *
3767 gnm_solver_db_get (void)
3768 {
3769 return solvers;
3770 }
3771
3772 /**
3773 * gnm_solver_factory_new:
3774 * @id: Unique identifier
3775 * @name: Translated name for UI purposes
3776 * @type: Model type created by factory
3777 * @creator: (scope notified): callback for creating a solver
3778 * @functional: (scope notified): callback for checking if factory is functional
3779 * @data: User pointer for @creator and @functional
3780 * @notify: Destroy notification for @data.
3781 *
3782 * Returns: (transfer full): a new #GnmSolverFactory
3783 */
3784 GnmSolverFactory *
3785 gnm_solver_factory_new (const char *id,
3786 const char *name,
3787 GnmSolverModelType type,
3788 GnmSolverCreator creator,
3789 GnmSolverFactoryFunctional functional,
3790 gpointer data,
3791 GDestroyNotify notify)
3792 {
3793 GnmSolverFactory *res;
3794
3795 g_return_val_if_fail (id != NULL, NULL);
3796 g_return_val_if_fail (name != NULL, NULL);
3797 g_return_val_if_fail (creator != NULL, NULL);
3798
3799 res = g_object_new (GNM_SOLVER_FACTORY_TYPE, NULL);
3800 res->id = g_strdup (id);
3801 res->name = g_strdup (name);
3802 res->type = type;
3803 res->creator = creator;
3804 res->functional = functional;
3805 res->data = data;
3806 res->notify = notify;
3807 return res;
3808 }
3809
3810 /**
3811 * gnm_solver_factory_create:
3812 * @factory: #GnmSolverFactory
3813 * @param: #GnmSolverParameters
3814 *
3815 * Returns: (transfer full): a new #GnmSolver
3816 */
3817 GnmSolver *
3818 gnm_solver_factory_create (GnmSolverFactory *factory,
3819 GnmSolverParameters *param)
3820 {
3821 g_return_val_if_fail (GNM_IS_SOLVER_FACTORY (factory), NULL);
3822 return factory->creator (factory, param, factory->data);
3823 }
3824
3825 gboolean
3826 gnm_solver_factory_functional (GnmSolverFactory *factory,
3827 WBCGtk *wbcg)
3828 {
3829 if (factory == NULL)
3830 return FALSE;
3831
3832 return (factory->functional == NULL ||
3833 factory->functional (factory, wbcg, factory->data));
3834 }
3835
3836 static int
3837 cb_compare_factories (GnmSolverFactory *a, GnmSolverFactory *b)
3838 {
3839 return go_utf8_collate_casefold (a->name, b->name);
3840 }
3841
3842 void
3843 gnm_solver_db_register (GnmSolverFactory *factory)
3844 {
3845 if (gnm_solver_debug ())
3846 g_printerr ("Registering %s\n", factory->id);
3847 g_object_ref (factory);
3848 solvers = g_slist_insert_sorted (solvers, factory,
3849 (GCompareFunc)cb_compare_factories);
3850 }
3851
3852 void
3853 gnm_solver_db_unregister (GnmSolverFactory *factory)
3854 {
3855 if (gnm_solver_debug ())
3856 g_printerr ("Unregistering %s\n", factory->id);
3857 solvers = g_slist_remove (solvers, factory);
3858 g_object_unref (factory);
3859 }
3860
3861 /* ------------------------------------------------------------------------- */
gnumeric mirror