ARTS 2.5.10 (git: 2f1c442c)
covariance_matrix.h
Go to the documentation of this file.
1
2/* Copyright (C) 2017
3 Simon Pfreundschuh <simonpf@chalmers.se>
4
5 This program is free software; you can redistribute it and/or modify it
6 under the terms of the GNU General Public License as published by the
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18 USA. */
19
34#ifndef covariance_matrix_h
35#define covariance_matrix_h
36
37#include <memory>
38#include <utility>
39
40#include "jacobian.h"
41#include "matpackI.h"
42#include "matpackII.h"
43
45
46//------------------------------------------------------------------------------
47// Type Aliases
48//------------------------------------------------------------------------------
49
50using IndexPair = std::pair<Index, Index>;
51
52//------------------------------------------------------------------------------
53// Block
54//------------------------------------------------------------------------------
64class Block {
65 public:
66 /*
67 * Enum class representing the type of the matrix representing the block, i.e.
68 * whether its of type Matrix (dense) or of type Sparse(sparse).
69 */
70 enum class MatrixType { dense, sparse };
71
72 /*
73 * Create a correlation block from given row_range, column_range and shared_ptr
74 * to a dense matrix.
75 *
76 * @param row_range The range of element-rows covered by the block
77 * @param column_range The range of element-column covered by the block
78 * @param indices The pair of block-indices that identifies the block in the
79 * covariance matrix.
80 * @param dense A shared pointer to a den matrix of type Matrix
81 */
82 Block(Range row_range,
83 Range column_range,
84 IndexPair indices,
85 std::shared_ptr<Matrix> dense)
86 : row_range_(row_range),
87 column_range_(column_range),
88 indices_(std::move(indices)),
90 dense_(std::move(dense)),
91 sparse_(nullptr) {
92 // Nothing to do here.
93 }
94
95 /*
96 * Create a correlation block from given row_range, column_range and shared_ptr
97 * to a sparse matrix.
98 *
99 * @param row_range The range of element-rows covered by the block
100 * @param column_range The range of element-column covered by the block
101 * @param indices The pair of block-indices that identifies the block in the
102 * covariance matrix.
103 * @param dense A shared pointer to a den matrix of type Sparse
104 */
105 Block(Range row_range,
106 Range column_range,
107 IndexPair indices,
108 std::shared_ptr<Sparse> sparse)
109 : row_range_(row_range),
110 column_range_(column_range),
111 indices_(std::move(indices)),
113 dense_(nullptr),
114 sparse_(std::move(sparse)) {
115 // Nothing to do here.
116 }
117
118 Block(const Block &) = default;
119 Block(Block &&) = default;
120 Block &operator=(const Block &) = default;
121 Block &operator=(Block &&) = default;
122
123 ~Block() = default;
124
126 [[nodiscard]] Index nrows() const { return row_range_.get_extent(); }
128 [[nodiscard]] Index ncols() const { return column_range_.get_extent(); }
129
131 [[nodiscard]] Range get_row_range() const { return row_range_; }
133 [[nodiscard]] Range get_column_range() const { return column_range_; }
134
139
140 void set_matrix(std::shared_ptr<Sparse> sparse) { sparse_ = std::move(sparse); dense_ = nullptr; matrix_type_=MatrixType::sparse; }
141 void set_matrix(std::shared_ptr<Matrix> dense) { dense_ = std::move(dense); sparse_ = nullptr; matrix_type_=MatrixType::dense; }
142
144 [[nodiscard]] Vector diagonal() const {
145 if (dense_) {
146 return dense_->diagonal();
147 }
148 return sparse_->diagonal();
149
150 }
151
153 [[nodiscard]] IndexPair get_indices() const { return indices_; }
154
156 void set_indices(Index f, Index s) { indices_ = {f, s}; }
157
159 [[nodiscard]] MatrixType get_matrix_type() const { return matrix_type_; }
160
161 [[nodiscard]] const Matrix &get_dense() const {
163 return *dense_;
164 }
165
168 return *dense_;
169 }
170
171 [[nodiscard]] const Sparse &get_sparse() const {
173 return *sparse_;
174 }
175
178 return *sparse_;
179 }
180
181 // Friend declarations.
182 friend void mult(MatrixView, ConstMatrixView, const Block &);
183 friend void mult(MatrixView, const Block &, ConstMatrixView);
184 friend void mult(VectorView, const Block &, ConstVectorView);
185
186 friend MatrixView &operator+=(MatrixView &, const Block &);
187
188 private:
191
193
194 std::shared_ptr<Matrix> dense_;
195 std::shared_ptr<Sparse> sparse_;
196};
197
198void mult(MatrixView, ConstMatrixView, const Block &);
199void mult(MatrixView, const Block &, ConstMatrixView);
200void mult(VectorView, const Block &, ConstVectorView);
201
203void add_inv(MatrixView A, const Block &);
204
205//------------------------------------------------------------------------------
206// Covariance Matrices
207//------------------------------------------------------------------------------
229 public:
230 CovarianceMatrix() = default;
235
236 ~CovarianceMatrix() = default;
237
238 explicit operator Matrix() const;
239 Matrix get_inverse() const;
240
241 Index nrows() const;
242 Index ncols() const;
243
245 Index ndiagblocks() const;
246
248 Index ninvdiagblocks() const;
249
251 Index nblocks() const;
252
260 bool has_block(Index i, Index j);
261
274 const Block *get_block(Index i = -1, Index j = -1);
275
281 std::vector<Block> &get_blocks() { return correlations_; };
282
288 std::vector<Block> &get_inverse_blocks() { return inverses_; };
289
301 bool has_diagonal_blocks(const ArrayOfArrayOfIndex &jis) const;
302
315 bool is_consistent(const ArrayOfArrayOfIndex &jis) const;
316
327 bool is_consistent(const Block &block) const;
328
334 void compute_inverse() const;
335
344 void add_correlation(Block c);
345
355
363 Vector diagonal() const;
364
373 Vector inverse_diagonal() const;
374
375 // Friend declarations.
376 friend void mult(MatrixView, ConstMatrixView, const CovarianceMatrix &);
377 friend void mult(MatrixView, const CovarianceMatrix &, ConstMatrixView);
378 friend void mult(VectorView, const CovarianceMatrix &, ConstVectorView);
379
382 friend void solve(VectorView, const CovarianceMatrix &, ConstVectorView);
383
385 friend void add_inv(MatrixView, const CovarianceMatrix &);
386
387 friend void xml_read_from_stream(istream &,
389 bifstream *,
390 const Verbosity &);
391 friend void xml_write_to_stream(ostream &,
392 const CovarianceMatrix &,
393 bofstream *,
394 const String &,
395 const Verbosity &);
396 friend std::ostream &operator<<(std::ostream &os, const CovarianceMatrix &v);
397
398 private:
399 void generate_blocks(std::vector<std::vector<const Block *>> &) const;
400 void invert_correlation_block(std::vector<Block> &inverses,
401 std::vector<const Block *> &blocks) const;
402 bool has_inverse(IndexPair indices) const;
403
404 std::vector<Block> correlations_;
405 mutable std::vector<Block> inverses_;
406};
407
411
415
417void add_inv(MatrixView, const CovarianceMatrix &);
418
419std::ostream &operator<<(std::ostream &os, const ConstVectorView &v);
420
421#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()
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)
friend MatrixView & operator+=(MatrixView &, const Block &)
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()
A constant view of a Matrix.
Definition: matpackI.h:1065
A constant view of a Vector.
Definition: matpackI.h:521
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...
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.
friend MatrixView & operator+=(MatrixView &, const CovarianceMatrix &)
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
The VectorView class.
Definition: matpackI.h:674
The Vector class.
Definition: matpackI.h:910
Binary output file stream class.
Definition: bifstream.h:43
Binary output file stream class.
Definition: bofstream.h:42
void solve(VectorView, const CovarianceMatrix &, ConstVectorView)
MatrixView & operator+=(MatrixView &, const Block &)
std::ostream & operator<<(std::ostream &os, const ConstVectorView &v)
Definition: matpackI.cc:108
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:102
Routines for setting up the jacobian.
Header file for sparse matrices.
Implementation of Matrix, Vector, and such stuff.
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
The Sparse class.
Definition: matpackII.h:75
#define v
#define c