Mercurial > hg > qm-dsp
diff maths/Polyfit.h @ 477:fa407c1d9923
Untabify + indent
author | Chris Cannam <cannam@all-day-breakfast.com> |
---|---|
date | Thu, 30 May 2019 18:36:48 +0100 |
parents | 7b7691b49197 |
children | 7132d952b9a6 |
line wrap: on
line diff
--- a/maths/Polyfit.h Thu May 30 18:30:58 2019 +0100 +++ b/maths/Polyfit.h Thu May 30 18:36:48 2019 +0100 @@ -50,8 +50,8 @@ public: static double PolyFit2 (const vector<double> &x, // does the work - const vector<double> &y, - vector<double> &coef); + const vector<double> &y, + vector<double> &coef); private: @@ -61,21 +61,21 @@ static void Square (const Matrix &x, // Matrix multiplication routine - const vector<double> &y, - Matrix &a, // A = transpose X times X - vector<double> &g, // G = Y times X - const int nrow, const int ncol); + const vector<double> &y, + Matrix &a, // A = transpose X times X + vector<double> &g, // G = Y times X + const int nrow, const int ncol); // Forms square coefficient matrix static bool GaussJordan (Matrix &b, // square matrix of coefficients - const vector<double> &y, // constant vector - vector<double> &coef); // solution vector + const vector<double> &y, // constant vector + vector<double> &coef); // solution vector // returns false if matrix singular static bool GaussJordan2(Matrix &b, - const vector<double> &y, - Matrix &w, - vector<vector<int> > &index); + const vector<double> &y, + Matrix &w, + vector<vector<int> > &index); }; // some utility functions @@ -114,8 +114,8 @@ // main PolyFit routine double TPolyFit::PolyFit2 (const vector<double> &x, - const vector<double> &y, - vector<double> &coefs) + const vector<double> &y, + vector<double> &coefs) // nterms = coefs.size() // npoints = x.size() { @@ -144,38 +144,38 @@ } for(i = 0; i < npoints; ++i) { - // { setup x matrix } - xi = x[i]; - xmatr[i][0] = 1.0; // { first column } - for(j = 1; j < nterms; ++j) - xmatr[i][j] = xmatr [i][j - 1] * xi; + // { setup x matrix } + xi = x[i]; + xmatr[i][0] = 1.0; // { first column } + for(j = 1; j < nterms; ++j) + xmatr[i][j] = xmatr [i][j - 1] * xi; } Square (xmatr, y, a, g, npoints, nterms); if(!GaussJordan (a, g, coefs)) - return -1; + return -1; sum_y = 0.0; sum_y2 = 0.0; srs = 0.0; for(i = 0; i < npoints; ++i) { - yi = y[i]; - yc = 0.0; - for(j = 0; j < nterms; ++j) - yc += coefs [j] * xmatr [i][j]; - srs += NSUtility::sqr (yc - yi); - sum_y += yi; - sum_y2 += yi * yi; + yi = y[i]; + yc = 0.0; + for(j = 0; j < nterms; ++j) + yc += coefs [j] * xmatr [i][j]; + srs += NSUtility::sqr (yc - yi); + sum_y += yi; + sum_y2 += yi * yi; } // If all Y values are the same, avoid dividing by zero correl_coef = sum_y2 - NSUtility::sqr (sum_y) / npoints; // Either return 0 or the correct value of correlation coefficient if (correl_coef != 0) - correl_coef = srs / correl_coef; + correl_coef = srs / correl_coef; if (correl_coef >= 1) - correl_coef = 0.0; + correl_coef = 0.0; else - correl_coef = sqrt (1.0 - correl_coef); + correl_coef = sqrt (1.0 - correl_coef); return correl_coef; } @@ -189,36 +189,36 @@ // Form square coefficient matrix void TPolyFit::Square (const Matrix &x, - const vector<double> &y, - Matrix &a, - vector<double> &g, - const int nrow, - const int ncol) + const vector<double> &y, + Matrix &a, + vector<double> &g, + const int nrow, + const int ncol) { int i, k, l; for(k = 0; k < ncol; ++k) { - for(l = 0; l < k + 1; ++l) - { - a [k][l] = 0.0; - for(i = 0; i < nrow; ++i) - { - a[k][l] += x[i][l] * x [i][k]; - if(k != l) - a[l][k] = a[k][l]; - } - } - g[k] = 0.0; - for(i = 0; i < nrow; ++i) - g[k] += y[i] * x[i][k]; + for(l = 0; l < k + 1; ++l) + { + a [k][l] = 0.0; + for(i = 0; i < nrow; ++i) + { + a[k][l] += x[i][l] * x [i][k]; + if(k != l) + a[l][k] = a[k][l]; + } + } + g[k] = 0.0; + for(i = 0; i < nrow; ++i) + g[k] += y[i] * x[i][k]; } } //--------------------------------------------------------------------------------- bool TPolyFit::GaussJordan (Matrix &b, - const vector<double> &y, - vector<double> &coef) + const vector<double> &y, + vector<double> &coef) //b square matrix of coefficients //y constant vector //coef solution vector @@ -247,44 +247,44 @@ NSUtility::zeroise(index, ncol, 3); if(!GaussJordan2(b, y, w, index)) - return false; + return false; // Interchange columns int m; for (int i = 0; i < ncol; ++i) { - m = ncol - i - 1; - if(index [m][0] != index [m][1]) - { - irow = index [m][0]; - icol = index [m][1]; - for(int k = 0; k < ncol; ++k) - swap (b[k][irow], b[k][icol]); - } + m = ncol - i - 1; + if(index [m][0] != index [m][1]) + { + irow = index [m][0]; + icol = index [m][1]; + for(int k = 0; k < ncol; ++k) + swap (b[k][irow], b[k][icol]); + } } for(int k = 0; k < ncol; ++k) { - if(index [k][2] != 0) - { + if(index [k][2] != 0) + { std::cerr << "ERROR: Error in PolyFit::GaussJordan: matrix is singular" << std::endl; return false; - } + } } for( int i = 0; i < ncol; ++i) - coef[i] = w[i][0]; + coef[i] = w[i][0]; return true; -} // end; { procedure GaussJordan } +} // end; { procedure GaussJordan } //---------------------------------------------------------------------------------------------- bool TPolyFit::GaussJordan2(Matrix &b, - const vector<double> &y, - Matrix &w, - vector<vector<int> > &index) + const vector<double> &y, + Matrix &w, + vector<vector<int> > &index) { //GaussJordan2; // first half of GaussJordan // actual start of gaussj @@ -297,75 +297,75 @@ int nv = 1; // single constant vector for(int i = 0; i < ncol; ++i) { - w[i][0] = y[i]; // copy constant vector - index[i][2] = -1; + w[i][0] = y[i]; // copy constant vector + index[i][2] = -1; } determ = 1.0; for(int i = 0; i < ncol; ++i) { - // Search for largest element - big = 0.0; - for(int j = 0; j < ncol; ++j) - { - if(index[j][2] != 0) - { - for(int k = 0; k < ncol; ++k) - { - if(index[k][2] > 0) { + // Search for largest element + big = 0.0; + for(int j = 0; j < ncol; ++j) + { + if(index[j][2] != 0) + { + for(int k = 0; k < ncol; ++k) + { + if(index[k][2] > 0) { std::cerr << "ERROR: Error in PolyFit::GaussJordan2: matrix is singular" << std::endl; return false; } - if(index[k][2] < 0 && fabs(b[j][k]) > big) - { - irow = j; - icol = k; - big = fabs(b[j][k]); - } - } // { k-loop } - } - } // { j-loop } - index [icol][2] = index [icol][2] + 1; - index [i][0] = irow; - index [i][1] = icol; + if(index[k][2] < 0 && fabs(b[j][k]) > big) + { + irow = j; + icol = k; + big = fabs(b[j][k]); + } + } // { k-loop } + } + } // { j-loop } + index [icol][2] = index [icol][2] + 1; + index [i][0] = irow; + index [i][1] = icol; - // Interchange rows to put pivot on diagonal - // GJ3 - if(irow != icol) - { - determ = -determ; - for(int m = 0; m < ncol; ++m) - swap (b [irow][m], b[icol][m]); - if (nv > 0) - for (int m = 0; m < nv; ++m) - swap (w[irow][m], w[icol][m]); - } // end GJ3 + // Interchange rows to put pivot on diagonal + // GJ3 + if(irow != icol) + { + determ = -determ; + for(int m = 0; m < ncol; ++m) + swap (b [irow][m], b[icol][m]); + if (nv > 0) + for (int m = 0; m < nv; ++m) + swap (w[irow][m], w[icol][m]); + } // end GJ3 - // divide pivot row by pivot column - pivot = b[icol][icol]; - determ *= pivot; - b[icol][icol] = 1.0; + // divide pivot row by pivot column + pivot = b[icol][icol]; + determ *= pivot; + b[icol][icol] = 1.0; - for(int m = 0; m < ncol; ++m) - b[icol][m] /= pivot; - if(nv > 0) - for(int m = 0; m < nv; ++m) - w[icol][m] /= pivot; + for(int m = 0; m < ncol; ++m) + b[icol][m] /= pivot; + if(nv > 0) + for(int m = 0; m < nv; ++m) + w[icol][m] /= pivot; - // Reduce nonpivot rows - for(int n = 0; n < ncol; ++n) - { - if(n != icol) - { - t = b[n][icol]; - b[n][icol] = 0.0; - for(int m = 0; m < ncol; ++m) - b[n][m] -= b[icol][m] * t; - if(nv > 0) - for(int m = 0; m < nv; ++m) - w[n][m] -= w[icol][m] * t; - } - } + // Reduce nonpivot rows + for(int n = 0; n < ncol; ++n) + { + if(n != icol) + { + t = b[n][icol]; + b[n][icol] = 0.0; + for(int m = 0; m < ncol; ++m) + b[n][m] -= b[icol][m] * t; + if(nv > 0) + for(int m = 0; m < nv; ++m) + w[n][m] -= w[icol][m] * t; + } + } } // { i-loop } return true; }