source: trunk/yat/utility/vector.h @ 754

#ifndef _theplu_yat_utility_vector_
#define _theplu_yat_utility_vector_

// $Id: vector.h 754 2007-02-17 22:33:44Z 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, 2007 Jari Häkkinen

  This file is part of the yat library, http://lev.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 "Exception.h"

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

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

namespace theplu {
namespace yat {
namespace utility {

  class matrix;

  /**
     This is the yat 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 utility::vectors. A support function is added,
     utility::vector::isview(), that can be used to check if a vector
     object is a view. Note that view vectors do not own the
     underlying data, and a view is not valid if the vector owing the
     data is deallocated.

     \par
     Currently there is no restriction on how a vector is used in
     cases when the vector is a const view into another vector or
     matrix. To avoid unexpected runtime errors, the programmer must
     declare const view vectors as const in order to get compile time
     diagnostics about improper use of a const view vector object. An
     example of improper use of a const view vector is assignment of
     values to an element in the object. If the const view object was
     declared const the assigment will be caught by the compiler
     whereas it will cause a (un-catchable) runtime error. Example:
     @code
  const vector vm(13,1.0);
  vector v1(vm,2,4);       // bad code! not const!
  v1(0)=-123;              // accepted by compiler, runtime abort will occur
  const vector v2(vm,2,4); // proper code
  v2(0)=-123;              // not acceptable for the compiler
     @endcode

     @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:

    /**
       \brief The default constructor.
    */
    vector(void);

    /**
       \brief Allocates memory space for \a n elements, and sets all
       elements to \a init_value.
       
       \throw GSL_error if memory allocation fails.
    */
    vector(size_t n, double init_value=0);

    /**
       \brief 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.

       \throw A GSL_error is indirectly thrown if memory allocation
       fails.
    */
    vector(const vector& other);

    /**
       \brief 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.

       \throw GSL_error if a view cannot be set up.
    */
    vector(vector& v, size_t offset, size_t n, size_t stride=1);

    /**
       \brief 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.

       \throw GSL_error if a view cannot be set up.
    */
    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);

    /**
       \brief 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.

       \throw GSL_error if memory allocation fails.
    */
    explicit vector(std::istream &, char sep='\0')
      throw(utility::IO_error, std::exception);

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

    /**
       \brief Vector addition, \f$ this_i = this_i + other_i \;
       \forall i \f$.

       \throw GSL_error if dimensions mis-match.
    */
    void add(const vector& other);

    /**
       \brief Add a constant to a vector, \f$ this_i = this_i + term \;
       \forall i \f$.
    */
    void add(double 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
    int div(const vector& other);

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

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

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

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

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

    ///
    /// @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
    size_t max_index(void) const;

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

    ///
    /// @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
    size_t min_index(void) const;

    ///
    /// @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
    int mul(const vector& other);

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

    /**
       \brief Rescale vector, \f$ this_i = this_i * factor \; \forall i \f$.
    */
    void scale(double factor);

    /**
       \brief Set element values to values in \a vec.

       This function is needed for assignment of viewed elements.

       \see const vector& operator=(const vector&)

       \throw GSL_error if dimensions mis-match.
    */
    void set(const vector& vec);

    ///
    /// Set all elements to \a value.
    ///
    // Jari, doxygen group as Initializing vector elements
    void set_all(const double& 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
    void set_basis(const size_t i);
    
    ///
    /// Set all elements to zero.
    ///
    // Jari, doxygen group as Initializing vector elements
    void set_zero(void);

    ///
    /// @return the number of elements in the vector.
    ///
    size_t size(void) const;

    ///
    /// 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 ????
    void sort(void);

    /**
       \brief Vector subtraction, \f$ this_i = this_i - other_i \;
       \forall i \f$.

       \throw GSL_error if dimensions mis-match.
    */
    void sub(const vector& other);

    ///
    /// 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.
    ///
    int swap(vector& other);

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

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

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

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

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

    ///
    /// 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;

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

    /**
       \brief 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 void set(const vector&).

       \throw A GSL_error is indirectly thrown if memory allocation
       fails, or if dimensions mis-match.
    */
    const vector& operator=(const vector&);

    /**
       \brief Addition and assign operator.

       \return A const reference to the resulting vector.

       \throw GSL_error if dimensions mis-match.
    */
    const vector& operator+=(const vector&);

    /**
       \brief Subtract and assign operator.

       \return A const reference to the resulting vector.

       \throw GSL_error if dimensions mis-match.
    */
    const vector& operator-=(const vector&);

    ///
    /// Multiply with scalar and assign operator.
    ///
    /// @return A const reference to the resulting vector.
    ///
    const vector& operator*=(const double);


  private:

    /**
       \brief 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.

       \throw GSL_error if memory cannot be allocated for the new
       copy, or if dimensions mis-match.
    */
    gsl_vector* create_gsl_vector_copy(void) const;

    gsl_vector* v_;
    const gsl_vector* v_const_;
    gsl_vector_view* view_;
    gsl_vector_const_view* view_const_;
    // proxy_v_ is used to access the proper underlying v_ or v_const_
    // in all const member functions.
    const gsl_vector* proxy_v_;
  };

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


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

#endif