source: trunk/yat/utility/vector.cc @ 880

// $Id: vector.cc 880 2007-09-21 23:43:22Z peter $

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

  This file is part of the yat library, http://trac.thep.lu.se/trac/yat

  The yat 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 yat 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 "vector.h"
#include "matrix.h"
#include "utility.h"
#include "yat/random/random.h"

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

namespace theplu {
namespace yat {
namespace utility {


  vector::vector(void)
    : v_(NULL), view_(NULL), view_const_(NULL), proxy_v_(NULL)
  {
  }


  vector::vector(size_t n, double init_value)
    : v_(gsl_vector_alloc(n)), view_(NULL), view_const_(NULL),
      proxy_v_(v_)
  {
    if (!v_)
      throw utility::GSL_error("vector::vector failed to allocate memory");

    all(init_value);
  }


  vector::vector(const vector& other)
    : v_(other.create_gsl_vector_copy()), view_(NULL),
      view_const_(NULL), proxy_v_(v_)
  {
  }


  vector::vector(vector& v, size_t offset, size_t n, size_t stride)
    : view_const_(NULL)
  {
    view_ = new gsl_vector_view(gsl_vector_subvector_with_stride(v.v_,offset,
                                                                 stride,n));
    if (!view_)
      throw utility::GSL_error("vector::vector failed to setup view");
    proxy_v_ = v_ = &(view_->vector);
  }


  vector::vector(const vector& v, size_t offset, size_t n, size_t stride)
    : v_(NULL), view_(NULL)
  {
    view_const_ = new gsl_vector_const_view(
                   gsl_vector_const_subvector_with_stride(v.v_,offset,stride,n));
    if (!view_const_)
      throw utility::GSL_error("vector::vector failed to setup view");
    proxy_v_ = &(view_const_->vector);
  }


  vector::vector(matrix& m, size_t i, bool row)
    : view_const_(NULL)
  {
    view_=new gsl_vector_view(row ?
                              gsl_matrix_row   (m.gsl_matrix_p(),i) :
                              gsl_matrix_column(m.gsl_matrix_p(),i)  );
    if (!view_)
      throw utility::GSL_error("vector::vector failed to setup view");
    proxy_v_ = v_ = &(view_->vector);
  }


  vector::vector(const matrix& m, size_t i, bool row)
    : v_(NULL), view_(NULL)
  {
    view_const_ = new gsl_vector_const_view(row ?
                                   gsl_matrix_const_row   (m.gsl_matrix_p(),i) :
                                   gsl_matrix_const_column(m.gsl_matrix_p(),i) );
    if (!view_const_)
      throw utility::GSL_error("vector::vector failed to setup view");
    proxy_v_ = &(view_const_->vector);
  }


  vector::vector(std::istream& is, char sep)
    throw (utility::IO_error, std::exception)
    : view_(NULL), view_const_(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 {
          std::stringstream ss("Warning: '");
          ss << element << "' is not a double.";
          throw IO_error(ss.str());
        }
      } 
      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
    proxy_v_ = v_ = gsl_vector_alloc(nof_rows*nof_columns);
    if (!v_)
      throw utility::GSL_error("vector::vector failed to allocate memory");
    size_t n=0;
    // if gsl error handler disabled, out of bounds index will not
    // abort the program.
    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(void)
  {
    delete_allocated_memory();
  }


  vector::iterator vector::begin(void)
  {
    return vector::iterator(*this, 0);
  }


  vector::const_iterator vector::begin(void) const
  {
    return vector::const_iterator(*this, 0);
  }


  const vector& vector::clone(const vector& other)
  {
    if (this!=&other) {

      delete_allocated_memory();

      if (other.view_) {
        view_ = new gsl_vector_view(*other.view_);
        proxy_v_ = v_ = &(view_->vector);
      }
      else if (other.view_const_) {
        view_const_ = new gsl_vector_const_view(*other.view_const_);
        proxy_v_ = &(view_const_->vector);
        v_=NULL;
      }
      else if (other.v_)
        proxy_v_ = v_ = other.create_gsl_vector_copy();

    }
    return *this;
  } 


  gsl_vector* vector::create_gsl_vector_copy(void) const
  {
    gsl_vector* vec = gsl_vector_alloc(size());
    if (!vec)
      throw utility::GSL_error("vector::create_gsl_vector_copy failed to allocate memory");
    if (gsl_vector_memcpy(vec, proxy_v_))
      throw utility::GSL_error("vector::create_gsl_matrix_copy dimension mis-match");
    return vec;
  }


  void vector::delete_allocated_memory(void)
  {
    if (view_)
      delete view_;
    else if (view_const_)
      delete view_const_;
    else if (v_)
      gsl_vector_free(v_);
    proxy_v_=v_=NULL;
  }


  void vector::div(const vector& other)
  {
    assert(v_);
    int status=gsl_vector_div(v_,other.gsl_vector_p());
    if (status)
      throw utility::GSL_error(std::string("vector::div",status));
  }


  vector::iterator vector::end(void)
  {
    return vector::iterator(*this, size());
  }


  vector::const_iterator vector::end(void) const
  {
    return vector::const_iterator(*this, size());
  }


  bool vector::equal(const vector& other, const double d) const
  {
    if (this==&other)
      return true;
    if (size()!=other.size())
      return false;
    // if gsl error handler disabled, out of bounds index will not
    // abort the program.
    for (size_t i=0; i<size(); ++i)
      // The two last condition checks are needed for NaN detection
      if (fabs( (*this)(i)-other(i) ) > d ||
          (*this)(i)!=(*this)(i) || other(i)!=other(i))
        return false;
    return true;
  }


  const gsl_vector* vector::gsl_vector_p(void) const
  {
    return proxy_v_;
  }


  gsl_vector* vector::gsl_vector_p(void)
  {
    return v_;
  }


  bool vector::isview(void) const
  {
    return view_ || view_const_;
  }


  void vector::mul(const vector& other)
  {
    assert(v_);
    int status=gsl_vector_mul(v_,other.gsl_vector_p());
    if (status)
      throw utility::GSL_error(std::string("vector::mul",status));
  }


  void vector::resize(size_t n, double init_value)
  {
    delete_allocated_memory();
    proxy_v_ = v_ = gsl_vector_alloc(n);
    if (!v_)
      throw utility::GSL_error("vector::vector failed to allocate memory");
    all(init_value);
  } 


  void vector::reverse(void)
  {
    assert(v_);
    gsl_vector_reverse(v_);
  }


  void vector::all(const double& value)
  {
    assert(v_);
    gsl_vector_set_all(v_,value);
  }


  size_t vector::size(void) const
  {
    if (!proxy_v_)
      return 0;
    return proxy_v_->size;
  }


  void vector::swap(size_t i, size_t j)
  {
    assert(v_);
    int status=gsl_vector_swap_elements(v_, i, j);
    if (status)
      throw utility::GSL_error(std::string("vector::swap_elements",status));
  }


  double& vector::operator()(size_t i)
  {
    assert(v_);
    double* d=gsl_vector_ptr(v_, i);
    if (!d)
      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
    return *d;
  }


  const double& vector::operator()(size_t i) const
  {
    const double* d=gsl_vector_const_ptr(proxy_v_, i);
    if (!d)
      throw utility::GSL_error("vector::operator()",GSL_EINVAL);
    return *d;
  }


  double& vector::operator[](size_t i)
  {
    return this->operator()(i);
  }


  const double& vector::operator[](size_t i) const
  {
    return this->operator()(i);
  }


  bool vector::operator==(const vector& other) const
  {
    return equal(other);
  }


  bool vector::operator!=(const vector& other) const
  {
    return !equal(other);
  }


  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 )
  {
    assert(v_);
    if (gsl_vector_memcpy(v_,other.gsl_vector_p()))
      throw utility::GSL_error("vector::set dimension mis-match");
    return *this;
  } 


  const vector& vector::operator+=(const vector& other)
  {
    assert(v_);
    int status=gsl_vector_add(v_, other.gsl_vector_p());
    if (status)
      throw utility::GSL_error(std::string("vector::add", status));
    return *this;
  }


  const vector& vector::operator+=(double d)
  {
    assert(v_);
    gsl_vector_add_constant(v_, d);
    return *this;
  }


  const vector& vector::operator-=(const vector& other)
  {
    assert(v_);
    int status=gsl_vector_sub(v_, other.gsl_vector_p());
    if (status)
      throw utility::GSL_error(std::string("vector::sub", status));
    return *this;
  }


  const vector& vector::operator-=(const double d)
  {
    assert(v_);
    gsl_vector_add_constant(v_, -d);
    return *this;
  }


  const vector& vector::operator*=(const double d)
  {
    assert(v_);
    gsl_vector_scale(v_, d);
    return *this;
  }


  bool isnull(const vector& v)
  {
    return gsl_vector_isnull(v.gsl_vector_p());
  }


  double max(const vector& v)
  {
    return gsl_vector_max(v.gsl_vector_p());
  }


  size_t max_index(const vector& v)
  {
    return gsl_vector_max_index(v.gsl_vector_p());
  }


  double min(const vector& v)
  {
    return gsl_vector_min(v.gsl_vector_p());
  }


  size_t min_index(const vector& v)
  {
    return gsl_vector_min_index(v.gsl_vector_p());
  }


  bool nan(const vector& templat, vector& flag)
  {
    size_t vsize=templat.size();
    if (vsize!=flag.size())
      flag.clone(vector(vsize,1.0));
    else
      flag.all(1.0);
    bool nan=false;
    for (size_t i=0; i<vsize; i++)
      if (std::isnan(templat(i))) {
        flag(i)=0;
        nan=true;
      }
    return nan;
  }


  void set_basis(vector& v, size_t i)
  {
    assert(v.gsl_vector_p());
    gsl_vector_set_basis(v.gsl_vector_p(),i);
  }


  void shuffle(vector& invec)
  {
    vector vec(invec);
    random::DiscreteUniform rnd;
    for (size_t i=0; i<vec.size(); i++){
      size_t index = rnd(vec.size()-i);
      invec[i]=vec(index);
      vec(index)=vec(vec.size()-i-1);
    }
  }


  void sort(vector& v)
  {
    assert(v.gsl_vector_p());
    gsl_sort_vector(v.gsl_vector_p());
  }


  void sort_index(std::vector<size_t>& sort_index, const vector& invec)
  {
    assert(invec.gsl_vector_p());
    gsl_permutation* p = gsl_permutation_alloc(invec.size());
    int status=gsl_sort_vector_index(p,invec.gsl_vector_p());
    if (status) {
      gsl_permutation_free(p);
      throw utility::GSL_error(std::string("sort_index(vector&,const vector&)",status));      
    }
    sort_index=std::vector<size_t>(p->data,p->data+p->size);
    gsl_permutation_free(p);
  }


  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, 
                            const vector& invec)
  {
    assert(invec.gsl_vector_p());
    assert(k<=invec.size());
    sort_index.resize(k);
    gsl_sort_vector_smallest_index(&sort_index[0],k,invec.gsl_vector_p());
  }
  
  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, 
                            const vector& invec)
  {
    assert(invec.gsl_vector_p());
    assert(k<=invec.size());
    sort_index.resize(k);
    gsl_sort_vector_largest_index(&sort_index[0],k,invec.gsl_vector_p());
  }


  double sum(const vector& v)
  {
    double sum = 0;
    size_t vsize=v.size();
    for (size_t i=0; i<vsize; ++i)
      sum += v[i];
    return sum;
  }


  void swap(vector& v, vector& w)
  {
    assert(v.gsl_vector_p()); assert(w.gsl_vector_p());
    int status=gsl_vector_swap(v.gsl_vector_p(),w.gsl_vector_p());
    if (status)
      throw utility::GSL_error(std::string("swap(vector&,vector&)",status));
  }


  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 << s.fill();
    }
    return s;
  }

}}} // of namespace utility, yat, and thep