ARTS 2.5.11 (git: 725533f0)
covariance_matrix.h
Go to the documentation of this file.
1
15#ifndef covariance_matrix_h
16#define covariance_matrix_h
17
18#include <memory>
19#include <utility>
20
21#include "jacobian.h"
22#include "matpack_data.h"
23#include "matpack_math.h"
24#include "matpack_sparse.h"
25
27
28//------------------------------------------------------------------------------
29// Type Aliases
30//------------------------------------------------------------------------------
31
32using IndexPair = std::pair<Index, Index>;
33
34//------------------------------------------------------------------------------
35// Block
36//------------------------------------------------------------------------------
46class Block {
47 public:
48 /*
49 * Enum class representing the type of the matrix representing the block, i.e.
50 * whether its of type Matrix (dense) or of type Sparse(sparse).
51 */
52 enum class MatrixType { dense, sparse };
53
54 /*
55 * Create a correlation block from given row_range, column_range and shared_ptr
56 * to a dense matrix.
57 *
58 * @param row_range The range of element-rows covered by the block
59 * @param column_range The range of element-column covered by the block
60 * @param indices The pair of block-indices that identifies the block in the
61 * covariance matrix.
62 * @param dense A shared pointer to a den matrix of type Matrix
63 */
64 Block(Range row_range,
65 Range column_range,
66 IndexPair indices,
67 std::shared_ptr<Matrix> dense)
68 : row_range_(row_range),
69 column_range_(column_range),
70 indices_(std::move(indices)),
72 dense_(std::move(dense)),
73 sparse_(nullptr) {
74 // Nothing to do here.
75 }
76
77 /*
78 * Create a correlation block from given row_range, column_range and shared_ptr
79 * to a sparse matrix.
80 *
81 * @param row_range The range of element-rows covered by the block
82 * @param column_range The range of element-column covered by the block
83 * @param indices The pair of block-indices that identifies the block in the
84 * covariance matrix.
85 * @param dense A shared pointer to a den matrix of type Sparse
86 */
87 Block(Range row_range,
88 Range column_range,
89 IndexPair indices,
90 std::shared_ptr<Sparse> sparse)
91 : row_range_(row_range),
92 column_range_(column_range),
93 indices_(std::move(indices)),
95 dense_(nullptr),
96 sparse_(std::move(sparse)) {
97 // Nothing to do here.
98 }
99
100 Block(const Block &) = default;
101 Block(Block &&) = default;
102 Block &operator=(const Block &) = default;
103 Block &operator=(Block &&) = default;
104
105 ~Block() = default;
106
108 [[nodiscard]] Index nrows() const { return row_range_.extent; }
110 [[nodiscard]] Index ncols() const { return column_range_.extent; }
111
113 [[nodiscard]] Range get_row_range() const { return row_range_; }
115 [[nodiscard]] Range get_column_range() const { return column_range_; }
116
118 Range& get_row_range() { return row_range_; }
120 Range& get_column_range() { return column_range_; }
121
122 void set_matrix(std::shared_ptr<Sparse> sparse) { sparse_ = std::move(sparse); dense_ = nullptr; matrix_type_=MatrixType::sparse; }
123 void set_matrix(std::shared_ptr<Matrix> dense) { dense_ = std::move(dense); sparse_ = nullptr; matrix_type_=MatrixType::dense; }
124
126 [[nodiscard]] Vector diagonal() const {
127 if (dense_) {
128 return ::diagonal(*dense_);
129 }
130 return sparse_->diagonal();
131
132 }
133
135 [[nodiscard]] IndexPair get_indices() const { return indices_; }
136
138 void set_indices(Index f, Index s) { indices_ = {f, s}; }
139
141 [[nodiscard]] MatrixType get_matrix_type() const { return matrix_type_; }
142
143 [[nodiscard]] const Matrix &get_dense() const {
145 return *dense_;
146 }
147
148 Matrix &get_dense() {
150 return *dense_;
151 }
152
153 [[nodiscard]] const Sparse &get_sparse() const {
155 return *sparse_;
156 }
157
158 Sparse &get_sparse() {
160 return *sparse_;
161 }
162
163 // Friend declarations.
164 friend void mult(MatrixView, ConstMatrixView, const Block &);
165 friend void mult(MatrixView, const Block &, ConstMatrixView);
166 friend void mult(VectorView, const Block &, ConstVectorView);
167
168 friend MatrixView operator+=(MatrixView, const Block &);
169
170 private:
173
175
176 std::shared_ptr<Matrix> dense_;
177 std::shared_ptr<Sparse> sparse_;
178};
179
180void mult(MatrixView, ConstMatrixView, const Block &);
181void mult(MatrixView, const Block &, ConstMatrixView);
182void mult(VectorView, const Block &, ConstVectorView);
183
184MatrixView operator+=(MatrixView, const Block &);
185void add_inv(MatrixView A, const Block &);
186
187//------------------------------------------------------------------------------
188// Covariance Matrices
189//------------------------------------------------------------------------------
211 public:
212 CovarianceMatrix() = default;
217
218 ~CovarianceMatrix() = default;
219
220 explicit operator Matrix() const;
221 Matrix get_inverse() const;
222
223 Index nrows() const;
224 Index ncols() const;
225
227 Index ndiagblocks() const;
228
230 Index ninvdiagblocks() const;
231
233 Index nblocks() const;
234
242 bool has_block(Index i, Index j);
243
256 const Block *get_block(Index i = -1, Index j = -1);
257
263 std::vector<Block> &get_blocks() { return correlations_; };
264
270 std::vector<Block> &get_inverse_blocks() { return inverses_; };
271
283 bool has_diagonal_blocks(const ArrayOfArrayOfIndex &jis) const;
284
297 bool is_consistent(const ArrayOfArrayOfIndex &jis) const;
298
309 bool is_consistent(const Block &block) const;
310
316 void compute_inverse() const;
317
326 void add_correlation(Block c);
327
337
345 Vector diagonal() const;
346
355 Vector inverse_diagonal() const;
356
357 // Friend declarations.
358 friend void mult(MatrixView, ConstMatrixView, const CovarianceMatrix &);
359 friend void mult(MatrixView, const CovarianceMatrix &, ConstMatrixView);
360 friend void mult(VectorView, const CovarianceMatrix &, ConstVectorView);
361
362 friend void mult_inv(MatrixView, ConstMatrixView, const CovarianceMatrix &);
363 friend void mult_inv(MatrixView, const CovarianceMatrix &, ConstMatrixView);
364 friend void solve(VectorView, const CovarianceMatrix &, ConstVectorView);
365
366 friend MatrixView operator+=(MatrixView, const CovarianceMatrix &);
367 friend void add_inv(MatrixView, const CovarianceMatrix &);
368
369 friend void xml_read_from_stream(istream &,
371 bifstream *,
372 const Verbosity &);
373 friend void xml_write_to_stream(ostream &,
374 const CovarianceMatrix &,
375 bofstream *,
376 const String &,
377 const Verbosity &);
378 friend std::ostream &operator<<(std::ostream &os, const CovarianceMatrix &v);
379
380 private:
381 void generate_blocks(std::vector<std::vector<const Block *>> &) const;
382 void invert_correlation_block(std::vector<Block> &inverses,
383 std::vector<const Block *> &blocks) const;
384 bool has_inverse(IndexPair indices) const;
385
386 std::vector<Block> correlations_;
387 mutable std::vector<Block> inverses_;
388};
389
390void mult(MatrixView, ConstMatrixView, const CovarianceMatrix &);
391void mult(MatrixView, const CovarianceMatrix &, ConstMatrixView);
392void mult(VectorView, const CovarianceMatrix &, ConstVectorView);
393
394void mult_inv(MatrixView, ConstMatrixView, const CovarianceMatrix &);
395void mult_inv(MatrixView, const CovarianceMatrix &, ConstMatrixView);
396void solve(VectorView, const CovarianceMatrix &, ConstVectorView);
397
398MatrixView operator+=(MatrixView, const CovarianceMatrix &);
399void add_inv(MatrixView, const CovarianceMatrix &);
400
401#endif // covariance_matrix_h
MatrixType get_matrix_type() const
Block(Range row_range, Range column_range, IndexPair indices, std::shared_ptr< Sparse > sparse)
Sparse & get_sparse()
Block & operator=(const Block &)=default
const Sparse & get_sparse() const
const Matrix & get_dense() const
Block(Block &&)=default
~Block()=default
Block(Range row_range, Range column_range, IndexPair indices, std::shared_ptr< Matrix > dense)
Range & get_row_range()
friend MatrixView operator+=(MatrixView, const Block &)
Range row_range_
Range column_range_
IndexPair get_indices() const
Index nrows() const
void set_indices(Index f, Index s)
Vector diagonal() const
void set_matrix(std::shared_ptr< Sparse > sparse)
MatrixType matrix_type_
Range get_row_range() const
Range get_column_range() const
void set_matrix(std::shared_ptr< Matrix > dense)
IndexPair indices_
Matrix & get_dense()
Block & operator=(Block &&)=default
std::shared_ptr< Matrix > dense_
std::shared_ptr< Sparse > sparse_
Index ncols() const
friend void mult(MatrixView, ConstMatrixView, const Block &)
Block(const Block &)=default
Range & get_column_range()
friend void solve(VectorView, const CovarianceMatrix &, ConstVectorView)
std::vector< Block > correlations_
~CovarianceMatrix()=default
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).
friend void xml_read_from_stream(istream &, CovarianceMatrix &, bifstream *, const Verbosity &)
Reads CovarianceMatrix from XML input stream.
std::vector< Block > & get_inverse_blocks()
Blocks of the inverse covariance matrix.
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...
friend MatrixView operator+=(MatrixView, const CovarianceMatrix &)
CovarianceMatrix & operator=(CovarianceMatrix &&)=default
friend std::ostream & operator<<(std::ostream &os, const CovarianceMatrix &v)
std::vector< Block > inverses_
void generate_blocks(std::vector< std::vector< const Block * > > &) const
CovarianceMatrix & operator=(const CovarianceMatrix &)=default
CovarianceMatrix(CovarianceMatrix &&)=default
friend void mult_inv(MatrixView, ConstMatrixView, const CovarianceMatrix &)
void invert_correlation_block(std::vector< Block > &inverses, std::vector< const Block * > &blocks) const
friend void xml_write_to_stream(ostream &, const CovarianceMatrix &, bofstream *, const String &, const Verbosity &)
Write CovarianceMatrix to XML output stream.
CovarianceMatrix()=default
std::vector< Block > & get_blocks()
Block in the covariance matrix.
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.
friend void mult(MatrixView, ConstMatrixView, const CovarianceMatrix &)
friend void add_inv(MatrixView, const CovarianceMatrix &)
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.
CovarianceMatrix(const CovarianceMatrix &)=default
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.
Binary output file stream class.
Definition bifstream.h:26
Binary output file stream class.
Definition bofstream.h:25
void solve(VectorView, const CovarianceMatrix &, ConstVectorView)
MatrixView operator+=(MatrixView, const Block &)
void mult_inv(MatrixView, ConstMatrixView, const CovarianceMatrix &)
void add_inv(MatrixView A, const Block &)
std::pair< Index, Index > IndexPair
void mult(MatrixView, ConstMatrixView, const Block &)
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
Routines for setting up the jacobian.
#define v
#define c