array_cwise.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) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #include "main.h"
11 
12 
13 // Test the corner cases of pow(x, y) for real types.
14 template<typename Scalar>
15 void pow_test() {
16  const Scalar zero = Scalar(0);
18  const Scalar one = Scalar(1);
19  const Scalar two = Scalar(2);
20  const Scalar three = Scalar(3);
21  const Scalar sqrt_half = Scalar(std::sqrt(0.5));
22  const Scalar sqrt2 = Scalar(std::sqrt(2));
25  const Scalar denorm_min = std::numeric_limits<Scalar>::denorm_min();
28  const Scalar max_exp = (static_cast<Scalar>(int(Eigen::NumTraits<Scalar>::max_exponent())) * Scalar(EIGEN_LN2)) / eps;
29 
30  const static Scalar abs_vals[] = {zero,
31  denorm_min,
32  min,
33  eps,
34  sqrt_half,
35  one,
36  sqrt2,
37  two,
38  three,
39  max_exp,
40  max,
41  inf,
42  nan};
43  const int abs_cases = 13;
44  const int num_cases = 2*abs_cases * 2*abs_cases;
45  // Repeat the same value to make sure we hit the vectorized path.
46  const int num_repeats = 32;
47  Array<Scalar, Dynamic, Dynamic> x(num_repeats, num_cases);
48  Array<Scalar, Dynamic, Dynamic> y(num_repeats, num_cases);
49  int count = 0;
50  for (int i = 0; i < abs_cases; ++i) {
51  const Scalar abs_x = abs_vals[i];
52  for (int sign_x = 0; sign_x < 2; ++sign_x) {
53  Scalar x_case = sign_x == 0 ? -abs_x : abs_x;
54  for (int j = 0; j < abs_cases; ++j) {
55  const Scalar abs_y = abs_vals[j];
56  for (int sign_y = 0; sign_y < 2; ++sign_y) {
57  Scalar y_case = sign_y == 0 ? -abs_y : abs_y;
58  for (int repeat = 0; repeat < num_repeats; ++repeat) {
59  x(repeat, count) = x_case;
60  y(repeat, count) = y_case;
61  }
62  ++count;
63  }
64  }
65  }
66  }
67 
68  Array<Scalar, Dynamic, Dynamic> actual = x.pow(y);
69  const Scalar tol = test_precision<Scalar>();
70  bool all_pass = true;
71  for (int i = 0; i < 1; ++i) {
72  for (int j = 0; j < num_cases; ++j) {
73  Scalar e = static_cast<Scalar>(std::pow(x(i,j), y(i,j)));
74  Scalar a = actual(i, j);
75  bool fail = !(a==e) && !internal::isApprox(a, e, tol) && !((numext::isnan)(a) && (numext::isnan)(e));
76  all_pass &= !fail;
77  if (fail) {
78  std::cout << "pow(" << x(i,j) << "," << y(i,j) << ") = " << a << " != " << e << std::endl;
79  }
80  }
81  }
82  VERIFY(all_pass);
83 }
84 
85 template<typename ArrayType> void array(const ArrayType& m)
86 {
87  typedef typename ArrayType::Scalar Scalar;
88  typedef typename ArrayType::RealScalar RealScalar;
91 
92  Index rows = m.rows();
93  Index cols = m.cols();
94 
95  ArrayType m1 = ArrayType::Random(rows, cols),
96  m2 = ArrayType::Random(rows, cols),
97  m3(rows, cols);
98  ArrayType m4 = m1; // copy constructor
99  VERIFY_IS_APPROX(m1, m4);
100 
101  ColVectorType cv1 = ColVectorType::Random(rows);
102  RowVectorType rv1 = RowVectorType::Random(cols);
103 
104  Scalar s1 = internal::random<Scalar>(),
105  s2 = internal::random<Scalar>();
106 
107  // scalar addition
108  VERIFY_IS_APPROX(m1 + s1, s1 + m1);
109  VERIFY_IS_APPROX(m1 + s1, ArrayType::Constant(rows,cols,s1) + m1);
110  VERIFY_IS_APPROX(s1 - m1, (-m1)+s1 );
111  VERIFY_IS_APPROX(m1 - s1, m1 - ArrayType::Constant(rows,cols,s1));
112  VERIFY_IS_APPROX(s1 - m1, ArrayType::Constant(rows,cols,s1) - m1);
113  VERIFY_IS_APPROX((m1*Scalar(2)) - s2, (m1+m1) - ArrayType::Constant(rows,cols,s2) );
114  m3 = m1;
115  m3 += s2;
116  VERIFY_IS_APPROX(m3, m1 + s2);
117  m3 = m1;
118  m3 -= s1;
119  VERIFY_IS_APPROX(m3, m1 - s1);
120 
121  // scalar operators via Maps
122  m3 = m1;
123  ArrayType::Map(m1.data(), m1.rows(), m1.cols()) -= ArrayType::Map(m2.data(), m2.rows(), m2.cols());
124  VERIFY_IS_APPROX(m1, m3 - m2);
125 
126  m3 = m1;
127  ArrayType::Map(m1.data(), m1.rows(), m1.cols()) += ArrayType::Map(m2.data(), m2.rows(), m2.cols());
128  VERIFY_IS_APPROX(m1, m3 + m2);
129 
130  m3 = m1;
131  ArrayType::Map(m1.data(), m1.rows(), m1.cols()) *= ArrayType::Map(m2.data(), m2.rows(), m2.cols());
132  VERIFY_IS_APPROX(m1, m3 * m2);
133 
134  m3 = m1;
135  m2 = ArrayType::Random(rows,cols);
136  m2 = (m2==0).select(1,m2);
137  ArrayType::Map(m1.data(), m1.rows(), m1.cols()) /= ArrayType::Map(m2.data(), m2.rows(), m2.cols());
138  VERIFY_IS_APPROX(m1, m3 / m2);
139 
140  // reductions
141  VERIFY_IS_APPROX(m1.abs().colwise().sum().sum(), m1.abs().sum());
142  VERIFY_IS_APPROX(m1.abs().rowwise().sum().sum(), m1.abs().sum());
143  using std::abs;
144  VERIFY_IS_MUCH_SMALLER_THAN(abs(m1.colwise().sum().sum() - m1.sum()), m1.abs().sum());
145  VERIFY_IS_MUCH_SMALLER_THAN(abs(m1.rowwise().sum().sum() - m1.sum()), m1.abs().sum());
146  if (!internal::isMuchSmallerThan(abs(m1.sum() - (m1+m2).sum()), m1.abs().sum(), test_precision<Scalar>()))
147  VERIFY_IS_NOT_APPROX(((m1+m2).rowwise().sum()).sum(), m1.sum());
148  VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op<Scalar,Scalar>()));
149 
150  // vector-wise ops
151  m3 = m1;
152  VERIFY_IS_APPROX(m3.colwise() += cv1, m1.colwise() + cv1);
153  m3 = m1;
154  VERIFY_IS_APPROX(m3.colwise() -= cv1, m1.colwise() - cv1);
155  m3 = m1;
156  VERIFY_IS_APPROX(m3.rowwise() += rv1, m1.rowwise() + rv1);
157  m3 = m1;
158  VERIFY_IS_APPROX(m3.rowwise() -= rv1, m1.rowwise() - rv1);
159 
160  // Conversion from scalar
161  VERIFY_IS_APPROX((m3 = s1), ArrayType::Constant(rows,cols,s1));
162  VERIFY_IS_APPROX((m3 = 1), ArrayType::Constant(rows,cols,1));
163  VERIFY_IS_APPROX((m3.topLeftCorner(rows,cols) = 1), ArrayType::Constant(rows,cols,1));
164  typedef Array<Scalar,
165  ArrayType::RowsAtCompileTime==Dynamic?2:ArrayType::RowsAtCompileTime,
166  ArrayType::ColsAtCompileTime==Dynamic?2:ArrayType::ColsAtCompileTime,
167  ArrayType::Options> FixedArrayType;
168  {
169  FixedArrayType f1(s1);
170  VERIFY_IS_APPROX(f1, FixedArrayType::Constant(s1));
171  FixedArrayType f2(numext::real(s1));
172  VERIFY_IS_APPROX(f2, FixedArrayType::Constant(numext::real(s1)));
173  FixedArrayType f3((int)100*numext::real(s1));
174  VERIFY_IS_APPROX(f3, FixedArrayType::Constant((int)100*numext::real(s1)));
175  f1.setRandom();
176  FixedArrayType f4(f1.data());
177  VERIFY_IS_APPROX(f4, f1);
178  }
179  #if EIGEN_HAS_CXX11
180  {
181  FixedArrayType f1{s1};
182  VERIFY_IS_APPROX(f1, FixedArrayType::Constant(s1));
183  FixedArrayType f2{numext::real(s1)};
184  VERIFY_IS_APPROX(f2, FixedArrayType::Constant(numext::real(s1)));
185  FixedArrayType f3{(int)100*numext::real(s1)};
186  VERIFY_IS_APPROX(f3, FixedArrayType::Constant((int)100*numext::real(s1)));
187  f1.setRandom();
188  FixedArrayType f4{f1.data()};
190  }
191  #endif
192 
193  // pow
194  VERIFY_IS_APPROX(m1.pow(2), m1.square());
195  VERIFY_IS_APPROX(pow(m1,2), m1.square());
196  VERIFY_IS_APPROX(m1.pow(3), m1.cube());
197  VERIFY_IS_APPROX(pow(m1,3), m1.cube());
198  VERIFY_IS_APPROX((-m1).pow(3), -m1.cube());
199  VERIFY_IS_APPROX(pow(2*m1,3), 8*m1.cube());
200  ArrayType exponents = ArrayType::Constant(rows, cols, RealScalar(2));
201  VERIFY_IS_APPROX(Eigen::pow(m1,exponents), m1.square());
202  VERIFY_IS_APPROX(m1.pow(exponents), m1.square());
203  VERIFY_IS_APPROX(Eigen::pow(2*m1,exponents), 4*m1.square());
204  VERIFY_IS_APPROX((2*m1).pow(exponents), 4*m1.square());
205  VERIFY_IS_APPROX(Eigen::pow(m1,2*exponents), m1.square().square());
206  VERIFY_IS_APPROX(m1.pow(2*exponents), m1.square().square());
207  VERIFY_IS_APPROX(Eigen::pow(m1(0,0), exponents), ArrayType::Constant(rows,cols,m1(0,0)*m1(0,0)));
208 
209  // Check possible conflicts with 1D ctor
210  typedef Array<Scalar, Dynamic, 1> OneDArrayType;
211  {
212  OneDArrayType o1(rows);
213  VERIFY(o1.size()==rows);
214  OneDArrayType o2(static_cast<int>(rows));
215  VERIFY(o2.size()==rows);
216  }
217  #if EIGEN_HAS_CXX11
218  {
219  OneDArrayType o1{rows};
220  VERIFY(o1.size()==rows);
221  OneDArrayType o4{int(rows)};
222  VERIFY(o4.size()==rows);
223  }
224  #endif
225  // Check possible conflicts with 2D ctor
226  typedef Array<Scalar, Dynamic, Dynamic> TwoDArrayType;
227  typedef Array<Scalar, 2, 1> ArrayType2;
228  {
229  TwoDArrayType o1(rows,cols);
230  VERIFY(o1.rows()==rows);
231  VERIFY(o1.cols()==cols);
232  TwoDArrayType o2(static_cast<int>(rows),static_cast<int>(cols));
233  VERIFY(o2.rows()==rows);
234  VERIFY(o2.cols()==cols);
235 
236  ArrayType2 o3(rows,cols);
237  VERIFY(o3(0)==Scalar(rows) && o3(1)==Scalar(cols));
238  ArrayType2 o4(static_cast<int>(rows),static_cast<int>(cols));
239  VERIFY(o4(0)==Scalar(rows) && o4(1)==Scalar(cols));
240  }
241  #if EIGEN_HAS_CXX11
242  {
243  TwoDArrayType o1{rows,cols};
244  VERIFY(o1.rows()==rows);
245  VERIFY(o1.cols()==cols);
246  TwoDArrayType o2{int(rows),int(cols)};
247  VERIFY(o2.rows()==rows);
248  VERIFY(o2.cols()==cols);
249 
250  ArrayType2 o3{rows,cols};
251  VERIFY(o3(0)==Scalar(rows) && o3(1)==Scalar(cols));
252  ArrayType2 o4{int(rows),int(cols)};
253  VERIFY(o4(0)==Scalar(rows) && o4(1)==Scalar(cols));
254  }
255  #endif
256 }
257 
258 template<typename ArrayType> void comparisons(const ArrayType& m)
259 {
260  using std::abs;
261  typedef typename ArrayType::Scalar Scalar;
262  typedef typename NumTraits<Scalar>::Real RealScalar;
263 
264  Index rows = m.rows();
265  Index cols = m.cols();
266 
267  Index r = internal::random<Index>(0, rows-1),
268  c = internal::random<Index>(0, cols-1);
269 
270  ArrayType m1 = ArrayType::Random(rows, cols),
271  m2 = ArrayType::Random(rows, cols),
272  m3(rows, cols),
273  m4 = m1;
274 
275  m4 = (m4.abs()==Scalar(0)).select(1,m4);
276 
277  VERIFY(((m1 + Scalar(1)) > m1).all());
278  VERIFY(((m1 - Scalar(1)) < m1).all());
279  if (rows*cols>1)
280  {
281  m3 = m1;
282  m3(r,c) += 1;
283  VERIFY(! (m1 < m3).all() );
284  VERIFY(! (m1 > m3).all() );
285  }
286  VERIFY(!(m1 > m2 && m1 < m2).any());
287  VERIFY((m1 <= m2 || m1 >= m2).all());
288 
289  // comparisons array to scalar
290  VERIFY( (m1 != (m1(r,c)+1) ).any() );
291  VERIFY( (m1 > (m1(r,c)-1) ).any() );
292  VERIFY( (m1 < (m1(r,c)+1) ).any() );
293  VERIFY( (m1 == m1(r,c) ).any() );
294 
295  // comparisons scalar to array
296  VERIFY( ( (m1(r,c)+1) != m1).any() );
297  VERIFY( ( (m1(r,c)-1) < m1).any() );
298  VERIFY( ( (m1(r,c)+1) > m1).any() );
299  VERIFY( ( m1(r,c) == m1).any() );
300 
301  // test Select
302  VERIFY_IS_APPROX( (m1<m2).select(m1,m2), m1.cwiseMin(m2) );
303  VERIFY_IS_APPROX( (m1>m2).select(m1,m2), m1.cwiseMax(m2) );
304  Scalar mid = (m1.cwiseAbs().minCoeff() + m1.cwiseAbs().maxCoeff())/Scalar(2);
305  for (int j=0; j<cols; ++j)
306  for (int i=0; i<rows; ++i)
307  m3(i,j) = abs(m1(i,j))<mid ? 0 : m1(i,j);
308  VERIFY_IS_APPROX( (m1.abs()<ArrayType::Constant(rows,cols,mid))
309  .select(ArrayType::Zero(rows,cols),m1), m3);
310  // shorter versions:
311  VERIFY_IS_APPROX( (m1.abs()<ArrayType::Constant(rows,cols,mid))
312  .select(0,m1), m3);
313  VERIFY_IS_APPROX( (m1.abs()>=ArrayType::Constant(rows,cols,mid))
314  .select(m1,0), m3);
315  // even shorter version:
316  VERIFY_IS_APPROX( (m1.abs()<mid).select(0,m1), m3);
317 
318  // count
319  VERIFY(((m1.abs()+1)>RealScalar(0.1)).count() == rows*cols);
320 
321  // and/or
322  VERIFY( (m1<RealScalar(0) && m1>RealScalar(0)).count() == 0);
323  VERIFY( (m1<RealScalar(0) || m1>=RealScalar(0)).count() == rows*cols);
324  RealScalar a = m1.abs().mean();
325  VERIFY( (m1<-a || m1>a).count() == (m1.abs()>a).count());
326 
327  typedef Array<Index, Dynamic, 1> ArrayOfIndices;
328 
329  // TODO allows colwise/rowwise for array
330  VERIFY_IS_APPROX(((m1.abs()+1)>RealScalar(0.1)).colwise().count(), ArrayOfIndices::Constant(cols,rows).transpose());
331  VERIFY_IS_APPROX(((m1.abs()+1)>RealScalar(0.1)).rowwise().count(), ArrayOfIndices::Constant(rows, cols));
332 }
333 
334 template<typename ArrayType> void array_real(const ArrayType& m)
335 {
336  using std::abs;
337  using std::sqrt;
338  typedef typename ArrayType::Scalar Scalar;
339  typedef typename NumTraits<Scalar>::Real RealScalar;
340 
341  Index rows = m.rows();
342  Index cols = m.cols();
343 
344  ArrayType m1 = ArrayType::Random(rows, cols),
345  m2 = ArrayType::Random(rows, cols),
346  m3(rows, cols),
347  m4 = m1;
348 
349  m4 = (m4.abs()==Scalar(0)).select(Scalar(1),m4);
350 
351  Scalar s1 = internal::random<Scalar>();
352 
353  // these tests are mostly to check possible compilation issues with free-functions.
354  VERIFY_IS_APPROX(m1.sin(), sin(m1));
355  VERIFY_IS_APPROX(m1.cos(), cos(m1));
356  VERIFY_IS_APPROX(m1.tan(), tan(m1));
357  VERIFY_IS_APPROX(m1.asin(), asin(m1));
358  VERIFY_IS_APPROX(m1.acos(), acos(m1));
359  VERIFY_IS_APPROX(m1.atan(), atan(m1));
360  VERIFY_IS_APPROX(m1.sinh(), sinh(m1));
361  VERIFY_IS_APPROX(m1.cosh(), cosh(m1));
362  VERIFY_IS_APPROX(m1.tanh(), tanh(m1));
363 #if EIGEN_HAS_CXX11_MATH
364  VERIFY_IS_APPROX(m1.tanh().atanh(), atanh(tanh(m1)));
365  VERIFY_IS_APPROX(m1.sinh().asinh(), asinh(sinh(m1)));
366  VERIFY_IS_APPROX(m1.cosh().acosh(), acosh(cosh(m1)));
367 #endif
368  VERIFY_IS_APPROX(m1.logistic(), logistic(m1));
369 
370  VERIFY_IS_APPROX(m1.arg(), arg(m1));
371  VERIFY_IS_APPROX(m1.round(), round(m1));
372  VERIFY_IS_APPROX(m1.rint(), rint(m1));
373  VERIFY_IS_APPROX(m1.floor(), floor(m1));
374  VERIFY_IS_APPROX(m1.ceil(), ceil(m1));
375  VERIFY((m1.isNaN() == (Eigen::isnan)(m1)).all());
376  VERIFY((m1.isInf() == (Eigen::isinf)(m1)).all());
377  VERIFY((m1.isFinite() == (Eigen::isfinite)(m1)).all());
378  VERIFY_IS_APPROX(m4.inverse(), inverse(m4));
379  VERIFY_IS_APPROX(m1.abs(), abs(m1));
380  VERIFY_IS_APPROX(m1.abs2(), abs2(m1));
381  VERIFY_IS_APPROX(m1.square(), square(m1));
382  VERIFY_IS_APPROX(m1.cube(), cube(m1));
383  VERIFY_IS_APPROX(cos(m1+RealScalar(3)*m2), cos((m1+RealScalar(3)*m2).eval()));
384  VERIFY_IS_APPROX(m1.sign(), sign(m1));
385  VERIFY((m1.sqrt().sign().isNaN() == (Eigen::isnan)(sign(sqrt(m1)))).all());
386 
387  // avoid inf and NaNs so verification doesn't fail
388  m3 = m4.abs();
389  VERIFY_IS_APPROX(m3.sqrt(), sqrt(abs(m3)));
390  VERIFY_IS_APPROX(m3.rsqrt(), Scalar(1)/sqrt(abs(m3)));
391  VERIFY_IS_APPROX(rsqrt(m3), Scalar(1)/sqrt(abs(m3)));
392  VERIFY_IS_APPROX(m3.log(), log(m3));
393  VERIFY_IS_APPROX(m3.log1p(), log1p(m3));
394  VERIFY_IS_APPROX(m3.log10(), log10(m3));
395  VERIFY_IS_APPROX(m3.log2(), log2(m3));
396 
397 
398  VERIFY((!(m1>m2) == (m1<=m2)).all());
399 
400  VERIFY_IS_APPROX(sin(m1.asin()), m1);
401  VERIFY_IS_APPROX(cos(m1.acos()), m1);
402  VERIFY_IS_APPROX(tan(m1.atan()), m1);
403  VERIFY_IS_APPROX(sinh(m1), Scalar(0.5)*(exp(m1)-exp(-m1)));
404  VERIFY_IS_APPROX(cosh(m1), Scalar(0.5)*(exp(m1)+exp(-m1)));
405  VERIFY_IS_APPROX(tanh(m1), (Scalar(0.5)*(exp(m1)-exp(-m1)))/(Scalar(0.5)*(exp(m1)+exp(-m1))));
406  VERIFY_IS_APPROX(logistic(m1), (Scalar(1)/(Scalar(1)+exp(-m1))));
407  VERIFY_IS_APPROX(arg(m1), ((m1<Scalar(0)).template cast<Scalar>())*Scalar(std::acos(Scalar(-1))));
408  VERIFY((round(m1) <= ceil(m1) && round(m1) >= floor(m1)).all());
409  VERIFY((rint(m1) <= ceil(m1) && rint(m1) >= floor(m1)).all());
410  VERIFY(((ceil(m1) - round(m1)) <= Scalar(0.5) || (round(m1) - floor(m1)) <= Scalar(0.5)).all());
411  VERIFY(((ceil(m1) - round(m1)) <= Scalar(1.0) && (round(m1) - floor(m1)) <= Scalar(1.0)).all());
412  VERIFY(((ceil(m1) - rint(m1)) <= Scalar(0.5) || (rint(m1) - floor(m1)) <= Scalar(0.5)).all());
413  VERIFY(((ceil(m1) - rint(m1)) <= Scalar(1.0) && (rint(m1) - floor(m1)) <= Scalar(1.0)).all());
414  VERIFY((Eigen::isnan)((m1*Scalar(0))/Scalar(0)).all());
415  VERIFY((Eigen::isinf)(m4/Scalar(0)).all());
416  VERIFY(((Eigen::isfinite)(m1) && (!(Eigen::isfinite)(m1*Scalar(0)/Scalar(0))) && (!(Eigen::isfinite)(m4/Scalar(0)))).all());
418  VERIFY((abs(m1) == m1 || abs(m1) == -m1).all());
419  VERIFY_IS_APPROX(m3, sqrt(abs2(m3)));
420  VERIFY_IS_APPROX(m1.absolute_difference(m2), (m1 > m2).select(m1 - m2, m2 - m1));
421  VERIFY_IS_APPROX( m1.sign(), -(-m1).sign() );
422  VERIFY_IS_APPROX( m1*m1.sign(),m1.abs());
423  VERIFY_IS_APPROX(m1.sign() * m1.abs(), m1);
424 
429 
430  // shift argument of logarithm so that it is not zero
431  Scalar smallNumber = NumTraits<Scalar>::dummy_precision();
432  VERIFY_IS_APPROX((m3 + smallNumber).log() , log(abs(m3) + smallNumber));
433  VERIFY_IS_APPROX((m3 + smallNumber + Scalar(1)).log() , log1p(abs(m3) + smallNumber));
434 
435  VERIFY_IS_APPROX(m1.exp() * m2.exp(), exp(m1+m2));
436  VERIFY_IS_APPROX(m1.exp(), exp(m1));
437  VERIFY_IS_APPROX(m1.exp() / m2.exp(),(m1-m2).exp());
438 
439  VERIFY_IS_APPROX(m1.expm1(), expm1(m1));
440  VERIFY_IS_APPROX((m3 + smallNumber).exp() - Scalar(1), expm1(abs(m3) + smallNumber));
441 
442  VERIFY_IS_APPROX(m3.pow(RealScalar(0.5)), m3.sqrt());
443  VERIFY_IS_APPROX(pow(m3,RealScalar(0.5)), m3.sqrt());
444 
445  VERIFY_IS_APPROX(m3.pow(RealScalar(-0.5)), m3.rsqrt());
446  VERIFY_IS_APPROX(pow(m3,RealScalar(-0.5)), m3.rsqrt());
447 
448  // Avoid inf and NaN.
449  m3 = (m1.square()<NumTraits<Scalar>::epsilon()).select(Scalar(1),m3);
450  VERIFY_IS_APPROX(m3.pow(RealScalar(-2)), m3.square().inverse());
451  pow_test<Scalar>();
452 
455 
456  // scalar by array division
457  const RealScalar tiny = sqrt(std::numeric_limits<RealScalar>::epsilon());
458  s1 += Scalar(tiny);
459  m1 += ArrayType::Constant(rows,cols,Scalar(tiny));
460  VERIFY_IS_APPROX(s1/m1, s1 * m1.inverse());
461 
462  // check inplace transpose
463  m3 = m1;
464  m3.transposeInPlace();
465  VERIFY_IS_APPROX(m3, m1.transpose());
466  m3.transposeInPlace();
467  VERIFY_IS_APPROX(m3, m1);
468 }
469 
470 template<typename ArrayType> void array_complex(const ArrayType& m)
471 {
472  typedef typename ArrayType::Scalar Scalar;
473  typedef typename NumTraits<Scalar>::Real RealScalar;
474 
475  Index rows = m.rows();
476  Index cols = m.cols();
477 
478  ArrayType m1 = ArrayType::Random(rows, cols),
479  m2(rows, cols),
480  m4 = m1;
481 
482  m4.real() = (m4.real().abs()==RealScalar(0)).select(RealScalar(1),m4.real());
483  m4.imag() = (m4.imag().abs()==RealScalar(0)).select(RealScalar(1),m4.imag());
484 
485  Array<RealScalar, -1, -1> m3(rows, cols);
486 
487  for (Index i = 0; i < m.rows(); ++i)
488  for (Index j = 0; j < m.cols(); ++j)
489  m2(i,j) = sqrt(m1(i,j));
490 
491  // these tests are mostly to check possible compilation issues with free-functions.
492  VERIFY_IS_APPROX(m1.sin(), sin(m1));
493  VERIFY_IS_APPROX(m1.cos(), cos(m1));
494  VERIFY_IS_APPROX(m1.tan(), tan(m1));
495  VERIFY_IS_APPROX(m1.sinh(), sinh(m1));
496  VERIFY_IS_APPROX(m1.cosh(), cosh(m1));
497  VERIFY_IS_APPROX(m1.tanh(), tanh(m1));
498  VERIFY_IS_APPROX(m1.logistic(), logistic(m1));
499  VERIFY_IS_APPROX(m1.arg(), arg(m1));
500  VERIFY((m1.isNaN() == (Eigen::isnan)(m1)).all());
501  VERIFY((m1.isInf() == (Eigen::isinf)(m1)).all());
502  VERIFY((m1.isFinite() == (Eigen::isfinite)(m1)).all());
503  VERIFY_IS_APPROX(m4.inverse(), inverse(m4));
504  VERIFY_IS_APPROX(m1.log(), log(m1));
505  VERIFY_IS_APPROX(m1.log10(), log10(m1));
506  VERIFY_IS_APPROX(m1.log2(), log2(m1));
507  VERIFY_IS_APPROX(m1.abs(), abs(m1));
508  VERIFY_IS_APPROX(m1.abs2(), abs2(m1));
509  VERIFY_IS_APPROX(m1.sqrt(), sqrt(m1));
510  VERIFY_IS_APPROX(m1.square(), square(m1));
511  VERIFY_IS_APPROX(m1.cube(), cube(m1));
512  VERIFY_IS_APPROX(cos(m1+RealScalar(3)*m2), cos((m1+RealScalar(3)*m2).eval()));
513  VERIFY_IS_APPROX(m1.sign(), sign(m1));
514 
515 
516  VERIFY_IS_APPROX(m1.exp() * m2.exp(), exp(m1+m2));
517  VERIFY_IS_APPROX(m1.exp(), exp(m1));
518  VERIFY_IS_APPROX(m1.exp() / m2.exp(),(m1-m2).exp());
519 
520  VERIFY_IS_APPROX(m1.expm1(), expm1(m1));
521  VERIFY_IS_APPROX(expm1(m1), exp(m1) - 1.);
522  // Check for larger magnitude complex numbers that expm1 matches exp - 1.
523  VERIFY_IS_APPROX(expm1(10. * m1), exp(10. * m1) - 1.);
524 
525  VERIFY_IS_APPROX(sinh(m1), 0.5*(exp(m1)-exp(-m1)));
526  VERIFY_IS_APPROX(cosh(m1), 0.5*(exp(m1)+exp(-m1)));
527  VERIFY_IS_APPROX(tanh(m1), (0.5*(exp(m1)-exp(-m1)))/(0.5*(exp(m1)+exp(-m1))));
528  VERIFY_IS_APPROX(logistic(m1), (1.0/(1.0 + exp(-m1))));
529 
530  for (Index i = 0; i < m.rows(); ++i)
531  for (Index j = 0; j < m.cols(); ++j)
532  m3(i,j) = std::atan2(m1(i,j).imag(), m1(i,j).real());
533  VERIFY_IS_APPROX(arg(m1), m3);
534 
535  std::complex<RealScalar> zero(0.0,0.0);
536  VERIFY((Eigen::isnan)(m1*zero/zero).all());
537 #if EIGEN_COMP_MSVC
538  // msvc complex division is not robust
539  VERIFY((Eigen::isinf)(m4/RealScalar(0)).all());
540 #else
541 #if EIGEN_COMP_CLANG
542  // clang's complex division is notoriously broken too
543  if((numext::isinf)(m4(0,0)/RealScalar(0))) {
544 #endif
545  VERIFY((Eigen::isinf)(m4/zero).all());
546 #if EIGEN_COMP_CLANG
547  }
548  else
549  {
550  VERIFY((Eigen::isinf)(m4.real()/zero.real()).all());
551  }
552 #endif
553 #endif // MSVC
554 
555  VERIFY(((Eigen::isfinite)(m1) && (!(Eigen::isfinite)(m1*zero/zero)) && (!(Eigen::isfinite)(m1/zero))).all());
556 
558  VERIFY_IS_APPROX(conj(m1.conjugate()), m1);
559  VERIFY_IS_APPROX(abs(m1), sqrt(square(m1.real())+square(m1.imag())));
560  VERIFY_IS_APPROX(abs(m1), sqrt(abs2(m1)));
561  VERIFY_IS_APPROX(log10(m1), log(m1)/log(10));
562  VERIFY_IS_APPROX(log2(m1), log(m1)/log(2));
563 
564  VERIFY_IS_APPROX( m1.sign(), -(-m1).sign() );
565  VERIFY_IS_APPROX( m1.sign() * m1.abs(), m1);
566 
567  // scalar by array division
568  Scalar s1 = internal::random<Scalar>();
569  const RealScalar tiny = std::sqrt(std::numeric_limits<RealScalar>::epsilon());
570  s1 += Scalar(tiny);
571  m1 += ArrayType::Constant(rows,cols,Scalar(tiny));
572  VERIFY_IS_APPROX(s1/m1, s1 * m1.inverse());
573 
574  // check inplace transpose
575  m2 = m1;
576  m2.transposeInPlace();
577  VERIFY_IS_APPROX(m2, m1.transpose());
578  m2.transposeInPlace();
579  VERIFY_IS_APPROX(m2, m1);
580  // Check vectorized inplace transpose.
581  ArrayType m5 = ArrayType::Random(131, 131);
582  ArrayType m6 = m5;
583  m6.transposeInPlace();
584  VERIFY_IS_APPROX(m6, m5.transpose());
585 }
586 
587 template<typename ArrayType> void min_max(const ArrayType& m)
588 {
589  typedef typename ArrayType::Scalar Scalar;
590 
591  Index rows = m.rows();
592  Index cols = m.cols();
593 
594  ArrayType m1 = ArrayType::Random(rows, cols);
595 
596  // min/max with array
597  Scalar maxM1 = m1.maxCoeff();
598  Scalar minM1 = m1.minCoeff();
599 
600  VERIFY_IS_APPROX(ArrayType::Constant(rows,cols, minM1), (m1.min)(ArrayType::Constant(rows,cols, minM1)));
601  VERIFY_IS_APPROX(m1, (m1.min)(ArrayType::Constant(rows,cols, maxM1)));
602 
603  VERIFY_IS_APPROX(ArrayType::Constant(rows,cols, maxM1), (m1.max)(ArrayType::Constant(rows,cols, maxM1)));
604  VERIFY_IS_APPROX(m1, (m1.max)(ArrayType::Constant(rows,cols, minM1)));
605 
606  // min/max with scalar input
607  VERIFY_IS_APPROX(ArrayType::Constant(rows,cols, minM1), (m1.min)( minM1));
608  VERIFY_IS_APPROX(m1, (m1.min)( maxM1));
609 
610  VERIFY_IS_APPROX(ArrayType::Constant(rows,cols, maxM1), (m1.max)( maxM1));
611  VERIFY_IS_APPROX(m1, (m1.max)( minM1));
612 
613 
614  // min/max with various NaN propagation options.
615  if (m1.size() > 1 && !NumTraits<Scalar>::IsInteger) {
617  maxM1 = m1.template maxCoeff<PropagateNaN>();
618  minM1 = m1.template minCoeff<PropagateNaN>();
619  VERIFY((numext::isnan)(maxM1));
620  VERIFY((numext::isnan)(minM1));
621 
622  maxM1 = m1.template maxCoeff<PropagateNumbers>();
623  minM1 = m1.template minCoeff<PropagateNumbers>();
624  VERIFY(!(numext::isnan)(maxM1));
625  VERIFY(!(numext::isnan)(minM1));
626  }
627 }
628 
629 template<int N>
630 struct shift_left {
631  template<typename Scalar>
632  Scalar operator()(const Scalar& v) const {
633  return v << N;
634  }
635 };
636 
637 template<int N>
639  template<typename Scalar>
640  Scalar operator()(const Scalar& v) const {
641  return v >> N;
642  }
643 };
644 
645 template<typename ArrayType> void array_integer(const ArrayType& m)
646 {
647  Index rows = m.rows();
648  Index cols = m.cols();
649 
650  ArrayType m1 = ArrayType::Random(rows, cols),
651  m2(rows, cols);
652 
653  m2 = m1.template shiftLeft<2>();
654  VERIFY( (m2 == m1.unaryExpr(shift_left<2>())).all() );
655  m2 = m1.template shiftLeft<9>();
656  VERIFY( (m2 == m1.unaryExpr(shift_left<9>())).all() );
657 
658  m2 = m1.template shiftRight<2>();
659  VERIFY( (m2 == m1.unaryExpr(arithmetic_shift_right<2>())).all() );
660  m2 = m1.template shiftRight<9>();
661  VERIFY( (m2 == m1.unaryExpr(arithmetic_shift_right<9>())).all() );
662 }
663 
664 EIGEN_DECLARE_TEST(array_cwise)
665 {
666  for(int i = 0; i < g_repeat; i++) {
668  CALL_SUBTEST_2( array(Array22f()) );
669  CALL_SUBTEST_3( array(Array44d()) );
670  CALL_SUBTEST_4( array(ArrayXXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
671  CALL_SUBTEST_5( array(ArrayXXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
672  CALL_SUBTEST_6( array(ArrayXXi(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
673  CALL_SUBTEST_6( array(Array<Index,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
674  CALL_SUBTEST_6( array_integer(ArrayXXi(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
675  CALL_SUBTEST_6( array_integer(Array<Index,Dynamic,Dynamic>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
676  }
677  for(int i = 0; i < g_repeat; i++) {
679  CALL_SUBTEST_2( comparisons(Array22f()) );
680  CALL_SUBTEST_3( comparisons(Array44d()) );
681  CALL_SUBTEST_5( comparisons(ArrayXXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
682  CALL_SUBTEST_6( comparisons(ArrayXXi(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
683  }
684  for(int i = 0; i < g_repeat; i++) {
686  CALL_SUBTEST_2( min_max(Array22f()) );
687  CALL_SUBTEST_3( min_max(Array44d()) );
688  CALL_SUBTEST_5( min_max(ArrayXXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
689  CALL_SUBTEST_6( min_max(ArrayXXi(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
690  }
691  for(int i = 0; i < g_repeat; i++) {
693  CALL_SUBTEST_2( array_real(Array22f()) );
694  CALL_SUBTEST_3( array_real(Array44d()) );
695  CALL_SUBTEST_5( array_real(ArrayXXf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
698  }
699  for(int i = 0; i < g_repeat; i++) {
700  CALL_SUBTEST_4( array_complex(ArrayXXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
701  }
702 
706  typedef CwiseUnaryOp<internal::scalar_abs_op<double>, ArrayXd > Xpr;
709  >::value));
710 }
EIGEN_DEVICE_FUNC const Log1pReturnType log1p() const
Matrix3f m
SCALAR Scalar
Definition: bench_gemm.cpp:46
Scalar operator()(const Scalar &v) const
#define max(a, b)
Definition: datatypes.h:20
EIGEN_DEVICE_FUNC bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
#define CALL_SUBTEST_6(FUNC)
#define CALL_SUBTEST_4(FUNC)
Jet< T, N > cos(const Jet< T, N > &f)
Definition: jet.h:426
float real
Definition: datatypes.h:10
void array_real(const ArrayType &m)
static const Eigen::internal::all_t all
Scalar * y
EIGEN_DECLARE_TEST(array_cwise)
void array_integer(const ArrayType &m)
#define min(a, b)
Definition: datatypes.h:19
#define CALL_SUBTEST_3(FUNC)
MatrixType m2(n_dims)
EIGEN_DONT_INLINE Scalar zero()
Definition: svd_common.h:296
#define CALL_SUBTEST_7(FUNC)
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
void comparisons(const ArrayType &m)
Jet< T, N > acos(const Jet< T, N > &f)
Definition: jet.h:432
EIGEN_DEVICE_FUNC const TanhReturnType tanh() const
Jet< T, N > sin(const Jet< T, N > &f)
Definition: jet.h:439
DerType::Scalar imag(const AutoDiffScalar< DerType > &)
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
double f2(const Vector2 &x)
#define isfinite(X)
Definition: main.h:95
void array(const ArrayType &m)
Definition: array_cwise.cpp:85
#define N
Definition: gksort.c:12
#define isinf(X)
Definition: main.h:94
EIGEN_DEVICE_FUNC const SinhReturnType sinh() const
EIGEN_DEVICE_FUNC const LogReturnType log() const
EIGEN_DEVICE_FUNC const InverseReturnType inverse() const
EIGEN_DEVICE_FUNC const CubeReturnType cube() const
AnnoyingScalar conj(const AnnoyingScalar &x)
static double epsilon
Definition: testRot3.cpp:37
EIGEN_DEVICE_FUNC const AtanReturnType atan() const
EIGEN_DEVICE_FUNC const ExpReturnType exp() const
void array_complex(const ArrayType &m)
EIGEN_DEVICE_FUNC const FloorReturnType floor() const
#define VERIFY_IS_APPROX(a, b)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Abs2ReturnType abs2() const
EIGEN_DEVICE_FUNC const SquareReturnType square() const
void pow_test()
Definition: array_cwise.cpp:15
EIGEN_DEVICE_FUNC const Expm1ReturnType expm1() const
#define CALL_SUBTEST_1(FUNC)
EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC bfloat16 log2(const bfloat16 &a)
Definition: BFloat16.h:508
Matrix3d m1
Definition: IOFormat.cpp:2
EIGEN_DEVICE_FUNC const RintReturnType rint() const
static int g_repeat
Definition: main.h:169
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const ArgReturnType arg() const
Array< int, Dynamic, 1 > v
#define EIGEN_LN2
EIGEN_DEVICE_FUNC const SignReturnType sign() const
#define CALL_SUBTEST_8(FUNC)
Array< double, 1, 3 > e(1./3., 0.5, 2.)
Base class for all 1D and 2D array, and related expressions.
Definition: ArrayBase.h:39
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
EIGEN_DEVICE_FUNC const RsqrtReturnType rsqrt() const
EIGEN_DEVICE_FUNC const Eigen::ArrayBase< Derived > & exponents
AnnoyingScalar atan2(const AnnoyingScalar &y, const AnnoyingScalar &x)
#define VERIFY_IS_MUCH_SMALLER_THAN(a, b)
Definition: main.h:390
void min_max(const ArrayType &m)
#define CALL_SUBTEST_5(FUNC)
Scalar operator()(const Scalar &v) const
static double inf
Definition: testMatrix.cpp:31
#define VERIFY(a)
Definition: main.h:380
#define EIGEN_TEST_MAX_SIZE
EIGEN_DEVICE_FUNC const Log10ReturnType log10() const
Point2 f1(const Point3 &p, OptionalJacobian< 2, 3 > H)
double f4(double x, double y, double z)
EIGEN_DEVICE_FUNC const TanReturnType tan() const
General-purpose arrays with easy API for coefficient-wise operations.
Definition: Array.h:45
double f3(double x1, double x2)
EIGEN_DEVICE_FUNC const ImagReturnType imag() const
#define VERIFY_IS_NOT_APPROX(a, b)
#define CALL_SUBTEST_2(FUNC)
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
const G double tol
Definition: Group.h:86
const int Dynamic
Definition: Constants.h:22
Jet< T, N > pow(const Jet< T, N > &f, double g)
Definition: jet.h:570
Generic expression where a coefficient-wise unary operator is applied to an expression.
Definition: CwiseUnaryOp.h:55
EIGEN_DEVICE_FUNC const LogisticReturnType logistic() const
EIGEN_DEVICE_FUNC bool isApprox(const Scalar &x, const Scalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
EIGEN_DEVICE_FUNC const CoshReturnType cosh() const
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
const AutoDiffScalar< DerType > & real(const AutoDiffScalar< DerType > &x)
#define abs(x)
Definition: datatypes.h:17
EIGEN_DEVICE_FUNC const AsinReturnType asin() const
EIGEN_DEVICE_FUNC const CeilReturnType ceil() const
std::ptrdiff_t j
#define isnan(X)
Definition: main.h:93
constexpr array< t, n > repeat(t v)
Definition: CXX11Meta.h:483
EIGEN_DEVICE_FUNC const RoundReturnType round() const


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:33:54