ARTS 2.5.9 (git: 825fa5f2)
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
24using std::abs;
25using std::cout;
26using std::endl;
27using std::min;
28using std::ofstream;
29using std::setw;
30using std::string;
31
32string source_dir = SOURCEDIR;
33string atmlab_dir = ATMLABDIR;
34
35// Forward declarations.
36
37void write_matrix(ConstMatrixView A, const char* fname);
38
40
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
170void 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
192void 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
258Index 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
285Index 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
321void 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
347void 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
360void tidy_up_test_environment(Engine* eng) {
361 int out = system("rm *_t.txt");
362 (void)out;
363
364 engEvalString(eng, "close()");
365}
366
368
381void 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
441void 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
499void 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
575void 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
686void 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
790void 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
987int 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}
A constant view of a Matrix.
Definition: matpackI.h:1065
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
A constant view of a Vector.
Definition: matpackI.h:521
The Joker class.
Definition: matpackI.h:126
Linear Forward model.
Definition: test_oem.cc:51
const Index n
Definition: test_oem.cc:77
OEMVector y0
Definition: test_oem.cc:81
LinearModel(Index m_, Index n_)
Default Constructor.
Definition: test_oem.cc:54
OEMVector evaluate(const ConstVectorView &xi)
See ForwardModel class.
Definition: test_oem.cc:69
LinearModel(ConstMatrixView J_, ConstVectorView y0_)
Construct a linear model from a given Jacobian J and offset vector y0.
Definition: test_oem.cc:62
OEMMatrix J
Definition: test_oem.cc:80
const Index m
Definition: test_oem.cc:77
const OEMMatrix & Jacobian(const ConstVectorView &xi)
See ForwardModel class.
Definition: test_oem.cc:66
The MatrixView class.
Definition: matpackI.h:1188
The Matrix class.
Definition: matpackI.h:1285
Quadratic Forward Model.
Definition: test_oem.cc:91
~QuadraticModel()
Destructor.
Definition: test_oem.cc:122
OEMMatrix * Hessians
Definition: test_oem.cc:159
const unsigned int n
Definition: test_oem.cc:155
OEMMatrix J
Definition: test_oem.cc:158
OEMVector evaluate(const OEMVector &xi)
Virtual function of the FowardModel class.
Definition: test_oem.cc:139
QuadraticModel(Index m_, Index n_)
Constructor.
Definition: test_oem.cc:103
OEMMatrix Jacobian(const ConstVectorView &xi)
Virtual function of the FowardModel class.
Definition: test_oem.cc:125
const unsigned int m
Definition: test_oem.cc:155
The VectorView class.
Definition: matpackI.h:674
The Vector class.
Definition: matpackI.h:910
void mult(MatrixView C, ConstMatrixView A, const Block &B)
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:488
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:168
Linear algebra functions.
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
Header file for sparse matrices.
void mult_general(VectorView y, const ConstMatrixView &M, const ConstVectorView &x) ARTS_NOEXCEPT
Matrix Vector multiplication.
Definition: matpackI.cc:1138
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
Defines the ARTS interface to the invlib library.
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:705
The Sparse class.
Definition: matpackII.h:75
void test_oem_levenberg_marquardt_sparse(Index m, Index n, Index ntests)
Test sparse non-linear OEM.
Definition: test_oem.cc:985
void setup_test_environment(Engine *&eng)
Setup the test environment.
Definition: test_oem.cc:321
void tidy_up_test_environment(Engine *eng)
Tidy up test environment.
Definition: test_oem.cc:360
void test_oem_gauss_newton(Engine *eng, Index m, Index n, Index ntests)
Test non-linear OEM.
Definition: test_oem.cc:686
void benchmark_inv(Engine *eng, Index n0, Index n1, Index ntests)
Matrix inversion benchmark.
Definition: test_oem.cc:381
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
string atmlab_dir
Definition: test_oem.cc:33
void generate_linear_model(MatrixView K)
Generate linear forward model.
Definition: test_oem.cc:245
void test_oem_levenberg_marquardt(Engine *eng, Index m, Index n, Index ntests)
Test non-linear OEM.
Definition: test_oem.cc:790
void test_oem_gauss_newton_sparse(Index m, Index n, Index ntests)
Test sparse non-linear OEM.
Definition: test_oem.cc:897
void test_oem_linear(Engine *eng, Index m, Index n, Index ntests)
Test linear_oem.
Definition: test_oem.cc:575
void run_plot_script(Engine *eng, string filename, string title)
Plot benchmark results.
Definition: test_oem.cc:347
string source_dir
Definition: test_oem.cc:32
Index run_test_matlab(Engine *eng, string filename)
Run test script in matlab.
Definition: test_oem.cc:258
void generate_test_data(VectorView y, VectorView xa, MatrixView Se, MatrixView Sx)
Generate test data for linear OEM retrieval.
Definition: test_oem.cc:220
void write_matrix(ConstMatrixView A, const char *fname)
Write matrix to text file.
Definition: test_oem.cc:170
void benchmark_oem_linear(Engine *eng, Index n0, Index n1, Index ntests)
Benchmark linear oem.
Definition: test_oem.cc:499
int main()
Definition: test_oem.cc:987
void write_vector(ConstVectorView v, const char *filename)
Write vector to text file.
Definition: test_oem.cc:192
void benchmark_mult(Engine *eng, Index n0, Index n1, Index ntests)
Matrix multiplication benchmark.
Definition: test_oem.cc:441
void random_fill_matrix(MatrixView A, Numeric range, bool positive)
Fill matrix with random values.
Definition: test_utils.cc:59
void random_fill_matrix_symmetric(MatrixView A, Numeric range, bool positive)
Generate random, symmetric matrix.
Definition: test_utils.cc:156
void add_noise(VectorView v, Numeric scale)
Add noise to vector.
Definition: test_utils.cc:38
Numeric get_maximum_error(ConstVectorView v1, ConstVectorView v2, bool relative)
Maximum element-wise error of two vectors.
Definition: test_utils.cc:297
void random_fill_vector(VectorView v, Numeric range, bool positive)
Fill vector with random values.
Definition: test_utils.cc:230
Utility functions for testing.
#define v