ublasmatrix.hpp

良好的代码实现

#ifdef BAYES_FILTER_GAPPY
typedef FMVec<detail::BaseSparseVector> SparseVec;
typedef FMMatrix<detail::BaseDenseRowMatrix> SparseRowMatrix;
typedef SparseRowMatrix SparseMatrix;
typedef FMMatrix<detail::BaseSparseColMatrix> SparseColMatrix;
typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseSparseRowMatrix> > SparseSymMatrix;
#endif


/*
 * Matrix Adaptors, simply hide the uBLAS details
 */
template <class M>
const ublas::triangular_adaptor<const M, ublas::upper>
 UpperTri(const M& m)
/*
 * View Upper triangle of m
 * ISSUE VC7 cannot cope with UTriMatrix::functor1_type
 */
{
	return ublas::triangular_adaptor<const M, ublas::upper>(m);
}

template <class M>
const ublas::triangular_adaptor<const M, ublas::lower>
 LowerTri(const M& m)
/*
 * View Lower triangle of m
 */
{
	return ublas::triangular_adaptor<const M, ublas::lower>(m);
}



/*
 * Matrix Support Operations
 */
template <class Base>
ublas::matrix_vector_range<FMMatrix<Base> >
 diag(FMMatrix<Base>& M, std::size_t n)
{	// Return a vector proxy to the first n diagonal elements of M
	return ublas::matrix_vector_range<FMMatrix<Base> >(M, ublas::range(0,n), ublas::range(0,n));
}

template <class Base>
const ublas::matrix_vector_range<const FMMatrix<Base> >
 diag(const FMMatrix<Base>& M, std::size_t n)
{	// Return a vector proxy to the first n diagonal elements of M
	return ublas::matrix_vector_range<const FMMatrix<Base> >(M, ublas::range(0,n), ublas::range(0,n));
}

template <class Base>
ublas::matrix_vector_range<FMMatrix<Base> >
 diag(FMMatrix<Base>& M)
{	// Return a vector proxy to the diagonal elements of M
	const std::size_t common_size = std::min(M.size1(),M.size2());
	return ublas::matrix_vector_range<FMMatrix<Base> >(M, ublas::range(0,common_size), ublas::range(0,common_size));
}

template <class Base>
const ublas::matrix_vector_range<const FMMatrix<Base> >
 diag(const FMMatrix<Base>& M)
{	// Return a vector proxy to the diagonal elements of M
	const std::size_t common_size = std::min(M.size1(),M.size2());
	return ublas::matrix_vector_range<const FMMatrix<Base> >(M, ublas::range(0,common_size), ublas::range(0,common_size));
}

template <class Base>
void identity(FMMatrix<Base>& I)
{	// Set I to generalised Identity matrix. Clear I and set diag(I) to one
	I.clear();
							// Set common diagonal elements
	std::size_t common_size = std::min(I.size1(),I.size2());
	typedef typename Base::value_type Base_value_type;
	diag(I).assign (ublas::scalar_vector<Base_value_type>(common_size, Base_value_type(1)));
}



/*
 * Symmetric Positive (Semi) Definate multiplication: X*S*X' and X'*S*X
 *  The name is slightly misleading. The result is actually only PD if S is PD
 *  Algorithms are intended to exploit the symmerty of the result
 *  and also where possible the row by row multiplication inherent in X*X'
 */

namespace detail {		// mult_SPD now an implementation detail
	
template <class MatrixX>
void mult_SPD (const MatrixX& X, const Vec& s, SymMatrix& P)
/*
 * Symmetric Positive (Semi) Definate multiply:	P += X*diag_matrix(s)*X'
 */
{
	Vec::const_iterator si, send = s.end();
	typename MatrixX::const_iterator1 Xa = X.begin1();
	const typename MatrixX::const_iterator1 Xend = X.end1();
	typename MatrixX::const_iterator1 Xb;

	// P(a,b) = X.row(a) * X.row(b)
	for (; Xa != Xend; ++Xa)				// Iterate Rows
	{
		typename MatrixX::const_Row Xav = MatrixX::rowi(Xa);
		Xb = Xa;							// Start at the row Xa only one triangle of symetric result required
		for (; Xb != Xend; ++Xb)
		{
			SymMatrix::value_type p = 0;	// Tripple vector inner product
			typename MatrixX::const_Row Xbv = MatrixX::rowi(Xb);
			for (si = s.begin(); si != send; ++si) {
				Vec::size_type i = si.index();
				p += Xav[i] * (*si) * Xbv[i];
			}
			P(Xa.index1(),Xb.index1()) += p;
		}
	}
}

template <class MatrixX>
void mult_SPDT (const MatrixX& X, const Vec& s, SymMatrix& P)
/*
 * Symmetric Positive (Semi) Definate multiply:	P += X'*diag_matrix(s)*X
 */
{
	Vec::const_iterator si, send = s.end();
	typename MatrixX::const_iterator2 Xa = X.begin2();
	const typename MatrixX::const_iterator2 Xend = X.end2();
	typename MatrixX::const_iterator2 Xb;

	// P(a,b) = X.col(a) * X.col(b)
	for (; Xa != Xend; ++Xa)				// Iterate vectors
	{
		typename MatrixX::const_Column Xav = MatrixX::columni(Xa);
		Xb = Xa;							// Start at the row Xa only one triangle of symertric result required
		for (; Xb != Xend; ++Xb)
		{									// Tripple vector inner product
			SymMatrix::value_type p = 0;
			typename MatrixX::const_Column Xbv = MatrixX::columni(Xb);
			for (si = s.begin(); si != send; ++si) {
				Vec::size_type i = si.index();
				p += Xav[i] * (*si) * Xbv[i];
			}
			P(Xa.index2(),Xb.index2()) += p;
		}
	}
}

}//namespace detail


/*
 * prod_SPD - uBLAS Expression templates for X*X' and X*X'
 * Functions are only defined for the type for which the operation is efficient
 * ISSUE Although numerically symmetric, uBlas has no expression type to represent this property
 *  The result must be assigned to a symmetric container to exploit the symmetry
 */

template <class E1, class E2>
struct prod_expression_result
{	// Provide ET result E1E2T_type of prod(matrix_expression<E1>,trans(matrix_expression<E2>)
	typedef BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits<E2, ublas::scalar_identity<BOOST_UBLAS_TYPENAME E2::value_type> >::result_type  E2T_type;
	typedef BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits<BOOST_UBLAS_TYPENAME E1::value_type, E1,
                                        BOOST_UBLAS_TYPENAME E2T_type::value_type, E2T_type>::result_type  E1E2T_type;

	// Provide ET result E1TE2_type of prod(trans(matrix_expression<E1>),matrix_expression<E2>)
	typedef BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits<E1, ublas::scalar_identity<BOOST_UBLAS_TYPENAME E1::value_type> >::result_type  E1T_type;
	typedef BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits<BOOST_UBLAS_TYPENAME E1T_type::value_type, E1T_type,
                                        BOOST_UBLAS_TYPENAME E2::value_type, E2>::result_type  E1TE2_type;
};

 
template<class E> inline
typename prod_expression_result<E,E>::E1E2T_type
 prod_SPD (const ublas::matrix_expression<E>& X)
/*
 * Symmetric Positive (Semi) Definate product: X*X'
 */
{
	// ISSUE: uBLAS post Boost 1_31_0 introduces a trans(const matrix_expression<E>& e) which propogates the const expression type
	// Bypass this to avoid having to detect the Boost version
	return prod( X, trans(const_cast<ublas::matrix_expression<E>&>(X) ));
}

template<class EX, class ES, class ET> inline
typename prod_expression_result<EX,ET>::E1E2T_type
 prod_SPD (const ublas::matrix_expression<EX>& X, const ublas::matrix_expression<ES>& S, ublas::matrix_expression<ET>& XStemp)
/*
 * Symmetric Positive (Semi) Definate product: X*(X*S)', XStemp = X*S
 *  Result symmetric if S is symmetric
 */
{
	return prod( X, trans(prod(X,S,XStemp())) );
}

template<class EX, class ES, class ET> inline
ET prod_SPD (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s, ublas::matrix_expression<ET>& Ptemp)
/*
 * Symmetric Positive (Semi) Definate product: X*diag_matrix(s)*X'
 * Precond: Ptemp must be size conformant with the product
 *  DEPRECATED With explicit Ptemp, do not rely on it's value
 */
{
	Vec::const_iterator si, send = s().end();
	const EX& XX = X();
	typename EX::const_iterator1 Xa = XX.begin1();
	const typename EX::const_iterator1 Xend = XX.end1();
	typename EX::const_iterator1 Xb;

	// P(a,b) = sum(X.row(a) * s * X.row(b))
	for (; Xa != Xend; ++Xa)		// Iterate rows
	{
		typedef const ublas::matrix_row<const EX> EX_row;
		EX_row Xav (ublas::row(XX, Xa.index1()));
		Xb = Xa;							// Start at the row Xa only one triangle of symetric result required
		for (; Xb != Xend; ++Xb)
		{
			typename EX::value_type p = 0;	// Triple vector inner product
			EX_row Xbv (ublas::row(XX, Xb.index1()));
			for (si = s().begin(); si != send; ++si) {
				Vec::size_type i = si.index();
				p += Xav[i] * (*si) * Xbv[i];
			}
			Ptemp()(Xa.index1(),Xb.index1()) = p;
			Ptemp()(Xb.index1(),Xa.index1()) = p;
		}
	}
	return Ptemp();
}

template<class EX, class ES> inline
SymMatrix prod_SPD (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s)
/*
 * Symmetric Positive (Semi) Definate product: X*diag_matrix(s)*X'
 * TODO requires a prod_diag(X,s)
 */
{
	SymMatrix Ptemp(X().size1(),X().size1());
	Ptemp.clear();

	(void) prod_SPD (X(),s(), Ptemp);
	return Ptemp;
}


template<class E> inline
typename prod_expression_result<E,E>::E1TE2_type
 prod_SPDT (const ublas::matrix_expression<E>& X)
/*
 * Symmetric Positive (Semi) Definate product: X'*X
 */
{
	// ISSUE: See prod_SPD
	return prod( trans(const_cast<ublas::matrix_expression<E>&>(X) ), X);
}

template<class EX, class ES, class ET> inline
typename prod_expression_result<ET,EX>::E1TE2_type
 prod_SPDT (const ublas::matrix_expression<EX>& X, const ublas::matrix_expression<ES>& S, ublas::matrix_expression<ET>& SXtemp)
/*
 * Symmetric Positive (Semi) Definate product: (S*X)'*X, SXtemp = S*X
 *  Result symmetric if S is symmetric
 */
{
	return prod( trans(prod(S,X,SXtemp())), X);
}

template<class EX, class ES, class ET> inline
ET prod_SPDT (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s, ublas::matrix_expression<ET>& Ptemp)
/*
 * Symmetric Positive (Semi) Definate product: X'*diag_matrix(s)*X, Ptemp = return value
 * Precond: Ptemp must be size conformant with the product
 *  DEPRECATED With explicit Ptemp, do not rely on it's value
 */
{
	const EX& XX = X();
	Vec::const_iterator si, send = s().end();
	typename EX::const_iterator2 Xa = X.begin2();
	const typename EX::const_iterator2 Xend = X.end2();
	typename EX::const_iterator2 Xb;

	// P(a,b) = sum(X.col(a) * s * X.col(b))
	for (; Xa != Xend; ++Xa)		// Iterate columns
	{
		typedef const ublas::matrix_column<const EX> EX_column;
		EX_column Xav = ublas::column(XX, Xa.index2());
		Xb = Xa;							// Start at the column Xa only one triangle of symetric result required
		for (; Xb != Xend; ++Xb)
		{
			typename EX::value_type p = 0;	// Triple vector inner product
			EX_column Xbv = ublas::column(XX, Xb.index2());
			for (si = s().begin(); si != send; ++si) {
				Vec::size_type i = si.index();
				p += Xav[i] * (*si) * Xbv[i];
			}
			Ptemp()(Xa.index2(),Xb.index2()) = p;
			Ptemp()(Xb.index2(),Xa.index2()) = p;
		}
	}
	return Ptemp();
}

template<class EX, class ES> inline
SymMatrix prod_SPDT (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s)
/*
 * Symmetric Positive (Semi) Definate product: X'*diag_matrix(s)*X
 * TODO requires a prod_diag(X,s)
 */
{
	SymMatrix Ptemp(X().size1(),X().size1());
	Ptemp.clear();

	(void) prod_SPDT (X(),s(), Ptemp);
	return Ptemp;
}

inline Vec::value_type prod_SPDT (const Vec& x, const Vec& s)
/*
 * Symmetric Positive (Semi) Definate product: x'*diag_matrix(s)*x
 */
{
	const Vec::const_iterator xi_end = x.end();
	Vec::const_iterator si = s.begin();
	Vec::const_iterator xi=x.begin();
	Vec::value_type p = Vec::value_type(0);
	while (xi != xi_end)
	{
		p += (*xi) * (*si) * (*xi);
		++xi; ++ si;
	}
	
	return p;
}

inline Vec::value_type prod_SPDT (const Vec& x, const SymMatrix& S)
/*
 * Symmetric Positive (Semi) Definate multiply:	p = x'*S*x
 */
{
	return inner_prod (x, prod(S,x) );
}

}//namespace