source: trunk/yat/utility/VectorBase.cc @ 1487

// $Id: VectorBase.cc 1487 2008-09-10 08:41:36Z jari $

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

  This file is part of the yat library, http://dev.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 3 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 yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "VectorBase.h"
#include "Vector.h"
#include "utility.h"
#include "yat/random/random.h"

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

namespace theplu {
namespace yat {
namespace utility {


  VectorBase::VectorBase(const gsl_vector* v)
    : const_vec_(v)
  {
  }


  VectorBase::~VectorBase(void)
  {
  }


  VectorBase::const_iterator VectorBase::begin(void) const
  {
    if (const_vec_)
      return const_iterator(&(this->operator()(0)), const_vec_->stride);
    return const_iterator(NULL, 1);
  }


  VectorBase::const_iterator VectorBase::end(void) const
  {
    if (const_vec_)
      return const_iterator(&(this->operator()(0))+const_vec_->stride*size(), 
                            const_vec_->stride);
    return const_iterator(NULL, 1);
  }


  bool VectorBase::equal(const VectorBase& 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)
      if (fabs( (*this)(i)-other(i) ) > d ||
          std::isnan((*this)(i)) || std::isnan(other(i)) )
        return false;
    return true;
  }


  const gsl_vector* VectorBase::gsl_vector_p(void) const
  {
    return const_vec_;
  }


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


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


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


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


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


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


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


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


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


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


  bool nan(const VectorBase& templat, Vector& flag)
  {
    size_t vsize(templat.size());
    flag = Vector(vsize, 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 sort_index(std::vector<size_t>& sort_index, const VectorBase& 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 VectorBase&)",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 VectorBase& 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 VectorBase& 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 VectorBase& v)
  {
    double sum = 0;
    size_t vsize=v.size();
    for (size_t i=0; i<vsize; ++i)
      sum += v(i);
    return sum;
  }


  std::ostream& operator<<(std::ostream& s, const VectorBase& 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