[gnumeric] Solver: Add sensitivity report.



commit 5177f07bdbd40bd17632625a4caece0835bb51e4
Author: Morten Welinder <terra gnome org>
Date:   Wed Feb 10 21:31:13 2016 -0500

    Solver: Add sensitivity report.
    
    This only works for lpsolve.  Glpk is doable, but it is unlikely
    to ever work for a non-linear model.

 NEWS                            |    1 +
 plugins/lpsolve/ChangeLog       |    7 +
 plugins/lpsolve/gnm-lpsolve.c   |  157 ++++++++++++++++-
 plugins/lpsolve/lpsolve-write.c |    9 +-
 src/dialogs/dialog-solver.c     |   11 +-
 src/dialogs/solver.ui           |   17 ++-
 src/gnumeric-fwd.h              |    1 +
 src/tools/ChangeLog             |    4 +
 src/tools/gnm-solver.c          |  381 +++++++++++++++++++++++++++++++++++----
 src/tools/gnm-solver.h          |   44 +++++-
 src/xml-sax-read.c              |    3 +-
 src/xml-sax-write.c             |    2 +
 12 files changed, 585 insertions(+), 52 deletions(-)
---
diff --git a/NEWS b/NEWS
index 167b1af..3fec7e5 100644
--- a/NEWS
+++ b/NEWS
@@ -4,6 +4,7 @@ Morten:
        * Fuzzed file fixes.  [#761663] [#761727]
        * Plug leak.
        * Fix problems with ssconvert to lp/cplex formats.
+       * Add sensitivity report to solver.  (lpsolve only for now.)
 
 --------------------------------------------------------------------------
 Gnumeric 1.12.27
diff --git a/plugins/lpsolve/ChangeLog b/plugins/lpsolve/ChangeLog
index e5f8e1f..cefc5c2 100644
--- a/plugins/lpsolve/ChangeLog
+++ b/plugins/lpsolve/ChangeLog
@@ -1,5 +1,12 @@
 2016-02-10  Morten Welinder  <terra gnome org>
 
+       * lpsolve-write.c (lpsolve_create_program): Name constraints to
+       avoid a clash between default names ("Rn") and cell names.
+
+       * gnm-lpsolve.c (cb_read_stdout): Plug leak.  Read sensitivity
+       results.
+       (gnm_lpsolve_start): Run with -S6 get get sensitivity results.
+
        * lpsolve-write.c (lpsolve_file_save): Handle the situation where
        there is no assigned solver -- i.e., plain ssconvert -- by
        creating a temporary.
diff --git a/plugins/lpsolve/gnm-lpsolve.c b/plugins/lpsolve/gnm-lpsolve.c
index 3210d02..6ebcfdc 100644
--- a/plugins/lpsolve/gnm-lpsolve.c
+++ b/plugins/lpsolve/gnm-lpsolve.c
@@ -13,11 +13,14 @@
 #define SOLVER_PROGRAM "lp_solve"
 #define SOLVER_URL "http://sourceforge.net/projects/lpsolve/";
 #define PRIVATE_KEY "::lpsolve::"
+#define SOLVER_INF 1e30
 
 typedef struct {
        GnmSubSolver *parent;
        GnmSolverResult *result;
-       enum { SEC_UNKNOWN, SEC_VALUES } section;
+       GnmSolverSensitivity *sensitivity;
+       enum { SEC_UNKNOWN, SEC_VALUES,
+              SEC_LIMITS, SEC_DUAL_LIMITS } section;
 } GnmLPSolve;
 
 static void
@@ -29,6 +32,11 @@ gnm_lpsolve_cleanup (GnmLPSolve *lp)
                g_object_unref (lp->result);
                lp->result = NULL;
        }
+
+       if (lp->sensitivity) {
+               g_object_unref (lp->sensitivity);
+               lp->sensitivity = NULL;
+       }
 }
 
 static void
@@ -61,14 +69,18 @@ static GnmSolverResult *
 gnm_lpsolve_start_solution (GnmLPSolve *lp)
 {
        int n;
+       GnmSolver *sol;
 
        g_return_val_if_fail (lp->result == NULL, NULL);
 
-       n = GNM_SOLVER (lp->parent)->input_cells->len;
+       sol = GNM_SOLVER (lp->parent);
+       n = sol->input_cells->len;
 
        lp->result = g_object_new (GNM_SOLVER_RESULT_TYPE, NULL);
        lp->result->solution = g_new0 (gnm_float, n);
 
+       lp->sensitivity = gnm_solver_sensitivity_new (sol);
+
        return lp->result;
 }
 
@@ -80,6 +92,44 @@ gnm_lpsolve_flush_solution (GnmLPSolve *lp)
                g_object_unref (lp->result);
                lp->result = NULL;
        }
+       g_clear_object (&lp->sensitivity);
+}
+
+
+static char **
+my_strsplit (const char *line)
+{
+       GPtrArray *res = g_ptr_array_new ();
+
+       while (1) {
+               const char *end;
+
+               while (g_ascii_isspace (*line))
+                       line++;
+
+               if (!*line)
+                       break;
+
+               end = line;
+               while (*end && !g_ascii_isspace (*end))
+                       end++;
+
+               g_ptr_array_add (res, g_strndup (line, end - line));
+               line = end;
+       }
+       g_ptr_array_add (res, NULL);
+
+       return (char **)g_ptr_array_free (res, FALSE);
+}
+
+static double
+fixup_inf (double v)
+{
+       if (v <= -SOLVER_INF)
+               return go_ninf;
+       if (v >= +SOLVER_INF)
+               return go_pinf;
+       return v;
 }
 
 
@@ -91,12 +141,19 @@ cb_read_stdout (GIOChannel *channel, GIOCondition cond, GnmLPSolve *lp)
        size_t obj_line_len = sizeof (obj_line_prefix) - 1;
        const char val_header_line[] = "Actual values of the variables:";
        size_t val_header_len = sizeof (val_header_line) - 1;
+       const char limit_header_line[] = "Objective function limits:";
+       size_t limit_header_len = sizeof (limit_header_line) - 1;
+       const char dual_limit_header_line[] = "Dual values with from - till limits:";
+       size_t dual_limit_header_len = sizeof (dual_limit_header_line) - 1;
+       gchar *line = NULL;
 
        do {
                GIOStatus status;
-               gchar *line = NULL;
                gsize tpos;
 
+               g_free (line);
+               line = NULL;
+
                status = g_io_channel_read_line (channel,
                                                 &line, NULL, &tpos,
                                                 NULL);
@@ -105,7 +162,7 @@ cb_read_stdout (GIOChannel *channel, GIOCondition cond, GnmLPSolve *lp)
 
                line[tpos] = 0;
 
-               if (line[0] == 0 || g_ascii_isspace (line[0]))
+               if (line[0] == 0)
                        lp->section = SEC_UNKNOWN;
                else if (lp->section == SEC_UNKNOWN &&
                         !strncmp (line, obj_line_prefix, obj_line_len)) {
@@ -117,6 +174,12 @@ cb_read_stdout (GIOChannel *channel, GIOCondition cond, GnmLPSolve *lp)
                } else if (lp->section == SEC_UNKNOWN &&
                           !strncmp (line, val_header_line, val_header_len)) {
                        lp->section = SEC_VALUES;
+               } else if (lp->section == SEC_UNKNOWN &&
+                          !strncmp (line, limit_header_line, limit_header_len)) {
+                       lp->section = SEC_LIMITS;
+               } else if (lp->section == SEC_UNKNOWN &&
+                          !strncmp (line, dual_limit_header_line, dual_limit_header_len)) {
+                       lp->section = SEC_DUAL_LIMITS;
                } else if (lp->section == SEC_VALUES && lp->result) {
                        GnmSolverResult *r = lp->result;
                        double v;
@@ -140,10 +203,88 @@ cb_read_stdout (GIOChannel *channel, GIOCondition cond, GnmLPSolve *lp)
 
                        v = g_ascii_strtod (space + 1, NULL);
                        r->solution[idx] = v;
+               } else if (lp->section == SEC_LIMITS) {
+                       double low, high;
+                       GnmCell *cell;
+                       int idx;
+                       gchar **items;
+
+                       if (g_ascii_isspace (line[0]))
+                               continue;
+
+                       items = my_strsplit (line);
+
+                       if (g_strv_length (items) != 4)
+                               goto bad_limit;
+
+                       cell = gnm_sub_solver_find_cell (lp->parent, items[0]);
+                       idx = gnm_solver_cell_index (sol, cell);
+                       if (idx < 0)
+                               goto bad_limit;
+
+                       low = fixup_inf (g_ascii_strtod (items[1], NULL));
+                       high = fixup_inf (g_ascii_strtod (items[2], NULL));
+
+                       lp->sensitivity->vars[idx].low = low;
+                       lp->sensitivity->vars[idx].high = high;
+
+                       g_strfreev (items);
+
+                       continue;
+
+               bad_limit:
+                       g_printerr ("Strange limit line in output: %s\n",
+                                   line);
+                       lp->section = SEC_UNKNOWN;
+                       g_strfreev (items);
+               } else if (lp->section == SEC_DUAL_LIMITS) {
+                       double dual, low, high;
+                       GnmCell *cell;
+                       int idx, cidx;
+                       gchar **items;
+
+                       if (g_ascii_isspace (line[0]))
+                               continue;
+
+                       items = my_strsplit (line);
+
+                       if (g_strv_length (items) != 4)
+                               goto bad_dual;
+
+                       cell = gnm_sub_solver_find_cell (lp->parent, items[0]);
+                       idx = gnm_solver_cell_index (sol, cell);
+
+                       cidx = (idx == -1)
+                               ? gnm_sub_solver_find_constraint (lp->parent, items[0])
+                               : -1;
+
+                       dual = fixup_inf (g_ascii_strtod (items[1], NULL));
+                       low = fixup_inf (g_ascii_strtod (items[2], NULL));
+                       high = fixup_inf (g_ascii_strtod (items[3], NULL));
+
+                       if (idx >= 0) {
+                               lp->sensitivity->vars[idx].reduced_cost = dual;
+                       } else if (cidx >= 0) {
+                               lp->sensitivity->constraints[cidx].low = low;
+                               lp->sensitivity->constraints[cidx].high = high;
+                               lp->sensitivity->constraints[cidx].shadow_price = dual;
+                       } else {
+                               // Ignore
+                       }
+
+                       g_strfreev (items);
+                       continue;
+
+               bad_dual:
+                       g_printerr ("Strange dual limit line in output: %s\n",
+                                   line);
+                       lp->section = SEC_UNKNOWN;
+                       g_strfreev (items);
                }
-               g_free (line);
        } while (1);
 
+       g_free (line);
+
        return TRUE;
 }
 
@@ -166,6 +307,9 @@ gnm_lpsolve_child_exit (GnmSubSolver *subsol, gboolean normal, int code,
                        gnm_sub_solver_flush (subsol);
                        if (lp->result)
                                lp->result->quality = GNM_SOLVER_RESULT_OPTIMAL;
+                       g_object_set (lp->parent,
+                                     "sensitivity", lp->sensitivity,
+                                     NULL);
                        gnm_lpsolve_flush_solution (lp);
                        break;
 
@@ -249,7 +393,7 @@ gnm_lpsolve_start (GnmSolver *sol, WorkbookControl *wbc, GError **err,
 {
        GnmSubSolver *subsol = GNM_SUB_SOLVER (sol);
        gboolean ok;
-       gchar *argv[5];
+       gchar *argv[6];
        int argc = 0;
        GnmSolverParameters *param = sol->params;
        const char *binary;
@@ -265,6 +409,7 @@ gnm_lpsolve_start (GnmSolver *sol, WorkbookControl *wbc, GError **err,
        argv[argc++] = (gchar *)(param->options.automatic_scaling
                                 ? "-s1"
                                 : "-s0");
+       argv[argc++] = (gchar *)"-S6";
        argv[argc++] = subsol->program_filename;
        argv[argc] = NULL;
        g_assert (argc < (int)G_N_ELEMENTS (argv));
diff --git a/plugins/lpsolve/lpsolve-write.c b/plugins/lpsolve/lpsolve-write.c
index 06f8681..a36fff5 100644
--- a/plugins/lpsolve/lpsolve-write.c
+++ b/plugins/lpsolve/lpsolve-write.c
@@ -213,6 +213,7 @@ lpsolve_create_program (GnmSubSolver *ssol, GOIOContext *io_context, GError **er
                const char *op = NULL;
                const char *type = NULL;
                int i;
+               int cidx = 0;
                gnm_float cl, cr;
                GnmCell *lhs, *rhs;
 
@@ -240,7 +241,7 @@ lpsolve_create_program (GnmSubSolver *ssol, GOIOContext *io_context, GError **er
                     gnm_solver_constraint_get_part (c, sp, i,
                                                     &lhs, &cl,
                                                     &rhs, &cr);
-                    i++) {
+                    i++, cidx++) {
                        if (type) {
                                g_string_append (declarations, type);
                                g_string_append_c (declarations, ' ');
@@ -249,6 +250,12 @@ lpsolve_create_program (GnmSubSolver *ssol, GOIOContext *io_context, GError **er
                        } else {
                                gboolean ok;
 
+                               char *name = g_strdup_printf ("CONSTR_%d", cidx);
+                               g_string_append (constraints, name);
+                               g_string_append (constraints, ": ");
+                               gnm_sub_solver_name_constraint (ssol, cidx, name);
+                               g_free (name);
+
                                ok = lpsolve_affine_func
                                        (constraints, lhs, ssol,
                                         x1, x2,
diff --git a/src/dialogs/dialog-solver.c b/src/dialogs/dialog-solver.c
index 8226f75..ae69c0e 100644
--- a/src/dialogs/dialog-solver.c
+++ b/src/dialogs/dialog-solver.c
@@ -431,6 +431,7 @@ extract_settings (SolverState *state)
        GET_BOOL_ENTRY ("non_neg_button", options.assume_non_negative);
        GET_BOOL_ENTRY ("all_int_button", options.assume_discrete);
        GET_BOOL_ENTRY ("program", options.program_report);
+       GET_BOOL_ENTRY ("sensitivity", options.sensitivity_report);
 
        g_free (param->options.scenario_name);
        param->options.scenario_name = g_strdup
@@ -682,7 +683,7 @@ static void
 create_report (GnmSolver *sol, SolverState *state)
 {
        Sheet *sheet = state->sheet;
-       char *base = g_strdup_printf ("%s Report", sheet->name_unquoted);
+       char *base = g_strdup_printf (_("%s %%s Report"), sheet->name_unquoted);
        gnm_solver_create_report (sol, base);
        g_free (base);
 }
@@ -769,14 +770,11 @@ run_solver (SolverState *state, GnmSolverParameters *param)
                gnm_solver_store_result (sol);
                redo = clipboard_copy_range_undo (sr.sheet, &sr.range);
 
-               if (param->options.program_report) {
+               if (param->options.program_report ||
+                   param->options.sensitivity_report) {
                        Workbook *wb = param->sheet->workbook;
                        GOUndo *undo_report, *redo_report;
 
-                       /* This is a bit of overkill -- it just removes the
-                          sheet that create_report will add.  However, if
-                          in the future we add multiple sheets then this
-                          should still be good.  */
                        undo_report = go_undo_binary_new
                                (wb,
                                 workbook_sheet_state_new (wb),
@@ -1147,6 +1145,7 @@ dialog_init (SolverState *state)
        INIT_BOOL_ENTRY ("non_neg_button", options.assume_non_negative);
        INIT_BOOL_ENTRY ("all_int_button", options.assume_discrete);
        INIT_BOOL_ENTRY ("program", options.program_report);
+       INIT_BOOL_ENTRY ("sensitivity", options.sensitivity_report);
 
        input = gnm_solver_param_get_input (param);
        if (input != NULL)
diff --git a/src/dialogs/solver.ui b/src/dialogs/solver.ui
index 9d5e07a..9f26964 100644
--- a/src/dialogs/solver.ui
+++ b/src/dialogs/solver.ui
@@ -738,7 +738,6 @@
                     <child>
                       <object class="GtkCheckButton" id="program">
                         <property name="label" translatable="yes">P_rogram</property>
-                        <property name="use_action_appearance">False</property>
                         <property name="visible">True</property>
                         <property name="can_focus">True</property>
                         <property name="receives_default">False</property>
@@ -752,6 +751,22 @@
                         <property name="position">0</property>
                       </packing>
                     </child>
+                    <child>
+                      <object class="GtkCheckButton" id="sensitivity">
+                        <property name="label" translatable="yes">_Sensitivity</property>
+                        <property name="visible">True</property>
+                        <property name="can_focus">True</property>
+                        <property name="receives_default">False</property>
+                        <property name="use_underline">True</property>
+                        <property name="xalign">0</property>
+                        <property name="draw_indicator">True</property>
+                      </object>
+                      <packing>
+                        <property name="expand">False</property>
+                        <property name="fill">False</property>
+                        <property name="position">1</property>
+                      </packing>
+                    </child>
                   </object>
                   <packing>
                     <property name="expand">True</property>
diff --git a/src/gnumeric-fwd.h b/src/gnumeric-fwd.h
index bba2876..c4eb4c8 100644
--- a/src/gnumeric-fwd.h
+++ b/src/gnumeric-fwd.h
@@ -7,6 +7,7 @@ G_BEGIN_DECLS
 
 typedef struct GnmComplete_             GnmComplete;
 typedef struct GnmScenario_             GnmScenario;
+typedef struct GnmSolver_              GnmSolver;
 typedef struct GnmSolverConstraint_     GnmSolverConstraint;
 typedef struct GnmSolverFactory_        GnmSolverFactory;
 typedef struct GnmSolverParameters_    GnmSolverParameters;
diff --git a/src/tools/ChangeLog b/src/tools/ChangeLog
index 4dc1b6a..f9ceeac 100644
--- a/src/tools/ChangeLog
+++ b/src/tools/ChangeLog
@@ -1,5 +1,9 @@
 2016-02-10  Morten Welinder  <terra gnome org>
 
+       * gnm-solver.c (gnm_solver_create_program_report): Split from
+       gnm_solver_create_report and fix naming of constraints.
+       (gnm_solver_create_sensitivity_report): New function.
+
        * dao.c (dao_autofit_rows): New function.
        (dao_autofit_these_columns): Simply use colrow_autofit.  Don't
        shrink any columns.
diff --git a/src/tools/gnm-solver.c b/src/tools/gnm-solver.c
index 3726670..af357fe 100644
--- a/src/tools/gnm-solver.c
+++ b/src/tools/gnm-solver.c
@@ -401,6 +401,39 @@ gnm_solver_constraint_as_str (GnmSolverConstraint const *c, Sheet *sheet)
        return g_string_free (buf, FALSE);
 }
 
+char *
+gnm_solver_constraint_part_as_str (GnmSolverConstraint const *c, int i,
+                                  GnmSolverParameters *sp)
+{
+       const char * const type_str[] = {
+               "\xe2\x89\xa4" /* "<=" */,
+               "\xe2\x89\xa5" /* ">=" */,
+               "=",
+               N_("Int"),
+               N_("Bool")
+       };
+       const char *type = type_str[c->type];
+       gboolean translate = (c->type >= GNM_SOLVER_INTEGER);
+       GString *buf;
+       gnm_float cl, cr;
+       GnmCell *lhs, *rhs;
+
+       if (!gnm_solver_constraint_get_part (c, sp, i, &lhs, &cl, &rhs, &cr))
+               return NULL;
+
+       buf = g_string_new (NULL);
+
+       g_string_append (buf, cell_name (lhs));
+       g_string_append_c (buf, ' ');
+       g_string_append (buf, translate ? _(type) : type);
+       if (gnm_solver_constraint_has_rhs (c)) {
+               g_string_append_c (buf, ' ');
+               g_string_append (buf, cell_name (rhs));
+       }
+
+       return g_string_free (buf, FALSE);
+}
+
 /* ------------------------------------------------------------------------- */
 
 enum {
@@ -466,6 +499,7 @@ gnm_solver_param_equal (GnmSolverParameters const *a,
             a->options.assume_discrete != b->options.assume_discrete ||
             a->options.automatic_scaling != b->options.automatic_scaling ||
             a->options.program_report != b->options.program_report ||
+            a->options.sensitivity_report != b->options.sensitivity_report ||
             a->options.add_scenario != b->options.add_scenario ||
            strcmp (a->options.scenario_name, b->options.scenario_name) ||
            a->options.gradient_order != b->options.gradient_order)
@@ -775,6 +809,7 @@ enum {
        SOL_PROP_REASON,
        SOL_PROP_PARAMS,
        SOL_PROP_RESULT,
+       SOL_PROP_SENSITIVITY,
        SOL_PROP_STARTTIME,
        SOL_PROP_ENDTIME,
        SOL_PROP_FLIP_SIGN
@@ -836,6 +871,10 @@ gnm_solver_get_property (GObject *object, guint property_id,
                g_value_set_object (value, sol->result);
                break;
 
+       case SOL_PROP_SENSITIVITY:
+               g_value_set_object (value, sol->sensitivity);
+               break;
+
        case SOL_PROP_STARTTIME:
                g_value_set_double (value, sol->starttime);
                break;
@@ -884,6 +923,13 @@ gnm_solver_set_property (GObject *object, guint property_id,
                break;
        }
 
+       case SOL_PROP_SENSITIVITY: {
+               GnmSolverSensitivity *s = g_value_dup_object (value);
+               if (sol->sensitivity) g_object_unref (sol->sensitivity);
+               sol->sensitivity = s;
+               break;
+       }
+
        case SOL_PROP_STARTTIME:
                sol->starttime = g_value_get_double (value);
                break;
@@ -1488,8 +1534,8 @@ print_vector (const char *name, const gnm_float *v, int n)
        g_printerr ("\n");
 }
 
-void
-gnm_solver_create_report (GnmSolver *solver, const char *name)
+static void
+gnm_solver_create_program_report (GnmSolver *solver, const char *name)
 {
        GnmSolverParameters *params = solver->params;
        int R = 0;
@@ -1616,7 +1662,7 @@ gnm_solver_create_report (GnmSolver *solver, const char *name)
                                                     &rhs, &cr);
                     i++) {
                        gnm_float slack = 0;
-                       char *ctxt = gnm_solver_constraint_as_str (c, params->sheet);
+                       char *ctxt = gnm_solver_constraint_part_as_str (c, i, params);
                        dao_set_cell (dao, 1, R, ctxt);
                        g_free (ctxt);
 
@@ -1667,61 +1713,153 @@ gnm_solver_create_report (GnmSolver *solver, const char *name)
 
        /* ---------------------------------------- */
 
-       if (gnm_solver_has_solution (solver) &&
-           gnm_debug_flag ("solver-sensitivity")) {
+       dao_autofit_columns (dao);
+       dao_redraw_respan (dao);
+
+       dao_free (dao);
+}
+
+static void
+gnm_solver_create_sensitivity_report (GnmSolver *solver, const char *name)
+{
+       GnmSolverParameters *params = solver->params;
+       GnmSolverSensitivity *sols = solver->sensitivity;
+       int R = 0;
+       data_analysis_output_t *dao;
+       GSList *l;
+
+       if (!sols)
+               return;
+
+       dao = dao_init_new_sheet (NULL);
+       dao->sheet = params->sheet;
+       dao_prepare_output (NULL, dao, name);
+
+       /* ---------------------------------------- */
+
+       if (solver->input_cells->len > 0) {
                unsigned ui;
-               const int N = 500;
-               gnm_float const *xs0 = solver->result->solution;
 
-               if (gnm_solver_debug ()) {
-                       gnm_float *g = gnm_solver_compute_gradient (solver, xs0);
-                       print_vector ("Computed gradient", g, solver->input_cells->len);
-                       g_free (g);
-               }
+               ADD_HEADER (_("Variables"));
 
-               gnm_solver_set_vars (solver, xs0);
+               dao_set_cell (dao, 1, R, _("Cell"));
+               dao_set_cell (dao, 2, R, _("Final\nValue"));
+               dao_set_cell (dao, 3, R, _("Reduced\nCost"));
+               dao_set_cell (dao, 4, R, _("Lower\nLimit"));
+               dao_set_cell (dao, 5, R, _("Upper\nLimit"));
+               dao_set_align (dao, 1, R, 5, R, GNM_HALIGN_CENTER, GNM_VALIGN_BOTTOM);
+               dao_autofit_these_rows (dao, R, R);
+               R++;
 
                for (ui = 0; ui < solver->input_cells->len; ui++) {
-                       char *txt;
-                       int i, j;
-                       gnm_float x0, y0, x, y;
                        GnmCell *cell = g_ptr_array_index (solver->input_cells, ui);
+                       gnm_float L = sols->vars[ui].low;
+                       gnm_float H = sols->vars[ui].high;
+                       gnm_float red = sols->vars[ui].reduced_cost;
+                       gnm_float s = solver->result->solution[ui];
+
+                       char *cname = gnm_solver_cell_name (cell, params->sheet);
+                       dao_set_cell (dao, 1, R, cname);
+                       g_free (cname);
+                       dao_set_cell_value (dao, 2, R, value_new_float (s));
+                       add_value_or_special (dao, 3, R, red);
+                       add_value_or_special (dao, 4, R, L);
+                       add_value_or_special (dao, 5, R, H);
 
                        R++;
-                       txt = g_strdup_printf (_("Neighborhood for %s\n"),
-                                              cell_name (cell));
-                       ADD_HEADER (txt);
-                       g_free (txt);
+               }
 
-                       x0 = xs0[ui];
-                       y0 = solver->result->value;
+               R++;
+       }
 
-                       x = x0;
-                       for (i = 0; i < N; i++)
-                               for (j = 0; j < 10; j++)
-                                       x = nextafter (x, gnm_ninf);
+       /* ---------------------------------------- */
 
-                       for (i = -N; i <= +N; i++) {
-                               gnm_solver_set_var (solver, ui, x);
-                               y = gnm_solver_get_target_value (solver);
+       ADD_HEADER (_("Constraints"));
 
-                               add_value_or_special (dao, 1, R, x - x0);
-                               add_value_or_special (dao, 2, R, y - y0);
-                               R++;
+       if (params->constraints) {
+               dao_set_cell (dao, 1, R, _("Constraint"));
+               dao_set_cell (dao, 2, R, _("Shadow\nPrice"));
+               dao_set_cell (dao, 3, R, _("Constraint\nLHS"));
+               dao_set_cell (dao, 4, R, _("Constraint\nRHS"));
+               dao_set_cell (dao, 5, R, _("Lower\nLimit"));
+               dao_set_cell (dao, 6, R, _("Upper\nLimit"));
+               dao_set_align (dao, 1, R, 6, R, GNM_HALIGN_CENTER, GNM_VALIGN_BOTTOM);
+               dao_autofit_these_rows (dao, R, R);
+       } else {
+               dao_set_cell (dao, 1, R, _("No constraints"));
+       }
+       R++;
 
-                               for (j = 0; j < 10; j++)
-                                       x = nextafter (x, gnm_pinf);
+       for (l = params->constraints; l; l = l->next) {
+               GnmSolverConstraint *c = l->data;
+               int i, cidx = 0;
+               gnm_float cl, cr;
+               GnmCell *lhs, *rhs;
+
+               for (i = 0;
+                    gnm_solver_constraint_get_part (c, params, i,
+                                                    &lhs, &cl,
+                                                    &rhs, &cr);
+                    i++, cidx++) {
+                       char *ctxt;
+
+                       switch (c->type) {
+                       case GNM_SOLVER_INTEGER:
+                       case GNM_SOLVER_BOOLEAN:
+                               continue;
+                       default:
+                               ; // Nothing
+                       }
+
+                       ctxt = gnm_solver_constraint_part_as_str (c, i, params);
+                       dao_set_cell (dao, 1, R, ctxt);
+                       g_free (ctxt);
+
+                       if (lhs) {
+                               gnm_cell_eval (lhs);
+                               cl = value_get_as_float (lhs->value);
                        }
-                       gnm_solver_set_var (solver, ui, x0);
+                       if (rhs) {
+                               gnm_cell_eval (rhs);
+                               cr = value_get_as_float (rhs->value);
+                       }
+
+                       add_value_or_special (dao, 2, R, sols->constraints[cidx].shadow_price);
+                       add_value_or_special (dao, 3, R, cl);
+                       add_value_or_special (dao, 4, R, cr);
+                       add_value_or_special (dao, 5, R, sols->constraints[cidx].low);
+                       add_value_or_special (dao, 6, R, sols->constraints[cidx].high);
+
+                       R++;
                }
        }
 
        /* ---------------------------------------- */
+
+       dao_autofit_columns (dao);
        dao_redraw_respan (dao);
 
        dao_free (dao);
 }
 
+void
+gnm_solver_create_report (GnmSolver *solver, const char *base)
+{
+       GnmSolverParameters *params = solver->params;
+
+       if (params->options.program_report) {
+               char *name = g_strdup_printf (base, _("Program"));
+               gnm_solver_create_program_report (solver, name);
+               g_free (name);
+       }
+
+       if (params->options.sensitivity_report) {
+               char *name = g_strdup_printf (base, _("Sensitivity"));
+               gnm_solver_create_sensitivity_report (solver, name);
+               g_free (name);
+       }
+}
+
 #undef AT_LIMIT
 #undef ADD_HEADER
 #undef MARK_BAD
@@ -2211,6 +2349,14 @@ gnm_solver_class_init (GObjectClass *object_class)
                                     GSF_PARAM_STATIC |
                                     G_PARAM_READWRITE));
 
+       g_object_class_install_property (object_class, SOL_PROP_SENSITIVITY,
+               g_param_spec_object ("sensitivity",
+                                    P_("Sensitivity"),
+                                    P_("Sensitivity results"),
+                                    GNM_SOLVER_SENSITIVITY_TYPE,
+                                    GSF_PARAM_STATIC |
+                                    G_PARAM_READWRITE));
+
        g_object_class_install_property (object_class, SOL_PROP_STARTTIME,
                g_param_spec_double ("starttime",
                                     P_("Start Time"),
@@ -2296,6 +2442,133 @@ GSF_CLASS (GnmSolverResult, gnm_solver_result,
 
 /* ------------------------------------------------------------------------- */
 
+static GObjectClass *gnm_solver_sensitivity_parent_class;
+
+enum {
+       SOLS_PROP_0,
+       SOLS_PROP_SOLVER
+};
+
+static void
+gnm_solver_sensitivity_constructed (GObject *obj)
+{
+       GnmSolverSensitivity *sols = GNM_SOLVER_SENSITIVITY (obj);
+       GnmSolver *sol = sols->solver;
+       GnmSolverParameters *sp = sol->params;
+       const int n = sol->input_cells->len;
+       int i, cn;
+       GSList *l;
+
+       /* Chain to parent first */
+       gnm_solver_sensitivity_parent_class->constructed (obj);
+
+       sols->vars = g_new (struct GnmSolverSensitivityVars_, n);
+       for (i = 0; i < n; i++) {
+               sols->vars[i].low = gnm_nan;
+               sols->vars[i].high = gnm_nan;
+               sols->vars[i].reduced_cost = gnm_nan;
+       }
+
+       cn = 0;
+       for (l = sp->constraints; l; l = l->next) {
+               GnmSolverConstraint *c = l->data;
+               int i;
+               gnm_float cl, cr;
+               GnmCell *lhs, *rhs;
+
+               for (i = 0;
+                    gnm_solver_constraint_get_part (c, sp, i,
+                                                    &lhs, &cl,
+                                                    &rhs, &cr);
+                    i++) {
+                       cn++;
+               }
+       }
+       sols->constraints = g_new (struct GnmSolverSensitivityConstraints_, cn);
+       for (i = 0; i < cn; i++) {
+               sols->constraints[i].low = gnm_nan;
+               sols->constraints[i].high = gnm_nan;
+               sols->constraints[i].shadow_price = gnm_nan;
+       }
+}
+
+static void
+gnm_solver_sensitivity_finalize (GObject *obj)
+{
+       GnmSolverSensitivity *r = GNM_SOLVER_SENSITIVITY (obj);
+       g_free (r->vars);
+       g_free (r->constraints);
+       gnm_solver_sensitivity_parent_class->finalize (obj);
+}
+
+static void
+gnm_solver_sensitivity_get_property (GObject *object, guint property_id,
+                                    GValue *value, GParamSpec *pspec)
+{
+       GnmSolverSensitivity *sols = (GnmSolverSensitivity *)object;
+
+       switch (property_id) {
+       case SOLS_PROP_SOLVER:
+               g_value_set_object (value, sols->solver);
+               break;
+
+       default:
+               G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
+               break;
+       }
+}
+
+static void
+gnm_solver_sensitivity_set_property (GObject *object, guint property_id,
+                                    GValue const *value, GParamSpec *pspec)
+{
+       GnmSolverSensitivity *sols = (GnmSolverSensitivity *)object;
+
+       switch (property_id) {
+       case SOLS_PROP_SOLVER:
+               /* We hold no ref.  */
+               sols->solver = g_value_get_object (value);
+               break;
+
+       default:
+               G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
+               break;
+       }
+}
+
+static void
+gnm_solver_sensitivity_class_init (GObjectClass *object_class)
+{
+       gnm_solver_sensitivity_parent_class =
+               g_type_class_peek_parent (object_class);
+
+       object_class->finalize = gnm_solver_sensitivity_finalize;
+       object_class->constructed = gnm_solver_sensitivity_constructed;
+       object_class->set_property = gnm_solver_sensitivity_set_property;
+       object_class->get_property = gnm_solver_sensitivity_get_property;
+
+       g_object_class_install_property
+               (object_class, SOLS_PROP_SOLVER,
+                g_param_spec_object ("solver",
+                                     P_("Solver"),
+                                     P_("Solver"),
+                                     GNM_SOLVER_TYPE,
+                                     GSF_PARAM_STATIC |
+                                     G_PARAM_CONSTRUCT_ONLY |
+                                     G_PARAM_READWRITE));
+}
+
+GSF_CLASS (GnmSolverSensitivity, gnm_solver_sensitivity,
+          gnm_solver_sensitivity_class_init, NULL, G_TYPE_OBJECT)
+
+GnmSolverSensitivity *
+gnm_solver_sensitivity_new (GnmSolver *sol)
+{
+       return g_object_new (GNM_SOLVER_SENSITIVITY_TYPE, "solver", sol, NULL);
+}
+
+/* ------------------------------------------------------------------------- */
+
 static GObjectClass *gnm_sub_solver_parent_class;
 
 enum {
@@ -2351,6 +2624,9 @@ gnm_sub_solver_clear (GnmSubSolver *subsol)
 
        if (subsol->name_from_cell)
                g_hash_table_remove_all (subsol->name_from_cell);
+
+       if (subsol->constraint_from_name)
+               g_hash_table_remove_all (subsol->constraint_from_name);
 }
 
 static void
@@ -2380,6 +2656,9 @@ gnm_sub_solver_finalize (GObject *obj)
        g_hash_table_destroy (subsol->name_from_cell);
        subsol->name_from_cell = NULL;
 
+       g_hash_table_destroy (subsol->constraint_from_name);
+       subsol->constraint_from_name = NULL;
+
        gnm_sub_solver_parent_class->finalize (obj);
 }
 
@@ -2396,6 +2675,10 @@ gnm_sub_solver_init (GnmSubSolver *subsol)
                                       g_free, NULL);
        subsol->name_from_cell =
                g_hash_table_new (g_direct_hash, g_direct_equal);
+
+       subsol->constraint_from_name =
+               g_hash_table_new_full (g_str_hash, g_str_equal,
+                                      g_free, NULL);
 }
 
 static void
@@ -2555,6 +2838,34 @@ gnm_sub_solver_get_cell_name (GnmSubSolver *subsol,
        return g_hash_table_lookup (subsol->name_from_cell, (gpointer)cell);
 }
 
+const char *
+gnm_sub_solver_name_constraint (GnmSubSolver *subsol,
+                               int cidx,
+                               const char *name)
+{
+       char *name_copy = g_strdup (name);
+
+       g_hash_table_insert (subsol->constraint_from_name,
+                            name_copy,
+                            GINT_TO_POINTER (cidx));
+
+       return name_copy;
+}
+
+int
+gnm_sub_solver_find_constraint (GnmSubSolver *subsol, const char *name)
+{
+       gpointer idx;
+
+       if (g_hash_table_lookup_extended (subsol->constraint_from_name,
+                                         (gpointer)name,
+                                         NULL, &idx))
+               return GPOINTER_TO_INT (idx);
+       else
+               return -1;
+}
+
+
 char *
 gnm_sub_solver_locate_binary (const char *binary, const char *solver,
                              const char *url,
diff --git a/src/tools/gnm-solver.h b/src/tools/gnm-solver.h
index 11da950..8697575 100644
--- a/src/tools/gnm-solver.h
+++ b/src/tools/gnm-solver.h
@@ -98,6 +98,8 @@ void gnm_solver_constraint_side_as_str (GnmSolverConstraint const *c,
                                        Sheet const *sheet,
                                        GString *buf, gboolean lhs);
 char *gnm_solver_constraint_as_str (GnmSolverConstraint const *c, Sheet *sheet);
+char *gnm_solver_constraint_part_as_str (GnmSolverConstraint const *c, int i,
+                                        GnmSolverParameters *sp);
 
 /* ------------------------------------------------------------------------- */
 
@@ -110,6 +112,7 @@ typedef struct {
        gboolean            assume_discrete;
        gboolean            automatic_scaling;
        gboolean            program_report;
+       gboolean            sensitivity_report;
        gboolean            add_scenario;
        gchar               *scenario_name;
        unsigned            gradient_order;
@@ -186,13 +189,42 @@ typedef struct {
 GType gnm_solver_result_get_type (void);
 
 /* ------------------------------------------------------------------------- */
+
+#define GNM_SOLVER_SENSITIVITY_TYPE   (gnm_solver_sensitivity_get_type ())
+#define GNM_SOLVER_SENSITIVITY(o)     (G_TYPE_CHECK_INSTANCE_CAST ((o), GNM_SOLVER_SENSITIVITY_TYPE, 
GnmSolverSensitivity))
+
+typedef struct {
+       GObject parent;
+
+       GnmSolver *solver;
+
+       struct GnmSolverSensitivityVars_ {
+               gnm_float low, high;      //  Allowable range
+               gnm_float reduced_cost;
+       } *vars;
+
+       struct GnmSolverSensitivityConstraints_ {
+               gnm_float low, high;      //  Allowable range
+               gnm_float shadow_price;
+       } *constraints;
+} GnmSolverSensitivity;
+
+typedef struct {
+       GObjectClass parent_class;
+} GnmSolverSensitivityClass;
+
+GType gnm_solver_sensitivity_get_type (void);
+
+GnmSolverSensitivity *gnm_solver_sensitivity_new (GnmSolver *sol);
+
+/* ------------------------------------------------------------------------- */
 /* Generic Solver class. */
 
 #define GNM_SOLVER_TYPE        (gnm_solver_get_type ())
 #define GNM_SOLVER(o)          (G_TYPE_CHECK_INSTANCE_CAST ((o), GNM_SOLVER_TYPE, GnmSolver))
 #define GNM_IS_SOLVER(o)       (G_TYPE_CHECK_INSTANCE_TYPE ((o), GNM_SOLVER_TYPE))
 
-typedef struct {
+struct GnmSolver_ {
        GObject parent;
 
        GnmSolverStatus status;
@@ -200,6 +232,7 @@ typedef struct {
 
        GnmSolverParameters *params;
        GnmSolverResult *result;
+       GnmSolverSensitivity *sensitivity;
        double starttime, endtime;
        gboolean flip_sign;
 
@@ -210,7 +243,7 @@ typedef struct {
        gnm_float *min;
        gnm_float *max;
        guint8 *discrete;
-} GnmSolver;
+};
 
 typedef struct {
        GObjectClass parent_class;
@@ -295,6 +328,8 @@ typedef struct {
        GHashTable *cell_from_name;
        GHashTable *name_from_cell;
 
+       GHashTable *constraint_from_name;
+
        GPid child_pid;
        guint child_watch;
 
@@ -332,6 +367,11 @@ GnmCell *gnm_sub_solver_find_cell (GnmSubSolver *subsol, const char *name);
 const char *gnm_sub_solver_get_cell_name (GnmSubSolver *subsol,
                                          GnmCell const *cell);
 
+const char *gnm_sub_solver_name_constraint (GnmSubSolver *subsol,
+                                           int cidx,
+                                           const char *name);
+int gnm_sub_solver_find_constraint (GnmSubSolver *subsol, const char *name);
+
 char *gnm_sub_solver_locate_binary (const char *binary, const char *solver,
                                    const char *url,
                                    WBCGtk *wbcg);
diff --git a/src/xml-sax-read.c b/src/xml-sax-read.c
index fa7b419..7730659 100644
--- a/src/xml-sax-read.c
+++ b/src/xml-sax-read.c
@@ -2660,7 +2660,8 @@ xml_sax_solver_start (GsfXMLIn *xin, xmlChar const **attrs)
                           gnm_xml_attr_bool (attrs, "NonNeg", &(sp->options.assume_non_negative)) ||
                           gnm_xml_attr_bool (attrs, "Discr", &(sp->options.assume_discrete)) ||
                           gnm_xml_attr_bool (attrs, "AutoScale", &(sp->options.automatic_scaling)) ||
-                          gnm_xml_attr_bool (attrs, "ProgramR", &(sp->options.program_report)))
+                          gnm_xml_attr_bool (attrs, "ProgramR", &(sp->options.program_report)) ||
+                          gnm_xml_attr_bool (attrs, "SensitivityR", &(sp->options.sensitivity_report)))
                        ; /* Nothing */
        }
 
diff --git a/src/xml-sax-write.c b/src/xml-sax-write.c
index 13e9cc5..975ab2f 100644
--- a/src/xml-sax-write.c
+++ b/src/xml-sax-write.c
@@ -1087,6 +1087,8 @@ xml_write_solver (GnmOutputXML *state)
                param->options.automatic_scaling);
        gsf_xml_out_add_bool (state->output, "ProgramR",
                param->options.program_report);
+       gsf_xml_out_add_bool (state->output, "SensitivityR",
+               param->options.sensitivity_report);
 
        for (ptr = param->constraints; ptr != NULL ; ptr = ptr->next) {
                GnmSolverConstraint const *c = ptr->data;


[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]