matrix_d.h

Go to the documentation of this file.

/* Copyright (C) 2000-2013 CEA
 *
 * This software and supporting documentation were developed by
 *     bioPICSEL
 *     CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
 *     18, route du Panorama
 *     92265 Fontenay-aux-Roses
 *     France
 */

#ifndef PRIMATOLOGIST_MATH_MATRIX_D_H
#define PRIMATOLOGIST_MATH_MATRIX_D_H

#include <primatologist/math/matrix.h>
#include <aims/data/data_g.h>                             // AimsData
#include <aims/math/gausslu.h>                            // AimsInversionLU

namespace aims {
namespace math {

  //==========================================================================
  //   BASE CLASS FOR MATRICES
  //==========================================================================

  //--------------------------------------------------------------------------
  //   CONSTRUCTOR / DESTRUCTOR / COPY
  //--------------------------------------------------------------------------

  template <typename T>
  MatrixBase<T>::MatrixBase( int nrow, int ncol ):
    carto::VolumeRef<T>( nrow, ncol )
  {
    this->fill(0);
  }

  template <typename T>
  template <typename U>
  MatrixBase<T>::MatrixBase( const carto::VolumeRef<U> & volume ):
    carto::VolumeRef<T>( volume.template copy<T>() )
  {}


  template <typename T>
  template <typename U>
  MatrixBase<T>::MatrixBase( const std::vector<U> & vector ):
    carto::VolumeRef<T>( vector.size() )
  {
    for( int i = 0; i < vector.size(); ++i )
      this->at(i) = vector[i];
  }

  template <typename T>
  MatrixBase<T>::~MatrixBase()
  {}

  // Utility structure to compare adresses of volumes of possibly
  // different types
  namespace internal {

    template <typename T, typename U>
    struct CompareVolumeRefAdress
    {
      static bool equal( const carto::VolumeRef<T> *,
                         const carto::VolumeRef<U> * )
      {
        return false;
      }
    };

    template <typename T>
    struct CompareVolumeRefAdress<T, T>
    {
      static bool equal( const carto::VolumeRef<T> * a,
                         const carto::VolumeRef<T> * b )
      {
        return( a == b );
      }
    };

  } // namespace internal

  template <typename T>
  template <typename U>
  MatrixBase<T> & MatrixBase<T>::operator= ( const carto::VolumeRef<U> & volume )
  {
    if( !internal::CompareVolumeRefAdress<T, U>::equal( this, &volume ) )
    {
      *this = volume.copy();
    }
    return *this;
  }

  template <typename T>
  template <typename U>
  MatrixBase<T> & MatrixBase<T>::operator= ( const std::vector<U> & vector )
  {
    this->reallocate( vector.size() );
    for( int i = 0; i < vector.size(); ++i )
      this->at(i) = vector[i];
    return *this;
  }

  //--------------------------------------------------------------------------
  //   PRIVATE SHARED MEMORY CONSTRUCTOR
  //--------------------------------------------------------------------------

  template <typename T>
  MatrixBase<T>::MatrixBase( carto::VolumeRef<T> & volume ):
    carto::VolumeRef<T>( volume )
  {}

  //--------------------------------------------------------------------------
  //   CONVERSIONS
  //--------------------------------------------------------------------------

  template <typename T>
  template <typename U>
  MatrixBase<T>::operator MatrixBase<U>() const
  {
    return( MatrixBase<U>( *this ) );
  }

  //--------------------------------------------------------------------------
  //   MATRIX OPERATIONS
  //--------------------------------------------------------------------------

  template <typename T>
  template <typename U>
  MatrixBase<T> & MatrixBase<T>::operator*= ( const MatrixBase<U> & m )
  {
    *this = *this * m;;
    return *this;
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::multiplies<T,U>::result_type>
  operator* ( const MatrixBase<T> & a, const MatrixBase<U> & b )
  {
    typedef typename carto::volumeutil::multiplies<T,U>::result_type
            result_type;
    if( a.ncol() != b.nrow() )
      throw incompatible_matrix_exception( a.ncol(), b.nrow() );

    MatrixBase<result_type> result( a.nrow(), b.ncol() );
    double value;
    for( int i = 0; i < a.nrow(); ++i )
      for( int j = 0; j < b.ncol(); ++j )
      {
        value = 0.;
        for( int k = 0; k < a.ncol(); ++k )
          value += a(i, k) * b(k, j);
        result(i, j) = value;
      }
    return result;
  }

  template  <typename T>
  MatrixBase<T> & MatrixBase<T>::transpose()
  {
    *this = aims::math::transpose(*this);
    return *this;
  }

  template  <typename T>
  MatrixBase<T> transpose( const MatrixBase<T> & m )
  {
    MatrixBase<T> result( m.ncol(), m.nrow() );
    for( int i = 0; i < m.nrow(); ++i )
      for( int j = 0; j < m.ncol(); ++j )
        result( j, i ) = m( i, j );
    return result;
  }

  template  <typename T>
  MatrixBase<T> & MatrixBase<T>::invert()
  {
    *this = aims::math::invert(*this);
    return *this;
  }

  template  <typename T>
  MatrixBase<T> invert( const MatrixBase<T> & matrix )
  {
    // We use AIMS tools
    AimsData<float> aimsFloat( matrix.template copy<float>() );
    AimsData<float> invertedAimsFloat = AimsInversionLU( aimsFloat );
    carto::VolumeRef<T> volume = carto::VolumeRef<float>( invertedAimsFloat.volume() ).copy<T>();
    return asMatrix( volume );
  }

  //==========================================================================
  //   OVERLOAD OPERATORS
  //==========================================================================

  template <typename T, typename U>
  bool operator== ( const MatrixBase<T> & vol, const MatrixBase<U> & other )
  {
    if( vol.nrow() != other.nrow() || vol.ncol() != other.ncol() )
      return false;
    for( int i = 0; i < vol.nrow(); ++i )
      for( int j = 0; j < vol.ncol(); ++j )
        if( vol(i, j) != other(i, j) )
          return false;
    return true;
  }

  template <typename T, typename U>
  bool operator!= ( const MatrixBase<T> & vol, const MatrixBase<U> & other )
  {
    return !( vol == other );
  }

  template <typename T>
  MatrixBase<T> operator- ( const MatrixBase<T> & vol )
  {
    carto::VolumeRef<T> volresult = -(carto::VolumeRef<T> &)vol;
    return asMatrix(volresult);
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::plus<T,U>::result_type >
  operator+ ( const MatrixBase<T> & vol, const U & value )
  {
    carto::VolumeRef<typename carto::volumeutil::plus<T,U>::result_type> volresult = (carto::VolumeRef<T>&)vol + value;
    return asMatrix(volresult);
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::minus<T,U>::result_type >
  operator- ( const MatrixBase<T> & vol, const U & value )
  {
    carto::VolumeRef<typename carto::volumeutil::minus<T,U>::result_type> volresult = (carto::VolumeRef<T>&)vol - value;
    return asMatrix(volresult);
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::multiplies<T,U>::result_type >
  operator* ( const MatrixBase<T> & vol, const U & value )
  {
    carto::VolumeRef<typename carto::volumeutil::multiplies<T,U>::result_type> volresult = (carto::VolumeRef<T>&)vol * value;
    return asMatrix(volresult);
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::divides<T,U>::result_type >
  operator/ ( const MatrixBase<T> & vol, const U & value )
  {
    carto::VolumeRef<typename carto::volumeutil::divides<T,U>::result_type> volresult = (carto::VolumeRef<T>&)vol / value;
    return asMatrix(volresult);
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::plus<U,T>::result_type >
  operator+ ( const U & value, const MatrixBase<T> & vol )
  {
    carto::VolumeRef<typename carto::volumeutil::plus<T,U>::result_type> volresult = value + (carto::VolumeRef<T>&)vol;
    return asMatrix(volresult);
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::minus<U,T>::result_type >
  operator- ( const U & value, const MatrixBase<T> & vol )
  {
    carto::VolumeRef<typename carto::volumeutil::minus<T,U>::result_type> volresult = value - (carto::VolumeRef<T>&)vol;
    return asMatrix(volresult);
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::multiplies<U,T>::result_type >
  operator* ( const U & value, const MatrixBase<T> & vol )
  {
    carto::VolumeRef<typename carto::volumeutil::multiplies<T,U>::result_type> volresult = value * (carto::VolumeRef<T>&)vol;
    return asMatrix(volresult);
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::divides<U,T>::result_type >
  operator/ ( const U & value, const MatrixBase<T> & vol )
  {
    carto::VolumeRef<typename carto::volumeutil::divides<T,U>::result_type> volresult = value / (carto::VolumeRef<T>&)vol;
    return asMatrix(volresult);
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::plus<T,U>::result_type >
  operator+ ( const MatrixBase<T> & vol, const MatrixBase<U> & other )
  {
    carto::VolumeRef<typename carto::volumeutil::plus<T,U>::result_type> volresult = (carto::VolumeRef<T>&)vol + (carto::VolumeRef<U>&)other;
    return asMatrix(volresult);
  }

  template <typename T, typename U>
  MatrixBase<typename carto::volumeutil::minus<T,U>::result_type >
  operator- ( const MatrixBase<T> & vol, const MatrixBase<U> & other )
  {
    carto::VolumeRef<typename carto::volumeutil::minus<T,U>::result_type> volresult = (carto::VolumeRef<T>&)vol - (carto::VolumeRef<U>&)other;
    return asMatrix(volresult);
  }

  //--------------------------------------------------------------------------
  //   SIZE
  //--------------------------------------------------------------------------

  template  <typename T>
  int MatrixBase<T>::nrow() const
  {
    return this->getSizeX();
  }

  template  <typename T>
  int MatrixBase<T>::ncol() const
  {
    return this->getSizeY();
  }

  template  <typename T>
  int MatrixBase<T>::size() const
  {
    return this->getSizeX() * this->getSizeY();
  }

  //==========================================================================
  //   STREAM / WRITE
  //==========================================================================

  template <typename T>
  std::ostream & operator<< ( std::ostream & out,
                              const MatrixBase<T> & matrix )
  {
    return out << (carto::VolumeRef<T> &)matrix;
  }

  // template <typename T>
  // void write( const MatrixBase<T> & matrix  );

  //==========================================================================
  //   UTILITY FUNCTIONS
  //==========================================================================

  template <typename T>
  MatrixBase<T> asMatrix( carto::VolumeRef<T> & volume )
  {
    return MatrixBase<T>( volume );
  }

} // namespace math
} // namespace aims

#endif // PRIMATOLOGIST_MATH_MATRIX_D_H