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