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