source: trunk/src/vector.cc @ 132

// $Id: vector.cc 132 2004-08-09 13:24:38Z peter $

#include "vector.h"
#include "random_singleton.h"

#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <utility>





namespace theplu {
namespace gslapi {



  vector::vector(void)
    : v_(NULL)
  {
  }



  vector::vector(const size_t n, const double init_value)
  {
    v_ = gsl_vector_alloc(n);
    set_all(init_value);
  }



  vector::vector(const vector& other, const std::vector<size_t>& index)
  {
    if (index.size()){
      v_ = gsl_vector_alloc(index.size());
      for (size_t i=0; i<index.size(); i++) 
        gsl_vector_set(v_, i, other(index[i]));
    }
    else
      v_ = other.gsl_vector_copy();
  }



  vector::vector(gsl_vector* vec)
    : v_(vec)
  {
  }



  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()) {
          // Jari. change to throw exception, remove unnecessary includes
          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);
      }

      // Ignoring empty lines
      if(!v.size()){
        nof_rows--;
        continue;
      }

      if(nof_columns==0)
        nof_columns=v.size();
      else if(v.size()!=nof_columns) {
        // Jari. change to throw exception, remove unnecessary includes
        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 );
      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 );
      for(u_int i=0;i<nof_columns;i++) 
        gsl_vector_set( v_, i, data_matrix[0][i] );
    }
    else {
      // Jari. change to throw exception, remove unnecessary includes
      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(void)
  {
    if (v_) {
      gsl_vector_free(v_);
      v_=NULL;
    }
  }



  gsl_vector* vector::gsl_vector_copy(void) const
  {
    gsl_vector* vec = gsl_vector_alloc(size());
    gsl_vector_memcpy(vec,v_);
    return vec;
  }





  std::pair<size_t,size_t> vector::minmax_index(void) const
  {
    size_t min_index=0;
    size_t max_index=0;
    void gsl_vector_minmax_index (const gsl_vector * v_, size_t * 
                                  min_index, size_t * max_index);
    return std::pair<size_t,size_t>(min_index, max_index);
  }




  std::pair<size_t,size_t> vector::minmax_index(const std::vector<size_t>& subset ) const
  {
    size_t min_index = subset[0];
    size_t max_index = subset[0];
    for (unsigned int i=1; i<subset.size(); i++) {
      if (gsl_vector_get( v_, subset[i]) < gsl_vector_get( v_, min_index))
        min_index = subset[i];
      else if (gsl_vector_get( v_, subset[i]) > gsl_vector_get( v_, max_index))
        max_index = subset[i];
    }
    return std::pair<size_t,size_t>(min_index, max_index);
  }



  vector vector::mul_elements( const vector& other ) const
  {
    vector res( *this );
    gsl_vector_mul( res.v_, other.v_ );
    return res;
  }

  void vector::shuffle(void) const
  {
    vector vec(*this);
    theplu::cpptools::random_singleton* 
      rnd=theplu::cpptools::random_singleton::get_instance();
    for (size_t i=0; i<vec.size(); i++){
      size_t index = rnd->get_uniform_int(vec.size()-i);
      gsl_vector_set(v_,i,vec(index));
      vec(index)=vec(vec.size()-i-1);
    }
  }

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



  vector vector::operator-(void) const
  {
    vector res( *this );
    gsl_vector_scale (res.v_,-1.);
    return res;
  }



  double vector::operator*( const vector &other ) const
  {
    // Jari, check for gsl_support
    double res = 0.0;;
    for ( size_t i = 0; i < size(); ++i ) 
      res += other[i] * (*this)[i];
    return res;
  }



  vector vector::operator+(const vector &other) const
  {
    vector res(*this);
    gsl_vector_add(res.v_,other.v_);
    return res;
  }



  vector vector::operator-( const vector &other ) const
  {
    vector res( *this );
    gsl_vector_sub(res.v_,other.v_);
    return res;
  }



  vector& vector::operator=( const vector& other )
  {
    if( this != &other ) {
      if ( v_ )
        gsl_vector_free( v_ );
      v_ = other.gsl_vector_copy();
    }
    return *this;
  } 



  std::ostream& theplu::gslapi::operator<<(std::ostream& s, const vector& a)
  {
    s.setf(std::ios::fixed);

    for (size_t j = 0; j < a.size(); ++j) {
      s << a[j];
      if ( (j+1)<a.size() )
        s << " ";
    }

    return s;
  }


}} // of namespace gslapi and namespace thep