covariance_matrix.h

Go to the documentation of this file.

 
#ifndef covariance_matrix_h
#define covariance_matrix_h
 
#include <memory>
#include <utility>
 
#include "jacobian.h"
#include "matpack_data.h"
#include "matpack_math.h"
#include "matpack_sparse.h"
 
class CovarianceMatrix;
 
//------------------------------------------------------------------------------
// Type Aliases
//------------------------------------------------------------------------------
 
using IndexPair = std::pair<Index, Index>;
 
//------------------------------------------------------------------------------
// Block
//------------------------------------------------------------------------------
class Block {
 public:
  /*
    * Enum class representing the type of the matrix representing the block, i.e.
    * whether its of type Matrix (dense) or of type Sparse(sparse).
    */
  enum class MatrixType { dense, sparse };
 
  /*
     * Create a correlation block from given row_range, column_range and shared_ptr
     * to a dense matrix.
     *
     * @param row_range The range of element-rows covered by the block
     * @param column_range The range of element-column covered by the block
     * @param indices The pair of block-indices that identifies the block in the
     *        covariance matrix.
     * @param dense A shared pointer to a den matrix of type Matrix
     */
  Block(Range row_range,
        Range column_range,
        IndexPair indices,
        std::shared_ptr<Matrix> dense)
      : row_range_(row_range),
        column_range_(column_range),
        indices_(std::move(indices)),
        matrix_type_(MatrixType::dense),
        dense_(std::move(dense)),
        sparse_(nullptr) {
    // Nothing to do here.
  }
 
  /*
     * Create a correlation block from given row_range, column_range and shared_ptr
     * to a sparse matrix.
     *
     * @param row_range The range of element-rows covered by the block
     * @param column_range The range of element-column covered by the block
     * @param indices The pair of block-indices that identifies the block in the
     *        covariance matrix.
     * @param dense A shared pointer to a den matrix of type Sparse
     */
  Block(Range row_range,
        Range column_range,
        IndexPair indices,
        std::shared_ptr<Sparse> sparse)
      : row_range_(row_range),
        column_range_(column_range),
        indices_(std::move(indices)),
        matrix_type_(MatrixType::sparse),
        dense_(nullptr),
        sparse_(std::move(sparse)) {
    // Nothing to do here.
  }
 
  Block(const Block &) = default;
  Block(Block &&) = default;
  Block &operator=(const Block &) = default;
  Block &operator=(Block &&) = default;
 
  ~Block() = default;
 
  [[nodiscard]] Index nrows() const { return row_range_.extent; }
  [[nodiscard]] Index ncols() const { return column_range_.extent; }
 
  [[nodiscard]] Range get_row_range() const { return row_range_; }
  [[nodiscard]] Range get_column_range() const { return column_range_; }
 
  Range& get_row_range() { return row_range_; }
  Range& get_column_range() { return column_range_; }
 
  void set_matrix(std::shared_ptr<Sparse> sparse) { sparse_ = std::move(sparse); dense_ = nullptr; matrix_type_=MatrixType::sparse; }
  void set_matrix(std::shared_ptr<Matrix> dense) { dense_ = std::move(dense); sparse_ = nullptr; matrix_type_=MatrixType::dense; }
 
  [[nodiscard]] Vector diagonal() const {
    if (dense_) {
      return ::diagonal(*dense_);
    }        
    return sparse_->diagonal();
   
  }
 
  [[nodiscard]] IndexPair get_indices() const { return indices_; }
 
  void set_indices(Index f, Index s) { indices_ = {f, s}; }
 
  [[nodiscard]] MatrixType get_matrix_type() const { return matrix_type_; }
 
  [[nodiscard]] const Matrix &get_dense() const {
    ARTS_ASSERT(dense_);
    return *dense_;
  }
 
  Matrix &get_dense() {
    ARTS_ASSERT(dense_);
    return *dense_;
  }
 
  [[nodiscard]] const Sparse &get_sparse() const {
    ARTS_ASSERT(sparse_);
    return *sparse_;
  }
 
  Sparse &get_sparse() {
    ARTS_ASSERT(sparse_);
    return *sparse_;
  }
 
  // Friend declarations.
  friend void mult(MatrixView, ConstMatrixView, const Block &);
  friend void mult(MatrixView, const Block &, ConstMatrixView);
  friend void mult(VectorView, const Block &, ConstVectorView);
 
  friend MatrixView operator+=(MatrixView, const Block &);
 
 private:
  Range row_range_, column_range_;
  IndexPair indices_;
 
  MatrixType matrix_type_;
 
  std::shared_ptr<Matrix> dense_;
  std::shared_ptr<Sparse> sparse_;
};
 
void mult(MatrixView, ConstMatrixView, const Block &);
void mult(MatrixView, const Block &, ConstMatrixView);
void mult(VectorView, const Block &, ConstVectorView);
 
MatrixView operator+=(MatrixView, const Block &);
void add_inv(MatrixView A, const Block &);
 
//------------------------------------------------------------------------------
// Covariance Matrices
//------------------------------------------------------------------------------
class CovarianceMatrix {
 public:
  CovarianceMatrix() = default;
  CovarianceMatrix(const CovarianceMatrix &) = default;
  CovarianceMatrix(CovarianceMatrix &&) = default;
  CovarianceMatrix &operator=(const CovarianceMatrix &) = default;
  CovarianceMatrix &operator=(CovarianceMatrix &&) = default;
 
  ~CovarianceMatrix() = default;
 
  explicit operator Matrix() const;
  Matrix get_inverse() const;
 
  Index nrows() const;
  Index ncols() const;
 
  Index ndiagblocks() const;
 
  Index ninvdiagblocks() const;
 
  Index nblocks() const;
 
  bool has_block(Index i, Index j);
 
  const Block *get_block(Index i = -1, Index j = -1);
 
  std::vector<Block> &get_blocks() { return correlations_; };
 
  std::vector<Block> &get_inverse_blocks() { return inverses_; };
 
  bool has_diagonal_blocks(const ArrayOfArrayOfIndex &jis) const;
 
  bool is_consistent(const ArrayOfArrayOfIndex &jis) const;
 
  bool is_consistent(const Block &block) const;
 
  void compute_inverse() const;
 
  void add_correlation(Block c);
 
  void add_correlation_inverse(Block c);
 
  Vector diagonal() const;
 
  Vector inverse_diagonal() const;
 
  // Friend declarations.
  friend void mult(MatrixView, ConstMatrixView, const CovarianceMatrix &);
  friend void mult(MatrixView, const CovarianceMatrix &, ConstMatrixView);
  friend void mult(VectorView, const CovarianceMatrix &, ConstVectorView);
 
  friend void mult_inv(MatrixView, ConstMatrixView, const CovarianceMatrix &);
  friend void mult_inv(MatrixView, const CovarianceMatrix &, ConstMatrixView);
  friend void solve(VectorView, const CovarianceMatrix &, ConstVectorView);
 
  friend MatrixView operator+=(MatrixView, const CovarianceMatrix &);
  friend void add_inv(MatrixView, const CovarianceMatrix &);
 
  friend void xml_read_from_stream(istream &,
                                   CovarianceMatrix &,
                                   bifstream *,
                                   const Verbosity &);
  friend void xml_write_to_stream(ostream &,
                                  const CovarianceMatrix &,
                                  bofstream *,
                                  const String &,
                                  const Verbosity &);
  friend std::ostream &operator<<(std::ostream &os, const CovarianceMatrix &v);
 
 private:
  void generate_blocks(std::vector<std::vector<const Block *>> &) const;
  void invert_correlation_block(std::vector<Block> &inverses,
                                std::vector<const Block *> &blocks) const;
  bool has_inverse(IndexPair indices) const;
 
  std::vector<Block> correlations_;
  mutable std::vector<Block> inverses_;
};
 
void mult(MatrixView, ConstMatrixView, const CovarianceMatrix &);
void mult(MatrixView, const CovarianceMatrix &, ConstMatrixView);
void mult(VectorView, const CovarianceMatrix &, ConstVectorView);
 
void mult_inv(MatrixView, ConstMatrixView, const CovarianceMatrix &);
void mult_inv(MatrixView, const CovarianceMatrix &, ConstMatrixView);
void solve(VectorView, const CovarianceMatrix &, ConstVectorView);
 
MatrixView operator+=(MatrixView, const CovarianceMatrix &);
void add_inv(MatrixView, const CovarianceMatrix &);
 
#endif  // covariance_matrix_h