ARTS  2.4.0(git:4fb77825)
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 //------------------------------------------------------------------------------
38 void mult(MatrixView C, ConstMatrixView A, const Block &B) {
39  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()));
44  ConstMatrixView ATView(A(joker, B.get_column_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 
68 void mult(MatrixView C, const Block &A, ConstMatrixView B) {
69  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));
73  ConstMatrixView BView(B(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 
98 void mult(VectorView w, const Block &A, ConstVectorView v) {
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 //------------------------------------------------------------------------------
145 CovarianceMatrix::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());
183  }
184 
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 
222 Index CovarianceMatrix::ncols() const { return nrows(); }
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 
240  if (i > j) {
241  std::swap(i, j);
242  }
243 
244  bool result = false;
245  for (const Block &b : correlations_) {
246  result |= b.get_indices() == std::make_pair(i, j);
247  }
248 
249  return result;
250 }
251 
253  if (i > j) {
254  std::swap(i, j);
255  }
256  Index bi, bj;
257  for (const Block &b : correlations_) {
258  std::tie(bi, bj) = b.get_indices();
259  if (((i == bi) && (j == bj)) || ((i == -1) && (j == bj)) ||
260  ((i == -1) && (j == -1))) {
261  return &b;
262  }
263  }
264  return nullptr;
265 }
266 
268  const ArrayOfArrayOfIndex &jis) const {
269  for (Index i = 0; i < static_cast<Index>(jis.size()); ++i) {
270  Index n_blocks = 0;
271  for (const Block &b : correlations_) {
272  if (b.get_indices() == std::make_pair(i, i)) {
273  ++n_blocks;
274  }
275  }
276  if (n_blocks != 1) {
277  return false;
278  }
279  }
280  return true;
281 }
282 
284  auto pred = [&jis](const Block &b) {
285  Index i, j;
286  std::tie(i, j) = b.get_indices();
287 
288  Index row_start = jis[i][0];
289  Index row_extent = jis[i][1] - jis[i][0] + 1;
290  Range row_range = b.get_row_range();
291  if ((row_range.get_start() != row_start) ||
292  (row_range.get_extent() != row_extent)) {
293  return false;
294  }
295 
296  Index column_start = jis[j][0];
297  Index column_extent = jis[j][1] - jis[j][0] + 1;
298  Range column_range = b.get_column_range();
299  if ((column_range.get_start() != column_start) ||
300  (column_range.get_extent() != column_extent)) {
301  return false;
302  }
303  return true;
304  };
305 
306  if (!std::all_of(correlations_.begin(), correlations_.end(), pred)) {
307  return false;
308  }
309  if (!std::all_of(inverses_.begin(), inverses_.end(), pred)) {
310  return false;
311  }
312  return true;
313 }
314 
316  Index i, j;
317  std::tie(i, j) = b.get_indices();
318 
319  for (const Block &c : correlations_) {
320  Index ii, jj;
321  std::tie(ii, jj) = c.get_indices();
322 
323  if ((ii == i) && (c.nrows() != b.nrows())) {
324  return false;
325  }
326 
327  if ((jj == j) && (c.ncols() != b.ncols())) {
328  return false;
329  }
330 
331  if ((ii == i) && (jj == j)) {
332  return false;
333  }
334  }
335  return true;
336 }
337 
339  for (const Block &b : inverses_) {
340  if (indices == b.get_indices()) {
341  return true;
342  }
343  }
344  return false;
345 }
346 
348  std::vector<std::vector<const Block *>> &corr_blocks) const {
349  for (size_t i = 0; i < correlations_.size(); i++) {
350  Index ci, cj;
351  std::tie(ci, cj) = correlations_[i].get_indices();
352  }
353 
354  std::vector<bool> has_blocks(correlations_.size(), false);
355  std::queue<Index> rq_queue{};
356  for (size_t i = 0; i < correlations_.size(); ++i) {
357  if (!has_blocks[i]) {
358  Index ci, cj;
359  std::tie(ci, cj) = correlations_[i].get_indices();
360  rq_queue.push(ci);
361  if (ci != cj) {
362  rq_queue.push(cj);
363  }
364  has_blocks[i] = true;
365  corr_blocks.push_back(std::vector<const Block *>{&correlations_[i]});
366 
367  while (!rq_queue.empty()) {
368  Index rq_index = rq_queue.front();
369  rq_queue.pop();
370 
371  for (size_t j = 0; j < correlations_.size(); ++j) {
372  if (!has_blocks[j]) {
373  std::tie(ci, cj) = correlations_[j].get_indices();
374  if ((ci == rq_index) || (cj == rq_index)) {
375  if (ci != rq_index) {
376  rq_queue.push(ci);
377  }
378  if (cj != rq_index) {
379  rq_queue.push(cj);
380  }
381  corr_blocks.back().push_back(&correlations_[j]);
382  has_blocks[j] = true;
383  }
384  }
385  }
386  }
387  }
388  }
389 }
390 
392  std::vector<std::vector<const Block *>> correlation_blocks{};
393  generate_blocks(correlation_blocks);
394  for (std::vector<const Block *> &cb : correlation_blocks) {
396  }
397 }
398 
400  std::vector<Block> &inverses, std::vector<const Block *> &blocks) const {
401  // Can't compute inverse of empty block.
402  assert(blocks.size() > 0);
403 
404  // Sort blocks w.r.t. indices.
405  auto comp = [](const Block *a, const Block *b) {
406  Index a1, a2, b1, b2;
407  std::tie(a1, a2) = a->get_indices();
408  std::tie(b1, b2) = b->get_indices();
409  return ((a1 < b1) || ((a1 == b1) && (a2 < b2)));
410  };
411 
412  std::sort(blocks.begin(), blocks.end(), comp);
413 
414  auto block_has_inverse = [this](const Block *a) {
415  return has_inverse(a->get_indices());
416  };
417  if (std::all_of(blocks.begin(), blocks.end(), block_has_inverse)) return;
418 
419  // Otherwise go on to precompute the inverse of a block consisting
420  // of correlations between multiple retrieval quantities.
421 
422  // The single blocks corresponding to a set of correlated retrieval quantities
423  // can be distributed freely over the covariance matrix, so we need to establish
424  // a mapping to a continuous square matrix to compute the inverse. This is done by
425  // mapping the coordinates of each block to a start row and extent in the continuous
426  // matrix A. We also record which retrieval quantity indices belong to this
427  // set of retrieval quantities.
428  Index n = 0;
429  std::map<Index, Index> block_start{};
430  std::map<Index, Index> block_extent{};
431  std::map<Index, Index> block_start_cont{};
432  std::map<Index, Index> block_extent_cont{};
433  std::vector<Index> block_indices{};
434 
435  for (size_t i = 0; i < blocks.size(); ++i) {
436  Index ci, cj;
437  std::tie(ci, cj) = blocks[i]->get_indices();
438 
439  if (ci == cj) {
440  Index extent = blocks[i]->get_row_range().get_extent();
441  block_start.insert(
442  std::make_pair(ci, blocks[i]->get_row_range().get_start()));
443  block_extent.insert(
444  std::make_pair(ci, blocks[i]->get_row_range().get_extent()));
445  block_start_cont.insert(std::make_pair(ci, n));
446  block_extent_cont.insert(std::make_pair(ci, extent));
447  block_indices.push_back(ci);
448  n += extent;
449  }
450  }
451 
452  // Copy blocks into a single dense matrix.
453  Matrix A(n, n);
454  A = 0.0;
455 
456  for (size_t i = 0; i < blocks.size(); ++i) {
457  Index ci, cj;
458  std::tie(ci, cj) = blocks[i]->get_indices();
459  Range row_range(block_start_cont[ci], block_extent_cont[ci]);
460  Range column_range(block_start_cont[cj], block_extent_cont[cj]);
461  MatrixView A_view = A(row_range, column_range);
462 
463  if (blocks[i]->get_matrix_type() == Block::MatrixType::dense) {
464  A_view = blocks[i]->get_dense();
465  } else {
466  A_view = static_cast<const Matrix>(blocks[i]->get_sparse());
467  }
468  }
469 
470  for (Index i = 0; i < n; ++i) {
471  for (Index j = i + 1; j < n; ++j) {
472  A(j, i) = A(i, j);
473  }
474  }
475  inv(A, A);
476 
477  // // Invert matrix using LAPACK.
478  // char uplo = 'L';
479  // int ni, info1(0), info2(0);
480  // ni = static_cast<int>(n);
481  // lapack::dpotrf_(&uplo, &ni, A.get_raw_data(), &ni, &info1);
482  // lapack::dpotri_(&uplo, &ni, A.get_raw_data(), &ni, &info2);
483  // if ((info1 != 0) || info2 !=0) {
484  // throw std::runtime_error("Error inverting block of covariance matrix."
485  // "Make sure that it is symmetric, positive definite"
486  // "or provide the inverse manually.");
487  // }
488 
489  // Now we need to disassemble the matrix inverse bach to the separate block in the
490  // covariance matrix. Note, however, that blocks that previously were implicitly
491  // zero are now non-zero, i.e. the inverse may contain more blocks than the covariance
492  // matrix itself.
493  for (Index bi : block_indices) {
494  for (Index bj : block_indices) {
495  if (bi <= bj) {
496  Range row_range_A(block_start_cont[bi], block_extent_cont[bi]);
497  Range column_range_A(block_start_cont[bj], block_extent_cont[bj]);
498  Range row_range(block_start[bi], block_extent[bi]);
499  Range column_range(block_start[bj], block_extent[bj]);
500  MatrixView A_view = A(row_range_A, column_range_A);
501  inverses.push_back(Block(row_range,
502  column_range,
503  std::make_pair(bi, bj),
504  std::make_shared<Matrix>(A_view)));
505  }
506  }
507  }
508 }
509 
511 
513  inverses_.push_back(c);
514 }
515 
517  Vector diag(nrows());
518  for (const Block &b : correlations_) {
519  Index i, j;
520  tie(i, j) = b.get_indices();
521 
522  if (i == j) {
523  diag[b.get_row_range()] = b.diagonal();
524  }
525  }
526  return diag;
527 }
528 
530  compute_inverse();
531 
532  Vector diag(nrows());
533  for (const Block &b : inverses_) {
534  Index i, j;
535  tie(i, j) = b.get_indices();
536 
537  if (i == j) {
538  diag[b.get_row_range()] = b.diagonal();
539  }
540  }
541  return diag;
542 }
543 
545  C = 0.0;
546  Matrix T(C);
547  for (const Block &c : B.correlations_) {
548  T = 0.0;
549  mult(T, A, c);
550  C += T;
551  }
552 }
553 
555  C = 0.0;
556  Matrix T(C);
557  for (const Block &c : A.correlations_) {
558  T = 0.0;
559  mult(T, c, B);
560  C += T;
561  }
562 }
563 
565  w = 0.0;
566  Vector t(w);
567  for (const Block &c : A.correlations_) {
568  t = 0.0;
569  mult(t, c, v);
570  w += t;
571  }
572 }
573 
575  C = 0.0;
576  Matrix T(C);
577  for (const Block &c : B.inverses_) {
578  T = 0.0;
579  mult(T, A, c);
580  C += T;
581  }
582 }
583 
585  C = 0.0;
586  Matrix T(C);
587  for (const Block &c : A.inverses_) {
588  T = 0.0;
589  mult(T, c, B);
590  C += T;
591  }
592 }
593 
595  w = 0.0;
596  Vector t(w);
597  for (const Block &c : A.inverses_) {
598  t = 0.0;
599  mult(t, c, v);
600  w += t;
601  }
602 }
603 
605  for (const Block &c : B.correlations_) {
606  A += c;
607  }
608  return A;
609 }
610 
612  for (const Block &c : B.inverses_) {
613  A += c;
614  }
615 }
616 
617 std::ostream &operator<<(std::ostream &os, const CovarianceMatrix &covmat) {
618  os << "Covariance Matrix, ";
619  os << "\tDimensions: [" << covmat.nrows() << " x " << covmat.ncols() << "]"
620  << std::endl;
621  os << "Blocks:" << std::endl;
622  for (const Block &b : covmat.correlations_) {
623  Index i, j;
624  tie(i, j) = b.get_indices();
625  os << "\ti = " << i << ", j = " << j << ": "
626  << b.get_row_range().get_extent();
627  os << " x " << b.get_column_range().get_extent();
628  os << ", has inverse: "
629  << (covmat.has_inverse(std::make_pair(i, j)) ? "yes" : "no");
630  os << std::endl;
631  }
632  return os;
633 }
Block::dense_
std::shared_ptr< Matrix > dense_
Definition: covariance_matrix.h:192
Matrix
The Matrix class.
Definition: matpackI.h:1193
CovarianceMatrix::ncols
Index ncols() const
Definition: covariance_matrix.cc:222
Block::get_dense
const Matrix & get_dense() const
Definition: covariance_matrix.h:159
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
Block::get_sparse
const Sparse & get_sparse() const
Definition: covariance_matrix.h:169
lapack.h
Interface for the LAPACK library.
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
operator+=
MatrixView & operator+=(MatrixView &A, const Block &B)
Definition: covariance_matrix.cc:120
joker
const Joker joker
Block::get_column_range
Range get_column_range() const
Definition: covariance_matrix.h:131
CovarianceMatrix::correlations_
std::vector< Block > correlations_
Definition: covariance_matrix.h:399
inv
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
oem::Matrix
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:34
swap
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
Definition: complex.cc:731
CovarianceMatrix
Definition: covariance_matrix.h:226
Block::get_matrix_type
MatrixType get_matrix_type() const
Definition: covariance_matrix.h:157
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
CovarianceMatrix::diagonal
Vector diagonal() const
Diagonal elements as vector.
Definition: covariance_matrix.cc:516
CovarianceMatrix::has_diagonal_blocks
bool has_diagonal_blocks(const ArrayOfArrayOfIndex &jis) const
Checks that the covariance matrix contains one diagonal block per retrieval quantity.
Definition: covariance_matrix.cc:267
b2
#define b2
Definition: complex.h:59
covariance_matrix.h
Header files of CovarianceMatrix class.
CovarianceMatrix::compute_inverse
void compute_inverse() const
Compute the inverse of this correlation matrix.
Definition: covariance_matrix.cc:391
Range::get_start
constexpr Index get_start() const
Returns the start index of the range.
Definition: matpackI.h:327
Block::MatrixType::dense
@ dense
operator<<
std::ostream & operator<<(std::ostream &os, const CovarianceMatrix &covmat)
Definition: covariance_matrix.cc:617
Array< ArrayOfIndex >
CovarianceMatrix::inverses_
std::vector< Block > inverses_
Definition: covariance_matrix.h:400
CovarianceMatrix::nblocks
Index nblocks() const
! The number of blocks in the matrix excluding inverses.
Definition: covariance_matrix.cc:237
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
CovarianceMatrix::get_block
const Block * get_block(Index i=-1, Index j=-1)
Return a pointer to the block with indices (i,j).
Definition: covariance_matrix.cc:252
CovarianceMatrix::get_inverse
Matrix get_inverse() const
Definition: covariance_matrix.cc:172
VectorView
The VectorView class.
Definition: matpackI.h:610
CovarianceMatrix::add_correlation_inverse
void add_correlation_inverse(Block c)
Add block inverse to covariance matrix.
Definition: covariance_matrix.cc:512
CovarianceMatrix::ndiagblocks
Index ndiagblocks() const
! The number of diagonal blocks in the matrix excluding inverses.
Definition: covariance_matrix.cc:224
Block::get_row_range
Range get_row_range() const
Definition: covariance_matrix.h:129
CovarianceMatrix::add_correlation
void add_correlation(Block c)
Add block to covariance matrix.
Definition: covariance_matrix.cc:510
CovarianceMatrix::generate_blocks
void generate_blocks(std::vector< std::vector< const Block * >> &) const
Definition: covariance_matrix.cc:347
Block::matrix_type_
MatrixType matrix_type_
Definition: covariance_matrix.h:190
max
#define max(a, b)
Definition: legacy_continua.cc:20629
a1
#define a1
Definition: complex.h:56
mult_inv
void mult_inv(MatrixView C, ConstMatrixView A, const CovarianceMatrix &B)
Definition: covariance_matrix.cc:574
IndexPair
std::pair< Index, Index > IndexPair
Definition: covariance_matrix.h:48
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
Range::get_extent
constexpr Index get_extent() const
Returns the extent of the range.
Definition: matpackI.h:329
ARTS::Var::nrows
Index nrows(Workspace &ws) noexcept
Definition: autoarts.h:4682
Range
The range class.
Definition: matpackI.h:160
CovarianceMatrix::has_block
bool has_block(Index i, Index j)
Check if the block with indices (i, j) is contained in the covariance matrix.
Definition: covariance_matrix.cc:239
transpose_mult
void transpose_mult(VectorView y, const Sparse &M, ConstVectorView x)
Sparse matrix - Vector multiplication.
Definition: matpackII.cc:452
Block::ncols
Index ncols() const
Definition: covariance_matrix.h:126
mult
void mult(MatrixView C, ConstMatrixView A, const Block &B)
Definition: covariance_matrix.cc:38
Block::nrows
Index nrows() const
Definition: covariance_matrix.h:124
add_inv
void add_inv(MatrixView A, const CovarianceMatrix &B)
Definition: covariance_matrix.cc:611
Block
Definition: covariance_matrix.h:62
a2
#define a2
Definition: complex.h:58
CovarianceMatrix::has_inverse
bool has_inverse(IndexPair indices) const
Definition: covariance_matrix.cc:338
CovarianceMatrix::nrows
Index nrows() const
Definition: covariance_matrix.cc:199
Block::sparse_
std::shared_ptr< Sparse > sparse_
Definition: covariance_matrix.h:193
CovarianceMatrix::is_consistent
bool is_consistent(const ArrayOfArrayOfIndex &jis) const
Checks that the dimensions of the blocks in the covariance matrix are consistent with ranges occupied...
Definition: covariance_matrix.cc:283
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Block::get_indices
IndexPair get_indices() const
Definition: covariance_matrix.h:151
transpose
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
Definition: complex.cc:1509
Vector
The Vector class.
Definition: matpackI.h:860
CovarianceMatrix::inverse_diagonal
Vector inverse_diagonal() const
Diagonal of the inverse of the covariance matrix as vector.
Definition: covariance_matrix.cc:529
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
b1
#define b1
Definition: complex.h:57
solve
void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v)
Definition: covariance_matrix.cc:594
CovarianceMatrix::invert_correlation_block
void invert_correlation_block(std::vector< Block > &inverses, std::vector< const Block * > &blocks) const
Definition: covariance_matrix.cc:399