[goffice] Regression: add code for computing leverages.

[Morten Welinder <mortenw src gnome org>]

commit 121aab389332a25c8e0bf8f26f17de823abdb552
Author: Morten Welinder <terra gnome org>
Date:   Thu Jan 17 19:08:16 2013 -0500

    Regression: add code for computing leverages.
    
    ...without creating any m-times-m matrices.

 ChangeLog                    |   11 ++++--
 NEWS                         |    1 +
 goffice/math/go-regression.c |   73 +++++++++++++++++++++++++++++++++++++++++-
 goffice/math/go-regression.h |    5 +++
 4 files changed, 85 insertions(+), 5 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index b9b0222..5ca7d1b 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -3,12 +3,14 @@
 	* goffice/math/go-regression.c (QRH): Make qdet argument optional.
 	(back_solve, fwd_solve): New functions extracted from and used
 	throughout.
+	(go_linear_regression_leverage): New function.
 
 2013-01-16  Jean Brefort  <jean brefort normalesup org>
 
 	* goffice/graph/gog-smoothed-curve.c
-	(gog_smoothed_curve_parent_changed): call parent class implementation. Fix
-	crash after deleting a named smoothed curve. [#691796]
+	(gog_smoothed_curve_parent_changed): call parent class
+	implementation.  Fix crash after deleting a named smoothed
+	curve. [#691796]
 
 2013-01-15  Morten Welinder  <terra gnome org>
 
@@ -24,8 +26,9 @@
 
 2013-01-14  Jean Brefort  <jean brefort normalesup org>
 
-	* goffice/graph/gog-series-labels.c (used_selection_changed_cb): don't add
-	an extra separator at start of the format string. [#691704]
+	* goffice/graph/gog-series-labels.c (used_selection_changed_cb):
+	don't add an extra separator at start of the format string.
+	[#691704]
 
 2013-01-13  Morten Welinder  <terra gnome org>
 
diff --git a/NEWS b/NEWS
index 445beed..510dae2 100644
--- a/NEWS
+++ b/NEWS
@@ -7,6 +7,7 @@ Jean:
 Morten:
 	* Cleanup linear algebra.  [#691630]
 	* Work on new compiler warnings.
+	* Add efficient code for computing leverages.  [#691913]
 
 --------------------------------------------------------------------------
 goffice 0.10.0:
diff --git a/goffice/math/go-regression.c b/goffice/math/go-regression.c
index 7bbafea..f1519ff 100644
--- a/goffice/math/go-regression.c
+++ b/goffice/math/go-regression.c
@@ -523,7 +523,7 @@ SUFFIX(regres_from_condition) (CONSTQMATRIX R, int n)
 	   done with QUADs, we can afford to lose a lot.  */
 	lc = SUFFIX(log)(c) / SUFFIX(log)(FLT_RADIX);
 
-#if 0
+#if 1
 	if (lc > 0.95* DOUBLE_MANT_DIG)
 		return GO_REG_near_singular_bad;
 	if (lc > 0.75 * DOUBLE_MANT_DIG)
@@ -2026,3 +2026,74 @@ SUFFIX(go_non_linear_regression) (SUFFIX(GORegressionFunction) f,
 
 	return result;
 }
+
+/*
+ * Compute the diagonal of A (AT A)^-1 AT
+ *
+ * A is full-rank m-times-n.
+ * d is output for m-element diagonal.
+ */
+
+GORegressionResult
+SUFFIX(go_linear_regression_leverage) (MATRIX A, DOUBLE *d, int m, int n)
+{
+	QMATRIX V;
+	QMATRIX R;
+	gboolean has_result;
+	void *state = SUFFIX(go_quad_start) ();
+	GORegressionResult regres;
+
+	ALLOC_MATRIX2 (V, m, n, QUAD);
+	ALLOC_MATRIX2 (R, n, n, QUAD);
+
+	has_result = SUFFIX(QRH)(A, FALSE, V, R, m, n, NULL);
+
+	if (has_result) {
+		int k;
+		QUAD *b = g_new (QUAD, n);
+
+		regres = SUFFIX(regres_from_condition) (R, n);
+
+		for (k = 0; k < m; k++) {
+			int i;
+			QUAD acc;
+
+			/* b = AT e_k  */
+			for (i = 0; i < n; i++)
+				SUFFIX(go_quad_init) (&b[i], A[k][i]);
+
+			/* Solve R^T newb = b */
+			if (SUFFIX(fwd_solve) (R, b, b, n)) {
+				regres = GO_REG_singular;
+				break;
+			}
+
+			/* Solve R newb = b */
+			if (SUFFIX(back_solve) (R, b, b, n)) {
+				regres = GO_REG_singular;
+				break;
+			}
+
+			/* acc = (Ab)_k */
+			SUFFIX(go_quad_init) (&acc, 0);
+			for (i = 0; i < n; i++) {
+				QUAD p;
+				SUFFIX(go_quad_init) (&p, A[k][i]);
+				SUFFIX(go_quad_mul) (&p, &p, &b[i]);
+				SUFFIX(go_quad_add) (&acc, &acc, &p);
+			}
+
+			d[k] = SUFFIX(go_quad_value) (&acc);
+		}
+
+		g_free (b);
+	} else
+		regres = GO_REG_invalid_data;
+
+	FREE_MATRIX (V, n, n);
+	FREE_MATRIX (R, n, n);
+
+	SUFFIX(go_quad_end) (state);
+
+	return regres;
+}
diff --git a/goffice/math/go-regression.h b/goffice/math/go-regression.h
index 40c7aad..847cadb 100644
--- a/goffice/math/go-regression.h
+++ b/goffice/math/go-regression.h
@@ -48,6 +48,8 @@ GORegressionResult 	 go_linear_regression 		(double **xss, int dim,
 							 gboolean affine,
 							 double *res,
 							 go_regression_stat_t *stat_);
+GORegressionResult       go_linear_regression_leverage  (double **A, double *d,
+							 int m, int n);
 GORegressionResult 	 go_exponential_regression 	(double **xss, int dim,
 							 const double *ys, int n,
 							 gboolean affine,
@@ -143,6 +145,9 @@ GORegressionResult    go_linear_regressionl   	(long double **xss, int dim,
 						 gboolean affine,
 						 long double *res,
 						 go_regression_stat_tl *stat_);
+GORegressionResult    go_linear_regression_leveragel (long double **A,
+						      long double *d,
+						      int m, int n);
 GORegressionResult    go_exponential_regressionl	(long double **xss, int dim,
 						 const long double *ys, int n,
 						 gboolean affine,