source: trunk/c++_tools/utility/matrix.h @ 675

#ifndef _theplu_utility_matrix_
#define _theplu_utility_matrix_

// $Id: matrix.h 675 2006-10-10 12:08:45Z 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, 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 "yat/utility/vector.h"
#include "yat/utility/Exception.h"

#include <gsl/gsl_matrix.h>
#include <iostream>

namespace theplu {
namespace utility {

  class vector;

  ///
  /// This is the yat interface to GSL matrix. '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[Matrix views] GSL matrix views are supported and these are
  /// disguised as ordinary utility::matrix'es. A support function is
  /// added, utility::matrix::isview(), that can be used to check if a
  /// matrix object is a view. Note that view matricess do not own the
  /// undelying data, and a view is not valid if the matrix owning the
  /// data is deallocated.
  ///
  /// @note All GSL matrix related functions are not implement but
  /// most functionality defined for GSL matrices can be achieved with
  /// this interface class. Most notable GSL functionality not
  /// supported are no binary file support and views on arrays,
  /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and
  /// superdiagonals. If there is an interest from the user community,
  /// the omitted functionality could be included.
  ///
  /// @todo Maybe it would be smart to create temporary objects need
  /// for BLAS in constructors. As it is now temporary objects are
  /// called before BLAS functionality iss used, cf. const matrix&
  /// matrix::operator*=(const matrix& other)
  ///
  class matrix
  {
  public:

    /**
       @brief The default constructor.
    
       This contructor does not initialize underlying (essential)
       structures.
    */
    inline matrix(void) : m_(NULL), view_(NULL) {}

    ///
    /// Constructor. Allocates memory space for \a r times \a c
    /// elements, and sets all elements to \a init_value.
    ///
    inline matrix(const size_t& r, const size_t& c, double init_value=0)
      : view_(NULL) { m_ = gsl_matrix_alloc(r,c); set_all(init_value); }

    ///
    /// The copy constructor.
    ///
    /// @note If the object to be copied is a matrix view, the values
    /// of the view will be copied, i.e. the view is not copied.
    ///
    inline matrix(const matrix& o) : view_(NULL) {m_=o.create_gsl_matrix_copy();}

    ///
    /// The matrix view constructor.
    ///
    /// Create a view of matrix \a m, with starting row \a offset_row,
    /// starting column \a offset_column, row size \a n_row, and
    /// column size \a n_column.
    ///
    /// A matrix view can be used as any matrix 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
    /// matrix 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.
    ///
    matrix(matrix& m, size_t offset_row, size_t offset_column, size_t n_row,
           size_t n_column);

    ///
    /// The istream constructor.
    ///
    /// Matrix rows are sepearated with the new line character, and
    /// column elements in a row must be separated either with white space
    /// characters except the new line (\\n) character (default), or 
    /// by the delimiter \a sep.
    ///. Rows, and
    /// columns, are read consecutively and only rectangular matrices
    /// are supported. Empty lines are ignored. End of input is at end
    /// of file marker.
    ///
    //    explicit matrix(std::istream &) throw (utility::IO_error,std::exception);
    explicit matrix(std::istream &, char sep='\0') 
      throw (utility::IO_error,std::exception);


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

    ///
    /// Elementwise addition of the elements of matrix \a b to the
    /// elements of the calling matrix ,\f$ a_{ij} = a_{ij} + b_{ij} \;
    /// \forall i,j \f$. The result is stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int add(const matrix& b) { return gsl_matrix_add(m_,b.m_); }

    ///
    /// Add the scalar value \a d to the elements of the calling
    /// matrix, \f$ a_{ij} = a_{ij} + d \; \forall i,j \f$. The result
    /// is stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int
    add_constant(const double d) { return gsl_matrix_add_constant(m_,d); }

    ///
    /// @return The number of columns in the matrix.
    ///
    inline size_t columns(void) const { return m_->size2; }

    ///
    /// Elementwise division of the elemnts of the calling matrix by
    /// the elements of matrix \a b, \f$ a_{ij} = a_{ij} / b_{ij} \;
    /// \forall i,j \f$. The result is stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int
    div_elements(const matrix& b) { return gsl_matrix_div_elements(m_,b.m_); }

    ///
    /// Check whether matrices 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 matrix contain a NaN, false is always returned.
    ///
    bool equal(const matrix&, const double precision=0) const;

    ///
    /// @return A const pointer to the internal GSL matrix.
    ///
    inline const gsl_matrix* gsl_matrix_p(void) const { return m_; }

    ///
    /// @return A pointer to the internal GSL matrix.
    ///
    // Jari, is this needed?
    inline gsl_matrix* gsl_matrix_p(void) { return m_; };

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

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

    ///
    /// @return The maximum value of the matrix.
    ///
    inline double max(void) const { return gsl_matrix_max(m_); }

    ///
    /// @return The index to the element with the maximum value of the
    /// matrix. The lowest index has precedence (searching in
    /// row-major order).
    ///
    inline std::pair<size_t,size_t> max_index(void) const;

    ///
    /// @return The minimum value of the matrix.
    ///
    inline double min(void) const { return gsl_matrix_min(m_); }

    ///
    /// @return The index to the element with the minimum value of the
    /// matrix. The lowest index has precedence (searching in
    /// row-major order).
    ///
    inline std::pair<size_t,size_t> min_index(void) const;

    ///
    /// @return The indecies to the element with the minimum and
    /// maximum values of the matrix, respectively. The lowest index
    /// has precedence (searching in row-major order).
    ///
    inline void minmax_index(std::pair<size_t,size_t>& min,
                             std::pair<size_t,size_t>& max) const
      { gsl_matrix_minmax_index(m_,&min.first,&min.second,
                                &max.first,&max.second); }

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

    ///
    /// The argument is changed into a matrix containing 1's and 0's
    /// only. 1 means element in matrix is valid and a zero means
    /// element is a NaN.
    ///
    /// @return true if matrix contains at least one NaN.
    ///
    bool nan(matrix&) const;

    ///
    /// Multiply the elements of matrix \a b with the elements of the
    /// calling matrix ,\f$ a_{ij} = a_{ij} * b_{ij} \; \forall
    /// i,j \f$. The result is stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int
    mul_elements(const matrix& b) { return gsl_matrix_mul_elements(m_,b.m_); }

    ///
    /// @return The number of rows in the matrix.
    ///
    inline size_t rows(void) const { return m_->size1; }

    ///
    /// Multiply the elements of the calling matrix with a scalar \a
    /// d, \f$ a_{ij} = d * a_{ij} \; \forall i,j \f$. The result is
    /// stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int scale(const double d) { return gsl_matrix_scale(m_,d); }

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

    ///
    /// Set all elements to \a value.
    ///
    inline void set_all(const double value) { gsl_matrix_set_all(m_,value); }

    ///
    /// Set \a column values to values in \a vec.
    ///
    /// @return Whatever GSL returns.
    ///
    /// @note No check on size is done.
    ///
    inline int set_column(const size_t column, const vector& vec)
      { return gsl_matrix_set_col(m_, column, vec.gsl_vector_p()); }

    ///
    /// Set \a row values to values in \a vec.
    ///
    /// @return Whatever GSL returns.
    ///
    /// @note No check on size is done.
    ///
    inline int set_row(const size_t row, const vector& vec)
      { return gsl_matrix_set_row(m_, row, vec.gsl_vector_p()); }

    ///
    /// Subtract the elements of matrix \a b from the elements of the
    /// calling matrix ,\f$ a_{ij} = a_{ij} - b_{ij} \; \forall
    /// i,j \f$. The result is stored into the calling matrix.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int sub(const matrix& b) { return gsl_matrix_sub(m_,b.m_); }

    ///
    /// Exchange the elements of the this and \a other by copying. The
    /// two matrices must have the same size.
    ///
    /// @return Whatever GSL returns.
    ///
    inline int swap(matrix& other) { return gsl_matrix_swap(m_,other.m_); }

    ///
    /// Swap columns \a i and \a j.
    ///
    inline int swap_columns(const size_t i,const size_t j)
      { return gsl_matrix_swap_columns(m_,i,j); }

    ///
    /// Swap row \a i and column \a j.
    ///
    inline int swap_rowcol(const size_t i,const size_t j)
      { return gsl_matrix_swap_rowcol(m_,i,j); }

    ///
    /// Swap rows \a i and \a j.
    ///
    inline int swap_rows(const size_t i, const size_t j)
      { return gsl_matrix_swap_rows(m_,i,j); }

    ///
    /// Transpose the matrix.
    ///
    void transpose(void);

    ///
    /// @return Reference to the element position (\a row, \a column).
    ///
    inline double& operator()(size_t row,size_t column)
    { return (*gsl_matrix_ptr(m_,row,column)); }

    ///
    /// @return Const reference to the element position (\a row, \a
    /// column).
    ///
    inline const double& operator()(size_t row,size_t column) const
    { return (*gsl_matrix_const_ptr(m_,row,column)); }

    ///
    /// Matrix-vector multiplication.
    ///
    /// @return The resulting vector.
    ///
    // Jari, where should this go?
    const vector operator*(const vector&) const;

    ///
    /// Comparison operator.
    ///
    /// @return True if all elements are equal otherwise False.
    ///
    /// @see equal
    ///
    inline bool operator==(const matrix& other) const { return equal(other); }

    ///
    /// Comparison operator.
    ///
    /// @return False if all elements are equal otherwise True.
    ///
    /// @see equal
    ///
    inline bool operator!=(const matrix& other) const { return !equal(other); }

    ///
    /// The assignment operator. There is no requirements on
    /// dimensions, i.e. the matrix is remapped in memory if
    /// necessary. This implies that in general views cannot be
    /// assigned using this operator. Views will be mutated into
    /// normal matrices. 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 matrix.
    ///
    /// @see int set(const matrix&)
    ///
    const matrix& operator=(const matrix& other);

    ///
    /// Add and assign operator.
    ///
    inline const matrix& operator+=(const matrix& m) { add(m); return *this; }

    ///
    /// Add and assign operator.
    ///
    inline const matrix&
    operator+=(const double d) { add_constant(d); return *this; }

    ///
    /// Subtract and assign operator.
    ///
    inline const matrix& operator-=(const matrix& m) { sub(m); return *this; }

    ///
    /// Multiply and assigment operator.
    ///
    /// @return Const reference to the resulting matrix.
    ///
    const matrix& operator*=(const matrix&);

    ///
    /// Multiply and assign operator.
    ///
    inline const matrix& operator*=(const double d) { scale(d); return *this; }


  private:

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

    gsl_matrix* m_;
    gsl_matrix_view* view_;
  };

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


}} // of namespace utility and namespace theplu

#endif