From 7b22b8cd929d1dded5b678076875d393c27e9d3d Mon Sep 17 00:00:00 2001
From: curt <curt>
Date: Sat, 13 Mar 1999 17:34:44 +0000
Subject: [PATCH] Moved to math subdirectory.

---
 Math/Makefile.am  |   1 +
 Math/leastsqs.cxx | 152 ++++++++++++++++++++++++++++++++++++++++++++++
 Math/leastsqs.hxx |  90 +++++++++++++++++++++++++++
 3 files changed, 243 insertions(+)
 create mode 100644 Math/leastsqs.cxx
 create mode 100644 Math/leastsqs.hxx

diff --git a/Math/Makefile.am b/Math/Makefile.am
index e6654827d..eda2bc005 100644
--- a/Math/Makefile.am
+++ b/Math/Makefile.am
@@ -8,6 +8,7 @@ libMath_a_SOURCES = \
 	fg_geodesy.cxx fg_geodesy.hxx \
 	fg_random.c fg_random.h \
 	interpolater.cxx interpolater.hxx \
+	leastsqs.cxx leastsqs.hxx \
 	mat3.h mat3defs.h mat3err.h \
 	point3d.hxx \
 	polar3d.cxx polar3d.hxx \
diff --git a/Math/leastsqs.cxx b/Math/leastsqs.cxx
new file mode 100644
index 000000000..28e2bd35a
--- /dev/null
+++ b/Math/leastsqs.cxx
@@ -0,0 +1,152 @@
+// leastsqs.c -- Implements a simple linear least squares best fit routine
+//
+// Written by Curtis Olson, started September 1997.
+//
+// Copyright (C) 1997  Curtis L. Olson  - curt@infoplane.com
+//
+// 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.
+//
+// $Id$
+// (Log is kept at end of this file)
+//
+
+
+#include <stdio.h>
+
+#include "leastsqs.hxx"
+
+
+/* 
+Least squares fit:
+
+y = b0 + b1x
+
+     n*sum(xi*yi) - (sum(xi)*sum(yi))
+b1 = --------------------------------
+     n*sum(xi^2) - (sum(xi))^2
+
+
+b0 = sum(yi)/n - b1*(sum(xi)/n)
+*/
+
+double sum_xi, sum_yi, sum_xi_2, sum_xi_yi;
+int sum_n;
+
+void least_squares(double *x, double *y, int n, double *m, double *b) {
+    int i;
+
+    sum_xi = sum_yi = sum_xi_2 = sum_xi_yi = 0.0;
+    sum_n = n;
+
+    for ( i = 0; i < n; i++ ) {
+	sum_xi += x[i];
+	sum_yi += y[i];
+	sum_xi_2 += x[i] * x[i];
+	sum_xi_yi += x[i] * y[i];
+    }
+
+    /* printf("sum(xi)=%.2f  sum(yi)=%.2f  sum(xi^2)=%.2f  sum(xi*yi)=%.2f\n",
+	   sum_xi, sum_yi, sum_xi_2, sum_xi_yi); */
+
+    *m = ( (double)sum_n * sum_xi_yi - sum_xi * sum_yi ) / 
+	( (double)sum_n * sum_xi_2 - sum_xi * sum_xi );
+    *b = (sum_yi / (double)sum_n) - (*m) * (sum_xi / (double)sum_n);
+
+    /* printf("slope = %.2f  intercept = %.2f\n", *m, *b); */
+}
+
+
+/* incrimentally update existing values with a new data point */
+void least_squares_update(double x, double y, double *m, double *b) {
+    ++sum_n;
+
+    sum_xi += x;
+    sum_yi += y;
+    sum_xi_2 += x * x;
+    sum_xi_yi += x * y;
+
+    /* printf("sum(xi)=%.2f  sum(yi)=%.2f  sum(xi^2)=%.2f  sum(xi*yi)=%.2f\n",
+	   sum_xi, sum_yi, sum_xi_2, sum_xi_yi); */
+
+    *m = ( (double)sum_n * sum_xi_yi - sum_xi * sum_yi ) / 
+	( (double)sum_n * sum_xi_2 - sum_xi * sum_xi );
+    *b = (sum_yi / (double)sum_n) - (*m) * (sum_xi / (double)sum_n);
+
+    /* printf("slope = %.2f  intercept = %.2f\n", *m, *b); */
+}
+
+
+/* 
+  return the least squares error:
+
+              (y[i] - y_hat[i])^2
+              -------------------
+                      n
+ */
+double least_squares_error(double *x, double *y, int n, double m, double b) {
+    int i;
+    double error, sum;
+
+    sum = 0.0;
+
+    for ( i = 0; i < n; i++ ) {
+	error = y[i] - (m * x[i] + b);
+	sum += error * error;
+	// printf("%.2f %.2f\n", error, sum);
+    }
+
+    return ( sum / (double)n );
+}
+
+
+/* 
+  return the maximum least squares error:
+
+              (y[i] - y_hat[i])^2
+ */
+double least_squares_max_error(double *x, double *y, int n, double m, double b){
+    int i;
+    double error, max_error;
+
+    max_error = 0.0;
+
+    for ( i = 0; i < n; i++ ) {
+	error = y[i] - (m * x[i] + b);
+	error = error * error;
+	if ( error > max_error ) {
+	    max_error = error;
+	}
+    }
+
+    return ( max_error );
+}
+
+
+// $Log$
+// Revision 1.1  1999/03/13 17:34:45  curt
+// Moved to math subdirectory.
+//
+// Revision 1.2  1998/04/21 17:03:41  curt
+// Prepairing for C++ integration.
+//
+// Revision 1.1  1998/04/08 22:57:24  curt
+// Adopted Gnu automake/autoconf system.
+//
+// Revision 1.1  1998/03/19 02:54:47  curt
+// Reorganized into a class lib called fgDEM.
+//
+// Revision 1.1  1997/10/13 17:02:35  curt
+// Initial revision.
+//
diff --git a/Math/leastsqs.hxx b/Math/leastsqs.hxx
new file mode 100644
index 000000000..d8b40c80e
--- /dev/null
+++ b/Math/leastsqs.hxx
@@ -0,0 +1,90 @@
+// leastsqs.h -- Implements a simple linear least squares best fit routine
+//
+// Written by Curtis Olson, started September 1997.
+//
+// Copyright (C) 1997  Curtis L. Olson  - curt@infoplane.com
+//
+// 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.
+//
+// $Id$
+// (Log is kept at end of this file)
+///
+
+
+#ifndef _LEASTSQS_H
+#define _LEASTSQS_H
+
+
+#ifndef __cplusplus                                                          
+# error This library requires C++
+#endif                                   
+
+
+/* 
+Least squares fit:
+
+y = b0 + b1x
+
+     n*sum(xi*yi) - (sum(xi)*sum(yi))
+b1 = --------------------------------
+     n*sum(xi^2) - (sum(xi))^2
+
+
+b0 = sum(yi)/n - b1*(sum(xi)/n)
+*/
+
+void least_squares(double *x, double *y, int n, double *m, double *b);
+
+/* incrimentally update existing values with a new data point */
+void least_squares_update(double x, double y, double *m, double *b);
+
+
+/* 
+  return the least squares error:
+
+              (y[i] - y_hat[i])^2
+              -------------------
+                      n
+*/
+double least_squares_error(double *x, double *y, int n, double m, double b);
+
+
+/* 
+  return the maximum least squares error:
+
+              (y[i] - y_hat[i])^2
+*/
+double least_squares_max_error(double *x, double *y, int n, double m, double b);
+
+
+#endif // _LEASTSQS_H
+
+
+// $Log$
+// Revision 1.1  1999/03/13 17:34:45  curt
+// Moved to math subdirectory.
+//
+// Revision 1.2  1998/04/21 17:03:42  curt
+// Prepairing for C++ integration.
+//
+// Revision 1.1  1998/04/08 22:57:25  curt
+// Adopted Gnu automake/autoconf system.
+//
+// Revision 1.1  1998/03/19 02:54:48  curt
+// Reorganized into a class lib called fgDEM.
+//
+// Revision 1.1  1997/10/13 17:02:35  curt
+// Initial revision.
+//