SHOGUN  6.0.0
LinalgBackendEigen.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2016, Shogun-Toolbox e.V. <shogun-team@shogun-toolbox.org>
3  * All rights reserved.
4  * Redistribution and use in source and binary forms, with or without
5  * modification, are permitted provided that the following conditions are met:
6  *
7  * 1. Redistributions of source code must retain the above copyright notice,
8  * this list of conditions and the following disclaimer.
9  *
10  * 2. Redistributions in binary form must reproduce the above copyright
11  * notice, this list of conditions and the following disclaimer in the
12  * documentation and/or other materials provided with the distribution.
13  *
14  * 3. Neither the name of the copyright holder nor the names of its
15  * contributors may be used to endorse or promote products derived from
16  * this software without specific prior written permission.
17  *
18  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
22  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
28  * POSSIBILITY OF SUCH DAMAGE.
29  *
30  * Authors: 2016 Pan Deng, Soumyajit De, Heiko Strathmann, Viktor Gal
31  */
32 
33 #ifndef LINALG_BACKEND_EIGEN_H__
34 #define LINALG_BACKEND_EIGEN_H__
35 
36 #include <shogun/lib/SGVector.h>
39 #include <numeric>
40 
41 namespace shogun
42 {
43 
46 {
47 public:
48  #define DEFINE_FOR_ALL_PTYPE(METHODNAME, Container) \
49  METHODNAME(bool, Container); \
50  METHODNAME(char, Container); \
51  METHODNAME(int8_t, Container); \
52  METHODNAME(uint8_t, Container); \
53  METHODNAME(int16_t, Container); \
54  METHODNAME(uint16_t, Container); \
55  METHODNAME(int32_t, Container); \
56  METHODNAME(uint32_t, Container); \
57  METHODNAME(int64_t, Container); \
58  METHODNAME(uint64_t, Container); \
59  METHODNAME(float32_t, Container); \
60  METHODNAME(float64_t, Container); \
61  METHODNAME(floatmax_t, Container); \
62  METHODNAME(complex128_t, Container); \
63 
64  #define DEFINE_FOR_REAL_PTYPE(METHODNAME, Container) \
65  METHODNAME(bool, Container); \
66  METHODNAME(char, Container); \
67  METHODNAME(int8_t, Container); \
68  METHODNAME(uint8_t, Container); \
69  METHODNAME(int16_t, Container); \
70  METHODNAME(uint16_t, Container); \
71  METHODNAME(int32_t, Container); \
72  METHODNAME(uint32_t, Container); \
73  METHODNAME(int64_t, Container); \
74  METHODNAME(uint64_t, Container); \
75  METHODNAME(float32_t, Container); \
76  METHODNAME(float64_t, Container); \
77  METHODNAME(floatmax_t, Container);
78 
79  #define DEFINE_FOR_NON_INTEGER_PTYPE(METHODNAME, Container) \
80  METHODNAME(float32_t, Container); \
81  METHODNAME(float64_t, Container); \
82  METHODNAME(floatmax_t, Container); \
83  METHODNAME(complex128_t, Container);
84 
85  #define DEFINE_FOR_NUMERIC_PTYPE(METHODNAME, Container) \
86  METHODNAME(char, Container); \
87  METHODNAME(int8_t, Container); \
88  METHODNAME(uint8_t, Container); \
89  METHODNAME(int16_t, Container); \
90  METHODNAME(uint16_t, Container); \
91  METHODNAME(int32_t, Container); \
92  METHODNAME(uint32_t, Container); \
93  METHODNAME(int64_t, Container); \
94  METHODNAME(uint64_t, Container); \
95  METHODNAME(float32_t, Container); \
96  METHODNAME(float64_t, Container); \
97  METHODNAME(floatmax_t, Container);
98 
100  #define BACKEND_GENERIC_IN_PLACE_ADD(Type, Container) \
101  virtual void add(Container<Type>& a, Container<Type>& b, Type alpha, \
102  Type beta, Container<Type>& result) const \
103  { \
104  add_impl(a, b, alpha, beta, result); \
105  }
108  #undef BACKEND_GENERIC_IN_PLACE_ADD
109 
111  #define BACKEND_GENERIC_CHOLESKY_FACTOR(Type, Container) \
112  virtual Container<Type> cholesky_factor(const Container<Type>& A, \
113  const bool lower) const \
114  { \
115  return cholesky_factor_impl(A, lower); \
116  }
118  #undef BACKEND_GENERIC_CHOLESKY_FACTOR
119 
121  #define BACKEND_GENERIC_CHOLESKY_SOLVER(Type, Container) \
122  virtual SGVector<Type> cholesky_solver(const Container<Type>& L, \
123  const SGVector<Type>& b, const bool lower) const \
124  { \
125  return cholesky_solver_impl(L, b, lower); \
126  }
128  #undef BACKEND_GENERIC_CHOLESKY_SOLVER
129 
131  #define BACKEND_GENERIC_DOT(Type, Container) \
132  virtual Type dot(const Container<Type>& a, const Container<Type>& b) const \
133  { \
134  return dot_impl(a, b); \
135  }
137  #undef BACKEND_GENERIC_DOT
138 
140  #define BACKEND_GENERIC_IN_PLACE_ELEMENT_PROD(Type, Container) \
141  virtual void element_prod(Container<Type>& a, Container<Type>& b,\
142  Container<Type>& result) const \
143  { \
144  element_prod_impl(a, b, result); \
145  }
147  #undef BACKEND_GENERIC_IN_PLACE_ELEMENT_PROD
148 
150  #define BACKEND_GENERIC_IN_PLACE_BLOCK_ELEMENT_PROD(Type, Container) \
151  virtual void element_prod(linalg::Block<Container<Type>>& a, \
152  linalg::Block<Container<Type>>& b, Container<Type>& result) const \
153  { \
154  element_prod_impl(a, b, result); \
155  }
157  #undef BACKEND_GENERIC_IN_PLACE_BLOCK_ELEMENT_PROD
158 
160  #define BACKEND_GENERIC_LOGISTIC(Type, Container) \
161  virtual void logistic(Container<Type>& a, Container<Type>& result) const \
162  { \
163  logistic_impl(a, result); \
164  }
166  #undef BACKEND_GENERIC_LOGISTIC
167 
169  #define BACKEND_GENERIC_IN_PLACE_MATRIX_PROD(Type, Container) \
170  virtual void matrix_prod(SGMatrix<Type>& a, Container<Type>& b,\
171  Container<Type>& result, bool transpose_A, bool transpose_B) const \
172  { \
173  matrix_prod_impl(a, b, result, transpose_A, transpose_B); \
174  }
177  #undef BACKEND_GENERIC_IN_PLACE_MATRIX_PROD
178 
180  #define BACKEND_GENERIC_MAX(Type, Container) \
181  virtual Type max(const Container<Type>& a) const \
182  { \
183  return max_impl(a); \
184  }
187  #undef BACKEND_GENERIC_MAX
188 
190  #define BACKEND_GENERIC_REAL_MEAN(Type, Container) \
191  virtual float64_t mean(const Container<Type>& a) const \
192  { \
193  return mean_impl(a); \
194  }
197  #undef BACKEND_GENERIC_REAL_MEAN
198 
200  #define BACKEND_GENERIC_COMPLEX_MEAN(Container) \
201  virtual complex128_t mean(const Container<complex128_t>& a) const \
202  { \
203  return mean_impl(a); \
204  }
207  #undef BACKEND_GENERIC_COMPLEX_MEAN
208 
210  #define BACKEND_GENERIC_RANGE_FILL(Type, Container) \
211  virtual void range_fill(Container<Type>& a, const Type start) const \
212  { \
213  range_fill_impl(a, start); \
214  }
217  #undef BACKEND_GENERIC_RANGE_FILL
218 
220  #define BACKEND_GENERIC_IN_PLACE_SCALE(Type, Container) \
221  virtual void scale(Container<Type>& a, Type alpha, Container<Type>& result) const \
222  { \
223  scale_impl(a, alpha, result); \
224  }
227  #undef BACKEND_GENERIC_IN_PLACE_SCALE
228 
230  #define BACKEND_GENERIC_SET_CONST(Type, Container) \
231  virtual void set_const(Container<Type>& a, const Type value) const \
232  { \
233  set_const_impl(a, value); \
234  }
237  #undef BACKEND_GENERIC_SET_CONST
238 
240  #define BACKEND_GENERIC_SUM(Type, Container) \
241  virtual Type sum(const Container<Type>& a, bool no_diag) const \
242  { \
243  return sum_impl(a, no_diag); \
244  }
247  #undef BACKEND_GENERIC_SUM
248 
250  #define BACKEND_GENERIC_BLOCK_SUM(Type, Container) \
251  virtual Type sum(const linalg::Block<Container<Type>>& a, bool no_diag) const \
252  { \
253  return sum_impl(a, no_diag); \
254  }
256  #undef BACKEND_GENERIC_BLOCK_SUM
257 
259  #define BACKEND_GENERIC_SYMMETRIC_SUM(Type, Container) \
260  virtual Type sum_symmetric(const Container<Type>& a, bool no_diag) const \
261  { \
262  return sum_symmetric_impl(a, no_diag); \
263  }
265  #undef BACKEND_GENERIC_SYMMETRIC_SUM
266 
268  #define BACKEND_GENERIC_SYMMETRIC_BLOCK_SUM(Type, Container) \
269  virtual Type sum_symmetric(const linalg::Block<Container<Type>>& a, bool no_diag) const \
270  { \
271  return sum_symmetric_impl(a, no_diag); \
272  }
274  #undef BACKEND_GENERIC_SYMMETRIC_BLOCK_SUM
275 
277  #define BACKEND_GENERIC_COLWISE_SUM(Type, Container) \
278  virtual SGVector<Type> colwise_sum(const Container<Type>& a, bool no_diag) const \
279  { \
280  return colwise_sum_impl(a, no_diag); \
281  }
283  #undef BACKEND_GENERIC_COLWISE_SUM
284 
286  #define BACKEND_GENERIC_BLOCK_COLWISE_SUM(Type, Container) \
287  virtual SGVector<Type> colwise_sum(const linalg::Block<Container<Type>>& a, bool no_diag) const \
288  { \
289  return colwise_sum_impl(a, no_diag); \
290  }
292  #undef BACKEND_GENERIC_BLOCK_COLWISE_SUM
293 
295  #define BACKEND_GENERIC_ROWWISE_SUM(Type, Container) \
296  virtual SGVector<Type> rowwise_sum(const Container<Type>& a, bool no_diag) const \
297  { \
298  return rowwise_sum_impl(a, no_diag); \
299  }
301  #undef BACKEND_GENERIC_ROWWISE_SUM
302 
304  #define BACKEND_GENERIC_BLOCK_ROWWISE_SUM(Type, Container) \
305  virtual SGVector<Type> rowwise_sum(const linalg::Block<Container<Type>>& a, bool no_diag) const \
306  { \
307  return rowwise_sum_impl(a, no_diag); \
308  }
310  #undef BACKEND_GENERIC_BLOCK_ROWWISE_SUM
311 
312  #undef DEFINE_FOR_ALL_PTYPE
313 
314 private:
316  template <typename T>
317  void add_impl(SGVector<T>& a, SGVector<T>& b, T alpha, T beta, SGVector<T>& result) const
318  {
319  typename SGVector<T>::EigenVectorXtMap a_eig = a;
320  typename SGVector<T>::EigenVectorXtMap b_eig = b;
321  typename SGVector<T>::EigenVectorXtMap result_eig = result;
322 
323  result_eig = alpha * a_eig + beta * b_eig;
324  }
325 
327  template <typename T>
328  void add_impl(SGMatrix<T>& a, SGMatrix<T>& b, T alpha, T beta, SGMatrix<T>& result) const
329  {
330  typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
331  typename SGMatrix<T>::EigenMatrixXtMap b_eig = b;
332  typename SGMatrix<T>::EigenMatrixXtMap result_eig = result;
333 
334  result_eig = alpha * a_eig + beta * b_eig;
335  }
336 
338  template <typename T>
339  SGMatrix<T> cholesky_factor_impl(const SGMatrix<T>& A, const bool lower) const
340  {
341  SGMatrix<T> c(A.num_rows, A.num_cols);
342  set_const_impl<T>(c, 0);
343  typename SGMatrix<T>::EigenMatrixXtMap A_eig = A;
344  typename SGMatrix<T>::EigenMatrixXtMap c_eig = c;
345 
346  Eigen::LLT<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> > llt(A_eig);
347 
348  //compute matrix L or U
349  if(lower==false)
350  c_eig = llt.matrixU();
351  else
352  c_eig = llt.matrixL();
353 
354  /*
355  * checking for success
356  *
357  * 0: Eigen::Success. Decomposition was successful
358  * 1: Eigen::NumericalIssue. The provided data did not satisfy the prerequisites.
359  */
360  REQUIRE(llt.info()!=Eigen::NumericalIssue, "Matrix is not Hermitian positive definite!\n");
361 
362  return c;
363  }
364 
366  template <typename T>
367  SGVector<T> cholesky_solver_impl(const SGMatrix<T>& L, const SGVector<T>& b,
368  const bool lower) const
369  {
370  SGVector<T> x(b.size());
371  set_const_impl<T>(x, 0);
372  typename SGMatrix<T>::EigenMatrixXtMap L_eig = L;
373  typename SGVector<T>::EigenVectorXtMap b_eig = b;
374  typename SGVector<T>::EigenVectorXtMap x_eig = x;
375 
376  if (lower == false)
377  {
378  Eigen::TriangularView<Eigen::Map<typename SGMatrix<T>::EigenMatrixXt,
379  0, Eigen::Stride<0,0> >, Eigen::Upper> tlv(L_eig);
380 
381  x_eig = (tlv.transpose()).solve(tlv.solve(b_eig));
382  }
383  else
384  {
385  Eigen::TriangularView<Eigen::Map<typename SGMatrix<T>::EigenMatrixXt,
386  0, Eigen::Stride<0,0> >, Eigen::Lower> tlv(L_eig);
387  x_eig = (tlv.transpose()).solve(tlv.solve(b_eig));
388  }
389 
390  return x;
391  }
392 
394  template <typename T>
395  T dot_impl(const SGVector<T>& a, const SGVector<T>& b) const
396  {
397  return (typename SGVector<T>::EigenVectorXtMap(a)).dot(typename SGVector<T>::EigenVectorXtMap(b));
398  }
399 
401  template <typename T>
402  void element_prod_impl(SGMatrix<T>& a, SGMatrix<T>& b, SGMatrix<T>& result) const
403  {
404  typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
405  typename SGMatrix<T>::EigenMatrixXtMap b_eig = b;
406  typename SGMatrix<T>::EigenMatrixXtMap result_eig = result;
407 
408  result_eig = a_eig.array() * b_eig.array();
409  }
410 
412  template <typename T>
413  void logistic_impl(SGMatrix<T>& a, SGMatrix<T>& result) const
414  {
415  typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
416  typename SGMatrix<T>::EigenMatrixXtMap result_eig = result;
417 
418  result_eig = (T)1 / (1 + ((-1 * a_eig).array()).exp());
419  }
420 
422  template <typename T>
423  void element_prod_impl(linalg::Block<SGMatrix<T>>& a,
424  linalg::Block<SGMatrix<T>>& b, SGMatrix<T>& result) const
425  {
426  typename SGMatrix<T>::EigenMatrixXtMap a_eig = a.m_matrix;
427  typename SGMatrix<T>::EigenMatrixXtMap b_eig = b.m_matrix;
428  typename SGMatrix<T>::EigenMatrixXtMap result_eig = result;
429 
430  Eigen::Block<typename SGMatrix<T>::EigenMatrixXtMap> a_block =
431  a_eig.block(a.m_row_begin, a.m_col_begin, a.m_row_size, a.m_col_size);
432  Eigen::Block<typename SGMatrix<T>::EigenMatrixXtMap> b_block =
433  b_eig.block(b.m_row_begin, b.m_col_begin, b.m_row_size, b.m_col_size);
434 
435  result_eig = a_block.array() * b_block.array();
436  }
437 
439  template <typename T>
440  void matrix_prod_impl(SGMatrix<T>& a, SGVector<T>& b, SGVector<T>& result,
441  bool transpose, bool transpose_B=false) const
442  {
443  typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
444  typename SGVector<T>::EigenVectorXtMap b_eig = b;
445  typename SGVector<T>::EigenVectorXtMap result_eig = result;
446 
447  if (transpose)
448  result_eig = a_eig.transpose() * b_eig;
449  else
450  result_eig = a_eig * b_eig;
451  }
452 
454  template <typename T>
455  void matrix_prod_impl(SGMatrix<T>& a, SGMatrix<T>& b, SGMatrix<T>& result,
456  bool transpose_A, bool transpose_B) const
457  {
458  typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
459  typename SGMatrix<T>::EigenMatrixXtMap b_eig = b;
460  typename SGMatrix<T>::EigenMatrixXtMap result_eig = result;
461 
462  if (transpose_A && transpose_B)
463  result_eig = a_eig.transpose() * b_eig.transpose();
464 
465  else if (transpose_A)
466  result_eig = a_eig.transpose() * b_eig;
467 
468  else if (transpose_B)
469  result_eig = a_eig * b_eig.transpose();
470 
471  else
472  result_eig = a_eig * b_eig;
473  }
474 
476  template <typename T>
477  T max_impl(const SGVector<T>& vec) const
478  {
479  return (typename SGVector<T>::EigenVectorXtMap(vec)).maxCoeff();
480  }
481 
483  template <typename T>
484  T max_impl(const SGMatrix<T>& mat) const
485  {
486  return (typename SGMatrix<T>::EigenMatrixXtMap(mat)).maxCoeff();
487  }
488 
490  template <typename T, template <typename> class Container>
491  typename std::enable_if<!std::is_same<T, complex128_t>::value, float64_t>::type
492  mean_impl(const Container<T>& a) const
493  {
494  return sum_impl(a)/(float64_t(a.size()));
495  }
496 
498  template<template <typename> class Container>
499  complex128_t mean_impl(const Container<complex128_t>& a) const
500  {
501  return sum_impl(a)/(complex128_t(a.size()));
502  }
503 
505  template <typename T, template <typename> class Container>
506  void range_fill_impl(Container<T>& a, const T start) const
507  {
508  for (index_t i = 0; i < a.size(); ++i)
509  a[i] = start + T(i);
510  }
511 
513  template <typename T>
514  void scale_impl(SGVector<T>& a, T alpha, SGVector<T>& result) const
515  {
516  typename SGVector<T>::EigenVectorXtMap a_eig = a;
517  typename SGVector<T>::EigenVectorXtMap result_eig = result;
518 
519  result_eig = alpha * a_eig;
520  }
521 
523  template <typename T>
524  void scale_impl(SGMatrix<T>& a, T alpha, SGMatrix<T>& result) const
525  {
526  typename SGMatrix<T>::EigenMatrixXtMap a_eig = a;
527  typename SGMatrix<T>::EigenMatrixXtMap result_eig = result;
528 
529  result_eig = alpha * a_eig;
530  }
531 
533  template <typename T, template <typename> class Container>
534  void set_const_impl(Container<T>& a, T value) const
535  {
536  for (index_t i = 0; i < a.size(); ++i)
537  a[i] = value;
538  }
539 
541  template <typename T>
542  T sum_impl(const SGVector<T>& vec, bool no_diag=false) const
543  {
544  return (typename SGVector<T>::EigenVectorXtMap(vec)).sum();
545  }
546 
548  template <typename T>
549  T sum_impl(const SGMatrix<T>& mat, bool no_diag=false) const
550  {
551  typename SGMatrix<T>::EigenMatrixXtMap m = mat;
552  T result = m.sum();
553  if (no_diag)
554  result -= m.diagonal().sum();
555 
556  return result;
557  }
558 
560  template <typename T>
561  T sum_impl(const linalg::Block<SGMatrix<T>>& mat, bool no_diag=false) const
562  {
563  typename SGMatrix<T>::EigenMatrixXtMap m = mat.m_matrix;
564  Eigen::Block<typename SGMatrix<T>::EigenMatrixXtMap> m_block = m.block(
565  mat.m_row_begin, mat.m_col_begin, mat.m_row_size, mat.m_col_size);
566 
567  T result = m_block.sum();
568  if (no_diag)
569  result -= m_block.diagonal().sum();
570 
571  return result;
572  }
573 
575  template <typename T>
576  T sum_symmetric_impl(const SGMatrix<T>& mat, bool no_diag=false) const
577  {
578  typename SGMatrix<T>::EigenMatrixXtMap m = mat;
579 
580  // since the matrix is symmetric with main diagonal inside, we can save half
581  // the computation with using only the upper triangular part.
582  typename SGMatrix<T>::EigenMatrixXt m_upper =
583  m.template triangularView<Eigen::StrictlyUpper>();
584  T result = m_upper.sum();
585  result += result;
586 
587  if (!no_diag)
588  result += m.diagonal().sum();
589  return result;
590  }
591 
593  template <typename T>
594  T sum_symmetric_impl(const linalg::Block<SGMatrix<T>>& mat, bool no_diag=false) const
595  {
596  typename SGMatrix<T>::EigenMatrixXtMap m = mat.m_matrix;
597  Eigen::Block<typename SGMatrix<T>::EigenMatrixXtMap> m_block = m.block(
598  mat.m_row_begin, mat.m_col_begin, mat.m_row_size, mat.m_col_size);
599 
600  // since the matrix is symmetric with main diagonal inside, we can save half
601  // the computation with using only the upper triangular part.
602  typename SGMatrix<T>::EigenMatrixXt m_upper =
603  m_block.template triangularView<Eigen::StrictlyUpper>();
604  T result = m_upper.sum();
605  result += result;
606 
607  if (!no_diag)
608  result += m_block.diagonal().sum();
609  return result;
610  }
611 
613  template <typename T>
614  SGVector<T> colwise_sum_impl(const SGMatrix<T>& mat, bool no_diag) const
615  {
616  SGVector<T> result(mat.num_cols);
617 
618  typename SGMatrix<T>::EigenMatrixXtMap mat_eig = mat;
619  typename SGVector<T>::EigenVectorXtMap result_eig = result;
620 
621  result_eig = mat_eig.colwise().sum();
622 
623  // remove the main diagonal elements if required
624  if (no_diag)
625  {
626  index_t len_major_diag = mat_eig.rows() < mat_eig.cols()
627  ? mat_eig.rows() : mat_eig.cols();
628  for (index_t i = 0; i < len_major_diag; ++i)
629  result_eig[i] -= mat_eig(i,i);
630  }
631 
632  return result;
633  }
634 
636  template <typename T>
637  SGVector<T> colwise_sum_impl(const linalg::Block<SGMatrix<T>>& mat, bool no_diag) const
638  {
639  SGVector<T> result(mat.m_col_size);
640 
641  typename SGMatrix<T>::EigenMatrixXtMap m = mat.m_matrix;
642  Eigen::Block<typename SGMatrix<T>::EigenMatrixXtMap> m_block = m.block(
643  mat.m_row_begin, mat.m_col_begin, mat.m_row_size, mat.m_col_size);
644  typename SGVector<T>::EigenVectorXtMap result_eig = result;
645 
646  result_eig = m_block.colwise().sum();
647 
648  // remove the main diagonal elements if required
649  if (no_diag)
650  {
651  index_t len_major_diag = m_block.rows() < m_block.cols()
652  ? m_block.rows() : m_block.cols();
653  for (index_t i = 0; i < len_major_diag; ++i)
654  result_eig[i] -= m_block(i,i);
655  }
656 
657  return result;
658  }
659 
661  template <typename T>
662  SGVector<T> rowwise_sum_impl(const SGMatrix<T>& mat, bool no_diag) const
663  {
664  SGVector<T> result(mat.num_rows);
665 
666  typename SGMatrix<T>::EigenMatrixXtMap mat_eig = mat;
667  typename SGVector<T>::EigenVectorXtMap result_eig = result;
668 
669  result_eig = mat_eig.rowwise().sum();
670 
671  // remove the main diagonal elements if required
672  if (no_diag)
673  {
674  index_t len_major_diag = mat_eig.rows() < mat_eig.cols()
675  ? mat_eig.rows() : mat_eig.cols();
676  for (index_t i = 0; i < len_major_diag; ++i)
677  result_eig[i] -= mat_eig(i,i);
678  }
679 
680  return result;
681  }
682 
684  template <typename T>
685  SGVector<T> rowwise_sum_impl(const linalg::Block<SGMatrix<T>>& mat, bool no_diag) const
686  {
687  SGVector<T> result(mat.m_row_size);
688 
689  typename SGMatrix<T>::EigenMatrixXtMap m = mat.m_matrix;
690  Eigen::Block<typename SGMatrix<T>::EigenMatrixXtMap> m_block = m.block(
691  mat.m_row_begin, mat.m_col_begin, mat.m_row_size, mat.m_col_size);
692  typename SGVector<T>::EigenVectorXtMap result_eig = result;
693 
694  result_eig = m_block.rowwise().sum();
695 
696  // remove the main diagonal elements if required
697  if (no_diag)
698  {
699  index_t len_major_diag = m_block.rows() < m_block.cols()
700  ? m_block.rows() : m_block.cols();
701  for (index_t i = 0; i < len_major_diag; ++i)
702  result_eig[i] -= m_block(i,i);
703  }
704 
705  return result;
706  }
707 
708 #undef DEFINE_FOR_ALL_PTYPE
709 #undef DEFINE_FOR_REAL_PTYPE
710 #undef DEFINE_FOR_NON_INTEGER_PTYPE
711 #undef DEFINE_FOR_NUMERIC_PTYPE
712 };
713 
714 }
715 
716 #endif //LINALG_BACKEND_EIGEN_H__
#define BACKEND_GENERIC_BLOCK_COLWISE_SUM(Type, Container)
#define BACKEND_GENERIC_CHOLESKY_FACTOR(Type, Container)
std::complex< float64_t > complex128_t
Definition: common.h:77
#define BACKEND_GENERIC_COMPLEX_MEAN(Container)
#define DEFINE_FOR_ALL_PTYPE(METHODNAME, Container)
Eigen::Map< EigenMatrixXt, 0, Eigen::Stride< 0, 0 > > EigenMatrixXtMap
Definition: SGMatrix.h:44
#define BACKEND_GENERIC_ROWWISE_SUM(Type, Container)
int32_t index_t
Definition: common.h:72
#define BACKEND_GENERIC_IN_PLACE_BLOCK_ELEMENT_PROD(Type, Container)
Eigen::Matrix< T,-1,-1, 0,-1,-1 > EigenMatrixXt
Definition: SGMatrix.h:43
#define REQUIRE(x,...)
Definition: SGIO.h:205
#define BACKEND_GENERIC_DOT(Type, Container)
Eigen::Map< EigenVectorXt, 0, Eigen::Stride< 0, 0 > > EigenVectorXtMap
Definition: SGVector.h:50
#define BACKEND_GENERIC_BLOCK_SUM(Type, Container)
#define BACKEND_GENERIC_REAL_MEAN(Type, Container)
#define BACKEND_GENERIC_COLWISE_SUM(Type, Container)
#define BACKEND_GENERIC_BLOCK_ROWWISE_SUM(Type, Container)
#define BACKEND_GENERIC_MAX(Type, Container)
#define BACKEND_GENERIC_IN_PLACE_MATRIX_PROD(Type, Container)
#define BACKEND_GENERIC_SUM(Type, Container)
#define BACKEND_GENERIC_CHOLESKY_SOLVER(Type, Container)
#define BACKEND_GENERIC_IN_PLACE_ELEMENT_PROD(Type, Container)
double float64_t
Definition: common.h:60
#define DEFINE_FOR_REAL_PTYPE(METHODNAME, Container)
#define BACKEND_GENERIC_LOGISTIC(Type, Container)
#define BACKEND_GENERIC_IN_PLACE_SCALE(Type, Container)
#define BACKEND_GENERIC_SYMMETRIC_SUM(Type, Container)
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
Base interface of generic linalg methods and generic memory transfer methods.
#define BACKEND_GENERIC_IN_PLACE_ADD(Type, Container)
Linalg methods with Eigen3 backend.
#define DEFINE_FOR_NON_INTEGER_PTYPE(METHODNAME, Container)
#define BACKEND_GENERIC_SYMMETRIC_BLOCK_SUM(Type, Container)
#define BACKEND_GENERIC_SET_CONST(Type, Container)
#define DEFINE_FOR_NUMERIC_PTYPE(METHODNAME, Container)
#define BACKEND_GENERIC_RANGE_FILL(Type, Container)

SHOGUN Machine Learning Toolbox - Documentation