GNOME.org
 

[gnumeric] solver: initial cut at non-linear solver.

commit d9174845e9728aaa7b833bea6db451d779e1bfdd
Author: Morten Welinder <terra gnome org>
Date:   Sat May 22 21:19:11 2010 -0400

    solver: initial cut at non-linear solver.

 NEWS                          |    1 +
 configure.in                  |    3 +-
 plugins/Makefile.am           |    2 +-
 plugins/nlsolve/.gitignore    |    9 +
 plugins/nlsolve/ChangeLog     |    3 +
 plugins/nlsolve/Makefile.am   |   20 ++
 plugins/nlsolve/boot.c        |   24 ++
 plugins/nlsolve/boot.h        |   12 +
 plugins/nlsolve/gnm-nlsolve.c |  482 +++++++++++++++++++++++++++++++++++++++++
 plugins/nlsolve/plugin.xml.in |   19 ++
 src/regression.h              |    2 +
 11 files changed, 575 insertions(+), 2 deletions(-)
---
diff --git a/NEWS b/NEWS
index 655b45c..6099ee8 100644
--- a/NEWS
+++ b/NEWS
@@ -3,6 +3,7 @@ Gnumeric 1.10.5
 Morten:
 	* Fix stf crash.  [#619283]
 	* Persist solver model type.
+	* Add non-linear solver.
 
 --------------------------------------------------------------------------
 Gnumeric 1.10.4
diff --git a/configure.in b/configure.in
index 60248fc..367d3ae 100644
--- a/configure.in
+++ b/configure.in
@@ -146,7 +146,7 @@ PKG_PROG_PKG_CONFIG(0.18)
 
 dnl *****************************
 libspreadsheet_reqs="
-	libgoffice-${GOFFICE_API_VER}	>= 0.8.3
+	libgoffice-${GOFFICE_API_VER}	>= 0.8.5
 	libgsf-1		>= 1.14.15
 	libxml-2.0		>= 2.4.12
 "
@@ -1088,6 +1088,7 @@ plugins/gnome-glossary/Makefile
 plugins/html/Makefile
 plugins/lotus-123/Makefile
 plugins/lpsolve/Makefile
+plugins/nlsolve/Makefile
 plugins/glpk/Makefile
 plugins/mps/Makefile
 plugins/oleo/Makefile
diff --git a/plugins/Makefile.am b/plugins/Makefile.am
index cb20731..370871f 100644
--- a/plugins/Makefile.am
+++ b/plugins/Makefile.am
@@ -1,5 +1,5 @@
 SUBDIRS_FILE_FORMATS = excel lotus-123 oleo sc sylk xbase html dif qpro \
-	plan-perfect applix openoffice lpsolve glpk
+	plan-perfect applix openoffice lpsolve nlsolve glpk
 
 if ENABLE_SOLVER
   SUBDIRS_FILE_FORMATS += mps
diff --git a/plugins/nlsolve/.gitignore b/plugins/nlsolve/.gitignore
new file mode 100644
index 0000000..6429675
--- /dev/null
+++ b/plugins/nlsolve/.gitignore
@@ -0,0 +1,9 @@
+Makefile
+Makefile.in
+.deps
+.libs
+*.lo
+*.la
+plugin.xml
+*.loT
+*.sw*
diff --git a/plugins/nlsolve/ChangeLog b/plugins/nlsolve/ChangeLog
new file mode 100644
index 0000000..45e20de
--- /dev/null
+++ b/plugins/nlsolve/ChangeLog
@@ -0,0 +1,3 @@
+2010-05-21  Morten Welinder  <terra gnome org>
+
+	* Initial Implementation
diff --git a/plugins/nlsolve/Makefile.am b/plugins/nlsolve/Makefile.am
new file mode 100644
index 0000000..66e2521
--- /dev/null
+++ b/plugins/nlsolve/Makefile.am
@@ -0,0 +1,20 @@
+AM_CPPFLAGS = \
+    -DGNOMELOCALEDIR=\""$(datadir)/locale"\" 	\
+    -I$(top_srcdir)/src	-I$(top_builddir)/src	\
+    $(GNUMERIC_CFLAGS)
+
+gnumeric_plugin_nlsolvedir = $(gnumeric_plugindir)/nlsolve
+xmldir = $(gnumeric_plugin_nlsolvedir)
+gnumeric_plugin_nlsolve_LTLIBRARIES = nlsolve.la
+nlsolve_la_LDFLAGS = -module $(GNUMERIC_PLUGIN_LDFLAGS)
+nlsolve_la_SOURCES = \
+	boot.h boot.c \
+	gnm-nlsolve.c
+
+xml_in_files = plugin.xml.in
+xml_DATA = $(xml_in_files:.xml.in=.xml)
+
+ INTLTOOL_XML_RULE@
+
+EXTRA_DIST = ChangeLog $(xml_in_files)
+DISTCLEANFILES = $(xml_DATA)
diff --git a/plugins/nlsolve/boot.c b/plugins/nlsolve/boot.c
new file mode 100644
index 0000000..35f61ba
--- /dev/null
+++ b/plugins/nlsolve/boot.c
@@ -0,0 +1,24 @@
+/*
+ * Copyright (C) 2010 Morten Welinder (terra gnome org)
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ */
+
+#include <gnumeric-config.h>
+#include <gnumeric.h>
+#include "boot.h"
+#include <gnm-plugin.h>
+
+GNM_PLUGIN_MODULE_HEADER;
diff --git a/plugins/nlsolve/boot.h b/plugins/nlsolve/boot.h
new file mode 100644
index 0000000..392bb53
--- /dev/null
+++ b/plugins/nlsolve/boot.h
@@ -0,0 +1,12 @@
+#ifndef GNUMERIC_NLSOLVE_BOOT_H
+#define GNUMERIC_NLSOLVE_BOOT_H
+
+#include "gnumeric.h"
+#include <goffice/goffice.h>
+#include <gsf/gsf-output.h>
+
+void
+nlsolve_file_save (GOFileSaver const *fs, GOIOContext *io_context,
+		   WorkbookView const *wb_view, GsfOutput *output);
+
+#endif
diff --git a/plugins/nlsolve/gnm-nlsolve.c b/plugins/nlsolve/gnm-nlsolve.c
new file mode 100644
index 0000000..8ef492a
--- /dev/null
+++ b/plugins/nlsolve/gnm-nlsolve.c
@@ -0,0 +1,482 @@
+#include <gnumeric-config.h>
+#include "gnumeric.h"
+#include <tools/gnm-solver.h>
+#include "cell.h"
+#include "sheet.h"
+#include "value.h"
+#include "ranges.h"
+#include "regression.h"
+#include <gsf/gsf-impl-utils.h>
+#include <glib/gi18n-lib.h>
+#include <string.h>
+
+#define PRIVATE_KEY "::nlsolve::"
+
+typedef struct {
+	GnmSolver *parent;
+	GnmSheetRange srinput;
+
+	/* Input cells.  */
+	GPtrArray *vars;
+
+	/* Target cell. */
+	GnmCell *target;
+
+	gboolean debug;
+
+	int iterations;
+
+	/* Next axis direction to try.  */
+	int dim;
+
+	/* Parameters: */
+	gnm_float eps;
+	int max_iter;
+	gnm_float min_factor;
+
+	guint idle_tag;
+} GnmNlsolve;
+
+static void
+gnm_nlsolve_cleanup (GnmNlsolve *nl)
+{
+	if (nl->idle_tag) {
+		g_source_remove (nl->idle_tag);
+		nl->idle_tag = 0;
+	}
+
+	if (nl->vars) {
+		g_ptr_array_free (nl->vars, TRUE);
+		nl->vars = NULL;
+	}
+}
+
+static void
+gnm_nlsolve_final (GnmNlsolve *nl)
+{
+	gnm_nlsolve_cleanup (nl);
+	g_free (nl);
+}
+
+static gboolean
+check_program (const GnmSolverParameters *params, GError **err)
+{
+	GSList *l;
+
+	if (params->options.assume_discrete)
+		goto no_discrete;
+
+	for (l = params->constraints; l; l = l->next) {
+		GnmSolverConstraint *c  = l->data;
+		if (c->type == GNM_SOLVER_INTEGER ||
+		    c->type == GNM_SOLVER_BOOLEAN)
+			goto no_discrete;
+	}
+
+	if (params->problem_type != GNM_SOLVER_MINIMIZE) {
+		/* Not a fundamental problem, just not done.  */
+		g_set_error (err,
+			     go_error_invalid (),
+			     0,
+			     _("This solver handles only minimization."));
+		return FALSE;
+	}
+
+	return TRUE;
+
+no_discrete:
+	g_set_error (err,
+		     go_error_invalid (),
+		     0,
+		     _("This solver does not handle discrete variables."));
+	return FALSE;
+}
+
+static void
+set_value (GnmNlsolve *nl, int i, gnm_float x)
+{
+	GnmCell *cell = g_ptr_array_index (nl->vars, i);
+	gnm_cell_set_value (cell, value_new_float (x));
+	cell_queue_recalc (cell);
+}
+
+static gnm_float
+get_value (GnmNlsolve *nl)
+{
+	GnmValue const *v;
+
+	gnm_cell_eval (nl->target);
+	v = nl->target->value;
+
+	if (VALUE_IS_NUMBER (v) || VALUE_IS_EMPTY (v))
+		return value_get_as_float (v);
+	else
+		return gnm_nan;
+}
+
+static void
+gnm_nlsolve_set_solution (GnmNlsolve *nl)
+{
+	GnmSolver *sol = nl->parent;
+	GnmSolverResult *result = g_object_new (GNM_SOLVER_RESULT_TYPE, NULL);
+
+	result->solution = gnm_solver_get_current_values (sol);
+	result->quality = GNM_SOLVER_RESULT_FEASIBLE;
+	result->value = get_value (nl);
+
+	g_object_set (sol, "result", result, NULL);
+	g_object_unref (result);
+}
+
+static gboolean
+gnm_nlsolve_get_initial_solution (GnmNlsolve *nl, GError **err)
+{
+	GnmSolver *sol = nl->parent;
+
+	if (gnm_solver_check_constraints (sol))
+		goto got_it;
+
+	/* More? */
+
+	g_set_error (err,
+		     go_error_invalid (),
+		     0,
+		     _("The initial values do not satisfy the constraints."));
+	return FALSE;
+
+got_it:
+	gnm_nlsolve_set_solution (nl);
+	return TRUE;
+}
+
+static gboolean
+gnm_nlsolve_prepare (GnmSolver *sol, WorkbookControl *wbc, GError **err,
+		     GnmNlsolve *nl)
+{
+	gboolean ok;
+
+	g_return_val_if_fail (sol->status == GNM_SOLVER_STATUS_READY, FALSE);
+
+	gnm_solver_set_status (sol, GNM_SOLVER_STATUS_PREPARING);
+
+	ok = check_program (sol->params, err);
+	if (ok)
+		ok = gnm_nlsolve_get_initial_solution (nl, err);
+
+	if (ok) {
+		if (nl->debug)
+			g_printerr ("Initial value:\n%15.8" GNM_FORMAT_f "\n",
+				    sol->result->value);
+		gnm_solver_set_status (sol, GNM_SOLVER_STATUS_PREPARED);
+	} else {
+		gnm_nlsolve_cleanup (nl);
+		gnm_solver_set_status (sol, GNM_SOLVER_STATUS_ERROR);
+	}
+
+	return ok;
+}
+
+static gnm_float *
+compute_gradient (GnmNlsolve *nl)
+{
+	gnm_float *g;
+	gnm_float y0;
+	const int n = nl->vars->len;
+	int i;
+
+	y0 = get_value (nl);
+	
+	g = g_new (gnm_float, n);
+	for (i = 0; i < n; i++) {
+		GnmCell *cell = g_ptr_array_index (nl->vars, i);
+		GnmValue *old = value_dup (cell->value);
+		gnm_float x0 = value_get_as_float (old);
+		gnm_float dx;
+		gnm_float y1;
+
+		if (x0 == 0)
+			dx = nl->eps;
+		else
+			dx = gnm_abs (x0) * nl->eps;
+
+		set_value (nl, i, x0 + dx);
+#if 0
+		g_printerr ("  g[%d]: trying %.10" GNM_FORMAT_g "\n",
+			    i, x0 + nl->eps);
+#endif
+
+		y1 = get_value (nl);
+#if 0
+		g_printerr ("  y0=%g   y1=%g\n", y0, y1);
+#endif
+		g[i] = (y1 - y0) / dx;
+
+		gnm_cell_set_value (cell, old);
+		cell_queue_recalc (cell);
+	}
+
+	return g;
+}
+
+static gnm_float **
+compute_hessian (GnmNlsolve *nl, const gnm_float *g0)
+{
+	gnm_float **H;
+	const int n = nl->vars->len;
+	int i, j;
+
+	H = g_new (gnm_float *, n);
+
+	for (i = 0; i < n; i++) {
+		GnmCell *cell = g_ptr_array_index (nl->vars, i);
+		GnmValue *old = value_dup (cell->value);
+		gnm_float x0 = value_get_as_float (old);
+		gnm_float dx;
+		const gnm_float *g;
+
+		if (x0 == 0)
+			dx = nl->eps;
+		else
+			dx = gnm_abs (x0) * nl->eps;
+
+		set_value (nl, i, x0 + dx);
+
+		g = H[i] = compute_gradient (nl);
+		if (nl->debug) {
+			int j;
+
+			g_printerr ("  Gradient %d ", i);
+			for (j = 0; j < n; j++)
+				g_printerr ("%15.8" GNM_FORMAT_f " ", g[j]);
+			g_printerr ("\n");
+		}
+		for (j = 0; j < n; j++)
+			H[i][j] = (g[j] - g0[j]) / dx;
+
+		gnm_cell_set_value (cell, old);
+		cell_queue_recalc (cell);
+	}
+
+	return H;
+}
+
+static void
+print_vector (const char *name, const gnm_float *v, int n)
+{
+	int i;
+
+	if (name)
+		g_printerr ("%s:\n", name);
+	for (i = 0; i < n; i++)
+		g_printerr ("%15.8" GNM_FORMAT_f " ", v[i]);
+	g_printerr ("\n");
+}
+
+static gboolean
+try_direction (GnmNlsolve *nl, const gnm_float *x0, const gnm_float *d)
+{
+	GnmSolver *sol = nl->parent;
+	const int n = nl->vars->len;
+	gnm_float factor;
+
+	for (factor = 1; factor >= nl->min_factor; factor /= 2) {
+		gnm_float y;
+		int i;
+
+		for (i = 0; i < n; i++) {
+			gnm_float x = x0[i] - factor * d[i];
+			set_value (nl, i, x);
+		}
+
+		y = get_value (nl);
+		if (nl->debug)
+			g_printerr ("Factor=%-10" GNM_FORMAT_g
+				    "  value=%.8" GNM_FORMAT_f "\n",
+				    factor, y);
+
+		if (y < sol->result->value) {
+			if (gnm_solver_check_constraints (sol)) {
+				gnm_nlsolve_set_solution (nl);
+				return TRUE;
+			}
+			if (nl->debug)
+				g_printerr ("Would-be solution does not satisfy constraints.\n");
+		}
+	}
+
+	return FALSE;
+}
+
+static gint
+gnm_nlsolve_idle (gpointer data)
+{
+	GnmNlsolve *nl = data;
+	GnmSolver *sol = nl->parent;
+	gnm_float **H, *g, *d, *x0;
+	int i;
+	const int n = nl->vars->len;
+	gboolean ok;
+	gnm_float y0;
+	gboolean call_again = TRUE;
+
+	nl->iterations++;
+
+	y0 = get_value (nl);
+
+	x0 = g_new (gnm_float, n);
+	for (i = 0; i < n; i++) {
+		GnmCell *cell = g_ptr_array_index (nl->vars, i);
+		x0[i] = value_get_as_float (cell->value);
+	}
+	g = compute_gradient (nl);
+	H = compute_hessian (nl, g);
+
+	if (nl->debug) {
+		int i;
+
+		print_vector ("x0", x0, n);
+		print_vector ("Gradient", g, n);
+
+		g_printerr ("Hessian:\n");
+		for (i = 0; i < n; i++) {
+			print_vector (NULL, H[i], n);
+		}
+	}
+
+	d = g_new (gnm_float, n);
+	ok = (gnm_linear_solve (H, g, n, d) == 0);
+	if (!ok)
+		goto stop;
+
+	if (nl->debug)
+		print_vector ("Delta", d, n);
+
+	ok = try_direction (nl, x0, d);
+
+	if (!ok) {
+		int i, j;
+		gnm_float *d1 = g_new (gnm_float, n);
+
+		/*
+		 * The Newton step failed.  This is likely because the
+		 * surface is seriously non-linear.  Try one axis
+		 * direction at a time.
+		 */
+
+		for (j = 0; !ok && j < n; j++) {
+			int c = nl->dim++ % n;
+			if (g[c] == 0)
+				continue;
+
+			for (i = 0; i < n; i++)
+				d1[i] = (i == c)
+					? y0 / g[c]
+					: 0;
+
+			ok = try_direction (nl, x0, d1);
+		}
+
+		g_free (d1);
+	}
+
+	if (!ok) {
+		gnm_solver_set_status (sol, GNM_SOLVER_STATUS_DONE);
+		call_again = FALSE;
+	}
+
+	if (call_again && nl->iterations >= nl->max_iter) {
+		gnm_solver_set_status (sol, GNM_SOLVER_STATUS_DONE);
+		call_again = FALSE;
+	}
+
+out:
+	g_free (d);
+	g_free (x0);
+	g_free (g);
+	for (i = 0; i < n; i++)
+		g_free (H[i]);
+	g_free (H);
+
+	return call_again;
+
+stop:
+	gnm_solver_stop (sol, NULL);
+	call_again = FALSE;
+	goto out;
+}
+
+static gboolean
+gnm_nlsolve_start (GnmSolver *sol, WorkbookControl *wbc, GError **err,
+		   GnmNlsolve *nl)
+{
+	gboolean ok = TRUE;
+
+	g_return_val_if_fail (sol->status == GNM_SOLVER_STATUS_PREPARED, FALSE);
+
+	nl->idle_tag = g_idle_add (gnm_nlsolve_idle, nl);
+	gnm_solver_set_status (sol, GNM_SOLVER_STATUS_RUNNING);
+
+	return ok;
+}
+
+static gboolean
+gnm_nlsolve_stop (GnmSolver *sol, GError *err, GnmNlsolve *nl)
+{
+	g_return_val_if_fail (sol->status == GNM_SOLVER_STATUS_RUNNING, FALSE);
+
+	gnm_nlsolve_cleanup (nl);
+
+	gnm_solver_set_status (sol, GNM_SOLVER_STATUS_CANCELLED);
+
+	return TRUE;
+}
+
+gboolean
+nlsolve_solver_factory_functional (GnmSolverFactory *factory);
+
+gboolean
+nlsolve_solver_factory_functional (GnmSolverFactory *factory)
+{
+	return TRUE;
+}
+
+
+GnmSolver *
+nlsolve_solver_factory (GnmSolverFactory *factory, GnmSolverParameters *params);
+
+GnmSolver *
+nlsolve_solver_factory (GnmSolverFactory *factory, GnmSolverParameters *params)
+{
+	GnmSolver *res = g_object_new (GNM_SOLVER_TYPE,
+				       "params", params,
+				       NULL);
+	GnmNlsolve *nl = g_new0 (GnmNlsolve, 1);
+	GSList *input_cells, *l;
+
+	nl->parent = GNM_SOLVER (res);
+	gnm_sheet_range_from_value (&nl->srinput,
+				    gnm_solver_param_get_input (params));
+	if (nl->srinput.sheet) nl->srinput.sheet = params->sheet;
+
+	nl->debug = gnm_solver_debug ();
+	nl->eps = gnm_pow2 (-25);
+	nl->max_iter = 100;
+	nl->min_factor = 1e-10;
+
+	nl->target = gnm_solver_param_get_target_cell (params);
+
+	nl->vars = g_ptr_array_new ();
+	input_cells = gnm_solver_param_get_input_cells (params);
+	for (l = input_cells; l; l = l->next)
+		g_ptr_array_add (nl->vars, l->data);
+	g_slist_free (input_cells);
+
+	g_signal_connect (res, "prepare", G_CALLBACK (gnm_nlsolve_prepare), nl);
+	g_signal_connect (res, "start", G_CALLBACK (gnm_nlsolve_start), nl);
+	g_signal_connect (res, "stop", G_CALLBACK (gnm_nlsolve_stop), nl);
+
+	g_object_set_data_full (G_OBJECT (res), PRIVATE_KEY, nl,
+				(GDestroyNotify)gnm_nlsolve_final);
+
+	return res;
+}
diff --git a/plugins/nlsolve/plugin.xml.in b/plugins/nlsolve/plugin.xml.in
new file mode 100644
index 0000000..65283d1
--- /dev/null
+++ b/plugins/nlsolve/plugin.xml.in
@@ -0,0 +1,19 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<plugin id="Gnumeric_nlsolve">
+	<information>
+		<_name>Non-Linear Program Solver</_name>
+		<_description>Non-Linear Program Solver</_description>
+	</information>
+	<loader type="Gnumeric_Builtin:module">
+		<attribute name="module_file" value="nlsolve"/>
+	</loader>
+	<services>
+		<service type="solver"
+			 id="nlsolve"
+			 model_type="nlp">
+			<information>
+			  <_description>Nlsolve</_description>
+			</information>
+		</service>
+	</services>
+</plugin>
diff --git a/src/regression.h b/src/regression.h
index 660a8fd..1627b2c 100644
--- a/src/regression.h
+++ b/src/regression.h
@@ -19,6 +19,7 @@ G_BEGIN_DECLS
 #	define gnm_non_linear_regression go_non_linear_regressionl
 #	define gnm_matrix_invert go_matrix_invertl
 #	define gnm_matrix_determinant go_matrix_determinantl
+#	define gnm_linear_solve go_linear_solvel
 #else
 #	define gnm_regression_stat_t go_regression_stat_t
 #	define gnm_regression_stat_new go_regression_stat_new
@@ -31,6 +32,7 @@ G_BEGIN_DECLS
 #	define gnm_non_linear_regression go_non_linear_regression
 #	define gnm_matrix_invert go_matrix_invert
 #	define gnm_matrix_determinant go_matrix_determinant
+#	define gnm_linear_solve go_linear_solve
 #endif
 
 G_END_DECLS
[Date Prev][Date Next]   [Thread Prev][Thread Next]   [Thread Index] [Date Index] [Author Index]