1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
2 //---------------------------------------------------------------------------
6 //---------------------------------------------------------------------------
7 // Least-squares curve fitting class for arbitrary data types
10 { ******************************************
11 **** Scientific Subroutine Library ****
12 **** for C++ Builder ****
13 ******************************************
15 The following programs were written by Allen Miller and appear in the
16 book entitled "Pascal Programs For Scientists And Engineers" which is
17 published by Sybex, 1981.
18 They were originally typed and submitted to MTPUG in Oct. 1982
22 They have had minor corrections and adaptations for Turbo Pascal by
28 2000 Oct 28 Updated for Delphi 4, open array parameters.
29 This allows the routine to be generalised so that it is no longer
30 hard-coded to make a specific order of best fit or work with a
31 specific number of points.
32 2001 Jan 07 Update Web site address
34 Copyright � David J Taylor, Edinburgh and others listed above
35 Web site: www.satsignal.net
36 E-mail: davidtaylor@writeme.com
39 ///////////////////////////////////////////////////////////////////////////////
40 // Modified by CLandone for VC6 Aug 2004
41 ///////////////////////////////////////////////////////////////////////////////
49 typedef vector<vector<double> > Matrix;
52 static double PolyFit2 (const vector<double> &x, // does the work
53 const vector<double> &y,
54 vector<double> &coef);
58 TPolyFit &operator = (const TPolyFit &); // disable assignment
59 TPolyFit(); // and instantiation
60 TPolyFit(const TPolyFit&); // and copying
63 static void Square (const Matrix &x, // Matrix multiplication routine
64 const vector<double> &y,
65 Matrix &a, // A = transpose X times X
66 vector<double> &g, // G = Y times X
67 const int nrow, const int ncol);
68 // Forms square coefficient matrix
70 static bool GaussJordan (Matrix &b, // square matrix of coefficients
71 const vector<double> &y, // constant vector
72 vector<double> &coef); // solution vector
73 // returns false if matrix singular
75 static bool GaussJordan2(Matrix &b,
76 const vector<double> &y,
78 vector<vector<int> > &index);
81 // some utility functions
85 static void swap(double &a, double &b) {double t = a; a = b; b = t;}
86 // fills a vector with zeros.
87 static void zeroise(vector<double> &array, int n) {
89 for(int j = 0; j < n; ++j) array.push_back(0);
91 // fills a vector with zeros.
92 static void zeroise(vector<int> &array, int n) {
94 for(int j = 0; j < n; ++j) array.push_back(0);
96 // fills a (m by n) matrix with zeros.
97 static void zeroise(vector<vector<double> > &matrix, int m, int n) {
101 for(int j = 0; j < m; ++j) matrix.push_back(zero);
103 // fills a (m by n) matrix with zeros.
104 static void zeroise(vector<vector<int> > &matrix, int m, int n) {
108 for(int j = 0; j < m; ++j) matrix.push_back(zero);
110 static double sqr(const double &x) {return x * x;}
114 // main PolyFit routine
116 double TPolyFit::PolyFit2 (const vector<double> &x,
117 const vector<double> &y,
118 vector<double> &coefs)
119 // nterms = coefs.size()
120 // npoints = x.size()
123 double xi, yi, yc, srs, sum_y, sum_y2;
124 Matrix xmatr; // Data matrix
126 vector<double> g; // Constant vector
127 const int npoints(x.size());
128 const int nterms(coefs.size());
130 NSUtility::zeroise(g, nterms);
131 NSUtility::zeroise(a, nterms, nterms);
132 NSUtility::zeroise(xmatr, npoints, nterms);
134 std::cerr << "ERROR: PolyFit called with less than one term" << std::endl;
138 std::cerr << "ERROR: PolyFit called with less than two points" << std::endl;
141 if(npoints != (int)y.size()) {
142 std::cerr << "ERROR: PolyFit called with x and y of unequal size" << std::endl;
145 for(i = 0; i < npoints; ++i)
147 // { setup x matrix }
149 xmatr[i][0] = 1.0; // { first column }
150 for(j = 1; j < nterms; ++j)
151 xmatr[i][j] = xmatr [i][j - 1] * xi;
153 Square (xmatr, y, a, g, npoints, nterms);
154 if(!GaussJordan (a, g, coefs))
159 for(i = 0; i < npoints; ++i)
163 for(j = 0; j < nterms; ++j)
164 yc += coefs [j] * xmatr [i][j];
165 srs += NSUtility::sqr (yc - yi);
170 // If all Y values are the same, avoid dividing by zero
171 correl_coef = sum_y2 - NSUtility::sqr (sum_y) / npoints;
172 // Either return 0 or the correct value of correlation coefficient
173 if (correl_coef != 0)
174 correl_coef = srs / correl_coef;
175 if (correl_coef >= 1)
178 correl_coef = sqrt (1.0 - correl_coef);
183 //------------------------------------------------------------------------
185 // Matrix multiplication routine
186 // A = transpose X times X
189 // Form square coefficient matrix
191 void TPolyFit::Square (const Matrix &x,
192 const vector<double> &y,
199 for(k = 0; k < ncol; ++k)
201 for(l = 0; l < k + 1; ++l)
204 for(i = 0; i < nrow; ++i)
206 a[k][l] += x[i][l] * x [i][k];
212 for(i = 0; i < nrow; ++i)
213 g[k] += y[i] * x[i][k];
216 //---------------------------------------------------------------------------------
219 bool TPolyFit::GaussJordan (Matrix &b,
220 const vector<double> &y,
221 vector<double> &coef)
222 //b square matrix of coefficients
224 //coef solution vector
225 //ncol order of matrix got from b.size()
230 { Gauss Jordan matrix inversion and solution }
232 { B (n, n) coefficient matrix becomes inverse }
233 { Y (n) original constant vector }
234 { W (n, m) constant vector(s) become solution vector }
235 { DETERM is the determinant }
236 { ERROR = 1 if singular }
238 { NV is number of constant vectors }
243 vector<vector<int> >index;
246 NSUtility::zeroise(w, ncol, ncol);
247 NSUtility::zeroise(index, ncol, 3);
249 if(!GaussJordan2(b, y, w, index))
252 // Interchange columns
254 for (int i = 0; i < ncol; ++i)
257 if(index [m][0] != index [m][1])
261 for(int k = 0; k < ncol; ++k)
262 swap (b[k][irow], b[k][icol]);
266 for(int k = 0; k < ncol; ++k)
268 if(index [k][2] != 0)
270 std::cerr << "ERROR: Error in PolyFit::GaussJordan: matrix is singular" << std::endl;
275 for( int i = 0; i < ncol; ++i)
280 } // end; { procedure GaussJordan }
281 //----------------------------------------------------------------------------------------------
284 bool TPolyFit::GaussJordan2(Matrix &b,
285 const vector<double> &y,
287 vector<vector<int> > &index)
289 //GaussJordan2; // first half of GaussJordan
290 // actual start of gaussj
295 int irow = 0, icol = 0;
297 int nv = 1; // single constant vector
298 for(int i = 0; i < ncol; ++i)
300 w[i][0] = y[i]; // copy constant vector
304 for(int i = 0; i < ncol; ++i)
306 // Search for largest element
308 for(int j = 0; j < ncol; ++j)
312 for(int k = 0; k < ncol; ++k)
314 if(index[k][2] > 0) {
315 std::cerr << "ERROR: Error in PolyFit::GaussJordan2: matrix is singular" << std::endl;
319 if(index[k][2] < 0 && fabs(b[j][k]) > big)
328 index [icol][2] = index [icol][2] + 1;
332 // Interchange rows to put pivot on diagonal
337 for(int m = 0; m < ncol; ++m)
338 swap (b [irow][m], b[icol][m]);
340 for (int m = 0; m < nv; ++m)
341 swap (w[irow][m], w[icol][m]);
344 // divide pivot row by pivot column
345 pivot = b[icol][icol];
349 for(int m = 0; m < ncol; ++m)
352 for(int m = 0; m < nv; ++m)
355 // Reduce nonpivot rows
356 for(int n = 0; n < ncol; ++n)
362 for(int m = 0; m < ncol; ++m)
363 b[n][m] -= b[icol][m] * t;
365 for(int m = 0; m < nv; ++m)
366 w[n][m] -= w[icol][m] * t;