Matrix.h

Go to the documentation of this file.

#ifndef __MATRIX_H
#define __MATRIX_H

#include <valarray>
#include<numeric>

namespace Teem
{
    template<typename T> class SliceIter;
    
    template<typename T> bool operator==(const SliceIter<T>& p, const SliceIter<T>& q)
    {
        return p.curr == q.curr && p.s.stride() == q.s.stride() && p.s.start() == q.s.start();
    }
    
    template<typename T> bool operator!=(const SliceIter<T>& p, const SliceIter<T>& q)
    {
        return !(p==q);
    }
    
    template<typename T> bool operator<(const SliceIter<T>& p, const SliceIter<T>& q)
    {
        return p.curr < q.curr && p.s.stride() == q.s.stride() && p.s.start() == q.s.start();
    }
    
    template<typename T> class SliceIter
    {
    public:
        SliceIter(std::valarray<T> *vv, const std::slice &ss) : v(vv), s(ss), curr(0) { }
        
        SliceIter end() const
        {
            SliceIter t = *this;
            t.curr = s.size();
            return t;
        }
        
        SliceIter& operator++() { curr++; return *this; }
        SliceIter operator++(int) { SliceIter 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); }
        
        friend bool operator== <>(const SliceIter& p, const SliceIter& q);
        friend bool operator!= <>(const SliceIter& p, const SliceIter& q);
        friend bool operator< <>(const SliceIter& p, const SliceIter& q);
        
    protected:
        std::valarray<T> *v; 
        const std::slice s; 
        size_t curr; 
        
        T& ref(size_t i) const { return (*v)[s.start() + i * s.stride()]; }
    };
    
    
    template<typename T> class ConstSliceIter;
    
    template<typename T> bool operator==(const ConstSliceIter<T>& p, const ConstSliceIter<T>& q)
    {
        return p.curr == q.curr && p.s.stride() == q.s.stride() && p.s.start() == q.s.start();
    }
    
    template<typename T> bool operator!=(const ConstSliceIter<T>& p, const ConstSliceIter<T>& q)
    {
        return !(p==q);
    }
    
    template<typename T> bool operator<(const ConstSliceIter<T>& p, const ConstSliceIter<T>& q)
    {
        return p.curr < q.curr && p.s.stride() == q.s.stride() && p.s.start() == q.s.start();
    }
    
    template<typename T> class ConstSliceIter
    {
    public:
        ConstSliceIter(const std::valarray<T> *vv, const std::slice &ss) : v(vv), s(ss), curr(0) { }
        
        ConstSliceIter end() const
        {
            ConstSliceIter t = *this;
            t.curr = s.size();
            return t;
        }
        
        ConstSliceIter& operator++() { curr++; return *this; }
        ConstSliceIter operator++(int) { ConstSliceIter 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); }
        
        friend bool operator== <>(const ConstSliceIter& p, const ConstSliceIter& q);
        friend bool operator!= <>(const ConstSliceIter& p, const ConstSliceIter& q);
        friend bool operator< <>(const ConstSliceIter& p, const ConstSliceIter& q);
        
    protected:
        const std::valarray<T> *v; 
        const std::slice s; 
        size_t curr; 
        
        const T& ref(size_t i) const { return (*v)[s.start() + i * s.stride()]; }
    };
    
    template<typename T> class Matrix
    {
    public:
        Matrix(size_t nx = 0, size_t ny = 0, T c = T()) : items(c, nx*ny), xDim(nx), yDim(ny) { }
        ~Matrix() { }
        
        size_t size() const { return xDim * yDim; }
        size_t columnNum() const { return xDim; }
        size_t rowNum() const { return yDim; }
        
        std::valarray<T>& flat() { return items; }
        
        const std::valarray<T>& const_flat() const { return items; }
        
        void resize(size_t nx, size_t ny, T c = T()) { items.resize(nx*ny, c); xDim = nx; yDim = ny; }
        
        SliceIter<T> row(size_t i) { return SliceIter<T>(&items, std::slice(i, xDim, yDim)); }
        ConstSliceIter<T> row(size_t i) const { return ConstSliceIter<T>(&items, std::slice(i, xDim, yDim)); }
        
        SliceIter<T> column(size_t i) { return SliceIter<T>(&items, std::slice(i*yDim, yDim, 1)); }
        ConstSliceIter<T> column(size_t i) const { return ConstSliceIter<T>(&items, std::slice(i*yDim, yDim, 1)); }
        
        T& operator() (size_t x, size_t y) { return column(x)[y]; }
        const T& operator() (size_t x, size_t y) const { return column(x)[y]; }
        
        SliceIter<T> operator() (size_t i) { return column(i); }
        ConstSliceIter<T> operator() (size_t i) const { return column(i); }
        
    protected:
        std::valarray<T> items; 
        size_t xDim; 
        size_t yDim; 
    };
    
    template<typename T> T operator*(const ConstSliceIter<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<typename T> std::valarray<T> operator*(const Matrix<T>& m, const std::valarray<T>& v)
    {
        assert(m.columnNum() == v.size());
        
        std::valarray<T> res(m.rowNum());
        for(size_t i = 0; i < m.rowNum(); i++)
            res[i] = m.row(i) * v;
        return res;
    }
    
    template<typename T> std::valarray<T> operator*(const std::valarray<T>& v, const Matrix<T>& m)
    {
        assert(m.rowNum() == v.size());
        
        std::valarray<T> res(m.columnNum());
        for(size_t i = 0; i < m.columnNum(); i++)
            res[i] = m.column(i) * v;
        return res;
    }
}

#endif

---