49 typedef vector<vector<double> >
Matrix;
52 static double PolyFit2 (
const vector<double> &x,
53 const vector<double> &y,
54 vector<double> &coef);
63 static void Square (
const Matrix &x,
64 const vector<double> &y,
67 const int nrow,
const int ncol);
71 const vector<double> &y,
72 vector<double> &coef);
76 const vector<double> &y,
78 vector<vector<int> > &index);
85 inline void swap(
double &a,
double &b) {
double t = a; a = b; b = t;}
86 void zeroise(vector<double> &array,
int n);
87 void zeroise(vector<int> &array,
int n);
88 void zeroise(vector<vector<double> > &matrix,
int m,
int n);
89 void zeroise(vector<vector<int> > &matrix,
int m,
int n);
90 inline double sqr(
const double &x) {
return x * x;}
103 const vector<double> &y,
104 vector<double> &coefs)
109 double xi, yi, yc, srs, sum_y, sum_y2;
113 const int npoints(x.size());
114 const int nterms(coefs.size());
118 zeroise(xmatr, npoints, nterms);
120 std::cerr <<
"ERROR: PolyFit called with less than one term" << std::endl;
124 std::cerr <<
"ERROR: PolyFit called with less than two points" << std::endl;
127 if(npoints != (
int)y.size()) {
128 std::cerr <<
"ERROR: PolyFit called with x and y of unequal size" << std::endl;
131 for(i = 0; i < npoints; ++i)
136 for(j = 1; j < nterms; ++j)
137 xmatr[i][j] = xmatr [i][j - 1] * xi;
139 Square (xmatr, y, a, g, npoints, nterms);
145 for(i = 0; i < npoints; ++i)
149 for(j = 0; j < nterms; ++j)
150 yc += coefs [j] * xmatr [i][j];
151 srs +=
sqr (yc - yi);
157 correl_coef = sum_y2 -
sqr (sum_y) / npoints;
159 if (correl_coef != 0)
160 correl_coef = srs / correl_coef;
161 if (correl_coef >= 1)
164 correl_coef = sqrt (1.0 - correl_coef);
178 const vector<double> &y,
185 for(k = 0; k < ncol; ++k)
187 for(l = 0; l < k + 1; ++l)
190 for(i = 0; i < nrow; ++i)
192 a[k][l] += x[i][l] * x [i][k];
198 for(i = 0; i < nrow; ++i)
199 g[k] += y[i] * x[i][k];
206 const vector<double> &y,
207 vector<double> &coef)
229 vector<vector<int> >index;
240 for (
int i = 0; i < ncol; ++i)
243 if(index [m][0] != index [m][1])
247 for(
int k = 0; k < ncol; ++k)
248 swap (b[k][irow], b[k][icol]);
252 for(
int k = 0; k < ncol; ++k)
254 if(index [k][2] != 0)
256 std::cerr <<
"ERROR: Error in PolyFit::GaussJordan: matrix is singular" << std::endl;
261 for(
int i = 0; i < ncol; ++i)
271 const vector<double> &y,
273 vector<vector<int> > &index)
284 for(
int i = 0; i < ncol; ++i)
290 for(
int i = 0; i < ncol; ++i)
294 for(
int j = 0; j < ncol; ++j)
298 for(
int k = 0; k < ncol; ++k)
300 if(index[k][2] > 0) {
301 std::cerr <<
"ERROR: Error in PolyFit::GaussJordan2: matrix is singular" << std::endl;
305 if(index[k][2] < 0 && fabs(b[j][k]) > big)
314 index [icol][2] = index [icol][2] + 1;
323 for(
int m = 0; m < ncol; ++m)
324 swap (b [irow][m], b[icol][m]);
326 for (
int m = 0; m < nv; ++m)
327 swap (w[irow][m], w[icol][m]);
331 pivot = b[icol][icol];
335 for(
int m = 0; m < ncol; ++m)
338 for(
int m = 0; m < nv; ++m)
342 for(
int n = 0; n < ncol; ++n)
348 for(
int m = 0; m < ncol; ++m)
349 b[n][m] -= b[icol][m] * t;
351 for(
int m = 0; m < nv; ++m)
352 w[n][m] -= w[icol][m] * t;
369 for(
int j = 0; j < n; ++j)
378 for(
int j = 0; j < n; ++j)
389 for(
int j = 0; j < m; ++j)
390 matrix.push_back(zero);
400 for(
int j = 0; j < m; ++j)
401 matrix.push_back(zero);
static double PolyFit2(const vector< double > &x, const vector< double > &y, vector< double > &coef)
double sqr(const double &x)
void swap(double &a, double &b)
vector< vector< double > > Matrix
void zeroise(vector< double > &array, int n)
static void Square(const Matrix &x, const vector< double > &y, Matrix &a, vector< double > &g, const int nrow, const int ncol)
TPolyFit & operator=(const TPolyFit &)
static bool GaussJordan2(Matrix &b, const vector< double > &y, Matrix &w, vector< vector< int > > &index)
static bool GaussJordan(Matrix &b, const vector< double > &y, vector< double > &coef)
void zeroise(vector< vector< int > > &matrix, int m, int n)