[gnumeric] Solver: check linearity of model and constraints.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] Solver: check linearity of model and constraints.
- Date: Tue, 5 May 2015 20:02:02 +0000 (UTC)
commit 3fa0a8e49b60e2a0a3d393bec1b13963ec765703
Author: Morten Welinder <terra gnome org>
Date: Tue May 5 16:01:11 2015 -0400
Solver: check linearity of model and constraints.
This is only implemented for lpsolve for now.
NEWS | 1 +
plugins/lpsolve/lpsolve-write.c | 98 ++++++++------------
src/tools/gnm-solver.c | 195 ++++++++++++++++++++++++++++++++++++++-
src/tools/gnm-solver.h | 12 +++
4 files changed, 243 insertions(+), 63 deletions(-)
---
diff --git a/NEWS b/NEWS
index c3632a6..4743392 100644
--- a/NEWS
+++ b/NEWS
@@ -14,6 +14,7 @@ Morten:
* Solver code refactoring.
* Plug leaks.
* Fuzzed file fixes. [#748595] [#748597]
+ * Make solver check linearity of model.
--------------------------------------------------------------------------
Gnumeric 1.12.22
diff --git a/plugins/lpsolve/lpsolve-write.c b/plugins/lpsolve/lpsolve-write.c
index f3555e2..72dd5b8 100644
--- a/plugins/lpsolve/lpsolve-write.c
+++ b/plugins/lpsolve/lpsolve-write.c
@@ -36,44 +36,6 @@
#include <string.h>
-static gboolean
-gnm_solver_get_lp_coeff (GnmCell *target, GnmCell *cell,
- gnm_float *x, GError **err)
-{
- gnm_float x0, x1;
- gboolean res = FALSE;
-
- gnm_cell_eval (target);
- if (!VALUE_IS_NUMBER (target->value))
- goto fail;
- x0 = value_get_as_float (target->value);
-
- gnm_cell_set_value (cell, value_new_float (1));
- cell_queue_recalc (cell);
- gnm_cell_eval (target);
- if (!VALUE_IS_NUMBER (target->value))
- goto fail;
- x1 = value_get_as_float (target->value);
-
- *x = x1 - x0;
- res = TRUE;
- goto out;
-
-fail:
- g_set_error (err,
- go_error_invalid (),
- 0,
- _("Target cell did not evaluate to a number."));
- *x = 0;
-
-out:
- gnm_cell_set_value (cell, value_new_int (0));
- cell_queue_recalc (cell);
- gnm_cell_eval (target);
-
- return res;
-}
-
static const char *
lpsolve_var_name (GnmSubSolver *ssol, GnmCell const *cell)
{
@@ -88,38 +50,37 @@ lpsolve_var_name (GnmSubSolver *ssol, GnmCell const *cell)
static gboolean
lpsolve_affine_func (GString *dst, GnmCell *target, GnmSubSolver *ssol,
+ gnm_float const *x1, gnm_float const *x2, guint8 const *pdisc,
gnm_float cst, GError **err)
{
GnmSolver *sol = GNM_SOLVER (ssol);
unsigned ui;
gboolean any = FALSE;
gnm_float y;
- GPtrArray *old_values;
gboolean ok = TRUE;
GPtrArray *input_cells = sol->input_cells;
+ gnm_float *cs;
if (!target) {
gnm_string_add_number (dst, cst);
return TRUE;
}
- old_values = g_ptr_array_new ();
- for (ui = 0; ui < input_cells->len; ui++) {
- GnmCell *cell = g_ptr_array_index (input_cells, ui);
- g_ptr_array_add (old_values, value_dup (cell->value));
- gnm_cell_set_value (cell, value_new_int (0));
- cell_queue_recalc (cell);
- }
-
+ gnm_solver_set_vars (sol, x1);
gnm_cell_eval (target);
y = cst + value_get_as_float (target->value);
+ cs = gnm_solver_get_lp_coeffs (sol, target, x1, x2, pdisc, err);
+ if (!cs)
+ goto fail;
+
+ /* Adjust constant for choice of x1. */
+ for (ui = 0; ui < input_cells->len; ui++)
+ y -= x1[ui] * cs[ui];
+
for (ui = 0; ui < input_cells->len; ui++) {
GnmCell *cell = g_ptr_array_index (input_cells, ui);
- gnm_float x;
- ok = gnm_solver_get_lp_coeff (target, cell, &x, err);
- if (!ok)
- goto fail;
+ gnm_float x = cs[ui];
if (x == 0)
continue;
@@ -154,13 +115,7 @@ lpsolve_affine_func (GString *dst, GnmCell *target, GnmSubSolver *ssol,
}
fail:
- for (ui = 0; ui < input_cells->len; ui++) {
- GnmCell *cell = g_ptr_array_index (input_cells, ui);
- GnmValue *old = g_ptr_array_index (old_values, ui);
- gnm_cell_set_value (cell, old);
- cell_queue_recalc (cell);
- }
- g_ptr_array_free (old_values, TRUE);
+ g_free (cs);
return ok;
}
@@ -178,6 +133,9 @@ lpsolve_create_program (GnmSubSolver *ssol, GOIOContext *io_context, GError **er
GnmCell *target_cell = gnm_solver_param_get_target_cell (sp);
GPtrArray *input_cells = sol->input_cells;
gsize progress;
+ GPtrArray *old = NULL;
+ gnm_float *x1 = NULL, *x2 = NULL;
+ guint8 *pdisc = NULL;
/* ---------------------------------------- */
@@ -191,7 +149,7 @@ lpsolve_create_program (GnmSubSolver *ssol, GOIOContext *io_context, GError **er
/* ---------------------------------------- */
- progress = 2;
+ progress = 3;
if (sp->options.assume_non_negative) progress++;
if (sp->options.assume_discrete) progress++;
progress += g_slist_length (sp->constraints);
@@ -200,6 +158,13 @@ lpsolve_create_program (GnmSubSolver *ssol, GOIOContext *io_context, GError **er
/* ---------------------------------------- */
+ old = gnm_solver_save_vars (sol);
+
+ gnm_solver_pick_lp_coords (sol, &x1, &x2, &pdisc);
+ go_io_count_progress_update (io_context, 1);
+
+ /* ---------------------------------------- */
+
switch (sp->problem_type) {
case GNM_SOLVER_MINIMIZE:
g_string_append (objfunc, "min: ");
@@ -213,6 +178,7 @@ lpsolve_create_program (GnmSubSolver *ssol, GOIOContext *io_context, GError **er
go_io_count_progress_update (io_context, 1);
if (!lpsolve_affine_func (objfunc, target_cell, ssol,
+ x1, x2, pdisc,
0, err))
goto fail;
g_string_append (objfunc, ";\n");
@@ -285,7 +251,9 @@ lpsolve_create_program (GnmSubSolver *ssol, GOIOContext *io_context, GError **er
gboolean ok;
ok = lpsolve_affine_func
- (constraints, lhs, ssol, cl, err);
+ (constraints, lhs, ssol,
+ x1, x2, pdisc,
+ cl, err);
if (!ok)
goto fail;
@@ -294,7 +262,9 @@ lpsolve_create_program (GnmSubSolver *ssol, GOIOContext *io_context, GError **er
g_string_append_c (constraints, ' ');
ok = lpsolve_affine_func
- (constraints, rhs, ssol, cr, err);
+ (constraints, rhs, ssol,
+ x1, x2, pdisc,
+ cr, err);
if (!ok)
goto fail;
@@ -323,6 +293,12 @@ fail:
g_string_free (objfunc, TRUE);
g_string_free (constraints, TRUE);
g_string_free (declarations, TRUE);
+ g_free (x1);
+ g_free (x2);
+ g_free (pdisc);
+
+ if (old)
+ gnm_solver_restore_vars (sol, old);
return prg;
}
diff --git a/src/tools/gnm-solver.c b/src/tools/gnm-solver.c
index b190eec..e6c0f25 100644
--- a/src/tools/gnm-solver.c
+++ b/src/tools/gnm-solver.c
@@ -1370,7 +1370,9 @@ cell_is_constant (GnmCell *cell, gnm_float *pc)
static void
-gnm_solver_get_limits (GnmSolver *solver, gnm_float **pmin, gnm_float **pmax)
+gnm_solver_get_limits (GnmSolver *solver,
+ gnm_float **pmin, gnm_float **pmax,
+ guint8 **pdisc)
{
GnmValue const *vinput;
GnmSolverParameters *params = solver->params;
@@ -1379,15 +1381,18 @@ gnm_solver_get_limits (GnmSolver *solver, gnm_float **pmin, gnm_float **pmax)
unsigned ui;
*pmin = *pmax = NULL;
+ *pdisc = NULL;
vinput = gnm_solver_param_get_input (params);
if (!vinput) return;
*pmin = g_new (gnm_float, n);
*pmax = g_new (gnm_float, n);
+ *pdisc = g_new (guint8, n);
for (ui = 0; ui < n; ui++) {
(*pmin)[ui] = params->options.assume_non_negative ? 0 : gnm_ninf;
(*pmax)[ui] = gnm_pinf;
+ (*pdisc)[ui] = params->options.assume_discrete;
}
for (l = params->constraints; l; l = l->next) {
@@ -1410,8 +1415,10 @@ gnm_solver_get_limits (GnmSolver *solver, gnm_float **pmin, gnm_float **pmax)
switch (c->type) {
case GNM_SOLVER_INTEGER:
+ (*pdisc)[idx] = TRUE;
break;
case GNM_SOLVER_BOOLEAN:
+ (*pdisc)[idx] = TRUE;
SET_LOWER (0.0);
SET_UPPER (1.0);
break;
@@ -1432,6 +1439,17 @@ gnm_solver_get_limits (GnmSolver *solver, gnm_float **pmin, gnm_float **pmax)
}
}
}
+
+ /*
+ * If parameters are discrete, narrow the range by eliminating
+ * the fractional part of the limits.
+ */
+ for (ui = 0; ui < n; ui++) {
+ if ((*pdisc)[ui]) {
+ (*pmin)[ui] = gnm_ceil ((*pmin)[ui]);
+ (*pmax)[ui] = gnm_floor ((*pmax)[ui]);
+ }
+ }
}
#undef SET_LOWER
@@ -1530,6 +1548,7 @@ gnm_solver_create_report (GnmSolver *solver, const char *name)
vinput = gnm_solver_param_get_input (params);
if (vinput) {
gnm_float *pmin, *pmax;
+ guint8 *pdisc;
unsigned ui;
ADD_HEADER (_("Variables"));
@@ -1541,7 +1560,7 @@ gnm_solver_create_report (GnmSolver *solver, const char *name)
dao_set_cell (dao, 5, R, _("Slack"));
R++;
- gnm_solver_get_limits (solver, &pmin, &pmax);
+ gnm_solver_get_limits (solver, &pmin, &pmax, &pdisc);
for (ui = 0; ui < solver->input_cells->len; ui++) {
GnmCell *cell = g_ptr_array_index (solver->input_cells, ui);
@@ -1574,6 +1593,7 @@ gnm_solver_create_report (GnmSolver *solver, const char *name)
g_free (pmin);
g_free (pmax);
+ g_free (pdisc);
R++;
}
@@ -1709,6 +1729,35 @@ gnm_solver_set_vars (GnmSolver *sol, gnm_float const *xs)
gnm_solver_set_var (sol, i, xs[i]);
}
+GPtrArray *
+gnm_solver_save_vars (GnmSolver *sol)
+{
+ GPtrArray *vals = g_ptr_array_new ();
+ unsigned ui;
+
+ for (ui = 0; ui < sol->input_cells->len; ui++) {
+ GnmCell *cell = g_ptr_array_index (sol->input_cells, ui);
+ g_ptr_array_add (vals, value_dup (cell->value));
+ }
+
+ return vals;
+}
+
+void
+gnm_solver_restore_vars (GnmSolver *sol, GPtrArray *vals)
+{
+ unsigned ui;
+
+ for (ui = 0; ui < sol->input_cells->len; ui++) {
+ GnmCell *cell = g_ptr_array_index (sol->input_cells, ui);
+ gnm_cell_set_value (cell, g_ptr_array_index (vals, ui));
+ cell_queue_recalc (cell);
+ }
+
+ g_ptr_array_free (vals, TRUE);
+}
+
+
/**
* gnm_solver_compute_gradient:
* @sol: Solver
@@ -1923,6 +1972,148 @@ gnm_solver_line_search (GnmSolver *sol,
return s1;
}
+/**
+ * gnm_solver_pick_lp_coords:
+ * @sol: Solver
+ * @px1: (out): first coordinate value
+ * @px2: (out): second coordinate value
+ * @pdisc: (out): indicator for discrete coordinate
+ *
+ * Pick two good values for each coordinate. We prefer 0 and 1
+ * when they are valid.
+ */
+void
+gnm_solver_pick_lp_coords (GnmSolver *sol,
+ gnm_float **px1, gnm_float **px2,
+ guint8 **pdisc)
+{
+ const unsigned n = sol->input_cells->len;
+ gnm_float *pmin, *pmax;
+ gnm_float *x1 = g_new (gnm_float, n);
+ gnm_float *x2 = g_new (gnm_float, n);
+ unsigned ui;
+
+ gnm_solver_get_limits (sol, &pmin, &pmax, pdisc);
+
+ for (ui = 0; ui < sol->input_cells->len; ui++) {
+ const gnm_float L = pmin[ui], H = pmax[ui];
+
+ if (L == H) {
+ x1[ui] = x2[ui] = L;
+ } else if ((*pdisc)[ui] && H - L == 1) {
+ x1[ui] = L;
+ x2[ui] = H;
+ } else {
+ if (L <= 0 && H >= 0)
+ x1[ui] = 0;
+ else if (gnm_finite (L))
+ x1[ui] = L;
+ else
+ x1[ui] = H;
+
+ if (x1[ui] + 1 <= H)
+ x2[ui] = x1[ui] + 1;
+ else if (x1[ui] - 1 >= H)
+ x2[ui] = x1[ui] - 1;
+ else if (x1[ui] != H)
+ x2[ui] = (x1[ui] + H) / 2;
+ else
+ x2[ui] = (x1[ui] + L) / 2;
+ }
+ }
+
+ g_free (pmin);
+ g_free (pmax);
+ *px1 = x1;
+ *px2 = x2;
+}
+
+/**
+ * gnm_solver_get_lp_coeffs:
+ * @sol: Solver
+ * @ycell: Cell for which to compute coefficients
+ * @x1: first coordinate value
+ * @x2: second coordinate value
+ * @pdisc: indicator for discrete coordinate
+ * @err: error location
+ *
+ * Returns: (transfer full): coordinates, or NULL in case of error.
+ * Note: this function is not affected by the flip-sign property, even
+ * if @ycell happens to coindice with the solver target cell.
+ */
+gnm_float *
+gnm_solver_get_lp_coeffs (GnmSolver *sol, GnmCell *ycell,
+ gnm_float const *x1, gnm_float const *x2,
+ guint8 const *pdisc,
+ GError **err)
+{
+ const unsigned n = sol->input_cells->len;
+ unsigned ui;
+ gnm_float *res = g_new (gnm_float, n);
+ gnm_float y0;
+
+ gnm_solver_set_vars (sol, x1);
+ gnm_cell_eval (ycell);
+ y0 = VALUE_IS_NUMBER (ycell->value) ? value_get_as_float (ycell->value) : gnm_nan;
+ if (!gnm_finite (y0))
+ goto fail_calc;
+
+ for (ui = 0; ui < sol->input_cells->len; ui++) {
+ gnm_float dx = x2[ui] - x1[ui], dy, y1;
+
+ if (dx <= 0) {
+ res[ui] = 0;
+ continue;
+ }
+
+ gnm_solver_set_var (sol, ui, x2[ui]);
+ gnm_cell_eval (ycell);
+ y1 = VALUE_IS_NUMBER (ycell->value) ? value_get_as_float (ycell->value) : gnm_nan;
+
+ dy = y1 - y0;
+ res[ui] = dy / dx;
+ if (!gnm_finite (res[ui]))
+ goto fail_calc;
+
+ if (!pdisc[ui] || dx != 1) {
+ gnm_float x01, y01, e, emax;
+
+ x01 = (x1[ui] + x2[ui]) / 2;
+ if (pdisc[ui]) x01 = gnm_floor (x01);
+ gnm_solver_set_var (sol, ui, x01);
+ gnm_cell_eval (ycell);
+ y01 = VALUE_IS_NUMBER (ycell->value) ? value_get_as_float (ycell->value) : gnm_nan;
+ if (!gnm_finite (y01))
+ goto fail_calc;
+
+ emax = dy == 0 ? 1e-10 : gnm_abs (dy) / 1e-10;
+ e = dy - 2 * (y01 - y0);
+ if (gnm_abs (e) > emax)
+ goto fail_linear;
+ }
+
+ gnm_solver_set_var (sol, ui, x1[ui]);
+ }
+
+ return res;
+
+fail_calc:
+ g_set_error (err,
+ go_error_invalid (),
+ 0,
+ _("Target cell did not evaluate to a number."));
+ g_free (res);
+ return NULL;
+
+fail_linear:
+ g_set_error (err,
+ go_error_invalid (),
+ 0,
+ _("Target cell does not appear to depend linearly on input cells."));
+ g_free (res);
+ return NULL;
+}
+
static void
gnm_solver_class_init (GObjectClass *object_class)
diff --git a/src/tools/gnm-solver.h b/src/tools/gnm-solver.h
index bd454e4..e3e9c9c 100644
--- a/src/tools/gnm-solver.h
+++ b/src/tools/gnm-solver.h
@@ -257,6 +257,9 @@ gnm_float gnm_solver_get_target_value (GnmSolver *solver);
void gnm_solver_set_var (GnmSolver *sol, int i, gnm_float x);
void gnm_solver_set_vars (GnmSolver *sol, gnm_float const *xs);
+GPtrArray *gnm_solver_save_vars (GnmSolver *sol);
+void gnm_solver_restore_vars (GnmSolver *sol, GPtrArray *vals);
+
gnm_float *gnm_solver_compute_gradient (GnmSolver *sol, gnm_float const *xs);
gnm_float gnm_solver_line_search (GnmSolver *sol,
@@ -265,6 +268,15 @@ gnm_float gnm_solver_line_search (GnmSolver *sol,
gnm_float step, gnm_float max_step, gnm_float eps,
gnm_float *py);
+void gnm_solver_pick_lp_coords (GnmSolver *sol,
+ gnm_float **px1, gnm_float **px2,
+ guint8 **pdisc);
+
+gnm_float *gnm_solver_get_lp_coeffs (GnmSolver *sol, GnmCell *ycell,
+ gnm_float const *x1, gnm_float const *x2,
+ guint8 const *pdisc,
+ GError **err);
+
/* ------------------------------------------------------------------------- */
/* Solver subclass for subprocesses. */
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
