source: trunk/yat/utility/VectorBase.h @ 1027

#ifndef _theplu_yat_utility_vector_base_
#define _theplu_yat_utility_vector_base_

// $Id: VectorBase.h 1027 2008-02-02 21:29:29Z peter $

/*
  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
  Copyright (C) 2005 Jari Häkkinen, Markus Ringnér, Peter Johansson
  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér
  Copyright (C) 2007 Jari Häkkinen, Peter Johansson
  Copyright (C) 2008 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 "Exception.h"
#include "Iterator.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;
  class vector;

  /**
     @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 VectorBase
  {
  public:
    /// \brief VectorBase::const_iterator
    typedef Iterator<const double&, const VectorBase> const_iterator;

    /**
       \brief Constructor.
    */
    VectorBase(const gsl_vector* v=NULL);

    ///
    /// The destructor.
    ///
    virtual ~VectorBase(void);

    /**
       \return read-only iterator to start of VectorBase
     */
    const_iterator begin(void) const;

    /**
       \return read-only iterator to end of VectorBase
     */
    const_iterator end(void) const;

    /**
       \brief Check whether VectorBases 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 VectorBase contain a NaN, false is always returned.

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

    ///
    /// @return A const pointer to the internal GSL vector,
    ///
    const gsl_vector* gsl_vector_p(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.
     */
    virtual bool isview(void) const=0; 

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

    /**
       \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;
    // Peter, remove this one
    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 VectorBase&, const double precision=0)
    */
    bool operator==(const VectorBase&) 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 VectorBase&, const double precision=0)
    */
    bool operator!=(const VectorBase&) const;

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

  protected:
    const gsl_vector* const_vec_;

  private:
    // copy assignment no allowed
    const VectorBase& operator=(const VectorBase&);
  };

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  /**
     Create a vector \a sort_index containing the indeces of
     elements in a another VectorBase \a invec.  The elements of \a
     sort_index give the index of the VectorBase element which would
     have been stored in that position if the VectorBase had been sorted
     in place. The first element of \a sort_index gives the index of the least
     element in \a invec, and the last element of \a sort_index gives the 
     index of the greatest element in \a invec . The VectorBase \a invec 
     is not changed.
  */
  void sort_index(std::vector<size_t>& sort_index, const VectorBase& invec);

  /** 
      Similar to sort_index but creates a VectorBase with indices to the \a k
  smallest elements in \a invec.  
  */
  void sort_smallest_index(std::vector<size_t>& sort_index, size_t k, 
                           const VectorBase& invec);

  /** 
      Similar to sort_index but creates a VectorBase with indices to the \a k
  largest elements in \a invec.  
  */
  void sort_largest_index(std::vector<size_t>& sort_index, size_t k, 
                          const VectorBase& invec);

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

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

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

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

#endif