63 :
m(J_.nrows()),
n(J_.ncols()),
J(J_),
y0(y0_) {}
109 sprintf(fname,
"J_t.txt");
112 for (
Index i = 0; i <
m; i++) {
116 sprintf(fname,
"H_%d_t.txt", (
int)i);
129 for (
Index i = 0; i <
m; i++) {
145 for (
Index i = 0; i <
n; i++) {
155 const unsigned int m,
n;
174 ofstream ofs(filename, ofstream::out);
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) <<
" ";
195 ofstream ofs(filename, ofstream::out);
197 for (
Index i = 0; i < n; i++) {
198 ofs << std::setprecision(20) <<
v[i] << endl;
263 string cmd =
"run('" + filename +
"');";
264 engEvalString(eng, cmd.c_str());
267 t = engGetVariable(eng,
"t");
289 mxArray *x_m, *G_m, *t;
292 string cmd =
"run('" + filename +
"');";
293 engEvalString(eng, cmd.c_str());
296 x_m = engGetVariable(eng,
"x");
297 G_m = engGetVariable(eng,
"G");
299 for (
Index i = 0; i < n; i++) {
300 x[i] = ((
Numeric*)mxGetData(x_m))[i];
302 for (
Index j = 0; j < m; j++) {
303 G(i, j) = ((
Numeric*)mxGetData(G_m))[j * n + i];
308 t = engGetVariable(eng,
"t");
325 int out = chdir(cmd.c_str());
329 string atmlab_init =
"run('" +
atmlab_dir +
"/atmlab/atmlab_init.m');";
333 engEvalString(eng, atmlab_init.c_str());
334 cmd =
"cd('" +
source_dir +
"/test_oem_files');";
335 engEvalString(eng, cmd.c_str());
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');");
361 int out = system(
"rm *_t.txt");
364 engEvalString(eng,
"close()");
382 Index step = (n1 - n0) / (ntests - 1);
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;
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;
393 for (
Index i = 0; i < ntests; i++) {
399 Index t, t1, t2, t_blas, t1_blas, t2_blas, t_m;
404 t = (t2 - t1) * 1000 / CLOCKS_PER_SEC;
409 t_blas = (t2_blas - t1_blas) * 1000 / CLOCKS_PER_SEC;
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;
420 cout << endl << endl;
442 Index step = (n1 - n0) / (ntests - 1);
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;
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;
453 for (
Index i = 0; i < ntests; i++) {
459 Index t, t1, t2, t_blas, t1_blas, t2_blas, t_m;
464 t = (t2 - t1) * 1000 / CLOCKS_PER_SEC;
469 t_blas = (t2_blas - t1_blas) * 1000 / CLOCKS_PER_SEC;
473 ofs << setw(5) << n << setw(10) << t_blas << setw(10) << t << setw(10);
475 cout << setw(5) << n << setw(10) << t_blas << setw(10) << t << setw(10);
480 cout << endl << endl;
500 Index step = (n1 - n0) / (ntests - 1);
503 ofstream ofs(
"times_linear.txt", ofstream::out);
505 ofs <<
"#" << setw(4) <<
"n" << setw(10) <<
"Matlab";
506 ofs << setw(10) <<
"C++" << endl;
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;
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),
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);
531 Index t, t1, t2, t_m;
542 t = (t2 - t1) * 1000 / CLOCKS_PER_SEC;
551 ofs << setw(5) << n << setw(10) << t << setw(10) << 42;
552 ofs << setw(10) << t_m << endl;
553 cout << setw(5) << n << setw(10) << t << setw(10) << t_m;
554 cout << endl << endl;
558 cout << endl << endl;
576 cout <<
"Testing linear OEM: m = " << m <<
", n = ";
577 cout << n <<
", ntests = " << ntests << endl << endl;
579 cout <<
"Test No. " << setw(15) <<
"Gauss-Newton";
580 cout << setw(20) <<
"Levenberg-Marquardt" << setw(15) <<
"Gain Matrix"
584 for (
Index i = 0; i < ntests; i++) {
604 LM_Pre_D lm_pre(Sa, Std_Pre(
std, pre));
605 GN_Pre gn_pre(1e-12, 100, Std_Pre(
std, pre));
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);
614 for (
Index j = 0; j < n; j++) {
615 x_norm[j] =
sqrt(Sa(j, j));
624 OEMVector x_std, x_nform, x_mform;
629 oem_std.compute(x_std, y, gn);
630 oem_nform.compute(x_nform, y, gn);
631 oem_mform.compute(x_mform, y, gn);
635 oem_std.compute(x_std_lm, y, lm);
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);
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);
652 Numeric err, err_G, err_lm, err_norm;
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;
687 cout <<
"Testing Gauss-Newton OEM: m = " << m <<
", n = ";
688 cout << n <<
", ntests = " << ntests << endl << endl;
690 cout <<
"Test No. " << setw(15) <<
"Standard" << setw(15) <<
"Normalized";
691 cout << setw(15) <<
"CG Solver" << setw(15) <<
"CG Normalized" << endl;
693 for (
Index i = 0; i < ntests; i++) {
696 OEMMatrix Se, Sa, SeInv, SaInv;
721 GN_Pre gn_pre(Std_Pre(
std, pre));
723 GN_CG_Pre gn_cg_pre(CG_Pre(cg, pre));
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);
731 for (
Index j = 0; j < n; j++) {
732 x_norm[j] =
sqrt(
abs(Sa(j, j)));
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);
756 cout << setw(9) << i + 1;
757 cout << setw(15) << e1 << setw(15) << e2 << setw(15) << e3;
758 cout << setw(15) << e4 << std::endl;
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);
770 cout << setw(9) << i + 1;
771 cout << setw(15) << e1 << setw(15) << e2 << setw(15) << e3;
772 cout << setw(15) << e4 << std::endl;
791 cout <<
"Testing Levenberg-Marquardt OEM: m = " << m <<
", n = ";
792 cout << n <<
", ntests = " << ntests << endl << endl;
794 cout <<
"Test No. " << setw(15) <<
"Standard" << setw(15) <<
"Normalized";
795 cout << setw(15) <<
"CG Solver" << setw(15) <<
"CG Normalized" << endl;
797 for (
Index i = 0; i < ntests; i++) {
800 OEMMatrix Se, Sa, SeInv, SaInv;
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));
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);
835 for (
Index j = 0; j < n; j++) {
836 x_norm[j] =
sqrt(
abs(Sa(j, j)));
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);
863 cout << setw(9) << i + 1;
864 cout << setw(15) << e1 << setw(15) << e2 << setw(15) << e3;
865 cout << setw(15) << e4 << std::endl;
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);
877 cout << setw(9) << i + 1;
878 cout << setw(15) << e1 << setw(15) << e2 << setw(15) << e3;
879 cout << setw(15) << e4 << std::endl;
898 cout <<
"Testing Gauss-Newton OEM: m = " << m <<
", n = ";
899 cout << n <<
", ntests = " << ntests << endl << endl;
901 cout <<
"Test No. " << setw(15) <<
"Standard" << setw(15) <<
"Normalized";
902 cout << setw(15) <<
"CG Solver" << setw(15) <<
"CG Normalized" << endl;
904 for (
Index i = 0; i < ntests; i++) {
918 Sparse Se_sparse(m, m), Sa_sparse(n, n);
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);
936 GN_Pre gn_pre(Std_Pre(
std, pre));
938 GN_CG_Pre gn_cg_pre(CG_Pre(cg, pre));
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);
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);
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);
965 cout << setw(9) << i + 1;
966 cout << setw(15) << e1 << setw(15) << e2 << setw(15) << e3;
967 cout << setw(15) << e4 << std::endl;
A constant view of a Matrix.
Index nrows() const noexcept
Index ncols() const noexcept
A constant view of a Vector.
LinearModel(Index m_, Index n_)
Default Constructor.
OEMVector evaluate(const ConstVectorView &xi)
See ForwardModel class.
LinearModel(ConstMatrixView J_, ConstVectorView y0_)
Construct a linear model from a given Jacobian J and offset vector y0.
const OEMMatrix & Jacobian(const ConstVectorView &xi)
See ForwardModel class.
~QuadraticModel()
Destructor.
OEMVector evaluate(const OEMVector &xi)
Virtual function of the FowardModel class.
QuadraticModel(Index m_, Index n_)
Constructor.
OEMMatrix Jacobian(const ConstVectorView &xi)
Virtual function of the FowardModel class.
void mult(MatrixView C, ConstMatrixView A, const Block &B)
void id_mat(MatrixView I)
Identity Matrix.
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Linear algebra functions.
Header file for sparse matrices.
void mult_general(VectorView y, const ConstMatrixView &M, const ConstVectorView &x) ARTS_NOEXCEPT
Matrix Vector multiplication.
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
VectorView std(VectorView std, const Vector &y, const ArrayOfVector &ys, const Index start, const Index end_tmp)
Compute the standard deviation of the ranged ys.
constexpr std::string_view Joker
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
invlib::GaussNewton< Numeric, CG > GN_CG
Gauss-Newton (GN) optimization using normed CG solver.
invlib::GaussNewton< Numeric, Std > GN
OEM Gauss-Newton optimization using normed ARTS QR solver.
NormalizingSolver< Matrix, invlib::Standard > Std
The invlib standard solver.
NormalizingSolver< Matrix, invlib::ConjugateGradient<> > CG
The invlib CG solver.
Defines the ARTS interface to the invlib library.
Numeric sqrt(const Rational r)
Square root.
void test_oem_levenberg_marquardt_sparse(Index m, Index n, Index ntests)
Test sparse non-linear OEM.
void setup_test_environment(Engine *&eng)
Setup the test environment.
void tidy_up_test_environment(Engine *eng)
Tidy up test environment.
void test_oem_gauss_newton(Engine *eng, Index m, Index n, Index ntests)
Test non-linear OEM.
void benchmark_inv(Engine *eng, Index n0, Index n1, Index ntests)
Matrix inversion benchmark.
Index run_oem_matlab(VectorView x, MatrixView G, Engine *eng, string filename)
Run test script in matlab and return result vector.
void generate_linear_model(MatrixView K)
Generate linear forward model.
void test_oem_levenberg_marquardt(Engine *eng, Index m, Index n, Index ntests)
Test non-linear OEM.
void test_oem_gauss_newton_sparse(Index m, Index n, Index ntests)
Test sparse non-linear OEM.
void test_oem_linear(Engine *eng, Index m, Index n, Index ntests)
Test linear_oem.
void run_plot_script(Engine *eng, string filename, string title)
Plot benchmark results.
Index run_test_matlab(Engine *eng, string filename)
Run test script in matlab.
void generate_test_data(VectorView y, VectorView xa, MatrixView Se, MatrixView Sx)
Generate test data for linear OEM retrieval.
void write_matrix(ConstMatrixView A, const char *fname)
Write matrix to text file.
void benchmark_oem_linear(Engine *eng, Index n0, Index n1, Index ntests)
Benchmark linear oem.
void write_vector(ConstVectorView v, const char *filename)
Write vector to text file.
void benchmark_mult(Engine *eng, Index n0, Index n1, Index ntests)
Matrix multiplication benchmark.
void random_fill_matrix(MatrixView A, Numeric range, bool positive)
Fill matrix with random values.
void random_fill_matrix_symmetric(MatrixView A, Numeric range, bool positive)
Generate random, symmetric matrix.
void add_noise(VectorView v, Numeric scale)
Add noise to vector.
Numeric get_maximum_error(ConstVectorView v1, ConstVectorView v2, bool relative)
Maximum element-wise error of two vectors.
void random_fill_vector(VectorView v, Numeric range, bool positive)
Fill vector with random values.
Utility functions for testing.