|
Re: Combining Vectors and References
Hi David and Lucian,
Appreciate the thoughts. That's interesting that the segment Lucian wrote compiled under VC. This code was actually written for Borland, so you'd think it would work here! I do have the polyfit.cpp file linked in, under the "Add to Project" menu, and there is a polyfit.obj file being generated. And if I mess up the parameters, the compiler tells me that it cannot initialize the object. I didn't want to include the whole set of Polyfit files in my posting above as it is fairly large. That was a segment of the header file I showed above, thinking that the problem is in that little bit of code. The whole project is quite large including libraries to drive some National Instruments hardware and an Intel Signal Processing system (just so you know the underlying concepts are not new to me). Here are the complete Polyfit precompile header, header file and .cpp file, with thanks to the people mentioned in the .cpp file for posting this. I tried to contact Mr Kavanagh directly but the email address he gave here is no longer in service. So then I tried Mr. Taylor, but he is mostly a PASCAL guy.
Precompile Header FIle:
#ifndef PRECHH
#define PRECHH
#include <vcl.h>
#include <fstream.h>
#include <vector>
#include <string>
#include <algorithm>
#include <math.h>
#endif
Header File:
//---------------------------------------------------------------------------
#ifndef PolyfitH
#define PolyfitH
#include "PreCH.h"
//---------------------------------------------------------------------------
using std::vector;
using std::string;
template <class FloatType>
class CPolyFit
{
typedef vector<vector<FloatType> > Matrix;
public:
static FloatType __fastcall PolyFit(const vector<FloatType> &x, // independent variable
const vector<FloatType> &y, // dependent variable
vector<FloatType> &coef); // coefficients
// returns correlation coef, -1 on error
static FloatType __fastcall PolyFit2 (const vector<FloatType> &x, // does the work
const vector<FloatType> &y,
vector<FloatType> &coef);
// Least-squares fit to nrow sets of x and y pairs of points
static void __fastcall LinFit(const vector<FloatType> &x, // independent variable
vector<FloatType> &y, // dependent variable
vector<FloatType> &y_calc, // calculated dep. variable
vector<FloatType> &resid, // array of residuals
vector<FloatType> &coef, // coefficients
vector<FloatType> &sig, // error on coefficients
const size_t nrow, // length of variable array
size_t ncol); // number of terms
private:
CPolyFit &operator = (const CPolyFit &); // disable assignment
static void __fastcall Normalise(vector<FloatType> &s, // normalised x y values
vector<FloatType> &t,
FloatType &x0, // min x
FloatType &xm, // max x
FloatType &y0,
FloatType &ym);
static void __fastcall RestoreCoeffs(vector<FloatType> &coef,
FloatType x0,
FloatType xm,
FloatType y0,
FloatType ym);
static void __fastcall GeneratePascalsTriangle(vector<vector<int> > &Pascal, size_t n);
static void __fastcall Square (const Matrix &x, // Matrix multiplication routine
const vector<FloatType> &y,
Matrix &a, // A = transpose X times X
vector<FloatType> &g, // G = Y times X
const size_t nrow, const size_t ncol);
// Forms square coefficient matrix
static bool __fastcall GaussJordan (Matrix &b, // square matrix of coefficients
const vector<FloatType> &y, // constant vector
vector<FloatType> &coef); // solution vector
// returns false if matrix singular
static bool __fastcall GaussJordan2(Matrix &b,
const vector<FloatType> &y,
Matrix &w,
vector<vector<int> > &index);
};
// some utility functions
namespace NSUtility
{
template <typename X> inline void Swap(X &a, X &b) {X t = a; a = b; b = t;}
template <typename X> void Zeroise(vector<X> &array, size_t n);
template <typename X> void Zeroise(vector<vector<X> > &matrix, size_t m, size_t n);
template <typename X> inline X Sqr(const X &x) {return x * x;}
template <typename X> X Pow(const X &x, int i);
};
#endif
CPP File
//---------------------------------------------------------------------------
#pragma hdrstop
#include "Polyfit.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
// Purpose - Least-squares curve fit of arbitrary order
/*
{ ******************************************
**** Scientific Subroutine Library ****
**** for C++ Builder ****
******************************************
The following programs were written by Allen Miller and appear in the
book entitled "Pascal Programs For Scientists And Engineers" which is
published by Sybex, 1981.
They were originally typed and submitted to MTPUG in Oct. 1982
Juergen Loewner
Hoher Heckenweg 3
D-4400 Muenster
They have had minor corrections and adaptations for Turbo Pascal by
Jeff Weiss
1572 Peacock Ave.
Sunnyvale, CA 94087.
2000 Oct 28 Updated for Delphi 4, open array parameters.
This allows the routine to be generalised so that it is no longer
hard-coded to make a specific order of best fit or work with a
specific number of points.
2001 Jan 07 Update Web site address
Copyright © David J Taylor, Edinburgh and others listed above
Web site: [url]www.satsignal.net[/url]
E-mail: [email]davidtaylor@writeme.com[/email]
}
*/
// 2003 Feb 26 Updated for C++ Builder 5 as a template class,
// using vector<FloatType> parameters.
// Added a method to handle some EMathError exceptions.
// If you do not want to use this just call PolyFit2 directly.
//
// R H Kavanagh, Helwith Bridge, North Yorkshire,
// E-mail: [email]r.kavanagh@daelnet.net[/email]
//
//
// 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];
}
}
//---------------------------------------------------------------------------------
template <class FloatType>
bool __fastcall CPolyFit<FloatType>::GaussJordan (Matrix &b,
const vector<FloatType> &y,
vector<FloatType> &coef)
//b square matrix of coefficients
//y constant vector
//coef solution vector
//ncol order of matrix got from b.size()
{
/*
{ Gauss Jordan matrix inversion and solution }
{ B (n, n) coefficient matrix becomes inverse }
{ Y (n) original constant vector }
{ W (n, m) constant vector(s) become solution vector }
{ DETERM is the determinant }
{ ERROR = 1 if singular }
{ INDEX (n, 3) }
{ NV is number of constant vectors }
*/
try
{
size_t ncol(b.size());
size_t irow, icol;
vector<vector<int> >index;
Matrix w;
Zeroise(w, ncol, ncol);
Zeroise(index, ncol, 3);
if(!GaussJordan2(b, y, w, index))
return false;
// Interchange columns
size_t m;
for (size_t 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(size_t k = 0; k < ncol; ++k)
Swap (b[k][irow], b[k][icol]);
}
}
for(size_t k = 0; k < ncol; ++k)
{
if(index [k][2] != 0)
{
throw "Error in GaussJordan: matrix is singular";
}
}
for(size_t i = 0; i < ncol; ++i)
coef[i] = w[i][0];
}
catch(char *b)
{
Application->MessageBox(b, "PolyFit Exception",
MB_OK | MB_ICONEXCLAMATION | MB_DEFBUTTON1 | MB_APPLMODAL | MB_SETFOREGROUND | MB_TOPMOST);
return false;
}
return true;
} // end; { procedure GaussJordan }
//----------------------------------------------------------------------------------------------
template <class FloatType>
bool __fastcall CPolyFit<FloatType>::GaussJordan2(Matrix &b,
const vector<FloatType> &y,
Matrix &w,
vector<vector<int> > &index)
{
//GaussJordan2; // first half of GaussJordan
// actual start of gaussj
try
{
FloatType big, t;
FloatType pivot;
FloatType determ;
size_t irow, icol;
size_t ncol(b.size());
size_t nv = 1; // single constant vector
for(size_t i = 0; i < ncol; ++i)
{
w[i][0] = y[i]; // copy constant vector
index[i][2] = -1;
}
determ = 1.0;
for(size_t i = 0; i < ncol; ++i)
{
// Search for largest element
big = 0.0;
for(size_t j = 0; j < ncol; ++j)
{
if(index[j][2] != 0)
{
for(size_t k = 0; k < ncol; ++k)
{
if(index[k][2] > 0)
throw "Error in GaussJordan2: matrix is singular";
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(size_t m = 0; m < ncol; ++m)
Swap (b [irow][m], b[icol][m]);
if (nv > 0)
for (size_t 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;
for(size_t m = 0; m < ncol; ++m)
b[icol][m] /= pivot;
if(nv > 0)
for(size_t m = 0; m < nv; ++m)
w[icol][m] /= pivot;
// Reduce nonpivot rows
for(size_t n = 0; n < ncol; ++n)
{
if(n != icol)
{
t = b[n][icol];
b[n][icol] = 0.0;
for(size_t m = 0; m < ncol; ++m)
b[n][m] -= b[icol][m] * t;
if(nv > 0)
for(size_t m = 0; m < nv; ++m)
w[n][m] -= w[icol][m] * t;
}
}
} // { i-loop }
}
catch(char *b)
{
Application->MessageBox(b, "PolyFit Exception",
MB_OK | MB_ICONEXCLAMATION | MB_DEFBUTTON1 | MB_APPLMODAL | MB_SETFOREGROUND | MB_TOPMOST);
return false;
}
return true;
//end; { GaussJordan2 }
}
//----------------------------------------------------------------------------------------------
// Least-squares fit to nrow sets of x and y pairs of points
template <class FloatType>
void __fastcall CPolyFit<FloatType>::LinFit
(const vector<FloatType> &x, // independent variable
vector<FloatType> &y, // dependent variable
vector<FloatType> &y_calc, // calculated dep. variable
vector<FloatType> &resid, // array of residuals
vector<FloatType> &coef, // coefficients
vector<FloatType> &sig, // error on coefficients
const size_t nrow, // length of variable array
size_t ncol) // number of terms
{
// bool error;
size_t i, j, nm;
FloatType xi, yi, yc, srs, see, sum_y, sum_y2;
Matrix xmatr; // Data matrix
Matrix a;
vector<FloatType> g; // Constant vector
Zeroise(g, ncol);
Zeroise(a, ncol, ncol);
Zeroise(xmatr, nrow, ncol);
for(i = 0; i < nrow; ++i)
{
// { setup x matrix }
xi = x[i];
xmatr[i][0] = 1.0; // { first column }
for(j = 1; j < ncol; ++j)
xmatr [i][j] = xmatr [i][j - 1] * xi;
}
Square (xmatr, y, a, g, nrow, ncol);
GaussJordan (a, g, coef);
sum_y = 0.0;
sum_y2 = 0.0;
srs = 0.0;
for(i = 0; i < nrow; ++i)
{
yi = y[i];
yc = 0.0;
for(j = 0; j < ncol; ++j)
yc += coef [j] * xmatr [i][j];
y_calc[i] = yc;
resid[i] = yc - yi;
srs += Sqr (resid [i]);
sum_y += yi;
sum_y2 += yi * yi;
}
if(nrow == ncol)
nm = 1;
else nm = nrow - ncol;
see = sqrt (srs / nm);
for(i = 0; i < ncol; ++i) // { errors on solution }
sig[i] = see * sqrt (a [i][i]);
}
//------------------------------------------------------------------------------------
// Utility functions
// fills a vector with zeros.
template <typename X>
void NSUtility::Zeroise(vector<X> &array, size_t n)
{
array.clear();
for(size_t j = 0; j < n; ++j)
array.push_back(0);
}
//--------------------------------------------------------------------------
// fills a (m by n) matrix with zeros.
template <typename X>
void NSUtility::Zeroise(vector<vector<X> > &matrix, size_t m, size_t n)
{
vector<X> zero;
Zeroise(zero, n);
matrix.clear();
for(size_t j = 0; j < m; ++j)
matrix.push_back(zero);
}
//--------------------------------------------------------------------------
// calcs x^i
template <typename X>
X NSUtility::Pow(const X &x, int i)
{
X ret = 1.0;
X t = i < 0 ? 1.0 / x : x;
int n = abs(i);
for(i = 0; i < n; ++i)
ret *= t;
return ret;
}
//-------------------------------------------------------------------------
One thing that may or may not be of note is that I had to move the precompile header include from the .cpp file over to the header file to get it to compile, or it had problems with the "using" statements at the top of the header file.
Appreciate the help! This is about to drive me nuts. I can't seem to locate any other good parabolic curve fitting code.
Regards,
Frank |
13-Jan-2006, 07:55
A Week in Kuwait by crystalattice
Similar Threads
