/* * Bayes++ the Bayesian Filtering Library * Copyright (c) 2002 Michael Stevens * See accompanying Bayes++.htm for terms and conditions of use. * * $Id: uBLASmatrix.hpp 562 2006-04-05 20:46:23 +0200 (Wed, 05 Apr 2006) mistevens $ *//* * 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> funtionsusing 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;#elsenamespace 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_expressiontemplate <class E>BOOST_UBLAS_INLINEdetail::noalias_proxy<E> noalias (ublas::matrix_expression<E>& lvalue) {    return detail::noalias_proxy<E> (lvalue() );}template <class E>BOOST_UBLAS_INLINEdetail::noalias_proxy<E> noalias (ublas::vector_expression<E>& lvalue) {    return detail::noalias_proxy<E> (lvalue() );}#endifnamespace 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		VecBase::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		MatrixBase::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 matchingusing detail::FMMatrix;							// Default typestypedef 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 typestypedef 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)#ifdef BAYES_FILTER_GAPPY