ublasmatrix.hpp

良好的代码实现

/*
 * Bayes++ the Bayesian Filtering Library
 * Copyright (c) 2002 Michael Stevens
 * See accompanying Bayes++.htm for terms and conditions of use.
 *
 * $Header: /cvsroot/bayesclasses/Bayes++/BayesFilter/uBLASmatrix.hpp,v 1.24.2.10 2005/07/16 07:29:27 mistevens Exp $
 * $NoKeywords: $
 */

/*
 * Common type independant uBlas interface
 *  Should be include after base types have been defined
 *
 * Everything in namespace Bayes_filter_matrix is intended to support the matrix storage
 * and algebra requirements of the library. Therefore the interfaces and implementation is
 * not intended to be stable. Nor is this a general purpose adapator for uBLAS
 *
 * Note on row_major matrices
 *  The implementation uses row major extensively. The computation of symetric products P*P' is 
 *  most efficient with row operations. These products are used extensively so the default
 *  is to use row_major matrices
 */

/* Filter Matrix Namespace */
namespace Bayesian_filter_matrix
{
						// Allow use a few functions in own namespace (particular useful for compilers with Konig lookup)
using ublas::row;
using ublas::column;
using ublas::trans;
using ublas::prod;		// These do not apply to the templated prod<temp> funtions
using ublas::inner_prod;
using ublas::outer_prod;

enum EmptyTag {Empty};	// Tag type used for empty matrix constructor

						// Old compiler workaround removed from new uBLAS
#ifndef BOOST_UBLAS_TYPENAME
#define BOOST_UBLAS_TYPENAME typename
#endif

#if (BOOST_VERSION >= 103100)
using ublas::noalias;
#else
namespace detail
{
	// Assignment proxy.
    // Provides temporary free assigment when LHS has no alias on RHS
    template<class C>
    class noalias_proxy {
    public:
        BOOST_UBLAS_INLINE
        noalias_proxy (C& lval):
            lval_ (lval) {}

        template <class E>
        BOOST_UBLAS_INLINE
        void operator= (const E& e) {
            lval_.assign (e);
        }

        template <class E>
        BOOST_UBLAS_INLINE
        void operator+= (const E& e) {
            lval_.plus_assign (e);
        }

        template <class E>
        BOOST_UBLAS_INLINE
        void operator-= (const E& e) {
            lval_.minus_assign (e);
        }

    private:  // nonassignable
        void operator=( const noalias_proxy& );
    private:
        C& lval_;
    };

}//namespace detail

// Improve syntax of effcient assignment where no aliases of LHS appear on the RHS
//  noalias(lhs) = rhs_expression
template <class E>
BOOST_UBLAS_INLINE
detail::noalias_proxy<E> noalias (ublas::matrix_expression<E>& lvalue) {
    return detail::noalias_proxy<E> (lvalue() );
}
template <class E>
BOOST_UBLAS_INLINE
detail::noalias_proxy<E> noalias (ublas::vector_expression<E>& lvalue) {
    return detail::noalias_proxy<E> (lvalue() );
}

#endif


namespace detail		// Lots of implementation detail
{

/*
 * Filter Vec type
 */
template <class VecBase>
class FMVec : public VecBase
{
public:
	typedef typename VecBase::value_type value_type;
	typedef typename VecBase::vector_temporary_type vector_temporary_type;

	// No Default Constructor. Empty creation is very error prone
	explicit FMVec(EmptyTag) : VecBase()
	{}	// Empty constructor
	explicit FMVec(std::size_t size) : VecBase(size)
	{}	// Normal sized constructor
	FMVec(const FMVec& c) : VecBase(static_cast<const VecBase&>(c))
	{}	// Copy constructor
	template <class E>
	explicit FMVec(const ublas::vector_expression<E>& e) : VecBase(e)
	{}	// vector_expression copy constructor
	template <class E>
	FMVec(const ublas::matrix_column<E>& e) : VecBase(e)
	{}	// conversion copy constructor, hides the implict copy required for matrix column access

	template <class E>
	FMVec& operator= (const ublas::vector_expression<E>& r)
	{	// Expression assignment; may be dependant on r
		VecBase::operator=(r);
		return *this;
	}
	FMVec& operator= (const FMVec& r)
	{	// Vector assignment; independant
		assign(r);
		return *this;
	}

	// Sub-range selection operators
	const ublas::vector_range<const VecBase> sub_range(std::size_t b, std::size_t e) const
	{
		return ublas::vector_range<const VecBase>(*this, ublas::range(b,e));
	}
	ublas::vector_range<VecBase> sub_range(std::size_t b, std::size_t e)
	{
		return ublas::vector_range<VecBase>(*this, ublas::range(b,e));
	}
};


/*
 * Filter Matrix class template. Augmentation for uBlas MatrixBase
 */
template <class MatrixBase>
class FMMatrix : public MatrixBase
{
public:
	typedef typename MatrixBase::value_type value_type;
	typedef typename MatrixBase::vector_temporary_type vector_temporary_type;
	typedef typename MatrixBase::matrix_temporary_type matrix_temporary_type;

	// No Default Constructor. Empty creation is very error prone
	explicit FMMatrix(EmptyTag) : MatrixBase()
	{}	// Empty constructor
	FMMatrix(std::size_t size1, std::size_t size2) : MatrixBase(size1,size2)
	{}	// Normal sized constructor
	FMMatrix(const FMMatrix& c) : MatrixBase(static_cast<const MatrixBase&>(c))
	{}	// Copy constructor
	template <class E>
	explicit FMMatrix(const ublas::matrix_expression<E>& e) : MatrixBase(e)
	{}	// matrix_expression copy constructor

	template <class E>
	FMMatrix& operator= (const ublas::matrix_expression<E>& r)
	{	// Expression assignment; may be dependant on r
		MatrixBase::operator=(r);
		return *this;
	}
	FMMatrix& operator= (const FMMatrix& r)
	{	// Matrix assignment; independant
		assign (r);
		return *this;
	}

	// Row,Column vector proxies
	typedef ublas::matrix_row<FMMatrix> Row;
	typedef const ublas::matrix_row<const FMMatrix> const_Row;
	typedef ublas::matrix_column<FMMatrix> Column;
	typedef const ublas::matrix_column<const FMMatrix> const_Column;

	// Vector proxies from iterators - static members dependant on MatrixBase type
	// ri() returns container associated with iterator. static_cast required as typeof(ri()) may not be typeof(MM)
	static Row rowi(const typename MatrixBase::iterator1& ri)
	{
		return Row(static_cast<FMMatrix&>(ri()), ri.index1());
	}
	static const_Row rowi(const typename MatrixBase::const_iterator1& ri)
	{
		return const_Row(static_cast<const FMMatrix&>(ri()), ri.index1());
	}
	static Column columni(const typename MatrixBase::iterator2& ci)
	{
		return Column(static_cast<FMMatrix&>(ci()), ci.index2());
	}
	static const_Column columni(const typename MatrixBase::const_iterator2& ci)
	{
		return const_Column(static_cast<const FMMatrix&>(ci()), ci.index2());
	}

	// Sub-range selection operators
	ublas::matrix_range<const MatrixBase>
	sub_matrix(std::size_t s1, std::size_t e1, std::size_t s2, std::size_t e2) const
	{
		return ublas::matrix_range<const MatrixBase> (*this, ublas::range(s1,e1), ublas::range(s2,e2));
	}
	ublas::matrix_range<MatrixBase>
	sub_matrix(std::size_t s1, std::size_t e1, std::size_t s2, std::size_t e2)
	{
		return ublas::matrix_range<MatrixBase> (*this, ublas::range(s1,e1), ublas::range(s2,e2));
	}

	// Requires boost_1.30.0 which has a generalised matrix_vector_slice
	ublas::matrix_vector_slice<const MatrixBase>
	sub_column(std::size_t s1, std::size_t e1, std::size_t s2) const 
	// Column vector s2 with rows [s1,e1)
	{
		return ublas::matrix_vector_slice<const MatrixBase> (*this, ublas::slice(s1,1,e1-s1), ublas::slice(s2,0,e1-s1));
	}
	ublas::matrix_vector_slice<MatrixBase>
	sub_column(std::size_t s1, std::size_t e1, std::size_t s2)
	// Column vector s2 with rows [s1,e1)
	{
		return ublas::matrix_vector_slice<MatrixBase> (*this, ublas::slice(s1,1,e1-s1), ublas::slice(s2,0,e1-s1));
	}
};


/*
 * Helper template to allow member construction before base class
 *  Boost version does not work as it passes by value
 */
template <typename MemberType>
class BaseFromMember
{
protected:
	MemberType member;
	explicit BaseFromMember() : member()
	{}

	template <typename T1>
	explicit BaseFromMember( const T1& x1 ) : member( x1 )
	{}

	template <typename T1, typename T2>
	explicit BaseFromMember( const T1& x1, const T2& x2 ) : member( x1, x2 )
	{}
};


/*
 * We require static type conversion between Symmetric matrices and equivilent row major matrices
 * Therefore we create symmetric matrix types, using a MatrixBase for storage
 * and wraps this in a symmetric_adaptor
 */
template <class MatrixBase>
class SymMatrixWrapper :
	private BaseFromMember<MatrixBase>,  // allow construction of MatrixBase member before symmetric_adaptor
	public ublas::symmetric_adaptor<MatrixBase, ublas::upper>
{
	typedef BaseFromMember<MatrixBase> matrix_type;
	typedef ublas::symmetric_adaptor<MatrixBase, ublas::upper> symadaptor_type;
public:
	typedef typename MatrixBase::value_type value_type;
	typedef typename MatrixBase::vector_temporary_type vector_temporary_type;
	typedef typename MatrixBase::matrix_temporary_type matrix_temporary_type;

	SymMatrixWrapper () : matrix_type(), symadaptor_type(matrix_type::member)
	{}
	SymMatrixWrapper (std::size_t size1, std::size_t size2) : matrix_type(size1,size2), symadaptor_type(matrix_type::member)
	{}	// Normal sized constructor
	explicit SymMatrixWrapper (const SymMatrixWrapper& r) : matrix_type(reinterpret_cast<const MatrixBase&>(r)), symadaptor_type(matrix_type::member)
	{}	// Explict copy construction referencing the copy reinterpreted as a MatrixBase
	template <class E>
	explicit SymMatrixWrapper (const ublas::matrix_expression<E>& e) : matrix_type(e), symadaptor_type(matrix_type::member)
	{}	// Explict matrix_expression conversion constructor

	template <class E>
	SymMatrixWrapper& operator=(const ublas::matrix_expression<E>& r)
	{
		symadaptor_type::operator=(r);
		return *this;
	}

	// Conversions straight to a FMMatrix, equivilent to a RowMatrix types
	const FMMatrix<MatrixBase>& asRowMatrix() const
	{
		return static_cast<const FMMatrix<MatrixBase>& >(matrix_type::member);
	}
	FMMatrix<MatrixBase>& asRowMatrix()
	{
		return static_cast<FMMatrix<MatrixBase>& >(matrix_type::member);
	}

	// Matrix storage members
	void clear()
	{	matrix_type::member.clear();
	}
	void resize(std::size_t nsize1, std::size_t nsize2, bool preserve = true)
	{
		matrix_type::member.resize(nsize1, nsize2, preserve);
	}
};

}//namespace detail



/*
 * Vector / Matrix types
 *  Finally the definitions !
 */
using detail::FMVec;		// Template class for template parameter matching
using detail::FMMatrix;

							// Default types
typedef FMVec<detail::BaseVector> Vec;
typedef FMMatrix<detail::BaseRowMatrix> RowMatrix;
typedef RowMatrix Matrix;
typedef FMMatrix<detail::BaseColMatrix> ColMatrix;
typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseRowMatrix> > SymMatrix;
typedef FMMatrix<detail::BaseUpperTriMatrix> UTriMatrix;
typedef FMMatrix<detail::BaseLowerTriMatrix> LTriMatrix;
typedef FMMatrix<detail::BaseDiagMatrix> DiagMatrix;

							// Explicitly dense types
typedef FMVec<detail::BaseDenseVector> DenseVec;
typedef FMMatrix<detail::BaseDenseRowMatrix> DenseRowMatrix;
typedef DenseRowMatrix DenseMatrix;
typedef FMMatrix<detail::BaseDenseColMatrix> DenseColMatrix;
typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseDenseRowMatrix> > DenseSymMatrix;
typedef FMMatrix<detail::BaseDenseUpperTriMatrix> DenseUTriMatrix;
typedef FMMatrix<detail::BaseDenseLowerTriMatrix> DenseLTriMatrix;
typedef FMMatrix<detail::BaseDenseDiagMatrix> DenseDiagMatrix;

							// Explicitly sparse types (any of the gappy types)