cpolyfit.cpp

多项式曲线拟合 任意介数 Purpose - Least-squares curve fit of arbitrary order working in C++ Builder 2007 as

//---------------------------------------------------------------------------

#pragma hdrstop

#include "CPolyfit.h"

//---------------------------------------------------------------------------

#pragma package(smart_init)

/*
=============================================================================
 Purpose - Least-squares curve fit of arbitrary order

 working in C++ Builder 2007 as a template class,
 using vector<FloatType> parameters.
 Added a method to handle some EMathError exceptions.
 If do NOT want to use this just call PolyFit2 directly.

 usage:      Call PolyFit by something like this.

			 CPolyFit<double> PolyFitObj;
			 double correlation_coefficiant = PolyFitObj.PolyFit(X, Y, A);

			 where X and Y are vectors of doubles which must have the same size and
			 A is a vector of doubles which must be the same size as the number of
			 coefficients required.

 returns:    The correlation coefficient or -1 on failure.
 produces:   A vector (A) which holds the coefficients.
=============================================================================
*/

using namespace NSUtility;

template <class FloatType>
FloatType __fastcall CPolyFit<FloatType>::PolyFit (const vector<FloatType> &x,
                                                   const vector<FloatType> &y,
                                                   vector<FloatType> &coefs)
                                                   // npoints = x.size()
                                                   // nterms = coefs.size()
{
 try
 {
  return PolyFit2(x, y, coefs);   // the main polyfit function. Returns correlation
                                  // coefficient on success.
 }
 catch(EMathError &exception)
 {
  // A maths error (EOverFlow, EUnderFlow or EZeroDivide has been thrown.
  // we will normalise the data and try again. This quite often solves
  // the problem.
  vector<FloatType> s(x);      // make a copy of the data
  vector<FloatType> t(y);
  FloatType x0, xm, y0, ym;    // x0 will be min(x) xm = max(x) etc.
  Normalise(s, t, x0, xm, y0, ym);
  try                       // try again using the normalised data
  {
   FloatType coef = PolyFit2(s, t, coefs);
   if(coef > 0)
    RestoreCoeffs(coefs, x0, xm, y0, ym);
        // the coefficients now refer to the normalised data.
        // so we have to convert them to refer to the origonal
        // data. - This was more awkward than I envisaged. RHK
   return coef;
  }
  catch(...){return -1;}    // give up
 }
 catch(...) {return -1;}    // will arrive here if an unhandled exception is thrown

 // Shouldn't get here.
 // return -1;   // remove the '//' if your compiler gives you gyp.
}
//------------------------------------------------------------------------------------------

template <class FloatType>
void __fastcall CPolyFit<FloatType>::Normalise(vector<FloatType> &s,
                                     vector<FloatType> &t,
                                     FloatType &x0,
                                     FloatType &xm,
                                     FloatType &y0,
                                     FloatType &ym)
{
 x0 = *min_element(s.begin(), s.end());
 xm = *max_element(s.begin(), s.end());
 y0 = *min_element(t.begin(), t.end());
 ym = *max_element(t.begin(), t.end());
 size_t n(s.size());
 for(size_t i = 0; i < n; ++i)
 {
  s[i] = (s[i] - x0) / (xm - x0);
  t[i] = (t[i] - y0) / (ym - y0);
 }
}
//-----------------------------------------------------------------------------------------

template <class FloatType>
void __fastcall CPolyFit<FloatType>::RestoreCoeffs(vector<FloatType> &b,  // coefficients
                                                   FloatType x0,  // min x
                                                   FloatType xm,  // max x
                                                   FloatType y0,
                                                   FloatType ym)
// we have coefficients bi such that
// t = b0 + b1 * s + b2 * s^2 +  .... + bn * s^n
// where t = (y - y0) / (ym - y0) and s = (x - x0) / (xm - x0)
//
// what we want is ai such that
// y = a0 + a1 * x + a2 * x^2 + .... + an * x^n
//
// substituting for s and t and looking at each power of x
// we get a horrible expression involving only known numbers.
// Fortunately it is easily solved by the following algorithm.
//
{
 vector<FloatType> a;
 size_t n = b.size();
 vector<vector<int> > p;
 GeneratePascalsTriangle(p, n);   // (1, 1 1, 1 2 1, 1 3 3 1, 1 4 6 4 1, etc. )

 FloatType yr = ym - y0;
 FloatType xr = x0 - xm;
 FloatType sum;
 for(size_t i = 0; i < n; ++i)
 {
  sum = 0.0;
  for(size_t j = 0; j < n - i; ++j)
  {
   sum += b[i + j] * p[i + j][j] * Pow(x0 / xr, j);
  }
  sum *= yr / Pow(-xr , i);
  a.push_back(sum);
 }
 a[0] += y0;
 b.assign(a.begin(), a.end());
}
//------------------------------------------------------------------------------------------

template <class FloatType>
void __fastcall CPolyFit<FloatType>::GeneratePascalsTriangle
                                     (vector<vector<int> > &Pascal, size_t n)
{
 vector<int> v;
 v.push_back(1);

 Pascal.clear();
 Pascal.push_back(v);

 int sum;
 for(size_t i = 0; i < n; ++i)
 {
  v.clear();
  for(size_t j = 0; j <= i + 1; ++j)
  {
   if(j == 0)
    sum = Pascal[i][j];
   else if(j == i + 1)
    sum = Pascal[i][j - 1];
   else
    sum = Pascal[i][j] + Pascal[i][j - 1];
   v.push_back(sum);
  }
  v.push_back(1);
  Pascal.push_back(v);
 }
}
//----------------------------------------------------------------------------------------

// main PolyFit routine
template <class FloatType>
FloatType __fastcall CPolyFit<FloatType>::PolyFit2 (const vector<FloatType> &x,
                                                    const vector<FloatType> &y,
                                                    vector<FloatType> &coefs)
                                                    // nterms = coefs.size()
                                                    // npoints = x.size()
{
 try
 {
  size_t i, j;
  FloatType xi, yi, yc, srs, sum_y, sum_y2;
  Matrix xmatr;        // Data matrix
  Matrix a;
  vector<FloatType> g;      // Constant vector
  const size_t npoints(x.size());
  const size_t nterms(coefs.size());
  FloatType correl_coef;
  Zeroise(g, nterms);
  Zeroise(a, nterms, nterms);
  Zeroise(xmatr, npoints, nterms);
  if(nterms < 1)
   throw "PolyFit called with less than one term";
  if(npoints < 2)
   throw "PolyFit called with less than two points";
  if(npoints != y.size())
   throw "PolyFit called with x and y of unequal size";
  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;
  }
  Square (xmatr, y, a, g, npoints, nterms);
  if(!GaussJordan (a, g, coefs))
   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 += 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 - Sqr (sum_y) / npoints;
  // Either return 0 or the correct value of correlation coefficient
  if (correl_coef != 0)
   correl_coef = srs / correl_coef;
  if (correl_coef >= 1)
   correl_coef = 0.0;
  else
   correl_coef = sqrt (1.0 - correl_coef);
  return correl_coef;
 }
 catch(char *b)
 {
  Application->MessageBox(b, "PolyFit Exception",
   MB_OK | MB_ICONEXCLAMATION | MB_DEFBUTTON1 | MB_APPLMODAL | MB_SETFOREGROUND | MB_TOPMOST);
  return -1;
 }
}
//------------------------------------------------------------------------

 // Matrix multiplication routine
 // A = transpose X times X
 // G = Y times X

 // Form square coefficient matrix
template <class FloatType>
void __fastcall CPolyFit<FloatType>::Square (const Matrix &x,
                                             const vector<FloatType> &y,
                                             Matrix &a,
                                             vector<FloatType> &g,
                                             const size_t nrow,
                                             const size_t ncol)
{
 size_t 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];
 }
}