nullary.cpp 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341
  1. // This file is part of Eigen, a lightweight C++ template library
  2. // for linear algebra.
  3. //
  4. // Copyright (C) 2010-2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
  5. // Copyright (C) 2016 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. #include "main.h"
  11. template<typename MatrixType>
  12. bool equalsIdentity(const MatrixType& A)
  13. {
  14. typedef typename MatrixType::Scalar Scalar;
  15. Scalar zero = static_cast<Scalar>(0);
  16. bool offDiagOK = true;
  17. for (Index i = 0; i < A.rows(); ++i) {
  18. for (Index j = i+1; j < A.cols(); ++j) {
  19. offDiagOK = offDiagOK && (A(i,j) == zero);
  20. }
  21. }
  22. for (Index i = 0; i < A.rows(); ++i) {
  23. for (Index j = 0; j < (std::min)(i, A.cols()); ++j) {
  24. offDiagOK = offDiagOK && (A(i,j) == zero);
  25. }
  26. }
  27. bool diagOK = (A.diagonal().array() == 1).all();
  28. return offDiagOK && diagOK;
  29. }
  30. template<typename VectorType>
  31. void check_extremity_accuracy(const VectorType &v, const typename VectorType::Scalar &low, const typename VectorType::Scalar &high)
  32. {
  33. typedef typename VectorType::Scalar Scalar;
  34. typedef typename VectorType::RealScalar RealScalar;
  35. RealScalar prec = internal::is_same<RealScalar,float>::value ? NumTraits<RealScalar>::dummy_precision()*10 : NumTraits<RealScalar>::dummy_precision()/10;
  36. Index size = v.size();
  37. if(size<20)
  38. return;
  39. for (int i=0; i<size; ++i)
  40. {
  41. if(i<5 || i>size-6)
  42. {
  43. Scalar ref = (low*RealScalar(size-i-1))/RealScalar(size-1) + (high*RealScalar(i))/RealScalar(size-1);
  44. if(std::abs(ref)>1)
  45. {
  46. if(!internal::isApprox(v(i), ref, prec))
  47. std::cout << v(i) << " != " << ref << " ; relative error: " << std::abs((v(i)-ref)/ref) << " ; required precision: " << prec << " ; range: " << low << "," << high << " ; i: " << i << "\n";
  48. VERIFY(internal::isApprox(v(i), (low*RealScalar(size-i-1))/RealScalar(size-1) + (high*RealScalar(i))/RealScalar(size-1), prec));
  49. }
  50. }
  51. }
  52. }
  53. template<typename VectorType>
  54. void testVectorType(const VectorType& base)
  55. {
  56. typedef typename VectorType::Scalar Scalar;
  57. typedef typename VectorType::RealScalar RealScalar;
  58. const Index size = base.size();
  59. Scalar high = internal::random<Scalar>(-500,500);
  60. Scalar low = (size == 1 ? high : internal::random<Scalar>(-500,500));
  61. if (numext::real(low)>numext::real(high)) std::swap(low,high);
  62. // check low==high
  63. if(internal::random<float>(0.f,1.f)<0.05f)
  64. low = high;
  65. // check abs(low) >> abs(high)
  66. else if(size>2 && std::numeric_limits<RealScalar>::max_exponent10>0 && internal::random<float>(0.f,1.f)<0.1f)
  67. low = -internal::random<Scalar>(1,2) * RealScalar(std::pow(RealScalar(10),std::numeric_limits<RealScalar>::max_exponent10/2));
  68. const Scalar step = ((size == 1) ? 1 : (high-low)/RealScalar(size-1));
  69. // check whether the result yields what we expect it to do
  70. VectorType m(base);
  71. m.setLinSpaced(size,low,high);
  72. if(!NumTraits<Scalar>::IsInteger)
  73. {
  74. VectorType n(size);
  75. for (int i=0; i<size; ++i)
  76. n(i) = low+RealScalar(i)*step;
  77. VERIFY_IS_APPROX(m,n);
  78. CALL_SUBTEST( check_extremity_accuracy(m, low, high) );
  79. }
  80. RealScalar range_length = numext::real(high-low);
  81. if((!NumTraits<Scalar>::IsInteger) || (range_length>=size && (Index(range_length)%(size-1))==0) || (Index(range_length+1)<size && (size%Index(range_length+1))==0))
  82. {
  83. VectorType n(size);
  84. if((!NumTraits<Scalar>::IsInteger) || (range_length>=size))
  85. for (int i=0; i<size; ++i)
  86. n(i) = size==1 ? low : (low + ((high-low)*Scalar(i))/RealScalar(size-1));
  87. else
  88. for (int i=0; i<size; ++i)
  89. n(i) = size==1 ? low : low + Scalar((double(range_length+1)*double(i))/double(size));
  90. VERIFY_IS_APPROX(m,n);
  91. // random access version
  92. m = VectorType::LinSpaced(size,low,high);
  93. VERIFY_IS_APPROX(m,n);
  94. VERIFY( internal::isApprox(m(m.size()-1),high) );
  95. VERIFY( size==1 || internal::isApprox(m(0),low) );
  96. VERIFY_IS_EQUAL(m(m.size()-1) , high);
  97. if(!NumTraits<Scalar>::IsInteger)
  98. CALL_SUBTEST( check_extremity_accuracy(m, low, high) );
  99. }
  100. VERIFY( numext::real(m(m.size()-1)) <= numext::real(high) );
  101. VERIFY( (m.array().real() <= numext::real(high)).all() );
  102. VERIFY( (m.array().real() >= numext::real(low)).all() );
  103. VERIFY( numext::real(m(m.size()-1)) >= numext::real(low) );
  104. if(size>=1)
  105. {
  106. VERIFY( internal::isApprox(m(0),low) );
  107. VERIFY_IS_EQUAL(m(0) , low);
  108. }
  109. // check whether everything works with row and col major vectors
  110. Matrix<Scalar,Dynamic,1> row_vector(size);
  111. Matrix<Scalar,1,Dynamic> col_vector(size);
  112. row_vector.setLinSpaced(size,low,high);
  113. col_vector.setLinSpaced(size,low,high);
  114. // when using the extended precision (e.g., FPU) the relative error might exceed 1 bit
  115. // when computing the squared sum in isApprox, thus the 2x factor.
  116. VERIFY( row_vector.isApprox(col_vector.transpose(), RealScalar(2)*NumTraits<Scalar>::epsilon()));
  117. Matrix<Scalar,Dynamic,1> size_changer(size+50);
  118. size_changer.setLinSpaced(size,low,high);
  119. VERIFY( size_changer.size() == size );
  120. typedef Matrix<Scalar,1,1> ScalarMatrix;
  121. ScalarMatrix scalar;
  122. scalar.setLinSpaced(1,low,high);
  123. VERIFY_IS_APPROX( scalar, ScalarMatrix::Constant(high) );
  124. VERIFY_IS_APPROX( ScalarMatrix::LinSpaced(1,low,high), ScalarMatrix::Constant(high) );
  125. // regression test for bug 526 (linear vectorized transversal)
  126. if (size > 1 && (!NumTraits<Scalar>::IsInteger)) {
  127. m.tail(size-1).setLinSpaced(low, high);
  128. VERIFY_IS_APPROX(m(size-1), high);
  129. }
  130. // regression test for bug 1383 (LinSpaced with empty size/range)
  131. {
  132. Index n0 = VectorType::SizeAtCompileTime==Dynamic ? 0 : VectorType::SizeAtCompileTime;
  133. low = internal::random<Scalar>();
  134. m = VectorType::LinSpaced(n0,low,low-RealScalar(1));
  135. VERIFY(m.size()==n0);
  136. if(VectorType::SizeAtCompileTime==Dynamic)
  137. {
  138. VERIFY_IS_EQUAL(VectorType::LinSpaced(n0,0,Scalar(n0-1)).sum(),Scalar(0));
  139. VERIFY_IS_EQUAL(VectorType::LinSpaced(n0,low,low-RealScalar(1)).sum(),Scalar(0));
  140. }
  141. m.setLinSpaced(n0,0,Scalar(n0-1));
  142. VERIFY(m.size()==n0);
  143. m.setLinSpaced(n0,low,low-RealScalar(1));
  144. VERIFY(m.size()==n0);
  145. // empty range only:
  146. VERIFY_IS_APPROX(VectorType::LinSpaced(size,low,low),VectorType::Constant(size,low));
  147. m.setLinSpaced(size,low,low);
  148. VERIFY_IS_APPROX(m,VectorType::Constant(size,low));
  149. if(NumTraits<Scalar>::IsInteger)
  150. {
  151. VERIFY_IS_APPROX( VectorType::LinSpaced(size,low,low+Scalar(size-1)), VectorType::LinSpaced(size,low+Scalar(size-1),low).reverse() );
  152. if(VectorType::SizeAtCompileTime==Dynamic)
  153. {
  154. // Check negative multiplicator path:
  155. for(Index k=1; k<5; ++k)
  156. VERIFY_IS_APPROX( VectorType::LinSpaced(size,low,low+Scalar((size-1)*k)), VectorType::LinSpaced(size,low+Scalar((size-1)*k),low).reverse() );
  157. // Check negative divisor path:
  158. for(Index k=1; k<5; ++k)
  159. VERIFY_IS_APPROX( VectorType::LinSpaced(size*k,low,low+Scalar(size-1)), VectorType::LinSpaced(size*k,low+Scalar(size-1),low).reverse() );
  160. }
  161. }
  162. }
  163. // test setUnit()
  164. if(m.size()>0)
  165. {
  166. for(Index k=0; k<10; ++k)
  167. {
  168. Index i = internal::random<Index>(0,m.size()-1);
  169. m.setUnit(i);
  170. VERIFY_IS_APPROX( m, VectorType::Unit(m.size(), i) );
  171. }
  172. if(VectorType::SizeAtCompileTime==Dynamic)
  173. {
  174. Index i = internal::random<Index>(0,2*m.size()-1);
  175. m.setUnit(2*m.size(),i);
  176. VERIFY_IS_APPROX( m, VectorType::Unit(m.size(),i) );
  177. }
  178. }
  179. }
  180. template<typename MatrixType>
  181. void testMatrixType(const MatrixType& m)
  182. {
  183. using std::abs;
  184. const Index rows = m.rows();
  185. const Index cols = m.cols();
  186. typedef typename MatrixType::Scalar Scalar;
  187. typedef typename MatrixType::RealScalar RealScalar;
  188. Scalar s1;
  189. do {
  190. s1 = internal::random<Scalar>();
  191. } while(abs(s1)<RealScalar(1e-5) && (!NumTraits<Scalar>::IsInteger));
  192. MatrixType A;
  193. A.setIdentity(rows, cols);
  194. VERIFY(equalsIdentity(A));
  195. VERIFY(equalsIdentity(MatrixType::Identity(rows, cols)));
  196. A = MatrixType::Constant(rows,cols,s1);
  197. Index i = internal::random<Index>(0,rows-1);
  198. Index j = internal::random<Index>(0,cols-1);
  199. VERIFY_IS_APPROX( MatrixType::Constant(rows,cols,s1)(i,j), s1 );
  200. VERIFY_IS_APPROX( MatrixType::Constant(rows,cols,s1).coeff(i,j), s1 );
  201. VERIFY_IS_APPROX( A(i,j), s1 );
  202. }
  203. template<int>
  204. void bug79()
  205. {
  206. // Assignment of a RowVectorXd to a MatrixXd (regression test for bug #79).
  207. VERIFY( (MatrixXd(RowVectorXd::LinSpaced(3, 0, 1)) - RowVector3d(0, 0.5, 1)).norm() < std::numeric_limits<double>::epsilon() );
  208. }
  209. template<int>
  210. void bug1630()
  211. {
  212. Array4d x4 = Array4d::LinSpaced(0.0, 1.0);
  213. Array3d x3(Array4d::LinSpaced(0.0, 1.0).head(3));
  214. VERIFY_IS_APPROX(x4.head(3), x3);
  215. }
  216. template<int>
  217. void nullary_overflow()
  218. {
  219. // Check possible overflow issue
  220. int n = 60000;
  221. ArrayXi a1(n), a2(n);
  222. a1.setLinSpaced(n, 0, n-1);
  223. for(int i=0; i<n; ++i)
  224. a2(i) = i;
  225. VERIFY_IS_APPROX(a1,a2);
  226. }
  227. template<int>
  228. void nullary_internal_logic()
  229. {
  230. // check some internal logic
  231. VERIFY(( internal::has_nullary_operator<internal::scalar_constant_op<double> >::value ));
  232. VERIFY(( !internal::has_unary_operator<internal::scalar_constant_op<double> >::value ));
  233. VERIFY(( !internal::has_binary_operator<internal::scalar_constant_op<double> >::value ));
  234. VERIFY(( internal::functor_has_linear_access<internal::scalar_constant_op<double> >::ret ));
  235. VERIFY(( !internal::has_nullary_operator<internal::scalar_identity_op<double> >::value ));
  236. VERIFY(( !internal::has_unary_operator<internal::scalar_identity_op<double> >::value ));
  237. VERIFY(( internal::has_binary_operator<internal::scalar_identity_op<double> >::value ));
  238. VERIFY(( !internal::functor_has_linear_access<internal::scalar_identity_op<double> >::ret ));
  239. VERIFY(( !internal::has_nullary_operator<internal::linspaced_op<float> >::value ));
  240. VERIFY(( internal::has_unary_operator<internal::linspaced_op<float> >::value ));
  241. VERIFY(( !internal::has_binary_operator<internal::linspaced_op<float> >::value ));
  242. VERIFY(( internal::functor_has_linear_access<internal::linspaced_op<float> >::ret ));
  243. // Regression unit test for a weird MSVC bug.
  244. // Search "nullary_wrapper_workaround_msvc" in CoreEvaluators.h for the details.
  245. // See also traits<Ref>::match.
  246. {
  247. MatrixXf A = MatrixXf::Random(3,3);
  248. Ref<const MatrixXf> R = 2.0*A;
  249. VERIFY_IS_APPROX(R, A+A);
  250. Ref<const MatrixXf> R1 = MatrixXf::Random(3,3)+A;
  251. VectorXi V = VectorXi::Random(3);
  252. Ref<const VectorXi> R2 = VectorXi::LinSpaced(3,1,3)+V;
  253. VERIFY_IS_APPROX(R2, V+Vector3i(1,2,3));
  254. VERIFY(( internal::has_nullary_operator<internal::scalar_constant_op<float> >::value ));
  255. VERIFY(( !internal::has_unary_operator<internal::scalar_constant_op<float> >::value ));
  256. VERIFY(( !internal::has_binary_operator<internal::scalar_constant_op<float> >::value ));
  257. VERIFY(( internal::functor_has_linear_access<internal::scalar_constant_op<float> >::ret ));
  258. VERIFY(( !internal::has_nullary_operator<internal::linspaced_op<int> >::value ));
  259. VERIFY(( internal::has_unary_operator<internal::linspaced_op<int> >::value ));
  260. VERIFY(( !internal::has_binary_operator<internal::linspaced_op<int> >::value ));
  261. VERIFY(( internal::functor_has_linear_access<internal::linspaced_op<int> >::ret ));
  262. }
  263. }
  264. EIGEN_DECLARE_TEST(nullary)
  265. {
  266. CALL_SUBTEST_1( testMatrixType(Matrix2d()) );
  267. CALL_SUBTEST_2( testMatrixType(MatrixXcf(internal::random<int>(1,300),internal::random<int>(1,300))) );
  268. CALL_SUBTEST_3( testMatrixType(MatrixXf(internal::random<int>(1,300),internal::random<int>(1,300))) );
  269. for(int i = 0; i < g_repeat*10; i++) {
  270. CALL_SUBTEST_3( testVectorType(VectorXcd(internal::random<int>(1,30000))) );
  271. CALL_SUBTEST_4( testVectorType(VectorXd(internal::random<int>(1,30000))) );
  272. CALL_SUBTEST_5( testVectorType(Vector4d()) ); // regression test for bug 232
  273. CALL_SUBTEST_6( testVectorType(Vector3d()) );
  274. CALL_SUBTEST_7( testVectorType(VectorXf(internal::random<int>(1,30000))) );
  275. CALL_SUBTEST_8( testVectorType(Vector3f()) );
  276. CALL_SUBTEST_8( testVectorType(Vector4f()) );
  277. CALL_SUBTEST_8( testVectorType(Matrix<float,8,1>()) );
  278. CALL_SUBTEST_8( testVectorType(Matrix<float,1,1>()) );
  279. CALL_SUBTEST_9( testVectorType(VectorXi(internal::random<int>(1,10))) );
  280. CALL_SUBTEST_9( testVectorType(VectorXi(internal::random<int>(9,300))) );
  281. CALL_SUBTEST_9( testVectorType(Matrix<int,1,1>()) );
  282. }
  283. CALL_SUBTEST_6( bug79<0>() );
  284. CALL_SUBTEST_6( bug1630<0>() );
  285. CALL_SUBTEST_9( nullary_overflow<0>() );
  286. CALL_SUBTEST_10( nullary_internal_logic<0>() );
  287. }