ARTS 2.5.9 (git: 825fa5f2)
test_covariance_matrix.cc
Go to the documentation of this file.
1/* Copyright (C) 2017
2 Simon Pfreundschuh <simon.pfreundschuh@chalmers.se>
3
4 This program is free software; you can redistribute it and/or modify it
5 under the terms of the GNU General Public License as published by the
6 Free Software Foundation; either version 2, or (at your option) any
7 later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 USA. */
18
27#include <cstdlib>
28#include <limits>
29#include <random>
30#include <utility>
31
32#include "array.h"
33#include "covariance_matrix.h"
34#include "jacobian.h"
35#include "lin_alg.h"
36#include "test_utils.h"
37#include "xml_io.h"
38
39// Type alias describing a retrieval setup for testing.
40using RetrievalData = std::tuple<ArrayOfRetrievalQuantity, ArrayOfArrayOfIndex>;
41
54template <typename F>
55std::shared_ptr<Matrix> create_covariance_matrix_1D(
56 Index i, Index j, ArrayOfRetrievalQuantity jqs, F f) {
57 if (i > j) {
58 std::swap(i, j);
59 }
60
61 const RetrievalQuantity& q1 = jqs[i];
62 const RetrievalQuantity& q2 = jqs[j];
63
64 const Vector& gv1 = q1.Grids()[0];
65 const Vector& gv2 = q2.Grids()[0];
66
67 std::shared_ptr<Matrix> c =
68 std::make_shared<Matrix>(gv1.nelem(), gv2.nelem());
69
70 for (Index k = 0; k < gv1.nelem(); k++) {
71 for (Index l = 0; l < gv2.nelem(); l++) {
72 (*c)(k, l) = f(gv1[k], gv2[l]);
73 }
74 }
75 return c;
76}
77
81template <typename F>
82std::shared_ptr<Sparse> create_sparse_covariance_matrix_1D(
83 Index i, Index j, ArrayOfRetrievalQuantity jqs, F f) {
84 if (i > j) {
85 std::swap(i, j);
86 }
87
88 const RetrievalQuantity& q1 = jqs[i];
89 const RetrievalQuantity& q2 = jqs[j];
90
91 const Vector& gv1 = q1.Grids()[0];
92 const Vector& gv2 = q2.Grids()[0];
93
94 std::shared_ptr<Sparse> s =
95 std::make_shared<Sparse>(gv1.nelem(), gv2.nelem());
96
97 Index nelem = gv1.nelem() * gv2.nelem();
98 ArrayOfIndex row_indices{}, col_indices{};
99 ArrayOfNumeric elements{};
100 row_indices.reserve(nelem);
101 col_indices.reserve(nelem);
102 elements.reserve(nelem);
103 for (Index k = 0; k < gv1.nelem(); k++) {
104 for (Index l = 0; l < gv2.nelem(); l++) {
105 row_indices.push_back(k);
106 col_indices.push_back(l);
107 elements.push_back(f(gv1[k], gv2[l]));
108 }
109 }
110 s->insert_elements(nelem, row_indices, col_indices, Vector(elements));
111 return s;
112}
113
121 std::random_device rd;
122 std::mt19937 gen(rd());
123 std::uniform_int_distribution<> n_rqs_dist(2, 10), n_gs_dist(10, 100);
124
125 Index n_rqs = n_rqs_dist(gen);
126 Numeric g_begin = 0.0;
127 Numeric g_end = 100.0;
129
130 for (Index i = 0; i < n_rqs; i++) {
131 Index n_gs = n_gs_dist(gen);
132 ArrayOfVector gvs(1);
133 gvs[0] = Vector(n_gs);
134 for (Index j = 0; j < n_gs; j++) {
135 Numeric frac = static_cast<Numeric>(j) / static_cast<Numeric>(n_gs);
136 gvs[0][j] = g_begin + (g_end - g_begin) * frac;
137 }
138 rqs.emplace_back(JacobianTarget(), "subtag", "subsubtag", "mode", 0.0, gvs);
139 }
140
141 ArrayOfArrayOfIndex jis(n_rqs);
142 Index ji = 0;
143 for (Index i = 0; i < n_rqs; i++) {
144 Index n_grid_points = 1;
145 const ArrayOfVector& gs = rqs[i].Grids();
146 for (auto& g : gs) {
147 n_grid_points *= g.nelem();
148 }
149 jis[i] = ArrayOfIndex{ji, ji + n_grid_points - 1};
150 ji += n_grid_points;
151 }
152
153 return std::make_tuple(rqs, jis);
154}
155
164 const ArrayOfArrayOfIndex& jis) {
165 std::random_device rd;
166 std::mt19937 gen(rd());
167 std::uniform_real_distribution<> dis(0.0, 1.0);
168 CovarianceMatrix covmat{};
169
170 auto g = [](Numeric x, Numeric y) -> Numeric {
171 if (x == y)
172 return 1.0;
173 else
174 return 0.0;
175 };
176 auto g2 = [](Numeric x, Numeric y) -> Numeric {
177 if (x == y)
178 return 0.1;
179 else
180 return 0.0;
181 };
182
183 for (size_t i = 0; i < rqs.size(); i++) {
184 Range range(jis[i][0], jis[i][1] - jis[i][0] + 1);
185 auto inds = std::make_pair(i, i);
186
187 Numeric p = dis(gen);
188 if (p < 0.0) {
189 std::shared_ptr<Matrix> m = create_covariance_matrix_1D(i, i, rqs, g);
190 covmat.add_correlation(Block(range, range, inds, m));
191 } else {
192 std::shared_ptr<Sparse> m =
194 covmat.add_correlation(Block(range, range, inds, m));
195 }
196 }
197
198 for (size_t i = 0; i < rqs.size(); i++) {
199 for (size_t j = i + 1; j < rqs.size(); j++) {
200 Range row_range(jis[i][0], jis[i][1] - jis[i][0] + 1);
201 Range col_range(jis[j][0], jis[j][1] - jis[j][0] + 1);
202 auto inds = std::make_pair(i, j);
203 Numeric p = dis(gen);
204 if (p > 0.8) {
205 p = dis(gen);
206 if (p < 0.0) {
207 std::shared_ptr<Matrix> m = create_covariance_matrix_1D(i, j, rqs, g);
208 covmat.add_correlation(Block(row_range, col_range, inds, m));
209 } else {
210 std::shared_ptr<Sparse> m =
212 covmat.add_correlation(Block(row_range, col_range, inds, m));
213 }
214 }
215 }
216 }
217 return covmat;
218}
219
228 Numeric e = 0.0;
229 for (Index i = 0; i < n_tests; i++) {
232 std::tie(rqs, jis) = setup_retrieval_1D();
234
235 Index n = covmat.ncols();
236 Matrix A(covmat);
237 Vector v(n), w(n), w_ref(n);
238 random_fill_vector(v, 10.0, false);
239
240 mult(w, covmat, v);
241 mult(w_ref, A, v);
242
243 e = std::max(e, get_maximum_error(w, w_ref, true));
244 }
245 return e;
246}
247
256 Numeric e = 0.0;
257 for (Index i = 0; i < n_tests; i++) {
260 std::tie(rqs, jis) = setup_retrieval_1D();
262
263 Index n = covmat.ncols();
264 Matrix A(covmat), B(n, n), C(n, n), C_ref(n, n);
265 random_fill_matrix(B, 10.0, false);
266
267 mult(C, covmat, B);
268 mult(C_ref, A, B);
269 e = std::max(e, get_maximum_error(C, C_ref, true));
270
271 mult(C, B, covmat);
272 mult(C_ref, B, A);
273 e = std::max(e, get_maximum_error(C, C_ref, true));
274 }
275 return e;
276}
277
289 Numeric e = 0.0;
290 for (Index i = 0; i < n_tests; i++) {
293 std::tie(rqs, jis) = setup_retrieval_1D();
295
296 Index n = covmat.ncols();
297 Matrix A(covmat), B(n, n), B_ref(n, n), C(n, n);
298 random_fill_matrix(B, 10.0, false);
299
300 covmat.compute_inverse();
301 mult_inv(B, covmat, A);
302 id_mat(B_ref);
303 mult(C, covmat, B_ref);
304 e = std::max(e, get_maximum_error(B, B_ref, true));
305 id_mat(B_ref);
306 mult_inv(B, A, covmat);
307 e = std::max(e, get_maximum_error(B, B_ref, true));
308 }
309 return e;
310}
311
320 Numeric e = 0.0;
321 for (Index i = 0; i < n_tests; i++) {
324 std::tie(rqs, jis) = setup_retrieval_1D();
326
327 Index n = covmat.ncols();
328 Matrix A(covmat), B(n, n), B_ref(n, n), C(n, n);
329 B = 0.0;
330 B_ref = B;
331
332 covmat.compute_inverse();
333
334 B += covmat;
335 B_ref += A;
336 e = std::max(e, get_maximum_error(B, B_ref, true));
337
338 add_inv(B, covmat);
339 inv(C, A);
340 B_ref += C;
341 e = std::max(e, get_maximum_error(B, B_ref, true));
342 }
343 return e;
344}
345
355 Numeric e = 0.0;
356 for (Index i = 0; i < n_tests; i++) {
359 std::tie(rqs, jis) = setup_retrieval_1D();
361
362 Vector diag_1 = covmat.diagonal();
363 Vector diag_2 = Matrix(covmat).diagonal();
364 e = std::max(e, get_maximum_error(diag_1, diag_2, true));
365 }
366 return e;
367}
368
377 Numeric e(0.0);
378 for (Index i = 0; i < n_tests; i++) {
381 std::tie(rqs, jis) = setup_retrieval_1D();
382 CovarianceMatrix covmat_1(random_covariance_matrix(rqs, jis)), covmat_2{};
383 covmat_1.compute_inverse();
384
385 xml_write_to_file("test.xml", covmat_1, FILE_TYPE_ASCII, 0, Verbosity());
386 xml_read_from_file("test.xml", covmat_2, Verbosity());
387
388 Matrix A(covmat_1), B(covmat_2), C(covmat_1);
389 e = std::max(e, get_maximum_error(A, B, true));
390
391 id_mat(C);
392 mult_inv(A, covmat_1, C);
393 mult_inv(B, covmat_2, C);
394 e = std::max(e, get_maximum_error(A, B, true));
395 }
396 return 0.0;
397}
398
399template <typename MatrixType>
400void covmat_seSet(CovarianceMatrix& covmat,
401 const MatrixType& block,
402 const Verbosity& /*v*/);
403
404template <typename MatrixType>
405void covmat_sxSet(CovarianceMatrix& covmat,
406 const MatrixType& block,
407 const Verbosity& /*v*/);
408
409template <typename MatrixType>
410void covmat_seAddBlock(CovarianceMatrix& covmat_se,
411 const MatrixType& block,
412 const Index& i,
413 const Index& j,
414 const Verbosity&);
415
416template <typename MatrixType>
418 const MatrixType& block,
419 const Index& i,
420 const Index& j,
421 const Verbosity&);
422
423template <typename MatrixType>
424void covmat_sxAddBlock(CovarianceMatrix& covmat_se,
425 const ArrayOfRetrievalQuantity& jq,
426 const MatrixType& block,
427 const Index& i,
428 const Index& j,
429 const Verbosity&);
430
431template <typename MatrixType>
433 const ArrayOfRetrievalQuantity& jq,
434 const MatrixType& block,
435 const Index& i,
436 const Index& j,
437 const Verbosity&);
438
443 CovarianceMatrix covmat{};
444 Matrix A(10, 10);
445 Sparse A_sparse(10, 10);
446 Matrix B(20, 20);
447 Matrix C(10, 20);
448
449 // covmat_seSet
450
451 covmat_seSet(covmat, A, Verbosity());
452 ARTS_ASSERT(covmat.ncols() == 10);
453
454 // covmat_seAddBlock
455
456 covmat_seAddBlock(covmat, A_sparse, -1, -1, Verbosity());
457 ARTS_ASSERT(covmat.ncols() == 20);
458
459 try {
460 covmat_seAddBlock(covmat, A, 3, 3, Verbosity());
461 // This should fail.
462 ARTS_ASSERT(false);
463 } catch (const std::runtime_error&) {
464 }
465
466 covmat_seAddBlock(covmat, A, 0, 1, Verbosity());
467
468 covmat_seAddBlock(covmat, A, 2, 2, Verbosity());
469
470 try {
471 covmat_seAddBlock(covmat, B, 1, 2, Verbosity());
472 // This should fail.
473 ARTS_ASSERT(false);
474 } catch (const std::runtime_error&) {
475 }
476
477 covmat_seAddInverseBlock(covmat, A, 0, 1, Verbosity());
478
479 try {
480 covmat_seAddInverseBlock(covmat, B, 3, 3, Verbosity());
481 // This should fail.
482 ARTS_ASSERT(false);
483 } catch (const std::runtime_error&) {
484 }
485
486 // covmat_sxSet
487 covmat_sxSet(covmat, A, Verbosity());
488 ARTS_ASSERT(covmat.ncols() == 10);
489
490 // covmat_sxAddBlock
491
492 A.resize(1000, 1000);
493 A_sparse.resize(8000, 8000);
494 C.resize(1000, 8000);
495
496 ArrayOfVector grids_1{Vector(10), Vector(10), Vector(10)};
497 ArrayOfVector grids_2{Vector(20), Vector(20), Vector(20)};
498
499 RetrievalQuantity rq_1(JacobianTarget(), "st", "sst", "m", 0.1, grids_1);
500 RetrievalQuantity rq_2(JacobianTarget(), "st", "sst", "m", 0.1, grids_2);
501
503 rqs.push_back(rq_1);
504 rqs.push_back(rq_2);
505
506 covmat_sxAddBlock(covmat, rqs, A_sparse, -1, -1, Verbosity());
507
508 try {
509 covmat_sxAddBlock(covmat, rqs, A_sparse, 0, 1, Verbosity());
510 // This should fail.
511 ARTS_ASSERT(false);
512 } catch (const std::runtime_error&) {
513 }
514
515 covmat_sxAddBlock(covmat, rqs, C, 0, 1, Verbosity());
516 covmat_sxAddInverseBlock(covmat, rqs, A, 0, 0, Verbosity());
517
518 try {
519 covmat_sxAddInverseBlock(covmat, rqs, C, 1, 1, Verbosity());
520 // This should fail.
521 ARTS_ASSERT(false);
522 } catch (const std::runtime_error&) {
523 }
524}
525
533namespace {
534#pragma GCC diagnostic push
535#pragma GCC diagnostic ignored "-Wunused-function"
536#pragma GCC diagnostic ignored "-Wconversion"
537
538#include "invlib/algebra.h"
539#include "invlib/interfaces/arts_wrapper.h"
540
541Numeric test_invlib_wrapper(Index n_tests) {
542 Numeric e = 0.0;
543 for (Index i = 0; i < n_tests; i++) {
546 std::tie(rqs, jis) = setup_retrieval_1D();
548 invlib::Matrix<ArtsCovarianceMatrixWrapper> wrapper(covmat);
549 Index n = covmat.ncols();
550
551 // Multiplication by Vector
552 invlib::Vector<ArtsVector> v{}, w{}, w_ref{};
553 v.resize(n);
554 w.resize(n);
555 w_ref.resize(n);
556 random_fill_vector(v, 10.0, false);
557
558 mult(w, covmat, v);
559 w_ref = wrapper * v;
560 e = std::max(e, get_maximum_error(w, w_ref, true));
561
562 solve(w, covmat, v);
563 w_ref = inv(wrapper) * v;
564 e = std::max(e, get_maximum_error(w, w_ref, true));
565 w_ref = inv(inv(wrapper)) * v;
566
567 // Multiplication by Matrix
568 invlib::Matrix<ArtsMatrix> A{static_cast<Matrix>(covmat)}, B{}, C{},
569 C_ref{};
570 B.resize(n, n);
571 C.resize(n, n);
572 C_ref.resize(n, n);
573
574 random_fill_matrix(B, 10.0, false);
575 mult(C, covmat, B);
576 C_ref = wrapper * B;
577 e = std::max(e, get_maximum_error(C, C_ref, true));
578
579 // Inverse
580 mult_inv(C, covmat, B);
581 C_ref = invlib::inv(wrapper) * B;
582 e = std::max(e, get_maximum_error(C, C_ref, true));
583
584 // Addition to Matrix
585 C += covmat;
586 C_ref = C_ref + wrapper;
587 e = std::max(e, get_maximum_error(C, C_ref, true));
588
589 // Addition of Inverse to Matrix
590 add_inv(C, covmat);
591 C_ref = C_ref + inv(wrapper);
592 e = std::max(e, get_maximum_error(C, C_ref, true));
593 }
594 return e;
595}
596
597#pragma GCC diagnostic pop
598} // namespace
599
600int main() {
601 Numeric e = 0.0, e_max = 0.0;
602
603 std::cout << "Testing covariance matrix: " << std::endl;
604 e = test_addition(10);
605 std::cout << "\tAddition: " << e << std::endl;
606 e_max = std::max(e, e_max);
607 if (e_max > 1e-5) {
608 return -1;
609 }
610
612 std::cout << "\tVector Multiplication: " << e << std::endl;
613 e_max = std::max(e, e_max);
614 if (e_max > 1e-5) {
615 return -1;
616 }
617
619 std::cout << "\tMatrix Multiplication: " << e << std::endl;
620 e_max = std::max(e, e_max);
621 if (e_max > 1e-5) {
622 return -1;
623 }
624
625 e = test_inverse(10);
626 std::cout << "\tInverse: " << e << std::endl;
627 e_max = std::max(e, e_max);
628 if (e_max > 1e-5) {
629 return -1;
630 }
631
632 e = test_io(10);
633 std::cout << "\tXML IO: " << e << std::endl;
634 e_max = std::max(e, e_max);
635 if (e_max > 1e-5) {
636 return -1;
637 }
638
639 e = test_invlib_wrapper(10);
640 std::cout << "\tinvlib Wrapper: " << e << std::endl;
641 e_max = std::max(e, e_max);
642 if (e_max > 1e-5) {
643 return -1;
644 }
645
646 e = test_diagonal(10);
647 std::cout << "\tdiagonal " << e << std::endl;
648 e_max = std::max(e, e_max);
649 if (e_max > 1e-5) {
650 return -1;
651 }
652
653 std::cout << std::endl << "\tTesting workspace functions ... ";
655 std::cout << " DONE." << std::endl;
656
657 return 0;
658}
This file contains the definition of Array.
This can be used to make arrays out of anything.
Definition: array.h:48
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
ConstVectorView diagonal() const ARTS_NOEXCEPT
Matrix diagonal as vector.
Definition: matpackI.cc:469
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
void compute_inverse() const
Compute the inverse of this correlation matrix.
Vector diagonal() const
Diagonal elements as vector.
The Matrix class.
Definition: matpackI.h:1285
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
The range class.
Definition: matpackI.h:160
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition: jacobian.h:325
const ArrayOfVector & Grids() const
Returns the grids of the retrieval.
Definition: jacobian.h:408
The Vector class.
Definition: matpackI.h:910
void mult(MatrixView C, ConstMatrixView A, const Block &B)
void add_inv(MatrixView A, const CovarianceMatrix &B)
void mult_inv(MatrixView C, ConstMatrixView A, const CovarianceMatrix &B)
void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v)
Header files of CovarianceMatrix class.
#define ARTS_ASSERT(condition,...)
Definition: debug.h:82
Routines for setting up the jacobian.
Jacobian::Target JacobianTarget
Definition: jacobian.h:320
Array< RetrievalQuantity > ArrayOfRetrievalQuantity
Definition: jacobian.h:544
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.
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
constexpr Numeric e
Elementary charge convenience name [C].
The Sparse class.
Definition: matpackII.h:75
void resize(Index r, Index c)
Resize function.
Definition: matpackII.cc:361
std::tuple< ArrayOfRetrievalQuantity, ArrayOfArrayOfIndex > RetrievalData
Numeric test_multiplication_by_vector(Index n_tests)
Tests the multiplication of covariance matrices with vectors.
void covmat_sxAddBlock(CovarianceMatrix &covmat_se, const ArrayOfRetrievalQuantity &jq, const MatrixType &block, const Index &i, const Index &j, const Verbosity &)
Definition: m_retrieval.cc:577
Numeric test_diagonal(Index n_tests)
Test extraction of diagonal.
void covmat_sxSet(CovarianceMatrix &covmat, const MatrixType &block, const Verbosity &)
Definition: m_retrieval.cc:558
void covmat_seSet(CovarianceMatrix &covmat, const MatrixType &block, const Verbosity &)
Definition: m_retrieval.cc:539
Numeric test_multiplication_by_matrix(Index n_tests)
Tests the multiplication of covariance matrices with matrices.
void test_workspace_methods()
Test general functionality of CovarianceMatrix class.
std::shared_ptr< Sparse > create_sparse_covariance_matrix_1D(Index i, Index j, ArrayOfRetrievalQuantity jqs, F f)
Same as create_covariance_matrix_1D but creates a matrix of type sparse.
Numeric test_addition(Index n_tests)
Test addition of covariance matrices and inverse covariance matrices.
RetrievalData setup_retrieval_1D()
Setup a random retrieval case for testing.
Numeric test_inverse(Index n_tests)
Tests the inversion of covariance matrices by computing the products of inverses of covariance matric...
void covmat_seAddBlock(CovarianceMatrix &covmat_se, const MatrixType &block, const Index &i, const Index &j, const Verbosity &)
Definition: m_retrieval.cc:390
CovarianceMatrix random_covariance_matrix(const ArrayOfRetrievalQuantity &rqs, const ArrayOfArrayOfIndex &jis)
Create a covariance matrix with correlations between retrieval quantitites given in rqs and jis.
std::shared_ptr< Matrix > create_covariance_matrix_1D(Index i, Index j, ArrayOfRetrievalQuantity jqs, F f)
Create a block of a covariance matrix representating a correlation between the retrieval quantities i...
void covmat_sxAddInverseBlock(CovarianceMatrix &covmat_se, const ArrayOfRetrievalQuantity &jq, const MatrixType &block, const Index &i, const Index &j, const Verbosity &)
Definition: m_retrieval.cc:649
Numeric test_io(Index n_tests)
Test input and output of covariance matrices.
int main()
void covmat_seAddInverseBlock(CovarianceMatrix &covmat_se, const MatrixType &block, const Index &i, const Index &j, const Verbosity &)
Definition: m_retrieval.cc:470
void random_fill_matrix(MatrixView A, Numeric range, bool positive)
Fill matrix with random values.
Definition: test_utils.cc:59
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
#define w
#define c
This file contains basic functions to handle XML data files.
void xml_write_to_file(const String &filename, const T &type, const FileType ftype, const Index no_clobber, const Verbosity &verbosity)
Write data to XML file.
Definition: xml_io.h:172
void xml_read_from_file(const String &filename, T &type, const Verbosity &verbosity)
Reads data from XML file.
Definition: xml_io.h:151
@ FILE_TYPE_ASCII
Definition: xml_io_base.h:43