ARTS  2.4.0(git:4fb77825)
test_oem.cc
Go to the documentation of this file.
1 
9 #include <stdlib.h>
10 #include <time.h>
11 #include <cmath>
12 #include <fstream>
13 #include <iomanip>
14 #include <iostream>
15 #include <string>
16 #include "engine.h"
17 #include "lin_alg.h"
18 #include "matpackII.h"
19 #include "matrix.h"
20 #include "oem.h"
21 #include "test_utils.h"
22 #include "unistd.h"
23 
24 using std::abs;
25 using std::cout;
26 using std::endl;
27 using std::min;
28 using std::ofstream;
29 using std::setw;
30 using std::string;
31 
32 string source_dir = SOURCEDIR;
33 string atmlab_dir = ATMLABDIR;
34 
35 // Forward declarations.
36 
37 void write_matrix(ConstMatrixView A, const char* fname);
38 
40 
51 class LinearModel {
52  public:
54  LinearModel(Index m_, Index n_) : m(m_), n(n_) {}
55 
63  : m(J_.nrows()), n(J_.ncols()), J(J_), y0(y0_) {}
64 
66  const OEMMatrix& Jacobian(const ConstVectorView& xi) { return J; }
67 
69  OEMVector evaluate(const ConstVectorView& xi) {
70  OEMVector y;
71  y.resize(m);
72  mult(y, J, xi);
73  y += y0;
74  return y;
75  }
76 
77  const Index m, n;
78 
79  private:
80  OEMMatrix J;
81  OEMVector y0;
82 };
83 
85 
92  public:
94 
103  QuadraticModel(Index m_, Index n_) : m(m_), n(n_) {
104  char fname[40];
105 
106  J = Matrix(m, n, 0);
107  Hessians = new OEMMatrix[m];
108  random_fill_matrix(J, 1.0, false);
109  sprintf(fname, "J_t.txt");
110  write_matrix(J, fname);
111 
112  for (Index i = 0; i < m; i++) {
113  Hessians[i] = OEMMatrix{};
114  Hessians[i].resize(n, n);
115  random_fill_matrix_symmetric(Hessians[i], 0.0001, false);
116  sprintf(fname, "H_%d_t.txt", (int)i);
117  write_matrix(Hessians[i], fname);
118  }
119  }
120 
122  ~QuadraticModel() { delete[] Hessians; }
123 
125  OEMMatrix Jacobian(const ConstVectorView& xi) {
126  OEMMatrix Ki{};
127  Ki.resize(m, n);
128 
129  for (Index i = 0; i < m; i++) {
130  mult(static_cast<MatrixView>(Ki)(i, Joker()), Hessians[i], xi);
131  }
132 
133  Ki *= 0.5;
134  Ki += J;
135  return Ki;
136  }
137 
139  OEMVector evaluate(const OEMVector& xi) {
140  OEMVector yi{};
141  yi.resize(m);
142  OEMMatrix Ki{};
143  Ki.resize(m, n);
144 
145  for (Index i = 0; i < n; i++) {
146  mult(static_cast<MatrixView>(Ki)(i, Joker()), Hessians[i], xi);
147  }
148 
149  Ki *= 0.5;
150  Ki += J;
151  mult(yi, Ki, xi);
152  return yi;
153  }
154 
155  const unsigned int m, n;
156 
157  private:
158  OEMMatrix J;
159  OEMMatrix* Hessians;
160 };
161 
163 
170 void write_matrix(ConstMatrixView A, const char* filename) {
171  Index m = A.nrows();
172  Index n = A.ncols();
173 
174  ofstream ofs(filename, ofstream::out);
175 
176  for (Index i = 0; i < m; i++) {
177  for (Index j = 0; j < (n - 1); j++) {
178  ofs << std::setprecision(40) << A(i, j) << " ";
179  }
180  ofs << A(i, n - 1);
181  ofs << endl;
182  }
183  ofs.close();
184 }
185 
187 
192 void write_vector(ConstVectorView v, const char* filename) {
193  Index n = v.nelem();
194 
195  ofstream ofs(filename, ofstream::out);
196 
197  for (Index i = 0; i < n; i++) {
198  ofs << std::setprecision(20) << v[i] << endl;
199  }
200  ofs.close();
201 }
202 
204 
221  VectorView xa,
222  MatrixView Se,
223  MatrixView Sx) {
224  random_fill_vector(y, 10, false);
225  random_fill_vector(xa, 10, false);
226 
227  random_fill_matrix(Se, 1.0, false);
228  Matrix tmp(Se);
229  // Make sure Se is positive semi-definite.
230  mult(Se, transpose(tmp), tmp);
231 
232  random_fill_matrix_symmetric(Sx, 1.0, false);
233  tmp = Sx;
234  // Make sure Sx is positive semi-definite.
235  mult(Sx, transpose(tmp), tmp);
236 }
237 
239 
246 
248 
258 Index run_test_matlab(Engine* eng, string filename) {
259  mxArray* t;
260  Index time;
261 
262  // Run test.
263  string cmd = "run('" + filename + "');";
264  engEvalString(eng, cmd.c_str());
265 
266  // Get execution time from matlab.
267  t = engGetVariable(eng, "t");
268  time = (Index)((Numeric*)mxGetData(t))[0];
269  return time;
270 }
271 
273 
285 Index run_oem_matlab(VectorView x, MatrixView G, Engine* eng, string filename) {
286  Index n = G.nrows();
287  Index m = G.ncols();
288  Index time;
289  mxArray *x_m, *G_m, *t;
290 
291  // Run test.
292  string cmd = "run('" + filename + "');";
293  engEvalString(eng, cmd.c_str());
294 
295  // Read out results.
296  x_m = engGetVariable(eng, "x");
297  G_m = engGetVariable(eng, "G");
298 
299  for (Index i = 0; i < n; i++) {
300  x[i] = ((Numeric*)mxGetData(x_m))[i];
301 
302  for (Index j = 0; j < m; j++) {
303  G(i, j) = ((Numeric*)mxGetData(G_m))[j * n + i];
304  }
305  }
306 
307  // Get execution time from matlab.
308  t = engGetVariable(eng, "t");
309  time = (Index)((Numeric*)mxGetData(t))[0];
310  return time;
311 }
312 
314 
321 void setup_test_environment(Engine*& eng) {
322  // swith to test folder
323  string cmd;
324  cmd = source_dir + "/test_oem_files";
325  int out = chdir(cmd.c_str());
326  (void)out;
327 
328  // Start MATLAB and try to initialize atmlab package.
329  string atmlab_init = "run('" + atmlab_dir + "/atmlab/atmlab_init.m');";
330 
331  eng = engOpen(NULL);
332 
333  engEvalString(eng, atmlab_init.c_str());
334  cmd = "cd('" + source_dir + "/test_oem_files');";
335  engEvalString(eng, cmd.c_str());
336 }
337 
339 
347 void run_plot_script(Engine* eng, string filename, string title) {
348  string cmd = "filename = '" + filename + "'";
349  engEvalString(eng, cmd.c_str());
350  cmd = "plot_title = '" + title + "'";
351  engEvalString(eng, cmd.c_str());
352  engEvalString(eng, "run('make_plot.m');");
353 }
354 
356 
360 void tidy_up_test_environment(Engine* eng) {
361  int out = system("rm *_t.txt");
362  (void)out;
363 
364  engEvalString(eng, "close()");
365 }
366 
368 
381 void benchmark_inv(Engine* eng, Index n0, Index n1, Index ntests) {
382  Index step = (n1 - n0) / (ntests - 1);
383  Index n = n0;
384 
385  ofstream ofs("times_inv.txt", ofstream::out);
386  ofs << "#" << setw(4) << "n" << setw(10) << "BLAS";
387  ofs << setw(10) << "arts" << setw(10) << "Matlab" << endl;
388 
389  cout << endl << "N TIMES N MATRIX INVERSION" << endl << endl;
390  cout << setw(5) << "n" << setw(10) << "BLAS" << setw(10);
391  cout << setw(10) << "arts" << setw(10) << "Matlab" << endl;
392 
393  for (Index i = 0; i < ntests; i++) {
394  Matrix A(n, n), B(n, n);
395 
396  random_fill_matrix(A, 100, false);
397  write_matrix(A, "A_t.txt");
398 
399  Index t, t1, t2, t_blas, t1_blas, t2_blas, t_m;
400 
401  t1 = clock();
402  inv(B, A);
403  t2 = clock();
404  t = (t2 - t1) * 1000 / CLOCKS_PER_SEC;
405 
406  t1_blas = clock();
407  inv(B, A);
408  t2_blas = clock();
409  t_blas = (t2_blas - t1_blas) * 1000 / CLOCKS_PER_SEC;
410 
411  t_m = run_test_matlab(eng, "test_inv.m");
412 
413  ofs << setw(5) << n << setw(10) << t_blas << setw(10);
414  ofs << t << setw(10) << t_m << endl;
415  cout << setw(5) << n << setw(10) << t_blas << setw(10);
416  cout << t << setw(10) << t_m << endl;
417 
418  n += step;
419  }
420  cout << endl << endl;
421 
422  // Tidy up
423  ofs.close();
424  run_plot_script(eng, "times_inv.txt", "Matrix Inversion");
425 }
426 
428 
441 void benchmark_mult(Engine* eng, Index n0, Index n1, Index ntests) {
442  Index step = (n1 - n0) / (ntests - 1);
443  Index n = n0;
444 
445  ofstream ofs("times_mult.txt", ofstream::out);
446  ofs << "#" << setw(4) << "n" << setw(10) << "BLAS";
447  ofs << setw(10) << "arts" << setw(10) << "Matlab" << endl;
448 
449  cout << endl << "N TIMES N MATRIX MULTIPLICATION" << endl << endl;
450  cout << setw(5) << "n" << setw(10) << "BLAS" << setw(10);
451  cout << setw(10) << "arts" << setw(10) << "Matlab" << endl;
452 
453  for (Index i = 0; i < ntests; i++) {
454  Matrix A(n, n), B(n, n);
455 
456  random_fill_matrix(A, 100, false);
457  write_matrix(A, "A_t.txt");
458 
459  Index t, t1, t2, t_blas, t1_blas, t2_blas, t_m;
460 
461  t1 = clock();
462  mult_general(B, A, A);
463  t2 = clock();
464  t = (t2 - t1) * 1000 / CLOCKS_PER_SEC;
465 
466  t1_blas = clock();
467  mult(B, A, A);
468  t2_blas = clock();
469  t_blas = (t2_blas - t1_blas) * 1000 / CLOCKS_PER_SEC;
470 
471  t_m = run_test_matlab(eng, "test_mult.m");
472 
473  ofs << setw(5) << n << setw(10) << t_blas << setw(10) << t << setw(10);
474  ofs << t_m << endl;
475  cout << setw(5) << n << setw(10) << t_blas << setw(10) << t << setw(10);
476  cout << t_m << endl;
477 
478  n += step;
479  }
480  cout << endl << endl;
481 
482  // Tidy up
483  ofs.close();
484  run_plot_script(eng, "times_mult.txt", "Matrix Multiplication");
485 }
486 
488 
499 void benchmark_oem_linear(Engine* eng, Index n0, Index n1, Index ntests) {
500  Index step = (n1 - n0) / (ntests - 1);
501  Index n = n0;
502 
503  ofstream ofs("times_linear.txt", ofstream::out);
504 
505  ofs << "#" << setw(4) << "n" << setw(10) << "Matlab";
506  ofs << setw(10) << "C++" << endl;
507 
508  cout << endl << "LINEAR OEM" << endl << endl;
509  cout << setw(5) << "n" << setw(10) << "C++" << setw(10);
510  cout << "Matlab" << setw(20) << "Max. Rel. Error" << endl;
511 
512  // Run tests.
513  for (Index i = 0; i < ntests; i++) {
514  Vector x_nform(n), x_mform(n), x_m(n), y(n), yf(n), xa(n), x_norm(n),
515  zero(n);
516  Matrix J(n, n), Se(n, n), Sa(n, n), SeInv(n, n), SxInv(n, n), G_nform(n, n),
517  G_mform(n, n), G_m(n, n);
518 
519  zero = 0.0;
520 
521  generate_test_data(y, xa, Se, Sa);
523  LinearModel K(J, zero);
524 
525  write_vector(xa, "xa_t.txt");
526  write_vector(y, "y_t.txt");
527  write_matrix(J, "J_t.txt");
528  write_matrix(Se, "Se_t.txt");
529  write_matrix(Sa, "Sa_t.txt");
530 
531  Index t, t1, t2, t_m;
532 
533  inv(SeInv, Se);
534  inv(SxInv, Sa);
535 
536  // n-form
537  t1 = clock();
538  mult(yf, J, xa);
539  // oem_linear_nform( x_nform, G_nform, J, yf, cost_x, cost_y, K,
540  // xa, x_norm, y, SeInv, SxInv, cost_x, false );
541  t2 = clock();
542  t = (t2 - t1) * 1000 / CLOCKS_PER_SEC;
543 
544  // m-form
545  mult(yf, J, xa);
546  // oem_linear_mform( x_mform, G_mform, xa, yf, y, J, Se, Sa );
547 
548  // Matlab
549  t_m = run_oem_matlab(x_m, G_m, eng, "test_oem");
550 
551  ofs << setw(5) << n << setw(10) << t << setw(10) << 42; // Dummy column
552  ofs << setw(10) << t_m << endl;
553  cout << setw(5) << n << setw(10) << t << setw(10) << t_m;
554  cout << endl << endl;
555 
556  n += step;
557  }
558  cout << endl << endl;
559 
560  // Tidy up
561  ofs.close();
562  run_plot_script(eng, "times_linear.txt", "Linear OEM");
563 }
564 
566 
575 void test_oem_linear(Engine* eng, Index m, Index n, Index ntests) {
576  cout << "Testing linear OEM: m = " << m << ", n = ";
577  cout << n << ", ntests = " << ntests << endl << endl;
578 
579  cout << "Test No. " << setw(15) << "Gauss-Newton";
580  cout << setw(20) << "Levenberg-Marquardt" << setw(15) << "Gain Matrix"
581  << endl;
582 
583  // Run tests.
584  for (Index i = 0; i < ntests; i++) {
585  // Create linear forward model and test data.
586  Matrix J(m, n);
588  OEMVector y, xa;
589  y.resize(m);
590  xa.resize(n);
591  OEMMatrix Se, Sa;
592  Se.resize(m, m);
593  Sa.resize(n, n);
594  generate_test_data(y, xa, Se, Sa);
595 
596  Vector zero(m);
597  zero = 0.0;
598  LinearModel K(J, zero);
599 
600  Pre pre(Sa);
601  Std std{};
602  LM_D lm(Sa);
603  GN gn(1e-12, 100);
604  LM_Pre_D lm_pre(Sa, Std_Pre(std, pre));
605  GN_Pre gn_pre(1e-12, 100, Std_Pre(std, pre));
606 
607  OEM_D_D<LinearModel> oem_std(K, xa, Sa, Se);
608  OEM_NFORM_D_D<LinearModel> oem_nform(K, xa, Sa, Se);
609  OEM_MFORM_D_D<LinearModel> oem_mform(K, xa, Sa, Se);
610 
611  // Write data for Matlab.
612 
613  Vector x_norm(n);
614  for (Index j = 0; j < n; j++) {
615  x_norm[j] = sqrt(Sa(j, j));
616  }
617 
618  write_vector(xa, "xa_t.txt");
619  write_vector(y, "y_t.txt");
620  write_matrix(J, "J_t.txt");
621  write_matrix(Se, "Se_t.txt");
622  write_matrix(Sa, "Sa_t.txt");
623 
624  OEMVector x_std, x_nform, x_mform;
625  x_std.resize(n);
626  x_nform.resize(n);
627  x_mform.resize(n);
628 
629  oem_std.compute(x_std, y, gn);
630  oem_nform.compute(x_nform, y, gn);
631  oem_mform.compute(x_mform, y, gn);
632 
633  OEMVector x_std_lm;
634  x_std_lm.resize(n);
635  oem_std.compute(x_std_lm, y, lm);
636 
637  OEMVector x_std_norm_gn, x_std_norm_lm;
638  x_std_norm_gn.resize(n);
639  x_std_norm_lm.resize(n);
640  oem_std.compute(x_std_norm_gn, y, gn_pre);
641  oem_std.compute(x_std_norm_lm, y, lm_pre);
642 
643  OEMMatrix G_std, G_nform, G_mform, G_norm;
644  G_std = oem_std.gain_matrix(x_std);
645  G_nform = oem_nform.gain_matrix(x_std);
646  G_mform = oem_mform.gain_matrix(x_std);
647 
648  Vector x_m(n);
649  Matrix G_m(m, n);
650  run_oem_matlab(x_m, G_m, eng, "test_oem");
651 
652  Numeric err, err_G, err_lm, err_norm;
653  err = get_maximum_error(x_std, x_m, true);
654  err = std::max(err, get_maximum_error(x_nform, x_m, true));
655  err = std::max(err, get_maximum_error(x_mform, x_m, true));
656 
657  err_lm = get_maximum_error(x_std_lm, x_m, true);
658 
659  err_norm = get_maximum_error(x_std_norm_lm, x_m, true);
660  err_norm = std::max(err_norm, get_maximum_error(x_std_norm_lm, x_m, true));
661 
662  err_G = get_maximum_error(G_std, G_m, true);
663  err_G = std::max(err, get_maximum_error(G_nform, G_m, true));
664  err_G = std::max(err, get_maximum_error(G_nform, G_m, true));
665 
666  cout << setw(8) << i + 1;
667  cout << setw(15) << err << setw(20) << err_lm << setw(15) << err_G;
668  cout << setw(15) << err_norm << endl;
669  }
670  cout << endl;
671 }
672 
674 
686 void test_oem_gauss_newton(Engine* eng, Index m, Index n, Index ntests) {
687  cout << "Testing Gauss-Newton OEM: m = " << m << ", n = ";
688  cout << n << ", ntests = " << ntests << endl << endl;
689 
690  cout << "Test No. " << setw(15) << "Standard" << setw(15) << "Normalized";
691  cout << setw(15) << "CG Solver" << setw(15) << "CG Normalized" << endl;
692 
693  for (Index i = 0; i < ntests; i++) {
694  // Generate random quadratic model.
695 
696  OEMMatrix Se, Sa, SeInv, SaInv;
697  Se.resize(m, m);
698  Sa.resize(n, n);
699  SeInv.resize(m, m);
700  SaInv.resize(n, n);
701  OEMVector xa;
702  xa.resize(n);
703  Vector y0(m), x0(n);
704 
705  QuadraticModel K(m, n);
706  generate_test_data(y0, xa, Se, Sa);
707  x0 = xa;
708  add_noise(x0, 0.1);
709  y0 = K.evaluate(x0);
710 
711  inv(SaInv, Sa);
712  inv(SeInv, Se);
713 
714  // Solvers.
715  Pre pre(Sa);
716  Std std{};
717  CG cg(1e-12);
718 
719  // Optimization Methods.
720  GN gn(std);
721  GN_Pre gn_pre(Std_Pre(std, pre));
722  GN_CG gn_cg(cg);
723  GN_CG_Pre gn_cg_pre(CG_Pre(cg, pre));
724 
725  PrecisionMatrix Pe(SeInv), Pa(SaInv);
726  OEM_D_D<QuadraticModel> map(K, xa, Sa, Se);
727  OEM_D_D<QuadraticModel> map_prec(K, xa, Pa, Pe);
728 
729  // Write data to disk.
730  Vector x_norm(n);
731  for (Index j = 0; j < n; j++) {
732  x_norm[j] = sqrt(abs(Sa(j, j)));
733  }
734 
735  write_vector(xa, "xa_t.txt");
736  write_vector(y0, "y_t.txt");
737  write_matrix(Se, "Se_t.txt");
738  write_matrix(Sa, "Sa_t.txt");
739 
740  OEMVector x, x_pre, x_cg, x_cg_pre;
741  map.compute(x, y0, gn);
742  map.compute(x_pre, y0, gn_pre);
743  map.compute(x_cg, y0, gn_cg);
744  map.compute(x_cg_pre, y0, gn_cg_pre);
745 
746  Vector x_m(n);
747  Matrix G_m(m, n);
748  run_oem_matlab(x_m, G_m, eng, "test_oem_gauss_newton");
749 
750  Numeric e1, e2, e3, e4;
751  e1 = get_maximum_error(x, x_m, true);
752  e2 = get_maximum_error(x_pre, x_m, true);
753  e3 = get_maximum_error(x_cg, x_m, true);
754  e4 = get_maximum_error(x_cg_pre, x_m, true);
755 
756  cout << setw(9) << i + 1;
757  cout << setw(15) << e1 << setw(15) << e2 << setw(15) << e3;
758  cout << setw(15) << e4 << std::endl;
759 
760  map_prec.compute(x, y0, gn);
761  map_prec.compute(x_pre, y0, gn_pre);
762  map_prec.compute(x_cg, y0, gn_cg);
763  map_prec.compute(x_cg_pre, y0, gn_cg_pre);
764 
765  e1 = get_maximum_error(x, x_m, true);
766  e2 = get_maximum_error(x_pre, x_m, true);
767  e3 = get_maximum_error(x_cg, x_m, true);
768  e4 = get_maximum_error(x_cg_pre, x_m, true);
769 
770  cout << setw(9) << i + 1;
771  cout << setw(15) << e1 << setw(15) << e2 << setw(15) << e3;
772  cout << setw(15) << e4 << std::endl;
773  }
774  cout << endl;
775 }
776 
778 
790 void test_oem_levenberg_marquardt(Engine* eng, Index m, Index n, Index ntests) {
791  cout << "Testing Levenberg-Marquardt OEM: m = " << m << ", n = ";
792  cout << n << ", ntests = " << ntests << endl << endl;
793 
794  cout << "Test No. " << setw(15) << "Standard" << setw(15) << "Normalized";
795  cout << setw(15) << "CG Solver" << setw(15) << "CG Normalized" << endl;
796 
797  for (Index i = 0; i < ntests; i++) {
798  // Generate random quadratic model.
799 
800  OEMMatrix Se, Sa, SeInv, SaInv;
801  Se.resize(m, m);
802  Sa.resize(n, n);
803  SeInv.resize(m, m);
804  SaInv.resize(n, n);
805  OEMVector xa;
806  xa.resize(n);
807  Vector y0(m), x0(n);
808 
809  QuadraticModel K(m, n);
810  generate_test_data(y0, xa, Se, Sa);
811  x0 = xa;
812  add_noise(x0, 0.1);
813  y0 = K.evaluate(x0);
814 
815  inv(SaInv, Sa);
816  inv(SeInv, Se);
817 
818  // Solvers.
819  Pre pre(Sa);
820  Std std{};
821  CG cg(1e-12);
822 
823  // Optimization Methods.
824  LM_D lm(SaInv, std);
825  LM_Pre_D lm_pre(SaInv, Std_Pre(std, pre));
826  LM_CG_D lm_cg(SaInv, cg);
827  LM_CG_Pre_D lm_cg_pre(SaInv, CG_Pre(cg, pre));
828 
829  PrecisionMatrix Pe(SeInv), Pa(SaInv);
830  OEM_D_D<QuadraticModel> map(K, xa, Sa, Se);
831  OEM_PD_PD<QuadraticModel> map_prec(K, xa, Pa, Pe);
832 
833  // Write data to disk.
834  Vector x_norm(n);
835  for (Index j = 0; j < n; j++) {
836  x_norm[j] = sqrt(abs(Sa(j, j)));
837  }
838 
839  write_vector(xa, "xa_t.txt");
840  write_vector(y0, "y_t.txt");
841  write_matrix(Se, "Se_t.txt");
842  write_matrix(Sa, "Sa_t.txt");
843 
844  OEMVector x, x_pre, x_cg, x_cg_pre;
845  map.compute(x, y0, lm);
846  map.compute(x_pre, y0, lm_pre);
847  map.compute(x_cg, y0, lm_cg);
848  map.compute(x_cg_pre, y0, lm_cg_pre);
849 
850  Vector x_m(n);
851  Matrix G_m(m, n);
852  run_oem_matlab(x_m, G_m, eng, "test_oem_gauss_newton");
853  Vector x_m2(n);
854  Matrix G_m2(m, n);
855  run_oem_matlab(x_m2, G_m2, eng, "test_oem_levenberg_marquardt");
856 
857  Numeric e1, e2, e3, e4;
858  e1 = get_maximum_error(x, x_m, true);
859  e2 = get_maximum_error(x_pre, x_m, true);
860  e3 = get_maximum_error(x_cg, x_m2, true);
861  e4 = get_maximum_error(x_cg_pre, x_m2, true);
862 
863  cout << setw(9) << i + 1;
864  cout << setw(15) << e1 << setw(15) << e2 << setw(15) << e3;
865  cout << setw(15) << e4 << std::endl;
866 
867  map_prec.compute(x, y0, lm);
868  map_prec.compute(x_pre, y0, lm_pre);
869  map_prec.compute(x_cg, y0, lm_cg);
870  map_prec.compute(x_cg_pre, y0, lm_cg_pre);
871 
872  e1 = get_maximum_error(x, x_m, true);
873  e2 = get_maximum_error(x_pre, x_m, true);
874  e3 = get_maximum_error(x_cg, x_m2, true);
875  e4 = get_maximum_error(x_cg_pre, x_m2, true);
876 
877  cout << setw(9) << i + 1;
878  cout << setw(15) << e1 << setw(15) << e2 << setw(15) << e3;
879  cout << setw(15) << e4 << std::endl;
880  }
881  cout << endl;
882 }
883 
885 
898  cout << "Testing Gauss-Newton OEM: m = " << m << ", n = ";
899  cout << n << ", ntests = " << ntests << endl << endl;
900 
901  cout << "Test No. " << setw(15) << "Standard" << setw(15) << "Normalized";
902  cout << setw(15) << "CG Solver" << setw(15) << "CG Normalized" << endl;
903 
904  for (Index i = 0; i < ntests; i++) {
905  // Generate random quadratic model.
906  OEMVector xa;
907  xa.resize(n);
908  Vector y0(m), x0(n);
909 
910  QuadraticModel K(m, n);
911  random_fill_vector(y0, 10, false);
912  random_fill_vector(xa, 10, false);
913  x0 = xa;
914  add_noise(x0, 0.1);
915  y0 = K.evaluate(x0);
916 
917  // Create sparse covariance matrices.
918  Sparse Se_sparse(m, m), Sa_sparse(n, n);
919  id_mat(Se_sparse);
920  //random_fill_matrix( Se_sparse, 1.0, false );
921  id_mat(Sa_sparse);
922  //random_fill_matrix( Sa_sparse, 1.0, false );
923  OEMMatrix Se = (Matrix)Se_sparse;
924  OEMMatrix Sa = (Matrix)Sa_sparse;
925  ArtsSparse Se_arts(Se_sparse);
926  ArtsSparse Sa_arts(Sa_sparse);
927  OEMSparse Se_map(Se_arts), Sa_map(Sa_arts);
928 
929  // Solvers.
930  Pre pre(Sa);
931  Std std{};
932  CG cg(1e-12);
933 
934  // Optimization Methods.
935  GN gn(std);
936  GN_Pre gn_pre(Std_Pre(std, pre));
937  GN_CG gn_cg(cg);
938  GN_CG_Pre gn_cg_pre(CG_Pre(cg, pre));
939 
940  PrecisionMatrix Pe(Se);
941  PrecisionMatrix Pa(Sa);
942  PrecisionSparse Pe_sparse(Se_map);
943  PrecisionSparse Pa_sparse(Sa_map);
944  OEM_PD_PD<QuadraticModel> map(K, xa, Pa, Pe);
945  OEM_PS_PS<QuadraticModel> map_sparse(K, xa, Pa_sparse, Pe_sparse);
946 
947  OEMVector x, x_pre, x_cg, x_cg_pre;
948  map.compute(x, y0, gn);
949  map.compute(x_pre, y0, gn_pre);
950  map.compute(x_cg, y0, gn_cg);
951  map.compute(x_cg_pre, y0, gn_cg_pre);
952 
953  OEMVector x_sparse, x_pre_sparse, x_cg_sparse, x_cg_pre_sparse;
954  map_sparse.compute(x_sparse, y0, gn);
955  map_sparse.compute(x_pre_sparse, y0, gn_pre);
956  map_sparse.compute(x_cg_sparse, y0, gn_cg);
957  map_sparse.compute(x_cg_pre_sparse, y0, gn_cg_pre);
958 
959  Numeric e1, e2, e3, e4;
960  e1 = get_maximum_error(x, x_sparse, true);
961  e2 = get_maximum_error(x_pre, x_pre_sparse, true);
962  e3 = get_maximum_error(x_cg, x_cg_sparse, true);
963  e4 = get_maximum_error(x_cg_pre, x_cg_pre_sparse, true);
964 
965  cout << setw(9) << i + 1;
966  cout << setw(15) << e1 << setw(15) << e2 << setw(15) << e3;
967  cout << setw(15) << e4 << std::endl;
968  }
969  cout << endl;
970 }
971 
973 
986 
987 int main() {
988  // Set up the test environment.
989  //Engine * eng;
990  //setup_test_environment( eng );
991 
992  //test_oem_linear(eng, 5, 5, 10);
993  //test_oem_gauss_newton(eng, 100, 100, 10);
994  //test_oem_levenberg_marquardt(eng, 100, 100, 10);
995  test_oem_gauss_newton_sparse(100, 50, 10);
996 
997  // Run tests and benchmarks.
998  //test_oem_linear( eng, 50, 50, 10 );
999  //test_oem_gauss_newton( eng, 50, 50, 10 );
1000  //test_oem_levenberg_marquardt( eng, 50, 50, 10 );
1001 
1002  //benchmark_inv( eng, 100, 2000, 16);
1003  //benchmark_mult( eng, 100, 2000, 16);
1004  //test_oem_linear( eng, 200, 200, 5 );
1005  //benchmark_oem_linear( eng, 100, 2000, 16);
1006 
1007  // Tidy up test environment.
1008  //tidy_up_test_environment(eng);
1009 
1010  // test_oem_gauss_newton_sparse( 100, 100, 10 );
1011  // test_oem_levenberg_marquardt_sparse( 100, 100, 10 );
1012 }
Matrix
The Matrix class.
Definition: matpackI.h:1193
tidy_up_test_environment
void tidy_up_test_environment(Engine *eng)
Tidy up test environment.
Definition: test_oem.cc:360
main
int main()
Definition: test_oem.cc:987
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
run_oem_matlab
Index run_oem_matlab(VectorView x, MatrixView G, Engine *eng, string filename)
Run test script in matlab and return result vector.
Definition: test_oem.cc:285
LinearModel::n
const Index n
Definition: test_oem.cc:77
random_fill_matrix_symmetric
void random_fill_matrix_symmetric(MatrixView A, Numeric range, bool positive)
Generate random, symmetric matrix.
Definition: test_utils.cc:156
test_utils.h
Utility functions for testing.
id_mat
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2240
setup_test_environment
void setup_test_environment(Engine *&eng)
Setup the test environment.
Definition: test_oem.cc:321
inv
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
oem::Matrix
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:34
write_matrix
void write_matrix(ConstMatrixView A, const char *fname)
Write matrix to text file.
Definition: test_oem.cc:170
abs
#define abs(x)
Definition: legacy_continua.cc:20626
LinearModel::m
const Index m
Definition: test_oem.cc:77
benchmark_inv
void benchmark_inv(Engine *eng, Index n0, Index n1, Index ntests)
Matrix inversion benchmark.
Definition: test_oem.cc:381
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
linalg::std
void std(VectorView std, const Vector &y, const ArrayOfVector &ys, const Index start=0, const Index end=-1)
Compute the standard deviation of the ranged ys.
Definition: raw.cc:66
ARTS::Var::xa
Vector xa(Workspace &ws) noexcept
Definition: autoarts.h:7359
Sparse
The Sparse class.
Definition: matpackII.h:60
random_fill_vector
void random_fill_vector(VectorView v, Numeric range, bool positive)
Fill vector with random values.
Definition: test_utils.cc:230
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
benchmark_mult
void benchmark_mult(Engine *eng, Index n0, Index n1, Index ntests)
Matrix multiplication benchmark.
Definition: test_oem.cc:441
QuadraticModel::Jacobian
OEMMatrix Jacobian(const ConstVectorView &xi)
Virtual function of the FowardModel class.
Definition: test_oem.cc:125
x0
#define x0
Definition: lineshapemodel.h:413
QuadraticModel::~QuadraticModel
~QuadraticModel()
Destructor.
Definition: test_oem.cc:122
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
LinearModel::LinearModel
LinearModel(ConstMatrixView J_, ConstVectorView y0_)
Construct a linear model from a given Jacobian J and offset vector y0.
Definition: test_oem.cc:62
QuadraticModel::m
const unsigned int m
Definition: test_oem.cc:155
LinearModel::LinearModel
LinearModel(Index m_, Index n_)
Default Constructor.
Definition: test_oem.cc:54
oem::GN
invlib::GaussNewton< Numeric, Std > GN
OEM Gauss-Newton optimization using normed ARTS QR solver.
Definition: oem.h:169
QuadraticModel
Quadratic Forward Model.
Definition: test_oem.cc:91
oem::CG
NormalizingSolver< Matrix, invlib::ConjugateGradient<> > CG
The invlib CG solver.
Definition: oem.h:166
oem::GN_CG
invlib::GaussNewton< Numeric, CG > GN_CG
Gauss-Newton (GN) optimization using normed CG solver.
Definition: oem.h:171
min
#define min(a, b)
Definition: legacy_continua.cc:20628
random_fill_matrix
void random_fill_matrix(MatrixView A, Numeric range, bool positive)
Fill matrix with random values.
Definition: test_utils.cc:59
test_oem_levenberg_marquardt_sparse
void test_oem_levenberg_marquardt_sparse(Index m, Index n, Index ntests)
Test sparse non-linear OEM.
Definition: test_oem.cc:985
LinearModel::evaluate
OEMVector evaluate(const ConstVectorView &xi)
See ForwardModel class.
Definition: test_oem.cc:69
QuadraticModel::evaluate
OEMVector evaluate(const OEMVector &xi)
Virtual function of the FowardModel class.
Definition: test_oem.cc:139
source_dir
string source_dir
Definition: test_oem.cc:32
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
QuadraticModel::J
OEMMatrix J
Definition: test_oem.cc:158
run_plot_script
void run_plot_script(Engine *eng, string filename, string title)
Plot benchmark results.
Definition: test_oem.cc:347
add_noise
void add_noise(VectorView v, Numeric scale)
Add noise to vector.
Definition: test_utils.cc:38
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
LinearModel::y0
OEMVector y0
Definition: test_oem.cc:81
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
lin_alg.h
Linear algebra functions.
LinearModel::Jacobian
const OEMMatrix & Jacobian(const ConstVectorView &xi)
See ForwardModel class.
Definition: test_oem.cc:66
ARTS::Var::yf
Vector yf(Workspace &ws) noexcept
Definition: autoarts.h:7658
oem.h
Defines the ARTS interface to the invlib library.
ARTS::Var::time
Time time(Workspace &ws) noexcept
Definition: autoarts.h:7043
max
#define max(a, b)
Definition: legacy_continua.cc:20629
LinearModel
Linear Forward model.
Definition: test_oem.cc:51
test_oem_gauss_newton_sparse
void test_oem_gauss_newton_sparse(Index m, Index n, Index ntests)
Test sparse non-linear OEM.
Definition: test_oem.cc:897
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
QuadraticModel::QuadraticModel
QuadraticModel(Index m_, Index n_)
Constructor.
Definition: test_oem.cc:103
mult_general
void mult_general(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1181
ARTS::Var::nrows
Index nrows(Workspace &ws) noexcept
Definition: autoarts.h:4682
LinearModel::J
OEMMatrix J
Definition: test_oem.cc:80
test_oem_gauss_newton
void test_oem_gauss_newton(Engine *eng, Index m, Index n, Index ntests)
Test non-linear OEM.
Definition: test_oem.cc:686
generate_test_data
void generate_test_data(VectorView y, VectorView xa, MatrixView Se, MatrixView Sx)
Generate test data for linear OEM retrieval.
Definition: test_oem.cc:220
Joker
The Joker class.
Definition: matpackI.h:126
write_vector
void write_vector(ConstVectorView v, const char *filename)
Write vector to text file.
Definition: test_oem.cc:192
QuadraticModel::n
const unsigned int n
Definition: test_oem.cc:155
oem::Std
NormalizingSolver< Matrix, invlib::Standard > Std
The invlib standard solver.
Definition: oem.h:158
matpackII.h
Header file for sparse matrices.
benchmark_oem_linear
void benchmark_oem_linear(Engine *eng, Index n0, Index n1, Index ntests)
Benchmark linear oem.
Definition: test_oem.cc:499
get_maximum_error
Numeric get_maximum_error(ConstVectorView v1, ConstVectorView v2, bool relative)
Maximum element-wise error of two vectors.
Definition: test_utils.cc:297
run_test_matlab
Index run_test_matlab(Engine *eng, string filename)
Run test script in matlab.
Definition: test_oem.cc:258
atmlab_dir
string atmlab_dir
Definition: test_oem.cc:33
test_oem_linear
void test_oem_linear(Engine *eng, Index m, Index n, Index ntests)
Test linear_oem.
Definition: test_oem.cc:575
QuadraticModel::Hessians
OEMMatrix * Hessians
Definition: test_oem.cc:159
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
transpose
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
Definition: complex.cc:1509
ARTS::Var::ncols
Index ncols(Workspace &ws) noexcept
Definition: autoarts.h:4546
Vector
The Vector class.
Definition: matpackI.h:860
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
generate_linear_model
void generate_linear_model(MatrixView K)
Generate linear forward model.
Definition: test_oem.cc:245
test_oem_levenberg_marquardt
void test_oem_levenberg_marquardt(Engine *eng, Index m, Index n, Index ntests)
Test non-linear OEM.
Definition: test_oem.cc:790