[goffice] regression: use insane precision.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [goffice] regression: use insane precision.
- Date: Mon, 17 May 2010 00:28:22 +0000 (UTC)
commit 3a989cd8f5171dffb58cb74037a963f2af9c2bf8
Author: Morten Welinder <terra gnome org>
Date: Sun May 16 20:28:06 2010 -0400
regression: use insane precision.
goffice/math/go-quad.c | 21 ++--
goffice/math/go-quad.h | 4 +-
goffice/math/go-regression.c | 208 +++++++++++++++++++++++++++---------------
3 files changed, 146 insertions(+), 87 deletions(-)
---
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
index eeb0aef..d81f4e9 100644
--- a/goffice/math/go-quad.c
+++ b/goffice/math/go-quad.c
@@ -189,15 +189,16 @@ SUFFIX(go_quad_div) (QUAD *res, const QUAD *a, const QUAD *b)
res->l = c.h - res->h + c.l;
}
-int
-SUFFIX(go_quad_sgn) (const QUAD *a)
+void
+SUFFIX(go_quad_sqrt) (QUAD *res, const QUAD *a)
{
- DOUBLE d = (a->h == 0 ? a->l : a->h);
-
- if (d == 0)
- return 0;
- else if (d > 0)
- return +1;
- else
- return -1;
+ if (a->h > 0) {
+ QUAD c, u;
+ c.h = SUFFIX(sqrt) (a->h);
+ SUFFIX(go_quad_mul12) (&u, c.h, c.h);
+ c.l = (a->h - u.h - u.l + a->l) * 0.5 / c.h;
+ res->h = c.h + c.l;
+ res->l = c.h - res->h + c.l;
+ } else
+ res->h = res->l = 0;
}
diff --git a/goffice/math/go-quad.h b/goffice/math/go-quad.h
index 53acf31..4dc6763 100644
--- a/goffice/math/go-quad.h
+++ b/goffice/math/go-quad.h
@@ -21,7 +21,7 @@ void go_quad_add (GOQuad *res, const GOQuad *a, const GOQuad *b);
void go_quad_sub (GOQuad *res, const GOQuad *a, const GOQuad *b);
void go_quad_mul (GOQuad *res, const GOQuad *a, const GOQuad *b);
void go_quad_div (GOQuad *res, const GOQuad *a, const GOQuad *b);
-int go_quad_sgn (const GOQuad *a);
+void go_quad_sqrt (GOQuad *res, const GOQuad *a);
void go_quad_mul12 (GOQuad *res, double x, double y);
@@ -42,7 +42,7 @@ void go_quad_addl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
void go_quad_subl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
void go_quad_mull (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
void go_quad_divl (GOQuadl *res, const GOQuadl *a, const GOQuadl *b);
-int go_quad_sgnl (const GOQuadl *a);
+void go_quad_sqrtl (GOQuadl *res, const GOQuadl *a);
void go_quad_mul12l (GOQuadl *res, long double x, long double y);
#endif
diff --git a/goffice/math/go-regression.c b/goffice/math/go-regression.c
index 45464f8..d46c253 100644
--- a/goffice/math/go-regression.c
+++ b/goffice/math/go-regression.c
@@ -71,6 +71,7 @@ general_linear_regressionl (long double *const *const xss, int xdim,
#define MATRIX DOUBLE **
#define CONSTMATRIX DOUBLE *const *const
#define QUAD SUFFIX(GOQuad)
+#define QMATRIX QUAD **
#undef DEBUG_NEAR_SINGULAR
#undef DEBUG_REFINEMENT
@@ -193,49 +194,58 @@ copy_stats (go_regression_stat_t *s1,
/* ------------------------------------------------------------------------- */
-static DOUBLE
-SUFFIX(dot_product) (const DOUBLE *x, const DOUBLE *y, int n)
+static void
+SUFFIX(dot_product) (const QUAD *x, const QUAD *y, int n, QUAD *dp)
{
- DOUBLE d = 0;
int i;
- for (i = 0; i < n; i++)
- d += x[i] * y[i];
- return d;
+ SUFFIX(go_quad_init) (dp, 0, 0);
+ for (i = 0; i < n; i++) {
+ QUAD d;
+ SUFFIX(go_quad_mul) (&d, &x[i], &y[i]);
+ SUFFIX(go_quad_add) (dp, dp, &d);
+ }
}
static GORegressionResult
-SUFFIX(QR) (CONSTMATRIX A, MATRIX Q, MATRIX R, int m, int n)
+SUFFIX(QR) (CONSTMATRIX A, QMATRIX Q, QMATRIX R, int m, int n)
{
- int k;
+ int i, j, k;
- COPY_MATRIX (Q, A, m, n);
+ for (i = 0; i < m; i++)
+ for (j = 0; j < n; j++)
+ SUFFIX(go_quad_init) (&Q[i][j], A[i][j], 0);
for (k = 0; k < m; k++) {
- DOUBLE L;
+ QUAD L;
int i;
- (void)SUFFIX(go_range_sumsq) (Q[k], n, &L);
- L = SUFFIX(sqrt) (L);
+ SUFFIX(dot_product) (Q[k], Q[k], n, &L);
+ SUFFIX(go_quad_sqrt) (&L, &L);
#if 0
PRINT_MATRIX (Q, m, n);
g_printerr ("L[%d] = %20.15" FORMAT_g "\n", k, L);
#endif
- if (L == 0)
+ if (SUFFIX(go_quad_value)(&L) == 0)
return GO_REG_singular;
for (i = 0; i < n; i++)
- Q[k][i] /= L;
+ SUFFIX(go_quad_div) (&Q[k][i], &Q[k][i], &L);
R[k][k] = L;
for (i = k + 1; i < m; i++) {
- DOUBLE d = SUFFIX(dot_product) (Q[k], Q[i], n);
+ QUAD ip;
int j;
- R[k][i] = 0;
- R[i][k] = d;
+ SUFFIX(go_quad_init) (&R[k][i], 0, 0);
- for (j = 0; j < n; j++)
- Q[i][j] -= d * Q[k][j];
+ SUFFIX(dot_product) (Q[k], Q[i], n, &ip);
+ R[i][k] = ip;
+
+ for (j = 0; j < n; j++) {
+ QUAD p;
+ SUFFIX(go_quad_mul) (&p, &ip, &Q[k][j]);
+ SUFFIX(go_quad_sub) (&Q[i][j], &Q[i][j], &p);
+ }
}
}
@@ -246,7 +256,7 @@ SUFFIX(QR) (CONSTMATRIX A, MATRIX Q, MATRIX R, int m, int n)
static void
SUFFIX(calc_residual) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
- const DOUBLE *res, QUAD *residual, QUAD *N2)
+ const QUAD *y, QUAD *residual, QUAD *N2)
{
int i, j;
@@ -258,8 +268,9 @@ SUFFIX(calc_residual) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
SUFFIX(go_quad_init) (&d, b[i], 0);
for (j = 0; j < dim; j++) {
- QUAD e;
- SUFFIX(go_quad_mul12) (&e, A[j][i], res[j]);
+ QUAD Aji, e;
+ SUFFIX(go_quad_init) (&Aji, A[j][i], 0);
+ SUFFIX(go_quad_mul) (&e, &Aji, &y[j]);
SUFFIX(go_quad_sub) (&d, &d, &e);
}
@@ -273,17 +284,14 @@ SUFFIX(calc_residual) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
static void
SUFFIX(refine) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
- CONSTMATRIX Q, CONSTMATRIX R, DOUBLE *result)
+ QMATRIX Q, QMATRIX R, QUAD *result)
{
QUAD *residual = g_new (QUAD, n);
QUAD *newresidual = g_new (QUAD, n);
- DOUBLE *delta = g_new (DOUBLE, dim);
- DOUBLE *newresult = g_new (DOUBLE, dim);
+ QUAD *delta = g_new (QUAD, dim);
+ QUAD *newresult = g_new (QUAD, dim);
int pass;
QUAD best_N2;
- void *state;
-
- state = SUFFIX(go_quad_start) ();
SUFFIX(calc_residual) (A, b, dim, n, result, residual, &best_N2);
#ifdef DEBUG_REFINEMENT
@@ -297,32 +305,34 @@ SUFFIX(refine) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
/* newresult = R^-1 Q^T residual */
for (i = dim - 1; i >= 0; i--) {
- QUAD acc, d, Rii;
+ QUAD acc;
SUFFIX(go_quad_init) (&acc, 0, 0);
for (j = 0; j < n; j++) {
- QUAD q, qr;
- SUFFIX(go_quad_init) (&q, Q[i][j], 0);
- SUFFIX(go_quad_mul) (&qr, &q, residual + j);
+ QUAD qr;
+ SUFFIX(go_quad_mul) (&qr, &Q[i][j], &residual[j]);
SUFFIX(go_quad_add) (&acc, &acc, &qr);
}
for (j = i + 1; j < dim; j++) {
QUAD Rd;
- SUFFIX(go_quad_mul12) (&Rd, R[j][i], delta[j]);
+ SUFFIX(go_quad_mul) (&Rd, &R[j][i], &delta[j]);
SUFFIX(go_quad_sub) (&acc, &acc, &Rd);
}
- SUFFIX(go_quad_init) (&Rii, R[i][i], 0);
- SUFFIX(go_quad_div) (&d, &acc, &Rii);
+ SUFFIX(go_quad_div) (&delta[i], &acc, &R[i][i]);
- delta[i] = SUFFIX(go_quad_value) (&d);
- newresult[i] = delta[i] + result[i];
+ SUFFIX(go_quad_add) (&newresult[i],
+ &delta[i],
+ &result[i]);
#ifdef DEBUG_REFINEMENT
g_printerr ("d[%2d] = %20.15" FORMAT_f
" %20.15" FORMAT_f
" %20.15" FORMAT_f "\n",
- i, delta[i], result[i], newresult[i]);
+ i,
+ SUFFIX(go_quad_value) (&delta[i]),
+ SUFFIX(go_quad_value) (&result[i]),
+ SUFFIX(go_quad_value) (&newresult[i]));
#endif
}
@@ -336,7 +346,7 @@ SUFFIX(refine) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
SUFFIX(go_quad_value) (&N2),
SUFFIX(go_quad_value) (&dres));
#endif
- if (SUFFIX(go_quad_sgn) (&dres) >= 0)
+ if (SUFFIX(go_quad_value) (&dres) >= 0)
break;
best_N2 = N2;
@@ -344,8 +354,6 @@ SUFFIX(refine) (CONSTMATRIX A, const DOUBLE *b, int dim, int n,
COPY_VECTOR (residual, newresidual, n);
}
- SUFFIX(go_quad_end) (state);
-
g_free (residual);
g_free (newresidual);
g_free (newresult);
@@ -575,13 +583,13 @@ SUFFIX(linear_solve) (CONSTMATRIX A, const DOUBLE *b, int n, DOUBLE *res)
gboolean
SUFFIX(go_matrix_invert) (MATRIX A, int n)
{
- MATRIX Q;
- MATRIX R;
+ QMATRIX Q;
+ QMATRIX R;
GORegressionResult err;
gboolean has_result;
- ALLOC_MATRIX (Q, n, n);
- ALLOC_MATRIX (R, n, n);
+ ALLOC_MATRIX2 (Q, n, n, QUAD);
+ ALLOC_MATRIX2 (R, n, n, QUAD);
err = SUFFIX(QR) (A, Q, R, n, n);
has_result = (err == GO_REG_ok || err == GO_REG_near_singular_good);
@@ -591,10 +599,15 @@ SUFFIX(go_matrix_invert) (MATRIX A, int n)
for (k = 0; k < n; k++) {
/* Solve R x = Q^T e_k */
for (i = n - 1; i >= 0; i--) {
- DOUBLE d = Q[i][k];
- for (j = i + 1; j < n; j++)
- d -= R[j][i] * A[k][j];
- A[k][i] = d / R[i][i];
+ QUAD d = Q[i][k];
+ for (j = i + 1; j < n; j++) {
+ QUAD p;
+ SUFFIX(go_quad_init) (&p, A[k][j], 0);
+ SUFFIX(go_quad_mul) (&p, &R[j][i], &p);
+ SUFFIX(go_quad_sub) (&d, &d, &p);
+ }
+ SUFFIX(go_quad_div) (&d, &d, &R[i][i]);
+ A[k][i] = SUFFIX(go_quad_value) (&d);
}
}
}
@@ -674,10 +687,12 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xss, int xdim,
g_free (ys2);
FREE_MATRIX (xss2, xdim, n);
#else
- MATRIX Q;
- MATRIX R;
+ QMATRIX Q;
+ QMATRIX R;
+ QUAD *qresult;
int i,j;
gboolean has_result;
+ void *state;
ZERO_VECTOR (result, xdim);
@@ -687,8 +702,11 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xss, int xdim,
if (xdim > n)
return GO_REG_not_enough_data;
- ALLOC_MATRIX (Q, xdim, n);
- ALLOC_MATRIX (R, xdim, xdim);
+ state = SUFFIX(go_quad_start) ();
+
+ ALLOC_MATRIX2 (Q, xdim, n, QUAD);
+ ALLOC_MATRIX2 (R, xdim, xdim, QUAD);
+ qresult = g_new0 (QUAD, xdim);
regerr = SUFFIX(QR) (xss, Q, R, xdim, n);
has_result = regerr == GO_REG_ok || regerr == GO_REG_near_singular_good;
@@ -696,19 +714,38 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xss, int xdim,
if (has_result) {
/* Solve R res = Q^T ys */
for (i = xdim - 1; i >= 0; i--) {
- DOUBLE d = SUFFIX(dot_product) (Q[i], ys, n);
- for (j = i + 1; j < xdim; j++)
- d -= R[j][i] * result[j];
- result[i] = d / R[i][i];
+ QUAD acc;
+
+ SUFFIX(go_quad_init) (&acc, 0, 0);
+ for (j = 0; j < n; j++) {
+ QUAD p;
+ SUFFIX(go_quad_init) (&p, ys[j], 0);
+ SUFFIX(go_quad_mul) (&p, &p, &Q[i][j]);
+ SUFFIX(go_quad_add) (&acc, &acc, &p);
+ }
+
+ for (j = i + 1; j < xdim; j++) {
+ QUAD d;
+ SUFFIX (go_quad_mul) (&d, &R[j][i], &qresult[j]);
+ SUFFIX (go_quad_sub) (&acc, &acc, &d);
+ }
+
+ SUFFIX(go_quad_div) (&qresult[i], &acc, &R[i][i]);
}
if (xdim > 2)
- SUFFIX(refine) (xss, ys, xdim, n, Q, R, result);
+ SUFFIX(refine) (xss, ys, xdim, n, Q, R, qresult);
+
+ /* Round to plain precision. */
+ for (i = 0; i < xdim; i++) {
+ result[i] = SUFFIX(go_quad_value) (&qresult[i]);
+ SUFFIX(go_quad_init) (&qresult[i], result[i], 0);
+ }
}
if (stat_ && has_result) {
GORegressionResult err2;
- DOUBLE *inv = g_new (DOUBLE, xdim);
+ QUAD *inv = g_new (QUAD, xdim);
int err;
int k;
QUAD N2;
@@ -731,7 +768,7 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xss, int xdim,
g_assert (err == 0);
}
- SUFFIX(calc_residual) (xss, ys, xdim, n, result, NULL, &N2);
+ SUFFIX(calc_residual) (xss, ys, xdim, n, qresult, NULL, &N2);
stat_->ss_resid = SUFFIX(go_quad_value) (&N2);
stat_->sqr_r = (stat_->ss_total == 0)
@@ -740,39 +777,57 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xss, int xdim,
/* FIXME: we want to guard against division by zero. */
stat_->adj_sqr_r = 1 - stat_->ss_resid * (n - 1) /
((n - xdim) * stat_->ss_total);
- stat_->var = (n == xdim)
- ? 0
- : stat_->ss_resid / (n - xdim);
+ if (n == xdim)
+ SUFFIX(go_quad_init) (&N2, 0, 0);
+ else {
+ QUAD d;
+ SUFFIX(go_quad_init) (&d, n - xdim, 0);
+ SUFFIX(go_quad_div) (&N2, &N2, &d);
+ }
+ stat_->var = SUFFIX(go_quad_value) (&N2);
err2 = GO_REG_ok;
stat_->se = g_new (DOUBLE, xdim);
for (k = 0; k < xdim; k++) {
/* Solve R^T z = e_k */
for (i = 0; i < xdim; i++) {
- DOUBLE d = (i == k ? 1 : 0);
- for (j = 0; j < i; j++)
- d -= R[i][j] * inv[j];
- inv[i] = d / R[i][i];
+ QUAD d;
+ SUFFIX(go_quad_init) (&d, i == k ? 1 : 0, 0);
+ for (j = 0; j < i; j++) {
+ QUAD p;
+ SUFFIX(go_quad_mul) (&p, &R[i][j], &inv[j]);
+ SUFFIX(go_quad_sub) (&d, &d, &p);
+ }
+ SUFFIX(go_quad_div) (&inv[i], &d, &R[i][i]);
}
/* Solve R inv = z */
for (i = xdim - 1; i >= 0; i--) {
- DOUBLE d = inv[i];
- for (j = i + 1; j < xdim; j++)
- d -= R[j][i] * inv[j];
- inv[i] = d / R[i][i];
+ QUAD d = inv[i];
+ for (j = i + 1; j < xdim; j++) {
+ QUAD p;
+ SUFFIX(go_quad_mul) (&p, &R[j][i], &inv[j]);
+ SUFFIX(go_quad_sub) (&d, &d, &p);
+ }
+ SUFFIX(go_quad_div) (&inv[i], &d, &R[i][i]);
}
- if (inv[k] < 0) {
+ if (SUFFIX(go_quad_value) (&inv[k]) < 0) {
/*
* If this happens, something is really
* wrong, numerically.
*/
g_printerr ("inv[%d]=%" FORMAT_g "\n",
- k, inv[k]);
- regerr = GO_REG_near_singular_bad;
+ k,
+ SUFFIX(go_quad_value) (&inv[k]));
+ //regerr = GO_REG_near_singular_bad;
+ }
+ {
+ QUAD p;
+ SUFFIX(go_quad_mul) (&p, &N2, &inv[k]);
+ SUFFIX(go_quad_sqrt) (&p, &p);
+ stat_->se[k] = SUFFIX(go_quad_value) (&p);
}
- stat_->se[k] = SUFFIX(sqrt) (stat_->var * inv[k]);
}
g_free (inv);
@@ -802,6 +857,9 @@ SUFFIX(general_linear_regression) (CONSTMATRIX xss, int xdim,
: stat_->ss_resid / stat_->df_resid;
}
+ SUFFIX(go_quad_end) (state);
+
+ g_free (qresult);
FREE_MATRIX (Q, xdim, n);
FREE_MATRIX (R, xdim, xdim);
#endif
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]