oem.h

Go to the documentation of this file.

 
#ifndef _ARTS_OEM_H_
#define _ARTS_OEM_H_
 
#include <type_traits>
 
#include "invlib/algebra.h"
#include "invlib/algebra/precision_matrix.h"
#include "invlib/algebra/solvers.h"
#include "invlib/interfaces/arts_wrapper.h"
#include "invlib/map.h"
#include "invlib/optimization.h"
 
//  Type Aliases
 
namespace oem {
 
using Vector = invlib::Vector<ArtsVector>;
using Matrix = invlib::Matrix<ArtsMatrix>;
using MatrixReference = invlib::Matrix<ArtsMatrixReference<::Matrix>>;
using CovarianceMatrix = invlib::Matrix<ArtsCovarianceMatrixWrapper>;
using Identity = invlib::MatrixIdentity<Matrix>;
 
// OEM Formulations
 
using invlib::Formulation;
 
template <typename ForwardModel>
using OEM_STANDARD = invlib::MAP<ForwardModel,
                                 Matrix,
                                 CovarianceMatrix,
                                 CovarianceMatrix,
                                 Vector,
                                 Formulation::STANDARD,
                                 invlib::Rodgers531>;
 
template <typename ForwardModel>
using OEM_NFORM = invlib::MAP<ForwardModel,
                              Matrix,
                              CovarianceMatrix,
                              CovarianceMatrix,
                              Vector,
                              Formulation::NFORM>;
 
template <typename ForwardModel>
using OEM_MFORM = invlib::MAP<ForwardModel,
                              Matrix,
                              CovarianceMatrix,
                              CovarianceMatrix,
                              Vector,
                              Formulation::MFORM>;
 
// Solvers
 
template <typename TransformationMatrixType,
          typename SolverType = invlib::Standard>
class NormalizingSolver : SolverType {
 public:
  template <typename... Params>
  NormalizingSolver(const TransformationMatrixType &trans,
                    bool apply,
                    Params... params)
      : SolverType(params...), apply_(apply), trans_(trans) {}
 
  template <typename MatrixType, typename VectorType>
  auto solve(const MatrixType &A, const VectorType &v) ->
      typename VectorType::ResultType {
    typename VectorType::ResultType w;
    if (apply_) {
      typename VectorType::ResultType vv = trans_ * v;
      auto &&ww = SolverType::solve(trans_ * A * trans_, vv);
      w = trans_ * ww;
    } else {
      w = SolverType::solve(A, v);
      VectorType u = v - A * w;
    }
    return w;
  }
 
 private:
  const bool apply_ = false;
  const TransformationMatrixType &trans_;
};
 
using Std = NormalizingSolver<Matrix, invlib::Standard>;
 
using CG = NormalizingSolver<Matrix, invlib::ConjugateGradient<>>;
 
using GN = invlib::GaussNewton<Numeric, Std>;
using GN_CG = invlib::GaussNewton<Numeric, CG>;
using LM = invlib::LevenbergMarquardt<Numeric, CovarianceMatrix, Std>;
using LM_CG = invlib::LevenbergMarquardt<Numeric, CovarianceMatrix, CG>;
 
//  Custom Log Class
 
template <typename T>
struct OptimizerLog;
 
template <typename RealType, typename DampingMatrix, typename Solver>
struct OptimizerLog<
    invlib::LevenbergMarquardt<RealType, DampingMatrix, Solver>> {
  static constexpr auto name = "Levenberg-Marquardt";
 
  static std::string header() {
    std::string out = "Gamma Factor";
    return out;
  }
 
  static std::string log(
      const invlib::LevenbergMarquardt<RealType, DampingMatrix, Solver> &g,
      Vector &gamma_history_,
      size_t i) {
    std::string lambda = std::to_string(g.get_lambda());
    std::string out(15 - std::min<size_t>(lambda.size(), 15), ' ');
    out += lambda;
    gamma_history_[i] = g.get_lambda();
    return out;
  }
};
 
template <typename RealType, typename Solver>
struct OptimizerLog<invlib::GaussNewton<RealType, Solver>> {
  static constexpr auto name = "Gauss-Newton";
 
  static std::string header() { return ""; }
 
  static std::string log(const invlib::GaussNewton<RealType, Solver> &,
                         Vector &,
                         size_t) {
    return "";
  }
};
 
template <invlib::LogType type>
class ArtsLog {
 public:
  ArtsLog(unsigned int v, ::Vector &g, bool l = false)
      : verbosity_(v), gamma_history_(g), linear_(l), finalized_(false) {}
 
  ~ArtsLog() {
    if ((verbosity_ >= 1) && (!finalized_)) {
      std::cout << invlib::separator() << std::endl << std::endl;
      std::cout << "Error during OEM computation." << std::endl;
      std::cout << std::endl;
      std::cout << invlib::center("----") << std::endl;
      std::cout << std::endl;
    }
  }
 
  template <typename... Params>
  void init(Params &... params) {
    if (verbosity_ >= 1) {
      std::tuple<Params &...> tuple(params...);
 
      auto &y = std::get<4>(tuple);
      scaling_factor_ = 1.0 / static_cast<Numeric>(y.nelem());
      std::cout << std::endl;
      std::cout << invlib::center("MAP Computation") << std::endl;
 
      // Print formulation.
      int formulation = static_cast<int>(std::get<6>(tuple));
      switch (formulation) {
        case 0:
          std::cout << "Formulation: Standard" << std::endl;
          break;
        case 1:
          std::cout << "Formulation: N-Form" << std::endl;
          break;
 
        case 2:
          std::cout << "Formulation: M-Form" << std::endl;
          break;
      }
 
      // Print optimization method.
      using OptimizationType = typename std::decay<
          typename std::tuple_element<5, decltype(tuple)>::type>::type;
      std::cout << "Method:      "
                << invlib::OptimizerLog<OptimizationType>::name;
      std::cout << std::endl;
 
      std::cout << std::endl;
      std::cout << std::setw(5) << "Step" << std::setw(15) << "Total Cost";
      std::cout << std::setw(15) << "x-Cost" << std::setw(15) << "y-Cost";
      std::cout << std::setw(15) << "Conv. Crit.";
      std::cout << std::setw(15) << OptimizerLog<OptimizationType>::header();
      std::cout << std::endl << invlib::separator() << std::endl;
    }
  }
 
  template <typename... Params>
  void step(const Params &... params) {
    if (verbosity_ >= 1) {
      std::tuple<const Params &...> tuple(params...);
      using OptimizationType = typename std::decay<
          typename std::tuple_element<5, decltype(tuple)>::type>::type;
 
      auto step_number = std::get<0>(tuple);
      std::cout << std::setw(5) << step_number;
      if (step_number == 0) {
        start_cost_ = std::get<1>(tuple);
      }
      std::cout << std::setw(15) << scaling_factor_ * std::get<1>(tuple);
      std::cout << std::setw(15) << scaling_factor_ * std::get<2>(tuple);
      std::cout << std::setw(15) << scaling_factor_ * std::get<3>(tuple);
 
      if (std::isnan(std::get<4>(tuple))) {
        std::cout << std::setw(15) << " ";
      } else {
        std::cout << std::setw(15) << std::get<4>(tuple);
      }
      std::cout << OptimizerLog<OptimizationType>::log(
          std::get<5>(tuple), gamma_history_, std::get<0>(tuple));
      std::cout << std::endl;
    }
  }
 
  template <typename... Params>
  void finalize(const Params &... params) {
    if (verbosity_ >= 1) {
      std::cout << invlib::separator() << std::endl;
 
      std::tuple<const Params &...> tuple(params...);
      std::cout << std::endl;
 
      std::cout << "Total number of steps:            ";
      std::cout << std::get<1>(tuple) << std::endl;
      std::cout << "Final scaled cost function value: ";
      std::cout << std::get<2>(tuple) * scaling_factor_ << std::endl;
 
      bool converged = std::get<0>(tuple);
      if (converged) {
        std::cout << "OEM computation converged." << std::endl;
      } else if (linear_) {
        std::cout << "Linear OEM computation finished." << std::endl;
      } else {
        std::cout << "OEM computation DID NOT converge!" << std::endl;
      }
    }
 
    finalized_ = true;
  }
 
  template <typename... Params>
  void time(const Params &... params) {
    if (verbosity_ >= 1) {
      std::tuple<const Params &...> tuple(params...);
      std::cout << std::endl;
      std::cout << "Elapsed Time for Retrieval:                       ";
      std::cout << std::get<0>(tuple) << std::endl;
      std::cout << "Time in inversion_iterate Agenda (No Jacobian):   ";
      std::cout << std::get<1>(tuple) << std::endl;
      std::cout << "Time in inversion_iterate Agenda (With Jacobian): ";
      std::cout << std::get<2>(tuple) << std::endl;
 
      std::cout << std::endl;
      std::cout << invlib::center("----") << std::endl;
      std::cout << std::endl;
    }
  }
 
 private:
  int verbosity_;
  Vector gamma_history_;
  Numeric scaling_factor_ = 0.0;
  Numeric start_cost_ = 0.0;
  bool linear_ = false;
  bool finalized_ = false;
};
 
//  Exception Handling
 
template <typename E>
std::vector<std::string> handle_nested_exception(const E &e, int level = 0) {
  const std::exception *re;
  std::vector<std::string> errors{};
 
  re = dynamic_cast<const std::exception *>(&e);
  if (re) {
    std::string s{};
 
    // If invlib level, extend error description.
    if (level == 0) {
      s = "Run-time error in oem computation: ";
    }
 
    s += re->what();
    errors.push_back(s);
  }
 
  try {
    std::rethrow_if_nested(e);
  } catch (const std::exception &ne) {
    std::vector<std::string> sv(handle_nested_exception(ne, level + 1));
    errors.insert(errors.end(), sv.begin(), sv.end());
  } catch (...) {
  }
  return errors;
}
 
// Forward model interface
 
class AgendaWrapper {
 public:
  const unsigned int m = 0;
  const unsigned int n = 0;
 
  AgendaWrapper(Workspace *ws,
                unsigned int measurement_space_dimension,
                unsigned int state_space_dimension,
                ::Matrix &arts_jacobian,
                ::Vector &arts_y,
                const Agenda *inversion_iterate_agenda)
      : m(measurement_space_dimension),
        n(state_space_dimension),
        inversion_iterate_agenda_(inversion_iterate_agenda),
        iteration_counter_(0),
        jacobian_(arts_jacobian),
        reuse_jacobian_((arts_jacobian.nrows() != 0) &&
                        (arts_jacobian.ncols() != 0) && (arts_y.nelem() != 0)),
        ws_(ws),
        yi_(arts_y) {}
 
  ArtsVector get_measurement_vector() { return yi_; }
 
  AgendaWrapper(const AgendaWrapper &) = delete;
  AgendaWrapper(AgendaWrapper &&) = delete;
  AgendaWrapper &operator=(const AgendaWrapper &) = delete;
  AgendaWrapper &operator=(AgendaWrapper &&) = delete;
 
  MatrixReference Jacobian(const Vector &xi, Vector &yi) {
    if (!reuse_jacobian_) {
      inversion_iterate_agendaExecute(
          *ws_, yi_, jacobian_, xi, 1, iteration_counter_, *inversion_iterate_agenda_);
      yi = yi_;
      iteration_counter_ += 1;
    } else {
      reuse_jacobian_ = false;
      yi = yi_;
    }
    return jacobian_;
  }
 
  Vector evaluate(const Vector &xi) {
    if (!reuse_jacobian_) {
      Matrix dummy;
      inversion_iterate_agendaExecute(*ws_,
                                      yi_,
                                      dummy,
                                      xi,
                                      0,
                                      iteration_counter_,
                                      *inversion_iterate_agenda_);
    } else {
      reuse_jacobian_ = false;
    }
    return yi_;
  }
 
 private:
  const Agenda *inversion_iterate_agenda_;
  unsigned int iteration_counter_;
  MatrixReference jacobian_;
  bool reuse_jacobian_;
  Workspace *ws_;
  Vector yi_;
};
}  // namespace oem
 
 
void Tensor4Clip(Tensor4& x,
                 const Index& iq,
                 const Numeric& limit_low,
                 const Numeric& limit_high) {
  // Sizes
  const Index nq = x.nbooks();
 
  ARTS_USER_ERROR_IF (iq < -1, "Argument *iq* must be >= -1.");
  ARTS_USER_ERROR_IF (iq >= nq,
      "Argument *iq* is too high.\n"
      "You have selected index: ", iq, "\n"
      "but the number of quantities is only: ", nq, "\n"
      "(Note that zero-based indexing is used)\n")
 
  Index ifirst = 0, ilast = nq - 1;
  if (iq > -1) {
    ifirst = iq;
    ilast = iq;
  }
 
  if (!std::isinf(limit_low)) {
    for (Index i = ifirst; i <= ilast; i++) {
      for (Index p = 0; p < x.npages(); p++) {
        for (Index r = 0; r < x.nrows(); r++) {
          for (Index c = 0; c < x.ncols(); c++) {
            if (x(i, p, r, c) < limit_low) x(i, p, r, c) = limit_low;
          }
        }
      }
    }
  }
 
  if (!std::isinf(limit_high)) {
    for (Index i = ifirst; i <= ilast; i++) {
      for (Index p = 0; p < x.npages(); p++) {
        for (Index r = 0; r < x.nrows(); r++) {
          for (Index c = 0; c < x.ncols(); c++) {
            if (x(i, p, r, c) > limit_high) x(i, p, r, c) = limit_high;
          }
        }
      }
    }
  }
}
 
// OEM error checking
 
void OEM_checks(Workspace& ws,
                Vector& x,
                Vector& yf,
                Matrix& jacobian,
                const Agenda& inversion_iterate_agenda,
                const Vector& xa,
                const CovarianceMatrix& covmat_sx,
                const Vector& y,
                const CovarianceMatrix& covmat_se,
                const ArrayOfRetrievalQuantity& jacobian_quantities,
                const String& method,
                const Vector& x_norm,
                const Index& max_iter,
                const Numeric& stop_dx,
                const Vector& lm_ga_settings,
                const Index& clear_matrices,
                const Index& display_progress) {
  const Index nq = jacobian_quantities.nelem();
  const Index n = xa.nelem();
  const Index m = y.nelem();
 
  ARTS_USER_ERROR_IF ((x.nelem() != n) && (x.nelem() != 0),
                      "The length of *x* must be either the same as *xa* or 0.");
  ARTS_USER_ERROR_IF (covmat_sx.ncols() != covmat_sx.nrows(),
                      "*covmat_sx* must be a square matrix.");
  ARTS_USER_ERROR_IF (covmat_sx.ncols() != n,
                      "Inconsistency in size between *x* and *covmat_sx*.");
  ARTS_USER_ERROR_IF ((yf.nelem() != m) && (yf.nelem() != 0),
                      "The length of *yf* must be either the same as *y* or 0.");
  ARTS_USER_ERROR_IF (covmat_se.ncols() != covmat_se.nrows(),
                      "*covmat_se* must be a square matrix.");
  ARTS_USER_ERROR_IF (covmat_se.ncols() != m,
                      "Inconsistency in size between *y* and *covmat_se*.");
  ARTS_USER_ERROR_IF ((jacobian.nrows() != m) && (!jacobian.empty()),
                      "The number of rows of the jacobian must be either the number of elements in *y* or 0.");
  ARTS_USER_ERROR_IF ((jacobian.ncols() != n) && (!jacobian.empty()),
                      "The number of cols of the jacobian must be either the number of elements in *xa* or 0.");
 
  ArrayOfArrayOfIndex jacobian_indices;
  bool any_affine;
  jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities);
  ARTS_USER_ERROR_IF (jacobian_indices.nelem() != nq,
        "Different number of elements in *jacobian_quantities* "
        "and *jacobian_indices*.");
  ARTS_USER_ERROR_IF (nq && jacobian_indices[nq - 1][1] + 1 != n,
        "Size of *covmat_sx* do not agree with Jacobian "
        "information (*jacobian_indices*).");
 
  // Check GINs
  ARTS_USER_ERROR_IF (!(method == "li" || method == "gn" || method == "li_m" ||
                        method == "gn_m" || method == "ml" || method == "lm" ||
                        method == "li_cg" || method == "gn_cg" || method == "li_cg_m" ||
                        method == "gn_cg_m" || method == "lm_cg" || method == "ml_cg"),
        "Valid options for *method* are \"nl\", \"gn\" and "
        "\"ml\" or \"lm\".");
 
  ARTS_USER_ERROR_IF (!(x_norm.nelem() == 0 || x_norm.nelem() == n),
        "The vector *x_norm* must have length 0 or match "
        "*covmat_sx*.");
 
  ARTS_USER_ERROR_IF (x_norm.nelem() > 0 && min(x_norm) <= 0,
                      "All values in *x_norm* must be > 0.");
 
  ARTS_USER_ERROR_IF (max_iter <= 0,
                      "The argument *max_iter* must be > 0.");
 
  ARTS_USER_ERROR_IF (stop_dx <= 0,
                      "The argument *stop_dx* must be > 0.");
 
  if ((method == "ml") || (method == "lm") || (method == "lm_cg") ||
      (method == "ml_cg")) {
    ARTS_USER_ERROR_IF (lm_ga_settings.nelem() != 6,
          "When using \"ml\", *lm_ga_setings* must be a "
          "vector of length 6.");
    ARTS_USER_ERROR_IF (min(lm_ga_settings) < 0,
          "The vector *lm_ga_setings* can not contain any "
          "negative value.");
  }
 
  ARTS_USER_ERROR_IF (clear_matrices < 0 || clear_matrices > 1,
                      "Valid options for *clear_matrices* are 0 and 1.");
  ARTS_USER_ERROR_IF (display_progress < 0 || display_progress > 1,
                      "Valid options for *display_progress* are 0 and 1.");
 
  // If necessary compute yf and jacobian.
  if (x.nelem() == 0) {
    x = xa;
    inversion_iterate_agendaExecute(
        ws, yf, jacobian, xa, 1, 0, inversion_iterate_agenda);
  }
  if ((yf.nelem() == 0) || (jacobian.empty())) {
    inversion_iterate_agendaExecute(
        ws, yf, jacobian, x, 1, 0, inversion_iterate_agenda);
  }
}
 
#endif  // _ARTS_OEM_H_