DGtal 2.1.0
Loading...
Searching...
No Matches
testSimpleMatrix.cpp
Go to the documentation of this file.
1
31#include <iostream>
32#include "DGtal/base/Common.h"
33#include "DGtal/math/linalg/SimpleMatrix.h"
34#include "DGtal/math/linalg/CStaticMatrix.h"
35#include "DGtal/math/linalg/CDenseMatrix.h"
36#include "DGtal/math/linalg/CStaticVector.h"
37#include "DGtal/math/linalg/CDenseVector.h"
38#include "DGtal/math/linalg/CLinearAlgebra.h"
39#include "DGtal/math/linalg/IntegerMatrixFunctions.h"
40#include "DGtal/helpers/StdDefs.h"
42
43using namespace std;
44using namespace DGtal;
45
47// Functions for testing class SimpleMatrix.
49
53bool testSimpleMatrix()
54{
55 unsigned int nbok = 0;
56 unsigned int nb = 0;
57
58 trace.beginBlock ( "Testing create ..." );
59
60 typedef SimpleMatrix<double,3,4> M34d;
61
62 M34d m34d;
63 trace.info() << m34d<<std::endl;
64
65 m34d.setComponent(1,2, 0.5);
66 trace.info() << m34d<<std::endl;
67
68 nbok += (m34d(1,2) == 0.5) ? 1 : 0;
69 nb++;
70 trace.info() << "(" << nbok << "/" << nb << ") "
71 << "true == true" << std::endl;
72
73 M34d matrix;
74 bool res=true;
75
76 matrix.constant(12.3);
77 trace.info() << matrix;
78 for(DGtal::Dimension i = 0; i< 3; ++i)
79 for(DGtal::Dimension j = 0; j< 4; ++j)
80 res = res && (matrix(i,j) == 12.3);
81 nbok += res ? 1 : 0;
82 nb++;
83 trace.info() << "(" << nbok << "/" << nb << ") "
84 << "all equals to 12.3" << std::endl;
85
86
87 trace.endBlock();
88
89 return nbok == nb;
90}
91
92bool testArithm()
93{
94 unsigned int nbok = 0;
95 unsigned int nb = 0;
96
97
98 typedef SimpleMatrix<double,3,4> M34d;
99 typedef SimpleMatrix<double,4,3> M43d;
100 typedef SimpleMatrix<double,3,3> M33d;
101
102 M34d m34d, two,four;
103 M34d m34dbis, resadd, ressub;
104
105 two.constant(2);
106 four.constant(4);
107
108 for(DGtal::Dimension i = 0; i< 3; ++i)
109 for(DGtal::Dimension j = 0; j< 4; ++j)
110 {
111 m34d.setComponent(i,j,i*j);
112 m34dbis.setComponent(i,j,i+j);
113 resadd.setComponent(i,j,i*j+i+j);
114 ressub.setComponent(i,j,i*j-(double)i-(double)j);
115 }
116
117
118 trace.info() << m34d <<std::endl;
119 trace.info() << m34dbis<<std::endl;
120
121 trace.beginBlock ( "Testing add ..." );
122 nbok += ((m34d + m34dbis) == resadd) ? 1 : 0;
123 nb++;
124 trace.info() << "(" << nbok << "/" << nb << ") "
125 << "ok" << std::endl;
126 nbok += ((m34dbis + m34d) == resadd) ? 1 : 0;
127 nb++;
128 trace.info() << "(" << nbok << "/" << nb << ") "
129 << "ok commutative" << std::endl;
130
131 M34d other;
132 other += m34d;
133 nbok += (other == m34d) ? 1 : 0;
134 nb++;
135 trace.info() << "(" << nbok << "/" << nb << ") "
136 << "ok +=" << std::endl;
137
138 trace.endBlock();
139
140 trace.beginBlock ( "Testing substraction ..." );
141 nbok += ((m34d - m34dbis) == ressub) ? 1 : 0;
142 nb++;
143 trace.info()<<ressub<<std::endl;
144 trace.info()<<m34d - m34dbis<<std::endl;
145
146 trace.info() << "(" << nbok << "/" << nb << ") "
147 << "ok simple" << std::endl;
148 trace.endBlock();
149
150 trace.beginBlock ( "Testing scalar product/divide ..." );
151 nbok += ( (two*2.0) == four) ? 1 : 0;
152 nb++;
153 trace.info()<<ressub<<std::endl;
154 trace.info() << "(" << nbok << "/" << nb << ") "
155 << " [2]*2 == [4]" << std::endl;
156
157 nbok += ( two == four/2.0) ? 1 : 0;
158 nb++;
159 trace.info()<<ressub<<std::endl;
160 trace.info() << "(" << nbok << "/" << nb << ") "
161 << " [2]= [4]/2" << std::endl;
162 trace.endBlock();
163
164
165 trace.beginBlock ( "Testing transpose ..." );
166 M43d transp = m34d.transpose();
167 nbok += (transp.transpose() == m34d) ? 1 : 0;
168 nb++;
169 trace.info() << "(" << nbok << "/" << nb << ") "
170 << "ok idem potent" << std::endl;
171 trace.endBlock();
172
173 trace.beginBlock ( "Testing product ..." );
174
175 M43d one;
176 M33d eight33;
177
178 one.constant(1);
179 eight33.constant(8);
180 trace.info() << two * one<<std::endl;
181 nbok += (two * one == eight33) ? 1 : 0;
182 nb++;
183 trace.info() << "(" << nbok << "/" << nb << ") "
184 << " [2]*[1] = [8]" << std::endl;
185 trace.endBlock();
186
187
188
189 return nbok == nb;
190
191}
192
193bool testColRow()
194{
195 unsigned int nbok = 0;
196 unsigned int nb = 0;
197
198 DGtal::SimpleMatrix<double,3,4> mat;
199 for(DGtal::Dimension i = 0; i< 3; ++i)
200 for(DGtal::Dimension j = 0; j< 4; ++j)
201 mat.setComponent(i,j,i+j);
202
203 trace.beginBlock("Get Row");
204 trace.info() << mat <<std::endl;
205 DGtal::SimpleMatrix<double,3,4>::RowVector row;
206 row = mat.row(1);
207 trace.info() << row << std::endl;
208 nbok += (row[1] == 2 ) ? 1 : 0;
209 nb++;
210 trace.info() << "(" << nbok << "/" << nb << ") "
211 << " row value" << std::endl;
212 trace.endBlock();
213
214 trace.beginBlock("Get Col");
215 DGtal::SimpleMatrix<double,3,4>::ColumnVector col;
216 col = mat.column(1);
217 trace.info() << row << std::endl;
218 nbok += (col[1] == 2 ) ? 1 : 0;
219 nb++;
220 trace.info() << "(" << nbok << "/" << nb << ") "
221 << " col value" << std::endl;
222 trace.endBlock();
223
224
225
226 trace.beginBlock("Prod Matrix x Row^t");
227 //Row vector is a dim 4 vector
228 DGtal::SimpleMatrix<double,3,4>::RowVector r(1,2,3,4);
229 DGtal::SimpleMatrix<double,3,4>::ColumnVector c;
230 c = mat*r;
231 DGtal::SimpleMatrix<double,3,4>::ColumnVector expected(20,30,40);
232
233 trace.info() << c << std::endl;
234 nbok += (c == expected) ? 1 : 0;
235 nb++;
236 trace.info() << "(" << nbok << "/" << nb << ") "
237 << " mat*row^t" << std::endl;
238 trace.endBlock();
239
240 return nbok == nb;
241}
242
243bool testDetCofactor()
244{
245 unsigned int nbok = 0;
246 unsigned int nb = 0;
247
248 typedef DGtal::SimpleMatrix<double,2,2> MAT2;
249 MAT2 mat2;
250 mat2.setComponent(0,0,1);
251 mat2.setComponent(1,1,2);
252
253 trace.beginBlock("det2x2 tests...");
254 trace.info() << mat2<<std::endl;
255 trace.info() << mat2.determinant() << std::endl;
256 nbok += (mat2.determinant() == 2) ? 1 : 0;
257 nb++;
258 trace.info() << "(" << nbok << "/" << nb << ") "
259 << " 2" << std::endl;
260 trace.endBlock();
261
262 typedef DGtal::SimpleMatrix<double,3,3> MAT;
263 MAT mat;
264 mat.setComponent(0,0,1);
265 mat.setComponent(1,1,2);
266 mat.setComponent(2,2,4);
267
268 trace.beginBlock("det3x3 tests...");
269 trace.info() << mat<<std::endl;
270 nbok += (mat.determinant() == 8) ? 1 : 0;
271 nb++;
272 trace.info() << "(" << nbok << "/" << nb << ") "
273 << " 8" << std::endl;
274 trace.endBlock();
275
276
277 typedef DGtal::SimpleMatrix<double,4,4> MAT44;
278 MAT44 mat44;
279 mat44.setComponent(0,0,1);
280 mat44.setComponent(1,1,2);
281 mat44.setComponent(2,2,4);
282 mat44.setComponent(3,3,4);
283
284 trace.beginBlock("det4x4 tests...");
285 trace.info() << mat44 <<std::endl;
286 trace.info() << mat44.determinant() << std::endl;
287 nbok += (mat44.determinant() == 32) ? 1 : 0;
288 nb++;
289 trace.info() << "(" << nbok << "/" << nb << ") "
290 << " 32" << std::endl;
291 trace.endBlock();
292
293
294 return nbok == nb;
295}
296
297bool testM1Matrix()
298{
299 trace.beginBlock("Mx1 matrix test");
300 SimpleMatrix<double, 3,1> mat;
301 trace.info() << mat<<std::endl;
302 trace.endBlock();
303 return true;
304}
305
306bool testInverse()
307{
308 unsigned int nbok = 0;
309 unsigned int nb = 0;
310
311 trace.beginBlock("Inverse tests 2X2...");
312
313 typedef DGtal::SimpleMatrix<double,2,2> MAT2;
314 MAT2 mat2;
315 mat2.setComponent(0,0,1);
316 mat2.setComponent(1,1,2);
317
318 MAT2 Id2;
319 Id2.identity();
320
321 trace.info() << mat2<<std::endl;
322 trace.info() << mat2.inverse() << std::endl;
323 nbok += (( mat2 * mat2.inverse() )== Id2 ) ? 1 : 0;
324 nb++;
325 trace.info() << "(" << nbok << "/" << nb << ") "
326 << " M*M^-1=Id" << std::endl;
327
328 trace.endBlock();
329
330 trace.beginBlock("Inverse tests 6x6 random...");
331
332 typedef DGtal::SimpleMatrix<double,6,6> MAT6;
333 MAT6 mat;
334
335 for(unsigned int i=0; i< 6; i++)
336 for(unsigned int j=0; j< 6; j++)
337 mat.setComponent(i,j, rand() % 10);
338
339 MAT6 Id6;
340 Id6.identity();
341
342 trace.info() << "M= "<<mat<<std::endl;
343 trace.info() << "M^-1=" <<mat.inverse() << std::endl;
344 trace.info() << "det(M)= "<<mat.determinant() <<std::endl;
345 trace.info() << "M*M^-1= "<<mat.inverse()*mat << std::endl;
346
347
348
349 trace.endBlock();
350
351
352 return nbok == nb;
353}
354
355bool testConstructor()
356{
357 unsigned int nbok = 0;
358 unsigned int nb = 0;
359
360 trace.beginBlock( "Initilizer-list constructor test" );
361 SimpleMatrix<double, 3, 3> mat = {1, 2, 3, 4, 5, 6, 7, 8, 9};
362 trace.info() << mat << std::endl;
363
364 trace.info() << "Testing values: ";
365 trace.info() << mat( 0, 0 );
366 nbok += ( mat( 0, 0 ) == 1 ) ? 1 : 0;
367 nb++;
368 trace.info() << "(" << nbok << "/" << nb << ") ";
369
370 trace.info() << mat( 0, 1 );
371 nbok += ( mat( 0, 1 ) == 2 ) ? 1 : 0;
372 nb++;
373 trace.info() << "(" << nbok << "/" << nb << ") ";
374
375 trace.info() << mat( 2, 2 );
376 nbok += ( mat( 2, 2 ) == 9 ) ? 1 : 0;
377 nb++;
378 trace.info() << "(" << nbok << "/" << nb << ") ";
379
380 trace.info() << std::endl;
381 trace.endBlock();
382 return nbok == nb;
383}
384
385bool testConcepts()
386{
387 typedef DGtal::SimpleMatrix<double,3,3> Matrix;
388 typedef Matrix::ColumnVector Vector;
389
390 BOOST_CONCEPT_ASSERT(( concepts::CStaticVector<Vector> ));
391 BOOST_CONCEPT_ASSERT(( concepts::CDenseVector<Vector> ));
392 BOOST_CONCEPT_ASSERT(( concepts::CStaticMatrix<Matrix> ));
393 BOOST_CONCEPT_ASSERT(( concepts::CDenseMatrix<Matrix> ));
394 BOOST_CONCEPT_ASSERT(( concepts::CLinearAlgebra<Vector, Matrix> ));
395 BOOST_CONCEPT_ASSERT(( concepts::CStaticVector<DGtal::Z3i::Space::Vector> ));
396 BOOST_CONCEPT_ASSERT(( concepts::CDenseVector<DGtal::Z3i::Space::Vector> ));
397
398 return true;
399}
400
401bool testBareissDeterminant()
402{
403 unsigned int nbok = 0;
404 unsigned int nb = 0;
405
406 trace.beginBlock( "Bareiss determinant test" );
407 {
408 typedef DGtal::SimpleMatrix<int,3,3> Matrix;
409
410 Matrix M = { 1, 2, 3, 4, 5, 7, 6, 8, 9 };
411 auto d = M.determinant();
412 int db;
413 DGtal::functions::getDeterminantBareiss( db, M );
414 trace.info() << "d=" << d << " db=" << db << "\n";
415 nbok += ( d == 7 ) ? 1 : 0;
416 nb++;
417 nbok += ( db == 7 ) ? 1 : 0;
418 nb++;
419 }
420
421 {
422 typedef DGtal::SimpleMatrix<int,4,4> Matrix;
423 Matrix M = { 1, 2, 3, -4, 13, 4, 5, 7, 6, 8, 17, 9, 21, 12, -5, 11 };
424 auto d = M.determinant();
425 int64_t db, dbv;
426 DGtal::functions::getDeterminantBareiss( db, M );
427 auto MV = DGtal::functions::matrixAsVectorVector( M );
428 DGtal::functions::getDeterminantBareiss( dbv, MV );
429 trace.info() << "d=" << d << " db=" << db << " dbv=" << dbv << "\n";
430 nbok += ( d == -12260 ) ? 1 : 0;
431 nb++;
432 nbok += ( db == -12260 ) ? 1 : 0;
433 nb++;
434 nbok += ( dbv == -12260 ) ? 1 : 0;
435 nb++;
436 }
437 {
438 std::vector<int> V = { 1, 2, 3, -4, 13, 4, 5, 7, 6, 8, 17, 9, 21, 12, -5, 11 };
439 auto MV = DGtal::functions::matrixAsVectorVector( 4, 4, V );
440 int64_t db;
441 DGtal::functions::getDeterminantBareiss( db, MV );
442 trace.info() << "db=" << db << "\n";
443 nbok += ( db == -12260 ) ? 1 : 0;
444 nb++;
445 }
446 {
447 typedef DGtal::SimpleMatrix<int,5,5> Matrix;
448 Matrix M = { 1311, 1, 2, 3, -4,
449 13, 457, 4, 5, 7,
450 6, 8, -535, 17, 9,
451 21, 12, -5, 243, 11,
452 123,-39,411,630,23 };
453 auto d = M.determinant();
454 int64_t db;
455 BigInteger big_db;
456 DGtal::functions::getDeterminantBareiss( db, M );
457 DGtal::functions::getDeterminantBareiss( big_db, M );
458 int64_t cdb = NumberTraits<BigInteger>::castToInt64_t( big_db );
459 trace.info() << "d=" << d << " (i64)db=" << db << " (big)db=" << cdb << "\n";
460 // The line below raises a warning in the macos compiler.
461 // nbok += ( int64_t(d) != -171492636038LL ) ? 1 : 0; // int overflow
462 nbok += ( int(d) != (-171492636038LL % 2147483648LL ) ) ? 1 : 0; // int overflow
463 nb++;
464 nbok += ( db != -171492636038LL ) ? 1 : 0; // int64 overflow (intermediate computation)
465 nb++;
466 nbok += ( cdb == -171492636038LL ) ? 1 : 0;
467 nb++;
468 }
469
470 {
471 typedef DGtal::SimpleMatrix<double,4,4> Matrix;
472 Matrix M = { 1.5, 2.2, 3.1, -4.6, 13.3, 4.2, 5.7, 7.3, 6.4, 8.0, 17.9, 9.3, 21.2, 12.2, -5.1, 11.8 };
473 auto d = M.determinant();
474 double db;
475 DGtal::functions::getDeterminantBareiss( db, M );
476 trace.info() << "d=" << d << " db=" << db << "\n";
477 nbok += ( std::fabs( d - db ) < 1e-10 ) ? 1 : 0;
478 nb++;
479 }
480
481 std::cout << "(" << nbok << "/" << nb << ")\n";
482 trace.endBlock();
483
484 return nbok == nb;
485}
486
487template <typename Number>
488std::ostream&
489operator<<( std::ostream& out, const std::vector< Number>& v )
490{
491 for ( auto i = 0; i < v.size(); i++ )
492 out << v[ i ] << " ";
493 out << "\n";
494 return out;
495}
496
497bool testLLL()
498{
499 unsigned int nbok = 0;
500 unsigned int nb = 0;
501
502 typedef int64_t Integer;
503 trace.beginBlock( "Test LLL-reduction on matrices." );
504 {
505 std::vector< std::vector< Integer > > B = {
506 {73,-127,63},
507 {99,-12,-14},
508 {17,-26,14}
509 };
510 auto L1 = DGtal::functions::computeLLLBasis( B, 0.75 );
511 auto L2 = DGtal::functions::computeLLLBasis( B, 0.99 );
512 std::cout << "Init base: \n" << B
513 << "\nLLL base: delta=0.75\n" << L1 << "\n"
514 << "\nLLL base: delta=0.99\n" << L2 << "\n";
515
516 std::vector< std::vector< Integer > > R = { {-12,3,-7}, {32,-17,-49}, {-7,-20,0} };
517 nbok += L1 == R ? 1 : 0;
518 nb++;
519 nbok += L2 == R ? 1 : 0;
520 nb++;
521 std::cout << "(" << nbok << "/" << nb << ") " << (nbok == nb ? "PASSED\n" : "ERROR\n");
522 }
523 {
524 std::vector< std::vector< Integer > > B = {
525 {118, 113, 90, 78},
526 {187, 229, 12, 109},
527 { 26, 163, 223, 21}
528 };
529 auto L1 = DGtal::functions::computeLLLBasis( B, 0.75 );
530 auto L2 = DGtal::functions::computeLLLBasis( B, 0.99 );
531 std::cout << "Init base: \n" << B
532 << "\nLLL base: delta=0.75\n" << L1 << "\n"
533 << "\nLLL base: delta=0.99\n" << L2 << "\n";
534
535 std::vector< std::vector< Integer > > R1 = {
536 { 69, 116, -78, 31},
537 {-23, 166, 55, -26},
538 { 49, -3, 168, 47}
539 };
540 std::vector< std::vector< Integer > > R2 = {
541 { 69, 116, -78, 31},
542 {-23, 166, 55, -26},
543 { 49, -3, 168, 47}
544 };
545 nbok += L1 == R1 ? 1 : 0;
546 nb++;
547 nbok += L2 == R2 ? 1 : 0;
548 nb++;
549 std::cout << "(" << nbok << "/" << nb << ") " << (nbok == nb ? "PASSED\n" : "ERROR\n");
550 }
551 {
552 std::vector< std::vector< Integer > > B = {
553 {16,16,11,13,31,25,21, 9},
554 {31, 3, 2,18,11,31,30, 2},
555 { 4,25,11,18, 8,20,13,10},
556 {16,21, 9, 2,10,23, 7,27},
557 {28,12,26, 1, 2,18, 4,19}
558 };
559 auto L1 = DGtal::functions::computeLLLBasis( B, 0.75 );
560 auto L2 = DGtal::functions::computeLLLBasis( B, 0.99 );
561 std::cout << "Init base: \n" << B
562 << "\nLLL base: delta=0.75\n" << L1 << "\n"
563 << "\nLLL base: delta=0.99\n" << L2 << "\n";
564
565 std::vector< std::vector< Integer > > R1 = {
566 { 15,-13, -9, 5,-20, 6, 9, -7},
567 {-12, 9, 0, 5,-23, -5, -8, 1},
568 { 12, -4, -2,-16, 2, 3, -6, 17},
569 { 12, -9, 17, -1, -8, -5, -3, -8},
570 { 4, 25, 11, 18, 8, 20, 13, 10}
571 };
572 std::vector< std::vector< Integer > > R2 = {
573 { 12, -9, 17, -1, -8, -5, -3, -8},
574 { 12, -4, -2,-16, 2, 3, -6, 17},
575 {-12, 9, 0, 5,-23, -5, -8, 1},
576 { 15,-13, -9, 5,-20, 6, 9, -7},
577 { 4, 25, 11, 18, 8, 20, 13, 10}
578 };
579 nbok += L1 == R1 ? 1 : 0;
580 nb++;
581 nbok += L2 == R2 ? 1 : 0;
582 nb++;
583 std::cout << "(" << nbok << "/" << nb << ") " << (nbok == nb ? "PASSED\n" : "ERROR\n");
584 }
585 {
586 // when the matrix is unimodular, outputs canonic vectors.
587 std::vector< std::vector< Integer > > B = {
588 { -3, 10, 47, 61, -53, -126, 713, 601,-1476, 1569},
589 { 2, -7, -33, -43, 37, 89, -502, -425, 1047,-1103},
590 { -3, 11, 53, 69, -59, -142, 800, 677,-1663, 1764},
591 { 1, -9, -48, -63, 52, 130, -727, -623, 1543,-1604},
592 { -2, 9, 48, 63, -49, -124, 680, 583,-1409, 1533},
593 { 5, -25, -118, -163, 113, 334,-1761,-1595, 4030,-3838},
594 { -3, 17, 85, 118, -84, -245, 1297, 1173,-2974, 2824},
595 { 5, -24, -119, -156, 126, 321,-1782,-1534, 3799,-3921},
596 { 2, -10, -44, -65, 42, 137, -699, -659, 1713,-1489},
597 { 1, -5, -23, -27, 33, 64, -405, -315, 784, -868}
598 };
599 auto L2 = DGtal::functions::computeLLLBasis( B, 0.99 );
600 Integer d, d2;
601 DGtal::functions::getDeterminantBareiss( d, B );
602 DGtal::functions::getDeterminantBareiss( d2, L2 );
603 std::cout << "Init base: \n" << B
604 << "\nLLL base: delta=0.99\n" << L2 << "\n"
605 << "det(B)=" << d << " det(L)=" << d2 << "\n";
606 nbok += d == 1 ? 1 : 0;
607 nb++;
608 nbok += d2 == 1 ? 1 : 0;
609 nb++;
610 std::cout << "(" << nbok << "/" << nb << ") "
611 << "When input matrix is unimodular, output matrix L is unimodular: "
612 << (nbok == nb ? "PASSED\n" : "ERROR\n");
613 for ( auto i = 0; i < B.size(); i++ )
614 {
615 nbok += DGtal::functions::normL1( L2[ i ] ) == 1 ? 1 : 0;
616 nb++;
617 }
618 std::cout << "(" << nbok << "/" << nb << ") "
619 << "The output matrix is then canonic: "
620 << (nbok == nb ? "PASSED\n" : "ERROR\n");
621 trace.beginBlock( "Shorten B" );
622 auto nbs = functions::shortenBasis( B );
623 std::cout << "Shorten base: \n" << B << "#nb shortening=" << nbs << "\n";
624 trace.endBlock();
625 }
626 trace.endBlock();
627 return nbok == nb;
628}
629
630bool testOrthogonalLattice()
631{
632 unsigned int nbok = 0;
633 unsigned int nb = 0;
634
635 typedef int64_t Integer;
636 trace.beginBlock( "Test orthogonal lattice computation." );
637 for ( auto i = 0; i < 10000; i++ )
638 {
639 vector<int64_t> n = { rand() % 30 - 15, rand() % 30 - 15, rand() % 30 - 15 };
640
641 auto g = functions::makePrimitive( n );
642 if ( g==0 ) continue;
643 vector<int64_t> no = n;
644 functions::negate( no );
645 auto L = DGtal::functions::computeOrthogonalLattice( n );
646 functions::shortenBasis( L );
647 Integer l0 = functions::dotProduct( L[ 0 ], n ); // vectors are orthogonal
648 Integer l1 = functions::dotProduct( L[ 1 ], n );
649 auto c = functions::crossProduct( L[ 0 ], L[ 1 ] ); // u x v = n
650 nbok += l0 == 0 ? 1 : 0;
651 nbok += l1 == 0 ? 1 : 0;
652 nbok += ( functions::equals( c, n ) || functions::equals( c, no ) ) ? 1 : 0;
653 nb += 3;
654 if ( nbok != nb )
655 {
656 std::cout << "----------- " << nbok << "/" << nb << " ----------------\n";
657 std::cout << "Error for vector n=" << n;
658 std::cout << "u.n=" << l0 << " v.n=" << l1 << " uxv=" << c << "\n";
659 std::cout << "u=" << L[0] << "v=" << L[1];
660 std::cout << "u.n=" << l0 << " v.n=" << l1 << " uxv=" << c << "\n";
661 std::cout << "------------------------------------\n";
662 break;
663 }
664 }
665 trace.info() << "(" << nbok << "/" << nb << ") "
666 << "Orthogonal lattice in 3D: "
667 << (nbok == nb ? "PASSED\n" : "ERROR\n");
668 for ( auto i = 0; i < 10000; i++ )
669 {
670 vector<int64_t> n = { rand() % 30 - 15, rand() % 30 - 15,
671 rand() % 30 - 15, rand() % 30 - 15 };
672 auto g = functions::makePrimitive( n );
673 if ( g==0 ) continue;
674 vector<int64_t> no = n;
675 functions::negate( no );
676 auto L = DGtal::functions::computeOrthogonalLattice( n );
677 functions::shortenBasis( L );
678 Integer l0 = functions::dotProduct( L[ 0 ], n );
679 Integer l1 = functions::dotProduct( L[ 1 ], n );
680 Integer l2 = functions::dotProduct( L[ 2 ], n );
681 Integer n0 = functions::dotProduct( L[ 0 ], L[ 0 ] );
682 Integer n1 = functions::dotProduct( L[ 1 ], L[ 1 ] );
683 Integer n2 = functions::dotProduct( L[ 2 ], L[ 2 ] );
684 nbok += l0 == 0 ? 1 : 0; // vectors are orthogonal
685 nbok += l1 == 0 ? 1 : 0;
686 nbok += l2 == 0 ? 1 : 0;
687 nbok += n0 > 0 ? 1 : 0; // vectors are non null
688 nbok += n1 > 0 ? 1 : 0;
689 nbok += n2 > 0 ? 1 : 0;
690 nb += 6;
691 if ( nbok != nb )
692 {
693 std::cout << "----------- " << nbok << "/" << nb << " ----------------\n";
694 std::cout << "Error for vector n=" << n;
695 std::cout << "u=" << L[0] << "v=" << L[1] << "w=" << L[2];
696 std::cout << "u.n=" << l0 << " v.n=" << l1 << " w.n=" << l2 << "\n";
697 std::cout << "u.u=" << n0 << " v.v=" << n1 << " w.w=" << n2 << "\n";
698 std::cout << "------------------------------------\n";
699 break;
700 }
701 }
702 trace.info() << "(" << nbok << "/" << nb << ") "
703 << "Orthogonal lattice in 4D: "
704 << (nbok == nb ? "PASSED\n" : "ERROR\n");
705 trace.endBlock();
706 return nbok == nb;
707}
708
710// Standard services - public :
711
712int main( int argc, char** argv )
713{
714 trace.beginBlock ( "Testing class SimpleMatrix" );
715 trace.info() << "Args:";
716 for ( int i = 0; i < argc; ++i )
717 trace.info() << " " << argv[ i ];
718 trace.info() << endl;
719
720 bool res = testSimpleMatrix() && testArithm() && testColRow()
721 && testDetCofactor() && testM1Matrix() && testInverse()
722 && testConcepts() && testConstructor()
723 && testBareissDeterminant()
724 && testLLL()
725 && testOrthogonalLattice();
726 trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
727 trace.endBlock();
728 return res ? 0 : 1;
729}
730// //
DGtal::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition PointVector.h:593
DGtal::SimpleMatrix
Aim: implements basic MxN Matrix services (M,N>=1).
Definition SimpleMatrix.h:76
DGtal::SimpleMatrix::determinant
Component determinant() const
DGtal::SimpleMatrix::column
ColumnVector column(const DGtal::Dimension j) const
DGtal::SimpleMatrix::row
RowVector row(const DGtal::Dimension i) const
DGtal::SimpleMatrix::setComponent
void setComponent(const DGtal::Dimension i, const DGtal::Dimension j, const Component &aValue)
DGtal::SimpleMatrix::constant
void constant(const Component &aScalar)
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Trace::emphase
std::ostream & emphase()
DGtal::Trace::info
std::ostream & info()
DGtal::Trace::endBlock
double endBlock()
Vector
DigitalPlane::Point Vector
Definition examplePlaneProbingParallelepipedEstimator.cpp:45
DGtal::functions::matrixAsVectorVector
std::vector< std::vector< TComponent > > matrixAsVectorVector(std::size_t m, std::size_t n, const std::vector< TComponent > &c)
DGtal::functions::shortenBasis
std::size_t shortenBasis(std::vector< std::vector< TComponent > > &B)
DGtal::functions::equals
bool equals(const std::vector< TComponent > &a, const std::vector< TComponent > &b)
DGtal::functions::crossProduct
std::vector< typename DGtal::ArithmeticConversionTraits< T, U >::type > crossProduct(const std::vector< T > &a, const std::vector< U > &b)
DGtal::functions::dotProduct
DGtal::ArithmeticConversionTraits< T, U >::type dotProduct(const std::vector< T > &a, const std::vector< U > &b)
DGtal::functions::computeLLLBasis
std::vector< std::vector< TComponent > > computeLLLBasis(const std::vector< std::vector< TComponent > > &B, TDouble delta=0.75)
DGtal::functions::makePrimitive
TComponent makePrimitive(std::vector< TComponent > &N)
DGtal::functions::computeOrthogonalLattice
std::vector< std::vector< TComponent > > computeOrthogonalLattice(std::vector< TComponent > N)
DGtal::functions::negate
void negate(std::vector< TComponent > &V)
DGtal::functions::getDeterminantBareiss
void getDeterminantBareiss(TInternalNumber &result, const SimpleMatrix< TComponent, TN, TN > &matrix)
DGtal::functions::normL1
T normL1(const std::vector< T > &a)
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
Definition ClosedIntegerHalfPlane.h:49
DGtal::int64_t
std::int64_t int64_t
signed 94-bit integer.
Definition BasicTypes.h:73
DGtal::operator<<
std::ostream & operator<<(std::ostream &out, const ClosedIntegerHalfPlane< TSpace > &object)
DGtal::Dimension
DGtal::uint32_t Dimension
Definition Common.h:119
DGtal::trace
Trace trace
DGtal::BigInteger
boost::multiprecision::number< boost::multiprecision::cpp_int_backend<>, boost::multiprecision::et_off > BigInteger
Definition BasicTypes.h:75
std
STL namespace.
DGtal::NumberTraits
Aim: The traits class for all models of Cinteger.
Definition NumberTraits.h:562
DGtal::concepts::CDenseMatrix
Aim: Represent any dynamic or static sized matrix having dense representation.
Definition CDenseMatrix.h:89
DGtal::concepts::CDenseVector
Aim: Represent any dynamic or static sized matrix having dense representation.
Definition CDenseVector.h:89
DGtal::concepts::CLinearAlgebra
Aim: Check right multiplication between matrix and vector and internal matrix multiplication....
Definition CLinearAlgebra.h:94
DGtal::concepts::CStaticMatrix
Aim: Represent any static sized matrix having sparse or dense representation.
Definition CStaticMatrix.h:90
DGtal::concepts::CStaticVector
Aim: Represent any static sized column vector having sparse or dense representation.
Definition CStaticVector.h:91
g
std::mt19937 g(rd())
main
int main()
Definition testBits.cpp:56
testInverse
bool testInverse()
Definition testSimpleMatrix.cpp:306
testBareissDeterminant
bool testBareissDeterminant()
Definition testSimpleMatrix.cpp:401
testM1Matrix
bool testM1Matrix()
Definition testSimpleMatrix.cpp:297
testDetCofactor
bool testDetCofactor()
Definition testSimpleMatrix.cpp:243
testColRow
bool testColRow()
Definition testSimpleMatrix.cpp:193
testOrthogonalLattice
bool testOrthogonalLattice()
Definition testSimpleMatrix.cpp:630
testLLL
bool testLLL()
Definition testSimpleMatrix.cpp:497
testConstructor
bool testConstructor()
Definition testSimpleMatrix.cpp:355
testConcepts
bool testConcepts()
Definition testSimpleMatrix.cpp:385
testArithm
bool testArithm()
Definition testSimpleMatrix.cpp:92
testSimpleMatrix
bool testSimpleMatrix()
Definition testSimpleMatrix.cpp:53