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