Bayes++ Bayesian Filtering Classes  Release 2014.5 - Copyright (c) 2003,2004,2005,2006,2011,2012,2014 Michael Stevens
uBLASmatrix.hpp
Go to the documentation of this file.
1 /*
2  * Bayes++ the Bayesian Filtering Library
3  * Copyright (c) 2002 Michael Stevens
4  * See accompanying Bayes++.htm for terms and conditions of use.
5  *
6  * $Id$
7  */
8 
9 /*
10  * Common type independent uBlas interface
11  * Should be include after base types have been defined
12  *
13  * Everything in namespace Bayes_filter_matrix is intended to support the matrix storage
14  * and algebra requirements of the library. Therefore the interfaces and implementation is
15  * not intended to be stable. Nor is this a general purpose adaptor for uBLAS
16  *
17  * Note on row_major matrices
18  * The implementation uses row major extensively. The computation of symmetric products P*P' is
19  * most efficient with row operations. These products are used extensively so the default
20  * is to use row_major matrices
21  */
22 
23 /* Filter Matrix Namespace */
24 namespace Bayesian_filter_matrix
25 {
26  // Allow use a few functions in own namespace (particular useful for compilers with Konig lookup)
27 using ublas::row;
28 using ublas::column;
29 using ublas::trans;
30 using ublas::prod; // These do not apply to the templated prod<temp> funtions
31 using ublas::inner_prod;
32 using ublas::outer_prod;
33 
34 enum EmptyTag {Empty}; // Tag type used for empty matrix constructor
35 
36  // Old compiler workaround removed from new uBLAS
37 #ifndef BOOST_UBLAS_TYPENAME
38 #define BOOST_UBLAS_TYPENAME typename
39 #endif
40 
41 #if (BOOST_VERSION >= 103100)
42 using ublas::noalias;
43 #else
44 namespace detail
45 {
46  // Assignment proxy.
47  // Provides temporary free assignment when LHS has no alias on RHS
48  template<class C>
49  class noalias_proxy {
50  public:
51  BOOST_UBLAS_INLINE
52  noalias_proxy (C& lval):
53  lval_ (lval) {}
54 
55  template <class E>
56  BOOST_UBLAS_INLINE
57  void operator= (const E& e) {
58  lval_.assign (e);
59  }
60 
61  template <class E>
62  BOOST_UBLAS_INLINE
63  void operator+= (const E& e) {
64  lval_.plus_assign (e);
65  }
66 
67  template <class E>
68  BOOST_UBLAS_INLINE
69  void operator-= (const E& e) {
70  lval_.minus_assign (e);
71  }
72 
73  private: // nonassignable
74  void operator=( const noalias_proxy& );
75  private:
76  C& lval_;
77  };
78 
79 }//namespace detail
80 
81 // Improve syntax of efficient assignment where no aliases of LHS appear on the RHS
82 // noalias(lhs) = rhs_expression
83 template <class E>
84 BOOST_UBLAS_INLINE
85 detail::noalias_proxy<E> noalias (ublas::matrix_expression<E>& lvalue) {
86  return detail::noalias_proxy<E> (lvalue() );
87 }
88 template <class E>
89 BOOST_UBLAS_INLINE
90 detail::noalias_proxy<E> noalias (ublas::vector_expression<E>& lvalue) {
91  return detail::noalias_proxy<E> (lvalue() );
92 }
93 
94 #endif
95 
96 
97 namespace detail // Lots of implementation detail
98 {
99 
100 /*
101  * Filter Vec type
102  */
103 template <class VecBase>
104 class FMVec : public VecBase
105 {
106 public:
107  typedef typename VecBase::value_type value_type;
108  typedef typename VecBase::vector_temporary_type vector_temporary_type;
109 
110  // No Default Constructor. Empty creation is very error prone
111  explicit FMVec(EmptyTag) : VecBase()
112  {} // Empty constructor
113  explicit FMVec(std::size_t size) : VecBase(size)
114  {} // Normal sized constructor
115  FMVec(const FMVec& c) : VecBase(static_cast<const VecBase&>(c))
116  {} // Copy constructor
117  template <class E>
118  explicit FMVec(const ublas::vector_expression<E>& e) : VecBase(e)
119  {} // vector_expression copy constructor
120  template <class E>
121  FMVec(const ublas::matrix_column<E>& e) : VecBase(e)
122  {} // conversion copy constructor, hides the implicit copy required for matrix column access
123 
124  template <class E>
125  FMVec& operator= (const ublas::vector_expression<E>& r)
126  { // Expression assignment; may be dependent on r
127  VecBase::operator=(r);
128  return *this;
129  }
130  FMVec& operator= (const FMVec& r)
131  { // Vector assignment; independent
132  VecBase::assign(r);
133  return *this;
134  }
135 
136  // Sub-range selection operators
137  const ublas::vector_range<const VecBase> sub_range(std::size_t b, std::size_t e) const
138  {
139  return ublas::vector_range<const VecBase>(*this, ublas::range(b,e));
140  }
141  ublas::vector_range<VecBase> sub_range(std::size_t b, std::size_t e)
142  {
143  return ublas::vector_range<VecBase>(*this, ublas::range(b,e));
144  }
145 };
146 
147 
148 /*
149  * Filter Matrix class template. Augmentation for uBlas MatrixBase
150  */
151 template <class MatrixBase>
152 class FMMatrix : public MatrixBase
153 {
154 public:
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;
158 
159  // No Default Constructor. Empty creation is very error prone
161  {} // Empty constructor
162  FMMatrix(std::size_t size1, std::size_t size2) : MatrixBase(size1,size2)
163  {} // Normal sized constructor
164  FMMatrix(const FMMatrix& c) : MatrixBase(static_cast<const MatrixBase&>(c))
165  {} // Copy constructor
166  template <class E>
167  explicit FMMatrix(const ublas::matrix_expression<E>& e) : MatrixBase(e)
168  {} // matrix_expression copy constructor
169 
170  template <class E>
171  FMMatrix& operator= (const ublas::matrix_expression<E>& r)
172  { // Expression assignment; may be dependent on r
173  MatrixBase::operator=(r);
174  return *this;
175  }
177  { // Matrix assignment; independent
178  MatrixBase::assign (r);
179  return *this;
180  }
181 
182  // Row,Column vector proxies
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;
187 
188  // Vector proxies from iterators - static members dependent on MatrixBase type
189  // ri() returns container associated with iterator. static_cast required as typeof(ri()) may not be typeof(MM)
190  static Row rowi(const typename MatrixBase::iterator1& ri)
191  {
192  return Row(static_cast<FMMatrix&>(ri()), ri.index1());
193  }
194  static const_Row rowi(const typename MatrixBase::const_iterator1& ri)
195  {
196  return const_Row(static_cast<const FMMatrix&>(ri()), ri.index1());
197  }
198  static Column columni(const typename MatrixBase::iterator2& ci)
199  {
200  return Column(static_cast<FMMatrix&>(ci()), ci.index2());
201  }
202  static const_Column columni(const typename MatrixBase::const_iterator2& ci)
203  {
204  return const_Column(static_cast<const FMMatrix&>(ci()), ci.index2());
205  }
206 
207  // Sub-range selection operators
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
210  {
211  return ublas::matrix_range<const MatrixBase> (*this, ublas::range(s1,e1), ublas::range(s2,e2));
212  }
213  ublas::matrix_range<MatrixBase>
214  sub_matrix(std::size_t s1, std::size_t e1, std::size_t s2, std::size_t e2)
215  {
216  return ublas::matrix_range<MatrixBase> (*this, ublas::range(s1,e1), ublas::range(s2,e2));
217  }
218 
219  // Requires boost_1.30.0 which has a generalised matrix_vector_slice
220  ublas::matrix_vector_slice<const MatrixBase>
221  sub_column(std::size_t s1, std::size_t e1, std::size_t s2) const
222  // Column vector s2 with rows [s1,e1)
223  {
224  return ublas::matrix_vector_slice<const MatrixBase> (*this, ublas::slice(s1,1,e1-s1), ublas::slice(s2,0,e1-s1));
225  }
226  ublas::matrix_vector_slice<MatrixBase>
227  sub_column(std::size_t s1, std::size_t e1, std::size_t s2)
228  // Column vector s2 with rows [s1,e1)
229  {
230  return ublas::matrix_vector_slice<MatrixBase> (*this, ublas::slice(s1,1,e1-s1), ublas::slice(s2,0,e1-s1));
231  }
232 };
233 
234 
235 /*
236  * Helper template to allow member construction before base class
237  * Boost version does not work as it passes by value
238  */
239 template <typename MemberType>
241 {
242 protected:
243  MemberType member;
244  explicit BaseFromMember() : member()
245  {}
246 
247  template <typename T1>
248  explicit BaseFromMember( const T1& x1 ) : member( x1 )
249  {}
250 
251  template <typename T1, typename T2>
252  explicit BaseFromMember( const T1& x1, const T2& x2 ) : member( x1, x2 )
253  {}
254 };
255 
256 
257 /*
258  * We require static type conversion between Symmetric matrices and equivalent row major matrices
259  * Therefore we create symmetric matrix types, using a MatrixBase for storage
260  * and wraps this in a symmetric_adaptor
261  */
262 template <class MatrixBase>
264  private BaseFromMember<MatrixBase>, // allow construction of MatrixBase member before symmetric_adaptor
265  public ublas::symmetric_adaptor<MatrixBase, ublas::upper>
266 {
268  typedef ublas::symmetric_adaptor<MatrixBase, ublas::upper> symadaptor_type;
269 public:
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;
273 
274  SymMatrixWrapper () : matrix_type(), symadaptor_type(matrix_type::member)
275  {}
276  SymMatrixWrapper (std::size_t size1, std::size_t size2) : matrix_type(size1,size2), symadaptor_type(matrix_type::member)
277  {} // Normal sized constructor
278  explicit SymMatrixWrapper (const SymMatrixWrapper& r) : matrix_type(reinterpret_cast<const MatrixBase&>(r)), symadaptor_type(matrix_type::member)
279  {} // Explicit copy construction referencing the copy reinterpreted as a MatrixBase
280  template <class E>
281  explicit SymMatrixWrapper (const ublas::matrix_expression<E>& e) : matrix_type(e), symadaptor_type(matrix_type::member)
282  {} // Explicit matrix_expression conversion constructor
283 
284  template <class E>
285  SymMatrixWrapper& operator=(const ublas::matrix_expression<E>& r)
286  {
287  symadaptor_type::operator=(r);
288  return *this;
289  }
290 
291  // Conversions straight to a FMMatrix, equivalent to a RowMatrix types
293  {
294  return static_cast<const FMMatrix<MatrixBase>& >(matrix_type::member);
295  }
297  {
298  return static_cast<FMMatrix<MatrixBase>& >(matrix_type::member);
299  }
300 
301  // Matrix storage members
302  void clear()
303  { matrix_type::member.clear();
304  }
305  void resize(std::size_t nsize1, std::size_t nsize2, bool preserve = true)
306  {
307  matrix_type::member.resize(nsize1, nsize2, preserve);
308  }
309 };
310 
311 }//namespace detail
312 
313 
314 
315 /*
316  * Vector / Matrix types
317  * Finally the definitions !
318  */
319 using detail::FMVec; // Template class for template parameter matching
320 using detail::FMMatrix;
321 
322  // Default types
323 typedef FMVec<detail::BaseVector> Vec;
324 typedef FMMatrix<detail::BaseRowMatrix> RowMatrix;
325 typedef RowMatrix Matrix;
326 typedef FMMatrix<detail::BaseColMatrix> ColMatrix;
327 typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseRowMatrix> > SymMatrix;
328 typedef FMMatrix<detail::BaseUpperTriMatrix> UTriMatrix;
329 typedef FMMatrix<detail::BaseLowerTriMatrix> LTriMatrix;
330 typedef FMMatrix<detail::BaseDiagMatrix> DiagMatrix;
331 
332  // Explicitly dense types
333 typedef FMVec<detail::BaseDenseVector> DenseVec;
334 typedef FMMatrix<detail::BaseDenseRowMatrix> DenseRowMatrix;
335 typedef DenseRowMatrix DenseMatrix;
336 typedef FMMatrix<detail::BaseDenseColMatrix> DenseColMatrix;
337 typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseDenseRowMatrix> > DenseSymMatrix;
338 typedef FMMatrix<detail::BaseDenseUpperTriMatrix> DenseUTriMatrix;
339 typedef FMMatrix<detail::BaseDenseLowerTriMatrix> DenseLTriMatrix;
340 typedef FMMatrix<detail::BaseDenseDiagMatrix> DenseDiagMatrix;
341 
342  // Explicitly sparse types (any of the gappy types)
343 #ifdef BAYES_FILTER_GAPPY
344 typedef FMVec<detail::BaseSparseVector> SparseVec;
345 typedef FMMatrix<detail::BaseDenseRowMatrix> SparseRowMatrix;
346 typedef SparseRowMatrix SparseMatrix;
347 typedef FMMatrix<detail::BaseSparseColMatrix> SparseColMatrix;
348 typedef FMMatrix<detail::SymMatrixWrapper<detail::BaseSparseRowMatrix> > SparseSymMatrix;
349 #endif
350 
351 
352 /*
353  * Matrix Adaptors, simply hide the uBLAS details
354  */
355 template <class M>
356 const ublas::triangular_adaptor<const M, ublas::upper>
357  UpperTri(const M& m)
358 /*
359  * View Upper triangle of m
360  * ISSUE VC7 cannot cope with UTriMatrix::functor1_type
361  */
362 {
363  return ublas::triangular_adaptor<const M, ublas::upper>(m);
364 }
365 
366 template <class M>
367 const ublas::triangular_adaptor<const M, ublas::lower>
368  LowerTri(const M& m)
369 /*
370  * View Lower triangle of m
371  */
372 {
373  return ublas::triangular_adaptor<const M, ublas::lower>(m);
374 }
375 
376 
377 
378 /*
379  * Matrix Support Operations
380  */
381 template <class Base>
382 ublas::matrix_vector_range<FMMatrix<Base> >
383  diag(FMMatrix<Base>& M, std::size_t n)
384 { // Return a vector proxy to the first n diagonal elements of M
385  return ublas::matrix_vector_range<FMMatrix<Base> >(M, ublas::range(0,n), ublas::range(0,n));
386 }
387 
388 template <class Base>
389 const ublas::matrix_vector_range<const FMMatrix<Base> >
390  diag(const FMMatrix<Base>& M, std::size_t n)
391 { // Return a vector proxy to the first n diagonal elements of M
392  return ublas::matrix_vector_range<const FMMatrix<Base> >(M, ublas::range(0,n), ublas::range(0,n));
393 }
394 
395 template <class Base>
396 ublas::matrix_vector_range<FMMatrix<Base> >
397  diag(FMMatrix<Base>& M)
398 { // Return a vector proxy to the diagonal elements of M
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));
401 }
402 
403 template <class Base>
404 const ublas::matrix_vector_range<const FMMatrix<Base> >
405  diag(const FMMatrix<Base>& M)
406 { // Return a vector proxy to the diagonal elements of M
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));
409 }
410 
411 template <class Base>
412 void identity(FMMatrix<Base>& I)
413 { // Set I to generalised Identity matrix. Clear I and set diag(I) to one
414  I.clear();
415  // Set common diagonal elements
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)));
419 }
420 
421 
422 
423 /*
424  * Symmetric Positive (Semi) Definite multiplication: X*S*X' and X'*S*X
425  * The name is slightly misleading. The result is actually only PD if S is PD
426  * Algorithms are intended to exploit the symmetry of the result
427  * and also where possible the row by row multiplication inherent in X*X'
428  */
429 
430 namespace detail { // mult_SPD now an implementation detail
431 
432 template <class MatrixX>
433 void mult_SPD (const MatrixX& X, const Vec& s, SymMatrix& P)
434 /*
435  * Symmetric Positive (Semi) Definite multiply: P += X*diag_matrix(s)*X'
436  */
437 {
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;
442 
443  // P(a,b) = X.row(a) * X.row(b)
444  for (; Xa != Xend; ++Xa) // Iterate Rows
445  {
446  typename MatrixX::const_Row Xav = MatrixX::rowi(Xa);
447  Xb = Xa; // Start at the row Xa only one triangle of symmetric result required
448  for (; Xb != Xend; ++Xb)
449  {
450  SymMatrix::value_type p = 0; // Triple vector inner product
451  typename MatrixX::const_Row Xbv = MatrixX::rowi(Xb);
452  for (si = s.begin(); si != send; ++si) {
453  Vec::size_type i = si.index();
454  p += Xav[i] * (*si) * Xbv[i];
455  }
456  P(Xa.index1(),Xb.index1()) += p;
457  }
458  }
459 }
460 
461 template <class MatrixX>
462 void mult_SPDT (const MatrixX& X, const Vec& s, SymMatrix& P)
463 /*
464  * Symmetric Positive (Semi) Definite multiply: P += X'*diag_matrix(s)*X
465  */
466 {
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;
471 
472  // P(a,b) = X.col(a) * X.col(b)
473  for (; Xa != Xend; ++Xa) // Iterate vectors
474  {
475  typename MatrixX::const_Column Xav = MatrixX::columni(Xa);
476  Xb = Xa; // Start at the row Xa only one triangle of symmetric result required
477  for (; Xb != Xend; ++Xb)
478  { // Triple vector inner product
479  SymMatrix::value_type p = 0;
480  typename MatrixX::const_Column Xbv = MatrixX::columni(Xb);
481  for (si = s.begin(); si != send; ++si) {
482  Vec::size_type i = si.index();
483  p += Xav[i] * (*si) * Xbv[i];
484  }
485  P(Xa.index2(),Xb.index2()) += p;
486  }
487  }
488 }
489 
490 }//namespace detail
491 
492 
493 /*
494  * prod_SPD - uBLAS Expression templates for X*X' and X*X'
495  * Functions are only defined for the type for which the operation is efficient
496  * ISSUE Although numerically symmetric, uBlas has no expression type to represent this property
497  * The result must be assigned to a symmetric container to exploit the symmetry
498  */
499 
500 template <class E1, class E2>
502 { // Provide ET result E1E2T_type of prod(matrix_expression<E1>,trans(matrix_expression<E2>)
503  typedef BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits<E2, ublas::scalar_identity<BOOST_UBLAS_TYPENAME E2::value_type> >::result_type E2T_type;
504  typedef BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits<BOOST_UBLAS_TYPENAME E1::value_type, E1,
505  BOOST_UBLAS_TYPENAME E2T_type::value_type, E2T_type>::result_type E1E2T_type;
506 
507  // Provide ET result E1TE2_type of prod(trans(matrix_expression<E1>),matrix_expression<E2>)
508  typedef BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits<E1, ublas::scalar_identity<BOOST_UBLAS_TYPENAME E1::value_type> >::result_type E1T_type;
509  typedef BOOST_UBLAS_TYPENAME ublas::matrix_matrix_binary_traits<BOOST_UBLAS_TYPENAME E1T_type::value_type, E1T_type,
510  BOOST_UBLAS_TYPENAME E2::value_type, E2>::result_type E1TE2_type;
511 };
512 
513 
514 template<class E> inline
516  prod_SPD (const ublas::matrix_expression<E>& X)
517 /*
518  * Symmetric Positive (Semi) Definite product: X*X'
519  */
520 {
521  // ISSUE: uBLAS post Boost 1_31_0 introduces a trans(const matrix_expression<E>& e) which propagates the const expression type
522  // Bypass this to avoid having to detect the Boost version
523  return prod( X, trans(const_cast<ublas::matrix_expression<E>&>(X) ));
524 }
525 
526 template<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)
529 /*
530  * Symmetric Positive (Semi) Definite product: X*(X*S)', XStemp = X*S
531  * Result symmetric if S is symmetric
532  */
533 {
534  return prod( X, trans(prod(X,S,XStemp())) );
535 }
536 
537 template<class EX, class ES, class ET> inline
538 ET prod_SPD (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s, ublas::matrix_expression<ET>& Ptemp)
539 /*
540  * Symmetric Positive (Semi) Definite product: X*diag_matrix(s)*X'
541  * Precond: Ptemp must be size conformant with the product
542  * DEPRECATED With explicit Ptemp, do not rely on it's value
543  */
544 {
545  Vec::const_iterator si, send = s().end();
546  const EX& XX = X();
547  typename EX::const_iterator1 Xa = XX.begin1();
548  const typename EX::const_iterator1 Xend = XX.end1();
549  typename EX::const_iterator1 Xb;
550 
551  // P(a,b) = sum(X.row(a) * s * X.row(b))
552  for (; Xa != Xend; ++Xa) // Iterate rows
553  {
554  typedef const ublas::matrix_row<const EX> EX_row;
555  EX_row Xav (ublas::row(XX, Xa.index1()));
556  Xb = Xa; // Start at the row Xa only one triangle of symmetric result required
557  for (; Xb != Xend; ++Xb)
558  {
559  typename EX::value_type p = 0; // Triple vector inner product
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];
564  }
565  Ptemp()(Xa.index1(),Xb.index1()) = p;
566  Ptemp()(Xb.index1(),Xa.index1()) = p;
567  }
568  }
569  return Ptemp();
570 }
571 
572 template<class EX, class ES> inline
573 SymMatrix prod_SPD (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s)
574 /*
575  * Symmetric Positive (Semi) Definite product: X*diag_matrix(s)*X'
576  * TODO requires a prod_diag(X,s)
577  */
578 {
579  SymMatrix Ptemp(X().size1(),X().size1());
580  Ptemp.clear();
581 
582  (void) prod_SPD (X(),s(), Ptemp);
583  return Ptemp;
584 }
585 
586 
587 template<class E> inline
589  prod_SPDT (const ublas::matrix_expression<E>& X)
590 /*
591  * Symmetric Positive (Semi) Definite product: X'*X
592  */
593 {
594  // ISSUE: See prod_SPD
595  return prod( trans(const_cast<ublas::matrix_expression<E>&>(X) ), X);
596 }
597 
598 template<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)
601 /*
602  * Symmetric Positive (Semi) Definite product: (S*X)'*X, SXtemp = S*X
603  * Result symmetric if S is symmetric
604  */
605 {
606  return prod( trans(prod(S,X,SXtemp())), X);
607 }
608 
609 template<class EX, class ES, class ET> inline
610 ET prod_SPDT (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s, ublas::matrix_expression<ET>& Ptemp)
611 /*
612  * Symmetric Positive (Semi) Definite product: X'*diag_matrix(s)*X, Ptemp = return value
613  * Precond: Ptemp must be size conformant with the product
614  * DEPRECATED With explicit Ptemp, do not rely on it's value
615  */
616 {
617  const EX& XX = X();
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;
622 
623  // P(a,b) = sum(X.col(a) * s * X.col(b))
624  for (; Xa != Xend; ++Xa) // Iterate columns
625  {
626  typedef const ublas::matrix_column<const EX> EX_column;
627  EX_column Xav = ublas::column(XX, Xa.index2());
628  Xb = Xa; // Start at the column Xa only one triangle of symmetric result required
629  for (; Xb != Xend; ++Xb)
630  {
631  typename EX::value_type p = 0; // Triple vector inner product
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];
636  }
637  Ptemp()(Xa.index2(),Xb.index2()) = p;
638  Ptemp()(Xb.index2(),Xa.index2()) = p;
639  }
640  }
641  return Ptemp();
642 }
643 
644 template<class EX, class ES> inline
645 SymMatrix prod_SPDT (const ublas::matrix_expression<EX>& X, const ublas::vector_expression<ES>& s)
646 /*
647  * Symmetric Positive (Semi) Definite product: X'*diag_matrix(s)*X
648  * TODO requires a prod_diag(X,s)
649  */
650 {
651  SymMatrix Ptemp(X().size1(),X().size1());
652  Ptemp.clear();
653 
654  (void) prod_SPDT (X(),s(), Ptemp);
655  return Ptemp;
656 }
657 
658 inline Vec::value_type prod_SPDT (const Vec& x, const Vec& s)
659 /*
660  * Symmetric Positive (Semi) Definite product: x'*diag_matrix(s)*x
661  */
662 {
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);
667  while (xi != xi_end)
668  {
669  p += (*xi) * (*si) * (*xi);
670  ++xi; ++ si;
671  }
672 
673  return p;
674 }
675 
676 inline Vec::value_type prod_SPDT (const Vec& x, const SymMatrix& S)
677 /*
678  * Symmetric Positive (Semi) Definite multiply: p = x'*S*x
679  */
680 {
681  return inner_prod (x, prod(S,x) );
682 }
683 
684 }//namespace
FMVec< detail::BaseDenseVector > DenseVec
Definition: uBLASmatrix.hpp:333
FMMatrix(const FMMatrix &c)
Definition: uBLASmatrix.hpp:164
void clear()
Definition: uBLASmatrix.hpp:302
FMVec(EmptyTag)
Definition: uBLASmatrix.hpp:111
MatrixBase::vector_temporary_type vector_temporary_type
Definition: uBLASmatrix.hpp:271
SymMatrixWrapper(const SymMatrixWrapper &r)
Definition: uBLASmatrix.hpp:278
ublas::matrix_range< const MatrixBase > sub_matrix(std::size_t s1, std::size_t e1, std::size_t s2, std::size_t e2) const
Definition: uBLASmatrix.hpp:209
static const_Row rowi(const typename MatrixBase::const_iterator1 &ri)
Definition: uBLASmatrix.hpp:194
void resize(std::size_t nsize1, std::size_t nsize2, bool preserve=true)
Definition: uBLASmatrix.hpp:305
FMMatrix(EmptyTag)
Definition: uBLASmatrix.hpp:160
BOOST_UBLAS_INLINE void operator=(const E &e)
Definition: uBLASmatrix.hpp:57
prod_expression_result< E, E >::E1E2T_type prod_SPD(const ublas::matrix_expression< E > &X)
Definition: uBLASmatrix.hpp:516
void mult_SPD(const MatrixX &X, const Vec &s, SymMatrix &P)
Definition: uBLASmatrix.hpp:433
FMMatrix< detail::BaseUpperTriMatrix > UTriMatrix
Definition: uBLASmatrix.hpp:328
FMMatrix< detail::BaseDenseDiagMatrix > DenseDiagMatrix
Definition: uBLASmatrix.hpp:340
BaseFromMember()
Definition: uBLASmatrix.hpp:244
void identity(FMMatrix< Base > &I)
Definition: uBLASmatrix.hpp:412
ublas::matrix_vector_range< FMMatrix< Base > > diag(FMMatrix< Base > &M, std::size_t n)
Definition: uBLASmatrix.hpp:383
MemberType member
Definition: uBLASmatrix.hpp:243
VecBase::value_type value_type
Definition: uBLASmatrix.hpp:107
Definition: uBLASmatrix.hpp:34
MatrixBase::matrix_temporary_type matrix_temporary_type
Definition: uBLASmatrix.hpp:157
FMMatrix< detail::BaseDiagMatrix > DiagMatrix
Definition: uBLASmatrix.hpp:330
BaseFromMember(const T1 &x1)
Definition: uBLASmatrix.hpp:248
Definition: uBLASmatrix.hpp:152
FMMatrix< MatrixBase > & asRowMatrix()
Definition: uBLASmatrix.hpp:296
MatrixBase::matrix_temporary_type matrix_temporary_type
Definition: uBLASmatrix.hpp:272
SymMatrixWrapper(std::size_t size1, std::size_t size2)
Definition: uBLASmatrix.hpp:276
Definition: uBLASmatrix.hpp:501
MatrixBase::value_type value_type
Definition: uBLASmatrix.hpp:155
FMMatrix(const ublas::matrix_expression< E > &e)
Definition: uBLASmatrix.hpp:167
static Row rowi(const typename MatrixBase::iterator1 &ri)
Definition: uBLASmatrix.hpp:190
BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits< E2, ublas::scalar_identity< BOOST_UBLAS_TYPENAME E2::value_type > >::result_type E2T_type
Definition: uBLASmatrix.hpp:503
ublas::matrix_column< FMMatrix > Column
Definition: uBLASmatrix.hpp:185
static const_Column columni(const typename MatrixBase::const_iterator2 &ci)
Definition: uBLASmatrix.hpp:202
BOOST_UBLAS_INLINE void operator+=(const E &e)
Definition: uBLASmatrix.hpp:63
DenseRowMatrix DenseMatrix
Definition: uBLASmatrix.hpp:335
BaseFromMember(const T1 &x1, const T2 &x2)
Definition: uBLASmatrix.hpp:252
const ublas::matrix_row< const FMMatrix > const_Row
Definition: uBLASmatrix.hpp:184
MatrixBase::vector_temporary_type vector_temporary_type
Definition: uBLASmatrix.hpp:156
FMVec(const ublas::matrix_column< E > &e)
Definition: uBLASmatrix.hpp:121
ublas::matrix_range< MatrixBase > sub_matrix(std::size_t s1, std::size_t e1, std::size_t s2, std::size_t e2)
Definition: uBLASmatrix.hpp:214
const ublas::triangular_adaptor< const M, ublas::upper > UpperTri(const M &m)
Definition: uBLASmatrix.hpp:357
FMVec< detail::BaseVector > Vec
Definition: uBLASmatrix.hpp:323
BOOST_UBLAS_INLINE detail::noalias_proxy< E > noalias(ublas::matrix_expression< E > &lvalue)
Definition: uBLASmatrix.hpp:85
ublas::vector_range< VecBase > sub_range(std::size_t b, std::size_t e)
Definition: uBLASmatrix.hpp:141
static Column columni(const typename MatrixBase::iterator2 &ci)
Definition: uBLASmatrix.hpp:198
const FMMatrix< MatrixBase > & asRowMatrix() const
Definition: uBLASmatrix.hpp:292
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
Definition: uBLASmatrix.hpp:510
void mult_SPDT(const MatrixX &X, const Vec &s, SymMatrix &P)
Definition: uBLASmatrix.hpp:462
Definition: uBLASmatrix.hpp:49
RowMatrix Matrix
Definition: uBLASmatrix.hpp:325
const ublas::matrix_column< const FMMatrix > const_Column
Definition: uBLASmatrix.hpp:186
SymMatrixWrapper & operator=(const ublas::matrix_expression< E > &r)
Definition: uBLASmatrix.hpp:285
const ublas::triangular_adaptor< const M, ublas::lower > LowerTri(const M &m)
Definition: uBLASmatrix.hpp:368
FMMatrix(std::size_t size1, std::size_t size2)
Definition: uBLASmatrix.hpp:162
BOOST_UBLAS_INLINE detail::noalias_proxy< E > noalias(ublas::vector_expression< E > &lvalue)
Definition: uBLASmatrix.hpp:90
FMVec(std::size_t size)
Definition: uBLASmatrix.hpp:113
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
Definition: uBLASmatrix.hpp:505
BOOST_UBLAS_INLINE noalias_proxy(C &lval)
Definition: uBLASmatrix.hpp:52
Definition: uBLASmatrix.hpp:104
prod_expression_result< E, E >::E1TE2_type prod_SPDT(const ublas::matrix_expression< E > &X)
Definition: uBLASmatrix.hpp:589
FMMatrix< detail::BaseDenseColMatrix > DenseColMatrix
Definition: uBLASmatrix.hpp:336
ublas::matrix_vector_slice< MatrixBase > sub_column(std::size_t s1, std::size_t e1, std::size_t s2)
Definition: uBLASmatrix.hpp:227
FMMatrix< detail::BaseDenseUpperTriMatrix > DenseUTriMatrix
Definition: uBLASmatrix.hpp:338
VecBase::vector_temporary_type vector_temporary_type
Definition: uBLASmatrix.hpp:108
FMMatrix< detail::SymMatrixWrapper< detail::BaseRowMatrix > > SymMatrix
Definition: uBLASmatrix.hpp:327
FMMatrix< detail::BaseDenseRowMatrix > DenseRowMatrix
Definition: uBLASmatrix.hpp:334
Definition: matSup.cpp:33
FMMatrix< detail::BaseDenseLowerTriMatrix > DenseLTriMatrix
Definition: uBLASmatrix.hpp:339
FMMatrix< detail::BaseLowerTriMatrix > LTriMatrix
Definition: uBLASmatrix.hpp:329
#define BOOST_UBLAS_TYPENAME
Definition: uBLASmatrix.hpp:38
SymMatrixWrapper()
Definition: uBLASmatrix.hpp:274
FMMatrix< detail::BaseColMatrix > ColMatrix
Definition: uBLASmatrix.hpp:326
MatrixBase::value_type value_type
Definition: uBLASmatrix.hpp:270
BOOST_UBLAS_INLINE void operator-=(const E &e)
Definition: uBLASmatrix.hpp:69
ublas::matrix_row< FMMatrix > Row
Definition: uBLASmatrix.hpp:183
EmptyTag
Definition: uBLASmatrix.hpp:34
FMMatrix< detail::BaseRowMatrix > RowMatrix
Definition: uBLASmatrix.hpp:324
Definition: uBLASmatrix.hpp:240
const ublas::vector_range< const VecBase > sub_range(std::size_t b, std::size_t e) const
Definition: uBLASmatrix.hpp:137
FMVec(const ublas::vector_expression< E > &e)
Definition: uBLASmatrix.hpp:118
BOOST_UBLAS_TYPENAME ublas::matrix_unary2_traits< E1, ublas::scalar_identity< BOOST_UBLAS_TYPENAME E1::value_type > >::result_type E1T_type
Definition: uBLASmatrix.hpp:508
FMVec(const FMVec &c)
Definition: uBLASmatrix.hpp:115
FMMatrix< detail::SymMatrixWrapper< detail::BaseDenseRowMatrix > > DenseSymMatrix
Definition: uBLASmatrix.hpp:337
SymMatrixWrapper(const ublas::matrix_expression< E > &e)
Definition: uBLASmatrix.hpp:281
ublas::matrix_vector_slice< const MatrixBase > sub_column(std::size_t s1, std::size_t e1, std::size_t s2) const
Definition: uBLASmatrix.hpp:221