17#include "matpack_math.h"
22void mult(MatrixView C, ConstMatrixView A,
const Block &B) {
45 Matrix D(CTView.ncols(), CTView.nrows());
47 CTView += transpose(D);
52void mult(MatrixView C,
const Block &A, ConstMatrixView B) {
75 Matrix D(CTView.ncols(), CTView.nrows());
77 CTView += transpose(D);
99 transpose_mult(wtview, *A.
sparse_, vtview);
110 Aview +=
static_cast<const Matrix
>(B.
get_sparse());
120 ATview += transpose(
static_cast<const Matrix
>(B.
get_sparse()));
129CovarianceMatrix::operator Matrix()
const {
134 for (
const Block &
c : correlations_) {
135 MatrixView Aview = A(
c.get_row_range(),
c.get_column_range());
137 Aview =
c.get_dense();
139 Aview =
static_cast<const Matrix
>(
c.get_sparse());
143 std::tie(ci, cj) =
c.get_indices();
145 MatrixView ATview = A(
c.get_column_range(),
c.get_row_range());
147 ATview = transpose(
c.get_dense());
149 ATview = transpose(
static_cast<const Matrix
>(
c.get_sparse()));
162 MatrixView Aview = A(
c.get_row_range(),
c.get_column_range());
164 Aview =
c.get_dense();
166 Aview =
static_cast<const Matrix
>(
c.get_sparse());
170 std::tie(ci, cj) =
c.get_indices();
172 MatrixView ATview = A(
c.get_column_range(),
c.get_row_range());
174 ATview = transpose(
c.get_dense());
176 ATview = transpose(
static_cast<const Matrix
>(
c.get_sparse()));
188 std::tie(i, j) =
c.get_indices();
197 std::tie(i, j) =
c.get_indices();
203 return std::max(m1, m2);
213 std::tie(i, j) =
c.get_indices();
226 std::tie(i, j) =
c.get_indices();
243 result |=
b.get_indices() == std::make_pair(i, j);
255 std::tie(bi, bj) =
b.get_indices();
256 if (((i == bi) && (j == bj)) || ((i == -1) && (j == bj)) ||
257 ((i == -1) && (j == -1))) {
266 for (Index i = 0; i < static_cast<Index>(jis.size()); ++i) {
269 if (
b.get_indices() == std::make_pair(i, i)) {
281 auto pred = [&jis](
const Block &
b) {
283 std::tie(i, j) =
b.get_indices();
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)) {
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)) {
314 std::tie(i, j) =
b.get_indices();
318 std::tie(ii, jj) =
c.get_indices();
320 if ((ii == i) && (
c.nrows() !=
b.nrows())) {
324 if ((jj == j) && (
c.ncols() !=
b.ncols())) {
328 if ((ii == i) && (jj == j)) {
337 if (indices ==
b.get_indices()) {
345 std::vector<std::vector<const Block *>> &corr_blocks)
const {
352 std::queue<Index> rq_queue{};
354 if (!has_blocks[i]) {
361 has_blocks[i] =
true;
364 while (!rq_queue.empty()) {
365 Index rq_index = rq_queue.front();
369 if (!has_blocks[j]) {
371 if ((ci == rq_index) || (cj == rq_index)) {
372 if (ci != rq_index) {
375 if (cj != rq_index) {
379 has_blocks[j] =
true;
389 std::vector<std::vector<const Block *>> correlation_blocks{};
391 for (std::vector<const Block *> &cb : correlation_blocks) {
397 std::vector<Block> &inverses, std::vector<const Block *> &blocks)
const {
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)));
409 std::sort(blocks.begin(), blocks.end(), comp);
411 auto block_has_inverse = [
this](
const Block *
a) {
414 if (std::all_of(blocks.begin(), blocks.end(), block_has_inverse))
return;
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{};
432 for (
size_t i = 0; i < blocks.size(); ++i) {
434 std::tie(ci, cj) = blocks[i]->get_indices();
437 Index extent = blocks[i]->get_row_range().extent;
439 std::make_pair(ci, blocks[i]->get_row_range().offset));
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);
453 for (
size_t i = 0; i < blocks.size(); ++i) {
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);
461 A_view = blocks[i]->get_dense();
463 A_view =
static_cast<const Matrix
>(blocks[i]->get_sparse());
467 for (Index i = 0; i < n; ++i) {
468 for (Index j = i + 1; j < n; ++j) {
490 for (Index bi : block_indices) {
491 for (Index bj : block_indices) {
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,
500 std::make_pair(bi, bj),
501 std::make_shared<Matrix>(A_view)));
514 Vector diag(
nrows());
517 tie(i, j) =
b.get_indices();
520 diag[
b.get_row_range()] =
b.diagonal();
529 Vector diag(
nrows());
532 tie(i, j) =
b.get_indices();
535 diag[
b.get_row_range()] =
b.diagonal();
615 os <<
"Covariance Matrix, ";
616 os <<
"\tDimensions: [" << covmat.
nrows() <<
" x " << covmat.
ncols() <<
"]"
618 os <<
"Blocks:" << std::endl;
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");
MatrixType get_matrix_type() const
const Sparse & get_sparse() const
const Matrix & get_dense() const
IndexPair get_indices() const
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,...)
Interface for the LAPACK library.