ARTS 2.5.0 (git: 9ee3ac6c)
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
39#include "jacobian.h"
40#include "matpackI.h"
41#include "matpackII.h"
42
44
45//------------------------------------------------------------------------------
46// Type Aliases
47//------------------------------------------------------------------------------
48
49using IndexPair = std::pair<Index, Index>;
50
51//------------------------------------------------------------------------------
52// Block
53//------------------------------------------------------------------------------
63class Block {
64 public:
65 /*
66 * Enum class representing the type of the matrix representing the block, i.e.
67 * whether its of type Matrix (dense) or of type Sparse(sparse).
68 */
69 enum class MatrixType { dense, sparse };
70
71 /*
72 * Create a correlation block from given row_range, column_range and shared_ptr
73 * to a dense matrix.
74 *
75 * @param row_range The range of element-rows covered by the block
76 * @param column_range The range of element-column covered by the block
77 * @param indices The pair of block-indices that identifies the block in the
78 * covariance matrix.
79 * @param dense A shared pointer to a den matrix of type Matrix
80 */
81 Block(Range row_range,
82 Range column_range,
83 IndexPair indices,
84 std::shared_ptr<Matrix> dense)
85 : row_range_(row_range),
86 column_range_(column_range),
87 indices_(indices),
89 dense_(dense),
90 sparse_(nullptr) {
91 // Nothing to do here.
92 }
93
94 /*
95 * Create a correlation block from given row_range, column_range and shared_ptr
96 * to a sparse matrix.
97 *
98 * @param row_range The range of element-rows covered by the block
99 * @param column_range The range of element-column covered by the block
100 * @param indices The pair of block-indices that identifies the block in the
101 * covariance matrix.
102 * @param dense A shared pointer to a den matrix of type Sparse
103 */
104 Block(Range row_range,
105 Range column_range,
106 IndexPair indices,
107 std::shared_ptr<Sparse> sparse)
108 : row_range_(row_range),
109 column_range_(column_range),
110 indices_(indices),
111 matrix_type_(MatrixType::sparse),
112 dense_(nullptr),
113 sparse_(sparse) {
114 // Nothing to do here.
115 }
116
117 Block(const Block &) = default;
118 Block(Block &&) = default;
119 Block &operator=(const Block &) = default;
120 Block &operator=(Block &&) = default;
121
122 ~Block() = default;
123
125 Index nrows() const { return row_range_.get_extent(); }
127 Index ncols() const { return column_range_.get_extent(); }
128
130 Range get_row_range() const { return row_range_; }
133
138
139 void set_matrix(std::shared_ptr<Sparse> sparse) { sparse_ = sparse; }
140 void set_matrix(std::shared_ptr<Matrix> dense) { dense_ = dense; }
141
143 Vector diagonal() const {
144 if (dense_) {
145 return dense_->diagonal();
146 } else {
147 return sparse_->diagonal();
148 }
149 }
150
152 IndexPair get_indices() const { return indices_; }
153
155 void set_indices(Index f, Index s) { indices_ = {f, s}; }
156
159
160 const Matrix &get_dense() const {
162 return *dense_;
163 }
164
167 return *dense_;
168 }
169
170 const Sparse &get_sparse() const {
172 return *sparse_;
173 }
174
177 return *sparse_;
178 }
179
180 // Friend declarations.
181 friend void mult(MatrixView, ConstMatrixView, const Block &);
182 friend void mult(MatrixView, const Block &, ConstMatrixView);
183 friend void mult(VectorView, const Block &, ConstVectorView);
184
185 friend MatrixView &operator+=(MatrixView &, const Block &);
186
187 private:
190
192
193 std::shared_ptr<Matrix> dense_;
194 std::shared_ptr<Sparse> sparse_;
195};
196
197void mult(MatrixView, ConstMatrixView, const Block &);
198void mult(MatrixView, const Block &, ConstMatrixView);
199void mult(VectorView, const Block &, ConstVectorView);
200
202void add_inv(MatrixView A, const Block &);
203
204//------------------------------------------------------------------------------
205// Covariance Matrices
206//------------------------------------------------------------------------------
228 public:
229 CovarianceMatrix() = default;
234
235 ~CovarianceMatrix() = default;
236
237 explicit operator Matrix() const;
238 Matrix get_inverse() const;
239
240 Index nrows() const;
241 Index ncols() const;
242
244 Index ndiagblocks() const;
245
247 Index ninvdiagblocks() const;
248
250 Index nblocks() const;
251
259 bool has_block(Index i, Index j);
260
273 const Block *get_block(Index i = -1, Index j = -1);
274
280 std::vector<Block> &get_blocks() { return correlations_; };
281
287 std::vector<Block> &get_inverse_blocks() { return inverses_; };
288
300 bool has_diagonal_blocks(const ArrayOfArrayOfIndex &jis) const;
301
314 bool is_consistent(const ArrayOfArrayOfIndex &jis) const;
315
326 bool is_consistent(const Block &block) const;
327
333 void compute_inverse() const;
334
343 void add_correlation(Block c);
344
354
362 Vector diagonal() const;
363
372 Vector inverse_diagonal() const;
373
374 // Friend declarations.
375 friend void mult(MatrixView, ConstMatrixView, const CovarianceMatrix &);
376 friend void mult(MatrixView, const CovarianceMatrix &, ConstMatrixView);
377 friend void mult(VectorView, const CovarianceMatrix &, ConstVectorView);
378
381 friend void solve(VectorView, const CovarianceMatrix &, ConstVectorView);
382
384 friend void add_inv(MatrixView, const CovarianceMatrix &);
385
386 friend void xml_read_from_stream(istream &,
388 bifstream *,
389 const Verbosity &);
390 friend void xml_write_to_stream(ostream &,
391 const CovarianceMatrix &,
392 bofstream *,
393 const String &,
394 const Verbosity &);
395 friend std::ostream &operator<<(std::ostream &os, const CovarianceMatrix &v);
396
397 private:
398 void generate_blocks(std::vector<std::vector<const Block *>> &) const;
399 void invert_correlation_block(std::vector<Block> &inverses,
400 std::vector<const Block *> &blocks) const;
401 bool has_inverse(IndexPair indices) const;
402
403 std::vector<Block> correlations_;
404 mutable std::vector<Block> inverses_;
405};
406
410
414
416void add_inv(MatrixView, const CovarianceMatrix &);
417
418std::ostream &operator<<(std::ostream &os, const ConstVectorView &v);
419
420#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:1014
A constant view of a Vector.
Definition: matpackI.h:489
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:1125
The Matrix class.
Definition: matpackI.h:1225
The range class.
Definition: matpackI.h:165
constexpr Index get_extent() const ARTS_NOEXCEPT
Returns the extent of the range.
Definition: matpackI.h:334
The Sparse class.
Definition: matpackII.h:67
The VectorView class.
Definition: matpackI.h:626
The Vector class.
Definition: matpackI.h:876
Binary output file stream class.
Definition: bifstream.h:42
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:106
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:83
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
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:33
#define v
#define c