source: trunk/lib/gslapi/vector.cc @ 527

// $Id: vector.cc 527 2006-03-01 11:23:53Z peter $

#include <c++_tools/gslapi/vector.h>
#include <c++_tools/gslapi/matrix.h>
#include <c++_tools/utility/stl_utility.h>
#include <c++_tools/utility/utility.h>
#include <c++_tools/classifier/DataLookup1D.h>

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


namespace theplu {
namespace gslapi {



  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
    : const_view_(NULL)
  {
    // Jari, exception handling needed here. Failure in setting up a
    // proper gsl_vector_view is communicated by NULL pointer in the
    // view structure (cf. GSL manual). How about GSL error state?
    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
                                                                 stride,n));
    v_ = &(view_->vector);
  }



  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
    : view_(NULL)
  {
    // Jari, exception handling needed here. Failure in setting up a
    // proper gsl_vector_view is communicated by NULL pointer in the
    // view structure (cf. GSL manual). How about GSL error state?
    const_view_ = new gsl_vector_const_view(
                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
  }



  vector::vector(matrix& m, size_t i, bool row)
    : const_view_(NULL)
  {
    view_=new gsl_vector_view(row ?
                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
    v_ = &(view_->vector);
  }



  vector::vector(const matrix& m, size_t i, bool row)
    : view_(NULL)
  {
    const_view_ = new gsl_vector_const_view(row ?
                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
    v_ = const_cast<gsl_vector*>(&(const_view_->vector));
  }



  vector::vector(std::istream& is, char sep) 
    throw (utility::IO_error,std::exception)
    : view_(NULL), const_view_(NULL)
  {
    // 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;
    std::string line;
    while(getline(is, line, '\n')) {
      // Empty lines
      if (!line.size())
          continue;
      nof_rows++;
      
      std::vector<double> v;
      std::string element;
      std::stringstream ss(line);
      bool ok=true;
      while(ok) {
        if(sep=='\0')
          ok=(ss>>element);
        else 
          ok=getline(ss, element, sep);
        if(!ok)
          break;        

        if(utility::is_double(element)) {
          v.push_back(atof(element.c_str()));
        }
        else if (!element.size() || utility::is_nan(element)) {
          v.push_back(std::numeric_limits<double>::quiet_NaN());
        }
        else {
          // Jari, this should be communicated with as an exception.
          // std::cerr << "Warning: '" << element 
          //           << "' is not an integer." << std::endl;
        }
      } 
      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
          v.push_back(std::numeric_limits<double>::quiet_NaN());
      if (!nof_columns)
        nof_columns=v.size();
      else if (nof_rows && (nof_columns>1)) {
        std::ostringstream s;
        s << "vector::vector(std::istream&) data file error:\n"
          << "    File has inconsistent number of rows (" << nof_rows
          << ") and columns (" << nof_columns
          << ").\n    Expected a row or a column vector.";
        throw utility::IO_error(s.str());
      }
      else if (v.size()!=nof_columns) {
        std::ostringstream s;
        s << "vector::vector(std::istream&) data file error:\n"
          << "    Line " << nof_rows << " has " << v.size()
          << " columns; expected " << nof_columns << " column.";
        throw utility::IO_error(s.str());
      }
      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
    v_ = gsl_vector_alloc(nof_rows*nof_columns);
    size_t n=0;
    for (size_t i=0; i<nof_rows; i++)
      for (size_t j=0; j<nof_columns; j++) 
        gsl_vector_set( v_, n++, data_matrix[i][j] );
  }



  vector::vector(const classifier::DataLookup1D& data)
    : view_(NULL), const_view_(NULL)
  {
    v_ = gsl_vector_alloc(data.size());
    size_t n=0;
    for (size_t i=0; i<data.size(); i++)
      gsl_vector_set( v_, n++, data(i) );
    
  }



  vector::~vector(void)
  {
    if (view_)
      delete view_;
    else if (const_view_)
      delete const_view_;
    else if (v_)
      gsl_vector_free(v_);
    v_=NULL;
  }



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



  std::pair<double,double> vector::minmax(void) const
  {
    double min, max;
    gsl_vector_minmax(v_,&min,&max);
    return std::pair<double,double>(min,max);
  }



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



  double vector::sum(void) const
  {
    double sum = 0;
    for (size_t i=0; i<size(); i++)
      sum += (*this)(i);
    return( sum );
  }  

  bool vector::operator==(const vector& a) const
  {
    if (size()!=a.size())
      return false;
    for (size_t i=0; i<size(); ++i)
      if (gsl_vector_get(v_,i)!=a(i))
        return false;
    return true;
  }



  double vector::operator*( const vector &other ) const
  {
    double res = 0.0;;
    for ( size_t i = 0; i < size(); ++i ) 
      res += other(i) * (*this)(i);
    return res;
  }



  const vector& vector::operator=( const vector& other )
  {
    if( this != &other ) {
      if (view_)
        delete view_;
      else if (const_view_)
        delete const_view_;
      else if ( v_ )
        gsl_vector_free( v_ );
      v_ = other.create_gsl_vector_copy();
    }
    return *this;
  } 



  std::ostream& operator<<(std::ostream& s, const vector& a)
  {
    s.setf(std::ios::dec);
    s.precision(12);
    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