?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
1.12.42
[gnumeric.git] / src / tools / gnm-solver.c
blobcca96d5cb3f6b5d0cb6976a94601eb00720b0888
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 <tools/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 sol->gradient_status = 0;
932 if (sol->gradient) {
933 g_ptr_array_unref (sol->gradient);
934 sol->gradient = NULL;
935 }
936
937 sol->hessian_status = 0;
938 if (sol->hessian) {
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 static gnm_float
1271 get_cell_value (GnmCell *cell)
1272 {
1273 GnmValue const *v;
1274 gnm_cell_eval (cell);
1275 v = cell->value;
1276 return VALUE_IS_NUMBER (v) || VALUE_IS_EMPTY (v)
1277 ? value_get_as_float (v)
1278 : gnm_nan;
1279 }
1280
1281 gboolean
1282 gnm_solver_check_constraints (GnmSolver *solver)
1283 {
1284 GSList *l;
1285 GnmSolverParameters *sp = solver->params;
1286 GnmCell *target_cell;
1287
1288 if (sp->options.assume_non_negative ||
1289 sp->options.assume_discrete) {
1290 unsigned ui;
1291 gboolean bad;
1292
1293 for (ui = 0; ui < solver->input_cells->len; ui++) {
1294 GnmCell *cell = g_ptr_array_index (solver->input_cells, ui);
1295 gnm_float val = get_cell_value (cell);
1296
1297 if (!gnm_finite (val))
1298 break;
1299 if (sp->options.assume_non_negative && val < 0)
1300 break;
1301 if (sp->options.assume_discrete &&
1302 val != gnm_floor (val))
1303 break;
1304 }
1305 bad = (ui < solver->input_cells->len);
1306
1307 if (bad)
1308 return FALSE;
1309 }
1310
1311 for (l = sp->constraints; l; l = l->next) {
1312 GnmSolverConstraint *c = l->data;
1313 int i;
1314 gnm_float cl, cr;
1315 GnmCell *lhs, *rhs;
1316
1317 for (i = 0;
1318 gnm_solver_constraint_get_part (c, sp, i,
1319 &lhs, &cl,
1320 &rhs, &cr);
1321 i++) {
1322 if (lhs)
1323 cl = get_cell_value (lhs);
1324 if (rhs)
1325 cr = get_cell_value (rhs);
1326
1327 switch (c->type) {
1328 case GNM_SOLVER_INTEGER:
1329 if (cl == gnm_floor (cl))
1330 continue;
1331 return FALSE;
1332 case GNM_SOLVER_BOOLEAN:
1333 if (cl == 0 || cl == 1)
1334 continue;
1335 return FALSE;
1336 case GNM_SOLVER_LE:
1337 if (cl <= cr)
1338 continue;
1339 return FALSE;
1340 case GNM_SOLVER_GE:
1341 if (cl >= cr)
1342 continue;
1343 return FALSE;
1344 case GNM_SOLVER_EQ:
1345 if (cl == cr)
1346 continue;
1347 return FALSE;
1348 default:
1349 g_assert_not_reached ();
1350 return FALSE;
1351 }
1352 }
1353 }
1354
1355 target_cell = gnm_solver_param_get_target_cell (sp);
1356 gnm_cell_eval (target_cell);
1357 if (!target_cell || !VALUE_IS_NUMBER (target_cell->value))
1358 return FALSE;
1359
1360 return TRUE;
1361 }
1362
1363 gboolean
1364 gnm_solver_saveas (GnmSolver *solver, WorkbookControl *wbc,
1365 GOFileSaver *fs,
1366 const char *templ, char **filename,
1367 GError **err)
1368 {
1369 int fd;
1370 GsfOutput *output;
1371 FILE *file;
1372 GOIOContext *io_context;
1373 gboolean ok;
1374 WorkbookView *wbv = wb_control_view (wbc);
1375
1376 fd = g_file_open_tmp (templ, filename, err);
1377 if (fd == -1) {
1378 g_set_error (err, G_FILE_ERROR, 0,
1379 _("Failed to create file for linear program"));
1380 return FALSE;
1381 }
1382
1383 file = fdopen (fd, "wb");
1384 if (!file) {
1385 /* This shouldn't really happen. */
1386 close (fd);
1387 g_set_error (err, G_FILE_ERROR, 0,
1388 _("Failed to create linear program file"));
1389 return FALSE;
1390 }
1391
1392 /* Give the saver a way to talk to the solver. */
1393 g_object_set_data_full (G_OBJECT (fs),
1394 "solver", g_object_ref (solver),
1395 (GDestroyNotify)g_object_unref);
1396
1397 output = gsf_output_stdio_new_FILE (*filename, file, TRUE);
1398 io_context = go_io_context_new (GO_CMD_CONTEXT (wbc));
1399 workbook_view_save_to_output (wbv, fs, output, io_context);
1400 ok = !go_io_error_occurred (io_context);
1401 g_object_unref (io_context);
1402 g_object_unref (output);
1403
1404 g_object_set_data (G_OBJECT (fs), "solver", NULL);
1405
1406 if (!ok) {
1407 g_set_error (err, G_FILE_ERROR, 0,
1408 _("Failed to save linear program"));
1409 return FALSE;
1410 }
1411
1412 return TRUE;
1413 }
1414
1415 int
1416 gnm_solver_cell_index (GnmSolver *solver, GnmCell const *cell)
1417 {
1418 gpointer idx;
1419
1420 if (g_hash_table_lookup_extended (solver->index_from_cell,
1421 (gpointer)cell,
1422 NULL, &idx))
1423 return GPOINTER_TO_INT (idx);
1424 else
1425 return -1;
1426 }
1427
1428 static int
1429 cell_in_cr (GnmSolver *sol, GnmCell const *cell, gboolean follow)
1430 {
1431 int idx;
1432
1433 if (!cell)
1434 return -1;
1435
1436 idx = gnm_solver_cell_index (sol, cell);
1437 if (idx < 0 && follow) {
1438 /* If the expression is just =X42 then look at X42 instead.
1439 This is because the mps loader uses such a level of
1440 indirection. Note: we follow only one such step. */
1441 GnmCellRef const *cr = gnm_expr_top_get_cellref (cell->base.texpr);
1442 GnmCellRef cr2;
1443 GnmCell const *new_cell;
1444 GnmEvalPos ep;
1445
1446 if (!cr)
1447 return -1;
1448
1449 eval_pos_init_cell (&ep, cell);
1450 gnm_cellref_make_abs (&cr2, cr, &ep);
1451 new_cell = sheet_cell_get (eval_sheet (cr2.sheet, cell->base.sheet),
1452 cr2.col, cr2.row);
1453 return cell_in_cr (sol, new_cell, FALSE);
1454 }
1455
1456 return idx;
1457 }
1458
1459 static gboolean
1460 cell_is_constant (GnmCell *cell, gnm_float *pc)
1461 {
1462 if (!cell) {
1463 *pc = 0;
1464 return TRUE;
1465 }
1466
1467 if (cell->base.texpr)
1468 return FALSE;
1469
1470 *pc = get_cell_value (cell);
1471 return gnm_finite (*pc);
1472 }
1473
1474 #define SET_LOWER(l_) \
1475 do { \
1476 sol->min[idx] = MAX (sol->min[idx], (l_)); \
1477 } while (0)
1478
1479 #define SET_UPPER(l_) \
1480 do { \
1481 sol->max[idx] = MIN (sol->max[idx], (l_)); \
1482 } while (0)
1483
1484
1485
1486 static void
1487 gnm_solver_update_derived (GnmSolver *sol)
1488 {
1489 GnmSolverParameters *params = sol->params;
1490
1491 if (sol->input_cells) {
1492 g_ptr_array_free (sol->input_cells, TRUE);
1493 sol->input_cells = NULL;
1494 }
1495
1496 if (sol->index_from_cell) {
1497 g_hash_table_destroy (sol->index_from_cell);
1498 sol->index_from_cell = NULL;
1499 }
1500 sol->target = NULL;
1501
1502 g_free (sol->min);
1503 sol->min = NULL;
1504
1505 g_free (sol->max);
1506 sol->max = NULL;
1507
1508 g_free (sol->discrete);
1509 sol->discrete = NULL;
1510
1511 if (params) {
1512 unsigned ui, n;
1513 GSList *l;
1514
1515 sol->target = gnm_solver_param_get_target_cell (params);
1516
1517 sol->input_cells = gnm_solver_param_get_input_cells (params);
1518
1519 n = sol->input_cells->len;
1520 sol->index_from_cell = g_hash_table_new (g_direct_hash, g_direct_equal);
1521 for (ui = 0; ui < n; ui++) {
1522 GnmCell *cell = g_ptr_array_index (sol->input_cells, ui);
1523 g_hash_table_insert (sol->index_from_cell, cell, GUINT_TO_POINTER (ui));
1524 }
1525
1526 sol->min = g_new (gnm_float, n);
1527 sol->max = g_new (gnm_float, n);
1528 sol->discrete = g_new (guint8, n);
1529 for (ui = 0; ui < n; ui++) {
1530 sol->min[ui] = params->options.assume_non_negative ? 0 : gnm_ninf;
1531 sol->max[ui] = gnm_pinf;
1532 sol->discrete[ui] = params->options.assume_discrete;
1533 }
1534
1535 for (l = params->constraints; l; l = l->next) {
1536 GnmSolverConstraint *c = l->data;
1537 int i;
1538 gnm_float cl, cr;
1539 GnmCell *lhs, *rhs;
1540
1541 for (i = 0;
1542 gnm_solver_constraint_get_part (c, params, i,
1543 &lhs, &cl,
1544 &rhs, &cr);
1545 i++) {
1546 int idx = cell_in_cr (sol, lhs, TRUE);
1547
1548 if (idx < 0)
1549 continue;
1550 if (!cell_is_constant (rhs, &cr))
1551 continue;
1552
1553 switch (c->type) {
1554 case GNM_SOLVER_INTEGER:
1555 sol->discrete[idx] = TRUE;
1556 break;
1557 case GNM_SOLVER_BOOLEAN:
1558 sol->discrete[idx] = TRUE;
1559 SET_LOWER (0.0);
1560 SET_UPPER (1.0);
1561 break;
1562 case GNM_SOLVER_LE:
1563 SET_UPPER (cr);
1564 break;
1565 case GNM_SOLVER_GE:
1566 SET_LOWER (cr);
1567 break;
1568 case GNM_SOLVER_EQ:
1569 SET_LOWER (cr);
1570 SET_UPPER (cr);
1571 break;
1572
1573 default:
1574 g_assert_not_reached ();
1575 break;
1576 }
1577 }
1578 }
1579
1580 /*
1581 * If parameters are discrete, narrow the range by eliminating
1582 * the fractional part of the limits.
1583 */
1584 for (ui = 0; ui < n; ui++) {
1585 if (sol->discrete[ui]) {
1586 sol->min[ui] = gnm_ceil (sol->min[ui]);
1587 sol->max[ui] = gnm_floor (sol->max[ui]);
1588 }
1589 }
1590 }
1591 }
1592
1593 #undef SET_LOWER
1594 #undef SET_UPPER
1595
1596
1597 #define ADD_HEADER(txt_) do { \
1598 dao_set_bold (dao, 0, R, 0, R); \
1599 dao_set_cell (dao, 0, R, (txt_)); \
1600 R++; \
1601 } while (0)
1602
1603 #define AT_LIMIT(s_,l_) \
1604 (gnm_finite (l_) ? gnm_abs ((s_) - (l_)) <= (gnm_abs ((s_)) + gnm_abs ((l_))) / 1e10 : (s_) == (l_))
1605
1606 #define MARK_BAD(col_) \
1607 do { \
1608 int c = (col_); \
1609 dao_set_colors (dao, c, R, c, R, \
1610 gnm_color_new_rgb8 (255, 0, 0), \
1611 NULL); \
1612 } while (0)
1613
1614
1615 static void
1616 add_value_or_special (data_analysis_output_t *dao, int col, int row,
1617 gnm_float x)
1618 {
1619 if (gnm_finite (x))
1620 dao_set_cell_float (dao, col, row, x);
1621 else {
1622 dao_set_cell (dao, col, row, "-");
1623 dao_set_align (dao, col, row, col, row,
1624 GNM_HALIGN_CENTER, GNM_VALIGN_TOP);
1625 }
1626 }
1627
1628 static void
1629 print_vector (const char *name, const gnm_float *v, int n)
1630 {
1631 int i;
1632
1633 if (name)
1634 g_printerr ("%s:\n", name);
1635 for (i = 0; i < n; i++)
1636 g_printerr ("%15.8" GNM_FORMAT_f " ", v[i]);
1637 g_printerr ("\n");
1638 }
1639
1640 static void
1641 gnm_solver_create_program_report (GnmSolver *solver, const char *name)
1642 {
1643 GnmSolverParameters *params = solver->params;
1644 int R = 0;
1645 data_analysis_output_t *dao;
1646 GSList *l;
1647
1648 dao = dao_init_new_sheet (NULL);
1649 dao->sheet = params->sheet;
1650 dao_prepare_output (NULL, dao, name);
1651
1652 /* ---------------------------------------- */
1653
1654 {
1655 char *tmp;
1656
1657 ADD_HEADER (_("Target"));
1658 dao_set_cell (dao, 1, R, _("Cell"));
1659 dao_set_cell (dao, 2, R, _("Value"));
1660 dao_set_cell (dao, 3, R, _("Type"));
1661 dao_set_cell (dao, 4, R, _("Status"));
1662 R++;
1663
1664 tmp = gnm_solver_cell_name
1665 (gnm_solver_param_get_target_cell (params),
1666 params->sheet);
1667 dao_set_cell (dao, 1, R, tmp);
1668 g_free (tmp);
1669
1670 dao_set_cell_float (dao, 2, R, solver->result->value);
1671
1672 switch (params->problem_type) {
1673 case GNM_SOLVER_MINIMIZE:
1674 dao_set_cell (dao, 3, R, _("Minimize"));
1675 break;
1676 case GNM_SOLVER_MAXIMIZE:
1677 dao_set_cell (dao, 3, R, _("Maximize"));
1678 break;
1679 }
1680
1681 switch (solver->result->quality) {
1682 default:
1683 break;
1684 case GNM_SOLVER_RESULT_FEASIBLE:
1685 dao_set_cell (dao, 4, R, _("Feasible"));
1686 break;
1687 case GNM_SOLVER_RESULT_OPTIMAL:
1688 dao_set_cell (dao, 4, R, _("Optimal"));
1689 break;
1690 }
1691
1692 R++;
1693 R++;
1694 }
1695
1696 /* ---------------------------------------- */
1697
1698 if (solver->input_cells->len > 0) {
1699 unsigned ui;
1700
1701 ADD_HEADER (_("Variables"));
1702
1703 dao_set_cell (dao, 1, R, _("Cell"));
1704 dao_set_cell (dao, 2, R, _("Value"));
1705 dao_set_cell (dao, 3, R, _("Lower"));
1706 dao_set_cell (dao, 4, R, _("Upper"));
1707 dao_set_cell (dao, 5, R, _("Slack"));
1708 R++;
1709
1710 for (ui = 0; ui < solver->input_cells->len; ui++) {
1711 GnmCell *cell = g_ptr_array_index (solver->input_cells, ui);
1712 gnm_float L = solver->min[ui];
1713 gnm_float H = solver->max[ui];
1714 gnm_float s = solver->result->solution[ui];
1715 gnm_float slack = MIN (s - L, H - s);
1716
1717 char *cname = gnm_solver_cell_name (cell, params->sheet);
1718 dao_set_cell (dao, 1, R, cname);
1719 g_free (cname);
1720 dao_set_cell_value (dao, 2, R, value_new_float (s));
1721 add_value_or_special (dao, 3, R, L);
1722 add_value_or_special (dao, 4, R, H);
1723
1724 add_value_or_special (dao, 5, R, slack);
1725 if (slack < 0)
1726 MARK_BAD (5);
1727
1728 if (AT_LIMIT (s, L) || AT_LIMIT (s, H))
1729 dao_set_cell (dao, 6, R, _("At limit"));
1730
1731 if (s < L || s > H) {
1732 dao_set_cell (dao, 7, R, _("Outside bounds"));
1733 MARK_BAD (7);
1734 }
1735
1736 R++;
1737 }
1738
1739 R++;
1740 }
1741
1742 /* ---------------------------------------- */
1743
1744 ADD_HEADER (_("Constraints"));
1745
1746 if (params->constraints) {
1747 dao_set_cell (dao, 1, R, _("Condition"));
1748 dao_set_cell (dao, 2, R, _("Value"));
1749 dao_set_cell (dao, 3, R, _("Limit"));
1750 dao_set_cell (dao, 4, R, _("Slack"));
1751 } else {
1752 dao_set_cell (dao, 1, R, _("No constraints"));
1753 }
1754 R++;
1755
1756 for (l = params->constraints; l; l = l->next) {
1757 GnmSolverConstraint *c = l->data;
1758 int i;
1759 gnm_float cl, cr;
1760 GnmCell *lhs, *rhs;
1761
1762 for (i = 0;
1763 gnm_solver_constraint_get_part (c, params, i,
1764 &lhs, &cl,
1765 &rhs, &cr);
1766 i++) {
1767 gnm_float slack = 0;
1768 char *ctxt = gnm_solver_constraint_part_as_str (c, i, params);
1769 dao_set_cell (dao, 1, R, ctxt);
1770 g_free (ctxt);
1771
1772 if (lhs)
1773 cl = get_cell_value (lhs);
1774 if (rhs)
1775 cr = get_cell_value (rhs);
1776
1777 switch (c->type) {
1778 case GNM_SOLVER_INTEGER: {
1779 gnm_float c = gnm_fake_round (cl);
1780 slack = 0 - gnm_abs (c - cl);
1781 break;
1782 }
1783 case GNM_SOLVER_BOOLEAN: {
1784 gnm_float c = (cl > 0.5 ? 1 : 0);
1785 slack = 0 - gnm_abs (c - cl);
1786 break;
1787 }
1788 case GNM_SOLVER_LE:
1789 slack = cr - cl;
1790 break;
1791 case GNM_SOLVER_GE:
1792 slack = cl - cr;
1793 break;
1794 case GNM_SOLVER_EQ:
1795 slack = 0 - gnm_abs (cl - cr);
1796 break;
1797 default:
1798 g_assert_not_reached ();
1799 }
1800
1801 add_value_or_special (dao, 2, R, cl);
1802 if (rhs)
1803 add_value_or_special (dao, 3, R, cr);
1804
1805 add_value_or_special (dao, 4, R, slack);
1806 if (slack < 0)
1807 MARK_BAD (4);
1808
1809 R++;
1810 }
1811 }
1812
1813 /* ---------------------------------------- */
1814
1815 dao_autofit_columns (dao);
1816 dao_redraw_respan (dao);
1817
1818 dao_free (dao);
1819 }
1820
1821 static void
1822 gnm_solver_create_sensitivity_report (GnmSolver *solver, const char *name)
1823 {
1824 GnmSolverParameters *params = solver->params;
1825 GnmSolverSensitivity *sols = solver->sensitivity;
1826 int R = 0;
1827 data_analysis_output_t *dao;
1828 GSList *l;
1829
1830 if (!sols)
1831 return;
1832
1833 dao = dao_init_new_sheet (NULL);
1834 dao->sheet = params->sheet;
1835 dao_prepare_output (NULL, dao, name);
1836
1837 /* ---------------------------------------- */
1838
1839 if (solver->input_cells->len > 0) {
1840 unsigned ui;
1841
1842 ADD_HEADER (_("Variables"));
1843
1844 dao_set_cell (dao, 1, R, _("Cell"));
1845 dao_set_cell (dao, 2, R, _("Final\nValue"));
1846 dao_set_cell (dao, 3, R, _("Reduced\nCost"));
1847 dao_set_cell (dao, 4, R, _("Lower\nLimit"));
1848 dao_set_cell (dao, 5, R, _("Upper\nLimit"));
1849 dao_set_align (dao, 1, R, 5, R, GNM_HALIGN_CENTER, GNM_VALIGN_BOTTOM);
1850 dao_autofit_these_rows (dao, R, R);
1851 R++;
1852
1853 for (ui = 0; ui < solver->input_cells->len; ui++) {
1854 GnmCell *cell = g_ptr_array_index (solver->input_cells, ui);
1855 gnm_float L = sols->vars[ui].low;
1856 gnm_float H = sols->vars[ui].high;
1857 gnm_float red = sols->vars[ui].reduced_cost;
1858 gnm_float s = solver->result->solution[ui];
1859
1860 char *cname = gnm_solver_cell_name (cell, params->sheet);
1861 dao_set_cell (dao, 1, R, cname);
1862 g_free (cname);
1863 dao_set_cell_value (dao, 2, R, value_new_float (s));
1864 add_value_or_special (dao, 3, R, red);
1865 add_value_or_special (dao, 4, R, L);
1866 add_value_or_special (dao, 5, R, H);
1867
1868 R++;
1869 }
1870
1871 R++;
1872 }
1873
1874 /* ---------------------------------------- */
1875
1876 ADD_HEADER (_("Constraints"));
1877
1878 if (params->constraints) {
1879 dao_set_cell (dao, 1, R, _("Constraint"));
1880 dao_set_cell (dao, 2, R, _("Shadow\nPrice"));
1881 dao_set_cell (dao, 3, R, _("Constraint\nLHS"));
1882 dao_set_cell (dao, 4, R, _("Constraint\nRHS"));
1883 dao_set_cell (dao, 5, R, _("Lower\nLimit"));
1884 dao_set_cell (dao, 6, R, _("Upper\nLimit"));
1885 dao_set_align (dao, 1, R, 6, R, GNM_HALIGN_CENTER, GNM_VALIGN_BOTTOM);
1886 dao_autofit_these_rows (dao, R, R);
1887 } else {
1888 dao_set_cell (dao, 1, R, _("No constraints"));
1889 }
1890 R++;
1891
1892 for (l = params->constraints; l; l = l->next) {
1893 GnmSolverConstraint *c = l->data;
1894 int i, cidx = 0;
1895 gnm_float cl, cr;
1896 GnmCell *lhs, *rhs;
1897
1898 for (i = 0;
1899 gnm_solver_constraint_get_part (c, params, i,
1900 &lhs, &cl,
1901 &rhs, &cr);
1902 i++, cidx++) {
1903 char *ctxt;
1904
1905 switch (c->type) {
1906 case GNM_SOLVER_INTEGER:
1907 case GNM_SOLVER_BOOLEAN:
1908 continue;
1909 default:
1910 ; // Nothing
1911 }
1912
1913 ctxt = gnm_solver_constraint_part_as_str (c, i, params);
1914 dao_set_cell (dao, 1, R, ctxt);
1915 g_free (ctxt);
1916
1917 if (lhs)
1918 cl = get_cell_value (lhs);
1919 if (rhs)
1920 cr = get_cell_value (rhs);
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 gnm_float y = get_cell_value (solver->target);
1973 return solver->flip_sign ? 0 - y : y;
1974 }
1975
1976 void
1977 gnm_solver_set_var (GnmSolver *sol, int i, gnm_float x)
1978 {
1979 GnmCell *cell = g_ptr_array_index (sol->input_cells, i);
1980
1981 if (cell->value &&
1982 VALUE_IS_FLOAT (cell->value) &&
1983 value_get_as_float (cell->value) == x)
1984 return;
1985
1986 gnm_cell_set_value (cell, value_new_float (x));
1987 cell_queue_recalc (cell);
1988 }
1989
1990 void
1991 gnm_solver_set_vars (GnmSolver *sol, gnm_float const *xs)
1992 {
1993 const int n = sol->input_cells->len;
1994 int i;
1995
1996 for (i = 0; i < n; i++)
1997 gnm_solver_set_var (sol, i, xs[i]);
1998 }
1999
2000 /**
2001 * gnm_solver_save_vars:
2002 * @sol: #GnmSolver
2003 *
2004 * Returns: (transfer full) (element-type GnmValue):
2005 */
2006 GPtrArray *
2007 gnm_solver_save_vars (GnmSolver *sol)
2008 {
2009 GPtrArray *vals = g_ptr_array_new ();
2010 unsigned ui;
2011
2012 for (ui = 0; ui < sol->input_cells->len; ui++) {
2013 GnmCell *cell = g_ptr_array_index (sol->input_cells, ui);
2014 g_ptr_array_add (vals, value_dup (cell->value));
2015 }
2016
2017 return vals;
2018 }
2019
2020 /**
2021 * gnm_solver_restore_vars:
2022 * @sol: #GnmSolver
2023 * @vals: (transfer full) (element-type GnmValue): values to restore
2024 */
2025 void
2026 gnm_solver_restore_vars (GnmSolver *sol, GPtrArray *vals)
2027 {
2028 unsigned ui;
2029
2030 for (ui = 0; ui < sol->input_cells->len; ui++) {
2031 GnmCell *cell = g_ptr_array_index (sol->input_cells, ui);
2032 gnm_cell_set_value (cell, g_ptr_array_index (vals, ui));
2033 cell_queue_recalc (cell);
2034 }
2035
2036 g_ptr_array_free (vals, TRUE);
2037 }
2038
2039 /**
2040 * gnm_solver_has_analytic_gradient:
2041 * @sol: the solver
2042 *
2043 * Returns: %TRUE if the gradient can be computed analytically.
2044 */
2045 gboolean
2046 gnm_solver_has_analytic_gradient (GnmSolver *sol)
2047 {
2048 const int n = sol->input_cells->len;
2049
2050 if (sol->gradient_status == 0) {
2051 int i;
2052
2053 sol->gradient_status++;
2054
2055 sol->gradient = g_ptr_array_new_with_free_func ((GDestroyNotify)gnm_expr_top_unref);
2056 for (i = 0; i < n; i++) {
2057 GnmCell *cell = g_ptr_array_index (sol->input_cells, i);
2058 GnmExprTop const *te =
2059 gnm_expr_cell_deriv (sol->target, cell);
2060 if (te)
2061 g_ptr_array_add (sol->gradient, (gpointer)te);
2062 else {
2063 if (gnm_solver_debug ())
2064 g_printerr ("Unable to compute analytic gradient\n");
2065 g_ptr_array_unref (sol->gradient);
2066 sol->gradient = NULL;
2067 sol->gradient_status++;
2068 break;
2069 }
2070 }
2071 }
2072
2073 return sol->gradient_status == 1;
2074 }
2075
2076 static gnm_float *
2077 gnm_solver_compute_gradient_analytically (GnmSolver *sol, gnm_float const *xs)
2078 {
2079 const int n = sol->input_cells->len;
2080 int i;
2081 gnm_float *g = g_new (gnm_float, n);
2082 GnmEvalPos ep;
2083
2084 eval_pos_init_cell (&ep, sol->target);
2085 for (i = 0; i < n; i++) {
2086 GnmExprTop const *te = g_ptr_array_index (sol->gradient, i);
2087 GnmValue *v = gnm_expr_top_eval
2088 (te, &ep, GNM_EXPR_EVAL_SCALAR_NON_EMPTY);
2089 g[i] = VALUE_IS_NUMBER (v) ? value_get_as_float (v) : gnm_nan;
2090 if (sol->flip_sign)
2091 g[i] = 0 - g[i];
2092 value_release (v);
2093 }
2094
2095 if (gnm_solver_debug ())
2096 print_vector ("Analytic gradient", g, n);
2097
2098 return g;
2099 }
2100
2101 /**
2102 * gnm_solver_compute_gradient: (skip)
2103 * @sol: Solver
2104 * @xs: Point to compute gradient at
2105 *
2106 * Returns: xxx(transfer full): A vector containing the gradient. This
2107 * function takes the flip-sign property into account. This will use
2108 * analytic gradient, if possible, and a numerical approximation otherwise.
2109 */
2110 gnm_float *
2111 gnm_solver_compute_gradient (GnmSolver *sol, gnm_float const *xs)
2112 {
2113 gnm_float *g;
2114 gnm_float y0;
2115 const int n = sol->input_cells->len;
2116 int i;
2117 const int order = sol->params->options.gradient_order;
2118
2119 gnm_solver_set_vars (sol, xs);
2120 y0 = gnm_solver_get_target_value (sol);
2121
2122 if (gnm_solver_has_analytic_gradient (sol))
2123 return gnm_solver_compute_gradient_analytically (sol, xs);
2124
2125 g = g_new (gnm_float, n);
2126 for (i = 0; i < n; i++) {
2127 gnm_float x0 = xs[i];
2128 gnm_float dx, dy;
2129 int j;
2130
2131 /*
2132 * This computes the least-squares fit of an affine function
2133 * based on 2*order equidistant points symmetrically around
2134 * the value we compute the derivative for.
2135 *
2136 * We use an even number of ULPs for the step size. This
2137 * ensures that the x values are computed without rounding
2138 * error except, potentially, a single step that crosses an
2139 * integer power of 2.
2140 */
2141 dx = 16 * (go_add_epsilon (x0) - x0);
2142 dy = 0;
2143 for (j = -order; j <= order; j++) {
2144 gnm_float y;
2145
2146 if (j == 0)
2147 continue;
2148
2149 gnm_solver_set_var (sol, i, x0 + j * dx);
2150 y = gnm_solver_get_target_value (sol);
2151 dy += j * (y - y0);
2152 }
2153 dy /= 2 * (order * (2 * order * order + 3 * order + 1) / 6);
2154 g[i] = dy / dx;
2155
2156 gnm_solver_set_var (sol, i, x0);
2157 }
2158
2159 if (gnm_solver_debug ())
2160 print_vector ("Numerical gradient", g, n);
2161
2162 return g;
2163 }
2164
2165 /**
2166 * gnm_solver_has_analytic_hessian:
2167 * @sol: the solver
2168 *
2169 * Returns: %TRUE if the Hessian can be computed analytically.
2170 */
2171 gboolean
2172 gnm_solver_has_analytic_hessian (GnmSolver *sol)
2173 {
2174 const int n = sol->input_cells->len;
2175 int i, j;
2176 GnmEvalPos ep;
2177 GnmExprDeriv *info;
2178
2179 if (!gnm_solver_has_analytic_gradient (sol))
2180 sol->hessian_status = sol->gradient_status;
2181
2182 if (sol->hessian_status)
2183 return sol->hessian_status == 1;
2184
2185 sol->hessian_status++;
2186 sol->hessian = g_ptr_array_new_with_free_func ((GDestroyNotify)gnm_expr_top_unref);
2187
2188 eval_pos_init_cell (&ep, sol->target);
2189 info = gnm_expr_deriv_info_new ();
2190 for (i = 0; i < n && sol->hessian_status == 1; i++) {
2191 GnmExprTop const *gi = g_ptr_array_index (sol->gradient, i);
2192 for (j = i; j < n; j++) {
2193 GnmCell *cell;
2194 GnmExprTop const *te;
2195 GnmEvalPos var;
2196
2197 cell = g_ptr_array_index (sol->input_cells, j);
2198 eval_pos_init_cell (&var, cell);
2199 gnm_expr_deriv_info_set_var (info, &var);
2200 te = gnm_expr_top_deriv (gi, &ep, info);
2201
2202 if (te)
2203 g_ptr_array_add (sol->hessian, (gpointer)te);
2204 else {
2205 if (gnm_solver_debug ())
2206 g_printerr ("Unable to compute analytic hessian\n");
2207 sol->hessian_status++;
2208 break;
2209 }
2210 }
2211 }
2212
2213 gnm_expr_deriv_info_unref (info);
2214
2215 return sol->hessian_status == 1;
2216 }
2217
2218 /**
2219 * gnm_solver_compute_hessian:
2220 * @sol: Solver
2221 * @xs: Point to compute Hessian at
2222 *
2223 * Returns: (transfer full): A matrix containing the Hessian. This
2224 * function takes the flip-sign property into account.
2225 */
2226 GnmMatrix *
2227 gnm_solver_compute_hessian (GnmSolver *sol, gnm_float const *xs)
2228 {
2229 int i, j, k;
2230 GnmMatrix *H;
2231 GnmEvalPos ep;
2232 int const n = sol->input_cells->len;
2233
2234 if (!gnm_solver_has_analytic_hessian (sol))
2235 return NULL;
2236
2237 gnm_solver_set_vars (sol, xs);
2238
2239 H = gnm_matrix_new (n, n);
2240 eval_pos_init_cell (&ep, sol->target);
2241 for (i = k = 0; i < n; i++) {
2242 for (j = i; j < n; j++, k++) {
2243 GnmExprTop const *te =
2244 g_ptr_array_index (sol->hessian, k);
2245 GnmValue *v = gnm_expr_top_eval
2246 (te, &ep, GNM_EXPR_EVAL_SCALAR_NON_EMPTY);
2247 gnm_float x = VALUE_IS_NUMBER (v)
2248 ? value_get_as_float (v)
2249 : gnm_nan;
2250 if (sol->flip_sign)
2251 x = 0 - x;
2252 value_release (v);
2253
2254 H->data[i][j] = x;
2255 H->data[j][i] = x;
2256 }
2257 }
2258
2259 return H;
2260 }
2261
2262
2263 static gnm_float
2264 try_step (GnmSolver *sol, gnm_float const *x0, gnm_float const *dir, gnm_float step)
2265 {
2266 int const n = sol->input_cells->len;
2267 gnm_float *x = g_new (gnm_float, n);
2268 int i;
2269 gnm_float y;
2270
2271 for (i = 0; i < n; i++)
2272 x[i] = x0[i] + step * dir[i];
2273 gnm_solver_set_vars (sol, x);
2274 y = gnm_solver_get_target_value (sol);
2275 g_free (x);
2276 return y;
2277 }
2278
2279 /**
2280 * gnm_solver_line_search:
2281 * @sol: Solver
2282 * @x0: Starting point
2283 * @dir: direction
2284 * @try_reverse: whether to try reverse direction at all
2285 * @step: initial step size
2286 * @max_step: largest allowed step
2287 * @eps: tolerance for optimal step
2288 * @py: (out): location to store resulting objective function value
2289 *
2290 * Returns: optimal step size.
2291 */
2292 gnm_float
2293 gnm_solver_line_search (GnmSolver *sol,
2294 gnm_float const *x0, gnm_float const *dir,
2295 gboolean try_reverse,
2296 gnm_float step, gnm_float max_step, gnm_float eps,
2297 gnm_float *py)
2298 {
2299 /*
2300 * 0: Initial
2301 * 1: Found improvement, but not far side of it
2302 * 2: Have (s0,s1,s2) with s1 lowest
2303 */
2304 int phase = 0;
2305 gnm_float s, s0, s1, s2;
2306 gnm_float y, y0, y1, y2;
2307 gnm_float const phi = (gnm_sqrt (5) + 1) / 2;
2308 gboolean rbig;
2309 gboolean debug = gnm_solver_debug ();
2310
2311 g_return_val_if_fail (eps >= 0, gnm_nan);
2312 g_return_val_if_fail (step > 0, gnm_nan);
2313 g_return_val_if_fail (max_step >= step, gnm_nan);
2314
2315 if (debug) {
2316 g_printerr ("LS: step=%" GNM_FORMAT_g ", max=%" GNM_FORMAT_g ", eps=%" GNM_FORMAT_g "\n", step, max_step, eps);
2317 print_vector (NULL, dir, sol->input_cells->len);
2318 }
2319
2320 gnm_solver_set_vars (sol, x0);
2321 y0 = gnm_solver_get_target_value (sol);
2322 s0 = 0;
2323
2324 s = step;
2325 while (phase == 0) {
2326 gboolean flat = TRUE;
2327
2328 y = try_step (sol, x0, dir, s);
2329 if (0 && debug)
2330 g_printerr ("LS0: s:%.6" GNM_FORMAT_g " y:%.6" GNM_FORMAT_g "\n",
2331 s, y);
2332 if (y < y0 && gnm_solver_check_constraints (sol)) {
2333 y1 = y;
2334 s1 = s;
2335 phase = 1;
2336 break;
2337 } else if (y != y0)
2338 flat = FALSE;
2339
2340 if (try_reverse) {
2341 y = try_step (sol, x0, dir, -s);
2342 if (0 && debug)
2343 g_printerr ("LS0: s:%.6" GNM_FORMAT_g " y:%.6" GNM_FORMAT_g "\n",
2344 -s, y);
2345 if (y < y0 && gnm_solver_check_constraints (sol)) {
2346 y1 = y;
2347 s1 = -s;
2348 phase = 1;
2349 break;
2350 } else if (y != y0)
2351 flat = FALSE;
2352 }
2353
2354 s /= 32;
2355
2356 if (s <= 0 || flat)
2357 return gnm_nan;
2358 }
2359
2360 while (phase == 1) {
2361 s = s1 * (phi + 1);
2362
2363 if (gnm_abs (s) >= max_step)
2364 goto bail;
2365
2366 y = try_step (sol, x0, dir, s);
2367 if (!gnm_finite (y) || !gnm_solver_check_constraints (sol))
2368 goto bail;
2369
2370 if (y < y1) {
2371 y1 = y;
2372 s1 = s;
2373 continue;
2374 }
2375
2376 y2 = y;
2377 s2 = s;
2378 phase = 2;
2379 }
2380
2381 /*
2382 * Phase 2: we have three steps, s0/s1/s2, in order (descending or ascending) such
2383 * that
2384 * 1. y1<=y0 (equality unlikely)
2385 * 2. y1<=y2 (equality unlikely)
2386 * 3a. (s2-s1)=phi*(s1-s0) or
2387 * 3b. (s2-s1)*phi=(s1-s0)
2388 */
2389 rbig = TRUE; /* Initially 3a holds. */
2390 while (phase == 2) {
2391 if (0 && debug) {
2392 gnm_float s01 = s1 - s0;
2393 gnm_float s12 = s2 - s1;
2394 g_printerr ("LS2: s0:%.6" GNM_FORMAT_g " s01=%.6" GNM_FORMAT_g " s12=%.6" GNM_FORMAT_g
2395 " r=%" GNM_FORMAT_g
2396 " y:%.10" GNM_FORMAT_g "/%.10" GNM_FORMAT_g "/%.10" GNM_FORMAT_g "\n",
2397 s0, s01, s12, s12 / s01, y0, y1, y2);
2398 }
2399
2400 s = rbig ? s1 + (s1 - s0) * (phi - 1) : s1 - (s2 - s1) * (phi - 1);
2401 if (s <= s0 || s >= s2 || gnm_abs (s - s1) <= eps)
2402 break;
2403
2404 y = try_step (sol, x0, dir, s);
2405 if (!gnm_finite (y) || !gnm_solver_check_constraints (sol))
2406 goto bail;
2407
2408 if (y < y1) {
2409 if (rbig) {
2410 y0 = y1;
2411 s0 = s1;
2412 } else {
2413 y2 = y1;
2414 s2 = s1;
2415 }
2416 y1 = y;
2417 s1 = s;
2418 } else {
2419 if (rbig) {
2420 y2 = y;
2421 s2 = s;
2422 } else {
2423 y0 = y;
2424 s0 = s;
2425 }
2426 rbig = !rbig;
2427
2428 if (y0 == y1 && y1 == y2)
2429 break;
2430 }
2431 }
2432
2433 bail:
2434 if (debug)
2435 g_printerr ("LS: step %.6" GNM_FORMAT_g "\n", s1);
2436
2437 *py = y1;
2438 return s1;
2439 }
2440
2441 /**
2442 * gnm_solver_pick_lp_coords:
2443 * @sol: Solver
2444 * @px1: (out): first coordinate value
2445 * @px2: (out): second coordinate value
2446 *
2447 * Pick two good values for each coordinate. We prefer 0 and 1
2448 * when they are valid.
2449 */
2450 void
2451 gnm_solver_pick_lp_coords (GnmSolver *sol,
2452 gnm_float **px1, gnm_float **px2)
2453 {
2454 const unsigned n = sol->input_cells->len;
2455 gnm_float *x1 = *px1 = g_new (gnm_float, n);
2456 gnm_float *x2 = *px2 = g_new (gnm_float, n);
2457 unsigned ui;
2458
2459 for (ui = 0; ui < n; ui++) {
2460 const gnm_float L = sol->min[ui], H = sol->max[ui];
2461
2462 if (L == H) {
2463 x1[ui] = x2[ui] = L;
2464 } else if (sol->discrete[ui] && H - L == 1) {
2465 x1[ui] = L;
2466 x2[ui] = H;
2467 } else {
2468 if (L <= 0 && H >= 0)
2469 x1[ui] = 0;
2470 else if (gnm_finite (L))
2471 x1[ui] = L;
2472 else
2473 x1[ui] = H;
2474
2475 if (x1[ui] + 1 <= H)
2476 x2[ui] = x1[ui] + 1;
2477 else if (x1[ui] - 1 >= H)
2478 x2[ui] = x1[ui] - 1;
2479 else if (x1[ui] != H)
2480 x2[ui] = (x1[ui] + H) / 2;
2481 else
2482 x2[ui] = (x1[ui] + L) / 2;
2483 }
2484 }
2485 }
2486
2487 /**
2488 * gnm_solver_get_lp_coeffs: (skip)
2489 * @sol: Solver
2490 * @ycell: Cell for which to compute coefficients
2491 * @x1: first coordinate value
2492 * @x2: second coordinate value
2493 * @err: error location
2494 *
2495 * Returns: xxx(transfer full) (nullable): coordinates, or %NULL in case of error.
2496 * Note: this function is not affected by the flip-sign property, even
2497 * if @ycell happens to coindice with the solver target cell.
2498 */
2499 gnm_float *
2500 gnm_solver_get_lp_coeffs (GnmSolver *sol, GnmCell *ycell,
2501 gnm_float const *x1, gnm_float const *x2,
2502 GError **err)
2503 {
2504 const unsigned n = sol->input_cells->len;
2505 unsigned ui;
2506 gnm_float *res = g_new (gnm_float, n);
2507 gnm_float y0;
2508
2509 gnm_solver_set_vars (sol, x1);
2510 y0 = get_cell_value (ycell);
2511 if (!gnm_finite (y0))
2512 goto fail_calc;
2513
2514 for (ui = 0; ui < n; ui++) {
2515 gnm_float dx = x2[ui] - x1[ui], dy, y1;
2516
2517 if (dx <= 0) {
2518 res[ui] = 0;
2519 continue;
2520 }
2521
2522 gnm_solver_set_var (sol, ui, x2[ui]);
2523 y1 = get_cell_value (ycell);
2524
2525 dy = y1 - y0;
2526 res[ui] = dy / dx;
2527 if (!gnm_finite (res[ui]))
2528 goto fail_calc;
2529
2530 if (!sol->discrete[ui] || dx != 1) {
2531 gnm_float x01, y01, e, emax;
2532
2533 x01 = (x1[ui] + x2[ui]) / 2;
2534 if (sol->discrete[ui]) x01 = gnm_floor (x01);
2535 gnm_solver_set_var (sol, ui, x01);
2536 y01 = get_cell_value (ycell);
2537 if (!gnm_finite (y01))
2538 goto fail_calc;
2539
2540 emax = dy == 0 ? 1e-10 : gnm_abs (dy) / 1e-10;
2541 e = dy - 2 * (y01 - y0);
2542 if (gnm_abs (e) > emax)
2543 goto fail_linear;
2544 }
2545
2546 gnm_solver_set_var (sol, ui, x1[ui]);
2547 }
2548
2549 return res;
2550
2551 fail_calc:
2552 g_set_error (err,
2553 go_error_invalid (),
2554 0,
2555 _("Target cell did not evaluate to a number."));
2556 g_free (res);
2557 return NULL;
2558
2559 fail_linear:
2560 g_set_error (err,
2561 go_error_invalid (),
2562 0,
2563 _("Target cell does not appear to depend linearly on input cells."));
2564 g_free (res);
2565 return NULL;
2566 }
2567
2568
2569 static void
2570 gnm_solver_class_init (GObjectClass *object_class)
2571 {
2572 gnm_solver_parent_class = g_type_class_peek_parent (object_class);
2573
2574 object_class->dispose = gnm_solver_dispose;
2575 object_class->set_property = gnm_solver_set_property;
2576 object_class->get_property = gnm_solver_get_property;
2577
2578 g_object_class_install_property (object_class, SOL_PROP_STATUS,
2579 g_param_spec_enum ("status",
2580 P_("status"),
2581 P_("The solver's current status"),
2582 GNM_SOLVER_STATUS_TYPE,
2583 GNM_SOLVER_STATUS_READY,
2584 GSF_PARAM_STATIC |
2585 G_PARAM_READWRITE));
2586
2587 g_object_class_install_property (object_class, SOL_PROP_REASON,
2588 g_param_spec_string ("reason",
2589 P_("reason"),
2590 P_("The reason behind the solver's status"),
2591 NULL,
2592 GSF_PARAM_STATIC |
2593 G_PARAM_READWRITE));
2594
2595 g_object_class_install_property (object_class, SOL_PROP_PARAMS,
2596 g_param_spec_object ("params",
2597 P_("Parameters"),
2598 P_("Solver parameters"),
2599 GNM_SOLVER_PARAMETERS_TYPE,
2600 GSF_PARAM_STATIC |
2601 G_PARAM_CONSTRUCT_ONLY |
2602 G_PARAM_READWRITE));
2603
2604 g_object_class_install_property (object_class, SOL_PROP_RESULT,
2605 g_param_spec_object ("result",
2606 P_("Result"),
2607 P_("Current best feasible result"),
2608 GNM_SOLVER_RESULT_TYPE,
2609 GSF_PARAM_STATIC |
2610 G_PARAM_READWRITE));
2611
2612 g_object_class_install_property (object_class, SOL_PROP_SENSITIVITY,
2613 g_param_spec_object ("sensitivity",
2614 P_("Sensitivity"),
2615 P_("Sensitivity results"),
2616 GNM_SOLVER_SENSITIVITY_TYPE,
2617 GSF_PARAM_STATIC |
2618 G_PARAM_READWRITE));
2619
2620 g_object_class_install_property (object_class, SOL_PROP_STARTTIME,
2621 g_param_spec_double ("starttime",
2622 P_("Start Time"),
2623 P_("Time the solver was started"),
2624 -1, 1e10, -1,
2625 GSF_PARAM_STATIC |
2626 G_PARAM_READWRITE));
2627
2628 g_object_class_install_property (object_class, SOL_PROP_ENDTIME,
2629 g_param_spec_double ("endtime",
2630 P_("End Time"),
2631 P_("Time the solver finished"),
2632 -1, 1e10, -1,
2633 GSF_PARAM_STATIC |
2634 G_PARAM_READWRITE));
2635
2636 g_object_class_install_property (object_class, SOL_PROP_FLIP_SIGN,
2637 g_param_spec_boolean ("flip-sign",
2638 P_("Flip Sign"),
2639 P_("Flip sign of target value"),
2640 FALSE,
2641 GSF_PARAM_STATIC | G_PARAM_READWRITE));
2642
2643 solver_signals[SOL_SIG_PREPARE] =
2644 g_signal_new ("prepare",
2645 G_OBJECT_CLASS_TYPE (object_class),
2646 G_SIGNAL_RUN_LAST,
2647 G_STRUCT_OFFSET (GnmSolverClass, prepare),
2648 NULL, NULL,
2649 gnm__BOOLEAN__OBJECT_POINTER,
2650 G_TYPE_BOOLEAN, 2,
2651 G_TYPE_OBJECT,
2652 G_TYPE_POINTER);
2653
2654 solver_signals[SOL_SIG_START] =
2655 g_signal_new ("start",
2656 G_OBJECT_CLASS_TYPE (object_class),
2657 G_SIGNAL_RUN_LAST,
2658 G_STRUCT_OFFSET (GnmSolverClass, start),
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_STOP] =
2666 g_signal_new ("stop",
2667 G_OBJECT_CLASS_TYPE (object_class),
2668 G_SIGNAL_RUN_LAST,
2669 G_STRUCT_OFFSET (GnmSolverClass, stop),
2670 NULL, NULL,
2671 gnm__BOOLEAN__POINTER,
2672 G_TYPE_BOOLEAN, 1,
2673 G_TYPE_POINTER);
2674 }
2675
2676 GSF_CLASS (GnmSolver, gnm_solver,
2677 &gnm_solver_class_init, NULL, G_TYPE_OBJECT)
2678
2679 /* ------------------------------------------------------------------------- */
2680
2681 static GObjectClass *gnm_solver_result_parent_class;
2682
2683 static void
2684 gnm_solver_result_finalize (GObject *obj)
2685 {
2686 GnmSolverResult *r = GNM_SOLVER_RESULT (obj);
2687 g_free (r->solution);
2688 gnm_solver_result_parent_class->finalize (obj);
2689 }
2690
2691 static void
2692 gnm_solver_result_class_init (GObjectClass *object_class)
2693 {
2694 gnm_solver_result_parent_class =
2695 g_type_class_peek_parent (object_class);
2696
2697 object_class->finalize = gnm_solver_result_finalize;
2698 }
2699
2700 GSF_CLASS (GnmSolverResult, gnm_solver_result,
2701 &gnm_solver_result_class_init, NULL, G_TYPE_OBJECT)
2702
2703 /* ------------------------------------------------------------------------- */
2704
2705 static GObjectClass *gnm_solver_sensitivity_parent_class;
2706
2707 enum {
2708 SOLS_PROP_0,
2709 SOLS_PROP_SOLVER
2710 };
2711
2712 static void
2713 gnm_solver_sensitivity_constructed (GObject *obj)
2714 {
2715 GnmSolverSensitivity *sols = GNM_SOLVER_SENSITIVITY (obj);
2716 GnmSolver *sol = sols->solver;
2717 GnmSolverParameters *sp = sol->params;
2718 const int n = sol->input_cells->len;
2719 int i, cn;
2720 GSList *l;
2721
2722 /* Chain to parent first */
2723 gnm_solver_sensitivity_parent_class->constructed (obj);
2724
2725 sols->vars = g_new (struct GnmSolverSensitivityVars_, n);
2726 for (i = 0; i < n; i++) {
2727 sols->vars[i].low = gnm_nan;
2728 sols->vars[i].high = gnm_nan;
2729 sols->vars[i].reduced_cost = gnm_nan;
2730 }
2731
2732 cn = 0;
2733 for (l = sp->constraints; l; l = l->next) {
2734 GnmSolverConstraint *c = l->data;
2735 int i;
2736 gnm_float cl, cr;
2737 GnmCell *lhs, *rhs;
2738
2739 for (i = 0;
2740 gnm_solver_constraint_get_part (c, sp, i,
2741 &lhs, &cl,
2742 &rhs, &cr);
2743 i++) {
2744 cn++;
2745 }
2746 }
2747 sols->constraints = g_new (struct GnmSolverSensitivityConstraints_, cn);
2748 for (i = 0; i < cn; i++) {
2749 sols->constraints[i].low = gnm_nan;
2750 sols->constraints[i].high = gnm_nan;
2751 sols->constraints[i].shadow_price = gnm_nan;
2752 }
2753 }
2754
2755 static void
2756 gnm_solver_sensitivity_finalize (GObject *obj)
2757 {
2758 GnmSolverSensitivity *r = GNM_SOLVER_SENSITIVITY (obj);
2759 g_free (r->vars);
2760 g_free (r->constraints);
2761 gnm_solver_sensitivity_parent_class->finalize (obj);
2762 }
2763
2764 static void
2765 gnm_solver_sensitivity_get_property (GObject *object, guint property_id,
2766 GValue *value, GParamSpec *pspec)
2767 {
2768 GnmSolverSensitivity *sols = (GnmSolverSensitivity *)object;
2769
2770 switch (property_id) {
2771 case SOLS_PROP_SOLVER:
2772 g_value_set_object (value, sols->solver);
2773 break;
2774
2775 default:
2776 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
2777 break;
2778 }
2779 }
2780
2781 static void
2782 gnm_solver_sensitivity_set_property (GObject *object, guint property_id,
2783 GValue const *value, GParamSpec *pspec)
2784 {
2785 GnmSolverSensitivity *sols = (GnmSolverSensitivity *)object;
2786
2787 switch (property_id) {
2788 case SOLS_PROP_SOLVER:
2789 /* We hold no ref. */
2790 sols->solver = g_value_get_object (value);
2791 break;
2792
2793 default:
2794 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
2795 break;
2796 }
2797 }
2798
2799 static void
2800 gnm_solver_sensitivity_class_init (GObjectClass *object_class)
2801 {
2802 gnm_solver_sensitivity_parent_class =
2803 g_type_class_peek_parent (object_class);
2804
2805 object_class->finalize = gnm_solver_sensitivity_finalize;
2806 object_class->constructed = gnm_solver_sensitivity_constructed;
2807 object_class->set_property = gnm_solver_sensitivity_set_property;
2808 object_class->get_property = gnm_solver_sensitivity_get_property;
2809
2810 g_object_class_install_property
2811 (object_class, SOLS_PROP_SOLVER,
2812 g_param_spec_object ("solver",
2813 P_("Solver"),
2814 P_("Solver"),
2815 GNM_SOLVER_TYPE,
2816 GSF_PARAM_STATIC |
2817 G_PARAM_CONSTRUCT_ONLY |
2818 G_PARAM_READWRITE));
2819 }
2820
2821 GSF_CLASS (GnmSolverSensitivity, gnm_solver_sensitivity,
2822 gnm_solver_sensitivity_class_init, NULL, G_TYPE_OBJECT)
2823
2824 GnmSolverSensitivity *
2825 gnm_solver_sensitivity_new (GnmSolver *sol)
2826 {
2827 return g_object_new (GNM_SOLVER_SENSITIVITY_TYPE, "solver", sol, NULL);
2828 }
2829
2830 /* ------------------------------------------------------------------------- */
2831
2832 static GObjectClass *gnm_sub_solver_parent_class;
2833
2834 enum {
2835 SUB_SOL_SIG_CHILD_EXIT,
2836 SUB_SOL_SIG_LAST
2837 };
2838
2839 static guint sub_solver_signals[SUB_SOL_SIG_LAST] = { 0 };
2840
2841 void
2842 gnm_sub_solver_clear (GnmSubSolver *subsol)
2843 {
2844 int i;
2845
2846 if (subsol->child_watch) {
2847 g_source_remove (subsol->child_watch);
2848 subsol->child_watch = 0;
2849 }
2850
2851 if (subsol->child_pid) {
2852 #ifdef G_OS_WIN32
2853 TerminateProcess (subsol->child_pid, 127);
2854 #else
2855 kill (subsol->child_pid, SIGKILL);
2856 #endif
2857 g_spawn_close_pid (subsol->child_pid);
2858 subsol->child_pid = (GPid)0;
2859 }
2860
2861 for (i = 0; i <= 2; i++) {
2862 if (subsol->channel_watches[i]) {
2863 g_source_remove (subsol->channel_watches[i]);
2864 subsol->channel_watches[i] = 0;
2865 }
2866 if (subsol->channels[i]) {
2867 g_io_channel_unref (subsol->channels[i]);
2868 subsol->channels[i] = NULL;
2869 }
2870 if (subsol->fd[i] != -1) {
2871 close (subsol->fd[i]);
2872 subsol->fd[i] = -1;
2873 }
2874 }
2875
2876 if (subsol->program_filename) {
2877 g_unlink (subsol->program_filename);
2878 g_free (subsol->program_filename);
2879 subsol->program_filename = NULL;
2880 }
2881
2882 if (subsol->cell_from_name)
2883 g_hash_table_remove_all (subsol->cell_from_name);
2884
2885 if (subsol->name_from_cell)
2886 g_hash_table_remove_all (subsol->name_from_cell);
2887
2888 if (subsol->constraint_from_name)
2889 g_hash_table_remove_all (subsol->constraint_from_name);
2890 }
2891
2892 static void
2893 gnm_sub_solver_dispose (GObject *obj)
2894 {
2895 GnmSubSolver *subsol = GNM_SUB_SOLVER (obj);
2896
2897 gnm_sub_solver_clear (subsol);
2898
2899 gnm_sub_solver_parent_class->dispose (obj);
2900 }
2901
2902 static void
2903 gnm_sub_solver_finalize (GObject *obj)
2904 {
2905 GnmSubSolver *subsol = GNM_SUB_SOLVER (obj);
2906
2907 /*
2908 * The weird finalization in gnm_lpsolve_final makes it important that
2909 * we leave the object in a state that gnm_sub_solver_clear is happy
2910 * with.
2911 */
2912
2913 g_hash_table_destroy (subsol->cell_from_name);
2914 subsol->cell_from_name = NULL;
2915
2916 g_hash_table_destroy (subsol->name_from_cell);
2917 subsol->name_from_cell = NULL;
2918
2919 g_hash_table_destroy (subsol->constraint_from_name);
2920 subsol->constraint_from_name = NULL;
2921
2922 gnm_sub_solver_parent_class->finalize (obj);
2923 }
2924
2925 static void
2926 gnm_sub_solver_init (GnmSubSolver *subsol)
2927 {
2928 int i;
2929
2930 for (i = 0; i <= 2; i++)
2931 subsol->fd[i] = -1;
2932
2933 subsol->cell_from_name =
2934 g_hash_table_new_full (g_str_hash, g_str_equal,
2935 g_free, NULL);
2936 subsol->name_from_cell =
2937 g_hash_table_new (g_direct_hash, g_direct_equal);
2938
2939 subsol->constraint_from_name =
2940 g_hash_table_new_full (g_str_hash, g_str_equal,
2941 g_free, NULL);
2942 }
2943
2944 static void
2945 cb_child_exit (G_GNUC_UNUSED GPid pid, gint status, GnmSubSolver *subsol)
2946 {
2947 gboolean normal = WIFEXITED (status);
2948 int code;
2949
2950 subsol->child_watch = 0;
2951
2952 if (normal) {
2953 code = WEXITSTATUS (status);
2954 if (gnm_solver_debug ())
2955 g_printerr ("Solver process exited with code %d\n",
2956 code);
2957 #ifndef G_OS_WIN32
2958 } else if (WIFSIGNALED (status)) {
2959 code = WTERMSIG (status);
2960 if (gnm_solver_debug ())
2961 g_printerr ("Solver process received signal %d\n",
2962 code);
2963 #endif
2964 } else {
2965 code = -1;
2966 g_printerr ("Solver process exited with status 0x%x\n",
2967 status);
2968 }
2969
2970 g_signal_emit (subsol, sub_solver_signals[SUB_SOL_SIG_CHILD_EXIT], 0,
2971 normal, code);
2972
2973 if (subsol->child_pid) {
2974 g_spawn_close_pid (subsol->child_pid);
2975 subsol->child_pid = (GPid)0;
2976 }
2977 }
2978
2979 /**
2980 * gnm_sub_solver_spawn: (skip)
2981 */
2982 gboolean
2983 gnm_sub_solver_spawn (GnmSubSolver *subsol,
2984 char **argv,
2985 GSpawnChildSetupFunc child_setup, gpointer setup_data,
2986 GIOFunc io_stdout, gpointer stdout_data,
2987 GIOFunc io_stderr, gpointer stderr_data,
2988 GError **err)
2989 {
2990 GnmSolver *sol = GNM_SOLVER (subsol);
2991 gboolean ok;
2992 GSpawnFlags spflags = G_SPAWN_DO_NOT_REAP_CHILD;
2993 int fd;
2994
2995 g_return_val_if_fail (subsol->child_watch == 0, FALSE);
2996 g_return_val_if_fail (sol->status == GNM_SOLVER_STATUS_PREPARED, FALSE);
2997
2998 if (!g_path_is_absolute (argv[0]))
2999 spflags |= G_SPAWN_SEARCH_PATH;
3000
3001 if (io_stdout == NULL && !gnm_solver_debug ())
3002 spflags |= G_SPAWN_STDOUT_TO_DEV_NULL;
3003
3004 if (gnm_solver_debug ()) {
3005 GString *msg = g_string_new ("Spawning");
3006 int i;
3007 for (i = 0; argv[i]; i++) {
3008 g_string_append_c (msg, ' ');
3009 g_string_append (msg, argv[i]);
3010 }
3011 g_printerr ("%s\n", msg->str);
3012 g_string_free (msg, TRUE);
3013 }
3014
3015 #ifdef G_OS_WIN32
3016 /* Hope for the best... */
3017 child_setup = NULL;
3018 setup_data = NULL;
3019 #endif
3020
3021 ok = g_spawn_async_with_pipes
3022 (g_get_home_dir (), /* PWD */
3023 argv,
3024 NULL, /* environment */
3025 spflags,
3026 child_setup, setup_data,
3027 &subsol->child_pid,
3028 NULL, /* stdin */
3029 io_stdout ? &subsol->fd[1] : NULL, /* stdout */
3030 io_stdout ? &subsol->fd[2] : NULL, /* stderr */
3031 err);
3032 if (!ok)
3033 goto fail;
3034
3035 subsol->child_watch =
3036 g_child_watch_add (subsol->child_pid,
3037 (GChildWatchFunc)cb_child_exit, subsol);
3038
3039 subsol->io_funcs[1] = io_stdout;
3040 subsol->io_funcs_data[1] = stdout_data;
3041 subsol->io_funcs[2] = io_stderr;
3042 subsol->io_funcs_data[2] = stderr_data;
3043
3044 for (fd = 1; fd <= 2; fd++) {
3045 GIOFlags ioflags;
3046
3047 if (subsol->io_funcs[fd] == NULL)
3048 continue;
3049
3050 /*
3051 * Despite the name these are documented to work on Win32.
3052 * Let us hope that is actually true.
3053 */
3054 subsol->channels[fd] = g_io_channel_unix_new (subsol->fd[fd]);
3055 ioflags = g_io_channel_get_flags (subsol->channels[fd]);
3056 g_io_channel_set_flags (subsol->channels[fd],
3057 ioflags | G_IO_FLAG_NONBLOCK,
3058 NULL);
3059 subsol->channel_watches[fd] =
3060 g_io_add_watch (subsol->channels[fd],
3061 G_IO_IN,
3062 subsol->io_funcs[fd],
3063 subsol->io_funcs_data[fd]);
3064 }
3065
3066 gnm_solver_set_status (sol, GNM_SOLVER_STATUS_RUNNING);
3067 return TRUE;
3068
3069 fail:
3070 gnm_sub_solver_clear (subsol);
3071 gnm_solver_set_status (sol, GNM_SOLVER_STATUS_ERROR);
3072 return FALSE;
3073 }
3074
3075 const char *
3076 gnm_sub_solver_name_cell (GnmSubSolver *subsol, GnmCell const *cell,
3077 const char *name)
3078 {
3079 char *name_copy = g_strdup (name);
3080
3081 g_hash_table_insert (subsol->cell_from_name,
3082 name_copy,
3083 (gpointer)cell);
3084 g_hash_table_insert (subsol->name_from_cell,
3085 (gpointer)cell,
3086 name_copy);
3087
3088 return name_copy;
3089 }
3090
3091 GnmCell *
3092 gnm_sub_solver_find_cell (GnmSubSolver *subsol, const char *name)
3093 {
3094 return g_hash_table_lookup (subsol->cell_from_name, name);
3095 }
3096
3097 const char *
3098 gnm_sub_solver_get_cell_name (GnmSubSolver *subsol,
3099 GnmCell const *cell)
3100 {
3101 return g_hash_table_lookup (subsol->name_from_cell, (gpointer)cell);
3102 }
3103
3104 const char *
3105 gnm_sub_solver_name_constraint (GnmSubSolver *subsol,
3106 int cidx,
3107 const char *name)
3108 {
3109 char *name_copy = g_strdup (name);
3110
3111 g_hash_table_insert (subsol->constraint_from_name,
3112 name_copy,
3113 GINT_TO_POINTER (cidx));
3114
3115 return name_copy;
3116 }
3117
3118 int
3119 gnm_sub_solver_find_constraint (GnmSubSolver *subsol, const char *name)
3120 {
3121 gpointer idx;
3122
3123 if (g_hash_table_lookup_extended (subsol->constraint_from_name,
3124 (gpointer)name,
3125 NULL, &idx))
3126 return GPOINTER_TO_INT (idx);
3127 else
3128 return -1;
3129 }
3130
3131
3132 char *
3133 gnm_sub_solver_locate_binary (const char *binary, const char *solver,
3134 const char *url,
3135 WBCGtk *wbcg)
3136 {
3137 GtkWindow *parent;
3138 GtkWidget *dialog;
3139 char *path = NULL;
3140 int res;
3141 GtkFileChooser *fsel;
3142 char *title;
3143
3144 parent = wbcg ? wbcg_toplevel (wbcg) : NULL;
3145 dialog = gtk_message_dialog_new_with_markup
3146 (parent,
3147 GTK_DIALOG_DESTROY_WITH_PARENT,
3148 GTK_MESSAGE_QUESTION,
3149 GTK_BUTTONS_YES_NO,
3150 _("Gnumeric is unable to locate the program <i>%s</i> needed "
3151 "for the <i>%s</i> solver. For more information see %s.\n\n"
3152 "Would you like to locate it yourself?"),
3153 binary, solver, url);
3154 title = g_strdup_printf (_("Unable to locate %s"), binary);
3155 g_object_set (G_OBJECT (dialog),
3156 "title", title,
3157 NULL);
3158 g_free (title);
3159
3160 res = go_gtk_dialog_run (GTK_DIALOG (dialog), parent);
3161 switch (res) {
3162 case GTK_RESPONSE_NO:
3163 case GTK_RESPONSE_DELETE_EVENT:
3164 default:
3165 return NULL;
3166 case GTK_RESPONSE_YES:
3167 break;
3168 }
3169
3170 title = g_strdup_printf (_("Locate the %s program"), binary);
3171 fsel = GTK_FILE_CHOOSER
3172 (g_object_new (GTK_TYPE_FILE_CHOOSER_DIALOG,
3173 "action", GTK_FILE_CHOOSER_ACTION_OPEN,
3174 "local-only", TRUE,
3175 "title", title,
3176 NULL));
3177 g_free (title);
3178 go_gtk_dialog_add_button (GTK_DIALOG (fsel), GNM_STOCK_CANCEL,
3179 "gtk-cancel", GTK_RESPONSE_CANCEL);
3180 go_gtk_dialog_add_button (GTK_DIALOG (fsel), GNM_STOCK_OK,
3181 "system-run", GTK_RESPONSE_OK);
3182 g_object_ref (fsel);
3183 if (go_gtk_file_sel_dialog (parent, GTK_WIDGET (fsel))) {
3184 path = gtk_file_chooser_get_filename (fsel);
3185 if (!g_file_test (path, G_FILE_TEST_IS_EXECUTABLE)) {
3186 g_free (path);
3187 path = NULL;
3188 }
3189 }
3190
3191 gtk_widget_destroy (GTK_WIDGET (fsel));
3192 g_object_unref (fsel);
3193
3194 return path;
3195 }
3196
3197
3198 void
3199 gnm_sub_solver_flush (GnmSubSolver *subsol)
3200 {
3201 int fd;
3202
3203 for (fd = 1; fd <= 2; fd++) {
3204 if (subsol->io_funcs[fd] == NULL)
3205 continue;
3206
3207 subsol->io_funcs[fd] (subsol->channels[fd],
3208 G_IO_IN,
3209 subsol->io_funcs_data[fd]);
3210 }
3211 }
3212
3213 static void
3214 gnm_sub_solver_class_init (GObjectClass *object_class)
3215 {
3216 gnm_sub_solver_parent_class = g_type_class_peek_parent (object_class);
3217
3218 object_class->dispose = gnm_sub_solver_dispose;
3219 object_class->finalize = gnm_sub_solver_finalize;
3220
3221 sub_solver_signals[SUB_SOL_SIG_CHILD_EXIT] =
3222 g_signal_new ("child-exit",
3223 G_OBJECT_CLASS_TYPE (object_class),
3224 G_SIGNAL_RUN_LAST,
3225 G_STRUCT_OFFSET (GnmSubSolverClass, child_exit),
3226 NULL, NULL,
3227 gnm__VOID__BOOLEAN_INT,
3228 G_TYPE_NONE, 2,
3229 G_TYPE_BOOLEAN, G_TYPE_INT);
3230 }
3231
3232 GSF_CLASS (GnmSubSolver, gnm_sub_solver,
3233 gnm_sub_solver_class_init, gnm_sub_solver_init, GNM_SOLVER_TYPE)
3234
3235 /* ------------------------------------------------------------------------- */
3236
3237 enum {
3238 SOL_ITER_SIG_ITERATE,
3239 SOL_ITER_SIG_LAST
3240 };
3241
3242 static guint solver_iterator_signals[SOL_ITER_SIG_LAST] = { 0 };
3243
3244 static void
3245 gnm_solver_iterator_class_init (GObjectClass *object_class)
3246 {
3247 solver_iterator_signals[SOL_ITER_SIG_ITERATE] =
3248 g_signal_new ("iterate",
3249 G_OBJECT_CLASS_TYPE (object_class),
3250 G_SIGNAL_RUN_LAST,
3251 G_STRUCT_OFFSET (GnmSolverIteratorClass, iterate),
3252 NULL, NULL,
3253 gnm__BOOLEAN__VOID,
3254 G_TYPE_BOOLEAN, 0);
3255 }
3256
3257 gboolean
3258 gnm_solver_iterator_iterate (GnmSolverIterator *iter)
3259 {
3260 gboolean progress = FALSE;
3261 g_signal_emit (iter, solver_iterator_signals[SOL_ITER_SIG_ITERATE], 0, &progress);
3262 return progress;
3263 }
3264
3265 /**
3266 * gnm_solver_iterator_new_func: (skip)
3267 */
3268 GnmSolverIterator *
3269 gnm_solver_iterator_new_func (GCallback iterate, gpointer user)
3270 {
3271 GnmSolverIterator *iter;
3272
3273 iter = g_object_new (GNM_SOLVER_ITERATOR_TYPE, NULL);
3274 g_signal_connect (iter, "iterate", G_CALLBACK (iterate), user);
3275 return iter;
3276 }
3277
3278 static gboolean
3279 cb_polish_iter (GnmSolverIterator *iter, GnmIterSolver *isol)
3280 {
3281 GnmSolver *sol = GNM_SOLVER (isol);
3282 const int n = sol->input_cells->len;
3283 gnm_float *dir;
3284 gboolean progress = FALSE;
3285 int c;
3286
3287 dir = g_new0 (gnm_float, n);
3288 for (c = 0; c < n; c++) {
3289 gnm_float s, y, s0, sm, xc = isol->xk[c];
3290
3291 if (xc == 0) {
3292 s0 = 0.5;
3293 sm = 1;
3294 } else {
3295 int e;
3296 (void)gnm_frexp (xc, &e);
3297 s0 = gnm_ldexp (1, e - 10);
3298 if (s0 == 0) s0 = GNM_MIN;
3299 sm = gnm_abs (xc);
3300 }
3301
3302 dir[c] = 1;
3303 s = gnm_solver_line_search (sol, isol->xk, dir, TRUE,
3304 s0, sm, 0.0, &y);
3305 dir[c] = 0;
3306
3307 if (gnm_finite (s) && s != 0) {
3308 isol->xk[c] += s;
3309 isol->yk = y;
3310 progress = TRUE;
3311 }
3312 }
3313 g_free (dir);
3314
3315 if (progress)
3316 gnm_iter_solver_set_solution (isol);
3317
3318 return progress;
3319 }
3320
3321 /**
3322 * gnm_solver_iterator_new_polish:
3323 * @isol: the solver to operate on
3324 *
3325 * Returns: (transfer full): an iterator object that can be used to polish
3326 * a solution by simple axis-parallel movement.
3327 */
3328 GnmSolverIterator *
3329 gnm_solver_iterator_new_polish (GnmIterSolver *isol)
3330 {
3331 return gnm_solver_iterator_new_func (G_CALLBACK (cb_polish_iter), isol);
3332 }
3333
3334
3335 static gboolean
3336 cb_gradient_iter (GnmSolverIterator *iter, GnmIterSolver *isol)
3337 {
3338 GnmSolver *sol = GNM_SOLVER (isol);
3339 const int n = sol->input_cells->len;
3340 gboolean progress = FALSE;
3341 gnm_float s, y;
3342 gnm_float *g;
3343 int i;
3344
3345 /* Search in opposite direction of gradient. */
3346 g = gnm_solver_compute_gradient (sol, isol->xk);
3347 for (i = 0; i < n; i++)
3348 g[i] = -g[i];
3349
3350 s = gnm_solver_line_search (sol, isol->xk, g, FALSE,
3351 1, gnm_pinf, 0.0, &y);
3352 if (s > 0) {
3353 for (i = 0; i < n; i++)
3354 isol->xk[i] += s * g[i];
3355 isol->yk = y;
3356 progress = TRUE;
3357 }
3358
3359 g_free (g);
3360
3361 if (progress)
3362 gnm_iter_solver_set_solution (isol);
3363
3364 return progress;
3365 }
3366
3367 /**
3368 * gnm_solver_iterator_new_gradient:
3369 * @isol: the solver to operate on
3370 *
3371 * Returns: (transfer full): an iterator object that can be used to perform
3372 * a gradient descent step.
3373 */
3374 GnmSolverIterator *
3375 gnm_solver_iterator_new_gradient (GnmIterSolver *isol)
3376 {
3377 return gnm_solver_iterator_new_func (G_CALLBACK (cb_gradient_iter), isol);
3378 }
3379
3380 GSF_CLASS (GnmSolverIterator, gnm_solver_iterator,
3381 gnm_solver_iterator_class_init, NULL, G_TYPE_OBJECT)
3382
3383 /* ------------------------------------------------------------------------- */
3384
3385 enum {
3386 SOLIC_PROP_0,
3387 SOLIC_PROP_CYCLES
3388 };
3389
3390 static GObjectClass *gnm_solver_iterator_compound_parent_class;
3391
3392 /**
3393 * gnm_solver_iterator_compound_add:
3394 * @ic: Compound iterator
3395 * @iter: (transfer full): sub-iterator
3396 * @count: repeat count
3397 *
3398 * Add an iterator to a compound iterator with a given repeat count. As a
3399 * special case, a repeat count of zero means to try the iterator once
3400 * in a cycle, but only if no other sub-iterator has shown any progress so far.
3401 */
3402 void
3403 gnm_solver_iterator_compound_add (GnmSolverIteratorCompound *ic,
3404 GnmSolverIterator *iter,
3405 unsigned count)
3406 {
3407 g_ptr_array_add (ic->iterators, iter);
3408 ic->counts = g_renew (unsigned, ic->counts, ic->iterators->len);
3409 ic->counts[ic->iterators->len - 1] = count;
3410 }
3411
3412 static gboolean
3413 gnm_solver_iterator_compound_iterate (GnmSolverIterator *iter)
3414 {
3415 GnmSolverIteratorCompound *ic = (GnmSolverIteratorCompound *)iter;
3416 gboolean progress;
3417
3418 while (TRUE) {
3419 if (ic->cycle >= ic->cycles)
3420 return FALSE;
3421
3422 if (ic->next >= ic->iterators->len) {
3423 /* We've been through all iterators. */
3424 if (!ic->cycle_progress)
3425 return FALSE;
3426 ic->cycle_progress = FALSE;
3427 ic->next = 0;
3428 ic->next_counter = 0;
3429 ic->cycle++;
3430 continue;
3431 }
3432
3433 if (ic->next_counter < ic->counts[ic->next])
3434 break;
3435
3436 /* Special case: when count==0, use only if no progress. */
3437 if (!ic->cycle_progress && ic->next_counter == 0)
3438 break;
3439
3440 ic->next++;
3441 ic->next_counter = 0;
3442 }
3443
3444 progress = gnm_solver_iterator_iterate (g_ptr_array_index (ic->iterators, ic->next));
3445 if (progress) {
3446 ic->cycle_progress = TRUE;
3447 ic->next_counter++;
3448 } else {
3449 /* No progress, so don't retry. */
3450 ic->next++;
3451 ic->next_counter = 0;
3452 }
3453
3454 /* Report progress as long as we have stuff to try. */
3455 return TRUE;
3456 }
3457
3458 static void
3459 gnm_solver_iterator_compound_init (GnmSolverIteratorCompound *ic)
3460 {
3461 ic->iterators = g_ptr_array_new ();
3462 ic->cycles = G_MAXUINT;
3463 }
3464
3465 static void
3466 gnm_solver_iterator_compound_finalize (GObject *obj)
3467 {
3468 GnmSolverIteratorCompound *ic = (GnmSolverIteratorCompound *)obj;
3469 g_ptr_array_foreach (ic->iterators, (GFunc)g_object_unref, NULL);
3470 g_ptr_array_free (ic->iterators, TRUE);
3471 g_free (ic->counts);
3472 gnm_solver_iterator_compound_parent_class->finalize (obj);
3473 }
3474
3475 static void
3476 gnm_solver_iterator_compound_get_property (GObject *object, guint property_id,
3477 GValue *value, GParamSpec *pspec)
3478 {
3479 GnmSolverIteratorCompound *it = (GnmSolverIteratorCompound *)object;
3480
3481 switch (property_id) {
3482 case SOLIC_PROP_CYCLES:
3483 g_value_set_uint (value, it->cycles);
3484 break;
3485
3486 default:
3487 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
3488 break;
3489 }
3490 }
3491
3492 static void
3493 gnm_solver_iterator_compound_set_property (GObject *object, guint property_id,
3494 GValue const *value, GParamSpec *pspec)
3495 {
3496 GnmSolverIteratorCompound *it = (GnmSolverIteratorCompound *)object;
3497
3498 switch (property_id) {
3499 case SOLIC_PROP_CYCLES:
3500 it->cycles = g_value_get_uint (value);
3501 break;
3502
3503 default:
3504 G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
3505 break;
3506 }
3507 }
3508
3509 static void
3510 gnm_solver_iterator_compound_class_init (GObjectClass *object_class)
3511 {
3512 GnmSolverIteratorClass *iclass = (GnmSolverIteratorClass *)object_class;
3513
3514 gnm_solver_iterator_compound_parent_class = g_type_class_peek_parent (object_class);
3515
3516 object_class->finalize = gnm_solver_iterator_compound_finalize;
3517 object_class->set_property = gnm_solver_iterator_compound_set_property;
3518 object_class->get_property = gnm_solver_iterator_compound_get_property;
3519 iclass->iterate = gnm_solver_iterator_compound_iterate;
3520
3521 g_object_class_install_property (object_class, SOLIC_PROP_CYCLES,
3522 g_param_spec_uint ("cycles",
3523 P_("Cycles"),
3524 P_("Maximum number of cycles"),
3525 0, G_MAXUINT, G_MAXUINT,
3526 GSF_PARAM_STATIC | G_PARAM_READWRITE));
3527 }
3528
3529 GSF_CLASS (GnmSolverIteratorCompound, gnm_solver_iterator_compound,
3530 gnm_solver_iterator_compound_class_init, gnm_solver_iterator_compound_init, GNM_SOLVER_ITERATOR_TYPE)
3531
3532 /* ------------------------------------------------------------------------- */
3533
3534 static GObjectClass *gnm_iter_solver_parent_class;
3535
3536 void
3537 gnm_iter_solver_set_iterator (GnmIterSolver *isol, GnmSolverIterator *iterator)
3538 {
3539 GnmSolverIterator *old_iterator;
3540
3541 g_return_if_fail (GNM_IS_ITER_SOLVER (isol));
3542
3543 old_iterator = isol->iterator;
3544 isol->iterator = iterator ? g_object_ref (iterator) : NULL;
3545 if (old_iterator)
3546 g_object_unref (old_iterator);
3547 }
3548
3549 gboolean
3550 gnm_iter_solver_get_initial_solution (GnmIterSolver *isol, GError **err)
3551 {
3552 GnmSolver *sol = GNM_SOLVER (isol);
3553 const int n = sol->input_cells->len;
3554 int i;
3555
3556 if (gnm_solver_check_constraints (sol))
3557 goto got_it;
3558
3559 /* More? */
3560
3561 g_set_error (err,
3562 go_error_invalid (),
3563 0,
3564 _("The initial values do not satisfy the constraints."));
3565 return FALSE;
3566
3567 got_it:
3568 for (i = 0; i < n; i++) {
3569 GnmCell *cell = g_ptr_array_index (sol->input_cells, i);
3570 isol->xk[i] = value_get_as_float (cell->value);
3571 }
3572 isol->yk = gnm_solver_get_target_value (sol);
3573
3574 gnm_iter_solver_set_solution (isol);
3575
3576 return TRUE;
3577 }
3578
3579 void
3580 gnm_iter_solver_set_solution (GnmIterSolver *isol)
3581 {
3582 GnmSolver *sol = GNM_SOLVER (isol);
3583 GnmSolverResult *result = g_object_new (GNM_SOLVER_RESULT_TYPE, NULL);
3584 const int n = sol->input_cells->len;
3585
3586 result->quality = GNM_SOLVER_RESULT_FEASIBLE;
3587 result->value = sol->flip_sign ? 0 - isol->yk : isol->yk;
3588 result->solution = g_memdup (isol->xk, n * sizeof (gnm_float));
3589 g_object_set (sol, "result", result, NULL);
3590 g_object_unref (result);
3591
3592 if (!gnm_solver_check_constraints (sol)) {
3593 g_printerr ("Infeasible solution set\n");
3594 }
3595 }
3596
3597 static void
3598 gnm_iter_solver_clear (GnmIterSolver *isol)
3599 {
3600 if (isol->idle_tag) {
3601 g_source_remove (isol->idle_tag);
3602 isol->idle_tag = 0;
3603 }
3604 }
3605
3606 static void
3607 gnm_iter_solver_dispose (GObject *obj)
3608 {
3609 GnmIterSolver *isol = GNM_ITER_SOLVER (obj);
3610 gnm_iter_solver_clear (isol);
3611 gnm_iter_solver_parent_class->dispose (obj);
3612 }
3613
3614 static void
3615 gnm_iter_solver_finalize (GObject *obj)
3616 {
3617 GnmIterSolver *isol = GNM_ITER_SOLVER (obj);
3618 g_free (isol->xk);
3619 gnm_iter_solver_parent_class->finalize (obj);
3620 }
3621
3622 static void
3623 gnm_iter_solver_constructed (GObject *obj)
3624 {
3625 GnmIterSolver *isol = GNM_ITER_SOLVER (obj);
3626 GnmSolver *sol = GNM_SOLVER (obj);
3627
3628 /* Chain to parent first */
3629 gnm_iter_solver_parent_class->constructed (obj);
3630
3631 isol->xk = g_new0 (gnm_float, sol->input_cells->len);
3632 }
3633
3634 static void
3635 gnm_iter_solver_init (GnmIterSolver *isol)
3636 {
3637 }
3638
3639 static gint
3640 gnm_iter_solver_idle (gpointer data)
3641 {
3642 GnmIterSolver *isol = data;
3643 GnmSolver *sol = &isol->parent;
3644 GnmSolverParameters *params = sol->params;
3645 gboolean progress;
3646
3647 progress = isol->iterator && gnm_solver_iterator_iterate (isol->iterator);
3648 isol->iterations++;
3649
3650 if (!gnm_solver_finished (sol)) {
3651 if (!progress) {
3652 gnm_solver_set_status (sol, GNM_SOLVER_STATUS_DONE);
3653 } else if (isol->iterations >= params->options.max_iter) {
3654 gnm_solver_stop (sol, NULL);
3655 gnm_solver_set_reason (sol, _("Iteration limit exceeded"));
3656 }
3657 }
3658
3659 if (gnm_solver_finished (sol)) {
3660 isol->idle_tag = 0;
3661
3662 gnm_app_recalc ();
3663
3664 return FALSE;
3665 } else {
3666 /* Call again. */
3667 return TRUE;
3668 }
3669 }
3670
3671 static gboolean
3672 gnm_iter_solver_start (GnmSolver *solver, WorkbookControl *wbc, GError **err)
3673 {
3674 GnmIterSolver *isol = GNM_ITER_SOLVER (solver);
3675
3676 g_return_val_if_fail (isol->idle_tag == 0, FALSE);
3677
3678 isol->idle_tag = g_idle_add (gnm_iter_solver_idle, solver);
3679 gnm_solver_set_status (solver, GNM_SOLVER_STATUS_RUNNING);
3680
3681 return TRUE;
3682 }
3683
3684 static gboolean
3685 gnm_iter_solver_stop (GnmSolver *solver, GError **err)
3686 {
3687 GnmIterSolver *isol = GNM_ITER_SOLVER (solver);
3688 GnmSolver *sol = &isol->parent;
3689
3690 gnm_iter_solver_clear (isol);
3691
3692 g_clear_object (&isol->iterator);
3693
3694 gnm_solver_set_status (sol, GNM_SOLVER_STATUS_CANCELLED);
3695
3696 return TRUE;
3697 }
3698
3699 static void
3700 gnm_iter_solver_class_init (GObjectClass *object_class)
3701 {
3702 GnmSolverClass *sclass = (GnmSolverClass *)object_class;
3703
3704 gnm_iter_solver_parent_class = g_type_class_peek_parent (object_class);
3705
3706 object_class->dispose = gnm_iter_solver_dispose;
3707 object_class->finalize = gnm_iter_solver_finalize;
3708 object_class->constructed = gnm_iter_solver_constructed;
3709 sclass->start = gnm_iter_solver_start;
3710 sclass->stop = gnm_iter_solver_stop;
3711 }
3712
3713 GSF_CLASS (GnmIterSolver, gnm_iter_solver,
3714 gnm_iter_solver_class_init, gnm_iter_solver_init, GNM_SOLVER_TYPE)
3715
3716 /* ------------------------------------------------------------------------- */
3717
3718 static GObjectClass *gnm_solver_factory_parent_class;
3719
3720 static void
3721 gnm_solver_factory_finalize (GObject *obj)
3722 {
3723 GnmSolverFactory *factory = GNM_SOLVER_FACTORY (obj);
3724
3725 if (factory->notify)
3726 factory->notify (factory->data);
3727
3728 g_free (factory->id);
3729 g_free (factory->name);
3730
3731 gnm_solver_factory_parent_class->finalize (obj);
3732 }
3733
3734 static void
3735 gnm_solver_factory_class_init (GObjectClass *object_class)
3736 {
3737 gnm_solver_factory_parent_class =
3738 g_type_class_peek_parent (object_class);
3739
3740 object_class->finalize = gnm_solver_factory_finalize;
3741 }
3742
3743 GSF_CLASS (GnmSolverFactory, gnm_solver_factory,
3744 gnm_solver_factory_class_init, NULL, G_TYPE_OBJECT)
3745
3746
3747 static GSList *solvers;
3748
3749 /**
3750 * gnm_solver_db_get:
3751 *
3752 * Returns: (transfer none) (element-type GnmSolverFactory): list of
3753 * registered solver factories.
3754 */
3755 GSList *
3756 gnm_solver_db_get (void)
3757 {
3758 return solvers;
3759 }
3760
3761 /**
3762 * gnm_solver_factory_new:
3763 * @id: Unique identifier
3764 * @name: Translated name for UI purposes
3765 * @type: Model type created by factory
3766 * @creator: (scope notified): callback for creating a solver
3767 * @functional: (scope notified): callback for checking if factory is functional
3768 * @data: User pointer for @creator and @functional
3769 * @notify: Destroy notification for @data.
3770 *
3771 * Returns: (transfer full): a new #GnmSolverFactory
3772 */
3773 GnmSolverFactory *
3774 gnm_solver_factory_new (const char *id,
3775 const char *name,
3776 GnmSolverModelType type,
3777 GnmSolverCreator creator,
3778 GnmSolverFactoryFunctional functional,
3779 gpointer data,
3780 GDestroyNotify notify)
3781 {
3782 GnmSolverFactory *res;
3783
3784 g_return_val_if_fail (id != NULL, NULL);
3785 g_return_val_if_fail (name != NULL, NULL);
3786 g_return_val_if_fail (creator != NULL, NULL);
3787
3788 res = g_object_new (GNM_SOLVER_FACTORY_TYPE, NULL);
3789 res->id = g_strdup (id);
3790 res->name = g_strdup (name);
3791 res->type = type;
3792 res->creator = creator;
3793 res->functional = functional;
3794 res->data = data;
3795 res->notify = notify;
3796 return res;
3797 }
3798
3799 /**
3800 * gnm_solver_factory_create:
3801 * @factory: #GnmSolverFactory
3802 * @param: #GnmSolverParameters
3803 *
3804 * Returns: (transfer full): a new #GnmSolver
3805 */
3806 GnmSolver *
3807 gnm_solver_factory_create (GnmSolverFactory *factory,
3808 GnmSolverParameters *param)
3809 {
3810 g_return_val_if_fail (GNM_IS_SOLVER_FACTORY (factory), NULL);
3811 return factory->creator (factory, param, factory->data);
3812 }
3813
3814 gboolean
3815 gnm_solver_factory_functional (GnmSolverFactory *factory,
3816 WBCGtk *wbcg)
3817 {
3818 if (factory == NULL)
3819 return FALSE;
3820
3821 return (factory->functional == NULL ||
3822 factory->functional (factory, wbcg, factory->data));
3823 }
3824
3825 static int
3826 cb_compare_factories (GnmSolverFactory *a, GnmSolverFactory *b)
3827 {
3828 return go_utf8_collate_casefold (a->name, b->name);
3829 }
3830
3831 void
3832 gnm_solver_db_register (GnmSolverFactory *factory)
3833 {
3834 if (gnm_solver_debug ())
3835 g_printerr ("Registering %s\n", factory->id);
3836 g_object_ref (factory);
3837 solvers = g_slist_insert_sorted (solvers, factory,
3838 (GCompareFunc)cb_compare_factories);
3839 }
3840
3841 void
3842 gnm_solver_db_unregister (GnmSolverFactory *factory)
3843 {
3844 if (gnm_solver_debug ())
3845 g_printerr ("Unregistering %s\n", factory->id);
3846 solvers = g_slist_remove (solvers, factory);
3847 g_object_unref (factory);
3848 }
3849
3850 /* ------------------------------------------------------------------------- */
gnumeric mirror