gsWeightMapper.h

Go to the documentation of this file.

 
#pragma once
 
#include <gsCore/gsExport.h>
#include <gsCore/gsLinearAlgebra.h>
 
namespace gismo
{
 
template<class T>
class gsWeightMapper
{
public:
    typedef T weightType;
    typedef index_t indexType; //indices of gsMatrix
 
    typedef gsSparseMatrix<weightType,gsEigen::RowMajor,indexType> LToGMatrix;
    typedef gsSparseMatrix<weightType,gsEigen::ColMajor,indexType> GToLMatrix;
 
    typedef std::vector<indexType> IndexContainer;
    typedef std::vector<indexType>::const_iterator CIndexIter;
    typedef std::vector<weightType> WeightContainer;
    typedef typename std::vector<weightType>::const_iterator CWeightIter;
    typedef typename LToGMatrix::InnerIterator InIterMat;
public:
    struct Iterator
    {
    friend class gsWeightMapper;
    private:
        const weightType  *m_beg;
        const weightType  *m_end;
        const weightType  *m_value;
        const indexType   *m_index;
        Iterator(const weightType  *beg,const weightType  *end,const indexType   *ind,const weightType  *value=NULL)
            : m_beg(beg), m_end(end), m_value(value?value:beg), m_index(ind)
        {}
        template<int ordering>
        Iterator( const gsSparseMatrix<weightType,ordering,indexType> &matrix, indexType outerId)
        {
            m_beg=matrix.valuePtr()+matrix.outerIndexPtr()[outerId];
            m_end=matrix.valuePtr()+matrix.outerIndexPtr()[outerId+1];
            m_value=m_beg;
            m_index=matrix.innerIndexPtr()+matrix.outerIndexPtr()[outerId];
        }
    public:
        Iterator()
            : m_beg(NULL),m_end(NULL),m_value(NULL),m_index(NULL)
            {
            }
 
        typedef std::bidirectional_iterator_tag iterator_category; //or another tag
 
        operator const weightType*() {return m_value;}
        operator const indexType*()  {return m_index;}
 
        // assignement shoud be implicitly generated by the compiler // Iterator& operator=(const Iterator& other) {m_value=other.m_value; m_index=other.m_index;  return *this;}
        bool operator==(const Iterator& other) const     {return m_value==other.m_value && m_index==other.m_index;}
        bool operator!=(const Iterator& other) const     {return !this->operator==(other);}
 
        Iterator& operator++()                     {++m_value;++m_index; return *this;}
        Iterator& operator--()                     {--m_value;--m_index; return *this;}
        Iterator& operator+=(ptrdiff_t a)          {m_value+=a; m_index+=a; return *this;}
        Iterator  operator+ (ptrdiff_t a) const    {return Iterator(m_beg,m_end, m_index+a, m_value+a);}
        Iterator& operator-=(ptrdiff_t a)          {m_value-=a; m_index-=a; return *this;}
        Iterator  operator- (ptrdiff_t a) const    {return Iterator(m_beg,m_end, m_index-a,m_value-a);}
 
        const weightType& operator*()    const     {return *m_value;}
 
        const weightType& weight()       const     {return *m_value;}
        const indexType&  index()        const     {return *m_index;}
        Iterator   end()        const     {return Iterator(m_beg,m_end,m_index,m_end);}
        Iterator  begin()      const     {return Iterator(m_beg,m_end,m_index,m_beg);}
 
        operator bool()        {return m_value<m_end && m_value>=m_beg;}
 
    };
 
    // operators and constructors
public:
    gsWeightMapper()
    { m_optimizationMatrix=NULL; }
 
    gsWeightMapper(indexType sourceSize,indexType targetSize)
    {
        m_matrix.resize(sourceSize,targetSize);
        m_optimizationMatrix=NULL;
    }
 
    virtual ~gsWeightMapper()
    {
        if(m_optimizationMatrix)
            delete m_optimizationMatrix;
    }
 
    gsWeightMapper(const gsWeightMapper& other)
        : m_matrix(other.m_matrix)
    {
        m_optimizationMatrix=NULL;
        optimize(other.getOptimizationFlags());
    }
 
    gsWeightMapper(const gsSparseMatrix<weightType,gsEigen::RowMajor,indexType> & other)
    {
        m_optimizationMatrix=NULL;
        *this=other;
    }
 
    gsWeightMapper(const gsSparseMatrix<weightType,gsEigen::ColMajor,indexType> & other)
    {
        m_optimizationMatrix=NULL;
        *this=other;
    }
 
    gsWeightMapper& operator=(const gsWeightMapper & other)
    {
        removeOptimization();
        m_matrix=other.m_matrix;
        optimize(other.getOptimizationFlags());
        return *this;
    }
 
    template<typename MatrixT>
    gsWeightMapper& operator*=(const MatrixT &other)
    {
        removeOptimization();
        // EIGEN problem
        // force a temporary because otherwise
        // permutation matrices do not work
        m_matrix=(m_matrix*other).eval();
        return *this;
    }
 
    template<typename MatrixT>
    gsWeightMapper operator*(const MatrixT &other)
    {
        gsWeightMapper result;
        result.m_matrix=m_matrix*other;
        return result;
    }
 
    template<typename MatrixT>
    gsWeightMapper& operator=(const MatrixT &other)
    {
        removeOptimization();
        m_matrix=other;
        optimize();
        return *this;
    }
 
    operator const LToGMatrix& ()
    { return m_matrix; }
 
    const LToGMatrix& asMatrix() const
    { return m_matrix; }
 
    LToGMatrix& asMatrix()
    { return m_matrix; }
    // functions for working with the mapper
public:
 
    void setEntry(indexType source, indexType target, weightType weight=1)
    {
        removeOptimization();
        m_matrix.at(source,target)=weight;
    }
 
    weightType getWeight(indexType source, indexType target) const
    { return m_matrix.at(source,target); }
 
    indexType getNrOfSources() const
    { return m_matrix.rows(); }
 
    indexType getNrOfTargets() const
    { return m_matrix.cols(); }
 
    bool sourceIsId(indexType source) const
    {
        IndexContainer indices;
        sourceToTarget(source,indices);
        return (indices.size()==1 && math::almostEqual<14>(m_matrix.at(source,indices[0]),T(1.0)));
    }
 
    bool targetIsId(indexType target) const
    {
        IndexContainer indices;
        targetToSource(target,indices);
        return (indices.size()==1 && math::almostEqual<14>(m_matrix.at(indices[0],target),T(1.0)));
    }
 
    // functions for transforming the coefficients
public:
 
    void mapToSourceCoefs(gsMatrix<weightType> const & targetCoefs,gsMatrix<weightType> & sourceCoefs) const
    {
        // from target to source it's just a multiplication
        GISMO_ASSERT(targetCoefs.rows()==m_matrix.cols(),"Wrong coefficient size!");
        sourceCoefs.noalias()=m_matrix * targetCoefs;
    }
 
    void mapToTargetCoefs(gsMatrix<weightType> const & sourceCoefs,gsMatrix<weightType> & targetCoefs) const;
 
    // functions for applying the map between target and source
public:
 
    void sourceToTarget(IndexContainer const & source,IndexContainer & target) const
    {
        target.reserve( source.size() );
        target.clear();
        IndexContainer res_i;
        for(size_t i = 0;i<source.size();i++)
        {
            //just call sourceToTarget for single bfs for all given sources
            sourceToTarget(source[i],res_i);
            target.insert(target.end(), res_i.begin(), res_i.end());
        }
        //get rid of duplicates and sort them
        sort( target.begin(), target.end() );
        target.erase( unique( target.begin(), target.end() ), target.end() );
    }
 
    void sourceToTarget(indexType source,IndexContainer & target) const
    {
        WeightContainer weights;
        sourceToTarget(source,target,weights);
    }
 
    void sourceToTarget(indexType source, IndexContainer & target, WeightContainer & weights) const;
 
    void targetToSource(IndexContainer const & target, IndexContainer & source) const
    {
        source.clear();
        IndexContainer res_i;
        for(size_t i = 0;i<target.size();i++)
        {
            //just call targetToSource for single bfs for all given targets
            targetToSource(target[i],res_i);
            source.insert(source.end(), res_i.begin(), res_i.end());
        }
        //get rid of duplicates and sort them
        sort( source.begin(), source.end() );
        source.erase( unique( source.begin(), source.end() ), source.end() );
    }
 
    void targetToSource(indexType target, IndexContainer & source) const
    {
        WeightContainer weights;
        targetToSource(target,source,weights);
    }
 
    void targetToSource(indexType target, IndexContainer & source, WeightContainer & weights) const;
 
    // functions for fast access to the mapping data
public:
    enum optimizeFlags
    {
        optSourceToTarget = 1U<<0,
        optTargetToSource = 1U<<1
    };
 
    void optimize (size_t flags=optSourceToTarget) const
    {
        flags|= flags & optTargetToSource ? optSourceToTarget : 0;
        flags = flags & ~getOptimizationFlags();
        if(flags & optSourceToTarget)
        {
            m_matrix.prune(1,T(10)*std::numeric_limits<weightType>::epsilon());
            m_matrix.makeCompressed();
        }
        if(flags & optTargetToSource)
        {
            this->m_optimizationMatrix=new GToLMatrix();
            *m_optimizationMatrix=m_matrix;
            m_optimizationMatrix->prune(1,T(10)*std::numeric_limits<weightType>::epsilon());
            m_optimizationMatrix->makeCompressed();
        }
    }
 
    size_t getOptimizationFlags() const
    {
        size_t flags=0;
        flags |= m_matrix.isCompressed() ? optSourceToTarget : 0;
        flags |= m_optimizationMatrix ? optTargetToSource : 0;
        return flags;
    }
 
    Iterator fastTargetToSource(indexType target) const
    {
        GISMO_ASSERT(m_optimizationMatrix,"optimize() must be called on the mapper with fastTargetToSource flag before using this function.");
        GISMO_ASSERT(target<m_matrix.cols() && 0<=target,"index out of bounds");
        return Iterator(*m_optimizationMatrix,target);
    }
    Iterator fastSourceToTarget(indexType source) const
    {
        GISMO_ASSERT(m_matrix.isCompressed(),"optimize() must be called on the mapper with fastSourceToTarget flag before using this function.");
        GISMO_ASSERT(source<m_matrix.rows() && 0<=source,"index "<<source<<" out of bounds (matrix.rows() = "<<m_matrix.rows());
        return Iterator(m_matrix,source);
    }
 
 
    void fastSourceToTarget(IndexContainer const & source,IndexContainer & target) const;
 
    void fastTargetToSource(IndexContainer const & target,IndexContainer & source) const;
 
 
     void getLocalMap (IndexContainer const & source, IndexContainer const & target, gsMatrix<T> &map) const;
 
     void getLocalMap (IndexContainer const & source, gsMatrix<T> &map) const;
 
 
private:
    void removeOptimization() {if(m_optimizationMatrix) delete m_optimizationMatrix; m_optimizationMatrix=NULL;}
    mutable GToLMatrix *m_optimizationMatrix;
    mutable LToGMatrix  m_matrix;
};//gsWeightMapper
 
}//namespace gismo
 
#ifndef GISMO_BUILD_LIB
#include GISMO_HPP_HEADER(gsWeightMapper.hpp)
#endif