ARTS 2.5.11 (git: 725533f0)
covariance_matrix.cc
Go to the documentation of this file.
1
9#include <queue>
10#include <tuple>
11#include <utility>
12#include <vector>
13
14#include "covariance_matrix.h"
15#include "lapack.h"
16#include "lin_alg.h"
17#include "matpack_math.h"
18
19//------------------------------------------------------------------------------
20// Correlations
21//------------------------------------------------------------------------------
22void mult(MatrixView C, ConstMatrixView A, const Block &B) {
23 ARTS_ASSERT((B.sparse_ != nullptr) || (B.dense_ != nullptr));
24
25 MatrixView CView(C(joker, B.get_column_range()));
26 MatrixView CTView(C(joker, B.get_row_range()));
27 ConstMatrixView AView(A(joker, B.get_row_range()));
28 ConstMatrixView ATView(A(joker, B.get_column_range()));
29
30 Index i, j;
31 std::tie(i, j) = B.get_indices();
32
34 mult(CView, AView, *B.dense_);
35 } else {
36 mult(CView, AView, *B.sparse_);
37 }
38
39 // Only upper blocks are stored, so if the correlation is between different RQs
40 // we also need to account for the implicit lower block.
41 if (i != j) {
43 mult(CTView, ATView, transpose(*B.dense_));
44 } else {
45 Matrix D(CTView.ncols(), CTView.nrows());
46 mult(D, *B.sparse_, transpose(ATView));
47 CTView += transpose(D);
48 }
49 }
50}
51
52void mult(MatrixView C, const Block &A, ConstMatrixView B) {
53 ARTS_ASSERT((A.sparse_ != nullptr) || (A.dense_ != nullptr));
54
55 MatrixView CView(C(A.get_row_range(), joker));
56 MatrixView CTView(C(A.get_column_range(), joker));
57 ConstMatrixView BView(B(A.get_column_range(), joker));
58 ConstMatrixView BTView(B(A.get_row_range(), joker));
59
61 mult(CView, *A.dense_, BView);
62 } else {
63 mult(CView, *A.sparse_, BView);
64 }
65
66 Index i, j;
67 std::tie(i, j) = A.get_indices();
68
69 // Only upper blocks are stored, so if the correlation is between different RQs
70 // we also need to account for the implicit lower block.
71 if (i != j) {
73 mult(CTView, transpose(*A.dense_), BTView);
74 } else {
75 Matrix D(CTView.ncols(), CTView.nrows());
76 mult(D, transpose(BTView), *A.sparse_);
77 CTView += transpose(D);
78 }
79 }
80}
81
82void mult(VectorView w, const Block &A, ConstVectorView v) {
83 VectorView wview(w[A.get_row_range()]), wtview(w[A.get_column_range()]);
84 ConstVectorView vview(v[A.get_column_range()]), vtview(v[A.get_row_range()]);
85
87 mult(wview, *A.dense_, vview);
88 } else {
89 mult(wview, *A.sparse_, vview);
90 }
91
92 Index i, j;
93 std::tie(i, j) = A.get_indices();
94
95 if (i != j) {
97 mult(wtview, transpose(*A.dense_), vtview);
98 } else {
99 transpose_mult(wtview, *A.sparse_, vtview);
100 }
101 }
102}
103
104MatrixView operator+=(MatrixView A, const Block &B) {
105 MatrixView Aview(A(B.get_row_range(), B.get_column_range()));
106 MatrixView ATview(A(B.get_column_range(), B.get_row_range()));
108 Aview += B.get_dense();
109 } else {
110 Aview += static_cast<const Matrix>(B.get_sparse());
111 }
112
113 Index i, j;
114 std::tie(i, j) = B.get_indices();
115
116 if (i != j) {
118 ATview += transpose(B.get_dense());
119 } else {
120 ATview += transpose(static_cast<const Matrix>(B.get_sparse()));
121 }
122 }
123 return A;
124}
125
126//------------------------------------------------------------------------------
127// Covariance Matrix
128//------------------------------------------------------------------------------
129CovarianceMatrix::operator Matrix() const {
130 Index n = nrows();
131 Matrix A(n, n);
132 A = 0.0;
133
134 for (const Block &c : correlations_) {
135 MatrixView Aview = A(c.get_row_range(), c.get_column_range());
136 if (c.get_matrix_type() == Block::MatrixType::dense) {
137 Aview = c.get_dense();
138 } else {
139 Aview = static_cast<const Matrix>(c.get_sparse());
140 }
141
142 Index ci, cj;
143 std::tie(ci, cj) = c.get_indices();
144 if (ci != cj) {
145 MatrixView ATview = A(c.get_column_range(), c.get_row_range());
146 if (c.get_matrix_type() == Block::MatrixType::dense) {
147 ATview = transpose(c.get_dense());
148 } else {
149 ATview = transpose(static_cast<const Matrix>(c.get_sparse()));
150 }
151 }
152 }
153 return A;
154}
155
157 Index n = nrows();
158 Matrix A(n, n);
159 A = 0.0;
160
161 for (const Block &c : inverses_) {
162 MatrixView Aview = A(c.get_row_range(), c.get_column_range());
163 if (c.get_matrix_type() == Block::MatrixType::dense) {
164 Aview = c.get_dense();
165 } else {
166 Aview = static_cast<const Matrix>(c.get_sparse());
167 }
169 Index ci, cj;
170 std::tie(ci, cj) = c.get_indices();
171 if (ci != cj) {
172 MatrixView ATview = A(c.get_column_range(), c.get_row_range());
173 if (c.get_matrix_type() == Block::MatrixType::dense) {
174 ATview = transpose(c.get_dense());
175 } else {
176 ATview = transpose(static_cast<const Matrix>(c.get_sparse()));
177 }
178 }
179 }
180 return A;
181}
182
184 Index m1 = 0;
185
186 for (const Block &c : correlations_) {
187 Index i, j;
188 std::tie(i, j) = c.get_indices();
189 if (i == j) {
190 m1 += c.nrows();
191 }
192 }
193
194 Index m2 = 0;
195 for (const Block &c : inverses_) {
196 Index i, j;
197 std::tie(i, j) = c.get_indices();
198 if (i == j) {
199 m2 += c.nrows();
200 }
201 }
202
203 return std::max(m1, m2);
204}
205
206Index CovarianceMatrix::ncols() const { return nrows(); }
207
209 Index m = 0;
210
211 for (const Block &c : correlations_) {
212 Index i, j;
213 std::tie(i, j) = c.get_indices();
214 if (i == j) {
215 ++m;
216 }
217 }
218 return m;
219}
220
222 Index m = 0;
223
224 for (const Block &c : inverses_) {
225 Index i, j;
226 std::tie(i, j) = c.get_indices();
227 if (i == j) {
228 ++m;
229 }
230 }
231 return m;
232}
233
234Index CovarianceMatrix::nblocks() const { return correlations_.size(); }
235
236bool CovarianceMatrix::has_block(Index i, Index j) {
237 if (i > j) {
238 std::swap(i, j);
239 }
240
241 bool result = false;
242 for (const Block &b : correlations_) {
243 result |= b.get_indices() == std::make_pair(i, j);
244 }
245
246 return result;
247}
248
249const Block *CovarianceMatrix::get_block(Index i, Index j) {
250 if (i > j) {
251 std::swap(i, j);
252 }
253 Index bi, bj;
254 for (const Block &b : correlations_) {
255 std::tie(bi, bj) = b.get_indices();
256 if (((i == bi) && (j == bj)) || ((i == -1) && (j == bj)) ||
257 ((i == -1) && (j == -1))) {
258 return &b;
259 }
260 }
261 return nullptr;
262}
263
265 const ArrayOfArrayOfIndex &jis) const {
266 for (Index i = 0; i < static_cast<Index>(jis.size()); ++i) {
267 Index n_blocks = 0;
268 for (const Block &b : correlations_) {
269 if (b.get_indices() == std::make_pair(i, i)) {
270 ++n_blocks;
271 }
272 }
273 if (n_blocks != 1) {
274 return false;
275 }
276 }
277 return true;
278}
279
281 auto pred = [&jis](const Block &b) {
282 Index i, j;
283 std::tie(i, j) = b.get_indices();
284
285 Index row_start = jis[i][0];
286 Index row_extent = jis[i][1] - jis[i][0] + 1;
287 Range row_range = b.get_row_range();
288 if ((row_range.offset != row_start) ||
289 (row_range.extent != row_extent)) {
290 return false;
291 }
292
293 Index column_start = jis[j][0];
294 Index column_extent = jis[j][1] - jis[j][0] + 1;
295 Range column_range = b.get_column_range();
296 if ((column_range.offset != column_start) ||
297 (column_range.extent != column_extent)) {
298 return false;
299 }
300 return true;
301 };
302
303 if (!std::all_of(correlations_.begin(), correlations_.end(), pred)) {
304 return false;
305 }
306 if (!std::all_of(inverses_.begin(), inverses_.end(), pred)) {
307 return false;
308 }
309 return true;
310}
311
313 Index i, j;
314 std::tie(i, j) = b.get_indices();
315
316 for (const Block &c : correlations_) {
317 Index ii, jj;
318 std::tie(ii, jj) = c.get_indices();
319
320 if ((ii == i) && (c.nrows() != b.nrows())) {
321 return false;
322 }
323
324 if ((jj == j) && (c.ncols() != b.ncols())) {
325 return false;
326 }
327
328 if ((ii == i) && (jj == j)) {
329 return false;
330 }
331 }
332 return true;
333}
334
336 for (const Block &b : inverses_) {
337 if (indices == b.get_indices()) {
338 return true;
339 }
340 }
341 return false;
342}
343
345 std::vector<std::vector<const Block *>> &corr_blocks) const {
346 for (size_t i = 0; i < correlations_.size(); i++) {
347 Index ci, cj;
348 std::tie(ci, cj) = correlations_[i].get_indices();
349 }
350
351 std::vector<bool> has_blocks(correlations_.size(), false);
352 std::queue<Index> rq_queue{};
353 for (size_t i = 0; i < correlations_.size(); ++i) {
354 if (!has_blocks[i]) {
355 Index ci, cj;
356 std::tie(ci, cj) = correlations_[i].get_indices();
357 rq_queue.push(ci);
358 if (ci != cj) {
359 rq_queue.push(cj);
361 has_blocks[i] = true;
362 corr_blocks.push_back(std::vector<const Block *>{&correlations_[i]});
364 while (!rq_queue.empty()) {
365 Index rq_index = rq_queue.front();
366 rq_queue.pop();
368 for (size_t j = 0; j < correlations_.size(); ++j) {
369 if (!has_blocks[j]) {
370 std::tie(ci, cj) = correlations_[j].get_indices();
371 if ((ci == rq_index) || (cj == rq_index)) {
372 if (ci != rq_index) {
373 rq_queue.push(ci);
374 }
375 if (cj != rq_index) {
376 rq_queue.push(cj);
377 }
378 corr_blocks.back().push_back(&correlations_[j]);
379 has_blocks[j] = true;
380 }
381 }
382 }
383 }
384 }
385 }
386}
387
389 std::vector<std::vector<const Block *>> correlation_blocks{};
390 generate_blocks(correlation_blocks);
391 for (std::vector<const Block *> &cb : correlation_blocks) {
393 }
394}
395
397 std::vector<Block> &inverses, std::vector<const Block *> &blocks) const {
398 // Can't compute inverse of empty block.
399 ARTS_ASSERT(blocks.size() > 0);
400
401 // Sort blocks w.r.t. indices.
402 auto comp = [](const Block *a, const Block *b) {
403 Index a1, a2, b1, b2;
404 std::tie(a1, a2) = a->get_indices();
405 std::tie(b1, b2) = b->get_indices();
406 return ((a1 < b1) || ((a1 == b1) && (a2 < b2)));
407 };
408
409 std::sort(blocks.begin(), blocks.end(), comp);
410
411 auto block_has_inverse = [this](const Block *a) {
412 return has_inverse(a->get_indices());
413 };
414 if (std::all_of(blocks.begin(), blocks.end(), block_has_inverse)) return;
415
416 // Otherwise go on to precompute the inverse of a block consisting
417 // of correlations between multiple retrieval quantities.
418
419 // The single blocks corresponding to a set of correlated retrieval quantities
420 // can be distributed freely over the covariance matrix, so we need to establish
421 // a mapping to a continuous square matrix to compute the inverse. This is done by
422 // mapping the coordinates of each block to a start row and extent in the continuous
423 // matrix A. We also record which retrieval quantity indices belong to this
424 // set of retrieval quantities.
425 Index n = 0;
426 std::map<Index, Index> block_start{};
427 std::map<Index, Index> block_extent{};
428 std::map<Index, Index> block_start_cont{};
429 std::map<Index, Index> block_extent_cont{};
430 std::vector<Index> block_indices{};
431
432 for (size_t i = 0; i < blocks.size(); ++i) {
433 Index ci, cj;
434 std::tie(ci, cj) = blocks[i]->get_indices();
435
436 if (ci == cj) {
437 Index extent = blocks[i]->get_row_range().extent;
438 block_start.insert(
439 std::make_pair(ci, blocks[i]->get_row_range().offset));
440 block_extent.insert(
441 std::make_pair(ci, blocks[i]->get_row_range().extent));
442 block_start_cont.insert(std::make_pair(ci, n));
443 block_extent_cont.insert(std::make_pair(ci, extent));
444 block_indices.push_back(ci);
445 n += extent;
446 }
447 }
448
449 // Copy blocks into a single dense matrix.
450 Matrix A(n, n);
451 A = 0.0;
452
453 for (size_t i = 0; i < blocks.size(); ++i) {
454 Index ci, cj;
455 std::tie(ci, cj) = blocks[i]->get_indices();
456 Range row_range(block_start_cont[ci], block_extent_cont[ci]);
457 Range column_range(block_start_cont[cj], block_extent_cont[cj]);
458 MatrixView A_view = A(row_range, column_range);
459
460 if (blocks[i]->get_matrix_type() == Block::MatrixType::dense) {
461 A_view = blocks[i]->get_dense();
462 } else {
463 A_view = static_cast<const Matrix>(blocks[i]->get_sparse());
464 }
465 }
466
467 for (Index i = 0; i < n; ++i) {
468 for (Index j = i + 1; j < n; ++j) {
469 A(j, i) = A(i, j);
470 }
471 }
472 inv(A, A);
473
474 // // Invert matrix using LAPACK.
475 // char uplo = 'L';
476 // int ni, info1(0), info2(0);
477 // ni = static_cast<int>(n);
478 // lapack::dpotrf_(&uplo, &ni, A.get_raw_data(), &ni, &info1);
479 // lapack::dpotri_(&uplo, &ni, A.get_raw_data(), &ni, &info2);
480 // if ((info1 != 0) || info2 !=0) {
481 // throw std::runtime_error("Error inverting block of covariance matrix."
482 // "Make sure that it is symmetric, positive definite"
483 // "or provide the inverse manually.");
484 // }
485
486 // Now we need to disassemble the matrix inverse bach to the separate block in the
487 // covariance matrix. Note, however, that blocks that previously were implicitly
488 // zero are now non-zero, i.e. the inverse may contain more blocks than the covariance
489 // matrix itself.
490 for (Index bi : block_indices) {
491 for (Index bj : block_indices) {
492 if (bi <= bj) {
493 Range row_range_A(block_start_cont[bi], block_extent_cont[bi]);
494 Range column_range_A(block_start_cont[bj], block_extent_cont[bj]);
495 Range row_range(block_start[bi], block_extent[bi]);
496 Range column_range(block_start[bj], block_extent[bj]);
497 MatrixView A_view = A(row_range_A, column_range_A);
498 inverses.push_back(Block(row_range,
499 column_range,
500 std::make_pair(bi, bj),
501 std::make_shared<Matrix>(A_view)));
502 }
503 }
504 }
505}
506
508
510 inverses_.push_back(c);
511}
512
514 Vector diag(nrows());
515 for (const Block &b : correlations_) {
516 Index i, j;
517 tie(i, j) = b.get_indices();
518
519 if (i == j) {
520 diag[b.get_row_range()] = b.diagonal();
521 }
522 }
523 return diag;
524}
525
528
529 Vector diag(nrows());
530 for (const Block &b : inverses_) {
531 Index i, j;
532 tie(i, j) = b.get_indices();
533
534 if (i == j) {
535 diag[b.get_row_range()] = b.diagonal();
536 }
537 }
538 return diag;
539}
540
541void mult(MatrixView C, ConstMatrixView A, const CovarianceMatrix &B) {
542 C = 0.0;
543 Matrix T(C);
544 for (const Block &c : B.correlations_) {
545 T = 0.0;
546 mult(T, A, c);
547 C += T;
548 }
549}
550
551void mult(MatrixView C, const CovarianceMatrix &A, ConstMatrixView B) {
552 C = 0.0;
553 Matrix T(C);
554 for (const Block &c : A.correlations_) {
555 T = 0.0;
556 mult(T, c, B);
557 C += T;
558 }
559}
560
561void mult(VectorView w, const CovarianceMatrix &A, ConstVectorView v) {
562 w = 0.0;
563 Vector t(w);
564 for (const Block &c : A.correlations_) {
565 t = 0.0;
566 mult(t, c, v);
567 w += t;
568 }
569}
570
571void mult_inv(MatrixView C, ConstMatrixView A, const CovarianceMatrix &B) {
572 C = 0.0;
573 Matrix T(C);
574 for (const Block &c : B.inverses_) {
575 T = 0.0;
576 mult(T, A, c);
577 C += T;
578 }
579}
580
581void mult_inv(MatrixView C, const CovarianceMatrix &A, ConstMatrixView B) {
582 C = 0.0;
583 Matrix T(C);
584 for (const Block &c : A.inverses_) {
585 T = 0.0;
586 mult(T, c, B);
587 C += T;
588 }
589}
590
591void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v) {
592 w = 0.0;
593 Vector t(w);
594 for (const Block &c : A.inverses_) {
595 t = 0.0;
596 mult(t, c, v);
597 w += t;
598 }
599}
600
601MatrixView operator+=(MatrixView A, const CovarianceMatrix &B) {
602 for (const Block &c : B.correlations_) {
603 A += c;
604 }
605 return A;
606}
607
608void add_inv(MatrixView A, const CovarianceMatrix &B) {
609 for (const Block &c : B.inverses_) {
610 A += c;
611 }
612}
613
614std::ostream &operator<<(std::ostream &os, const CovarianceMatrix &covmat) {
615 os << "Covariance Matrix, ";
616 os << "\tDimensions: [" << covmat.nrows() << " x " << covmat.ncols() << "]"
617 << std::endl;
618 os << "Blocks:" << std::endl;
619 for (const Block &b : covmat.correlations_) {
620 Index i, j;
621 tie(i, j) = b.get_indices();
622 os << "\ti = " << i << ", j = " << j << ": "
623 << b.get_row_range().extent;
624 os << " x " << b.get_column_range().extent;
625 os << ", has inverse: "
626 << (covmat.has_inverse(std::make_pair(i, j)) ? "yes" : "no");
627 os << std::endl;
628 }
629 return os;
630}
MatrixType get_matrix_type() const
const Sparse & get_sparse() const
const Matrix & get_dense() const
IndexPair get_indices() const
MatrixType matrix_type_
Range get_row_range() const
Range get_column_range() const
std::shared_ptr< Matrix > dense_
std::shared_ptr< Sparse > sparse_
std::vector< Block > correlations_
void add_correlation(Block c)
Add block to covariance matrix.
void compute_inverse() const
Compute the inverse of this correlation matrix.
bool has_block(Index i, Index j)
Check if the block with indices (i, j) is contained in the covariance matrix.
const Block * get_block(Index i=-1, Index j=-1)
Return a pointer to the block with indices (i,j).
bool has_inverse(IndexPair indices) const
Matrix get_inverse() const
bool is_consistent(const ArrayOfArrayOfIndex &jis) const
Checks that the dimensions of the blocks in the covariance matrix are consistent with ranges occupied...
std::vector< Block > inverses_
void generate_blocks(std::vector< std::vector< const Block * > > &) const
void invert_correlation_block(std::vector< Block > &inverses, std::vector< const Block * > &blocks) const
Vector diagonal() const
Diagonal elements as vector.
bool has_diagonal_blocks(const ArrayOfArrayOfIndex &jis) const
Checks that the covariance matrix contains one diagonal block per retrieval quantity.
void add_correlation_inverse(Block c)
Add block inverse to covariance matrix.
Index ndiagblocks() const
! The number of diagonal blocks in the matrix excluding inverses.
Vector inverse_diagonal() const
Diagonal of the inverse of the covariance matrix as vector.
Index nblocks() const
! The number of blocks in the matrix excluding inverses.
Index ninvdiagblocks() const
! The number of inverse diagonal blocks in the matrix.
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)
std::ostream & operator<<(std::ostream &os, const CovarianceMatrix &covmat)
void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v)
MatrixView operator+=(MatrixView A, const Block &B)
Header files of CovarianceMatrix class.
std::pair< Index, Index > IndexPair
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
Interface for the LAPACK library.
#define v
#define w
#define a
#define c
#define b