31using ublas::inner_prod;
32using ublas::outer_prod;
37#ifndef BOOST_UBLAS_TYPENAME
38#define BOOST_UBLAS_TYPENAME typename
41#if (BOOST_VERSION >= 103100)
64 lval_.plus_assign (e);
70 lval_.minus_assign (e);
103template <
class VecBase>
107 typedef typename VecBase::value_type value_type;
108 typedef typename VecBase::vector_temporary_type vector_temporary_type;
113 explicit FMVec(std::size_t size) :
VecBase(size)
118 explicit FMVec(
const ublas::vector_expression<E>& e) :
VecBase(e)
121 FMVec(
const ublas::matrix_column<E>& e) :
VecBase(e)
125 FMVec& operator= (
const ublas::vector_expression<E>& r)
127 VecBase::operator=(r);
130 FMVec& operator= (
const FMVec& r)
137 const ublas::vector_range<const VecBase> sub_range(std::size_t b, std::size_t e)
const
139 return ublas::vector_range<const VecBase>(*
this, ublas::range(b,e));
141 ublas::vector_range<VecBase> sub_range(std::size_t b, std::size_t e)
143 return ublas::vector_range<VecBase>(*
this, ublas::range(b,e));
151template <
class MatrixBase>
155 typedef typename MatrixBase::value_type value_type;
156 typedef typename MatrixBase::vector_temporary_type vector_temporary_type;
157 typedef typename MatrixBase::matrix_temporary_type matrix_temporary_type;
162 FMMatrix(std::size_t size1, std::size_t size2) :
MatrixBase(size1,size2)
167 explicit FMMatrix(
const ublas::matrix_expression<E>& e) :
MatrixBase(e)
171 FMMatrix& operator= (
const ublas::matrix_expression<E>& r)
173 MatrixBase::operator=(r);
176 FMMatrix& operator= (
const FMMatrix& r)
178 MatrixBase::assign (r);
183 typedef ublas::matrix_row<FMMatrix> Row;
184 typedef const ublas::matrix_row<const FMMatrix> const_Row;
185 typedef ublas::matrix_column<FMMatrix> Column;
186 typedef const ublas::matrix_column<const FMMatrix> const_Column;
190 static Row rowi(
const typename MatrixBase::iterator1& ri)
192 return Row(
static_cast<FMMatrix&
>(ri()), ri.index1());
194 static const_Row rowi(
const typename MatrixBase::const_iterator1& ri)
196 return const_Row(
static_cast<const FMMatrix&
>(ri()), ri.index1());
198 static Column columni(
const typename MatrixBase::iterator2& ci)
200 return Column(
static_cast<FMMatrix&
>(ci()), ci.index2());
202 static const_Column columni(
const typename MatrixBase::const_iterator2& ci)
204 return const_Column(
static_cast<const FMMatrix&
>(ci()), ci.index2());
208 ublas::matrix_range<const MatrixBase>
209 sub_matrix(std::size_t s1, std::size_t e1, std::size_t s2, std::size_t e2)
const
211 return ublas::matrix_range<const MatrixBase> (*
this, ublas::range(s1,e1), ublas::range(s2,e2));
213 ublas::matrix_range<MatrixBase>
214 sub_matrix(std::size_t s1, std::size_t e1, std::size_t s2, std::size_t e2)
216 return ublas::matrix_range<MatrixBase> (*
this, ublas::range(s1,e1), ublas::range(s2,e2));
220 ublas::matrix_vector_slice<const MatrixBase>
221 sub_column(std::size_t s1, std::size_t e1, std::size_t s2)
const
224 return ublas::matrix_vector_slice<const MatrixBase> (*
this, ublas::slice(s1,1,e1-s1), ublas::slice(s2,0,e1-s1));
226 ublas::matrix_vector_slice<MatrixBase>
227 sub_column(std::size_t s1, std::size_t e1, std::size_t s2)
230 return ublas::matrix_vector_slice<MatrixBase> (*
this, ublas::slice(s1,1,e1-s1), ublas::slice(s2,0,e1-s1));
239template <
typename MemberType>
244 explicit BaseFromMember() : member()
247 template <
typename T1>
248 explicit BaseFromMember(
const T1& x1 ) : member( x1 )
251 template <
typename T1,
typename T2>
252 explicit BaseFromMember(
const T1& x1,
const T2& x2 ) : member( x1, x2 )
262template <
class MatrixBase>
263class SymMatrixWrapper :
264 private BaseFromMember<MatrixBase>,
265 public ublas::symmetric_adaptor<MatrixBase, ublas::upper>
267 typedef BaseFromMember<MatrixBase> matrix_type;
268 typedef ublas::symmetric_adaptor<MatrixBase, ublas::upper> symadaptor_type;
270 typedef typename MatrixBase::value_type value_type;
271 typedef typename MatrixBase::vector_temporary_type vector_temporary_type;
272 typedef typename MatrixBase::matrix_temporary_type matrix_temporary_type;
274 SymMatrixWrapper () : matrix_type(), symadaptor_type(matrix_type::member)
276 SymMatrixWrapper (std::size_t size1, std::size_t size2) : matrix_type(size1,size2), symadaptor_type(matrix_type::member)
278 explicit SymMatrixWrapper (
const SymMatrixWrapper& r) : matrix_type(reinterpret_cast<const
MatrixBase&>(r)), symadaptor_type(matrix_type::member)
281 explicit SymMatrixWrapper (
const ublas::matrix_expression<E>& e) : matrix_type(e), symadaptor_type(matrix_type::member)
285 SymMatrixWrapper& operator=(
const ublas::matrix_expression<E>& r)
287 symadaptor_type::operator=(r);
292 const FMMatrix<MatrixBase>& asRowMatrix()
const
294 return static_cast<const FMMatrix<MatrixBase>&
>(matrix_type::member);
296 FMMatrix<MatrixBase>& asRowMatrix()
298 return static_cast<FMMatrix<MatrixBase>&
>(matrix_type::member);
303 { matrix_type::member.clear();
305 void resize(std::size_t nsize1, std::size_t nsize2,
bool preserve =
true)
307 matrix_type::member.resize(nsize1, nsize2, preserve);
320using detail::FMMatrix;
323typedef FMVec<detail::BaseVector>
Vec;
327typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseRowMatrix> >
SymMatrix;
337typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseDenseRowMatrix> >
DenseSymMatrix;
343#ifdef BAYES_FILTER_GAPPY
344typedef FMVec<detail::BaseSparseVector> SparseVec;
345typedef FMMatrix<detail::BaseDenseRowMatrix> SparseRowMatrix;
346typedef SparseRowMatrix SparseMatrix;
347typedef FMMatrix<detail::BaseSparseColMatrix> SparseColMatrix;
348typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseSparseRowMatrix> > SparseSymMatrix;
356const ublas::triangular_adaptor<const M, ublas::upper>
363 return ublas::triangular_adaptor<const M, ublas::upper>(m);
367const ublas::triangular_adaptor<const M, ublas::lower>
373 return ublas::triangular_adaptor<const M, ublas::lower>(m);
382ublas::matrix_vector_range<FMMatrix<Base> >
383 diag(FMMatrix<Base>& M, std::size_t n)
385 return ublas::matrix_vector_range<FMMatrix<Base> >(M, ublas::range(0,n), ublas::range(0,n));
389const ublas::matrix_vector_range<const FMMatrix<Base> >
390 diag(
const FMMatrix<Base>& M, std::size_t n)
392 return ublas::matrix_vector_range<const FMMatrix<Base> >(M, ublas::range(0,n), ublas::range(0,n));
396ublas::matrix_vector_range<FMMatrix<Base> >
399 const std::size_t common_size = std::min(M.size1(),M.size2());
400 return ublas::matrix_vector_range<FMMatrix<Base> >(M, ublas::range(0,common_size), ublas::range(0,common_size));
404const ublas::matrix_vector_range<const FMMatrix<Base> >
407 const std::size_t common_size = std::min(M.size1(),M.size2());
408 return ublas::matrix_vector_range<const FMMatrix<Base> >(M, ublas::range(0,common_size), ublas::range(0,common_size));
416 std::size_t common_size = std::min(I.size1(),I.size2());
417 typedef typename Base::value_type Base_value_type;
418 diag(I).assign (ublas::scalar_vector<Base_value_type>(common_size, Base_value_type(1)));
432template <
class MatrixX>
438 Vec::const_iterator
si,
send = s.end();
439 typename MatrixX::const_iterator1
Xa = X.begin1();
440 const typename MatrixX::const_iterator1
Xend = X.end1();
441 typename MatrixX::const_iterator1
Xb;
446 typename MatrixX::const_Row
Xav = MatrixX::rowi(
Xa);
450 SymMatrix::value_type
p = 0;
451 typename MatrixX::const_Row
Xbv = MatrixX::rowi(
Xb);
453 Vec::size_type i =
si.index();
456 P(
Xa.index1(),
Xb.index1()) +=
p;
461template <
class MatrixX>
467 Vec::const_iterator
si,
send = s.end();
468 typename MatrixX::const_iterator2
Xa = X.begin2();
469 const typename MatrixX::const_iterator2
Xend = X.end2();
470 typename MatrixX::const_iterator2
Xb;
475 typename MatrixX::const_Column
Xav = MatrixX::columni(
Xa);
479 SymMatrix::value_type
p = 0;
480 typename MatrixX::const_Column
Xbv = MatrixX::columni(
Xb);
482 Vec::size_type i =
si.index();
485 P(
Xa.index2(),
Xb.index2()) +=
p;
500template <
class E1,
class E2>
514template<
class E>
inline
523 return prod( X, trans(
const_cast<ublas::matrix_expression<E>&
>(X) ));
526template<
class EX,
class ES,
class ET>
inline
528 prod_SPD (
const ublas::matrix_expression<EX>& X,
const ublas::matrix_expression<ES>& S, ublas::matrix_expression<ET>& XStemp)
534 return prod( X, trans(prod(X,S,XStemp())) );
537template<
class EX,
class ES,
class ET>
inline
538ET
prod_SPD (
const ublas::matrix_expression<EX>& X,
const ublas::vector_expression<ES>& s, ublas::matrix_expression<ET>& Ptemp)
545 Vec::const_iterator si, send = s().end();
547 typename EX::const_iterator1 Xa = XX.begin1();
548 const typename EX::const_iterator1 Xend = XX.end1();
549 typename EX::const_iterator1 Xb;
552 for (; Xa != Xend; ++Xa)
554 typedef const ublas::matrix_row<const EX> EX_row;
555 EX_row Xav (ublas::row(XX, Xa.index1()));
557 for (; Xb != Xend; ++Xb)
559 typename EX::value_type p = 0;
560 EX_row Xbv (ublas::row(XX, Xb.index1()));
561 for (si = s().begin(); si != send; ++si) {
562 Vec::size_type i = si.index();
563 p += Xav[i] * (*si) * Xbv[i];
565 Ptemp()(Xa.index1(),Xb.index1()) = p;
566 Ptemp()(Xb.index1(),Xa.index1()) = p;
572template<
class EX,
class ES>
inline
573SymMatrix prod_SPD (
const ublas::matrix_expression<EX>& X,
const ublas::vector_expression<ES>& s)
579 SymMatrix Ptemp(X().size1(),X().size1());
587template<
class E>
inline
595 return prod( trans(
const_cast<ublas::matrix_expression<E>&
>(X) ), X);
598template<
class EX,
class ES,
class ET>
inline
600 prod_SPDT (
const ublas::matrix_expression<EX>& X,
const ublas::matrix_expression<ES>& S, ublas::matrix_expression<ET>& SXtemp)
606 return prod( trans(prod(S,X,SXtemp())), X);
609template<
class EX,
class ES,
class ET>
inline
610ET
prod_SPDT (
const ublas::matrix_expression<EX>& X,
const ublas::vector_expression<ES>& s, ublas::matrix_expression<ET>& Ptemp)
618 Vec::const_iterator si, send = s().end();
619 typename EX::const_iterator2 Xa = X.begin2();
620 const typename EX::const_iterator2 Xend = X.end2();
621 typename EX::const_iterator2 Xb;
624 for (; Xa != Xend; ++Xa)
626 typedef const ublas::matrix_column<const EX> EX_column;
627 EX_column Xav = ublas::column(XX, Xa.index2());
629 for (; Xb != Xend; ++Xb)
631 typename EX::value_type p = 0;
632 EX_column Xbv = ublas::column(XX, Xb.index2());
633 for (si = s().begin(); si != send; ++si) {
634 Vec::size_type i = si.index();
635 p += Xav[i] * (*si) * Xbv[i];
637 Ptemp()(Xa.index2(),Xb.index2()) = p;
638 Ptemp()(Xb.index2(),Xa.index2()) = p;
644template<
class EX,
class ES>
inline
651 SymMatrix Ptemp(X().size1(),X().size1());
663 const Vec::const_iterator xi_end = x.end();
664 Vec::const_iterator si = s.begin();
665 Vec::const_iterator xi=x.begin();
666 Vec::value_type p = Vec::value_type(0);
669 p += (*xi) * (*si) * (*xi);
681 return inner_prod (x, prod(S,x) );
Definition uBLASmatrix.hpp:49
BOOST_UBLAS_INLINE noalias_proxy(C &lval)
Definition uBLASmatrix.hpp:52
BOOST_UBLAS_INLINE void operator-=(const E &e)
Definition uBLASmatrix.hpp:69
BOOST_UBLAS_INLINE void operator=(const E &e)
Definition uBLASmatrix.hpp:57
BOOST_UBLAS_INLINE void operator+=(const E &e)
Definition uBLASmatrix.hpp:63
void mult_SPD(const MatrixX &X, const Vec &s, SymMatrix &P)
Definition uBLASmatrix.hpp:433
void mult_SPDT(const MatrixX &X, const Vec &s, SymMatrix &P)
Definition uBLASmatrix.hpp:462
FMMatrix< detail::BaseDiagMatrix > DiagMatrix
Definition uBLASmatrix.hpp:330
FMMatrix< detail::BaseColMatrix > ColMatrix
Definition uBLASmatrix.hpp:326
FMMatrix< detail::BaseRowMatrix > RowMatrix
Definition uBLASmatrix.hpp:324
FMVec< detail::BaseDenseVector > DenseVec
Definition uBLASmatrix.hpp:333
FMMatrix< detail::SymMatrixWrapper< detail::BaseRowMatrix > > SymMatrix
Definition uBLASmatrix.hpp:327
FMVec< detail::BaseVector > Vec
Definition uBLASmatrix.hpp:323
prod_expression_result< E, E >::E1TE2_type prod_SPDT(const ublas::matrix_expression< E > &X)
Definition uBLASmatrix.hpp:589
FMMatrix< detail::SymMatrixWrapper< detail::BaseDenseRowMatrix > > DenseSymMatrix
Definition uBLASmatrix.hpp:337
EmptyTag
Definition uBLASmatrix.hpp:34
@ Empty
Definition uBLASmatrix.hpp:34
void identity(FMMatrix< Base > &I)
Definition uBLASmatrix.hpp:412
FMMatrix< detail::BaseDenseLowerTriMatrix > DenseLTriMatrix
Definition uBLASmatrix.hpp:339
ublas::matrix_vector_range< FMMatrix< Base > > diag(FMMatrix< Base > &M, std::size_t n)
Definition uBLASmatrix.hpp:383
BOOST_UBLAS_INLINE detail::noalias_proxy< E > noalias(ublas::matrix_expression< E > &lvalue)
Definition uBLASmatrix.hpp:85
const ublas::triangular_adaptor< const M, ublas::upper > UpperTri(const M &m)
Definition uBLASmatrix.hpp:357
FMMatrix< detail::BaseDenseUpperTriMatrix > DenseUTriMatrix
Definition uBLASmatrix.hpp:338
FMMatrix< detail::BaseUpperTriMatrix > UTriMatrix
Definition uBLASmatrix.hpp:328
RowMatrix Matrix
Definition uBLASmatrix.hpp:325
FMMatrix< detail::BaseLowerTriMatrix > LTriMatrix
Definition uBLASmatrix.hpp:329
FMMatrix< detail::BaseDenseDiagMatrix > DenseDiagMatrix
Definition uBLASmatrix.hpp:340
FMMatrix< detail::BaseDenseRowMatrix > DenseRowMatrix
Definition uBLASmatrix.hpp:334
const ublas::triangular_adaptor< const M, ublas::lower > LowerTri(const M &m)
Definition uBLASmatrix.hpp:368
FMMatrix< detail::BaseDenseColMatrix > DenseColMatrix
Definition uBLASmatrix.hpp:336
prod_expression_result< E, E >::E1E2T_type prod_SPD(const ublas::matrix_expression< E > &X)
Definition uBLASmatrix.hpp:516
DenseRowMatrix DenseMatrix
Definition uBLASmatrix.hpp:335
Definition uBLASmatrix.hpp:502
BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits< BOOST_UBLAS_TYPENAMEE1T_type::value_type, E1T_type, BOOST_UBLAS_TYPENAMEE2::value_type, E2 >::result_type E1TE2_type
Definition uBLASmatrix.hpp:510
BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits< BOOST_UBLAS_TYPENAMEE1::value_type, E1, BOOST_UBLAS_TYPENAMEE2T_type::value_type, E2T_type >::result_type E1E2T_type
Definition uBLASmatrix.hpp:505
BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits< E1, ublas::scalar_identity< BOOST_UBLAS_TYPENAMEE1::value_type > >::result_type E1T_type
Definition uBLASmatrix.hpp:508
BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits< E2, ublas::scalar_identity< BOOST_UBLAS_TYPENAMEE2::value_type > >::result_type E2T_type
Definition uBLASmatrix.hpp:503
#define BOOST_UBLAS_TYPENAME
Definition uBLASmatrix.hpp:38