Code:
vec_a = LSPVecValued_GSL( deg, vec_x , vec_y);
for(int j=0; j< n ; j=j+20 )
{
cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) ( EVALUATED FROM REFERENCE) = "
<< evaluate_polynomial( vec_x[j], vec_a) << endl; // This version works without error
cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) (EVALUATED FROM VALUE) = "
<< evaluate_polynomial_ByValue( vec_x[j], LSPVecValued_GSL( deg, vec_x , vec_y) ) << endl; // This also works w/o error
/* cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) (EVALUATED FROM REFERENCE) = "
<< evaluate_polynomial( vec_x[j], LSPVecValued_GSL( deg, vec_x , vec_y)) << endl; */ // But this does not compile!!!
}
As you can see above, I am also able to call the second evaluation function (the one whose vector argument is passed by value)
Code:
/// COMPILATION SYNTAX: g++ -I/usr/local/include -L/usr/local/lib example4.cpp -lgsl -lgslcblas -lm
/// Here the necessary libraries are libgsl, libgslcblas, and libm
//
// One issue is that the double precision is very important for matrix LS computation, whereas the original data is float.
// So be careful about casting the data from float to double.
#include <iostream>
#include <stdio.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_math.h>
//#include<math>
#include<vector>
using namespace std;
float evaluate_polynomial(double, vector<double>& ) ;
float evaluate_polynomial_ByValue(double t, vector<double> vec_a) ;
void LeastSquarePolynomial_GSL( const int, const vector<float> &, const vector<float> &, vector<double> &);
vector<double> LSPVecValued_GSL( const int, const vector<float> &, const vector<float> &);
int main()
{
int i, n, deg, dim_a ;
double xi, yi;
float temp;
vector<double> vec_a;
//vector<double> vec_y;
//vector<double> vec_x;
vector<float> vec_y;
vector<float> vec_x;
deg = 10 ;
n = 250 ;
for(int j=0; j<n; j++)
{
//yi = 1 + j + j*j; // y[j]=polynomial p(x)=1+x+x^2 evaluated at x[i]=i
yi=0;
//xi = (double)j;
xi= ((double)j)/n; // rescale to prevent catastrophic cancellation
//xi = ((float)j)/(5.0 +j);
// temp = (double )xi;
// xi = 5.0*2*3.141592*j/n;
temp = (float)xi ;
vec_x.push_back(temp);
for(int k=0; k<=deg; k++)
{
yi= yi+ pow(xi,k);
// yi = cos(xi) ;
}
temp=(float)yi;
vec_y.push_back(temp);
}
// LeastSquarePolynomial_GSL( deg, vec_x , vec_y, vec_a );
vec_a = LSPVecValued_GSL( deg, vec_x , vec_y);
for(int k=0; k <= deg; k++)
{
// cout << "a["<<k<<"] = " << vec_a[k] << endl;
cout << "a["<<k<<"] = " << LSPVecValued_GSL( deg, vec_x , vec_y)[k] << endl; // k'th component of output vector of function
}
cout << "--------------------" << endl;
// vec_a = LSPVecValued_GSL( deg, vec_x , vec_y);
{
for(int j=0; j< n ; j=j+20 )
{
cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) ( EVALUATED FROM REFERENCE) = "
<< evaluate_polynomial( vec_x[j], vec_a) << endl; // This version works without error
cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) (EVALUATED FROM VALUE) = "
<< evaluate_polynomial_ByValue( vec_x[j], LSPVecValued_GSL( deg, vec_x , vec_y) ) << endl; // This also works w/o error
/* cout<<"x["<<j<<"] = " << vec_x[j] << " ,y["<<j<<"] = " << vec_y[j] <<" , p(x["<<j<<"]) (EVALUATED FROM REFERENCE) = "
<< evaluate_polynomial( vec_x[j], LSPVecValued_GSL( deg, vec_x , vec_y)) << endl; */ // But this does not compile!!!
}
}
vec_a.clear();
vec_x.clear();
vec_y.clear();
}
void LeastSquarePolynomial_GSL( const int deg, const vector<float>& vec_x , const vector<float>& vec_y, vector<double>& vec_a )
{ // this version takes float arguments for x and y values, but returns a double vec_a
//When using this float version instead of the original double precision version, numerical accuracy is lost, and therefore
// the maximum recommended polynomial degree = 10, Maximum number of data points = 500.
// (Always re-scaling 0.0 <= xi = xi/N <= 1.0)
int i, n, dim_a ;
double xi, yi, chisq;
gsl_matrix *X, *cov;
gsl_vector *y, *c;
float tempFloat;
n =vec_x.size(); // Of course vec_x.size() equals vec_y.size()
dim_a = deg+1;
double a[dim_a];
X = gsl_matrix_alloc(n,dim_a);
y = gsl_vector_alloc(n);
c = gsl_vector_alloc(dim_a);
cov = gsl_matrix_alloc(dim_a,dim_a);
for(int j=0; j<n; j++)
{
yi= (double)vec_y[j]; // cast into double the original floats..
xi= (double)vec_x[j];
gsl_vector_set(y, j, yi );
for(int k=0; k<= deg; k++)
{
gsl_matrix_set(X, j, k, pow( xi , k) );
}
}
{
gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n,dim_a);
gsl_multifit_linear(X, y, c, cov, &chisq, work);
gsl_multifit_linear_free(work);
}
for(int k=0; k< dim_a; k++)
{
a[k] = gsl_vector_get(c, (k) );
vec_a.push_back( a[k] );
}
gsl_matrix_free(X);
gsl_vector_free(y);
gsl_vector_free(c);
gsl_matrix_free(cov);
}
float evaluate_polynomial(double t, vector<double>& vec_a) // takes double arguments but returns a float value
{
double p=0.0; // accumulates in double precision
float temp;
int deg = vec_a.size() -1;
for(int k=0; k<=deg; k++)
{
p= p + vec_a[k]*pow(t,k) ;
// p= p + k*k*pow(t,k) + k ;
//cout << "Inside polynomial evaluation function vec_a["<<k<<"] = " << vec_a[k] << endl;
}
temp = (float)p;
return temp;
}
float evaluate_polynomial_ByValue(double t, vector<double> vec_a) // takes double arguments but returns a float value
{
double p=0.0; // accumulates in double precision
float temp;
int deg = vec_a.size() -1;
for(int k=0; k<=deg; k++)
{
p= p + vec_a[k]*pow(t,k) ;
// p= p + k*k*pow(t,k) + k ;
//cout << "Inside polynomial evaluation function vec_a["<<k<<"] = " << vec_a[k] << endl;
}
temp = (float)p;
return temp;
}
vector<double> LSPVecValued_GSL( const int deg, const vector<float>& vec_x , const vector<float>& vec_y)
{
vector<double> vec_a ; // note that this time vec_a is not reference but value.
// this version takes float arguments for x and y values, but returns a double vec_a
//When using this float version instead of the original double precision version, numerical accuracy is lost, and therefore
// the maximum recommended polynomial degree = 10, Maximum number of data points = 500.
// (Always re-scaling 0.0 <= xi = xi/N <= 1.0)
int i, n, dim_a ;
double xi, yi, chisq;
gsl_matrix *X, *cov;
gsl_vector *y, *c;
float tempFloat;
n =vec_x.size(); // Of course vec_x.size() equals vec_y.size()
dim_a = deg+1;
double a[dim_a];
X = gsl_matrix_alloc(n,dim_a);
y = gsl_vector_alloc(n);
c = gsl_vector_alloc(dim_a);
cov = gsl_matrix_alloc(dim_a,dim_a);
for(int j=0; j<n; j++)
{
yi= (double)vec_y[j]; // cast into double the original floats..
xi= (double)vec_x[j];
gsl_vector_set(y, j, yi );
for(int k=0; k<= deg; k++)
{
gsl_matrix_set(X, j, k, pow( xi , k) );
}
}
{
gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(n,dim_a);
gsl_multifit_linear(X, y, c, cov, &chisq, work);
gsl_multifit_linear_free(work);
}
for(int k=0; k< dim_a; k++)
{
a[k] = gsl_vector_get(c, (k) );
vec_a.push_back( a[k] );
}
gsl_matrix_free(X);
gsl_vector_free(y);
gsl_vector_free(c);
gsl_matrix_free(cov);
return vec_a;
}
Here is the output of the program, which verifies that when the polynomial is evaluated, it does reproduce the original test data within reasonable roundoff error:
Originally Posted by FortranLevelC++
