covariance_matrix.h

Go to the documentation of this file.

/* Copyright (C) 2017
   Simon Pfreundschuh <simonpf@chalmers.se>

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */

#ifndef covariance_matrix_h
#define covariance_matrix_h

#include <memory>

#include "jacobian.h"
#include "matpackI.h"
#include "matpackII.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_(indices),
        matrix_type_(MatrixType::dense),
        dense_(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_(indices),
        matrix_type_(MatrixType::sparse),
        dense_(nullptr),
        sparse_(sparse) {
    // Nothing to do here.
  }

  Block(const Block &) = default;
  Block(Block &&) = default;
  Block &operator=(const Block &) = default;
  Block &operator=(Block &&) = default;

  ~Block() = default;

  Index nrows() const { return row_range_.get_extent(); }
  Index ncols() const { return column_range_.get_extent(); }

  Range get_row_range() const { return row_range_; }
  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_ = sparse; }
  void set_matrix(std::shared_ptr<Matrix> dense) { dense_ = dense; }

  Vector diagonal() const {
    if (dense_) {
      return dense_->diagonal();
    } else {
      return sparse_->diagonal();
    }
  }

  IndexPair get_indices() const { return indices_; }
  
  void set_indices(Index f, Index s) { indices_ = {f, s}; }
  
  MatrixType get_matrix_type() const { return matrix_type_; }

  const Matrix &get_dense() const {
    assert(dense_);
    return *dense_;
  }

  Matrix &get_dense() {
    assert(dense_);
    return *dense_;
  }

  const Sparse &get_sparse() const {
    assert(sparse_);
    return *sparse_;
  }

  Sparse &get_sparse() {
    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 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 &);

std::ostream &operator<<(std::ostream &os, const ConstVectorView &v);

#endif  // covariance_matrix_h