source: trunk/src/vector.cc @ 36

// $Id: vector.cc 36 2004-02-12 10:29:21Z peter $

// #include <iostream>
#include <string>
#include <sstream>
#include <vector>
#include "vector.h"


using namespace thep_gsl_api;


// Constructors and Destructors
///////////////////////////////
vector::vector() : v_( NULL )
{
  is_col_ = true;
}


vector::vector( const size_t& N, bool is_col_vector, 
    bool init_to_zero ) :
  is_col_( is_col_vector )
{
  if( init_to_zero )
    { 
     v_ = gsl_vector_calloc( N );
    }
  else     
    {
     v_ = gsl_vector_alloc ( N );
    }
}


// Is this the way to do it? No copy ...
// internal data ...
vector::vector( gsl_vector* v, bool is_col ) : v_( v ), is_col_( is_col )
{
}


// Copy constructor
vector::vector( const vector& other )
{
  v_ = new_copy( other.get_gsl_vector() );
  is_col_ = other.is_column();
}


// Constructor that gets data from istream
vector::vector(std::istream& is)
{
  using namespace std;
  
  // read the data file and store in stl vectors (dynamically expandable)
  std::vector<std::vector<double> > data_matrix;
  u_int nof_columns=0;
  u_int nof_rows = 0;
  string s;
  for (nof_rows = 0; getline(is, s, '\n'); nof_rows++) {
   istringstream line(s);
   std::vector<double> v;
   string tmp_string;
   while (line >> tmp_string) {
    if(!is.good()) {
     cerr << "vector::vector(std::istream&): "
    << "error reading data file!" << endl;
     exit(1);
    }
    double t=atof(tmp_string.c_str());
    // Here we should check that it actually was a double!!!!!
    v.push_back(t);
   }
   if(nof_columns==0)
     nof_columns=v.size();
   else if(v.size()!=nof_columns) {
    cerr << "vector::vector(std::istream&) data file error: "
   << "line" << nof_rows+1 << " has " << v.size() 
   << " columns; expected " << nof_columns
   << " columns"
   << endl; 
    exit(1);
   }
   data_matrix.push_back(v);
  }
  
  // manipulate the state of the stream to be good
  is.clear(std::ios::goodbit);
  
  // convert the data to a gsl vector and check that data file is a
  // column vector or a row vector
  if(nof_columns==1) {
   v_ = gsl_vector_alloc ( nof_rows );
   is_col_ = true;
   for(u_int i=0;i<nof_rows;i++) 
     gsl_vector_set( v_, i, data_matrix[i][0] );
  }
  else if(nof_rows==1){
   v_ = gsl_vector_alloc ( nof_columns );
   is_col_ = false;
   for(u_int i=0;i<nof_columns;i++) 
     gsl_vector_set( v_, i, data_matrix[0][i] );
  }
  else {
    cerr << "vector::vector(std::istream&) data file error: "
   << "file has " << nof_rows << " rows and " << nof_columns 
   << " columns; expected a row vector or a column vector"
   << endl;
    exit(1);
  }
}

  
vector::~vector()
{
  if( v_ ) 
    {
     gsl_vector_free( v_ );
     v_ = NULL;
    }
}


gsl_vector* vector::new_copy( const gsl_vector* p_other )
{
  // Get dimenstions
  size_t N = p_other->size;
    
  // Create new empty vector
  gsl_vector* p_res =  gsl_vector_alloc( N );

  // Copy p_others elements into p_res 
  gsl_vector_memcpy( p_res, p_other );

  return p_res;
}


// Operators on vectors
////////////////////////
vector& vector::operator=( const vector& other )
{
  if( this != &other ) 
    {
     gsl_vector* v_new = new_copy( other.get_gsl_vector() );
     if( v_ ) gsl_vector_free( v_ );
     v_ = v_new;
    }
  return *this;
} 


vector vector::operator+( const vector &other ) const
{
  assert( size() == other.size() );
  vector res( *this );
  gsl_vector_add( res.get_gsl_vector(), other.get_gsl_vector() );
  return res;
}


vector vector::operator-( const vector &other ) const
{
  assert( size() == other.size() );
  vector res( *this );
  gsl_vector_sub( res.get_gsl_vector(), other.get_gsl_vector() );
  return res;
}

vector vector::operator-() const
{
  vector res( *this );
  gsl_vector_scale (res.get_gsl_vector() , -1.);
  return res;
}
double vector::operator*( const vector &other ) const
{
  assert( size() == other.size() );
  double res = 0.0;;
  for ( size_t i = 0; i < other.size(); ++i ) 
    {
     res += other.get( i ) * get( i );
    }
  return res;
}


int vector::operator/=( const vector& other ) const
{
  return gsl_vector_div( v_, other.get_gsl_vector() );
}

// Each element is divided by the "constant" 
int vector::operator/=( const double& constant ) 
{
  assert( constant != 0 );
  return gsl_vector_scale( v_, 1 / constant );
}

// This is not the dot product, hence a matrix
// is returned
// matrix vector::operator*( const vector &other ) const
// {
//   assert( rows() == other.cols() );
//   matrix res( rows(), other.cols() );
//   gsl_linalg_matmult( v_, other.get_gsl_vector(), 
//          res.get_gsl_vector() );
//   return res;
// }



vector vector::mul_elements( const vector& other ) const
{
  assert( size() == other.size() );
  vector res( *this );
  gsl_vector_mul( res.get_gsl_vector(), other.get_gsl_vector() );
  return res;
}

std::ostream& thep_gsl_api::operator<< ( std::ostream& s_out, 
           const vector& a )
{
  using namespace std;
  s_out.setf( ios::fixed );  

  if( a.is_column() )
    {
     for ( size_t i = 0; i < a.size(); ++i ) 
       {
  s_out << a.get( i  ) << endl;
       }
    }
  else
    {
     for ( size_t j = 0; j < a.size(); ++j ) 
       {
  s_out << a.get( j ) << " ";
       }
    }

  return s_out;
}

double vector::sum() const
{
  double sum = 0;
  for (int i=0; i<size(); i++)
   sum += gsl_vector_get( v_, i );
  return( sum );
}  
 
 

// // Public functions (A-Z)
// /////////////////////////