Like Tree3Likes
  • 1 Post By laserlight
  • 1 Post By laserlight
  • 1 Post By rcgldr

Passing C++ vectors to functions as arguments by reference and value

This is a discussion on Passing C++ vectors to functions as arguments by reference and value within the C++ Programming forums, part of the General Programming Boards category; Hello. I have a program that is working very well when I pass C++ vectors as arguments to my functions ...

  1. #1
    Registered User
    Join Date
    May 2013
    Location
    United States
    Posts
    58

    Passing C++ vectors to functions as arguments by reference and value

    Hello. I have a program that is working very well when I pass C++ vectors as arguments to my functions by reference, but I get some compilation errors when try to make a modification. I am also posting the entire program and its output below. so that you can see what is going on. I have commented out the line that causes an error.(I apologize for some of the indentation that got corrupted when I copied the code to the browser.)

    This program basically calculates the coefficients of a least square polynomial and then evaluates this polynomial at artificial data points and verifies that this actually reproduces the original data within reasonable floating point error.

    The function that computes the coefficients of the least square polynomial is
    Code:
     vector<double> LSPVecValued_GSL( const int,  const vector<float> &, const vector<float> &);
    and as you can see it returns a vector by value, and this vector contains the coefficients of the least square polynomial.

    There is also a function that evaluates this polynomial by accepting a vector argument by reference :
    Code:
    float evaluate_polynomial(double, vector<double>& ) ;
    I have also created another version of the evaluation function which accepts the same vector argument by value:
    Code:
    float evaluate_polynomial_ByValue(double t, vector<double> vec_a) ;
    In the program I call the first evaluation function (whose vector argument is passed by reference) by first using an intermediate vector variable containing the coefficients, and then I pass this vector as an argument to the evaluation function, as follows:
    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) directly by plugging in the function LSPVecValued_GSL"(...)" and this works without error, and this is a one step process, only one line of code is involved.

    However, I get a compilation error (line number 12 that I have commented out above) if I try to plug in the function "LSPVecValued_GSL(...)" into the first evaluation function that expects a vector argument by reference. I tried to put a "&" in front ofLSPVecValued_GSL but this did not fix the bug.

    What syntax is appropriate to use the first evaluation function (which accepts a vector argument by reference) if I want to plug in the vector-valued function LeastSquarePolynomial_GSL directly in the the first version of the evaluation function which expects a vector argument by reference?

    MANY THANKS FOR YOUR PRECIOUS TIME!!!



    HERE IS THE ENTIRE PROGRAM FOLLOWED BY ITS OUTPUT ( please see the line 96 (commented out) which causes the problem):
    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:
    > g++ -I/usr/local/include -L/usr/local/lib LSPExampleVecValued1_GSL.cpp -lgsl -lgslcblas -lm
    > ./a.out
    a[0] = 1
    a[1] = 0.999999
    a[2] = 1.00004
    a[3] = 0.999496
    a[4] = 1.00292
    a[5] = 0.99015
    a[6] = 1.02055
    a[7] = 0.972974
    a[8] = 1.0219
    a[9] = 0.98999
    a[10] = 1.00198
    --------------------
    x[0] = 0 ,y[0] = 1 , p(x[0]) ( EVALUATED FROM REFERENCE) = 1
    x[0] = 0 ,y[0] = 1 , p(x[0]) (EVALUATED FROM VALUE) = 1
    x[20] = 0.08 ,y[20] = 1.08696 , p(x[20]) ( EVALUATED FROM REFERENCE) = 1.08696
    x[20] = 0.08 ,y[20] = 1.08696 , p(x[20]) (EVALUATED FROM VALUE) = 1.08696
    x[40] = 0.16 ,y[40] = 1.19048 , p(x[40]) ( EVALUATED FROM REFERENCE) = 1.19048
    x[40] = 0.16 ,y[40] = 1.19048 , p(x[40]) (EVALUATED FROM VALUE) = 1.19048
    x[60] = 0.24 ,y[60] = 1.31579 , p(x[60]) ( EVALUATED FROM REFERENCE) = 1.31579
    x[60] = 0.24 ,y[60] = 1.31579 , p(x[60]) (EVALUATED FROM VALUE) = 1.31579
    x[80] = 0.32 ,y[80] = 1.47058 , p(x[80]) ( EVALUATED FROM REFERENCE) = 1.47058
    x[80] = 0.32 ,y[80] = 1.47058 , p(x[80]) (EVALUATED FROM VALUE) = 1.47058
    x[100] = 0.4 ,y[100] = 1.6666 , p(x[100]) ( EVALUATED FROM REFERENCE) = 1.6666
    x[100] = 0.4 ,y[100] = 1.6666 , p(x[100]) (EVALUATED FROM VALUE) = 1.6666
    x[120] = 0.48 ,y[120] = 1.92248 , p(x[120]) ( EVALUATED FROM REFERENCE) = 1.92248
    x[120] = 0.48 ,y[120] = 1.92248 , p(x[120]) (EVALUATED FROM VALUE) = 1.92248
    x[140] = 0.56 ,y[140] = 2.26887 , p(x[140]) ( EVALUATED FROM REFERENCE) = 2.26887
    x[140] = 0.56 ,y[140] = 2.26887 , p(x[140]) (EVALUATED FROM VALUE) = 2.26887
    x[160] = 0.64 ,y[160] = 2.75728 , p(x[160]) ( EVALUATED FROM REFERENCE) = 2.75728
    x[160] = 0.64 ,y[160] = 2.75728 , p(x[160]) (EVALUATED FROM VALUE) = 2.75728
    x[180] = 0.72 ,y[180] = 3.47516 , p(x[180]) ( EVALUATED FROM REFERENCE) = 3.47516
    x[180] = 0.72 ,y[180] = 3.47516 , p(x[180]) (EVALUATED FROM VALUE) = 3.47516
    x[200] = 0.8 ,y[200] = 4.5705 , p(x[200]) ( EVALUATED FROM REFERENCE) = 4.5705
    x[200] = 0.8 ,y[200] = 4.5705 , p(x[200]) (EVALUATED FROM VALUE) = 4.5705
    x[220] = 0.88 ,y[220] = 6.29099 , p(x[220]) ( EVALUATED FROM REFERENCE) = 6.29099
    x[220] = 0.88 ,y[220] = 6.29099 , p(x[220]) (EVALUATED FROM VALUE) = 6.29099
    x[240] = 0.96 ,y[240] = 9.04402 , p(x[240]) ( EVALUATED FROM REFERENCE) = 9.04402
    x[240] = 0.96 ,y[240] = 9.04402 , p(x[240]) (EVALUATED FROM VALUE) = 9.04402
    >
    Last edited by FortranLevelC++; 03-12-2014 at 10:31 PM.

  2. #2
    C++ Witch laserlight's Avatar
    Join Date
    Oct 2003
    Location
    Singapore
    Posts
    21,720
    Quote Originally Posted by FortranLevelC++
    What syntax is appropriate to use the first evaluation function (which accepts a vector argument by reference) if I want to plug in the vector-valued function LeastSquarePolynomial_GSL directly in the the first version of the evaluation function which expects a vector argument by reference?
    Well, since you passed the vector by (non-const) lvalue reference, you are presumably modifying the vector or its contents within evaluate_polynomial. Therefore, you should create the vector first, and then pass it as an argument:
    Code:
    vector<double> temp = LSPVecValued_GSL(deg, vec_x, vec_y);
    cout << "x[" << j << "] = " << vec_x[j] << " ,y[" << j << "] = " << vec_y[j]
         << " , p(x[" << j << "]) (EVALUATED FROM REFERENCE) = "
         << evaluate_polynomial(vec_x[j], temp) << endl;
    // do something with temp
    // ...
    If it turns out that you are not modifying the vector or its contents from within evaluate_polynomial, then you should change the parameter to be a const reference (if you do not make a copy of the vector). You will then be able to call evaluate_polynomial using the same syntax as if you used pass by value for that parameter.
    C + C++ Compiler: MinGW port of GCC
    Version Control System: Bazaar

    Look up a C++ Reference and learn How To Ask Questions The Smart Way

  3. #3
    Registered User
    Join Date
    May 2013
    Location
    United States
    Posts
    58
    Quote Originally Posted by laserlight View Post
    Well, since you passed the vector by (non-const) lvalue reference, you are presumably modifying the vector or its contents within evaluate_polynomial. Therefore, you should create the vector first, and then pass it as an argument:
    Code:
    vector<double> temp = LSPVecValued_GSL(deg, vec_x, vec_y);
    cout << "x[" << j << "] = " << vec_x[j] << " ,y[" << j << "] = " << vec_y[j]
         << " , p(x[" << j << "]) (EVALUATED FROM REFERENCE) = "
         << evaluate_polynomial(vec_x[j], temp) << endl;
    // do something with temp
    // ...
    If it turns out that you are not modifying the vector or its contents from within evaluate_polynomial, then you should change the parameter to be a const reference (if you do not make a copy of the vector). You will then be able to call evaluate_polynomial using the same syntax as if you used pass by value for that parameter.
    LaserLight,

    Many thanks! Actually, I was not intending to modify the arguments anyway. I have made the arguments of the function constants as you recommended, and this time it is working. However, this time I don't understand why the compiler is not giving an error message when the arguments are constants. Why would the error message go away when I make the argument constant?

    After all, I still cannot understand the difference between the two step evaluation like this:

    Code:
    vec_a = LSPVecValued_GSL(  deg,  vec_x, vec_y);
    cout<<  evaluate_polynomial(vec_x[j], vec_a)  << endl;  // which compiles and works without error
    and the one line version like this:

    Code:
    cout <<  evaluate_polynomial(vec_x[j],  LSPVecValued_GSL(  deg,  vec_x , vec_y))  << endl;  // compiler gives an error message
    I imagined that the only difference should be that within the for loop calling evaluate_function many times will be computationally inefficient, but this was just a way of testing the syntax of my code at this stage. Otherwise, I understand that if the function tries to modify constant arguments then there will be an error.

    Anyway, here is the program and its output with both direct and indirect function calls, and all versions are working:

    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( const double, const vector<double>& ) ;
    float evaluate_polynomial_ByValue(double , vector<double> ) ;
    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;
    
    
       
       
       {
         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 INDIRECTLY) = "
               <<  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 DIRECTLY) = "
               <<  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
    float evaluate_polynomial( const double t, const 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;
    
    
    }

    The output of the program:
    >
    > g++ -I/usr/local/include -L/usr/local/lib LSPExampleVecValued1_GSL.cpp -lgsl -lgslcblas -lm
    > ./a.out
    a[0] = 1
    a[1] = 0.999999
    a[2] = 1.00004
    a[3] = 0.999496
    a[4] = 1.00292
    a[5] = 0.99015
    a[6] = 1.02055
    a[7] = 0.972974
    a[8] = 1.0219
    a[9] = 0.98999
    a[10] = 1.00198
    --------------------
    x[0] = 0 ,y[0] = 1 , p(x[0]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1
    x[0] = 0 ,y[0] = 1 , p(x[0]) (EVALUATED FROM VALUE) = 1
    x[0] = 0 ,y[0] = 1 , p(x[0]) (EVALUATED FROM REFERENCE DIRECTLY) = 1
    x[20] = 0.08 ,y[20] = 1.08696 , p(x[20]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.08696
    x[20] = 0.08 ,y[20] = 1.08696 , p(x[20]) (EVALUATED FROM VALUE) = 1.08696
    x[20] = 0.08 ,y[20] = 1.08696 , p(x[20]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.08696
    x[40] = 0.16 ,y[40] = 1.19048 , p(x[40]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.19048
    x[40] = 0.16 ,y[40] = 1.19048 , p(x[40]) (EVALUATED FROM VALUE) = 1.19048
    x[40] = 0.16 ,y[40] = 1.19048 , p(x[40]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.19048
    x[60] = 0.24 ,y[60] = 1.31579 , p(x[60]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.31579
    x[60] = 0.24 ,y[60] = 1.31579 , p(x[60]) (EVALUATED FROM VALUE) = 1.31579
    x[60] = 0.24 ,y[60] = 1.31579 , p(x[60]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.31579
    x[80] = 0.32 ,y[80] = 1.47058 , p(x[80]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.47058
    x[80] = 0.32 ,y[80] = 1.47058 , p(x[80]) (EVALUATED FROM VALUE) = 1.47058
    x[80] = 0.32 ,y[80] = 1.47058 , p(x[80]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.47058
    x[100] = 0.4 ,y[100] = 1.6666 , p(x[100]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.6666
    x[100] = 0.4 ,y[100] = 1.6666 , p(x[100]) (EVALUATED FROM VALUE) = 1.6666
    x[100] = 0.4 ,y[100] = 1.6666 , p(x[100]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.6666
    x[120] = 0.48 ,y[120] = 1.92248 , p(x[120]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 1.92248
    x[120] = 0.48 ,y[120] = 1.92248 , p(x[120]) (EVALUATED FROM VALUE) = 1.92248
    x[120] = 0.48 ,y[120] = 1.92248 , p(x[120]) (EVALUATED FROM REFERENCE DIRECTLY) = 1.92248
    x[140] = 0.56 ,y[140] = 2.26887 , p(x[140]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 2.26887
    x[140] = 0.56 ,y[140] = 2.26887 , p(x[140]) (EVALUATED FROM VALUE) = 2.26887
    x[140] = 0.56 ,y[140] = 2.26887 , p(x[140]) (EVALUATED FROM REFERENCE DIRECTLY) = 2.26887
    x[160] = 0.64 ,y[160] = 2.75728 , p(x[160]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 2.75728
    x[160] = 0.64 ,y[160] = 2.75728 , p(x[160]) (EVALUATED FROM VALUE) = 2.75728
    x[160] = 0.64 ,y[160] = 2.75728 , p(x[160]) (EVALUATED FROM REFERENCE DIRECTLY) = 2.75728
    x[180] = 0.72 ,y[180] = 3.47516 , p(x[180]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 3.47516
    x[180] = 0.72 ,y[180] = 3.47516 , p(x[180]) (EVALUATED FROM VALUE) = 3.47516
    x[180] = 0.72 ,y[180] = 3.47516 , p(x[180]) (EVALUATED FROM REFERENCE DIRECTLY) = 3.47516
    x[200] = 0.8 ,y[200] = 4.5705 , p(x[200]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 4.5705
    x[200] = 0.8 ,y[200] = 4.5705 , p(x[200]) (EVALUATED FROM VALUE) = 4.5705
    x[200] = 0.8 ,y[200] = 4.5705 , p(x[200]) (EVALUATED FROM REFERENCE DIRECTLY) = 4.5705
    x[220] = 0.88 ,y[220] = 6.29099 , p(x[220]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 6.29099
    x[220] = 0.88 ,y[220] = 6.29099 , p(x[220]) (EVALUATED FROM VALUE) = 6.29099
    x[220] = 0.88 ,y[220] = 6.29099 , p(x[220]) (EVALUATED FROM REFERENCE DIRECTLY) = 6.29099
    x[240] = 0.96 ,y[240] = 9.04402 , p(x[240]) ( EVALUATED FROM REFERENCE INDIRECTLY) = 9.04402
    x[240] = 0.96 ,y[240] = 9.04402 , p(x[240]) (EVALUATED FROM VALUE) = 9.04402
    x[240] = 0.96 ,y[240] = 9.04402 , p(x[240]) (EVALUATED FROM REFERENCE DIRECTLY) = 9.04402
    >
    Last edited by FortranLevelC++; 03-13-2014 at 01:00 AM.

  4. #4
    C++ Witch laserlight's Avatar
    Join Date
    Oct 2003
    Location
    Singapore
    Posts
    21,720
    Quote Originally Posted by FortranLevelC++
    However, this time I don't understand why the compiler is not giving an error message when the arguments are constants. Why would the error message go away when I make the argument constant?
    A temporary, like that returned by LSPVecValued_GSL, cannot be bound to a non-const lvalue reference. However, it can be bound to a const lvalue reference, in which case its lifetime is extended to be the lifetime of the reference.

    By the way, I suggest writing evaluate_polynomial as something like:
    Code:
    float evaluate_polynomial(const double t, const vector<double>& vec_a)
    {
        double p = 0.0; // accumulates in double precision
        const vector<double>::size_type vec_a_size = vec_a.size();
        for (vector<double>::size_type k = 0; k < vec_a_size; ++k)
        {
            p += vec_a[k] * pow(t, k);
        }
        return static_cast<float>(p);
    }
    C + C++ Compiler: MinGW port of GCC
    Version Control System: Bazaar

    Look up a C++ Reference and learn How To Ask Questions The Smart Way

  5. #5
    Registered User
    Join Date
    May 2013
    Location
    United States
    Posts
    58
    Many thanks!

  6. #6
    Registered User
    Join Date
    Jun 2005
    Posts
    6,294
    Quote Originally Posted by laserlight View Post
    By the way, I suggest writing evaluate_polynomial as something like:
    Code:
    float evaluate_polynomial(const double t, const vector<double>& vec_a)
    {
        double p = 0.0; // accumulates in double precision
        const vector<double>::size_type vec_a_size = vec_a.size();
        for (vector<double>::size_type k = 0; k < vec_a_size; ++k)
        {
            p += vec_a[k] * pow(t, k);
        }
        return static_cast<float>(p);
    }
    Personally, I'd specify the return type as double, write a version that uses an templated iterator, avoid the use of pow(), myself. Something like
    Code:
    template<class Iter> double evaluate_polynomial(double t, Iter begin, Iter end)
    {
          double result = 0.0, tpower = 1.0;
          while (begin != end)
          {
                result += *begin * tpower;
                tpower *= t;
                ++begin;
          }
          return result;
    }
    
    double evaluate_polynomial(double t, const std::vector<double> & vec_a)
    {
        return evaluate_polynomial(t, vec_a.begin(), vec_a.end());
    }
    In C++11, I'd consider using a range-based for loop, as an alternative to template machinery (I'm still of two minds which is better - both have strengths and weaknesses).
    Right 98% of the time, and don't care about the other 3%.

  7. #7
    Master Apprentice phantomotap's Avatar
    Join Date
    Jan 2008
    Posts
    4,258
    Personally, I'd specify the return type as double, write a version that uses an templated iterator, avoid the use of pow(), myself.
    O_o

    If you are going that route, you should use templates and traits to determine the result and intermediate representation.

    With C++11, the mechanic should be an accumulator.

    Soma
    “Often out of periods of losing come the greatest strivings toward a new winning streak.” -- Fred Rogers
    “Salem Was Wrong!” -- Pedant Necromancer

  8. #8
    Registered User
    Join Date
    Jun 2005
    Posts
    6,294
    Quote Originally Posted by phantomotap View Post
    If you are going that route, you should use templates and traits to determine the result and intermediate representation.
    Yes, indeed. I wasn't going to write such a thing on the fly in a forum though.

    Quote Originally Posted by phantomotap View Post
    With C++11, the mechanic should be an accumulator.
    I would use the word "can" rather than "should". I'm not disagreeing with you, but I consider it best to make a decision on approach rather than blindly follow such a rule.
    Right 98% of the time, and don't care about the other 3%.

  9. #9
    Master Apprentice phantomotap's Avatar
    Join Date
    Jan 2008
    Posts
    4,258
    O_o

    Indeed. I was referring to this "flavor" of algorithm not attempting to make a rule.

    Soma
    “Often out of periods of losing come the greatest strivings toward a new winning streak.” -- Fred Rogers
    “Salem Was Wrong!” -- Pedant Necromancer

  10. #10
    Registered User
    Join Date
    Apr 2013
    Posts
    1,283
    If interested, here is a link to a document (.rtf) of an alternate method of doing a polynomial least squares fit that uses orthogonal polynomials, eliminating the need to invert a matrix (this improves the evaluation of higher order polynomials). It includes example C code.

    http://rcgldr.net/misc/opls.rtf
    Last edited by rcgldr; 03-13-2014 at 09:42 AM.

  11. #11
    Registered User
    Join Date
    Apr 2013
    Posts
    1,283
    Quote Originally Posted by rcgldr View Post
    If interested, here is a link to a document (.rtf) of an alternate method of doing a polynomial least squares fit that uses orthogonal polynomials, It includes example C code.
    The example code is really psuedo code, using an index of [-1][] on one of the variables. Here is example code, it's old. Note that f[] is the y values, where y = f(x).

    Code:
    static int m, n;                /* input values */
    static double x[nmax];
    static double f[nmax];
    static double w[nmax];
    static double b[mmax];          /* generated values */
    static double A[mmax];
    static double B[mmax];
    static double L[mmax];
    static double W[mmax];
    static double p2[nmax];
    static double p1[nmax];
    static double p0[nmax];
    static double c[mmax];          /* coefficients for y(x) */
    static double z[nmax];          /* calculated f[] */
    static double xi, zi, vr;
    static double D0, D1;           /* constants */
    
    
    static double *pp2;             /* pointers to logical p2, p1, p0 */
    static double *pp1;
    static double *pp0;
    static double *ppx;
    
    
    /*------------------------------------------------------*/
    /*      polyf           poly fit                        */
    /*      in:             x[], f[], w[], n, m             */
    /*      out:            b[], A[], B[], L[], W[]         */
    /*------------------------------------------------------*/
    static void polyf()
    {
    register int i, j;
    
        D0 = (double)0.;            /* init */
        D1 = (double)1.;
        pp2 = p2;
        pp1 = p1;
        pp0 = p0;
    
        j = 0;
        A[j] = D0;                  /* calc A, p[j], p[j-1], L, W */
        L[j] = D0;                  /* note A[0] not used */
        W[j] = D0;
        for(i = 0; i < n; i++){
            pp0[i] = D1;
            pp1[i] = D0;
            L[j] += w[i];
            W[j] += w[i]*f[i];}
        B[0] = D0;
        b[j] = W[j]/L[j];
    
        for(j = 1; j <= m; j++){
            ppx = pp2;              /* save old p[j], p[j-1] */
            pp2 = pp1;
            pp1 = pp0;
            pp0 = ppx;
            A[j] = D0;              /* calc A */
            for(i = 0; i < n; i++){
                A[j] += w[i]*x[i]*pp1[i]*pp1[i]/L[j-1];}
            L[j] = D0;              /* calc p[j], L, W */
            W[j] = D0;
            for(i = 0; i < n; i++){
                pp0[i] = (x[i]-A[j])*pp1[i]-B[j-1]*pp2[i];
                L[j] += w[i]*pp0[i]*pp0[i];
                W[j] += w[i]*f[i]*pp0[i];}
            B[j] = L[j]/L[j-1];     /* calc B[], b[] */
            b[j] = W[j]/L[j];}
    }
    
    /*------------------------------------------------------*/
    /*      gcoef   generate coefficients                   */
    /*      in:     b[], A[], B[]                           */
    /*      out:    c[]                                     */
    /*      uses:   p0[], p1[], p2[]                        */
    /*------------------------------------------------------*/
    static void gcoef()
    {
    register int i, j;
        for(i = 0; i <= m; i++){            /* init */
            c[i] = p2[i] = p1[i] = p0[i] = 0.;}
        p0[0] = D1;
        c[0] += b[0]*p0[0];
        for(j = 1; j <= m; j++){            /* generate coefs */
            p2[0] = p1[0];
            p1[0] = p0[0];
            p0[0] = -A[j]*p1[0]-B[j-1]*p2[0];
            c[0] += b[j]*p0[0];
            for(i = 1; i <= j; i++){
                p2[i] = p1[i];
                p1[i] = p0[i];
                p0[i] = p1[i-1]-A[j]*p1[i]-B[j-1]*p2[i];
                c[i] += b[j]*p0[i];}}
    }
    Last edited by rcgldr; 03-13-2014 at 10:58 AM.

  12. #12
    Registered User
    Join Date
    Apr 2006
    Posts
    2,032
    Quote Originally Posted by FortranLevelC++ View Post
    However, this time I don't understand why the compiler is not giving an error message when the arguments are constants. Why would the error message go away when I make the argument constant?
    Laserlight answered this partially, but there is a larger reason. Logically a parameter passed by non-const reference is an output parameter, who's later state reflects the result of the function. A temporary cannot be accessed after it's passed to the function, so it cannot be a good output parameter. A temporary can be a const ref or value parameter, serving as input. In C++11 it can also be a rvalue reference, in which case you're explicitly consuming the value, and that object is not intended to be accessed after the call.

    As a quirk, a temporary can be a non-const this parameter to functions and operators.
    It is too clear and so it is hard to see.
    A dunce once searched for fire with a lighted lantern.
    Had he known what fire was,
    He could have cooked his rice much sooner.

  13. #13
    Registered User
    Join Date
    May 2013
    Location
    United States
    Posts
    58
    Quote Originally Posted by rcgldr View Post
    If interested, here is a link to a document (.rtf) of an alternate method of doing a polynomial least squares fit that uses orthogonal polynomials, eliminating the need to invert a matrix (this improves the evaluation of higher order polynomials). It includes example C code.

    http://rcgldr.net/misc/opls.rtf

    Many thanks. I will investigate this under the microscope.

  14. #14
    Registered User
    Join Date
    Jun 2005
    Posts
    6,294
    Quote Originally Posted by King Mir View Post
    A temporary cannot be accessed after it's passed to the function, so it cannot be a good output parameter.
    Philosophically (albeit not in C++) a temporary can be a good output parameter if the caller wishes to discard that output.
    Right 98% of the time, and don't care about the other 3%.

  15. #15
    Registered User
    Join Date
    Apr 2006
    Posts
    2,032
    Quote Originally Posted by grumpy View Post
    Philosophically (albeit not in C++) a temporary can be a good output parameter if the caller wishes to discard that output.
    Arguably, yes. The counterpoint would be if when encountering a temporary used that way, is it more likely to be a garbage value or programmer error. If the latter, then more verbose ways of discarding output are better.
    It is too clear and so it is hard to see.
    A dunce once searched for fire with a lighted lantern.
    Had he known what fire was,
    He could have cooked his rice much sooner.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. int&, passing of arguments by reference, etc
    By jackson6612 in forum C++ Programming
    Replies: 6
    Last Post: 05-17-2011, 03:26 AM
  2. Replies: 10
    Last Post: 05-15-2011, 06:52 PM
  3. passing vectors to functions
    By igltorry in forum C Programming
    Replies: 5
    Last Post: 03-13-2011, 12:54 AM
  4. Managed C++, passing arguments by reference
    By jimzy in forum C++ Programming
    Replies: 8
    Last Post: 11-02-2007, 01:03 PM
  5. passing arguments through functions
    By cwd in forum C Programming
    Replies: 2
    Last Post: 09-30-2001, 06:07 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21