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

#ifndef _theplu_gslapi_vector_
#define _theplu_gslapi_vector_

// $Id: vector.h 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, 2008 Jari Häkkinen

  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/utility/Exception.h>

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

#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort_vector.h>

namespace theplu {
namespace gslapi {

  class matrix;

  ///
  /// This is the C++ tools interface to GSL vector. 'double' is the
  /// only type supported, maybe we should add a 'complex' type as
  /// well in the future.
  ///
  /// \par[File streams] Reading and writing vectors to file streams
  /// are of course supported. These are implemented without using GSL
  /// functionality, and thus binary read and write to streams are not
  /// supported.
  ///
  /// \par[Vector views] GSL vector views are supported and these are
  /// disguised as ordinary gslapi::vectors. A support function is
  /// added, gslapi::vector::isview(), that can be used to check if a
  /// vector object is a view. Note that view vectors do not own the
  /// undelying data, and a view is not valid if the vector owning the
  /// data is deallocated.
  ///
  /// @note Missing support to underlying GSL vector features can in
  /// principle be included to this class if requested by the user
  /// community.
  ///
  class vector
  {
  public:

    ///
    /// The default constructor.
    ///
    inline vector(void) : v_(NULL), view_(NULL), const_view_(NULL) {}

    ///
    /// Constructor. Allocates memory space for \a n elements, and
    /// sets all elements to \a init_value.
    ///
    inline vector(const size_t n,const double init_value=0)
      : view_(NULL), const_view_(NULL)
      { v_ = gsl_vector_alloc(n); set_all(init_value); }

    ///
    /// The copy constructor.
    ///
    /// @note If the object to be copied is a vector view, the values
    /// of the view will be copied, i.e. the view is not copied.
    ///
    inline vector(const vector& other) : view_(NULL), const_view_(NULL)
      { v_ = other.create_gsl_vector_copy(); }

    ///
    /// Vector view constructor.
    ///
    /// Create a view of vector \a v, with starting index \a offset,
    /// size \a n, and an optional \a stride.
    ///
    /// A vector view can be used as any vector with the difference
    /// that changes made to the view will also change the object that
    /// is viewed. Also, using the copy constructor will create a new
    /// vector object that is a copy of whatever is viewed. If a copy
    /// of the view is needed then you should use this constructor to
    /// obtain a copy.
    ///
    /// @note If the object viewed by the view goes out of scope or is
    /// deleted, the view becomes invalid and the result of further
    /// use is undefined.
    ///
    vector(vector& v, size_t offset, size_t n, size_t stride=1);

    ///
    /// Vector const view constructor.
    ///
    /// Create a view of vector \a v, with starting index \a offset,
    /// size \a n, and an optional \a stride.
    ///
    /// A vector view can be used as any const vector. Using the copy
    /// constructor will create a new vector object that is a copy of
    /// whatever is viewed. If a copy of the view is needed then you
    /// should use this constructor to obtain a copy.
    ///
    /// @note If the object viewed by the view goes out of scope or is
    /// deleted, the view becomes invalid and the result of further
    /// use is undefined.
    ///
    vector(const vector& v, size_t offset, size_t n, size_t stride=1);

    ///
    /// Matrix row/column view constructor.
    ///
    /// Create a row/column vector view of matrix \a m, pointing at
    /// row/column \a i. The parameter \a row is used to set whether
    /// the view should be a row or column view. If \a row is set to
    /// true, the view will be a row view (default behaviour), and,
    /// naturally, a column view otherwise.
    ///
    /// A vector view can be used as any vector with the difference
    /// that changes made to the view will also change the object that
    /// is viewed. Also, using the copy constructor will create a new
    /// vector object that is a copy of whatever is viewed. If a copy
    /// of the view is needed then you should use the vector view
    /// constructor to obtain a copy.
    ///
    /// @note If the object viewed by the view goes out of scope or is
    /// deleted, the view becomes invalid and the result of further
    /// use is undefined.
    ///
    vector(matrix& m, size_t i, bool row=true);

    ///
    /// Matrix row/column const view constructor.
    ///
    /// Create a row/column vector view of matrix \a m, pointing at
    /// row/column \a i. The parameter \a row is used to set whether
    /// the view should be a row or column view. If \a row is set to
    /// true, the view will be a row view (default behaviour), and,
    /// naturally, a column view otherwise.
    ///
    /// A const vector view can be used as any const vector. Using the
    /// copy constructor will create a new vector object that is a
    /// copy of whatever is viewed. If a copy of the view is needed
    /// then you should use the vector view constructor to obtain a
    /// copy.
    ///
    /// @note If the object viewed by the view goes out of scope or is
    /// deleted, the view becomes invalid and the result of further
    /// use is undefined.
    ///
    vector(const matrix& m, size_t i, bool row=true);

    ///
    /// The istream constructor.
    ///
    /// Either elements should be separated 
    /// with white space characters (default), or elements should be separated
    /// by the delimiter \a sep. When delimiter \a sep is used empty elements
    /// are stored as NaN's (except that empty lines are ignored). The
    /// end of input to the vector is at end of file marker.
    ///
    explicit vector(std::istream &, char sep='\0') throw (utility::IO_error,std::exception);

    ///
    /// The string constructor.
    ///
    /// Either elements should be separated 
    /// with white space characters (default), or elements should be separated
    /// by the delimiter \a sep. When delimiter \a sep is used empty elements
    /// are stored as NaN's (except that empty lines are ignored). The
    /// end of input to the vector is at end of file marker.
    ///
    explicit vector(std::string &, char sep='\0') throw (utility::IO_error,std::exception);


    ///
    /// The destructor.
    ///
    ~vector(void);

    ///
    /// Vector addition, \f$this_i = this_i + other_i \; \forall i\f$.
    ///
    /// @return GSL_SUCCESS on normal exit.
    ///
    // Jari, group as vector_operators
    inline int add(const vector& other) { return gsl_vector_add(v_,other.v_); }

    ///
    /// Add a constant to a vector, \f$this_i = this_i + term \;
    /// \forall i\f$.
    ///
    /// @return GSL_SUCCESS on normal exit.
    ///
    // Jari, group as vector_operators
    inline int add(double term) { return gsl_vector_add_constant(v_,term); }

    ///
    /// This function performs element-wise division, \f$this_i =
    /// this_i/other_i \; \forall i\f$.
    ///
    /// @return GSL_SUCCESS on normal exit.
    ///
    // Jari, doxygen group as Vector operators
    inline int div(const vector& other) { return gsl_vector_div(v_,other.v_); }

    ///
    /// @return A const pointer to the internal GSL vector,
    ///
    inline const gsl_vector* gsl_vector_p(void) const { return v_; }

    ///
    /// @return A pointer to the internal GSL vector,
    ///
    inline gsl_vector* gsl_vector_p(void) { return v_; }

    ///
    /// @return True if all elements in the vector is zero, false
    /// othwerwise;
    ///
    inline bool isnull(void) const { return gsl_vector_isnull(v_); }

    ///
    /// Check if the vector object is a view (sub-vector) to another
    /// vector.
    ///
    /// @return True if the object is a view, false othwerwise.
    ///
    inline bool isview(void) const { return view_ || const_view_; }

    ///
    /// @return The maximum value of the vector.
    ///
    // Jari, doxygen group as Finding maximum and minimum elements
    inline double max(void) const { return gsl_vector_max(v_); }

    ///
    /// @return The element index to the maximum value of the
    /// vector. The lowest index has precedence.
    ///
    // Jari, doxygen group as Finding maximum and minimum elements
    inline size_t max_index(void) const { return gsl_vector_max_index(v_); }

    ///
    /// @return The minimum value of the vector.
    ///
    // Jari, doxygen group as Finding maximum and minimum elements
    inline double min(void) const { return gsl_vector_min(v_); }

    ///
    /// @return The element index to the minimum value of the
    /// vector. The lowest index has precedence.
    ///
    // Jari, doxygen group as Finding maximum and minimum elements
    inline size_t min_index(void) const { return gsl_vector_min_index(v_); }

    ///
    /// @return The minimum and maximum values of the vector, as the
    /// \a first and \a second member of the returned \a pair,
    /// respectively.
    ///
    // Jari, doxygen group as Finding maximum and minimum elements
    std::pair<double,double> minmax(void) const;

    ///
    /// @return The indecies to the minimum and maximum values of the
    /// vector, as the \a first and \a second member of the returned
    /// \a pair, respectively. The lowest index has precedence.
    ///
    // Jari, doxygen group as Finding maximum and minimum elements
    std::pair<size_t,size_t> minmax_index(void) const;

    ///
    /// This function performs element-wise multiplication, \f$this_i =
    /// this_i * other_i \; \forall i\f$.
    ///
    /// @return GSL_SUCCESS on normal exit.
    ///
    // Jari, doxygen group as Vector operators
    inline int mul(const vector& other) { return gsl_vector_mul(v_,other.v_); }

    ///
    /// Reverse the order of elements in the vector.
    ///
    /// @return GSL_SUCCESS on normal exit.
    ///
    // Jari, doxygen group as Exchanging elements
    inline int reverse(void) { return gsl_vector_reverse(v_);}

    ///
    /// Rescale vector, \f$this_i = this_i * factor \; \forall i\f$.
    ///
    /// @return GSL_SUCCESS on normal exit.
    ///
    // Jari, doxygen group as Vector operators
    inline int scale(double factor) { return gsl_vector_scale(v_,factor); }

    ///
    /// Set element values to values in \a vec. This function is
    /// needed for assignment of viewed elements.
    ///
    /// @return Whatever GSL returns.
    ///
    /// @note No check on size is done.
    ///
    /// @see const vector& operator=(const vector&)
    ///
    inline int set(const vector& vec) { return gsl_vector_memcpy(v_,vec.v_); }

    ///
    /// Set all elements to \a value.
    ///
    // Jari, doxygen group as Initializing vector elements
    inline void set_all(const double& value) { gsl_vector_set_all(v_,value); }

    ///
    /// Makes a basis vector by setting all elements to
    /// zero except the \a i-th element which is set to
    /// one.
    ///
    // Jari, doxygen group as Initializing vector elements
    inline void set_basis(const size_t i) { gsl_vector_set_basis(v_,i); }
    
    ///
    /// Set all elements to zero.
    ///
    // Jari, doxygen group as Initializing vector elements
    inline void set_zero(void) { gsl_vector_set_zero(v_); }

    ///
    /// @return the number of elements in the vector.
    ///
    inline size_t size(void) const { return v_->size; }

    ///
    /// Sort the elements in the vector.
    ///
    /// Bug in gsl: if vector contains NaN an infinite loop is entered.
    ///
    // Markus to Jari, doxygen group as Exchanging elements ????
    inline void sort(void) { gsl_sort_vector(v_);}


    ///
    /// Vector subtraction, \f$this_i = this_i - other_i \; \forall i\f$.
    ///
    /// @return GSL_SUCCESS on normal exit.
    ///
    // Jari, doxygen group as Vector operators
    inline int sub(const vector& other) { return gsl_vector_sub(v_,other.v_); }

    ///
    /// Calculate the sum of all vector elements.
    ///
    /// @return The sum.
    ///
    double sum(void) const;

    ///
    /// Swap vector elements by copying. The two vectors must have the
    /// same length.
    ///
    /// @return GSL_SUCCESS on normal exit.
    ///
    inline int swap(vector& other) { return gsl_vector_swap(v_,other.v_); }

    ///
    /// Exchange elements \a i and \a j.
    ///
    /// @return GSL_SUCCESS on normal exit.
    ///
    // Jari, doxygen group as Exchanging elements
    inline int
    swap_elements(size_t i,size_t j) { return gsl_vector_swap_elements(v_,i,j);}

    ///
    /// Element access operator.
    ///
    /// @return Reference to element \a i.
    ///
    // Jari, doxygen group as Accessing vector elements
    inline double& operator()(size_t i) { return *gsl_vector_ptr(v_,i); }

    ///
    /// Const element access operator.
    ///
    /// @return The value of element \a i.
    ///
    // Jari, doxygen group as Accessing vector elements
    inline const double&
    operator()(size_t i) const { return *gsl_vector_const_ptr(v_,i); }

    ///
    /// Element access operator.
    ///
    /// @return Reference to element \a i.
    ///
    // Jari, doxygen group as Accessing vector elements
    inline double& operator[](size_t i) { return *gsl_vector_ptr(v_,i); }

    ///
    /// Const element access operator.
    ///
    /// @return The value of element \a i.
    ///
    // Jari, doxygen group as Accessing vector elements
    inline const double&
    operator[](size_t i) const { return *gsl_vector_const_ptr(v_,i); }

    ///
    /// Comparison operator. Takes linear time.
    ///
    /// @return True if the sequence of the elements in the vectors
    /// are element/wise equal.
    ///
    bool operator==(const vector&) const;

    ///
    /// The assignment operator. There is no requirements on
    /// dimensions, i.e. the vector is remapped in memory if
    /// necessary. This implies that in general views cannot be
    /// assigned using this operator. Views will be mutated into
    /// normal vectors. The only exception to this behaviour on views
    /// is when self-assignemnt is done, since self-assignment is
    /// ignored.
    ///
    /// @return A const reference to the resulting vector.
    ///
    /// @see int set(const vector&)
    ///
    const vector& operator=(const vector&);

    ///
    /// @return The dot product.
    ///
    double operator*( const vector &other ) const;

    ///
    /// Addition and assign operator.
    ///
    /// @return A const reference to the resulting vector.
    ///
    inline const vector&
    operator+=(const vector& other) { gsl_vector_add(v_,other.v_); return *this;}

    ///
    /// Subtract and assign operator.
    ///
    /// @return A const reference to the resulting vector.
    ///
    inline const vector&
    operator-=(const vector& other) { gsl_vector_sub(v_,other.v_); return *this;}

    ///
    /// Multiply with scalar and assign operator.
    ///
    /// @return A const reference to the resulting vector.
    ///
    inline const vector&
    operator*=(const double d) { gsl_vector_scale(v_,d); return *this; }


  private:

    ///
    /// Create a new copy of the internal GSL vector.
    ///
    /// Necessary memory for the new GSL vector is allocated and the
    /// caller is responsible for freeing the allocated memory.
    ///
    /// @return A pointer to a copy of the internal GSL vector.
    ///
    gsl_vector* create_gsl_vector_copy(void) const;

    gsl_vector* v_;
    gsl_vector_view* view_;
    gsl_vector_const_view* const_view_;
  };

  ///
  /// The output operator for the vector class.
  ///
  std::ostream& operator<<(std::ostream&, const vector& );


}} // of namespace gslapi and namespace theplu

#endif