ground/gcs/src/libs/eigen/test/nullary.cpp

   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 (low>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)/(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+i*step;
  93     VERIFY_IS_APPROX(m,n);
  94
  95     CALL_SUBTEST( check_extremity_accuracy(m, low, high) );
  96   }
  97
  98   if((!NumTraits<Scalar>::IsInteger) || ((high-low)>=size && (Index(high-low)%(size-1))==0) || (Index(high-low+1)<size && (size%Index(high-low+1))==0))
  99   {
 100     VectorType n(size);
 101     if((!NumTraits<Scalar>::IsInteger) || (high-low>=size))
 102       for (int i=0; i<size; ++i)
 103         n(i) = size==1 ? low : (low + ((high-low)*Scalar(i))/(size-1));
 104     else
 105       for (int i=0; i<size; ++i)
 106         n(i) = size==1 ? low : low + Scalar((double(high-low+1)*double(i))/double(size));
 107     VERIFY_IS_APPROX(m,n);
 108
 109     // random access version
 110     m = VectorType::LinSpaced(size,low,high);
 111     VERIFY_IS_APPROX(m,n);
 112     VERIFY( internal::isApprox(m(m.size()-1),high) );
 113     VERIFY( size==1 || internal::isApprox(m(0),low) );
 114     VERIFY_IS_EQUAL(m(m.size()-1) , high);
 115     if(!NumTraits<Scalar>::IsInteger)
 116       CALL_SUBTEST( check_extremity_accuracy(m, low, high) );
 117   }
 118
 119   VERIFY( m(m.size()-1) <= high );
 120   VERIFY( (m.array() <= high).all() );
 121   VERIFY( (m.array() >= low).all() );
 122
 123
 124   VERIFY( m(m.size()-1) >= low );
 125   if(size>=1)
 126   {
 127     VERIFY( internal::isApprox(m(0),low) );
 128     VERIFY_IS_EQUAL(m(0) , low);
 129   }
 130
 131   // check whether everything works with row and col major vectors
 132   Matrix<Scalar,Dynamic,1> row_vector(size);
 133   Matrix<Scalar,1,Dynamic> col_vector(size);
 134   row_vector.setLinSpaced(size,low,high);
 135   col_vector.setLinSpaced(size,low,high);
 136   // when using the extended precision (e.g., FPU) the relative error might exceed 1 bit
 137   // when computing the squared sum in isApprox, thus the 2x factor.
 138   VERIFY( row_vector.isApprox(col_vector.transpose(), Scalar(2)*NumTraits<Scalar>::epsilon()));
 139
 140   Matrix<Scalar,Dynamic,1> size_changer(size+50);
 141   size_changer.setLinSpaced(size,low,high);
 142   VERIFY( size_changer.size() == size );
 143
 144   typedef Matrix<Scalar,1,1> ScalarMatrix;
 145   ScalarMatrix scalar;
 146   scalar.setLinSpaced(1,low,high);
 147   VERIFY_IS_APPROX( scalar, ScalarMatrix::Constant(high) );
 148   VERIFY_IS_APPROX( ScalarMatrix::LinSpaced(1,low,high), ScalarMatrix::Constant(high) );
 149
 150   // regression test for bug 526 (linear vectorized transversal)
 151   if (size > 1 && (!NumTraits<Scalar>::IsInteger)) {
 152     m.tail(size-1).setLinSpaced(low, high);
 153     VERIFY_IS_APPROX(m(size-1), high);
 154   }
 155
 156   // regression test for bug 1383 (LinSpaced with empty size/range)
 157   {
 158     Index n0 = VectorType::SizeAtCompileTime==Dynamic ? 0 : VectorType::SizeAtCompileTime;
 159     low = internal::random<Scalar>();
 160     m = VectorType::LinSpaced(n0,low,low-1);
 161     VERIFY(m.size()==n0);
 162
 163     if(VectorType::SizeAtCompileTime==Dynamic)
 164     {
 165       VERIFY_IS_EQUAL(VectorType::LinSpaced(n0,0,Scalar(n0-1)).sum(),Scalar(0));
 166       VERIFY_IS_EQUAL(VectorType::LinSpaced(n0,low,low-1).sum(),Scalar(0));
 167     }
 168
 169     m.setLinSpaced(n0,0,Scalar(n0-1));
 170     VERIFY(m.size()==n0);
 171     m.setLinSpaced(n0,low,low-1);
 172     VERIFY(m.size()==n0);
 173
 174     // empty range only:
 175     VERIFY_IS_APPROX(VectorType::LinSpaced(size,low,low),VectorType::Constant(size,low));
 176     m.setLinSpaced(size,low,low);
 177     VERIFY_IS_APPROX(m,VectorType::Constant(size,low));
 178
 179     if(NumTraits<Scalar>::IsInteger)
 180     {
 181       VERIFY_IS_APPROX( VectorType::LinSpaced(size,low,Scalar(low+size-1)), VectorType::LinSpaced(size,Scalar(low+size-1),low).reverse() );
 182
 183       if(VectorType::SizeAtCompileTime==Dynamic)
 184       {
 185         // Check negative multiplicator path:
 186         for(Index k=1; k<5; ++k)
 187           VERIFY_IS_APPROX( VectorType::LinSpaced(size,low,Scalar(low+(size-1)*k)), VectorType::LinSpaced(size,Scalar(low+(size-1)*k),low).reverse() );
 188         // Check negative divisor path:
 189         for(Index k=1; k<5; ++k)
 190           VERIFY_IS_APPROX( VectorType::LinSpaced(size*k,low,Scalar(low+size-1)), VectorType::LinSpaced(size*k,Scalar(low+size-1),low).reverse() );
 191       }
 192     }
 193   }
 194 }
 195
 196 template<typename MatrixType>
 197 void testMatrixType(const MatrixType& m)
 198 {
 199   using std::abs;
 200   const Index rows = m.rows();
 201   const Index cols = m.cols();
 202   typedef typename MatrixType::Scalar Scalar;
 203   typedef typename MatrixType::RealScalar RealScalar;
 204
 205   Scalar s1;
 206   do {
 207     s1 = internal::random<Scalar>();
 208   } while(abs(s1)<RealScalar(1e-5) && (!NumTraits<Scalar>::IsInteger));
 209
 210   MatrixType A;
 211   A.setIdentity(rows, cols);
 212   VERIFY(equalsIdentity(A));
 213   VERIFY(equalsIdentity(MatrixType::Identity(rows, cols)));
 214
 215
 216   A = MatrixType::Constant(rows,cols,s1);
 217   Index i = internal::random<Index>(0,rows-1);
 218   Index j = internal::random<Index>(0,cols-1);
 219   VERIFY_IS_APPROX( MatrixType::Constant(rows,cols,s1)(i,j), s1 );
 220   VERIFY_IS_APPROX( MatrixType::Constant(rows,cols,s1).coeff(i,j), s1 );
 221   VERIFY_IS_APPROX( A(i,j), s1 );
 222 }
 223
 224 void test_nullary()
 225 {
 226   CALL_SUBTEST_1( testMatrixType(Matrix2d()) );
 227   CALL_SUBTEST_2( testMatrixType(MatrixXcf(internal::random<int>(1,300),internal::random<int>(1,300))) );
 228   CALL_SUBTEST_3( testMatrixType(MatrixXf(internal::random<int>(1,300),internal::random<int>(1,300))) );
 229
 230   for(int i = 0; i < g_repeat*10; i++) {
 231     CALL_SUBTEST_4( testVectorType(VectorXd(internal::random<int>(1,30000))) );
 232     CALL_SUBTEST_5( testVectorType(Vector4d()) );  // regression test for bug 232
 233     CALL_SUBTEST_6( testVectorType(Vector3d()) );
 234     CALL_SUBTEST_7( testVectorType(VectorXf(internal::random<int>(1,30000))) );
 235     CALL_SUBTEST_8( testVectorType(Vector3f()) );
 236     CALL_SUBTEST_8( testVectorType(Vector4f()) );
 237     CALL_SUBTEST_8( testVectorType(Matrix<float,8,1>()) );
 238     CALL_SUBTEST_8( testVectorType(Matrix<float,1,1>()) );
 239
 240     CALL_SUBTEST_9( testVectorType(VectorXi(internal::random<int>(1,10))) );
 241     CALL_SUBTEST_9( testVectorType(VectorXi(internal::random<int>(9,300))) );
 242     CALL_SUBTEST_9( testVectorType(Matrix<int,1,1>()) );
 243   }
 244
 245 #ifdef EIGEN_TEST_PART_6
 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 #endif
 249
 250 #ifdef EIGEN_TEST_PART_9
 251   // Check possible overflow issue
 252   {
 253     int n = 60000;
 254     ArrayXi a1(n), a2(n);
 255     a1.setLinSpaced(n, 0, n-1);
 256     for(int i=0; i<n; ++i)
 257       a2(i) = i;
 258     VERIFY_IS_APPROX(a1,a2);
 259   }
 260 #endif
 261
 262 #ifdef EIGEN_TEST_PART_10
 263   // check some internal logic
 264   VERIFY((  internal::has_nullary_operator<internal::scalar_constant_op<double> >::value ));
 265   VERIFY(( !internal::has_unary_operator<internal::scalar_constant_op<double> >::value ));
 266   VERIFY(( !internal::has_binary_operator<internal::scalar_constant_op<double> >::value ));
 267   VERIFY((  internal::functor_has_linear_access<internal::scalar_constant_op<double> >::ret ));
 268
 269   VERIFY(( !internal::has_nullary_operator<internal::scalar_identity_op<double> >::value ));
 270   VERIFY(( !internal::has_unary_operator<internal::scalar_identity_op<double> >::value ));
 271   VERIFY((  internal::has_binary_operator<internal::scalar_identity_op<double> >::value ));
 272   VERIFY(( !internal::functor_has_linear_access<internal::scalar_identity_op<double> >::ret ));
 273
 274   VERIFY(( !internal::has_nullary_operator<internal::linspaced_op<float,float> >::value ));
 275   VERIFY((  internal::has_unary_operator<internal::linspaced_op<float,float> >::value ));
 276   VERIFY(( !internal::has_binary_operator<internal::linspaced_op<float,float> >::value ));
 277   VERIFY((  internal::functor_has_linear_access<internal::linspaced_op<float,float> >::ret ));
 278
 279   // Regression unit test for a weird MSVC bug.
 280   // Search "nullary_wrapper_workaround_msvc" in CoreEvaluators.h for the details.
 281   // See also traits<Ref>::match.
 282   {
 283     MatrixXf A = MatrixXf::Random(3,3);
 284     Ref<const MatrixXf> R = 2.0*A;
 285     VERIFY_IS_APPROX(R, A+A);
 286
 287     Ref<const MatrixXf> R1 = MatrixXf::Random(3,3)+A;
 288
 289     VectorXi V = VectorXi::Random(3);
 290     Ref<const VectorXi> R2 = VectorXi::LinSpaced(3,1,3)+V;
 291     VERIFY_IS_APPROX(R2, V+Vector3i(1,2,3));
 292
 293     VERIFY((  internal::has_nullary_operator<internal::scalar_constant_op<float> >::value ));
 294     VERIFY(( !internal::has_unary_operator<internal::scalar_constant_op<float> >::value ));
 295     VERIFY(( !internal::has_binary_operator<internal::scalar_constant_op<float> >::value ));
 296     VERIFY((  internal::functor_has_linear_access<internal::scalar_constant_op<float> >::ret ));
 297
 298     VERIFY(( !internal::has_nullary_operator<internal::linspaced_op<int,int> >::value ));
 299     VERIFY((  internal::has_unary_operator<internal::linspaced_op<int,int> >::value ));
 300     VERIFY(( !internal::has_binary_operator<internal::linspaced_op<int,int> >::value ));
 301     VERIFY((  internal::functor_has_linear_access<internal::linspaced_op<int,int> >::ret ));
 302   }
 303 #endif
 304 }