mlMatrixTemplate.h

Go to the documentation of this file.

/*************************************************************************************
**
** Copyright 2012, MeVis Medical Solutions AG
**
** The user may use this file in accordance with the license agreement provided with
** the Software or, alternatively, in accordance with the terms contained in a
** written agreement between the user and MeVis Medical Solutions AG.
**
** For further information use the contact form at https://www.mevislab.de/contact
**
**************************************************************************************/
 
#ifndef ML_MATRIX_TEMPLATE_H
#define ML_MATRIX_TEMPLATE_H
 
 
 
#include <mlInitSystemML.h>
#include <mlPrintTemplateErrors.h>
 
#ifdef UNIX
#include "mlDebug.h"
#endif
 
#include <ThirdPartyWarningsDisable.h>
#include <iostream>
#include <valarray>
#include <numeric>
#include <iterator>
#include <ThirdPartyWarningsRestore.h>
 
//======================================================================
// namespaces
//======================================================================
ML_START_NAMESPACE
//======================================================================
// declaration of class *Slice_iter*
//======================================================================
template<class T>
class Slice_iter {
  std::valarray<T>* v;                            
  std::slice s;                                   
  size_t curr;                                    
 
  T& ref(size_t i) const { return (*v)[s.start()+i*s.stride()]; }
public:
  Slice_iter(std::valarray<T>* vv, std::slice ss) :v(vv), s(ss), curr(0) { }
 
  Slice_iter<T> end() const
  {
    Slice_iter<T> t = *this;
    t.curr = s.size();                            // index of one plus last element
    return t;
  }
 
  Slice_iter<T>& operator++() { curr++; return *this; }                         
  Slice_iter<T> operator++(int) { Slice_iter<T> t = *this; curr++; return t; }  
 
  T& operator[](size_t i) { return ref(i); }      
  T& operator()(size_t i) { return ref(i); }      
  T& operator*() { return ref(curr); }            
 
  // caution: friend function has to be implemented here inplace to support
  // explicit template specification according to ansi standart in gcc, i.e to
  // let the compiler know, that the friend functions are realy templated overloads
  // (note: the gcc conform syntax (friend bool operator== <> (...)) is not supported
  // on other compiler)
  friend bool operator== (const Slice_iter<T>& p, const Slice_iter<T>& q)   
  {
    return p.curr==q.curr
      && p.s.stride()==q.s.stride()
      && p.s.start()==q.s.start();
  }
 
  friend bool operator!= (const Slice_iter<T>& p, const Slice_iter<T>& q)   
  {
    return !(p==q);
  }
 
  friend bool operator< (const Slice_iter<T>& p, const Slice_iter<T>& q)   
  {
    return p.curr<q.curr
      && p.s.stride()==q.s.stride()
      && p.s.start()==q.s.start();
  }
};
 
 
//======================================================================
// declaration of class *Cslice_iter*
//======================================================================
template<class T>
class Cslice_iter
{
  std::valarray<T>* v;                            
  std::slice s;                                   
  size_t curr;                                    
 
  const T& ref(size_t i) const { return (*v)[s.start()+i*s.stride()]; }
public:
  Cslice_iter(std::valarray<T>* vv, std::slice ss): v(vv), s(ss), curr(0){}
 
  Cslice_iter<T> end() const
  {
    Cslice_iter<T> t = *this;
    t.curr = s.size();                            // index of one plus last element
    return t;
  }
 
  Cslice_iter<T>& operator++() { curr++; return *this; }                          
  Cslice_iter<T> operator++(int) { Cslice_iter<T> t = *this; curr++; return t; }  
 
  const T& operator[](size_t i) const { return ref(i); }  
  const T& operator()(size_t i) const { return ref(i); }  
  const T& operator*() const { return ref(curr); }        
 
  // caution: friend function has to be implemented here inplace to support
  // explicit template specification according to ansi standart in gcc, i.e to
  // let the compiler know, that the friend functions are realy templated overloads
  // (note: the gcc conform syntax (friend bool operator== <> (...)) is not supported
  // on other compiler)
  friend bool operator==(const Cslice_iter<T>& p, const Cslice_iter<T>& q) 
  {
    return p.curr==q.curr
      && p.s.stride()==q.s.stride()
      && p.s.start()==q.s.start();
  }
 
  friend bool operator!=(const Cslice_iter<T>& p, const Cslice_iter<T>& q) 
  {
    return !(p==q);
  }
 
  friend bool operator< (const Cslice_iter<T>& p, const Cslice_iter<T>& q) 
  {
    return p.curr<q.curr
      && p.s.stride()==q.s.stride()
      && p.s.start()==q.s.start();
  }
 
};
 
//======================================================================
// declaration of class *MatrixTemplate*
//======================================================================
template <class T>
class MatrixTemplate {
  std::valarray<T>* v;                            
  size_t size_x;                                  
  size_t size_y;                                  
public:
  MatrixTemplate(size_t x, size_t y);             
  MatrixTemplate(const MatrixTemplate<T>&);       
  MatrixTemplate<T>& operator=(const MatrixTemplate<T>&);   
  ~MatrixTemplate();                              
 
  void freeMatrix();                                   
  void resizeMatrix(size_t x, size_t y, const T& c = T());     
  void overrideMatrix(size_t x, size_t y, const T& c = T());   
 
  size_t size() const { return size_x*size_y; }   
  size_t sizeX() const { return size_x; }         
  size_t sizeY() const { return size_y; }         
 
  Slice_iter<T> row(size_t i);                    
  Cslice_iter<T> row(size_t i) const;             
 
  Slice_iter<T> column(size_t i);                 
  Cslice_iter<T> column(size_t i) const;          
 
  T& operator()(size_t x, size_t y);              
  T operator()(size_t x, size_t y) const;         
 
  Slice_iter<T> operator()(size_t i) { return column(i); }          
  Cslice_iter<T> operator()(size_t i) const { return column(i); }   
 
  Slice_iter<T> operator[](size_t i) { return column(i); }          
  Cslice_iter<T> operator[](size_t i) const { return column(i); }   
 
  MatrixTemplate<T>& operator*=(T);               
 
  std::valarray<T>& array() { return *v; }        
 
  operator T * () { return reinterpret_cast<T*>(&v); }           
};
 
 
//======================================================================
// declaration of class *MatrixSizedTemplate*
//======================================================================
template <class T, size_t SIZE_X, size_t SIZE_Y>
class MatrixSizedTemplate : public MatrixTemplate<T>
{
  public:
    MatrixSizedTemplate () : MatrixTemplate<T>(SIZE_X, SIZE_Y) {;}
};
 
 
 
 
 
 
//=======================================================================
// Implementation of MatrixTemplate - functions
//=======================================================================
 
template <class T>
inline Slice_iter<T> MatrixTemplate<T>::row(size_t i)
{
  return Slice_iter<T>(v,std::slice(i,size_x,size_y));
}
 
template <class T>
inline Cslice_iter<T> MatrixTemplate<T>::row(size_t i) const
{
  return Cslice_iter<T>(v,std::slice(i,size_x,size_y));
}
 
template <class T>
inline Slice_iter<T> MatrixTemplate<T>::column(size_t i)
{
  return Slice_iter<T>(v,std::slice(i*size_y,size_y,1));
}
 
template <class T>
inline Cslice_iter<T> MatrixTemplate<T>::column(size_t i) const
{
  return Cslice_iter<T>(v,std::slice(i*size_y,size_y,1));
}
 
template <class T>
MatrixTemplate<T>::MatrixTemplate(size_t x, size_t y)
{
  size_x = x;
  size_y = y;
  v = new std::valarray<T>(size_x*size_y);
}
 
template <class T>
MatrixTemplate<T>::MatrixTemplate(const MatrixTemplate<T>&mt)
{
  // set size according to transfered object
  size_x = mt.size_x;
  size_y = mt.size_y;
  // assign new valarray of appropriate size
  v = new std::valarray<T>(size_x*size_y);
 
  // copy memory content
  *v = *(mt.v);
}
 
template <class T>
MatrixTemplate<T>& MatrixTemplate<T>::operator=(const MatrixTemplate<T>&mt)
{
  if(this != &mt){
    // resize matrixTemplate
    this->resizeMatrix(mt.size_x, mt.size_y);
    // copy memory content
    *v = *(mt.v);
  }
  return *this;
}
 
template <class T>
MatrixTemplate<T>::~MatrixTemplate()
{
  delete v;
}
 
template <class T>
void MatrixTemplate<T>::resizeMatrix(size_t x, size_t y, const T& c)
{
  size_x = x;
  size_y = y;
  v->resize(size_x*size_y, c);
}
 
template <class T>
void MatrixTemplate<T>::overrideMatrix(size_t x, size_t y, const T& c)
{
  size_x = x;
  size_y = y;
  *v = c;
  v->resize(size_x*size_y, c);
}
 
 
template <class T>
void MatrixTemplate<T>::freeMatrix()
{
  size_x = 0;
  size_y = 0;
  v->free();
}
 
template <class T>
T& MatrixTemplate<T>::operator()(size_t x, size_t y)
{
  return column(x)[y];
}
 
 
 
//---------------------------------------------------------------------------------
//---------------------------------------------------------------------------------
template <class T>
T mul(const Cslice_iter<T>& v1, const std::valarray<T>& v2)
{
  T res = 0;
  for (size_t i = 0; i<v2.size(); i++) res+= v1[i]*v2[i];
  return res;
}
 
 
template <class T>
std::valarray<T> operator*(const MatrixTemplate<T>& m, const std::valarray<T>& v)
{
  if (m.sizeX()!=v.size()){
    printTemplateError("mlMatrixTemplate.h, operator*() # 1",
                   ML_BAD_PARAMETER,
                   "Wrong number of elements in m*v\n, ignoring problem");
  }
 
  std::valarray<T> res(m.sizeY());
  for (size_t i = 0; i<m.sizeY(); i++) res[i] = mul(m.row(i),v);
  return res;
}
 
 
//std::valarray<T> operator*(const MatrixTemplate<T>& m, std::valarray<T>& v)
template <class T>
std::valarray<T> mul_mv(const MatrixTemplate<T>& m, std::valarray<T>& v)
{
  if (m.sizeX()!=v.size()){
    printTemplateError("mlMatrixTemplate.h, mul_mv()",
                   ML_BAD_PARAMETER,
                   "Wrong number of elements in m*v\n, ignoring problem");
  }
 
  std::valarray<T> res(m.sizeY());
 
  for (size_t i = 0; i<m.sizeY(); i++) {
    const Cslice_iter<T>& ri = m.row(i);
    res[i] = inner_product(ri,ri.end(),&v[0],T(0));
  }
  return res;
}
 
 
 
template <class T>
std::valarray<T> operator*(std::valarray<T>& v, const MatrixTemplate<T>& m)
{
  if (v.size()!=m.sizeY()){
    printTemplateError("mlMatrixTemplate.h, operator*() # 2",
                   ML_BAD_PARAMETER,
                   "Wrong number of elements in v*m\n, ignoring problem");
  }
 
  std::valarray<T> res(m.sizeX());
 
  for (size_t i = 0; i<m.sizeX(); i++) {
    const Cslice_iter<T>& ci = m.column(i);
    res[i] = inner_product(ci,ci.end(),&v[0],T(0));
  }
  return res;
}
 
template <class T>
MatrixTemplate<T>& MatrixTemplate<T>::operator*=(T d)
{
  (*v) *= d;
  return *this;
}
 
template <class T>
std::ostream& operator<<(std::ostream& os, MatrixTemplate<T>& m)
{
  for(int y=0; y<m.sizeY(); y++)
  {
    for(int x=0; x<m.sizeX(); x++)
      os<<m[x][y]<<"\t";
    os << "\n";
  }
  return os;
}
 
 
 
 
 
} // namespace ml
 
#endif //__MatrixTemplate_h__