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

#ifndef _theplu_yat_utility_vector_
#define _theplu_yat_utility_vector_

// $Id: vector.h 830 2007-03-24 14:54:42Z 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
  Copyright (C) 2007 Jari Häkkinen, Peter Johansson

  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;

  /**
     @brief This is the yat interface to GSL vector. 

     For time being 'double' is the only type supported.

     \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/matrix
     owing the data is deallocated.

     \par
     Currently there is no restriction on how a vector is used 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. If
     'const' is not used and the const view vector is used erroneously
     (such as on the left hand side in assignments), the compiler will
     not catch the error and runtime error will occur. assert(3) is
     used to catch the runtime error during development. Example on
     bad and proper use of const view vectors:
     @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
                           // or catched by assert depending on compiler flags
  const vector v2(vm,2,4); // proper code
  v2(0)=-123;              // not acceptable for the compiler
     @endcode
  */

  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, IO_error if
       unexpected input is found in the input stream.
    */
    explicit vector(std::istream &, char sep='\0')
      throw(utility::IO_error, std::exception);

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

    /**
       \brief Make a copy of \a other.

       This function will make a deep copy of \a other. Memory is
       resized and view state is changed if needed.
    */
    const vector& clone(const vector& other);

    /**
       \brief This function performs element-wise division, \f$ this_i =
       this_i/other_i \; \forall i \f$.

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

    /**
       \brief Check whether vectors are equal within a user defined
       precision, set by \a precision.

       \return True if each element deviates less or equal than \a
       d. If any vector contain a NaN, false is always returned.

       \see operator== and operator!=
    */
    bool equal(const vector&, const double precision=0) const;

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

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

    /**
       \brief This function performs element-wise multiplication, \f$
       this_i = this_i * other_i \; \forall i \f$.

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

    /**
       @brief Resize vector
       
       All elements are set to @a init_value.

       \note Underlying GSL vector is destroyed and a view into this
       vector becomes invalid.
    */
    void resize(size_t, double init_value=0);

    /**
       \brief Reverse the order of elements in the vector.
    */
    void reverse(void);

    ///
    /// Set all elements to \a value.
    ///
    void all(const double& value);

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

    /**
       \brief Exchange elements \a i and \a j.

       \throw GSL_error if vector lengths differs.
    */
    void swap(size_t i, size_t j);

    /**
       \brief Element access operator.

       \return Reference to element \a i.

       \throw If GSL range checks are enabled in the underlying GSL
       library a GSL_error exception is thrown if either index is out
       of range.
    */
    double& operator()(size_t i);

    /**
       \brief Element access operator.

       \return Const reference to element \a i.

       \throw If GSL range checks are enabled in the underlying GSL
       library a GSL_error exception is thrown if either index is out
       of range.
    */
    const double& operator()(size_t i) const;

    ///
    /// Element access operator.
    ///
    /// @return Reference to element \a i.
    ///
    double& operator[](size_t i);

    ///
    /// Const element access operator.
    ///
    /// @return The value of element \a i.
    ///
    const double& operator[](size_t i) const;

    /**
       \brief Comparison operator. Takes linear time.

       Checks are performed with exact matching, i.e., rounding off
       effects may destroy comparison. Use the equal function for
       comparing elements within a user defined precision.

       \return True if all elements are equal otherwise false.

       \see equal(const vector&, const double precision=0)
    */
    bool operator==(const vector&) const;

    /**
       \brief Comparison operator. Takes linear time.

       Checks are performed with exact matching, i.e., rounding off
       effects may destroy comparison. Use the equal function for
       comparing elements within a user defined precision.

       \return False if all elements are equal otherwise true.

       \see equal(const vector&, const double precision=0)
    */
    bool operator!=(const vector&) const;

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

    /**
       \brief The assignment operator.

       Dimensions of the vectors must match. If the LHS vector is a
       view, the underlying data will be changed.

       \return A const reference to the resulting vector.

       \see void set(const vector&).

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

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

       \return A const reference to the resulting vector.

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

    /**
       \brief Add a constant to a vector, \f$ this_i = this_i + d \;
       \forall i \f$.

       \return A const reference to the resulting vector.
    */
    const vector& operator+=(double d);

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

       \return A const reference to the resulting vector.

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

    /**
       \brief Subtract a constant to a vector, \f$ this_i = this_i - d
       \; \forall i \f$.

       \return A const reference to the resulting vector.
    */
    const vector& operator-=(double d);

    /**
       \brief Multiply with scalar and assign operator, \f$ this_i =
       this_i * d \; \forall i \f$.

       \return A const reference to the resulting vector.
    */
    const vector& operator*=(double d);


  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;

    /**
       \brief Clear all dynamically allocated memory.

       Internal utility function.
    */
    void delete_allocated_memory(void);

    gsl_vector* v_;
    gsl_vector_view* view_;
    const gsl_vector_const_view* view_const_;
    // proxy_v_ is used to access the proper underlying gsl_vector
    // in all const member functions. It is not used by design for
    // non-const vector functions and operators. This is to make sure
    // that runtime errors occur if a const vector is used in an
    // inappropriate manner such as on left hand side in assignment
    // (remember, v_ is null for const vector views).
    const gsl_vector* proxy_v_;
  };

  /**
     \brief Check if all elements of the vector are zero.

     \return True if all elements in the vector is zero, false
     othwerwise.
  */
  bool isnull(const vector&);

  /**
     \brief Get the maximum value of the vector.

     \return The maximum value of the vector.
  */
  double max(const vector&);

  /**
     \brief Locate the maximum value in the vector.

     \return The index to the maximum value of the vector.

     \note Lower index has precedence.
  */
  size_t max_index(const vector&);

  /**
     \brief Get the minimum value of the vector.

     \return The minimum value of the vector.
  */
  double min(const vector&);

  /**
     \brief Locate the minimum value in the vector.

     \return The index to the minimum value of the vector.

     \note Lower index has precedence.
  */
  size_t min_index(const vector&);

  /**
     \brief Create a vector \a flag indicating NaN's in another vector
     \a templat.

     The \a flag vector is changed to contain 1's and 0's only. A 1
     means that the corresponding element in the \a templat vector is
     valid and a zero means that the corresponding element is a NaN.

     \note Space for vector \a flag is reallocated to fit the size of
     vector \a templat if sizes mismatch.

     \return True if the \a templat vector contains at least one NaN.
  */
  bool nan(const vector& templat, vector& flag);

  /**
     \brief Transforms a vector to a basis vector.

     All elements are set to zero except the \a i-th element which is
     set to one.
  */
  void set_basis(vector&, size_t i);

  /**
     Randomly shuffles the elements in vector \a invec
  */
  void shuffle(vector& invec); 

  /**
     Sort the elements in the vector.

     \note Bug in GSL: if the vector contains one ore more NaNs the
     call never returns.
  */
  void sort(vector&);

  /**
     \brief Calculate the sum of all vector elements.

     \return The sum.
  */
  double sum(const vector&);

  /**
     \brief Swap vector elements by copying.

     The two vectors must have the same length.

     \throw GSL_error if vector lengths differs.
  */
  void swap(vector&, vector&);

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

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

#endif