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>
equalsIdentity(const MatrixType & A)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>
check_extremity_accuracy(const VectorType & v,const typename VectorType::Scalar & low,const typename VectorType::Scalar & high)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>
testVectorType(const VectorType & base)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>
testMatrixType(const MatrixType & m)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>
bug79()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>
bug1630()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>
nullary_overflow()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>
nullary_internal_logic()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
EIGEN_DECLARE_TEST(nullary)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 }
342