source: trunk/se/lu/thep/wenni/lib/c++_tools/gslapi/vector.cc @ 597

// $Id: vector.cc 597 2008-02-27 23:36:46Z jari $

/*
  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005 Jari Häkkinen, Peter Johansson, Markus Ringnér
  Copyright (C) 2006 Jari Häkkinen, Peter Johansson

  This file is part of the thep c++ tools library,
                                http://lev.thep.lu.se/trac/c++_tools

  The c++ tools library is free software; you can redistribute it
  and/or modify it under the terms of the GNU General Public License
  as published by the Free Software Foundation; either version 2 of
  the License, or (at your option) any later version.

  The c++ tools library is distributed in the hope that it will be
  useful, but WITHOUT ANY WARRANTY; without even the implied warranty
  of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#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 <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(std::string& line, char sep)
    throw (utility::IO_error,std::exception)
    : view_(NULL), const_view_(NULL)
  {
    // Empty line
    if (!line.size())
      return;

    std::vector<double> v;
    v.reserve(line.length()/2);
    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());
      }
    }
    if (sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
      v.push_back(std::numeric_limits<double>::quiet_NaN());

    // convert the data to a gsl vector
    v_ = gsl_vector_alloc(v.size());
    size_t n=0;
    for (size_t i=0; i<v.size(); i++)
      gsl_vector_set(v_, n++, v[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