1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra. Eigen itself is part of the KDE project.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
|
---|
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 "sparse.h"
|
---|
11 |
|
---|
12 | template<typename SetterType,typename DenseType, typename Scalar, int Options>
|
---|
13 | bool test_random_setter(SparseMatrix<Scalar,Options>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords)
|
---|
14 | {
|
---|
15 | typedef SparseMatrix<Scalar,Options> SparseType;
|
---|
16 | {
|
---|
17 | sm.setZero();
|
---|
18 | SetterType w(sm);
|
---|
19 | std::vector<Vector2i> remaining = nonzeroCoords;
|
---|
20 | while(!remaining.empty())
|
---|
21 | {
|
---|
22 | int i = ei_random<int>(0,remaining.size()-1);
|
---|
23 | w(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y());
|
---|
24 | remaining[i] = remaining.back();
|
---|
25 | remaining.pop_back();
|
---|
26 | }
|
---|
27 | }
|
---|
28 | return sm.isApprox(ref);
|
---|
29 | }
|
---|
30 |
|
---|
31 | template<typename SetterType,typename DenseType, typename T>
|
---|
32 | bool test_random_setter(DynamicSparseMatrix<T>& sm, const DenseType& ref, const std::vector<Vector2i>& nonzeroCoords)
|
---|
33 | {
|
---|
34 | sm.setZero();
|
---|
35 | std::vector<Vector2i> remaining = nonzeroCoords;
|
---|
36 | while(!remaining.empty())
|
---|
37 | {
|
---|
38 | int i = ei_random<int>(0,remaining.size()-1);
|
---|
39 | sm.coeffRef(remaining[i].x(),remaining[i].y()) = ref.coeff(remaining[i].x(),remaining[i].y());
|
---|
40 | remaining[i] = remaining.back();
|
---|
41 | remaining.pop_back();
|
---|
42 | }
|
---|
43 | return sm.isApprox(ref);
|
---|
44 | }
|
---|
45 |
|
---|
46 | template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& ref)
|
---|
47 | {
|
---|
48 | const int rows = ref.rows();
|
---|
49 | const int cols = ref.cols();
|
---|
50 | typedef typename SparseMatrixType::Scalar Scalar;
|
---|
51 | enum { Flags = SparseMatrixType::Flags };
|
---|
52 |
|
---|
53 | double density = std::max(8./(rows*cols), 0.01);
|
---|
54 | typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
|
---|
55 | typedef Matrix<Scalar,Dynamic,1> DenseVector;
|
---|
56 | Scalar eps = 1e-6;
|
---|
57 |
|
---|
58 | SparseMatrixType m(rows, cols);
|
---|
59 | DenseMatrix refMat = DenseMatrix::Zero(rows, cols);
|
---|
60 | DenseVector vec1 = DenseVector::Random(rows);
|
---|
61 | Scalar s1 = ei_random<Scalar>();
|
---|
62 |
|
---|
63 | std::vector<Vector2i> zeroCoords;
|
---|
64 | std::vector<Vector2i> nonzeroCoords;
|
---|
65 | initSparse<Scalar>(density, refMat, m, 0, &zeroCoords, &nonzeroCoords);
|
---|
66 |
|
---|
67 | if (zeroCoords.size()==0 || nonzeroCoords.size()==0)
|
---|
68 | return;
|
---|
69 |
|
---|
70 | // test coeff and coeffRef
|
---|
71 | for (int i=0; i<(int)zeroCoords.size(); ++i)
|
---|
72 | {
|
---|
73 | VERIFY_IS_MUCH_SMALLER_THAN( m.coeff(zeroCoords[i].x(),zeroCoords[i].y()), eps );
|
---|
74 | if(ei_is_same_type<SparseMatrixType,SparseMatrix<Scalar,Flags> >::ret)
|
---|
75 | VERIFY_RAISES_ASSERT( m.coeffRef(zeroCoords[0].x(),zeroCoords[0].y()) = 5 );
|
---|
76 | }
|
---|
77 | VERIFY_IS_APPROX(m, refMat);
|
---|
78 |
|
---|
79 | m.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
|
---|
80 | refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
|
---|
81 |
|
---|
82 | VERIFY_IS_APPROX(m, refMat);
|
---|
83 | /*
|
---|
84 | // test InnerIterators and Block expressions
|
---|
85 | for (int t=0; t<10; ++t)
|
---|
86 | {
|
---|
87 | int j = ei_random<int>(0,cols-1);
|
---|
88 | int i = ei_random<int>(0,rows-1);
|
---|
89 | int w = ei_random<int>(1,cols-j-1);
|
---|
90 | int h = ei_random<int>(1,rows-i-1);
|
---|
91 |
|
---|
92 | // VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w));
|
---|
93 | for(int c=0; c<w; c++)
|
---|
94 | {
|
---|
95 | VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c));
|
---|
96 | for(int r=0; r<h; r++)
|
---|
97 | {
|
---|
98 | // VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r));
|
---|
99 | }
|
---|
100 | }
|
---|
101 | // for(int r=0; r<h; r++)
|
---|
102 | // {
|
---|
103 | // VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r));
|
---|
104 | // for(int c=0; c<w; c++)
|
---|
105 | // {
|
---|
106 | // VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c));
|
---|
107 | // }
|
---|
108 | // }
|
---|
109 | }
|
---|
110 |
|
---|
111 | for(int c=0; c<cols; c++)
|
---|
112 | {
|
---|
113 | VERIFY_IS_APPROX(m.col(c) + m.col(c), (m + m).col(c));
|
---|
114 | VERIFY_IS_APPROX(m.col(c) + m.col(c), refMat.col(c) + refMat.col(c));
|
---|
115 | }
|
---|
116 |
|
---|
117 | for(int r=0; r<rows; r++)
|
---|
118 | {
|
---|
119 | VERIFY_IS_APPROX(m.row(r) + m.row(r), (m + m).row(r));
|
---|
120 | VERIFY_IS_APPROX(m.row(r) + m.row(r), refMat.row(r) + refMat.row(r));
|
---|
121 | }
|
---|
122 | */
|
---|
123 |
|
---|
124 | // test SparseSetters
|
---|
125 | // coherent setter
|
---|
126 | // TODO extend the MatrixSetter
|
---|
127 | // {
|
---|
128 | // m.setZero();
|
---|
129 | // VERIFY_IS_NOT_APPROX(m, refMat);
|
---|
130 | // SparseSetter<SparseMatrixType, FullyCoherentAccessPattern> w(m);
|
---|
131 | // for (int i=0; i<nonzeroCoords.size(); ++i)
|
---|
132 | // {
|
---|
133 | // w->coeffRef(nonzeroCoords[i].x(),nonzeroCoords[i].y()) = refMat.coeff(nonzeroCoords[i].x(),nonzeroCoords[i].y());
|
---|
134 | // }
|
---|
135 | // }
|
---|
136 | // VERIFY_IS_APPROX(m, refMat);
|
---|
137 |
|
---|
138 | // random setter
|
---|
139 | // {
|
---|
140 | // m.setZero();
|
---|
141 | // VERIFY_IS_NOT_APPROX(m, refMat);
|
---|
142 | // SparseSetter<SparseMatrixType, RandomAccessPattern> w(m);
|
---|
143 | // std::vector<Vector2i> remaining = nonzeroCoords;
|
---|
144 | // while(!remaining.empty())
|
---|
145 | // {
|
---|
146 | // int i = ei_random<int>(0,remaining.size()-1);
|
---|
147 | // w->coeffRef(remaining[i].x(),remaining[i].y()) = refMat.coeff(remaining[i].x(),remaining[i].y());
|
---|
148 | // remaining[i] = remaining.back();
|
---|
149 | // remaining.pop_back();
|
---|
150 | // }
|
---|
151 | // }
|
---|
152 | // VERIFY_IS_APPROX(m, refMat);
|
---|
153 |
|
---|
154 | VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdMapTraits> >(m,refMat,nonzeroCoords) ));
|
---|
155 | #ifdef EIGEN_UNORDERED_MAP_SUPPORT
|
---|
156 | VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, StdUnorderedMapTraits> >(m,refMat,nonzeroCoords) ));
|
---|
157 | #endif
|
---|
158 | #ifdef _DENSE_HASH_MAP_H_
|
---|
159 | VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleDenseHashMapTraits> >(m,refMat,nonzeroCoords) ));
|
---|
160 | #endif
|
---|
161 | #ifdef _SPARSE_HASH_MAP_H_
|
---|
162 | VERIFY(( test_random_setter<RandomSetter<SparseMatrixType, GoogleSparseHashMapTraits> >(m,refMat,nonzeroCoords) ));
|
---|
163 | #endif
|
---|
164 |
|
---|
165 | // test fillrand
|
---|
166 | {
|
---|
167 | DenseMatrix m1(rows,cols);
|
---|
168 | m1.setZero();
|
---|
169 | SparseMatrixType m2(rows,cols);
|
---|
170 | m2.startFill();
|
---|
171 | for (int j=0; j<cols; ++j)
|
---|
172 | {
|
---|
173 | for (int k=0; k<rows/2; ++k)
|
---|
174 | {
|
---|
175 | int i = ei_random<int>(0,rows-1);
|
---|
176 | if (m1.coeff(i,j)==Scalar(0))
|
---|
177 | m2.fillrand(i,j) = m1(i,j) = ei_random<Scalar>();
|
---|
178 | }
|
---|
179 | }
|
---|
180 | m2.endFill();
|
---|
181 | VERIFY_IS_APPROX(m2,m1);
|
---|
182 | }
|
---|
183 |
|
---|
184 | // test RandomSetter
|
---|
185 | /*{
|
---|
186 | SparseMatrixType m1(rows,cols), m2(rows,cols);
|
---|
187 | DenseMatrix refM1 = DenseMatrix::Zero(rows, rows);
|
---|
188 | initSparse<Scalar>(density, refM1, m1);
|
---|
189 | {
|
---|
190 | Eigen::RandomSetter<SparseMatrixType > setter(m2);
|
---|
191 | for (int j=0; j<m1.outerSize(); ++j)
|
---|
192 | for (typename SparseMatrixType::InnerIterator i(m1,j); i; ++i)
|
---|
193 | setter(i.index(), j) = i.value();
|
---|
194 | }
|
---|
195 | VERIFY_IS_APPROX(m1, m2);
|
---|
196 | }*/
|
---|
197 | // std::cerr << m.transpose() << "\n\n" << refMat.transpose() << "\n\n";
|
---|
198 | // VERIFY_IS_APPROX(m, refMat);
|
---|
199 |
|
---|
200 | // test basic computations
|
---|
201 | {
|
---|
202 | DenseMatrix refM1 = DenseMatrix::Zero(rows, rows);
|
---|
203 | DenseMatrix refM2 = DenseMatrix::Zero(rows, rows);
|
---|
204 | DenseMatrix refM3 = DenseMatrix::Zero(rows, rows);
|
---|
205 | DenseMatrix refM4 = DenseMatrix::Zero(rows, rows);
|
---|
206 | SparseMatrixType m1(rows, rows);
|
---|
207 | SparseMatrixType m2(rows, rows);
|
---|
208 | SparseMatrixType m3(rows, rows);
|
---|
209 | SparseMatrixType m4(rows, rows);
|
---|
210 | initSparse<Scalar>(density, refM1, m1);
|
---|
211 | initSparse<Scalar>(density, refM2, m2);
|
---|
212 | initSparse<Scalar>(density, refM3, m3);
|
---|
213 | initSparse<Scalar>(density, refM4, m4);
|
---|
214 |
|
---|
215 | VERIFY_IS_APPROX(m1+m2, refM1+refM2);
|
---|
216 | VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3);
|
---|
217 | VERIFY_IS_APPROX(m3.cwise()*(m1+m2), refM3.cwise()*(refM1+refM2));
|
---|
218 | VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2);
|
---|
219 |
|
---|
220 | VERIFY_IS_APPROX(m1*=s1, refM1*=s1);
|
---|
221 | VERIFY_IS_APPROX(m1/=s1, refM1/=s1);
|
---|
222 |
|
---|
223 | VERIFY_IS_APPROX(m1+=m2, refM1+=refM2);
|
---|
224 | VERIFY_IS_APPROX(m1-=m2, refM1-=refM2);
|
---|
225 |
|
---|
226 | VERIFY_IS_APPROX(m1.col(0).eigen2_dot(refM2.row(0)), refM1.col(0).eigen2_dot(refM2.row(0)));
|
---|
227 |
|
---|
228 | refM4.setRandom();
|
---|
229 | // sparse cwise* dense
|
---|
230 | VERIFY_IS_APPROX(m3.cwise()*refM4, refM3.cwise()*refM4);
|
---|
231 | // VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4);
|
---|
232 | }
|
---|
233 |
|
---|
234 | // test innerVector()
|
---|
235 | {
|
---|
236 | DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
|
---|
237 | SparseMatrixType m2(rows, rows);
|
---|
238 | initSparse<Scalar>(density, refMat2, m2);
|
---|
239 | int j0 = ei_random(0,rows-1);
|
---|
240 | int j1 = ei_random(0,rows-1);
|
---|
241 | VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0));
|
---|
242 | VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1));
|
---|
243 | //m2.innerVector(j0) = 2*m2.innerVector(j1);
|
---|
244 | //refMat2.col(j0) = 2*refMat2.col(j1);
|
---|
245 | //VERIFY_IS_APPROX(m2, refMat2);
|
---|
246 | }
|
---|
247 |
|
---|
248 | // test innerVectors()
|
---|
249 | {
|
---|
250 | DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
|
---|
251 | SparseMatrixType m2(rows, rows);
|
---|
252 | initSparse<Scalar>(density, refMat2, m2);
|
---|
253 | int j0 = ei_random(0,rows-2);
|
---|
254 | int j1 = ei_random(0,rows-2);
|
---|
255 | int n0 = ei_random<int>(1,rows-std::max(j0,j1));
|
---|
256 | VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0));
|
---|
257 | VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0),
|
---|
258 | refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0));
|
---|
259 | //m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0);
|
---|
260 | //refMat2.block(0,j0,rows,n0) = refMat2.block(0,j0,rows,n0) + refMat2.block(0,j1,rows,n0);
|
---|
261 | }
|
---|
262 |
|
---|
263 | // test transpose
|
---|
264 | {
|
---|
265 | DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
|
---|
266 | SparseMatrixType m2(rows, rows);
|
---|
267 | initSparse<Scalar>(density, refMat2, m2);
|
---|
268 | VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval());
|
---|
269 | VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
|
---|
270 | }
|
---|
271 |
|
---|
272 | // test prune
|
---|
273 | {
|
---|
274 | SparseMatrixType m2(rows, rows);
|
---|
275 | DenseMatrix refM2(rows, rows);
|
---|
276 | refM2.setZero();
|
---|
277 | int countFalseNonZero = 0;
|
---|
278 | int countTrueNonZero = 0;
|
---|
279 | m2.startFill();
|
---|
280 | for (int j=0; j<m2.outerSize(); ++j)
|
---|
281 | for (int i=0; i<m2.innerSize(); ++i)
|
---|
282 | {
|
---|
283 | float x = ei_random<float>(0,1);
|
---|
284 | if (x<0.1)
|
---|
285 | {
|
---|
286 | // do nothing
|
---|
287 | }
|
---|
288 | else if (x<0.5)
|
---|
289 | {
|
---|
290 | countFalseNonZero++;
|
---|
291 | m2.fill(i,j) = Scalar(0);
|
---|
292 | }
|
---|
293 | else
|
---|
294 | {
|
---|
295 | countTrueNonZero++;
|
---|
296 | m2.fill(i,j) = refM2(i,j) = Scalar(1);
|
---|
297 | }
|
---|
298 | }
|
---|
299 | m2.endFill();
|
---|
300 | VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros());
|
---|
301 | VERIFY_IS_APPROX(m2, refM2);
|
---|
302 | m2.prune(1);
|
---|
303 | VERIFY(countTrueNonZero==m2.nonZeros());
|
---|
304 | VERIFY_IS_APPROX(m2, refM2);
|
---|
305 | }
|
---|
306 | }
|
---|
307 |
|
---|
308 | void test_eigen2_sparse_basic()
|
---|
309 | {
|
---|
310 | for(int i = 0; i < g_repeat; i++) {
|
---|
311 | CALL_SUBTEST_1( sparse_basic(SparseMatrix<double>(8, 8)) );
|
---|
312 | CALL_SUBTEST_2( sparse_basic(SparseMatrix<std::complex<double> >(16, 16)) );
|
---|
313 | CALL_SUBTEST_1( sparse_basic(SparseMatrix<double>(33, 33)) );
|
---|
314 |
|
---|
315 | CALL_SUBTEST_3( sparse_basic(DynamicSparseMatrix<double>(8, 8)) );
|
---|
316 | }
|
---|
317 | }
|
---|