Helios++
Helios software for LiDAR simulations
FluxionumTest.h
1 #pragma once
2 
3 #include <fluxionum/UnivariateNewtonRaphsonMinimizer.h>
4 #include <fluxionum/DesignMatrix.h>
5 #include <fluxionum/TemporalDesignMatrix.h>
6 #include <fluxionum/IndexedDesignMatrix.h>
7 #include <fluxionum/DiffDesignMatrix.h>
8 #include <fluxionum/DiffDesignMatrixInterpolator.h>
9 #include <fluxionum/LinearPiecesFunction.h>
10 #include <fluxionum/ParametricLinearPiecesFunction.h>
11 #include <fluxionum/FixedIterativeEulerMethod.h>
12 #include <fluxionum/FixedParametricIterativeEulerMethod.h>
13 #include <fluxionum/ClosestLesserSampleFunction.h>
14 #include <fluxionum/ParametricClosestLesserSampleFunction.h>
15 
16 #include <armadillo>
17 
18 #include <functional>
19 #include <cmath>
20 
21 namespace HeliosTests{
22 
23 using namespace fluxionum;
24 using std::function;
25 
31 class FluxionumTest : public BaseTest{
32 public:
33  // *** ATTRIBUTES *** //
34  // ******************** //
38  double const eps = 0.00001;
42  std::string testDir;
43 
44  // *** CONSTRUCTOR *** //
45  // ********************* //
49  FluxionumTest(std::string testDir) :
50  BaseTest("Fluxionum test"),
51  testDir(testDir)
52  {}
53 
54  // *** R U N *** //
55  // *************** //
59  bool run() override;
60 
61  // *** SUB-TESTS *** //
62  // ******************* //
67  bool testUnivariateNewtonRaphsonMinimization();
72  bool testDesignMatrixBuilding();
77  bool testDesignMatrixMethods();
82  bool testDiffDesignMatrix();
87  bool testDesignFunctions();
88 };
89 
90 // *** R U N *** //
91 // *************** //
93  // Run tests
94  if(!testUnivariateNewtonRaphsonMinimization()) return false;
95  if(!testDesignMatrixBuilding()) return false;
96  if(!testDesignMatrixMethods()) return false;
97  if(!testDiffDesignMatrix()) return false;
98  if(!testDesignFunctions()) return false;
99  return true;
100 }
101 
102 // *** SUB-TESTS *** //
103 // ******************* //
105  double expected = -3.0;
107  [] (double x) -> double {return std::pow(x, 2)/2.0 + 3.0*x;},
108  [] (double x) -> double {return x+3.0;},
109  [] (double x) -> double {return 1;}
110  );
111  double x = unrm.argmin(0.0);
112  return std::fabs(x-expected) <= eps;
113 }
114 
116  // Validation functions
117  std::function<bool( // Validate DesignMatrix
118  DesignMatrix<double> &, vector<string> const &,
119  arma::Mat<double> const &
120  )> validateDesignMatrix = [&] (
122  vector<string> const &colNames,
123  arma::Mat<double> const &X
124  ) -> bool
125  {
126  size_t const nRows = dm.getNumRows();
127  size_t const nCols = dm.getNumColumns();
128  if(X.n_rows != nRows || X.n_cols != nCols) return false;
129  if(dm.hasColumnNames()){
130  if(dm.getColumnNames().size() != colNames.size()) return false;
131  for(size_t j = 0 ; j < nCols ; ++j){ // Check names
132  if(dm.getColumnName(j) != colNames[j]) return false;
133  }
134  }
135  else if(!colNames.empty()) return false;
136  for(size_t i = 0 ; i < nRows ; ++i){ // Check values
137  for(size_t j = 0 ; j < nCols ; ++j){
138  if(std::fabs(dm(i, j)-X(i, j)) > eps) return false;
139  }
140  }
141  return true;
142  };
143  std::function<bool( // Validate TemporalDesignMatrix
144  TemporalDesignMatrix<double, double> &, string const &,
145  vector<string> const &,
146  arma::Col<double> const &, arma::Mat<double> const &
147  )> validateTemporalDesignMatrix = [&](
149  string const &timeName,
150  vector<string> const &colNames,
151  arma::Col<double> const &t,
152  arma::Mat<double> const &X
153  ) -> bool
154  {
155  if(!validateDesignMatrix(tdm, colNames, X)) return false;
156  size_t const nRows = tdm.getNumRows();
157  if(tdm.getTimeName() != timeName) return false;
158  for(size_t i = 0 ; i < nRows ; ++i){ // Check time
159  if(std::fabs(tdm[i]-t[i]) > eps) return false;
160  }
161  return true;
162  };
163  std::function<bool( // Validate IndexedDesignMatrix
164  IndexedDesignMatrix<int, double> &, string const&,
165  vector<string> const&,
166  vector<int> const&, arma::Mat<double> const&
167  )> validateIndexedDesignMatrix = [&](
169  string const &indexName,
170  vector<string> const &colNames,
171  vector<int> const &ids,
172  arma::Mat<double> const &X
173  ) -> bool
174  {
175  if(!validateDesignMatrix(idm, colNames, X)) return false;
176  size_t const nRows = idm.getNumRows();
177  if(ids.size() != nRows) return false;
178  if(idm.getIndexName() != indexName) return false;
179  for(size_t i = 0 ; i < nRows ; ++i){ // Check time
180  if(idm[i] != ids[i]) return false;
181  }
182  return true;
183  };
184 
185  // Basic design matrices
186  vector<string> colNamesR2({"x", "y"});
187  arma::Mat<double> X1 = arma::randn(5, 2);
188  DesignMatrix<double> dm1(X1, colNamesR2);
189  if(!validateDesignMatrix(dm1, colNamesR2, X1)) return false;
190  DesignMatrix<double> dm2(colNamesR2);
191  if(!validateDesignMatrix(dm2, colNamesR2, arma::Mat<double>()))
192  return false;
194  if(!validateDesignMatrix(dm3, vector<string>(0), arma::Mat<double>()))
195  return false;
196  DesignMatrix<double> dm4(dm1);
197  if(!validateDesignMatrix(dm4, colNamesR2, X1)) return false;
198  dm4 = dm1;
199  if(!validateDesignMatrix(dm4, colNamesR2, X1)) return false;
200 
201  // Temporal design matrices
202  vector<string> colNamesR2t({"x", "y", "t"});
203  arma::Col<double> t1 = arma::randn(5, 1);
204  arma::Mat<double> X2 = arma::randn(5, 3);
205  TemporalDesignMatrix<double, double> tdm1(X1, t1, "time", colNamesR2);
206  if(!validateTemporalDesignMatrix(tdm1, "time", colNamesR2, t1, X1))
207  return false;
209  if(!validateTemporalDesignMatrix(tdm2, "time", colNamesR2, t1, X1))
210  return false;
211  tdm2 = tdm1;
212  if(!validateTemporalDesignMatrix(tdm2, "time", colNamesR2, t1, X1))
213  return false;
214  TemporalDesignMatrix<double, double> tdm3(tdm1, t1, "TIME");
215  if(!validateTemporalDesignMatrix(tdm3, "TIME", colNamesR2, t1, X1))
216  return false;
218  if(!validateTemporalDesignMatrix(tdm4, "time", colNamesR2, t1, X1))
219  return false;
220  TemporalDesignMatrix<double, double> tdm5(X2, 2, "t", colNamesR2);
221  if(!validateTemporalDesignMatrix(
222  tdm5, "t", colNamesR2, X2.col(2), X2.cols(0, 1)
223  )) return false;
224 
225  // Indexed design matrices
226  vector<int> ids1({0, 1, 2, 3, 4});
227  IndexedDesignMatrix<int, double> idm1(X1, ids1, "index", colNamesR2);
228  if(!validateIndexedDesignMatrix(idm1, "index", colNamesR2, ids1, X1))
229  return false;
231  if(!validateIndexedDesignMatrix(idm2, "index", colNamesR2, ids1, X1))
232  return false;
233  idm2 = idm1;
234  if(!validateIndexedDesignMatrix(idm2, "index", colNamesR2, ids1, X1))
235  return false;
236  IndexedDesignMatrix<int, double> idm3(idm1, ids1, "INDEX");
237  if(!validateIndexedDesignMatrix(idm3, "INDEX", colNamesR2, ids1, X1))
238  return false;
239  IndexedDesignMatrix<int, double> idm4(idm1, ids1);
240  if(!validateIndexedDesignMatrix(idm4, "index", colNamesR2, ids1, X1))
241  return false;
242  IndexedDesignMatrix<int, double> idm5(X2, 2, "idx", colNamesR2);
243  if(!validateIndexedDesignMatrix(
244  idm5, "idx", colNamesR2,
245  vector<int>({
246  (int)X2(0, 2),
247  (int)X2(1, 2),
248  (int)X2(2, 2),
249  (int)X2(3, 2),
250  (int)X2(4, 2),
251  }),
252  X2.cols(0, 1)
253  )) return false;
254 
255  // Design matrix from file
256  vector<string> hf1({"x", "y", "z"}); // Header
257  arma::Mat<double> Xf1( // Matrix
258  "0 0 0;"
259  "0 0.1 0;"
260  "0.1 0.1 0;"
261  "0.1 0.2 0;"
262  "0.2 0.2 0;"
263  "0.2 0.2 0.1;"
264  "0.2 0.3 0.1;"
265  "0.3 0.3 0.1;"
266  "0.3 0.3 0.2"
267  );
268  string const dmf1Path = testDir + "design_matrix_header.txt";
269  DesignMatrix<double> dmf1(dmf1Path);
270  if(!validateDesignMatrix(dmf1, hf1, Xf1)) return false;
271  string const dmf2Path = testDir + "design_matrix_nonheader.txt";
272  DesignMatrix<double> dmf2(dmf2Path);
273  if(!validateDesignMatrix(dmf2, vector<string>(0), Xf1)) return false;
274 
275  // Temporal design matrix from file
276  vector<string> htf1({"x", "y"}); // Header
277  arma::Mat<double> Xtf1( // Matrix
278  "0 0;"
279  "0 0.1;"
280  "0.1 0.1;"
281  "0.1 0.2;"
282  "0.2 0.2;"
283  "0.2 0.2;"
284  "0.2 0.3;"
285  "0.3 0.3;"
286  "0.3 0.3"
287  );
288  arma::Col<double> ttf1("0 0.1 0.2 0.3 0.4 0.5 0.6 0.8 1.0");
289  string const tdmf1Path = testDir + "temporal_design_matrix_header.txt";
290  TemporalDesignMatrix<double, double> tdmf1(tdmf1Path);
291  if(!validateTemporalDesignMatrix(tdmf1, "t", htf1, ttf1, Xtf1))
292  return false;
293  string const tdmf2Path = testDir + "temporal_design_matrix_nonheader.txt";
294  TemporalDesignMatrix<double, double> tdmf2(tdmf2Path);
295  if(!validateTemporalDesignMatrix(
296  tdmf2, "time", vector<string>(0), ttf1, Xtf1
297  )) return false;
298 
299  // Indexed design matrix from file
300  vector<string> hif1({"x", "y"}); // Header
301  arma::Mat<double> Xif1( // Matrix
302  "0 0;"
303  "0 0.1;"
304  "0.1 0.1;"
305  "0.1 0.2;"
306  "0.2 0.2;"
307  "0.2 0.2;"
308  "0.2 0.3;"
309  "0.3 0.3;"
310  "0.3 0.3"
311  );
312  vector<int> iif1({0, 1, 2, 3, 4, 6, 8, 7, 5});
313  string const idmf1Path = testDir + "indexed_design_matrix_header.txt";
314  IndexedDesignMatrix<int, double> idmf1(idmf1Path);
315  if(!validateIndexedDesignMatrix(idmf1, "idx", hif1, iif1, Xif1))
316  return false;
317  string const idmf2Path = testDir + "indexed_design_matrix_nonheader.txt";
318  IndexedDesignMatrix<int, double> idmf2(idmf2Path);
319  if(!validateIndexedDesignMatrix(
320  idmf2, "index", vector<string>(0), iif1, Xif1
321  )) return false;
322 
323  // On passed return true
324  return true;
325 }
326 
328  // Prepare data for tests
330  arma::Mat<double>(
331  "0 1 2 3 4;"
332  "0 1 2 3 4;"
333  "0 1 2 3 4;"
334  )
335  );
337  arma::Mat<double>(
338  "5 6 7 8 9;"
339  "5 6 7 8 9;"
340  "5 6 7 8 9;"
341  )
342  );
344  arma::Mat<double>(
345  "1 2 3;"
346  "1 2 3;"
347  "1 2 3;"
348  ),
349  arma::Col<double>("0.0 0.1 0.2")
350  );
352  arma::Mat<double>(
353  "4 5 6;"
354  "4 5 6;"
355  "4 5 6;"
356  ),
357  arma::Col<double>("0.3 0.4 0.5")
358  );
360  arma::Mat<double>(
361  "4 5 6;"
362  "4 5 6;"
363  "4 5 6;"
364  ),
365  arma::Col<double>("0.4 0.5 0.3")
366  );
368  arma::Mat<double>(
369  "-2 -2 1;"
370  "0 0 1;"
371  "2 2 1;"
372  "4 4 1;"
373  "4 4 1;"
374  "4 4 1;"
375  "4 4 1;"
376  "4 4 1;"
377  "3 3 1;"
378  "2 2 1;"
379  "1 1 1;"
380  ),
381  arma::Col<double>("-5 -4 -3 -2 -1 0 1 2 3 4 5")
382  );
384  arma::Mat<double>(
385  "0.1 0.2 0.3;"
386  "0.2 0.3 0.4;"
387  "0.3 0.4 0.5;"
388  ),
389  vector<int>({1, 3, 5})
390  );
392  arma::Mat<double>(
393  "0.4 0.5 0.6;"
394  "0.5 0.6 0.7;"
395  "0.7 0.8 0.9;"
396  ),
397  vector<int>({6, 8, 10})
398  );
399 
400  // Validate mergeInPlace
401  DesignMatrix<double> dm(dm1);
402  dm.mergeInPlace(dm2);
403  arma::Mat<double> em( // Expected matrix
404  "0 1 2 3 4;"
405  "0 1 2 3 4;"
406  "0 1 2 3 4;"
407  "5 6 7 8 9;"
408  "5 6 7 8 9;"
409  "5 6 7 8 9;"
410  );
411  if(arma::any(arma::vectorise(arma::abs(dm.getX()-em)) > eps)) return false;
413  tdm.mergeInPlace(tdm2);
414  em = arma::Mat<double>(
415  "1 2 3;"
416  "1 2 3;"
417  "1 2 3;"
418  "4 5 6;"
419  "4 5 6;"
420  "4 5 6;"
421  );
422  arma::Col<double> et({ // Expected time
423  0.0, 0.1, 0.2, 0.3, 0.4, 0.5
424  });
425  if(arma::any(arma::vectorise(arma::abs(tdm.getX()-em)) > eps))
426  return false;
427  if(arma::any((tdm.getTimeVector()-et) >eps)) return false;
429  idm.mergeInPlace(idm2);
430  em = arma::Mat<double>(
431  "0.1 0.2 0.3;"
432  "0.2 0.3 0.4;"
433  "0.3 0.4 0.5;"
434  "0.4 0.5 0.6;"
435  "0.5 0.6 0.7;"
436  "0.7 0.8 0.9;"
437  );
438  vector<int> ei({ // Expected indices
439  1, 3, 5, 6, 8, 10
440  });
441  if(arma::any(arma::vectorise(arma::abs(idm.getX()-em)) > eps))
442  return false;
443  for(size_t i = 0 ; i < ei.size() ; ++i){
444  if(std::fabs(ei[i]-idm.getIndices()[i]) > eps) return false;
445  }
446 
447  // Validate swapColumns
448  dm = dm1;
449  dm.swapColumns(arma::uvec({0, 2, 4, 1, 3}));
450  em = arma::Mat<double>(
451  "0 2 4 1 3;"
452  "0 2 4 1 3;"
453  "0 2 4 1 3;"
454  );
455  if(arma::any(arma::vectorise(arma::abs(dm.getX()-em)) > eps)) return false;
456  tdm = tdm1;
457  tdm.swapColumns(arma::uvec({1, 2, 0}));
458  em = arma::Mat<double>(
459  "2 3 1;"
460  "2 3 1;"
461  "2 3 1;"
462  );
463  if(arma::any(arma::vectorise(arma::abs(tdm.getX()-em)) > eps))
464  return false;
465  idm = idm1;
466  idm.swapColumns(arma::uvec({2, 1, 0}));
467  em = arma::Mat<double>(
468  "0.3 0.2 0.1;"
469  "0.4 0.3 0.2;"
470  "0.5 0.4 0.3;"
471  );
472  if(arma::any(arma::vectorise(arma::abs(idm.getX()-em)) > eps))
473  return false;
474 
475  // Validate sortByTime
476  tdm3.sortByTime();
477  if(arma::any((tdm3.getTimeVector()-tdm2.getTimeVector()) > eps))
478  return false;
479 
480  // Validate shiftTime
481  et = tdm3.getTimeVector() + 13.37;
482  tdm3.shiftTime(13.37);
483  if(arma::any((tdm3.getTimeVector()-et) > eps)) return false;
484 
485  // Validate slopeFilter
486  et = arma::Col<double>("-5 -2 2 5");
487  em = arma::Mat<double>(
488  "-2 -2 1;"
489  "4 4 1;"
490  "4 4 1;"
491  "1 1 1;"
492  );
493  if(tdm4.slopeFilter(0.01) != 7) return false;
494  if(arma::any(arma::abs(tdm4.getTimeVector()-et) > eps)) return false;
495  if(arma::any(arma::abs(arma::vectorise(tdm4.getX()-em)) > eps))
496  return false;
497 
498  // On passed return true
499  return true;
500 }
501 
503  // Validation function
504  std::function<bool( // Validate DiffDesignMatrix
505  DiffDesignMatrix<double, double> &, string const&,
506  vector<string> const&,
507  arma::Col<double> const&, arma::Mat<double> const&,
508  DiffDesignMatrixType
509  )> validateDiffDesignMatrix = [&] (
511  string const &timeName,
512  vector<string> const &colNames,
513  arma::Col<double> const &t,
514  arma::Mat<double> const &A,
515  DiffDesignMatrixType diffType
516  ) -> bool
517  {
518  if(ddm.getA().n_rows != ddm.getTimeVector().n_rows) return false;
519  size_t const nRows = ddm.getNumRows();
520  size_t const nCols = ddm.getNumColumns();
521  if(A.n_rows != nRows || A.n_cols != nCols) return false;
522  if(ddm.hasColumnNames()){
523  if(ddm.getColumnNames().size() != colNames.size()) return false;
524  for(size_t j = 0 ; j < nCols ; ++j){ // Check names
525  if(ddm.getColumnName(j) != colNames[j]) return false;
526  }
527  }
528  else if(!colNames.empty()) return false;
529  for(size_t i = 0 ; i < nRows ; ++i){
530  for(size_t j = 0 ; j < nCols ; ++j){
531  if(std::fabs(ddm(i, j)-A(i, j)) > eps) return false;
532  }
533  }
534  if(diffType != ddm.getDiffType()) return false;
535  if(ddm.getTimeName() != timeName) return false;
536  for(size_t i = 0 ; i < nRows ; ++i){ // Check time
537  if(std::fabs(ddm[i]-t[i]) > eps) return false;
538  }
539  return true;
540  };
541 
542 
543  // Expected diff design matrices
544  arma::Col<double> Et1("-3 -2 -1 0 1 2");
545  arma::Col<double> Et3("-2.5 -1.5 -0.5 0.5 1.5");
546  arma::Mat<double> EA1(
547  "0.1 -0.5;"
548  "0.1 -0.3;"
549  "0.1 -0.1;"
550  "0.1 0.1;"
551  "0.1 0.3;"
552  "0.1 0.5"
553  );
554  arma::Mat<double> EA3(
555  "0.1 -0.4;"
556  "0.1 -0.2;"
557  "0.1 0.0;"
558  "0.1 0.2;"
559  "0.1 0.4;"
560  );
561  vector<string> EcolNames1({"f1", "f2"});
562  string EtimeName1 = "t";
563  string EtimeName2 = "time";
564 
565  // Build design matrix
566  vector<string> colNames1({"t", "f1", "f2"});
567  arma::Mat<double> X1 = arma::Mat<double>(
568  "-3 -0.3 0.9;"
569  "-2 -0.2 0.4;"
570  "-1 -0.1 0.1;"
571  "0 0 0;"
572  "1 0.1 0.1;"
573  "2 0.2 0.4;"
574  "3 0.3 0.9"
575  );
576  arma::Mat<double> X2 = arma::Mat<double>(
577  "3 0.3 0.9;"
578  "-1 -0.1 0.1;"
579  "0 0 0;"
580  "-3 -0.3 0.9;"
581  "1 0.1 0.1;"
582  "-2 -0.2 0.4;"
583  "2 0.2 0.4"
584  );
585  DesignMatrix<double> dm1(X1, colNames1); // time sorted
586  DesignMatrix<double> dm2(X2); // non time sorted
587 
588  // Build temporal design matrix
589  TemporalDesignMatrix<double, double> tdm1(dm1, 0); // time sorted
590  TemporalDesignMatrix<double, double> tdm2(dm2, 0); // non time sorted
591 
592  // Build diff design matrix from previous temporal design matrix
593  DiffDesignMatrix<double, double> ddm1( // built from time sorted
594  tdm1, DiffDesignMatrixType::FORWARD_FINITE_DIFFERENCES
595  );
596  DiffDesignMatrix<double, double> ddm2( // built from non time sorted
597  tdm2, DiffDesignMatrixType::FORWARD_FINITE_DIFFERENCES
598  );
599  DiffDesignMatrix<double, double> ddm3( // build from non time sorted
600  tdm2, DiffDesignMatrixType::CENTRAL_FINITE_DIFFERENCES
601  );
602 
603  // Validate built diff design matrix
604  if(!validateDiffDesignMatrix(
605  ddm1, EtimeName1, EcolNames1, Et1, EA1,
606  DiffDesignMatrixType::FORWARD_FINITE_DIFFERENCES
607  ))
608  return false;
609  if(!validateDiffDesignMatrix(
610  ddm2, EtimeName2, vector<string>(0), Et1, EA1,
611  DiffDesignMatrixType::FORWARD_FINITE_DIFFERENCES
612  )) return false;
613  if(!validateDiffDesignMatrix(
614  ddm3, EtimeName2, vector<string>(0), Et3, EA3,
615  DiffDesignMatrixType::CENTRAL_FINITE_DIFFERENCES
616  )) return false;
617 
618  // Load diff design matrix from file
619  string const ddmf1Path = testDir + "diff_design_matrix_header.txt";
620  DiffDesignMatrix<double, double> ddmf1(ddmf1Path);
621  if(!validateDiffDesignMatrix(
622  ddmf1, EtimeName1, EcolNames1, Et1, EA1,
623  DiffDesignMatrixType::FORWARD_FINITE_DIFFERENCES
624  )) return false;
625  string const ddmf2Path = testDir + "diff_design_matrix_nonheader.txt";
626  DiffDesignMatrix<double, double> ddmf2(ddmf2Path);
627  if(!validateDiffDesignMatrix(
628  ddmf2, EtimeName2, vector<string>(0), Et3, EA3,
629  DiffDesignMatrixType::CENTRAL_FINITE_DIFFERENCES
630  )) return false;
631 
632  // On passed return true
633  return true;
634 }
635 
637  // Prepare data
639  arma::Mat<double>(
640  "0 0;"
641  "2 2;"
642  "4 3;"
643  "6 3;"
644  "9 4;"
645  ),
646  0
647  );
649  arma::Mat<double>(
650  "0 0 3;"
651  "2 2 4;"
652  "4 3 3;"
653  "6 3 6;"
654  "9 4 5;"
655  ),
656  0
657  );
659  DiffDesignMatrixType::FORWARD_FINITE_DIFFERENCES
660  );
662  DiffDesignMatrixType::FORWARD_FINITE_DIFFERENCES
663  );
664 
665  // Validate LinearPiecesFunction
666  arma::Col<double> intercept(
667  tdm1.getColumnCopy(0).subvec(0,tdm1.getNumRows()-2)
668  );
669  arma::Col<double> slope = ffd1.getA().col(0);
671  DiffDesignMatrixInterpolator::makeLinearPiecesFunction(
672  ffd1,
673  slope,
674  intercept
675  );
677  DiffDesignMatrixInterpolator::makeLinearPiecesFunction(
678  ffd1, tdm1, 0, &intercept, &slope
679  );
680  arma::Col<double> lpf1t("-1 0 1 2 3 4 5 6 7.5 9");
681  arma::Col<double> lpf1E("-1 0 1 2 2.5 3 3 3 3.5 4");
682  for(size_t i = 0 ; i < lpf1t.n_elem ; ++i){
683  if(std::fabs(lpf1(lpf1t.at(i)) - lpf1E.at(i)) > eps) return false;
684  if(std::fabs(lpf2(lpf1t.at(i)) - lpf1E.at(i)) > eps) return false;
685  }
686 
687  // Validate ParametricLinearPiecesFunction
689  DiffDesignMatrixInterpolator::makeParametricLinearPiecesFunction(
690  ffd2,
691  tdm2.getX()
692  );
694  DiffDesignMatrixInterpolator::makeParametricLinearPiecesFunction(
695  ffd2, tdm2
696  );
697  arma::Mat<double> plpf1E(
698  "-1 2.5;"
699  "0 3;"
700  "1 3.5;"
701  "2 4;"
702  "2.5 3.5;"
703  "3 3;"
704  "3 4.5;"
705  "3 6;"
706  "3.5 5.5;"
707  "4 5"
708  );
709  for(size_t i = 0 ; i < lpf1t.n_elem ; ++i){
710  arma::Col<double> plpf1y = plpf1(lpf1t.at(i));
711  arma::Col<double> plpf2y = plpf2(lpf1t.at(i));
712  for(size_t j = 0 ; j < 2 ; ++j){
713  if(std::fabs(plpf1y.at(j) - plpf1E.at(i, j)) > eps) return false;
714  if(std::fabs(plpf2y.at(j) - plpf1E.at(i, j)) > eps) return false;
715  }
716  }
717 
718  // Validate FixedIterativeEulerMethod
719  arma::Col<double> fiem1t("0 1 2 3 4 5 6 7.5 9");
720  arma::Col<double> fiem1E("0 1 2 2.5 3 3 3 3.5 4");
721  arma::Col<double> ySamples1(
722  tdm1.getColumnCopy(0).subvec(0, tdm1.getNumRows()-2)
723  );
724  arma::Col<double> ySamples2;
725  arma::Col<double> dydtSamples1( // Vector of known derivative samples
726  ffd1.getA().col(0)
727  );
728  arma::Col<double> dydtSamples2;
730  ffd1.getTimeVector(),
731  dydtSamples1,
732  0
733  );
736  DiffDesignMatrixInterpolator::makeFixedIterativeEulerMethod(
737  ffd1,
738  ySamples1,
739  clsf1
740  );
742  DiffDesignMatrixInterpolator::makeFixedIterativeEulerMethod(
743  ffd1,
744  tdm1,
745  0,
746  &ySamples2,
747  &clsf2,
748  &dydtSamples2
749  );
750  for(size_t i = 0 ; i < fiem1t.n_elem ; ++i){
751  double const h = (i==0) ? 0 : fiem1t.at(i) - fiem1t.at(i-1);
752  if(std::fabs(fiem1(h) - fiem1E.at(i)) > eps) return false;
753  if(std::fabs(fiem2(h) - fiem1E.at(i)) > eps) return false;
754  }
755 
756  // Validate ParametricFixedIterativeEulerMethod
757  arma::Mat<double> pfiem1E(
758  "0 3;"
759  "1 3.5;"
760  "2 4;"
761  "2.5 3.5;"
762  "3 3;"
763  "3 4.5;"
764  "3 6;"
765  "3.5 5.5;"
766  "4 5"
767  );
768  arma::Mat<double> const &pySamples1 = tdm2.getX();
769  arma::Mat<double> const &pdydtSamples1 = ffd2.getA();
771  ffd2.getTimeVector(),
772  pdydtSamples1,
773  0
774  );
776  DiffDesignMatrixInterpolator::makeFixedParametricIterativeEulerMethod(
777  ffd2,
778  pySamples1,
779  pclsf1
780  );
783  DiffDesignMatrixInterpolator::makeFixedParametricIterativeEulerMethod(
784  ffd2,
785  tdm2,
786  &pclfs2
787  );
788 
789  for(size_t i = 0 ; i < fiem1t.n_elem ; ++i){
790  double const h = (i==0) ? 0 : fiem1t.at(i) - fiem1t.at(i-1);
791  if(arma::any(arma::abs(fpiem1(h) - pfiem1E.row(i).as_col()) > eps))
792  return false;
793  if(arma::any(arma::abs(fpiem2(h) - pfiem1E.row(i).as_col()) > eps))
794  return false;
795  }
796 
797  // On passed return true
798  return true;
799 }
800 
801 }
BaseTest class.
Definition: BaseTest.h:20
Fluxionum test.
Definition: FluxionumTest.h:31
bool run() override
Definition: FluxionumTest.h:92
std::string testDir
The directory where test files are located.
Definition: FluxionumTest.h:42
bool testDesignMatrixBuilding()
Test the building of different design matrices.
Definition: FluxionumTest.h:115
bool testDesignFunctions()
Test functions from design matrices.
Definition: FluxionumTest.h:636
bool testUnivariateNewtonRaphsonMinimization()
Test univariate Newton-Raphson minimization.
Definition: FluxionumTest.h:104
FluxionumTest(std::string testDir)
Fluxionum test constructor.
Definition: FluxionumTest.h:49
bool testDesignMatrixMethods()
Test the methods of DesignMatrix.
Definition: FluxionumTest.h:327
bool testDiffDesignMatrix()
Test the generation of differential design matrices.
Definition: FluxionumTest.h:502
string const & getColumnName(size_t const j) const
Obtain the name of the -th column.
Definition: AbstractDesignMatrix.h:111
bool hasColumnNames() const
Check whether there are available column names for the AbstractDesignMatrix (true) or not (false)
Definition: AbstractDesignMatrix.h:104
vector< string > const & getColumnNames() const
Obtain a constant/read reference to the column names.
Definition: AbstractDesignMatrix.h:126
Closest lesser sample function.
Definition: ClosestLesserSampleFunction.h:28
This class represents a DesignMatrix .
Definition: DesignMatrix.h:41
virtual void mergeInPlace(DesignMatrix const &dm)
Merge given DesignMatrix into this DesignMatrix.
Definition: DesignMatrix.h:119
size_t getNumColumns() const override
Obtain the number of columns of the DesignMatrix .
Definition: DesignMatrix.h:244
arma::Col< T > getColumnCopy(size_t const j) const
Like DesignMatrix::getColumn(size_t const) but returning a copy by value instead of a view-like refer...
Definition: DesignMatrix.h:279
size_t getNumRows() const override
Obtain the number of rows of the DesignMatrix .
Definition: DesignMatrix.h:237
arma::Mat< T > const & getX() const
Obtain a constant/read reference to the matrix.
Definition: DesignMatrix.h:219
void swapColumns(vector< long unsigned int > const &indices)
Definition: DesignMatrix.h:125
The heart of a differential design matrix is the idea that the columns are defining the values of var...
Definition: DiffDesignMatrix.h:169
size_t getNumRows() const override
Obtain the number of rows of the DiffDesignMatrix .
Definition: DiffDesignMatrix.h:443
string const & getTimeName() const
Obtain the name of the time attribute.
Definition: DiffDesignMatrix.h:469
size_t getNumColumns() const override
Obtain the number of columns of the DiffDesignMatrix .
Definition: DiffDesignMatrix.h:450
arma::Mat< VarType > const & getA() const
Obtain a constant/read reference to the matrix.
Definition: DiffDesignMatrix.h:481
arma::Col< TimeType > const & getTimeVector() const
Obtain the sorted series of time values DiffDesignMatrix::t.
Definition: DiffDesignMatrix.h:475
DiffDesignMatrixType getDiffType() const
Obtain the type of differential of the DiffDesignMatrix.
Definition: DiffDesignMatrix.h:463
Fixed iterative Euler method.
Definition: FixedIterativeEulerMethod.h:25
Definition: IndexedDesignMatrix.h:43
void mergeInPlace(DesignMatrix< VarType > const &dm) override
Extend DesignMatrix::mergeInPlace method so the indices are also merged.
Definition: IndexedDesignMatrix.h:269
vector< IndexType > const & getIndices() const
Obtain a constant/read reference to the vector of indices.
Definition: IndexedDesignMatrix.h:287
string const getIndexName() const
Obtain the name of the index attribute.
Definition: IndexedDesignMatrix.h:293
Linear pieces function.
Definition: LinearPiecesFunction.h:28
Parametric linear pieces function.
Definition: ParametricLinearPiecesFunction.h:35
This class represents a DesignMatrix where each row satisfy that its elements can be placed into a co...
Definition: TemporalDesignMatrix.h:74
virtual size_t slopeFilter(VarType const tau)
Apply a filter based on the first derivative understood as the slope.
Definition: TemporalDesignMatrix.h:391
string const & getTimeName() const
Obtain the name of the time attribute.
Definition: TemporalDesignMatrix.h:424
DiffDesignMatrix< TimeType, VarType > toDiffDesignMatrix(DiffDesignMatrixType diffType=DiffDesignMatrixType::FORWARD_FINITE_DIFFERENCES, bool const sort=true) const
Build a DiffDesignMatrix from the TemporalDesignMatrix.
arma::Col< TimeType > const & getTimeVector() const
Obtain a constant/read reference to the time vector .
Definition: TemporalDesignMatrix.h:418
virtual void sortByTime()
Sort all rows of matrix and all components of vector so , considering the row is associated with ...
Definition: TemporalDesignMatrix.h:351
void mergeInPlace(DesignMatrix< VarType > const &dm) override
Extend DesignMatrix::mergeInPlace method so the time values are also merged.
Definition: TemporalDesignMatrix.h:315
virtual void shiftTime(TimeType const x)
Add given time to all time values in vector .
Definition: TemporalDesignMatrix.h:343
Implementation of univariate Newton-Raphson minimizer.
Definition: UnivariateNewtonRaphsonMinimizer.h:23
IT argmin(IT x) override
Implementation of the univariate Newton-Raphson minimization.