[gnumeric] Solver: add analytic Hessian computation.

[Morten Welinder <mortenw src gnome org>]

commit c263307c22a955aa11fc3d43a8039bcb55195c40
Author: Morten Welinder <terra gnome org>
Date:   Sun Oct 2 15:23:32 2016 -0400

    Solver: add analytic Hessian computation.

 src/tools/ChangeLog    |    5 ++
 src/tools/gnm-solver.c |  106 ++++++++++++++++++++++++++++++++++++++++++++++++
 src/tools/gnm-solver.h |    7 +++
 3 files changed, 118 insertions(+), 0 deletions(-)
---
diff --git a/src/tools/ChangeLog b/src/tools/ChangeLog
index 3a9b616..2871e39 100644
--- a/src/tools/ChangeLog
+++ b/src/tools/ChangeLog
@@ -1,3 +1,8 @@
+2016-10-02  Morten Welinder  <terra gnome org>
+
+       * gnm-solver.c (gnm_solver_compute_hessian): New function to
+       compute analytic hessian of object function.
+
 2016-08-20  Morten Welinder <terra gnome org>
 
        * Release 1.12.32
diff --git a/src/tools/gnm-solver.c b/src/tools/gnm-solver.c
index c34b8a0..c162b88 100644
--- a/src/tools/gnm-solver.c
+++ b/src/tools/gnm-solver.c
@@ -857,6 +857,12 @@ gnm_solver_dispose (GObject *obj)
                sol->gradient = NULL;
        }
 
+       if (sol->hessian) {
+               sol->hessian_status = 0;
+               g_ptr_array_unref (sol->hessian);
+               sol->hessian = NULL;
+       }
+
        gnm_solver_parent_class->dispose (obj);
 }
 
@@ -1950,6 +1956,12 @@ gnm_solver_restore_vars (GnmSolver *sol, GPtrArray *vals)
        g_ptr_array_free (vals, TRUE);
 }
 
+/**
+ * gnm_solver_has_analytic_gradient:
+ * @sol: the solver
+ *
+ * Returns: %TRUE if the gradient can be computed analytically.
+ */
 gboolean
 gnm_solver_has_analytic_gradient (GnmSolver *sol)
 {
@@ -1993,6 +2005,8 @@ gnm_solver_compute_gradient_analytically (GnmSolver *sol, gnm_float const *xs)
                GnmValue *v = gnm_expr_top_eval
                        (te, &ep, GNM_EXPR_EVAL_SCALAR_NON_EMPTY);
                g[i] = VALUE_IS_NUMBER (v) ? value_get_as_float (v) : gnm_nan;
+               if (sol->flip_sign)
+                       g[i] = 0 - g[i];
                value_release (v);
        }
 
@@ -2066,6 +2080,98 @@ gnm_solver_compute_gradient (GnmSolver *sol, gnm_float const *xs)
        return g;
 }
 
+/**
+ * gnm_solver_has_analytic_hessian:
+ * @sol: the solver
+ *
+ * Returns: %TRUE if the Hessian can be computed analytically.
+ */
+gboolean
+gnm_solver_has_analytic_hessian (GnmSolver *sol)
+{
+       const int n = sol->input_cells->len;
+       int i, j;
+       GnmEvalPos ep;
+       GnmExprDeriv *info;
+
+       if (!gnm_solver_has_analytic_gradient (sol))
+               sol->hessian_status = sol->gradient_status;
+
+       if (sol->hessian_status)
+               return sol->hessian_status == 1;
+
+       sol->hessian_status++;
+       sol->hessian = g_ptr_array_new_with_free_func ((GDestroyNotify)gnm_expr_top_unref);
+
+       eval_pos_init_cell (&ep, sol->target);
+       info = gnm_expr_deriv_info_new ();
+       for (i = 0; i < n && sol->hessian_status == 1; i++) {
+               GnmExprTop const *gi = g_ptr_array_index (sol->gradient, i);
+               for (j = i; j < n; j++) {
+                       GnmCell *cell;
+                       GnmExprTop const *te;
+                       GnmEvalPos var;
+
+                       cell = g_ptr_array_index (sol->input_cells, j);
+                       eval_pos_init_cell (&var, cell);
+                       gnm_expr_deriv_info_set_var (info, &var);
+                       te = gnm_expr_top_deriv (gi, &ep, info);
+
+                       if (te)
+                               g_ptr_array_add (sol->hessian, (gpointer)te);
+                       else {
+                               if (gnm_solver_debug ())
+                                       g_printerr ("Unable to compute analytic hessian\n");
+                               sol->hessian_status++;
+                               break;
+                       }
+               }
+       }
+
+       gnm_expr_deriv_info_free (info);
+
+       return sol->hessian_status == 1;
+}
+
+/**
+ * gnm_solver_compute_hessian:
+ * @sol: Solver
+ * @xs: Point to compute Hessian at
+ *
+ * Returns: (transfer full): A vector containing the Hessian.  This
+ * function takes the flip-sign property into account.  The result vector
+ * will be n+(n-1)+...+2+1 elements long containing the triangular
+ * Hessian.  Use symmetry to obtain the full Hessian.
+ */
+gnm_float *
+gnm_solver_compute_hessian (GnmSolver *sol, gnm_float const *xs)
+{
+       int i, hlen;
+       gnm_float *H;
+       GnmEvalPos ep;
+
+       if (!gnm_solver_has_analytic_hessian (sol))
+               return NULL;
+
+       gnm_solver_set_vars (sol, xs);
+
+       hlen = sol->hessian->len;
+       H = g_new (gnm_float, hlen);
+       eval_pos_init_cell (&ep, sol->target);
+       for (i = 0; i < hlen; i++) {
+               GnmExprTop const *te = g_ptr_array_index (sol->hessian, i);
+               GnmValue *v = gnm_expr_top_eval
+                       (te, &ep, GNM_EXPR_EVAL_SCALAR_NON_EMPTY);
+               H[i] = VALUE_IS_NUMBER (v) ? value_get_as_float (v) : gnm_nan;
+               if (sol->flip_sign)
+                       H[i] = 0 - H[i];
+               value_release (v);
+       }
+
+       return H;
+}
+
+
 static gnm_float
 try_step (GnmSolver *sol, gnm_float const *x0, gnm_float const *dir, gnm_float step)
 {
diff --git a/src/tools/gnm-solver.h b/src/tools/gnm-solver.h
index db7fb22..63eea58 100644
--- a/src/tools/gnm-solver.h
+++ b/src/tools/gnm-solver.h
@@ -247,6 +247,10 @@ struct GnmSolver_ {
        // Analytic gradient
        int gradient_status; // 0: not tried; 1: ok; 2: fail
        GPtrArray *gradient;
+
+       // Analytic Hessian
+       int hessian_status; // 0: not tried; 1: ok; 2: fail
+       GPtrArray *hessian;
 };
 
 typedef struct {
@@ -304,6 +308,9 @@ void gnm_solver_restore_vars (GnmSolver *sol, GPtrArray *vals);
 gboolean gnm_solver_has_analytic_gradient (GnmSolver *sol);
 gnm_float *gnm_solver_compute_gradient (GnmSolver *sol, gnm_float const *xs);
 
+gboolean gnm_solver_has_analytic_hessian (GnmSolver *sol);
+gnm_float *gnm_solver_compute_hessian (GnmSolver *sol, gnm_float const *xs);
+
 gnm_float gnm_solver_line_search (GnmSolver *sol,
                                  gnm_float const *x0, gnm_float const *dir,
                                  gboolean try_reverse,