vectorwiseop.cpp
Go to the documentation of this file.
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 
11 #define TEST_ENABLE_TEMPORARY_TRACKING
12 #define EIGEN_NO_STATIC_ASSERT
13 
14 #include "main.h"
15 
16 template<typename ArrayType> void vectorwiseop_array(const ArrayType& m)
17 {
18  typedef typename ArrayType::Scalar Scalar;
19  typedef Array<Scalar, ArrayType::RowsAtCompileTime, 1> ColVectorType;
20  typedef Array<Scalar, 1, ArrayType::ColsAtCompileTime> RowVectorType;
21 
22  Index rows = m.rows();
23  Index cols = m.cols();
24  Index r = internal::random<Index>(0, rows-1),
25  c = internal::random<Index>(0, cols-1);
26 
27  ArrayType m1 = ArrayType::Random(rows, cols),
28  m2(rows, cols),
29  m3(rows, cols);
30 
31  ColVectorType colvec = ColVectorType::Random(rows);
32  RowVectorType rowvec = RowVectorType::Random(cols);
33 
34  // test addition
35 
36  m2 = m1;
37  m2.colwise() += colvec;
38  VERIFY_IS_APPROX(m2, m1.colwise() + colvec);
39  VERIFY_IS_APPROX(m2.col(c), m1.col(c) + colvec);
40 
41  VERIFY_RAISES_ASSERT(m2.colwise() += colvec.transpose());
42  VERIFY_RAISES_ASSERT(m1.colwise() + colvec.transpose());
43 
44  m2 = m1;
45  m2.rowwise() += rowvec;
46  VERIFY_IS_APPROX(m2, m1.rowwise() + rowvec);
47  VERIFY_IS_APPROX(m2.row(r), m1.row(r) + rowvec);
48 
49  VERIFY_RAISES_ASSERT(m2.rowwise() += rowvec.transpose());
50  VERIFY_RAISES_ASSERT(m1.rowwise() + rowvec.transpose());
51 
52  // test substraction
53 
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 
59  VERIFY_RAISES_ASSERT(m2.colwise() -= colvec.transpose());
60  VERIFY_RAISES_ASSERT(m1.colwise() - colvec.transpose());
61 
62  m2 = m1;
63  m2.rowwise() -= rowvec;
64  VERIFY_IS_APPROX(m2, m1.rowwise() - rowvec);
65  VERIFY_IS_APPROX(m2.row(r), m1.row(r) - rowvec);
66 
67  VERIFY_RAISES_ASSERT(m2.rowwise() -= rowvec.transpose());
68  VERIFY_RAISES_ASSERT(m1.rowwise() - rowvec.transpose());
69 
70  // test multiplication
71 
72  m2 = m1;
73  m2.colwise() *= colvec;
74  VERIFY_IS_APPROX(m2, m1.colwise() * colvec);
75  VERIFY_IS_APPROX(m2.col(c), m1.col(c) * colvec);
76 
77  VERIFY_RAISES_ASSERT(m2.colwise() *= colvec.transpose());
78  VERIFY_RAISES_ASSERT(m1.colwise() * colvec.transpose());
79 
80  m2 = m1;
81  m2.rowwise() *= rowvec;
82  VERIFY_IS_APPROX(m2, m1.rowwise() * rowvec);
83  VERIFY_IS_APPROX(m2.row(r), m1.row(r) * rowvec);
84 
85  VERIFY_RAISES_ASSERT(m2.rowwise() *= rowvec.transpose());
86  VERIFY_RAISES_ASSERT(m1.rowwise() * rowvec.transpose());
87 
88  // test quotient
89 
90  m2 = m1;
91  m2.colwise() /= colvec;
92  VERIFY_IS_APPROX(m2, m1.colwise() / colvec);
93  VERIFY_IS_APPROX(m2.col(c), m1.col(c) / colvec);
94 
95  VERIFY_RAISES_ASSERT(m2.colwise() /= colvec.transpose());
96  VERIFY_RAISES_ASSERT(m1.colwise() / colvec.transpose());
97 
98  m2 = m1;
99  m2.rowwise() /= rowvec;
100  VERIFY_IS_APPROX(m2, m1.rowwise() / rowvec);
101  VERIFY_IS_APPROX(m2.row(r), m1.row(r) / rowvec);
102 
103  VERIFY_RAISES_ASSERT(m2.rowwise() /= rowvec.transpose());
104  VERIFY_RAISES_ASSERT(m1.rowwise() / rowvec.transpose());
105 
106  m2 = m1;
107  // yes, there might be an aliasing issue there but ".rowwise() /="
108  // is supposed to evaluate " m2.colwise().sum()" into a temporary to avoid
109  // evaluating the reduction multiple times
110  if(ArrayType::RowsAtCompileTime>2 || ArrayType::RowsAtCompileTime==Dynamic)
111  {
112  m2.rowwise() /= m2.colwise().sum();
113  VERIFY_IS_APPROX(m2, m1.rowwise() / m1.colwise().sum());
114  }
115 
116  // all/any
117  Array<bool,Dynamic,Dynamic> mb(rows,cols);
118  mb = (m1.real()<=0.7).colwise().all();
119  VERIFY( (mb.col(c) == (m1.real().col(c)<=0.7).all()).all() );
120  mb = (m1.real()<=0.7).rowwise().all();
121  VERIFY( (mb.row(r) == (m1.real().row(r)<=0.7).all()).all() );
122 
123  mb = (m1.real()>=0.7).colwise().any();
124  VERIFY( (mb.col(c) == (m1.real().col(c)>=0.7).any()).all() );
125  mb = (m1.real()>=0.7).rowwise().any();
126  VERIFY( (mb.row(r) == (m1.real().row(r)>=0.7).any()).all() );
127 }
128 
129 template<typename MatrixType> void vectorwiseop_matrix(const MatrixType& m)
130 {
131  typedef typename MatrixType::Scalar Scalar;
132  typedef typename NumTraits<Scalar>::Real RealScalar;
133  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType;
134  typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
135  typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealColVectorType;
136  typedef Matrix<RealScalar, 1, MatrixType::ColsAtCompileTime> RealRowVectorType;
137  typedef Matrix<Scalar,Dynamic,Dynamic> MatrixX;
138 
139  Index rows = m.rows();
140  Index cols = m.cols();
141  Index r = internal::random<Index>(0, rows-1),
142  c = internal::random<Index>(0, cols-1);
143 
144  MatrixType m1 = MatrixType::Random(rows, cols),
145  m2(rows, cols),
146  m3(rows, cols);
147 
148  ColVectorType colvec = ColVectorType::Random(rows);
149  RowVectorType rowvec = RowVectorType::Random(cols);
150  RealColVectorType rcres;
151  RealRowVectorType rrres;
152 
153  // test broadcast assignment
154  m2 = m1;
155  m2.colwise() = colvec;
156  for(Index j=0; j<cols; ++j)
157  VERIFY_IS_APPROX(m2.col(j), colvec);
158  m2.rowwise() = rowvec;
159  for(Index i=0; i<rows; ++i)
160  VERIFY_IS_APPROX(m2.row(i), rowvec);
161  if(rows>1)
162  VERIFY_RAISES_ASSERT(m2.colwise() = colvec.transpose());
163  if(cols>1)
164  VERIFY_RAISES_ASSERT(m2.rowwise() = rowvec.transpose());
165 
166  // test addition
167 
168  m2 = m1;
169  m2.colwise() += colvec;
170  VERIFY_IS_APPROX(m2, m1.colwise() + colvec);
171  VERIFY_IS_APPROX(m2.col(c), m1.col(c) + colvec);
172 
173  if(rows>1)
174  {
175  VERIFY_RAISES_ASSERT(m2.colwise() += colvec.transpose());
176  VERIFY_RAISES_ASSERT(m1.colwise() + colvec.transpose());
177  }
178 
179  m2 = m1;
180  m2.rowwise() += rowvec;
181  VERIFY_IS_APPROX(m2, m1.rowwise() + rowvec);
182  VERIFY_IS_APPROX(m2.row(r), m1.row(r) + rowvec);
183 
184  if(cols>1)
185  {
186  VERIFY_RAISES_ASSERT(m2.rowwise() += rowvec.transpose());
187  VERIFY_RAISES_ASSERT(m1.rowwise() + rowvec.transpose());
188  }
189 
190  // test substraction
191 
192  m2 = m1;
193  m2.colwise() -= colvec;
194  VERIFY_IS_APPROX(m2, m1.colwise() - colvec);
195  VERIFY_IS_APPROX(m2.col(c), m1.col(c) - colvec);
196 
197  if(rows>1)
198  {
199  VERIFY_RAISES_ASSERT(m2.colwise() -= colvec.transpose());
200  VERIFY_RAISES_ASSERT(m1.colwise() - colvec.transpose());
201  }
202 
203  m2 = m1;
204  m2.rowwise() -= rowvec;
205  VERIFY_IS_APPROX(m2, m1.rowwise() - rowvec);
206  VERIFY_IS_APPROX(m2.row(r), m1.row(r) - rowvec);
207 
208  if(cols>1)
209  {
210  VERIFY_RAISES_ASSERT(m2.rowwise() -= rowvec.transpose());
211  VERIFY_RAISES_ASSERT(m1.rowwise() - rowvec.transpose());
212  }
213 
214  // ------ partial reductions ------
215 
216  #define TEST_PARTIAL_REDUX_BASIC(FUNC,ROW,COL,PREPROCESS) { \
217  ROW = m1 PREPROCESS .colwise().FUNC ; \
218  for(Index k=0; k<cols; ++k) VERIFY_IS_APPROX(ROW(k), m1.col(k) PREPROCESS .FUNC ); \
219  COL = m1 PREPROCESS .rowwise().FUNC ; \
220  for(Index k=0; k<rows; ++k) VERIFY_IS_APPROX(COL(k), m1.row(k) PREPROCESS .FUNC ); \
221  }
222 
223  TEST_PARTIAL_REDUX_BASIC(sum(), rowvec,colvec,EIGEN_EMPTY);
224  TEST_PARTIAL_REDUX_BASIC(prod(), rowvec,colvec,EIGEN_EMPTY);
225  TEST_PARTIAL_REDUX_BASIC(mean(), rowvec,colvec,EIGEN_EMPTY);
226  TEST_PARTIAL_REDUX_BASIC(minCoeff(), rrres, rcres, .real());
227  TEST_PARTIAL_REDUX_BASIC(maxCoeff(), rrres, rcres, .real());
228  TEST_PARTIAL_REDUX_BASIC(norm(), rrres, rcres, EIGEN_EMPTY);
229  TEST_PARTIAL_REDUX_BASIC(squaredNorm(),rrres, rcres, EIGEN_EMPTY);
230  TEST_PARTIAL_REDUX_BASIC(redux(internal::scalar_sum_op<Scalar,Scalar>()),rowvec,colvec,EIGEN_EMPTY);
231 
232  VERIFY_IS_APPROX(m1.cwiseAbs().colwise().sum(), m1.colwise().template lpNorm<1>());
233  VERIFY_IS_APPROX(m1.cwiseAbs().rowwise().sum(), m1.rowwise().template lpNorm<1>());
234  VERIFY_IS_APPROX(m1.cwiseAbs().colwise().maxCoeff(), m1.colwise().template lpNorm<Infinity>());
235  VERIFY_IS_APPROX(m1.cwiseAbs().rowwise().maxCoeff(), m1.rowwise().template lpNorm<Infinity>());
236 
237  // regression for bug 1158
238  VERIFY_IS_APPROX(m1.cwiseAbs().colwise().sum().x(), m1.col(0).cwiseAbs().sum());
239 
240  // test normalized
241  m2 = m1.colwise().normalized();
242  VERIFY_IS_APPROX(m2.col(c), m1.col(c).normalized());
243  m2 = m1.rowwise().normalized();
244  VERIFY_IS_APPROX(m2.row(r), m1.row(r).normalized());
245 
246  // test normalize
247  m2 = m1;
248  m2.colwise().normalize();
249  VERIFY_IS_APPROX(m2.col(c), m1.col(c).normalized());
250  m2 = m1;
251  m2.rowwise().normalize();
252  VERIFY_IS_APPROX(m2.row(r), m1.row(r).normalized());
253 
254  // test with partial reduction of products
255  Matrix<Scalar,MatrixType::RowsAtCompileTime,MatrixType::RowsAtCompileTime> m1m1 = m1 * m1.transpose();
256  VERIFY_IS_APPROX( (m1 * m1.transpose()).colwise().sum(), m1m1.colwise().sum());
257  Matrix<Scalar,1,MatrixType::RowsAtCompileTime> tmp(rows);
258  VERIFY_EVALUATION_COUNT( tmp = (m1 * m1.transpose()).colwise().sum(), 1);
259 
260  m2 = m1.rowwise() - (m1.colwise().sum()/RealScalar(m1.rows())).eval();
261  m1 = m1.rowwise() - (m1.colwise().sum()/RealScalar(m1.rows()));
262  VERIFY_IS_APPROX( m1, m2 );
263  VERIFY_EVALUATION_COUNT( m2 = (m1.rowwise() - m1.colwise().sum()/RealScalar(m1.rows())), (MatrixType::RowsAtCompileTime!=1 ? 1 : 0) );
264 
265  // test empty expressions
266  VERIFY_IS_APPROX(m1.matrix().middleCols(0,0).rowwise().sum().eval(), MatrixX::Zero(rows,1));
267  VERIFY_IS_APPROX(m1.matrix().middleRows(0,0).colwise().sum().eval(), MatrixX::Zero(1,cols));
268  VERIFY_IS_APPROX(m1.matrix().middleCols(0,fix<0>).rowwise().sum().eval(), MatrixX::Zero(rows,1));
269  VERIFY_IS_APPROX(m1.matrix().middleRows(0,fix<0>).colwise().sum().eval(), MatrixX::Zero(1,cols));
270 
271  VERIFY_IS_APPROX(m1.matrix().middleCols(0,0).rowwise().prod().eval(), MatrixX::Ones(rows,1));
272  VERIFY_IS_APPROX(m1.matrix().middleRows(0,0).colwise().prod().eval(), MatrixX::Ones(1,cols));
273  VERIFY_IS_APPROX(m1.matrix().middleCols(0,fix<0>).rowwise().prod().eval(), MatrixX::Ones(rows,1));
274  VERIFY_IS_APPROX(m1.matrix().middleRows(0,fix<0>).colwise().prod().eval(), MatrixX::Ones(1,cols));
275 
276  VERIFY_IS_APPROX(m1.matrix().middleCols(0,0).rowwise().squaredNorm().eval(), MatrixX::Zero(rows,1));
277 
278  VERIFY_RAISES_ASSERT(m1.real().middleCols(0,0).rowwise().minCoeff().eval());
279  VERIFY_RAISES_ASSERT(m1.real().middleRows(0,0).colwise().maxCoeff().eval());
280  VERIFY_IS_EQUAL(m1.real().middleRows(0,0).rowwise().maxCoeff().eval().rows(),0);
281  VERIFY_IS_EQUAL(m1.real().middleCols(0,0).colwise().maxCoeff().eval().cols(),0);
282  VERIFY_IS_EQUAL(m1.real().middleRows(0,fix<0>).rowwise().maxCoeff().eval().rows(),0);
283  VERIFY_IS_EQUAL(m1.real().middleCols(0,fix<0>).colwise().maxCoeff().eval().cols(),0);
284 }
285 
286 EIGEN_DECLARE_TEST(vectorwiseop)
287 {
288  CALL_SUBTEST_1( vectorwiseop_array(Array22cd()) );
289  CALL_SUBTEST_2( vectorwiseop_array(Array<double, 3, 2>()) );
290  CALL_SUBTEST_3( vectorwiseop_array(ArrayXXf(3, 4)) );
291  CALL_SUBTEST_4( vectorwiseop_matrix(Matrix4cf()) );
292  CALL_SUBTEST_5( vectorwiseop_matrix(Matrix4f()) );
293  CALL_SUBTEST_5( vectorwiseop_matrix(Vector4f()) );
294  CALL_SUBTEST_5( vectorwiseop_matrix(Matrix<float,4,5>()) );
295  CALL_SUBTEST_6( vectorwiseop_matrix(MatrixXd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
296  CALL_SUBTEST_7( vectorwiseop_matrix(VectorXd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
297  CALL_SUBTEST_7( vectorwiseop_matrix(RowVectorXd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
298 }
TEST_PARTIAL_REDUX_BASIC
#define TEST_PARTIAL_REDUX_BASIC(FUNC, ROW, COL, PREPROCESS)
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
EIGEN_EMPTY
#define EIGEN_EMPTY
Definition: Macros.h:1175
VERIFY_IS_EQUAL
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
c
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
m1
Matrix3d m1
Definition: IOFormat.cpp:2
Eigen::Array
General-purpose arrays with easy API for coefficient-wise operations.
Definition: Array.h:45
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(vectorwiseop)
Definition: vectorwiseop.cpp:286
real
float real
Definition: datatypes.h:10
vectorwiseop_array
void vectorwiseop_array(const ArrayType &m)
Definition: vectorwiseop.cpp:16
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
VERIFY_RAISES_ASSERT
#define VERIFY_RAISES_ASSERT(a)
Definition: main.h:340
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
m2
MatrixType m2(n_dims)
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
Eigen::Dynamic
const int Dynamic
Definition: Constants.h:22
Eigen::all
static const Eigen::internal::all_t all
Definition: IndexedViewHelper.h:171
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
m3
static const DiscreteKey m3(M(3), 2)
main.h
vectorwiseop_matrix
void vectorwiseop_matrix(const MatrixType &m)
Definition: vectorwiseop.cpp:129
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
VERIFY_EVALUATION_COUNT
#define VERIFY_EVALUATION_COUNT(XPR, N)
Definition: test/sparse_product.cpp:27
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
sampling::mean
static const Vector2 mean(20, 40)
CALL_SUBTEST_7
#define CALL_SUBTEST_7(FUNC)
Definition: split_test_helper.h:40
prod
EIGEN_DONT_INLINE void prod(const Lhs &a, const Rhs &b, Res &c)
Definition: product_threshold.cpp:39
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
VERIFY
#define VERIFY(a)
Definition: main.h:380
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74

gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:09:41