vectorwiseop.cpp 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298
  1. // This file is part of Eigen, a lightweight C++ template library
  2. // for linear algebra.
  3. //
  4. // Copyright (C) 2011 Benoit Jacob <jacob.benoit.1@gmail.com>
  5. // Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr>
  6. //
  7. // This Source Code Form is subject to the terms of the Mozilla
  8. // Public License v. 2.0. If a copy of the MPL was not distributed
  9. // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
  10. #define TEST_ENABLE_TEMPORARY_TRACKING
  11. #define EIGEN_NO_STATIC_ASSERT
  12. #include "main.h"
  13. template<typename ArrayType> void vectorwiseop_array(const ArrayType& m)
  14. {
  15. typedef typename ArrayType::Scalar Scalar;
  16. typedef Array<Scalar, ArrayType::RowsAtCompileTime, 1> ColVectorType;
  17. typedef Array<Scalar, 1, ArrayType::ColsAtCompileTime> RowVectorType;
  18. Index rows = m.rows();
  19. Index cols = m.cols();
  20. Index r = internal::random<Index>(0, rows-1),
  21. c = internal::random<Index>(0, cols-1);
  22. ArrayType m1 = ArrayType::Random(rows, cols),
  23. m2(rows, cols),
  24. m3(rows, cols);
  25. ColVectorType colvec = ColVectorType::Random(rows);
  26. RowVectorType rowvec = RowVectorType::Random(cols);
  27. // test addition
  28. m2 = m1;
  29. m2.colwise() += colvec;
  30. VERIFY_IS_APPROX(m2, m1.colwise() + colvec);
  31. VERIFY_IS_APPROX(m2.col(c), m1.col(c) + colvec);
  32. VERIFY_RAISES_ASSERT(m2.colwise() += colvec.transpose());
  33. VERIFY_RAISES_ASSERT(m1.colwise() + colvec.transpose());
  34. m2 = m1;
  35. m2.rowwise() += rowvec;
  36. VERIFY_IS_APPROX(m2, m1.rowwise() + rowvec);
  37. VERIFY_IS_APPROX(m2.row(r), m1.row(r) + rowvec);
  38. VERIFY_RAISES_ASSERT(m2.rowwise() += rowvec.transpose());
  39. VERIFY_RAISES_ASSERT(m1.rowwise() + rowvec.transpose());
  40. // test substraction
  41. m2 = m1;
  42. m2.colwise() -= colvec;
  43. VERIFY_IS_APPROX(m2, m1.colwise() - colvec);
  44. VERIFY_IS_APPROX(m2.col(c), m1.col(c) - colvec);
  45. VERIFY_RAISES_ASSERT(m2.colwise() -= colvec.transpose());
  46. VERIFY_RAISES_ASSERT(m1.colwise() - colvec.transpose());
  47. m2 = m1;
  48. m2.rowwise() -= rowvec;
  49. VERIFY_IS_APPROX(m2, m1.rowwise() - rowvec);
  50. VERIFY_IS_APPROX(m2.row(r), m1.row(r) - rowvec);
  51. VERIFY_RAISES_ASSERT(m2.rowwise() -= rowvec.transpose());
  52. VERIFY_RAISES_ASSERT(m1.rowwise() - rowvec.transpose());
  53. // test multiplication
  54. m2 = m1;
  55. m2.colwise() *= colvec;
  56. VERIFY_IS_APPROX(m2, m1.colwise() * colvec);
  57. VERIFY_IS_APPROX(m2.col(c), m1.col(c) * colvec);
  58. VERIFY_RAISES_ASSERT(m2.colwise() *= colvec.transpose());
  59. VERIFY_RAISES_ASSERT(m1.colwise() * colvec.transpose());
  60. m2 = m1;
  61. m2.rowwise() *= rowvec;
  62. VERIFY_IS_APPROX(m2, m1.rowwise() * rowvec);
  63. VERIFY_IS_APPROX(m2.row(r), m1.row(r) * rowvec);
  64. VERIFY_RAISES_ASSERT(m2.rowwise() *= rowvec.transpose());
  65. VERIFY_RAISES_ASSERT(m1.rowwise() * rowvec.transpose());
  66. // test quotient
  67. m2 = m1;
  68. m2.colwise() /= colvec;
  69. VERIFY_IS_APPROX(m2, m1.colwise() / colvec);
  70. VERIFY_IS_APPROX(m2.col(c), m1.col(c) / colvec);
  71. VERIFY_RAISES_ASSERT(m2.colwise() /= colvec.transpose());
  72. VERIFY_RAISES_ASSERT(m1.colwise() / colvec.transpose());
  73. m2 = m1;
  74. m2.rowwise() /= rowvec;
  75. VERIFY_IS_APPROX(m2, m1.rowwise() / rowvec);
  76. VERIFY_IS_APPROX(m2.row(r), m1.row(r) / rowvec);
  77. VERIFY_RAISES_ASSERT(m2.rowwise() /= rowvec.transpose());
  78. VERIFY_RAISES_ASSERT(m1.rowwise() / rowvec.transpose());
  79. m2 = m1;
  80. // yes, there might be an aliasing issue there but ".rowwise() /="
  81. // is supposed to evaluate " m2.colwise().sum()" into a temporary to avoid
  82. // evaluating the reduction multiple times
  83. if(ArrayType::RowsAtCompileTime>2 || ArrayType::RowsAtCompileTime==Dynamic)
  84. {
  85. m2.rowwise() /= m2.colwise().sum();
  86. VERIFY_IS_APPROX(m2, m1.rowwise() / m1.colwise().sum());
  87. }
  88. // all/any
  89. Array<bool,Dynamic,Dynamic> mb(rows,cols);
  90. mb = (m1.real()<=0.7).colwise().all();
  91. VERIFY( (mb.col(c) == (m1.real().col(c)<=0.7).all()).all() );
  92. mb = (m1.real()<=0.7).rowwise().all();
  93. VERIFY( (mb.row(r) == (m1.real().row(r)<=0.7).all()).all() );
  94. mb = (m1.real()>=0.7).colwise().any();
  95. VERIFY( (mb.col(c) == (m1.real().col(c)>=0.7).any()).all() );
  96. mb = (m1.real()>=0.7).rowwise().any();
  97. VERIFY( (mb.row(r) == (m1.real().row(r)>=0.7).any()).all() );
  98. }
  99. template<typename MatrixType> void vectorwiseop_matrix(const MatrixType& m)
  100. {
  101. typedef typename MatrixType::Scalar Scalar;
  102. typedef typename NumTraits<Scalar>::Real RealScalar;
  103. typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType;
  104. typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
  105. typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealColVectorType;
  106. typedef Matrix<RealScalar, 1, MatrixType::ColsAtCompileTime> RealRowVectorType;
  107. typedef Matrix<Scalar,Dynamic,Dynamic> MatrixX;
  108. Index rows = m.rows();
  109. Index cols = m.cols();
  110. Index r = internal::random<Index>(0, rows-1),
  111. c = internal::random<Index>(0, cols-1);
  112. MatrixType m1 = MatrixType::Random(rows, cols),
  113. m2(rows, cols),
  114. m3(rows, cols);
  115. ColVectorType colvec = ColVectorType::Random(rows);
  116. RowVectorType rowvec = RowVectorType::Random(cols);
  117. RealColVectorType rcres;
  118. RealRowVectorType rrres;
  119. // test broadcast assignment
  120. m2 = m1;
  121. m2.colwise() = colvec;
  122. for(Index j=0; j<cols; ++j)
  123. VERIFY_IS_APPROX(m2.col(j), colvec);
  124. m2.rowwise() = rowvec;
  125. for(Index i=0; i<rows; ++i)
  126. VERIFY_IS_APPROX(m2.row(i), rowvec);
  127. if(rows>1)
  128. VERIFY_RAISES_ASSERT(m2.colwise() = colvec.transpose());
  129. if(cols>1)
  130. VERIFY_RAISES_ASSERT(m2.rowwise() = rowvec.transpose());
  131. // test addition
  132. m2 = m1;
  133. m2.colwise() += colvec;
  134. VERIFY_IS_APPROX(m2, m1.colwise() + colvec);
  135. VERIFY_IS_APPROX(m2.col(c), m1.col(c) + colvec);
  136. if(rows>1)
  137. {
  138. VERIFY_RAISES_ASSERT(m2.colwise() += colvec.transpose());
  139. VERIFY_RAISES_ASSERT(m1.colwise() + colvec.transpose());
  140. }
  141. m2 = m1;
  142. m2.rowwise() += rowvec;
  143. VERIFY_IS_APPROX(m2, m1.rowwise() + rowvec);
  144. VERIFY_IS_APPROX(m2.row(r), m1.row(r) + rowvec);
  145. if(cols>1)
  146. {
  147. VERIFY_RAISES_ASSERT(m2.rowwise() += rowvec.transpose());
  148. VERIFY_RAISES_ASSERT(m1.rowwise() + rowvec.transpose());
  149. }
  150. // test substraction
  151. m2 = m1;
  152. m2.colwise() -= colvec;
  153. VERIFY_IS_APPROX(m2, m1.colwise() - colvec);
  154. VERIFY_IS_APPROX(m2.col(c), m1.col(c) - colvec);
  155. if(rows>1)
  156. {
  157. VERIFY_RAISES_ASSERT(m2.colwise() -= colvec.transpose());
  158. VERIFY_RAISES_ASSERT(m1.colwise() - colvec.transpose());
  159. }
  160. m2 = m1;
  161. m2.rowwise() -= rowvec;
  162. VERIFY_IS_APPROX(m2, m1.rowwise() - rowvec);
  163. VERIFY_IS_APPROX(m2.row(r), m1.row(r) - rowvec);
  164. if(cols>1)
  165. {
  166. VERIFY_RAISES_ASSERT(m2.rowwise() -= rowvec.transpose());
  167. VERIFY_RAISES_ASSERT(m1.rowwise() - rowvec.transpose());
  168. }
  169. // ------ partial reductions ------
  170. #define TEST_PARTIAL_REDUX_BASIC(FUNC,ROW,COL,PREPROCESS) { \
  171. ROW = m1 PREPROCESS .colwise().FUNC ; \
  172. for(Index k=0; k<cols; ++k) VERIFY_IS_APPROX(ROW(k), m1.col(k) PREPROCESS .FUNC ); \
  173. COL = m1 PREPROCESS .rowwise().FUNC ; \
  174. for(Index k=0; k<rows; ++k) VERIFY_IS_APPROX(COL(k), m1.row(k) PREPROCESS .FUNC ); \
  175. }
  176. TEST_PARTIAL_REDUX_BASIC(sum(), rowvec,colvec,EIGEN_EMPTY);
  177. TEST_PARTIAL_REDUX_BASIC(prod(), rowvec,colvec,EIGEN_EMPTY);
  178. TEST_PARTIAL_REDUX_BASIC(mean(), rowvec,colvec,EIGEN_EMPTY);
  179. TEST_PARTIAL_REDUX_BASIC(minCoeff(), rrres, rcres, .real());
  180. TEST_PARTIAL_REDUX_BASIC(maxCoeff(), rrres, rcres, .real());
  181. TEST_PARTIAL_REDUX_BASIC(norm(), rrres, rcres, EIGEN_EMPTY);
  182. TEST_PARTIAL_REDUX_BASIC(squaredNorm(),rrres, rcres, EIGEN_EMPTY);
  183. TEST_PARTIAL_REDUX_BASIC(redux(internal::scalar_sum_op<Scalar,Scalar>()),rowvec,colvec,EIGEN_EMPTY);
  184. VERIFY_IS_APPROX(m1.cwiseAbs().colwise().sum(), m1.colwise().template lpNorm<1>());
  185. VERIFY_IS_APPROX(m1.cwiseAbs().rowwise().sum(), m1.rowwise().template lpNorm<1>());
  186. VERIFY_IS_APPROX(m1.cwiseAbs().colwise().maxCoeff(), m1.colwise().template lpNorm<Infinity>());
  187. VERIFY_IS_APPROX(m1.cwiseAbs().rowwise().maxCoeff(), m1.rowwise().template lpNorm<Infinity>());
  188. // regression for bug 1158
  189. VERIFY_IS_APPROX(m1.cwiseAbs().colwise().sum().x(), m1.col(0).cwiseAbs().sum());
  190. // test normalized
  191. m2 = m1.colwise().normalized();
  192. VERIFY_IS_APPROX(m2.col(c), m1.col(c).normalized());
  193. m2 = m1.rowwise().normalized();
  194. VERIFY_IS_APPROX(m2.row(r), m1.row(r).normalized());
  195. // test normalize
  196. m2 = m1;
  197. m2.colwise().normalize();
  198. VERIFY_IS_APPROX(m2.col(c), m1.col(c).normalized());
  199. m2 = m1;
  200. m2.rowwise().normalize();
  201. VERIFY_IS_APPROX(m2.row(r), m1.row(r).normalized());
  202. // test with partial reduction of products
  203. Matrix<Scalar,MatrixType::RowsAtCompileTime,MatrixType::RowsAtCompileTime> m1m1 = m1 * m1.transpose();
  204. VERIFY_IS_APPROX( (m1 * m1.transpose()).colwise().sum(), m1m1.colwise().sum());
  205. Matrix<Scalar,1,MatrixType::RowsAtCompileTime> tmp(rows);
  206. VERIFY_EVALUATION_COUNT( tmp = (m1 * m1.transpose()).colwise().sum(), 1);
  207. m2 = m1.rowwise() - (m1.colwise().sum()/RealScalar(m1.rows())).eval();
  208. m1 = m1.rowwise() - (m1.colwise().sum()/RealScalar(m1.rows()));
  209. VERIFY_IS_APPROX( m1, m2 );
  210. VERIFY_EVALUATION_COUNT( m2 = (m1.rowwise() - m1.colwise().sum()/RealScalar(m1.rows())), (MatrixType::RowsAtCompileTime!=1 ? 1 : 0) );
  211. // test empty expressions
  212. VERIFY_IS_APPROX(m1.matrix().middleCols(0,0).rowwise().sum().eval(), MatrixX::Zero(rows,1));
  213. VERIFY_IS_APPROX(m1.matrix().middleRows(0,0).colwise().sum().eval(), MatrixX::Zero(1,cols));
  214. VERIFY_IS_APPROX(m1.matrix().middleCols(0,fix<0>).rowwise().sum().eval(), MatrixX::Zero(rows,1));
  215. VERIFY_IS_APPROX(m1.matrix().middleRows(0,fix<0>).colwise().sum().eval(), MatrixX::Zero(1,cols));
  216. VERIFY_IS_APPROX(m1.matrix().middleCols(0,0).rowwise().prod().eval(), MatrixX::Ones(rows,1));
  217. VERIFY_IS_APPROX(m1.matrix().middleRows(0,0).colwise().prod().eval(), MatrixX::Ones(1,cols));
  218. VERIFY_IS_APPROX(m1.matrix().middleCols(0,fix<0>).rowwise().prod().eval(), MatrixX::Ones(rows,1));
  219. VERIFY_IS_APPROX(m1.matrix().middleRows(0,fix<0>).colwise().prod().eval(), MatrixX::Ones(1,cols));
  220. VERIFY_IS_APPROX(m1.matrix().middleCols(0,0).rowwise().squaredNorm().eval(), MatrixX::Zero(rows,1));
  221. VERIFY_RAISES_ASSERT(m1.real().middleCols(0,0).rowwise().minCoeff().eval());
  222. VERIFY_RAISES_ASSERT(m1.real().middleRows(0,0).colwise().maxCoeff().eval());
  223. VERIFY_IS_EQUAL(m1.real().middleRows(0,0).rowwise().maxCoeff().eval().rows(),0);
  224. VERIFY_IS_EQUAL(m1.real().middleCols(0,0).colwise().maxCoeff().eval().cols(),0);
  225. VERIFY_IS_EQUAL(m1.real().middleRows(0,fix<0>).rowwise().maxCoeff().eval().rows(),0);
  226. VERIFY_IS_EQUAL(m1.real().middleCols(0,fix<0>).colwise().maxCoeff().eval().cols(),0);
  227. }
  228. EIGEN_DECLARE_TEST(vectorwiseop)
  229. {
  230. CALL_SUBTEST_1( vectorwiseop_array(Array22cd()) );
  231. CALL_SUBTEST_2( vectorwiseop_array(Array<double, 3, 2>()) );
  232. CALL_SUBTEST_3( vectorwiseop_array(ArrayXXf(3, 4)) );
  233. CALL_SUBTEST_4( vectorwiseop_matrix(Matrix4cf()) );
  234. CALL_SUBTEST_5( vectorwiseop_matrix(Matrix4f()) );
  235. CALL_SUBTEST_5( vectorwiseop_matrix(Vector4f()) );
  236. CALL_SUBTEST_5( vectorwiseop_matrix(Matrix<float,4,5>()) );
  237. CALL_SUBTEST_6( vectorwiseop_matrix(MatrixXd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
  238. CALL_SUBTEST_7( vectorwiseop_matrix(VectorXd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
  239. CALL_SUBTEST_7( vectorwiseop_matrix(RowVectorXd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
  240. }