[gnumeric] Solver: move axis-direction polisher generic and improve it.

[Morten Welinder <mortenw src gnome org>]

commit 7daa743d7b3cf997cf171a429884c01896b76efd
Author: Morten Welinder <terra gnome org>
Date:   Tue Apr 28 21:31:00 2015 -0400

    Solver: move axis-direction polisher generic and improve it.

 plugins/nlsolve/gnm-nlsolve.c |   74 ++++------------------------------
 src/tools/ChangeLog           |    5 ++
 src/tools/gnm-solver.c        |   87 +++++++++++++++++++++++++++++++++++++++++
 src/tools/gnm-solver.h        |    3 +
 4 files changed, 104 insertions(+), 65 deletions(-)
---
diff --git a/plugins/nlsolve/gnm-nlsolve.c b/plugins/nlsolve/gnm-nlsolve.c
index 4d7d0b6..022d3bf 100644
--- a/plugins/nlsolve/gnm-nlsolve.c
+++ b/plugins/nlsolve/gnm-nlsolve.c
@@ -525,55 +525,6 @@ rosenbrock_iter (GnmNlsolve *nl)
 }
 
 static gboolean
-polish_iter (GnmNlsolve *nl)
-{
-       GnmSolver *sol = nl->sol;
-       GnmIterSolver *isol = nl->isol;
-       const int n = nl->n;
-       gnm_float *x;
-       gnm_float step;
-       gboolean any_at_all = FALSE;
-
-       x = g_new (gnm_float, n);
-       for (step = gnm_pow2 (-10); step > GNM_EPSILON; step *= 0.75) {
-               int c, s;
-
-               for (c = 0; c < n; c++) {
-                       for (s = 0; s <= 1; s++) {
-                               gnm_float y;
-                               gnm_float dx = step * gnm_abs (isol->xk[c]);
-
-                               if (dx == 0) dx = step;
-                               if (s) dx = -dx;
-
-                               memcpy (x, isol->xk, n * sizeof (gnm_float));
-                               x[c] += dx;
-                               set_vector (nl, x);
-                               y = get_value (nl);
-
-                               if (y < isol->yk && gnm_solver_check_constraints (sol))  {
-                                       isol->yk = y;
-                                       memcpy (isol->xk, x, n * sizeof (gnm_float));
-                                       any_at_all = TRUE;
-                                       if (nl->debug)
-                                               g_printerr ("Polish step %.15" GNM_FORMAT_g
-                                                           " in direction %d\n",
-                                                           dx, c);
-                                       break;
-                               }
-                       }
-               }
-       }
-
-       g_free (x);
-
-       if (any_at_all)
-               set_solution (nl);
-
-       return any_at_all;
-}
-
-static gboolean
 gnm_nlsolve_iterate (GnmSolverIterator *iter, GnmNlsolve *nl)
 {
        GnmIterSolver *isol = nl->isol;
@@ -591,12 +542,6 @@ gnm_nlsolve_iterate (GnmSolverIterator *iter, GnmNlsolve *nl)
        return rosenbrock_iter (nl);
 }
 
-static gboolean
-gnm_nlsolve_polish (GnmSolverIterator *iter, GnmNlsolve *nl)
-{
-       return !nl->tentative && polish_iter (nl);
-}
-
 static void
 gnm_nlsolve_final (GnmNlsolve *nl)
 {
@@ -628,11 +573,12 @@ nlsolve_solver_factory_functional (GnmSolverFactory *factory)
 GnmSolver *
 nlsolve_solver_factory (GnmSolverFactory *factory, GnmSolverParameters *params)
 {
-       GnmIterSolver *res = g_object_new
+       GnmIterSolver *isol = g_object_new
                (GNM_ITER_SOLVER_TYPE,
                 "params", params,
                 "flip-sign", (params->problem_type == GNM_SOLVER_MAXIMIZE),
                 NULL);
+       GnmSolver *sol = GNM_SOLVER (isol);
        GnmNlsolve *nl = g_new0 (GnmNlsolve, 1);
        GnmSolverIteratorCompound *citer;
        GnmSolverIterator *iter;
@@ -642,24 +588,22 @@ nlsolve_solver_factory (GnmSolverFactory *factory, GnmSolverParameters *params)
        iter = gnm_solver_iterator_new_func (G_CALLBACK (gnm_nlsolve_iterate), nl);
        gnm_solver_iterator_compound_add (citer, iter, 1);
 
-       iter = gnm_solver_iterator_new_func (G_CALLBACK (gnm_nlsolve_polish), nl);
-       gnm_solver_iterator_compound_add (citer, iter, 0);
-
-       res->iterator = GNM_SOLVER_ITERATOR (citer);
+       gnm_solver_iterator_compound_add (citer, gnm_solver_iterator_new_polish (isol), 0);
 
-       nl->sol = GNM_SOLVER (res);
-       nl->isol = res;
+       isol->iterator = GNM_SOLVER_ITERATOR (citer);
 
+       nl->sol = sol;
+       nl->isol = isol;
        nl->debug = gnm_solver_debug ();
        nl->min_factor = 1e-10;
        nl->n = nl->sol->input_cells->len;
 
-       g_signal_connect (res, "prepare", G_CALLBACK (gnm_nlsolve_prepare), nl);
+       g_signal_connect (isol, "prepare", G_CALLBACK (gnm_nlsolve_prepare), nl);
 
-       g_object_set_data_full (G_OBJECT (res), PRIVATE_KEY, nl,
+       g_object_set_data_full (G_OBJECT (isol), PRIVATE_KEY, nl,
                                (GDestroyNotify)gnm_nlsolve_final);
 
-       return GNM_SOLVER (res);
+       return sol;
 }
 
 /* ------------------------------------------------------------------------- */
diff --git a/src/tools/ChangeLog b/src/tools/ChangeLog
index dee794e..f4a0047 100644
--- a/src/tools/ChangeLog
+++ b/src/tools/ChangeLog
@@ -1,3 +1,8 @@
+2015-04-28  Morten Welinder  <terra gnome org>
+
+       * gnm-solver.c (gnm_solver_iterator_new_polish): New function with
+       guts from nlsolve.
+
 2015-04-25  Morten Welinder  <terra gnome org>
 
        * gnm-solver.c (gnm_solver_param_get_input_cells): Return result
diff --git a/src/tools/gnm-solver.c b/src/tools/gnm-solver.c
index 70504fc..c0f3f34 100644
--- a/src/tools/gnm-solver.c
+++ b/src/tools/gnm-solver.c
@@ -2227,6 +2227,93 @@ gnm_solver_iterator_new_func (GCallback iterate, gpointer user)
        return iter;
 }
 
+static gboolean
+cb_polish_iter (GnmSolverIterator *iter, GnmIterSolver *isol)
+{
+       GnmSolver *sol = GNM_SOLVER (isol);
+       const int n = sol->input_cells->len;
+       gnm_float *x;
+       gboolean progress = FALSE;
+       gboolean debug = gnm_solver_debug ();
+       int c;
+
+       x = g_new (gnm_float, n);
+
+       for (c = 0; c < n; c++) {
+               gnm_float xc = isol->xk[c], dx;
+               int e;
+               gboolean try_reverse = TRUE;
+               gboolean try_bigger = TRUE;
+
+               (void)gnm_frexp (xc, &e);
+               dx = gnm_ldexp (1, e + 10);
+               if (dx == 0) dx = GNM_MIN;
+
+               memcpy (x, isol->xk, n * sizeof (gnm_float));
+               while (1) {
+                       gnm_float y;
+
+                       x[c] = xc + dx;
+                       if (x[c] == xc)
+                               break;
+                       gnm_solver_set_vars (sol, x);
+                       y = gnm_iter_solver_get_target_value (isol);
+
+                       if (y < isol->yk && gnm_solver_check_constraints (sol))  {
+                               /* Success!  */
+
+                               isol->yk = y;
+                               isol->xk[c] = xc = x[c];
+                               progress = TRUE;
+                               try_reverse = FALSE;
+                               if (debug)
+                                       g_printerr ("Polish step %.15" GNM_FORMAT_g
+                                                   " in direction %d\n",
+                                                   dx, c);
+                               if (try_bigger) {
+                                       dx *= 2;
+                                       if (gnm_finite (dx))
+                                               continue;
+                               }
+                               break;
+                       } else {
+                               /* Step didn't work.  Restore and find new step to try.  */
+                               x[c] = xc;
+
+                               if (try_reverse) {
+                                       dx = -dx;
+                                       if (dx < 0)
+                                               continue;
+                               }
+
+                               dx /= 2;
+                               try_bigger = FALSE;
+                       }
+               }
+       }
+
+       g_free (x);
+
+       if (progress)
+               gnm_iter_solver_set_solution (isol);
+
+       return progress;
+}
+
+/**
+ * gnm_solver_iterator_new_polish:
+ * @isol: the solver to operate on
+ *
+ * Returns: (transfer full): an iterator object that can be used to polish
+ * a solution by simple axis-parallel movement.
+ */
+GnmSolverIterator *
+gnm_solver_iterator_new_polish (GnmIterSolver *isol)
+{
+       return gnm_solver_iterator_new_func (G_CALLBACK (cb_polish_iter), isol);
+}
+
+
 GSF_CLASS (GnmSolverIterator, gnm_solver_iterator,
           gnm_solver_iterator_class_init, NULL, G_TYPE_OBJECT)
 
diff --git a/src/tools/gnm-solver.h b/src/tools/gnm-solver.h
index abdaa52..f457c70 100644
--- a/src/tools/gnm-solver.h
+++ b/src/tools/gnm-solver.h
@@ -334,7 +334,10 @@ typedef struct {
 } GnmSolverIteratorClass;
 
 GType gnm_solver_iterator_get_type (void);
+
 GnmSolverIterator *gnm_solver_iterator_new_func (GCallback iterate, gpointer user);
+GnmSolverIterator *gnm_solver_iterator_new_polish (GnmIterSolver *isol);
+
 gboolean gnm_solver_iterator_iterate (GnmSolverIterator *iter);