Visual Servoing Platform  version 3.4.1 under development (2021-10-17)
testMatrixDeterminant.cpp
1 /****************************************************************************
2  *
3  * ViSP, open source Visual Servoing Platform software.
4  * Copyright (C) 2005 - 2019 by Inria. All rights reserved.
5  *
6  * This software is free software; you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation; either version 2 of the License, or
9  * (at your option) any later version.
10  * See the file LICENSE.txt at the root directory of this source
11  * distribution for additional information about the GNU GPL.
12  *
13  * For using ViSP with software that can not be combined with the GNU
14  * GPL, please contact Inria about acquiring a ViSP Professional
15  * Edition License.
16  *
17  * See http://visp.inria.fr for more information.
18  *
19  * This software was developed at:
20  * Inria Rennes - Bretagne Atlantique
21  * Campus Universitaire de Beaulieu
22  * 35042 Rennes Cedex
23  * France
24  *
25  * If you have questions regarding the use of this file, please contact
26  * Inria at visp@inria.fr
27  *
28  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
29  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
30  *
31  * Description:
32  * Test various determinant computation methods.
33  *
34  * Authors:
35  * Fabien Spindler
36  *
37  *****************************************************************************/
38 
44 #include <visp3/core/vpMatrix.h>
45 #include <visp3/core/vpTime.h>
46 #include <visp3/io/vpParseArgv.h>
47 
48 // List of allowed command line options
49 #define GETOPTARGS "cdn:i:pf:R:C:vh"
50 
59 void usage(const char *name, const char *badparam)
60 {
61  fprintf(stdout, "\n\
62 Test matrix inversions\n\
63 using LU, QR and Cholesky methods as well as Pseudo-inverse.\n\
64 Outputs a comparison of these methods.\n\
65 \n\
66 SYNOPSIS\n\
67  %s [-n <number of matrices>] [-f <plot filename>]\n\
68  [-R <number of rows>] [-C <number of columns>]\n\
69  [-i <number of iterations>] [-p] [-h]\n", name);
70 
71  fprintf(stdout, "\n\
72 OPTIONS: Default\n\
73  -n <number of matrices> \n\
74  Number of matrices inverted during each test loop.\n\
75 \n\
76  -i <number of iterations> \n\
77  Number of iterations of the test.\n\
78 \n\
79  -f <plot filename> \n\
80  Set output path for plot output.\n\
81  The plot logs the times of \n\
82  the different inversion methods: \n\
83  QR,LU,Cholesky and Pseudo-inverse.\n\
84 \n\
85  -R <number of rows>\n\
86  Number of rows of the automatically generated matrices \n\
87  we test on.\n\
88 \n\
89  -C <number of columns>\n\
90  Number of colums of the automatically generated matrices \n\
91  we test on.\n\
92 \n\
93  -p \n\
94  Plot into filename in the gnuplot format. \n\
95  If this option is used, tests results will be logged \n\
96  into a filename specified with -f.\n\
97 \n\
98  -h\n\
99  Print the help.\n\n");
100 
101  if (badparam) {
102  fprintf(stderr, "ERROR: \n");
103  fprintf(stderr, "\nBad parameter [%s]\n", badparam);
104  }
105 }
106 
114 bool getOptions(int argc, const char **argv, unsigned int &nb_matrices, unsigned int &nb_iterations,
115  bool &use_plot_file, std::string &plotfile, unsigned int &nbrows, unsigned int &nbcols, bool &verbose)
116 {
117  const char *optarg_;
118  int c;
119  while ((c = vpParseArgv::parse(argc, argv, GETOPTARGS, &optarg_)) > 1) {
120 
121  switch (c) {
122  case 'h':
123  usage(argv[0], NULL);
124  return false;
125  break;
126  case 'n':
127  nb_matrices = (unsigned int)atoi(optarg_);
128  break;
129  case 'i':
130  nb_iterations = (unsigned int)atoi(optarg_);
131  break;
132  case 'f':
133  plotfile = optarg_;
134  use_plot_file = true;
135  break;
136  case 'p':
137  use_plot_file = true;
138  break;
139  case 'R':
140  nbrows = (unsigned int)atoi(optarg_);
141  break;
142  case 'C':
143  nbcols = (unsigned int)atoi(optarg_);
144  break;
145  case 'v':
146  verbose = true;
147  break;
148  // add default options -c -d
149  case 'c':
150  break;
151  case 'd':
152  break;
153  default:
154  usage(argv[0], optarg_);
155  return false;
156  break;
157  }
158  }
159 
160  if ((c == 1) || (c == -1)) {
161  // standalone param or error
162  usage(argv[0], NULL);
163  std::cerr << "ERROR: " << std::endl;
164  std::cerr << " Bad argument " << optarg_ << std::endl << std::endl;
165  return false;
166  }
167 
168  return true;
169 }
170 
171 vpMatrix make_random_matrix(unsigned int nbrows, unsigned int nbcols)
172 {
173  vpMatrix A;
174  A.resize(nbrows, nbcols);
175 
176  for (unsigned int i = 0; i < A.getRows(); i++)
177  for (unsigned int j = 0; j < A.getCols(); j++)
178  A[i][j] = (double)rand() / (double)RAND_MAX;
179  return A;
180 }
181 
182 void create_bench(unsigned int nb_matrices, unsigned int nb_rows, unsigned int nb_cols, bool verbose,
183  std::vector<vpMatrix> &bench)
184 {
185  if (verbose)
186  std::cout << "Create a bench of " << nb_matrices << " " << nb_rows << " by " << nb_cols << " matrices" << std::endl;
187  bench.clear();
188  for (unsigned int i = 0; i < nb_matrices; i++) {
189  vpMatrix M = make_random_matrix(nb_rows, nb_cols);
190  bench.push_back(M);
191  }
192 }
193 
194 void test_det_default(bool verbose, const std::vector<vpMatrix> &bench, double &time, std::vector<double> &result)
195 {
196  if (verbose)
197  std::cout << "Test determinant using default method" << std::endl;
198  // Compute inverse
199  if (verbose)
200  std::cout << " Matrix size: " << bench[0].AtA().getRows() << "x" << bench[0].AtA().getCols() << std::endl;
201 
202  result.resize(bench.size());
203  double t = vpTime::measureTimeMs();
204  for (unsigned int i = 0; i < bench.size(); i++) {
205  result[i] = bench[i].AtA().det();
206  }
207  time = vpTime::measureTimeMs() - t;
208 }
209 
210 #if defined(VISP_HAVE_EIGEN3)
211 void test_det_eigen3(bool verbose, const std::vector<vpMatrix> &bench, double &time, std::vector<double> &result)
212 {
213  if (verbose)
214  std::cout << "Test determinant using Eigen3 3rd party" << std::endl;
215  // Compute inverse
216  if (verbose)
217  std::cout << " Matrix size: " << bench[0].AtA().getRows() << "x" << bench[0].AtA().getCols() << std::endl;
218 
219  result.resize(bench.size());
220  double t = vpTime::measureTimeMs();
221  for (unsigned int i = 0; i < bench.size(); i++) {
222  result[i] = bench[i].AtA().detByLUEigen3();
223  }
224  time = vpTime::measureTimeMs() - t;
225 }
226 #endif
227 
228 #if defined(VISP_HAVE_LAPACK)
229 void test_det_lapack(bool verbose, const std::vector<vpMatrix> &bench, double &time, std::vector<double> &result)
230 {
231  if (verbose)
232  std::cout << "Test determinant using Lapack 3rd party" << std::endl;
233  // Compute inverse
234  if (verbose)
235  std::cout << " Matrix size: " << bench[0].AtA().getRows() << "x" << bench[0].AtA().getCols() << std::endl;
236 
237  result.resize(bench.size());
238  double t = vpTime::measureTimeMs();
239  for (unsigned int i = 0; i < bench.size(); i++) {
240  result[i] = bench[i].AtA().detByLULapack();
241  }
242  time = vpTime::measureTimeMs() - t;
243 }
244 #endif
245 
246 #if (VISP_HAVE_OPENCV_VERSION >= 0x020101)
247 void test_det_opencv(bool verbose, const std::vector<vpMatrix> &bench, double &time, std::vector<double> &result)
248 {
249  if (verbose)
250  std::cout << "Test determinant using OpenCV 3rd party" << std::endl;
251  // Compute inverse
252  if (verbose)
253  std::cout << " Matrix size: " << bench[0].AtA().getRows() << "x" << bench[0].AtA().getCols() << std::endl;
254 
255  result.resize(bench.size());
256  double t = vpTime::measureTimeMs();
257  for (unsigned int i = 0; i < bench.size(); i++) {
258  result[i] = bench[i].AtA().detByLUOpenCV();
259  }
260  time = vpTime::measureTimeMs() - t;
261 }
262 #endif
263 
264 void save_time(const std::string &method, bool verbose, bool use_plot_file, std::ofstream &of, double time)
265 {
266  if (use_plot_file)
267  of << time << "\t";
268  if (verbose || !use_plot_file) {
269  std::cout << method << time << std::endl;
270  }
271 }
272 
273 int main(int argc, const char *argv[])
274 {
275  try {
276 #if defined(VISP_HAVE_EIGEN3) || defined(VISP_HAVE_LAPACK) || (VISP_HAVE_OPENCV_VERSION >= 0x020101)
277  unsigned int nb_matrices = 1000;
278  unsigned int nb_iterations = 10;
279  unsigned int nb_rows = 6;
280  unsigned int nb_cols = 6;
281  bool verbose = false;
282  std::string plotfile("plot-det.csv");
283  bool use_plot_file = false;
284  std::ofstream of;
285 
286  // Read the command line options
287  if (getOptions(argc, argv, nb_matrices, nb_iterations, use_plot_file, plotfile, nb_rows, nb_cols, verbose) ==
288  false) {
289  exit(-1);
290  }
291 
292  if (use_plot_file) {
293  of.open(plotfile.c_str());
294  of << "iter"
295  << "\t";
296 
297  of << "\"Determinant default\""
298  << "\t";
299 
300 #if defined(VISP_HAVE_LAPACK)
301  of << "\"Determinant Lapack\""
302  << "\t";
303 #endif
304 #if defined(VISP_HAVE_EIGEN3)
305  of << "\"Determinant Eigen3\""
306  << "\t";
307 #endif
308 #if (VISP_HAVE_OPENCV_VERSION >= 0x020101)
309  of << "\"Determinant OpenCV\""
310  << "\t";
311 #endif
312  of << std::endl;
313  }
314 
315  int ret = EXIT_SUCCESS;
316  for (unsigned int iter = 0; iter < nb_iterations; iter++) {
317  std::vector<vpMatrix> bench;
318  create_bench(nb_matrices, nb_rows, nb_cols, verbose, bench);
319 
320  if (use_plot_file)
321  of << iter << "\t";
322 
323  double time;
324 
325  std::vector<double> result_default;
326  test_det_default(verbose, bench, time, result_default);
327  save_time("Determinant default: ", verbose, use_plot_file, of, time);
328 
329 #if defined(VISP_HAVE_LAPACK)
330  std::vector<double> result_lapack;
331  test_det_lapack(verbose, bench, time, result_lapack);
332  save_time("Determinant by Lapack: ", verbose, use_plot_file, of, time);
333 #endif
334 
335 #if defined(VISP_HAVE_EIGEN3)
336  std::vector<double> result_eigen3;
337  test_det_eigen3(verbose, bench, time, result_eigen3);
338  save_time("Determinant by Eigen3: ", verbose, use_plot_file, of, time);
339 #endif
340 
341 #if (VISP_HAVE_OPENCV_VERSION >= 0x020101)
342  std::vector<double> result_opencv;
343  test_det_opencv(verbose, bench, time, result_opencv);
344  save_time("Determinant by OpenCV: ", verbose, use_plot_file, of, time);
345 #endif
346 
347  if (use_plot_file)
348  of << std::endl;
349 
350 #if defined(VISP_HAVE_LAPACK) && (VISP_HAVE_OPENCV_VERSION >= 0x020101)
351  // Compare results
352  for (unsigned int i = 0; i < bench.size(); i++) {
353  if (std::fabs(result_lapack[i] - result_opencv[i]) > 1e-6) {
354  std::cout << "Determinant differ between Lapack and OpenCV: " << result_lapack[i] << " " << result_opencv[i]
355  << std::endl;
356  ret = EXIT_FAILURE;
357  }
358  }
359 #endif
360 #if defined(VISP_HAVE_EIGEN3) && (VISP_HAVE_OPENCV_VERSION >= 0x020101)
361  // Compare results
362  for (unsigned int i = 0; i < bench.size(); i++) {
363  if (std::fabs(result_eigen3[i] - result_opencv[i]) > 1e-6) {
364  std::cout << "Determinant differ between Eigen3 and OpenCV: " << result_eigen3[i] << " " << result_opencv[i]
365  << std::endl;
366  ret = EXIT_FAILURE;
367  }
368  }
369 #endif
370 #if defined(VISP_HAVE_EIGEN3) && defined(VISP_HAVE_LAPACK)
371  // Compare results
372  for (unsigned int i = 0; i < bench.size(); i++) {
373  if (std::fabs(result_eigen3[i] - result_lapack[i]) > 1e-6) {
374  std::cout << "Determinant differ between Eigen3 and Lapack: " << result_eigen3[i] << " " << result_lapack[i]
375  << std::endl;
376  ret = EXIT_FAILURE;
377  }
378  }
379 #endif
380  }
381  if (use_plot_file) {
382  of.close();
383  std::cout << "Result saved in " << plotfile << std::endl;
384  }
385 
386  if (ret == EXIT_SUCCESS) {
387  std::cout << "Test succeed" << std::endl;
388  } else {
389  std::cout << "Test failed" << std::endl;
390  }
391 
392  return ret;
393 #else
394  (void)argc;
395  (void)argv;
396  std::cout << "Test does nothing since you dont't have Lapack, Eigen3 or OpenCV 3rd party"
397  << std::endl;
398  return EXIT_SUCCESS;
399 #endif
400  } catch (const vpException &e) {
401  std::cout << "Catch an exception: " << e.getStringMessage() << std::endl;
402  return EXIT_FAILURE;
403  }
404 }
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:153
void resize(unsigned int nrows, unsigned int ncols, bool flagNullify=true, bool recopy_=true)
Definition: vpArray2D.h:304
error that can be emited by ViSP classes.
Definition: vpException.h:71
unsigned int getRows() const
Definition: vpArray2D.h:289
VISP_EXPORT double measureTimeMs()
Definition: vpTime.cpp:126
static bool parse(int *argcPtr, const char **argv, vpArgvInfo *argTable, int flags)
Definition: vpParseArgv.cpp:69
unsigned int getCols() const
Definition: vpArray2D.h:279
const std::string & getStringMessage() const
Send a reference (constant) related the error message (can be empty).
Definition: vpException.cpp:92