145CovarianceMatrix::operator
Matrix()
const {
150 for (
const Block &
c : correlations_) {
151 MatrixView Aview = A(
c.get_row_range(),
c.get_column_range());
153 Aview =
c.get_dense();
155 Aview =
static_cast<const Matrix>(
c.get_sparse());
159 std::tie(ci, cj) =
c.get_indices();
161 MatrixView ATview = A(
c.get_column_range(),
c.get_row_range());
178 MatrixView Aview = A(
c.get_row_range(),
c.get_column_range());
180 Aview =
c.get_dense();
182 Aview =
static_cast<const Matrix>(
c.get_sparse());
186 std::tie(ci, cj) =
c.get_indices();
188 MatrixView ATview = A(
c.get_column_range(),
c.get_row_range());
204 std::tie(i, j) =
c.get_indices();
213 std::tie(i, j) =
c.get_indices();
219 return std::max(m1, m2);
229 std::tie(i, j) =
c.get_indices();
242 std::tie(i, j) =
c.get_indices();
259 result |=
b.get_indices() == std::make_pair(i, j);
271 std::tie(bi, bj) =
b.get_indices();
272 if (((i == bi) && (j == bj)) || ((i == -1) && (j == bj)) ||
273 ((i == -1) && (j == -1))) {
282 for (
Index i = 0; i < static_cast<Index>(jis.size()); ++i) {
285 if (
b.get_indices() == std::make_pair(i, i)) {
297 auto pred = [&jis](
const Block &
b) {
299 std::tie(i, j) =
b.get_indices();
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) ||
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)) {
330 std::tie(i, j) =
b.get_indices();
334 std::tie(ii, jj) =
c.get_indices();
336 if ((ii == i) && (
c.nrows() !=
b.nrows())) {
340 if ((jj == j) && (
c.ncols() !=
b.ncols())) {
344 if ((ii == i) && (jj == j)) {
353 if (indices ==
b.get_indices()) {
361 std::vector<std::vector<const Block *>> &corr_blocks)
const {
368 std::queue<Index> rq_queue{};
370 if (!has_blocks[i]) {
377 has_blocks[i] =
true;
380 while (!rq_queue.empty()) {
385 if (!has_blocks[j]) {
387 if ((ci == rq_index) || (cj == rq_index)) {
388 if (ci != rq_index) {
391 if (cj != rq_index) {
395 has_blocks[j] =
true;
405 std::vector<std::vector<const Block *>> correlation_blocks{};
407 for (std::vector<const Block *> &cb : correlation_blocks) {
413 std::vector<Block> &inverses, std::vector<const Block *> &blocks)
const {
420 std::tie(
a1,
a2) =
a->get_indices();
421 std::tie(
b1,
b2) =
b->get_indices();
425 std::sort(blocks.begin(), blocks.end(), comp);
427 auto block_has_inverse = [
this](
const Block *
a) {
430 if (std::all_of(blocks.begin(), blocks.end(), block_has_inverse))
return;
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{};
448 for (
size_t i = 0; i < blocks.size(); ++i) {
450 std::tie(ci, cj) = blocks[i]->get_indices();
453 Index extent = blocks[i]->get_row_range().get_extent();
455 std::make_pair(ci, blocks[i]->get_row_range().get_start()));
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);
469 for (
size_t i = 0; i < blocks.size(); ++i) {
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);
477 A_view = blocks[i]->get_dense();
479 A_view =
static_cast<const Matrix>(blocks[i]->get_sparse());
483 for (
Index i = 0; i < n; ++i) {
484 for (
Index j = i + 1; j < n; ++j) {
506 for (
Index bi : block_indices) {
507 for (
Index bj : block_indices) {
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,
516 std::make_pair(bi, bj),
517 std::make_shared<Matrix>(A_view)));
533 tie(i, j) =
b.get_indices();
536 diag[
b.get_row_range()] =
b.diagonal();
548 tie(i, j) =
b.get_indices();
551 diag[
b.get_row_range()] =
b.diagonal();
631 os <<
"Covariance Matrix, ";
632 os <<
"\tDimensions: [" << covmat.
nrows() <<
" x " << covmat.
ncols() <<
"]"
634 os <<
"Blocks:" << std::endl;
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");
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_
A constant view of a Matrix.
Index nrows() const noexcept
Index ncols() const noexcept
A constant view of a Vector.
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.
constexpr Index get_extent() const noexcept
Returns the extent of the range.
constexpr Index get_start() const noexcept
Returns the start index of the range.
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,...)
Interface for the LAPACK library.
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
void transpose_mult(VectorView y, const Sparse &M, ConstVectorView x)
Sparse matrix - Vector multiplication.
INDEX Index
The type to use for all integer numbers and indices.
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.