JMatrixND.hh

Go to the documentation of this file.

#ifndef __JMATH__JMATRIXND__
#define __JMATH__JMATRIXND__

#include <limits>
#include <cmath>
#include <algorithm>
#include <ostream>
#include <iomanip>
#include <vector>
#include <utility>

#include "JIO/JSerialisable.hh"
#include "JLang/JEquals.hh"
#include "JLang/JException.hh"
#include "JLang/JSingleton.hh"
#include "Jeep/JPrint.hh"
#include "JMath/JMath.hh"


/**
 * \author mdejong
 */

namespace JMATH {}
namespace JPP { using namespace JMATH; }

namespace JMATH {

  using JIO::JReader;
  using JIO::JWriter;
  using JLANG::JEquals;
  using JLANG::JDivisionByZero;
  using JLANG::JException;
  using JLANG::JSingleton;


  /**
   * Type definition of a NxN matrix.
   */
  struct JMatrixND_t :
    public std::vector< std::vector<double> >
  {

    typedef std::vector< std::vector<double> >       matrix_type;

    typedef matrix_type::const_iterator              const_row_type;
    typedef std::vector<double>::const_iterator      const_col_type;

    typedef matrix_type::iterator                    row_type;
    typedef std::vector<double>::iterator            col_type;
  };


  /**
   * Type definition of permutation matrix.
   */
  typedef std::vector< std::pair<int,int> >          JPermutationMatrix;
  
    
  /**
   * N x N matrix.
   */
  class JMatrixND :
    public JMatrixND_t,
    public JMath      <JMatrixND>,
    public JEquals    <JMatrixND>
  {
  public:

    using JMath<JMatrixND>::mul;


    /**
     * Default constructor.
     */
    JMatrixND()
    {}


    /**
     * Constructor.
     *
     * \param  size     dimension
     */
    JMatrixND(const size_t size) 
    {
      resize(size);
      reset();
    }


    /**
     * Resize matrix.
     * Note that this method resets the matrix.
     *
     * \param  size     dimension
     */
    void resize(const size_t size)
    {
      matrix_type::resize(size);

      for (row_type row = this->begin(); row != this->end(); ++row) {

        row->resize(size);

        for (col_type col = row->begin(); col != row->end(); ++col) {
          *col = 0.0;
        }
      }
    }


    /**
     * Get matrix element.
     *
     * \param  row      row    number
     * \param  col      column number
     * \return          matrix element at (row,col)
     */
    double at(size_t row, size_t col) const
    {
      return (*this)[row][col];
    }


    /**
     * Get matrix element.
     *
     * \param  row      row    number
     * \param  col      column number
     * \return          matrix element at (row,col)
     */
    double& at(size_t row, size_t col)
    {
      return (*this)[row][col];
    }


    /**
     * Get matrix element.
     *
     * \param  row      row    number
     * \param  col      column number
     * \return          matrix element at (row,col)
     */
    double operator()(size_t row, size_t col) const
    {
      return at(row, col);
    }


    /**
     * Get matrix element.
     *
     * \param  row      row    number
     * \param  col      column number
     * \return          matrix element at (row,col)
     */
    double& operator()(size_t row, size_t col)
    {
      return at(row, col);
    }


    /**
     * Set matrix to the null matrix.
     *
     * \return          this matrix
     */
    JMatrixND& reset()
    {
      for (row_type row = this->begin(); row != this->end(); ++row) {
        for (col_type col = row->begin(); col != row->end(); ++col) {
          *col = 0.0;
        }
      }

      return *this;
    }


    /**
     * Set to identity matrix
     *
     * \return          this matrix
     */
    JMatrixND& setIdentity()
    {
      for (size_t i = 0; i != this->size(); ++i) {

        at(i,i) = 1.0;

        for (size_t j = 0; j != i; ++j) {
          at(i,j) = at(j,i) = 0.0;
        }
      }

      return *this;
    }


    /**
     * Transpose.
     *
     * \return          this matrix
     */
    JMatrixND& transpose()
    {
      for (size_t i = 0; i != this->size(); ++i) {
        for (size_t j = 0; j != i; ++j) {
          std::swap(at(i,j), at(j,i));
        }
      }

      return *this;
    }


    /**
     * Negate matrix.
     *
     * \return          this matrix
     */
    JMatrixND& negate()
    {
      for (row_type row = this->begin(); row != this->end(); ++row) {
        for (col_type col = row->begin(); col != row->end(); ++col) {
          *col = -*col;
        }
      }

      return *this;
    }


    /**
     * Matrix addition.
     *
     * \param  A        matrix
     * \return          this matrix
     */
    JMatrixND& add(const JMatrixND& A)
    {
      if (this->size() == A.size()) {

        for (size_t i = 0; i != this->size(); ++i) {
          for (size_t j = 0; j != this->size(); ++j) {
            at(i,j) += A.at(i,j);
          }
        }

      } else {
        THROW(JException, "JMatrixND::add() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
      }

      return *this;
    }


    /**
     * Matrix subtraction.
     *
     * \param  A        matrix
     * \return          this matrix
     */
    JMatrixND& sub(const JMatrixND& A)
    {
      if (this->size() == A.size()) {

        for (size_t i = 0; i != this->size(); ++i) {
          for (size_t j = 0; j != this->size(); ++j) {
            at(i,j) -= A.at(i,j);
          }
        }

      } else {
        THROW(JException, "JMatrixND::sub() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
      }

      return *this;
    }


    /**
     * Scale matrix.
     *
     * \param  factor   factor
     * \return          this matrix
     */
    JMatrixND& mul(const double factor)
    {
      for (row_type row = this->begin(); row != this->end(); ++row) {
        for (col_type col = row->begin(); col != row->end(); ++col) {
          *col *= factor;
        }
      }

      return *this;
    }


    /**
     * Scale matrix.
     *
     * \param  factor   factor
     * \return          this matrix
     */
    JMatrixND& div(const double factor)
    {
      for (row_type row = this->begin(); row != this->end(); ++row) {
        for (col_type col = row->begin(); col != row->end(); ++col) {
          *col /= factor;
        }
      }

      return *this;
    }


    /**
     * Matrix multiplication.
     *
     * \param  A        matrix
     * \param  B        matrix
     * \return          this matrix
     */
    JMatrixND& mul(const JMatrixND& A,
                   const JMatrixND& B)
    {
      if (A.size() == B.size()) {
        
        this->resize(A.size());
        
        for (size_t i = 0; i != A.size(); ++i) {
          for (size_t j = 0; j != A.size(); ++j) {
            
            this->at(i,j) = 0.0;
            
            for (size_t k = 0; k != A.size(); ++k) {
              this->at(i,j) += A[i][k] * B[k][j];
            }
          }
        }
        
      } else {

        THROW(JException, "JMatrixND::mul() inconsistent matrix dimensions " << A.size() << ' ' << B.size());
      }
      
      return *this;
    }


    /**
     * Equality.
     *
     * \param  A        matrix
     * \param  eps      numerical precision
     * \return          true if matrices identical; else false
     */
    bool equals(const JMatrixND& A,
                const double     eps = std::numeric_limits<double>::min()) const
    {
      if (this->size() == A.size()) {

        for (const_row_type row1 = this->begin(), row2 = A.begin(); row1 != this->end(); ++row1, ++row2) {
          for (const_col_type col1 = row1->begin(), col2 = row2->begin(); col1 != row1->end(); ++col1, ++col2) {
            if (fabs(*col1 - *col2) > eps) {
              return false;
            }
          }
        }

        return true;

      } else {
        THROW(JException, "JMatrixND::equals() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
      }
    }
      

    /**
     * Test identity.
     *
     * \param  eps      numerical precision
     * \return          true if identity matrix; else false
     */
    bool isIdentity(const double eps = std::numeric_limits<double>::min()) const
    {
      for (size_t i = 0; i != this->size(); ++i) {

        if (fabs(1.0 - at(i,i)) > eps) {
          return false;
        };
        
        for (size_t j = 0; j != i; ++j) {
          if (fabs(at(i,j)) > eps || fabs(at(j,i)) > eps) {
            return false;
          };
        }
      }

      return true;
    }


    /**
     * Get determinant of matrix.
     *
     * \return         determinant of matrix
     */
    double getDeterminant() const
    {
      JMatrixHelper& A = JMatrixHelper::getInstance();

      A.decompose(*this);
    
      double det = 1.0;
      
      for (size_t i = 0; i != A.size(); ++i) {
        det *= A[i][i];
      }
      
      return det * A.getSign();
    }


    /**
     * Test positive definiteness.
     *
     * \return         true if positive definite; else false
     */
    bool isPositiveDefinite() const
    {
      JMatrixHelper& A = JMatrixHelper::getInstance();

      A.decompose(*this);
    
      for (size_t i = 0; i != A.size(); ++i) {
        if (A[i][i] <= 0.0) {
          return false;
        }
      }

      return true;
    }


    /**
     * Test positive semi-definiteness.
     *
     * \return         true if positive semi-definite; else false
     */
    const bool isPositiveSemiDefinite() const
    {
      JMatrixHelper& A = JMatrixHelper::getInstance();

      A.decompose(*this);    

      for (size_t i = 0; i != A.size(); ++i) {
        if (A[i][i] < 0.0) {
          return false;
        }
      }

      return true;
    }


    /**
     * Read matrix from input.
     *
     * \param  in       reader
     * \param  matrix   matrix
     * \return          reader
     */
    friend inline JReader& operator>>(JReader& in, JMatrixND& matrix)
    {
      int size;

      in >> size;

      matrix.resize(size);

      for (row_type row = matrix.begin(); row != matrix.end(); ++row) {
        for (col_type col = row->begin(); col != row->end(); ++col) {
          in >> *col;
        }
      }

      return in;
    }


    /**
     * Write matrix to output.
     *
     * \param  out      writer
     * \param  matrix   matrix
     * \return          writer
     */
    friend inline JWriter& operator<<(JWriter& out, const JMatrixND& matrix)
    {
      int size = matrix.size();

      out << size;

      for (const_row_type row = matrix.begin(); row != matrix.end(); ++row) {
        for (const_col_type col = row->begin(); col != row->end(); ++col) {
          out << *col;
        }
      }

      return out;
    }


    /**
     * Print ASCII formatted output.
     *
     * \param  out      output stream
     * \param  A        matrix
     * \return          output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JMatrixND& A)
    {
      using namespace std;

      JFlags flags(out);

      for (const_row_type row = A.begin(); row != A.end(); ++row) {

        for (const_col_type col = row->begin(); col != row->end(); ++col) {
          out << FIXED(10,2) << *col << ' ';
        }

        out << endl;
      }

      return out;
    }

  protected:    
    /**
     * Auxiliary class for matrix decomposition.
     */
    struct JMatrixHelper :
      public JMatrixND_t,
      public JSingleton<JMatrixHelper>
    {
      /**
       * Set matrix helper.
       *
       * \param  A        matrix
       * \return          this matrix helper
       */
      const JMatrixHelper& set(const JMatrixND& A)
      {
        static_cast<JMatrixND_t&>(*this) = A;

        return *this;
      }
      
      
      /**
       * LU decomposition.
       *
       * \param  A        matrix
       */
      void decompose(const JMatrixND& A)
      {
        set(A);

        sign = +1;
        
        size_t   k;
        double   val, del;
        row_type row, root;
        col_type col, cursor;
      
        // LU decomposition

        for (root = this->begin(), k = 0; root != this->end(); ++root, ++k) {

          row_type pivot;

          val   = (*root)[k];
          pivot = root;

          for (row = root; ++row != this->end(); ) {
            if (fabs((*row)[k]) > fabs(val)) {
              val   = (*row)[k];
              pivot = row;
            }
          }

          if (val == 0) {
            throw JDivisionByZero("LU decomposition zero pivot");
          }

          if (pivot != root) {

            iter_swap(root, pivot);

            sign = -sign;
          }


          // decompose in LU

          for (row = root; ++row != this->end(); ) {
          
            del = (*row)[k] /= val;
          
            for (col = row->begin() + k, cursor = root->begin() + k; ++col != row->end(); )
              (*col) -= del * (*(++cursor));
          }
        }
      }


      /**
       * Get sign corresponding to permutation matrix.
       */
      int getSign() const
      {
        return sign;
      }
    
    private:
      int sign;
    };
  };


  /**
   * Determine dot product <tt>Y^T x V x Y</tt>, where
   *    - <tt>Y</tt> is a <tt>1xN</tt> matrix; and 
   *    - <tt>V</tt> is a <tt>NxN</tt> matrix. 
   *
   * \param  __begin          begin of data
   * \param  __end            end   of data
   * \param  V                matrix
   * \return                  dot product
   */
  template<class T>
  inline double getDot(T __begin,
                       T __end,
                       const JMatrixND& V)
  {
    double dot = 0.0;

    T yl = __begin;
    
    for (JMatrixND::const_row_type row = V.begin(); row != V.end(); ++row, ++yl) {

      T yr = __begin;

      for (JMatrixND::const_col_type col = row->begin(); col != row->end(); ++col, ++yr) {
        dot += (*yl) * (*col) * (*yr);
      }
    }
    
    return dot;
  }
}

#endif