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"
#include "invlib/profiling/timer.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, 0, *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();

  if (iq < -1) throw runtime_error("Argument *iq* must be >= -1.");
  if (iq >= nq) {
    ostringstream os;
    os << "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";
    throw runtime_error(os.str());
  }

  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();

  if ((x.nelem() != n) && (x.nelem() != 0))
    throw runtime_error(
        "The length of *x* must be either the same as *xa* or 0.");
  if (covmat_sx.ncols() != covmat_sx.nrows())
    throw runtime_error("*covmat_sx* must be a square matrix.");
  if (covmat_sx.ncols() != n)
    throw runtime_error("Inconsistency in size between *x* and *covmat_sx*.");
  if ((yf.nelem() != m) && (yf.nelem() != 0))
    throw runtime_error(
        "The length of *yf* must be either the same as *y* or 0.");
  if (covmat_se.ncols() != covmat_se.nrows())
    throw runtime_error("*covmat_se* must be a square matrix.");
  if (covmat_se.ncols() != m)
    throw runtime_error("Inconsistency in size between *y* and *covmat_se*.");
  if ((jacobian.nrows() != m) && (!jacobian.empty()))
    throw runtime_error(
        "The number of rows of the jacobian must be either the number of elements in *y* or 0.");
  if ((jacobian.ncols() != n) && (!jacobian.empty()))
    throw runtime_error(
        "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);
  if (jacobian_indices.nelem() != nq)
    throw runtime_error(
        "Different number of elements in *jacobian_quantities* "
        "and *jacobian_indices*.");
  if (nq && jacobian_indices[nq - 1][1] + 1 != n)
    throw runtime_error(
        "Size of *covmat_sx* do not agree with Jacobian "
        "information (*jacobian_indices*).");

  // Check GINs
  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")) {
    throw runtime_error(
        "Valid options for *method* are \"nl\", \"gn\" and "
        "\"ml\" or \"lm\".");
  }

  if (!(x_norm.nelem() == 0 || x_norm.nelem() == n)) {
    throw runtime_error(
        "The vector *x_norm* must have length 0 or match "
        "*covmat_sx*.");
  }

  if (x_norm.nelem() > 0 && min(x_norm) <= 0) {
    throw runtime_error("All values in *x_norm* must be > 0.");
  }

  if (max_iter <= 0) {
    throw runtime_error("The argument *max_iter* must be > 0.");
  }

  if (stop_dx <= 0) {
    throw runtime_error("The argument *stop_dx* must be > 0.");
  }

  if ((method == "ml") || (method == "lm") || (method == "lm_cg") ||
      (method == "ml_cg")) {
    if (lm_ga_settings.nelem() != 6) {
      throw runtime_error(
          "When using \"ml\", *lm_ga_setings* must be a "
          "vector of length 6.");
    }
    if (min(lm_ga_settings) < 0) {
      throw runtime_error(
          "The vector *lm_ga_setings* can not contain any "
          "negative value.");
    }
  }

  if (clear_matrices < 0 || clear_matrices > 1)
    throw runtime_error("Valid options for *clear_matrices* are 0 and 1.");
  if (display_progress < 0 || display_progress > 1)
    throw runtime_error("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_