[gnumeric] Mathfunc: api cleanup.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] Mathfunc: api cleanup.
- Date: Thu, 8 Dec 2016 23:43:54 +0000 (UTC)
commit 934bfe62b3a50146aea6e8057c33e52ff04ffed9
Author: Morten Welinder <terra gnome org>
Date: Thu Dec 8 18:43:26 2016 -0500
Mathfunc: api cleanup.
ChangeLog | 6 +
plugins/fn-math/functions.c | 2 +-
plugins/nlsolve/gnm-nlsolve.c | 182 ++++++++++++-----------------------
src/mathfunc.c | 214 +++++++++++++++++++++++++++++++++++++++++
src/mathfunc.h | 15 +++
src/regression.h | 4 -
6 files changed, 299 insertions(+), 124 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 9b3991c..3a0cd6d 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,9 @@
+2016-12-08 Morten Welinder <terra gnome org>
+
+ * src/mathfunc.c (gnm_linear_solve): Use proper matrix type. All
+ callers changed.
+ (gnm_linear_solve_multiple): Ditto.
+
2016-10-02 Morten Welinder <terra gnome org>
* src/libgnumeric.c (gnm_pre_parse_init): Don't pretend the
diff --git a/plugins/fn-math/functions.c b/plugins/fn-math/functions.c
index 0607bcb..7800812 100644
--- a/plugins/fn-math/functions.c
+++ b/plugins/fn-math/functions.c
@@ -3065,7 +3065,7 @@ gnumeric_linsolve (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
goto out;
}
- regres = gnm_linear_solve_multiple (A->data, B->data, A->rows, B->cols);
+ regres = gnm_linear_solve_multiple (A, B);
if (regres != GO_REG_ok && regres != GO_REG_near_singular_good) {
res = value_new_error_NUM (ei->pos);
diff --git a/plugins/nlsolve/gnm-nlsolve.c b/plugins/nlsolve/gnm-nlsolve.c
index 15d2232..794b343 100644
--- a/plugins/nlsolve/gnm-nlsolve.c
+++ b/plugins/nlsolve/gnm-nlsolve.c
@@ -7,6 +7,7 @@
#include "regression.h"
#include "rangefunc.h"
#include "workbook.h"
+#include "mathfunc.h"
#include <gsf/gsf-impl-utils.h>
#include <glib/gi18n-lib.h>
#include <string.h>
@@ -53,6 +54,9 @@ typedef struct {
int tentative;
gnm_float *tentative_xk, tentative_yk;
+ // Newton state
+ // (nothing right now)
+
/* Parameters: */
gboolean debug;
gnm_float min_factor;
@@ -120,11 +124,13 @@ print_vector (const char *name, const gnm_float *v, int n)
g_printerr ("\n");
}
+#if 0
static void
set_value (GnmNlsolve *nl, int i, gnm_float x)
{
gnm_solver_set_var (nl->sol, i, x);
}
+#endif
static void
set_vector (GnmNlsolve *nl, const gnm_float *xs)
@@ -185,126 +191,84 @@ compute_gradient (GnmNlsolve *nl, const gnm_float *xs)
return gnm_solver_compute_gradient (nl->sol, xs);
}
-static gnm_float **
-compute_hessian (GnmNlsolve *nl, const gnm_float *xs, const gnm_float *g0)
-{
- gnm_float **H, *xs2;
- const int n = nl->n;
- int i, j;
-
- H = g_new (gnm_float *, n);
-
- xs2 = g_memdup (xs, n * sizeof (gnm_float));
- for (i = 0; i < n; i++) {
- gnm_float x0 = xs[i];
- gnm_float dx;
- const gnm_float *g;
- gnm_float eps = gnm_pow2 (-25);
-
- if (x0 == 0)
- dx = eps;
- else
- dx = gnm_abs (x0) * eps;
-
- xs2[i] = x0 + dx;
- g = H[i] = compute_gradient (nl, xs2);
- xs2[i] = x0;
-
- 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;
-
- set_value (nl, i, x0);
- }
-
- g_free (xs2);
- return H;
-}
-
static gboolean
-newton_improve (GnmNlsolve *nl, gnm_float *xs, gnm_float *y, gnm_float ymax)
+newton_improve (GnmNlsolve *nl, gnm_float *xs)
{
GnmSolver *sol = nl->sol;
+ GnmIterSolver *isol = nl->isol;
const int n = nl->n;
- gnm_float *g, **H, *d;
+ int i;
+ gnm_float *g, *d, *xs2;
+ GnmMatrix *H;
gboolean ok;
+ xs2 = g_new (gnm_float, n);
g = compute_gradient (nl, xs);
- H = compute_hessian (nl, xs, g);
+ H = gnm_solver_compute_hessian (sol, xs);
d = g_new (gnm_float, n);
- ok = (gnm_linear_solve (H, g, n, d) == 0);
-
+ ok = (gnm_linear_solve_posdef (H, g, d) == GO_REG_ok);
if (ok) {
+ for (i = 0; i < n; i++)
+ d[i] = 0 - d[i];
+ }
+
+ if (nl->debug) {
int i;
- gnm_float y2, *xs2 = g_new (gnm_float, n);
- gnm_float best_f = 42;
- // We try these step sizes. We really should not need
- // negative, but if H isn't positive definite it might
- // work.
- static const gnm_float fs[] = {
- 1.0, 0.5, 1.0 / 16,
- -1.0, -1.0 / 16,
- };
- unsigned ui;
-
- if (nl->debug) {
- int i;
- for (i = 0; i < n; i++)
- print_vector (NULL, H[i], n);
+ g_printerr ("Hessian:\n");
+ for (i = 0; i < n; i++)
+ print_vector (NULL, H->data[i], n);
+ print_vector ("g", g, n);
+ if (ok)
print_vector ("d", d, n);
- print_vector ("g", g, n);
- }
+ else
+ g_printerr ("Failed to solve Newton step.\n");
+ }
- ok = FALSE;
- for (ui = 0 ; ui < G_N_ELEMENTS (fs); ui++) {
- gnm_float f = fs[ui];
- int i;
- for (i = 0; i < n; i++)
- xs2[i] = xs[i] - f * d[i];
- set_vector (nl, xs2);
- y2 = get_value (nl);
- if (nl->debug) {
- print_vector ("xs2", xs2, n);
- g_printerr ("Obj value %.15" GNM_FORMAT_g "\n",
- y2);
- }
+ if (ok) {
+ gnm_float y2, f;
- if (y2 < ymax && gnm_solver_check_constraints (sol)) {
- best_f = f;
- ymax = y2;
- break;
- }
- }
+ for (i = 0; i < n; i++)
+ xs2[i] = xs[i] + d[i];
+ set_vector (nl, xs2);
+ y2 = get_value (nl);
- if (best_f != 42) {
- for (i = 0; i < n; i++)
- xs[i] = xs[i] - best_f * d[i];
- *y = ymax;
+ ok = FALSE;
+ if (y2 < isol->yk && gnm_solver_check_constraints (sol)) {
+ if (nl->debug)
+ g_printerr ("Accepting newton step\n");
+ memcpy (isol->xk, xs2, n * sizeof (gnm_float));
+ isol->yk = y2;
+ set_solution (nl);
ok = TRUE;
- }
+ } else {
+ if (nl->debug)
+ g_printerr ("Full newton step would go to %g\n", y2);
- g_free (xs2);
- } else {
- if (nl->debug)
- g_printerr ("Failed to solve Newton step.\n");
+ f = gnm_solver_line_search
+ (sol, xs, d, FALSE, 0.75, 1, 0.01, &y2);
+
+ if (f > 0 && f < 1 && y2 < isol->yk) {
+ if (nl->debug)
+ g_printerr ("Accepting reduced newton step with f=%g\n", f);
+ for (i = 0; i < n; i++)
+ isol->xk[i] = xs[i] + f * d[i];
+ isol->yk = y2;
+ set_solution (nl);
+ ok = TRUE;
+ }
+ }
}
g_free (d);
g_free (g);
- free_matrix (H, n);
+ gnm_matrix_free (H);
+ g_free (xs2);
return ok;
}
static void
-rosenbrock_init (GnmNlsolve *nl)
+nlsolve_init (GnmNlsolve *nl)
{
const int n = nl->n;
int i, j;
@@ -366,9 +330,9 @@ rosenbrock_iter (GnmNlsolve *nl)
}
}
- if (isol->iterations % 100 == 0 &&
- gnm_solver_has_analytic_gradient (sol)) {
- if (newton_improve (nl, isol->xk, &isol->yk, isol->yk))
+ if ((isol->iterations < 20 || isol->iterations % 100 == 0) &&
+ gnm_solver_has_analytic_hessian (sol)) {
+ if (newton_improve (nl, isol->xk))
return TRUE;
}
@@ -512,26 +476,6 @@ rosenbrock_iter (GnmNlsolve *nl)
} else {
nl->smallsteps++;
}
-
- if (0 && !nl->tentative && nl->smallsteps > 50) {
- gnm_float yk = isol->yk;
-
- nl->tentative = 10;
- nl->tentative_xk = g_memdup (isol->xk, n * sizeof (gnm_float));
- nl->tentative_yk = yk;
-
- for (i = 0; i < 4; i++) {
- gnm_float ymax = yk +
- gnm_abs (yk) * (0.10 / (i + 1));
- if (i > 0)
- ymax = MIN (ymax, isol->yk);
- if (!newton_improve (nl, isol->xk, &isol->yk, ymax))
- break;
- }
-
- if (nl->debug)
- print_vector ("Tentative move to", isol->xk, n);
- }
}
g_free (x);
@@ -552,7 +496,7 @@ gnm_nlsolve_iterate (GnmSolverIterator *iter, GnmNlsolve *nl)
const int n = nl->n;
if (isol->iterations == 0)
- rosenbrock_init (nl);
+ nlsolve_init (nl);
if (nl->debug) {
g_printerr ("Iteration %ld at %.15" GNM_FORMAT_g "\n",
diff --git a/src/mathfunc.c b/src/mathfunc.c
index 0960245..10b4147 100644
--- a/src/mathfunc.c
+++ b/src/mathfunc.c
@@ -48,6 +48,7 @@
#include <goffice/goffice.h>
#include <glib/gstdio.h>
#include <value.h>
+#include <gutils.h>
/* R code wants this, so provide it. */
#ifndef IEEE_754
@@ -2493,6 +2494,10 @@ gnm_float qbinom(gnm_float p, gnm_float n, gnm_float pr, gboolean lower_tail, gb
}
}
+#if 0
+}
+#endif
+
/* ------------------------------------------------------------------------ */
/* Imported src/nmath/dnbinom.c from R. */
/*
@@ -2704,6 +2709,9 @@ gnm_float qnbinom(gnm_float p, gnm_float n, gnm_float pr, gboolean lower_tail, g
}
}
+#if 0
+}
+#endif
/* ------------------------------------------------------------------------ */
/* Imported src/nmath/dbeta.c from R. */
/*
@@ -2875,6 +2883,7 @@ gnm_float dhyper(gnm_float x, gnm_float r, gnm_float b, gnm_float n, gboolean gi
return( (give_log) ? p1 + p2 - p3 : p1*p2/p3 );
}
+
/* ------------------------------------------------------------------------ */
/* Imported src/nmath/phyper.c from R. */
/*
@@ -5336,6 +5345,211 @@ gnm_matrix_eigen (GnmMatrix const *m, GnmMatrix *EIG, gnm_float *eigenvalues)
/* ------------------------------------------------------------------------- */
+static void
+swap_row_and_col (GnmMatrix *M, int a, int b)
+{
+ gnm_float *r;
+ int i;
+
+ // Swap rows
+ r = M->data[a];
+ M->data[a] = M->data[b];
+ M->data[b] = r;
+
+ // Swap cols
+ for (i = 0; i < M->rows; i++) {
+ gnm_float d = M->data[i][a];
+ M->data[i][a] = M->data[i][b];
+ M->data[i][b] = d;
+ }
+}
+
+
+gboolean
+gnm_matrix_modified_cholesky (GnmMatrix const *A,
+ GnmMatrix *L,
+ gnm_float *D,
+ gnm_float *E,
+ int *P)
+{
+ int n = A->cols;
+ gnm_float nu, xsi, gam, bsqr, delta;
+ int i, j;
+ GnmMatrix *G, *C;
+
+ g_return_val_if_fail (A->rows == A->cols, FALSE);
+ g_return_val_if_fail (A->rows == L->rows, FALSE);
+ g_return_val_if_fail (A->cols == L->cols, FALSE);
+
+ // Copy A into L; Use G and C as aliases for L.
+ for (i = 0; i < n; i++)
+ for (j = 0; j < n; j++)
+ L->data[i][j] = A->data[i][j];
+ G = L;
+ C = G;
+
+ // Init permutation as identity.
+ for (i = 0; i < n; i++)
+ P[i] = i;
+
+ nu = n == 1 ? 1.0 : gnm_sqrt (n * n - 1);
+ gam = xsi = 0;
+ for (i = 0; i < n; i++) {
+ gnm_float aii = gnm_abs (G->data[i][i]);
+ gam = MAX (gam, aii);
+ for (j = i + 1; j < n; j++) {
+ gnm_float aij = gnm_abs (G->data[i][j]);
+ xsi = MAX (xsi, aij);
+ }
+ }
+ bsqr = MAX (MAX (gam, xsi / nu), GNM_EPSILON);
+ delta = MAX (gam + xsi, 1.0) * GNM_EPSILON;
+
+ for (j = 0; j < n; j++) {
+ int q, s;
+ gnm_float theta_j = 0, dj;
+
+ q = j;
+ for (i = j + 1; i < n; i++) {
+ if (gnm_abs (C->data[i][i]) > gnm_abs (C->data[q][q]))
+ q = i;
+ }
+
+ if (j != q) {
+ int a;
+ gnm_float b;
+
+ swap_row_and_col (L, j, q);
+ a = P[j]; P[j] = P[q]; P[q] = a;
+ b = D[j]; D[j] = D[q]; D[q] = b;
+ if (E) { b = E[j]; E[j] = E[q]; E[q] = b; }
+ }
+
+ for (s = 0; s < j; s++)
+ L->data[j][s] = C->data[j][s] / D[s];
+
+ for (i = j + 1; i < n; i++) {
+ int s;
+ gnm_float d = G->data[i][j];
+
+ for (s = 0; s < j; s++)
+ d -= L->data[j][s] * C->data[i][s];
+ C->data[i][j] = d;
+
+ theta_j = MAX (theta_j, gnm_abs (d));
+ }
+
+ dj = MAX (theta_j * theta_j / bsqr, delta);
+ dj = MAX (dj, gnm_abs (C->data[j][j]));
+ D[j] = dj;
+ if (E) E[j] = dj - C->data[j][j];
+
+ for (i = j + 1; i < n; i++) {
+ gnm_float cij = C->data[i][j];
+ C->data[i][i] -= cij * cij / D[j];
+ }
+ }
+
+ for (i = 0; i < n; i++) {
+ for (j = i + 1; j < n; j++)
+ L->data[i][j] = 0;
+ L->data[i][i] = 1;
+ }
+
+ return TRUE;
+}
+
+GORegressionResult
+gnm_linear_solve_posdef (GnmMatrix const *A, const gnm_float *b,
+ gnm_float *x)
+{
+ int i, j, n;
+ GnmMatrix *L;
+ gnm_float *D, *E;
+ int *P;
+ GORegressionResult res;
+ gboolean ok;
+
+ g_return_val_if_fail (A != NULL, GO_REG_invalid_dimensions);
+ g_return_val_if_fail (A->rows == A->cols, GO_REG_invalid_dimensions);
+ g_return_val_if_fail (b != NULL, GO_REG_invalid_dimensions);
+ g_return_val_if_fail (x != NULL, GO_REG_invalid_dimensions);
+
+ n = A->cols;
+ L = gnm_matrix_new (n, n);
+ D = g_new (gnm_float, n);
+ E = g_new (gnm_float, n);
+ P = g_new (int, n);
+
+ ok = gnm_matrix_modified_cholesky (A, L, D, E, P);
+ if (!ok) {
+ res = GO_REG_invalid_data;
+ goto done;
+ }
+
+ if (gnm_debug_flag ("posdef")) {
+ for (i = 0; i < n; i++)
+ g_printerr ("Posdef E[i] = %g\n", E[P[i]]);
+ }
+
+ // The only information from the above we use is E and P.
+ // However, we reuse the memory for L for A+E
+ for (i = 0; i < n; i++) {
+ for (j = 0; j < n; j++)
+ L->data[i][j] = A->data[i][j];
+ L->data[i][i] += E[P[i]];
+ }
+
+ res = gnm_linear_solve (L, b, x);
+
+done:
+ g_free (P);
+ g_free (E);
+ g_free (D);
+ gnm_matrix_free (L);
+
+ return res;
+}
+
+/* ------------------------------------------------------------------------- */
+
+GORegressionResult
+gnm_linear_solve (GnmMatrix const *A, const gnm_float *b,
+ gnm_float *x)
+{
+ g_return_val_if_fail (A != NULL, GO_REG_invalid_dimensions);
+ g_return_val_if_fail (A->rows == A->cols, GO_REG_invalid_dimensions);
+ g_return_val_if_fail (b != NULL, GO_REG_invalid_dimensions);
+ g_return_val_if_fail (x != NULL, GO_REG_invalid_dimensions);
+
+ return
+#ifdef GNM_WITH_LONG_DOUBLE
+ go_linear_solvel
+#else
+ go_linear_solve
+#endif
+ (A->data, b, A->rows, x);
+}
+
+GORegressionResult
+gnm_linear_solve_multiple (GnmMatrix const *A, GnmMatrix *B)
+{
+ g_return_val_if_fail (A != NULL, GO_REG_invalid_dimensions);
+ g_return_val_if_fail (B != NULL, GO_REG_invalid_dimensions);
+ g_return_val_if_fail (A->rows == A->cols, GO_REG_invalid_dimensions);
+ g_return_val_if_fail (A->rows == B->rows, GO_REG_invalid_dimensions);
+
+ return
+#ifdef GNM_WITH_LONG_DOUBLE
+ go_linear_solve_multiplel
+#else
+ go_linear_solve_multiple
+#endif
+ (A->data, B->data, A->rows, B->cols);
+}
+
+/* ------------------------------------------------------------------------- */
+
#ifdef GNM_SUPPLIES_ERFL
long double
erfl (long double x)
diff --git a/src/mathfunc.h b/src/mathfunc.h
index c6b7fb2..401b080 100644
--- a/src/mathfunc.h
+++ b/src/mathfunc.h
@@ -139,6 +139,21 @@ gboolean gnm_matrix_is_empty (GnmMatrix const *m);
void gnm_matrix_multiply (GnmMatrix *C, const GnmMatrix *A, const GnmMatrix *B);
gboolean gnm_matrix_eigen (GnmMatrix const *m, GnmMatrix *EIG, gnm_float *eigenvalues);
+
+gboolean gnm_matrix_modified_cholesky (GnmMatrix const *A,
+ GnmMatrix *L,
+ gnm_float *D,
+ gnm_float *E,
+ int *P);
+
+GORegressionResult gnm_linear_solve_posdef (GnmMatrix const *A, const gnm_float *b,
+ gnm_float *x);
+
+GORegressionResult gnm_linear_solve (GnmMatrix const *A, const gnm_float *b,
+ gnm_float *x);
+
+GORegressionResult gnm_linear_solve_multiple (GnmMatrix const *A, GnmMatrix *B);
+
/* ------------------------------------------------------------------------- */
void mathfunc_init (void);
diff --git a/src/regression.h b/src/regression.h
index 5aa8e4d..977a28e 100644
--- a/src/regression.h
+++ b/src/regression.h
@@ -21,8 +21,6 @@ G_BEGIN_DECLS
# define gnm_matrix_invert go_matrix_invertl
# define gnm_matrix_pseudo_inverse go_matrix_pseudo_inversel
# define gnm_matrix_determinant go_matrix_determinantl
-# define gnm_linear_solve go_linear_solvel
-# define gnm_linear_solve_multiple go_linear_solve_multiplel
#else
# define gnm_regression_stat_t go_regression_stat_t
# define gnm_regression_stat_new go_regression_stat_new
@@ -37,8 +35,6 @@ G_BEGIN_DECLS
# define gnm_matrix_invert go_matrix_invert
# define gnm_matrix_pseudo_inverse go_matrix_pseudo_inverse
# define gnm_matrix_determinant go_matrix_determinant
-# define gnm_linear_solve go_linear_solve
-# define gnm_linear_solve_multiple go_linear_solve_multiple
#endif
G_END_DECLS
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
