49 #include <visp3/core/vpColVector.h> 50 #include <visp3/core/vpMatrix.h> 51 #include <visp3/core/vpTime.h> 52 #include <visp3/io/vpParseArgv.h> 55 #define GETOPTARGS "cdn:i:pf:R:C:vh" 65 void usage(
const char *name,
const char *badparam)
68 Test matrix pseudo-inverse.\n\ 69 Outputs a comparison of the results obtained by supported 3rd parties.\n\ 72 %s [-n <number of matrices>] [-f <plot filename>]\n\ 73 [-R <number of rows>] [-C <number of columns>]\n\ 74 [-i <number of iterations>] [-p] [-h]\n", name);
78 -n <number of matrices> \n\ 79 Number of matrices inverted during each test loop.\n\ 81 -i <number of iterations> \n\ 82 Number of iterations of the test.\n\ 84 -f <plot filename> \n\ 85 Set output path for plot output.\n\ 86 The plot logs the times of \n\ 87 the different inversion methods: \n\ 88 QR,LU,Cholesky and Pseudo-inverse.\n\ 90 -R <number of rows>\n\ 91 Number of rows of the automatically generated matrices \n\ 94 -C <number of columns>\n\ 95 Number of colums of the automatically generated matrices \n\ 99 Plot into filename in the gnuplot format. \n\ 100 If this option is used, tests results will be logged \n\ 101 into a filename specified with -f.\n\ 104 Print the help.\n\n");
107 fprintf(stderr,
"ERROR: \n");
108 fprintf(stderr,
"\nBad parameter [%s]\n", badparam);
119 bool getOptions(
int argc,
const char **argv,
unsigned int &nb_matrices,
unsigned int &nb_iterations,
120 bool &use_plot_file, std::string &plotfile,
unsigned int &nbrows,
unsigned int &nbcols,
bool &verbose)
128 usage(argv[0], NULL);
132 nb_matrices = (
unsigned int)atoi(optarg_);
135 nb_iterations = (
unsigned int)atoi(optarg_);
139 use_plot_file =
true;
142 use_plot_file =
true;
145 nbrows = (
unsigned int)atoi(optarg_);
148 nbcols = (
unsigned int)atoi(optarg_);
159 usage(argv[0], optarg_);
165 if ((c == 1) || (c == -1)) {
167 usage(argv[0], NULL);
168 std::cerr <<
"ERROR: " << std::endl;
169 std::cerr <<
" Bad argument " << optarg_ << std::endl << std::endl;
176 vpMatrix make_random_matrix(
unsigned int nbrows,
unsigned int nbcols)
181 for (
unsigned int i = 0; i < A.
getRows(); i++) {
182 for (
unsigned int j = 0; j < A.
getCols(); j++) {
183 A[i][j] = (double)rand() / (double)RAND_MAX;
190 void create_bench_random_matrix(
unsigned int nb_matrices,
unsigned int nb_rows,
unsigned int nb_cols,
bool verbose,
191 std::vector<vpMatrix> &bench)
194 std::cout <<
"Create a bench of " << nb_matrices <<
" " << nb_rows <<
" by " << nb_cols <<
" matrices" << std::endl;
196 for (
unsigned int i = 0; i < nb_matrices; i++) {
197 vpMatrix M = make_random_matrix(nb_rows, nb_cols);
202 int test_pseudo_inverse(
const std::vector<vpMatrix> &A,
const std::vector<vpMatrix> &Api)
204 double allowed_error = 1e-3;
206 for (
unsigned int i = 0; i < A.size(); i++) {
207 double error = (A[i] * Api[i] * A[i] - A[i]).euclideanNorm();
208 if (error > allowed_error) {
209 std::cout <<
"Bad pseudo-inverse [" << i <<
"]: euclidean norm: " << error << std::endl;
216 int test_pseudo_inverse(
const std::vector<vpMatrix> &A,
const std::vector<vpMatrix> &Api,
217 const std::vector<vpColVector> &sv,
const std::vector<vpMatrix> &imA,
218 const std::vector<vpMatrix> &imAt,
const std::vector<vpMatrix> &kerAt)
220 double allowed_error = 1e-3;
222 for (
unsigned int i = 0; i < A.size(); i++) {
223 double error = (A[i] * Api[i] * A[i] - A[i]).euclideanNorm();
224 if (error > allowed_error) {
225 std::cout <<
"Bad pseudo-inverse [" << i <<
"]: euclidean norm: " << error << std::endl;
231 for (
unsigned int i = 0; i < kerAt.size(); i++) {
232 if (kerAt[i].size()) {
233 vpMatrix nullspace = A[i] * kerAt[i].
t();
236 if (error > allowed_error) {
237 std::cout <<
"Bad kernel [" << i <<
"]: euclidean norm: " << error << std::endl;
244 for (
unsigned int i = 0; i < kerAt.size(); i++) {
245 unsigned int rank = imA[i].getCols();
246 vpMatrix U, S(rank, A[i].getCols()), Vt(A[i].getCols(), A[i].getCols());
249 for (
unsigned int j = 0; j < rank; j++)
252 Vt.
insert(imAt[i].t(), 0, 0);
253 Vt.insert(kerAt[i], imAt[i].getCols(), 0);
255 double error = (U * S * Vt - A[i]).euclideanNorm();
257 if (error > allowed_error) {
258 std::cout <<
"Bad imA, imAt, sv, kerAt [" << i <<
"]: euclidean norm: " << error << std::endl;
266 int test_pseudo_inverse_default(
bool verbose,
const std::vector<vpMatrix> &bench, std::vector<double> &time)
269 std::cout <<
"Test pseudo-inverse using default 3rd party" << std::endl;
271 std::cout <<
" Pseudo-inverse on a " << bench[0].getRows() <<
"x" << bench[0].getCols() <<
" matrix" << std::endl;
273 size_t size = bench.size();
274 std::vector<vpMatrix> PI(size), imA(size), imAt(size), kerAt(size);
275 std::vector<vpColVector> sv(size);
276 int ret = EXIT_SUCCESS;
279 unsigned int test = 0;
281 for (
unsigned int i = 0; i < bench.size(); i++) {
282 PI[i] = bench[i].pseudoInverse();
285 for (
unsigned int i = 0; i < time.size(); i++) {
286 ret += test_pseudo_inverse(bench, PI);
292 for (
unsigned int i = 0; i < bench.size(); i++) {
293 bench[i].pseudoInverse(PI[i]);
296 for (
unsigned int i = 0; i < time.size(); i++) {
297 ret += test_pseudo_inverse(bench, PI);
303 for (
unsigned int i = 0; i < bench.size(); i++) {
304 bench[i].pseudoInverse(PI[i], sv[i]);
307 for (
unsigned int i = 0; i < time.size(); i++) {
308 ret += test_pseudo_inverse(bench, PI);
314 for (
unsigned int i = 0; i < bench.size(); i++) {
315 bench[i].pseudoInverse(PI[i], sv[i], 1e-6, imA[i], imAt[i], kerAt[i]);
319 for (
unsigned int i = 0; i < time.size(); i++) {
320 ret += test_pseudo_inverse(bench, PI, sv, imA, imAt, kerAt);
326 #if defined(VISP_HAVE_EIGEN3) 327 int test_pseudo_inverse_eigen3(
bool verbose,
const std::vector<vpMatrix> &bench, std::vector<double> &time)
330 std::cout <<
"Test pseudo-inverse using Eigen3 3rd party" << std::endl;
332 std::cout <<
" Pseudo-inverse on a " << bench[0].getRows() <<
"x" << bench[0].getCols() <<
" matrix" << std::endl;
334 size_t size = bench.size();
335 std::vector<vpMatrix> PI(size), imA(size), imAt(size), kerAt(size);
336 std::vector<vpColVector> sv(size);
337 int ret = EXIT_SUCCESS;
340 unsigned int test = 0;
342 for (
unsigned int i = 0; i < bench.size(); i++) {
343 PI[i] = bench[i].pseudoInverseEigen3();
346 for (
unsigned int i = 0; i < time.size(); i++) {
347 ret += test_pseudo_inverse(bench, PI);
353 for (
unsigned int i = 0; i < bench.size(); i++) {
354 bench[i].pseudoInverseEigen3(PI[i]);
357 for (
unsigned int i = 0; i < time.size(); i++) {
358 ret += test_pseudo_inverse(bench, PI);
364 for (
unsigned int i = 0; i < bench.size(); i++) {
365 bench[i].pseudoInverseEigen3(PI[i], sv[i]);
368 for (
unsigned int i = 0; i < time.size(); i++) {
369 ret += test_pseudo_inverse(bench, PI);
375 for (
unsigned int i = 0; i < bench.size(); i++) {
376 bench[i].pseudoInverseEigen3(PI[i], sv[i], 1e-6, imA[i], imAt[i], kerAt[i]);
380 for (
unsigned int i = 0; i < time.size(); i++) {
381 ret += test_pseudo_inverse(bench, PI, sv, imA, imAt, kerAt);
388 #if defined(VISP_HAVE_LAPACK) 389 int test_pseudo_inverse_lapack(
bool verbose,
const std::vector<vpMatrix> &bench, std::vector<double> &time)
392 std::cout <<
"Test pseudo-inverse using Lapack 3rd party" << std::endl;
394 std::cout <<
" Pseudo-inverse on a " << bench[0].getRows() <<
"x" << bench[0].getCols() <<
" matrix" << std::endl;
396 size_t size = bench.size();
397 std::vector<vpMatrix> PI(size), imA(size), imAt(size), kerAt(size);
398 std::vector<vpColVector> sv(size);
399 int ret = EXIT_SUCCESS;
402 unsigned int test = 0;
404 for (
unsigned int i = 0; i < bench.size(); i++) {
405 PI[i] = bench[i].pseudoInverseLapack();
408 for (
unsigned int i = 0; i < time.size(); i++) {
409 ret += test_pseudo_inverse(bench, PI);
415 for (
unsigned int i = 0; i < bench.size(); i++) {
416 bench[i].pseudoInverseLapack(PI[i]);
419 for (
unsigned int i = 0; i < time.size(); i++) {
420 ret += test_pseudo_inverse(bench, PI);
426 for (
unsigned int i = 0; i < bench.size(); i++) {
427 bench[i].pseudoInverseLapack(PI[i], sv[i]);
430 for (
unsigned int i = 0; i < time.size(); i++) {
431 ret += test_pseudo_inverse(bench, PI);
437 for (
unsigned int i = 0; i < bench.size(); i++) {
438 bench[i].pseudoInverseLapack(PI[i], sv[i], 1e-6, imA[i], imAt[i], kerAt[i]);
442 for (
unsigned int i = 0; i < time.size(); i++) {
443 ret += test_pseudo_inverse(bench, PI, sv, imA, imAt, kerAt);
450 #if defined(VISP_HAVE_GSL) 451 int test_pseudo_inverse_gsl(
bool verbose,
const std::vector<vpMatrix> &bench, std::vector<double> &time)
454 std::cout <<
"Test pseudo-inverse using Gsl 3rd party" << std::endl;
456 std::cout <<
" Pseudo-inverse on a " << bench[0].getRows() <<
"x" << bench[0].getCols() <<
" matrix" << std::endl;
458 size_t size = bench.size();
459 std::vector<vpMatrix> PI(size), imA(size), imAt(size), kerAt(size);
460 std::vector<vpColVector> sv(size);
461 int ret = EXIT_SUCCESS;
464 unsigned int test = 0;
466 for (
unsigned int i = 0; i < bench.size(); i++) {
467 PI[i] = bench[i].pseudoInverseGsl();
470 for (
unsigned int i = 0; i < time.size(); i++) {
471 ret += test_pseudo_inverse(bench, PI);
477 for (
unsigned int i = 0; i < bench.size(); i++) {
478 bench[i].pseudoInverseGsl(PI[i]);
481 for (
unsigned int i = 0; i < time.size(); i++) {
482 ret += test_pseudo_inverse(bench, PI);
488 for (
unsigned int i = 0; i < bench.size(); i++) {
489 bench[i].pseudoInverseGsl(PI[i], sv[i]);
492 for (
unsigned int i = 0; i < time.size(); i++) {
493 ret += test_pseudo_inverse(bench, PI);
499 for (
unsigned int i = 0; i < bench.size(); i++) {
500 bench[i].pseudoInverseGsl(PI[i], sv[i], 1e-6, imA[i], imAt[i], kerAt[i]);
504 for (
unsigned int i = 0; i < time.size(); i++) {
505 ret += test_pseudo_inverse(bench, PI, sv, imA, imAt, kerAt);
512 #if (VISP_HAVE_OPENCV_VERSION >= 0x020101) 513 int test_pseudo_inverse_opencv(
bool verbose,
const std::vector<vpMatrix> &bench, std::vector<double> &time)
516 std::cout <<
"Test pseudo-inverse using OpenCV 3rd party" << std::endl;
518 std::cout <<
" Pseudo-inverse on a " << bench[0].getRows() <<
"x" << bench[0].getCols() <<
" matrix" << std::endl;
520 size_t size = bench.size();
521 std::vector<vpMatrix> PI(size), imA(size), imAt(size), kerAt(size);
522 std::vector<vpColVector> sv(size);
523 int ret = EXIT_SUCCESS;
526 unsigned int test = 0;
528 for (
unsigned int i = 0; i < bench.size(); i++) {
529 PI[i] = bench[i].pseudoInverseOpenCV();
532 for (
unsigned int i = 0; i < time.size(); i++) {
533 ret += test_pseudo_inverse(bench, PI);
539 for (
unsigned int i = 0; i < bench.size(); i++) {
540 bench[i].pseudoInverseOpenCV(PI[i]);
543 for (
unsigned int i = 0; i < time.size(); i++) {
544 ret += test_pseudo_inverse(bench, PI);
550 for (
unsigned int i = 0; i < bench.size(); i++) {
551 bench[i].pseudoInverseOpenCV(PI[i], sv[i]);
554 for (
unsigned int i = 0; i < time.size(); i++) {
555 ret += test_pseudo_inverse(bench, PI);
561 for (
unsigned int i = 0; i < bench.size(); i++) {
562 bench[i].pseudoInverseOpenCV(PI[i], sv[i], 1e-6, imA[i], imAt[i], kerAt[i]);
566 for (
unsigned int i = 0; i < time.size(); i++) {
567 ret += test_pseudo_inverse(bench, PI, sv, imA, imAt, kerAt);
574 void save_time(
const std::string &method,
unsigned int nrows,
unsigned int ncols,
bool verbose,
bool use_plot_file,
575 std::ofstream &of,
const std::vector<double> &time)
577 for (
size_t i = 0; i < time.size(); i++) {
579 of << time[i] <<
"\t";
581 std::cout <<
" " << method <<
" svd(" << nrows <<
"x" << ncols <<
")" 582 <<
" test " << i <<
": " << time[i] << std::endl;
587 int main(
int argc,
const char *argv[])
590 #if defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_LAPACK) || (VISP_HAVE_OPENCV_VERSION >= 0x020101) || \ 591 defined(VISP_HAVE_GSL) 592 unsigned int nb_matrices = 10;
593 unsigned int nb_iterations = 10;
594 unsigned int nb_rows = 12;
595 unsigned int nb_cols = 6;
596 bool verbose =
false;
597 std::string plotfile(
"plot-pseudo-inv.csv");
598 bool use_plot_file =
false;
601 unsigned int nb_svd_functions = 4;
602 unsigned int nb_test_matrix_size = 3;
603 std::vector<double> time(nb_svd_functions);
604 std::vector<unsigned int> nrows(nb_test_matrix_size), ncols(nb_test_matrix_size);
607 if (getOptions(argc, argv, nb_matrices, nb_iterations, use_plot_file, plotfile, nb_rows, nb_cols, verbose) ==
612 for (
unsigned int s = 0; s < nb_test_matrix_size; s++) {
627 of.open(plotfile.c_str());
631 for (
unsigned int s = 0; s < nb_test_matrix_size; s++) {
632 for (
unsigned int i = 0; i < nb_svd_functions; i++)
633 of <<
"\"default " << nrows[s] <<
"x" << ncols[s] <<
" test " << i <<
"\"" 636 #if defined(VISP_HAVE_LAPACK) 637 for (
unsigned int i = 0; i < nb_svd_functions; i++)
638 of <<
"\"Lapack " << nrows[s] <<
"x" << ncols[s] <<
" test " << i <<
"\"" 641 #if defined(VISP_HAVE_EIGEN3) 642 for (
unsigned int i = 0; i < nb_svd_functions; i++)
643 of <<
"\"Eigen3 " << nrows[s] <<
"x" << ncols[s] <<
" test " << i <<
"\"" 646 #if (VISP_HAVE_OPENCV_VERSION >= 0x020101) 647 for (
unsigned int i = 0; i < nb_svd_functions; i++)
648 of <<
"\"OpenCV " << nrows[s] <<
"x" << ncols[s] <<
" test " << i <<
"\"" 651 #if defined(VISP_HAVE_GSL) 652 for (
unsigned int i = 0; i < nb_svd_functions; i++)
653 of <<
"\"GSL " << nrows[s] <<
"x" << ncols[s] <<
" test " << i <<
"\"" 660 int ret = EXIT_SUCCESS;
661 for (
unsigned int iter = 0; iter < nb_iterations; iter++) {
666 for (
unsigned int s = 0; s < nb_test_matrix_size; s++) {
667 std::vector<vpMatrix> bench_random_matrices;
668 create_bench_random_matrix(nb_matrices, nrows[s], ncols[s], verbose, bench_random_matrices);
670 ret += test_pseudo_inverse_default(verbose, bench_random_matrices, time);
671 save_time(
"default -", nrows[s], ncols[s], verbose, use_plot_file, of, time);
673 #if defined(VISP_HAVE_LAPACK) 674 ret += test_pseudo_inverse_lapack(verbose, bench_random_matrices, time);
675 save_time(
"Lapack -", nrows[s], ncols[s], verbose, use_plot_file, of, time);
678 #if defined(VISP_HAVE_EIGEN3) 679 ret += test_pseudo_inverse_eigen3(verbose, bench_random_matrices, time);
680 save_time(
"Eigen3 -", nrows[s], ncols[s], verbose, use_plot_file, of, time);
683 #if (VISP_HAVE_OPENCV_VERSION >= 0x020101) 684 ret += test_pseudo_inverse_opencv(verbose, bench_random_matrices, time);
685 save_time(
"OpenCV -", nrows[s], ncols[s], verbose, use_plot_file, of, time);
688 #if defined(VISP_HAVE_GSL) 689 ret += test_pseudo_inverse_gsl(verbose, bench_random_matrices, time);
690 save_time(
"GSL -", nrows[s], ncols[s], verbose, use_plot_file, of, time);
698 std::cout <<
"Result saved in " << plotfile << std::endl;
701 if (ret == EXIT_SUCCESS) {
702 std::cout <<
"Test succeed" << std::endl;
704 std::cout <<
"Test failed" << std::endl;
711 std::cout <<
"Test does nothing since you dont't have Eigen3, Lapack, " 712 "OpenCV or GSL 3rd party" Implementation of a matrix and operations on matrices.
void resize(const unsigned int nrows, const unsigned int ncols, const bool flagNullify=true, const bool recopy_=true)
error that can be emited by ViSP classes.
unsigned int getCols() const
VISP_EXPORT double measureTimeMs()
static bool parse(int *argcPtr, const char **argv, vpArgvInfo *argTable, int flags)
void insert(const vpMatrix &A, const unsigned int r, const unsigned int c)
unsigned int getRows() const
double euclideanNorm() const
const std::string & getStringMessage(void) const
Send a reference (constant) related the error message (can be empty).