17#include "invlib/algebra.h"
18#include "invlib/algebra/precision_matrix.h"
19#include "invlib/algebra/solvers.h"
20#include "invlib/interfaces/arts_wrapper.h"
21#include "invlib/map.h"
22#include "invlib/optimization.h"
31using Vector = invlib::Vector<ArtsVector>;
33using Matrix = invlib::Matrix<ArtsMatrix>;
38using Identity = invlib::MatrixIdentity<Matrix>;
45using invlib::Formulation;
55template <
typename ForwardModel>
61 Formulation::STANDARD,
72template <
typename ForwardModel>
88template <
typename ForwardModel>
110template <
typename TransformationMatrixType,
111 typename SolverType = invlib::Standard>
114 template <
typename... Params>
118 : SolverType(params...),
apply_(apply),
trans_(trans) {}
130 template <
typename MatrixType,
typename VectorType>
131 auto solve(
const MatrixType &A,
const VectorType &
v) ->
132 typename VectorType::ResultType {
133 typename VectorType::ResultType
w;
135 typename VectorType::ResultType vv =
trans_ *
v;
139 w = SolverType::solve(A,
v);
140 VectorType
u =
v - A *
w;
168using GN = invlib::GaussNewton<Numeric, Std>;
170using GN_CG = invlib::GaussNewton<Numeric, CG>;
172using LM = invlib::LevenbergMarquardt<Numeric, CovarianceMatrix, Std>;
174using LM_CG = invlib::LevenbergMarquardt<Numeric, CovarianceMatrix, CG>;
194template <
typename RealType,
typename DampingMatrix,
typename Solver>
196 invlib::LevenbergMarquardt<RealType, DampingMatrix, Solver>> {
198 static constexpr auto name =
"Levenberg-Marquardt";
202 std::string out =
"Gamma Factor";
208 const invlib::LevenbergMarquardt<RealType, DampingMatrix, Solver> &g,
211 std::string lambda = std::to_string(g.get_lambda());
212 std::string out(15 - std::min<size_t>(lambda.size(), 15),
' ');
214 gamma_history_[i] = g.get_lambda();
220template <
typename RealType,
typename Solver>
223 static constexpr auto name =
"Gauss-Newton";
226 static std::string
header() {
return ""; }
228 static std::string
log(
const invlib::GaussNewton<RealType, Solver> &,
242template <invlib::LogType type>
258 std::cout << invlib::separator() << std::endl << std::endl;
259 std::cout <<
"Error during OEM computation." << std::endl;
260 std::cout << std::endl;
261 std::cout << invlib::center(
"----") << std::endl;
262 std::cout << std::endl;
271 template <
typename... Params>
272 void init(Params &... params) {
274 std::tuple<Params &...> tuple(params...);
276 auto &y = std::get<4>(tuple);
278 std::cout << std::endl;
279 std::cout << invlib::center(
"MAP Computation") << std::endl;
282 int formulation =
static_cast<int>(std::get<6>(tuple));
283 switch (formulation) {
285 std::cout <<
"Formulation: Standard" << std::endl;
288 std::cout <<
"Formulation: N-Form" << std::endl;
292 std::cout <<
"Formulation: M-Form" << std::endl;
297 using OptimizationType =
typename std::decay<
298 typename std::tuple_element<5,
decltype(tuple)>::type>::type;
299 std::cout <<
"Method: "
300 << invlib::OptimizerLog<OptimizationType>::name;
301 std::cout << std::endl;
303 std::cout << std::endl;
304 std::cout << std::setw(5) <<
"Step" << std::setw(15) <<
"Total Cost";
305 std::cout << std::setw(15) <<
"x-Cost" << std::setw(15) <<
"y-Cost";
306 std::cout << std::setw(15) <<
"Conv. Crit.";
308 std::cout << std::endl << invlib::separator() << std::endl;
317 template <
typename... Params>
318 void step(
const Params &... params) {
320 std::tuple<
const Params &...> tuple(params...);
321 using OptimizationType =
typename std::decay<
322 typename std::tuple_element<5,
decltype(tuple)>::type>::type;
324 auto step_number = std::get<0>(tuple);
325 std::cout << std::setw(5) << step_number;
326 if (step_number == 0) {
329 std::cout << std::setw(15) << scaling_factor_ * std::get<1>(tuple);
330 std::cout << std::setw(15) << scaling_factor_ * std::get<2>(tuple);
331 std::cout << std::setw(15) << scaling_factor_ * std::get<3>(tuple);
333 if (std::isnan(std::get<4>(tuple))) {
334 std::cout << std::setw(15) <<
" ";
336 std::cout << std::setw(15) << std::get<4>(tuple);
338 std::cout << OptimizerLog<OptimizationType>::log(
340 std::cout << std::endl;
349 template <
typename... Params>
352 std::cout << invlib::separator() << std::endl;
354 std::tuple<
const Params &...> tuple(params...);
355 std::cout << std::endl;
357 std::cout <<
"Total number of steps: ";
358 std::cout << std::get<1>(tuple) << std::endl;
359 std::cout <<
"Final scaled cost function value: ";
362 bool converged = std::get<0>(tuple);
364 std::cout <<
"OEM computation converged." << std::endl;
366 std::cout <<
"Linear OEM computation finished." << std::endl;
368 std::cout <<
"OEM computation DID NOT converge!" << std::endl;
376 template <
typename... Params>
377 void time(
const Params &... params) {
379 std::tuple<
const Params &...> tuple(params...);
380 std::cout << std::endl;
381 std::cout <<
"Elapsed Time for Retrieval: ";
382 std::cout << std::get<0>(tuple) << std::endl;
383 std::cout <<
"Time in inversion_iterate Agenda (No Jacobian): ";
384 std::cout << std::get<1>(tuple) << std::endl;
385 std::cout <<
"Time in inversion_iterate Agenda (With Jacobian): ";
386 std::cout << std::get<2>(tuple) << std::endl;
388 std::cout << std::endl;
389 std::cout << invlib::center(
"----") << std::endl;
390 std::cout << std::endl;
425 const std::exception *re;
426 std::vector<std::string> errors{};
428 re =
dynamic_cast<const std::exception *
>(&e);
434 s =
"Run-time error in oem computation: ";
442 std::rethrow_if_nested(e);
443 }
catch (
const std::exception &ne) {
445 errors.insert(errors.end(), sv.begin(), sv.end());
464 const unsigned int m = 0;
466 const unsigned int n = 0;
483 unsigned int measurement_space_dimension,
484 unsigned int state_space_dimension,
487 const Agenda *inversion_iterate_agenda)
488 :
m(measurement_space_dimension),
489 n(state_space_dimension),
494 (arts_jacobian.ncols() != 0) && (arts_y.nelem() != 0)),
583 const Numeric& limit_low,
584 const Numeric& limit_high) {
586 const Index nq = x.nbooks();
589 ARTS_USER_ERROR_IF (iq >= nq,
590 "Argument *iq* is too high.\n
"
591 "You have selected index:
", iq, "\n
"
592 "but the number of quantities is only:
", nq, "\n
"
593 "(Note that zero-based indexing is used)\n
")
595 Index ifirst = 0, ilast = nq - 1;
601 if (!std::isinf(limit_low)) {
602 for (Index i = ifirst; i <= ilast; i++) {
603 for (Index p = 0; p < x.npages(); p++) {
604 for (Index r = 0; r < x.nrows(); r++) {
605 for (Index c = 0; c < x.ncols(); c++) {
606 if (x(i, p, r, c) < limit_low) x(i, p, r, c) = limit_low;
613 if (!std::isinf(limit_high)) {
614 for (Index i = ifirst; i <= ilast; i++) {
615 for (Index p = 0; p < x.npages(); p++) {
616 for (Index r = 0; r < x.nrows(); r++) {
617 for (Index c = 0; c < x.ncols(); c++) {
618 if (x(i, p, r, c) > limit_high) x(i, p, r, c) = limit_high;
662void OEM_checks(Workspace& ws,
666 const Agenda& inversion_iterate_agenda,
668 const CovarianceMatrix& covmat_sx,
670 const CovarianceMatrix& covmat_se,
671 const ArrayOfRetrievalQuantity& jacobian_quantities,
672 const String& method,
673 const Vector& x_norm,
674 const Index& max_iter,
675 const Numeric& stop_dx,
676 const Vector& lm_ga_settings,
677 const Index& clear_matrices,
678 const Index& display_progress) {
679 const Index nq = jacobian_quantities.nelem();
680 const Index n = xa.nelem();
681 const Index m = y.nelem();
683 ARTS_USER_ERROR_IF ((x.nelem() != n) && (x.nelem() != 0),
684 "The length of *x* must be either the same as *xa* or 0.
");
685 ARTS_USER_ERROR_IF (covmat_sx.ncols() != covmat_sx.nrows(),
686 "*covmat_sx* must be
a square matrix.
");
687 ARTS_USER_ERROR_IF (covmat_sx.ncols() != n,
688 "Inconsistency in size between *x* and *covmat_sx*.
");
689 ARTS_USER_ERROR_IF ((yf.nelem() != m) && (yf.nelem() != 0),
690 "The length of *yf* must be either the same as *y* or 0.
");
691 ARTS_USER_ERROR_IF (covmat_se.ncols() != covmat_se.nrows(),
692 "*covmat_se* must be
a square matrix.
");
693 ARTS_USER_ERROR_IF (covmat_se.ncols() != m,
694 "Inconsistency in size between *y* and *covmat_se*.
");
695 ARTS_USER_ERROR_IF ((jacobian.nrows() != m) && (!jacobian.empty()),
696 "The number of rows of the jacobian must be either the number of elements in *y* or 0.
");
697 ARTS_USER_ERROR_IF ((jacobian.ncols() != n) && (!jacobian.empty()),
698 "The number of cols of the jacobian must be either the number of elements in *xa* or 0.
");
700 ArrayOfArrayOfIndex jacobian_indices;
702 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities);
703 ARTS_USER_ERROR_IF (jacobian_indices.nelem() != nq,
704 "Different number of elements in *jacobian_quantities*
"
705 "and *jacobian_indices*.
");
706 ARTS_USER_ERROR_IF (nq && jacobian_indices[nq - 1][1] + 1 != n,
707 "Size of *covmat_sx*
do not agree with
Jacobian "
708 "information (*jacobian_indices*).
");
711 ARTS_USER_ERROR_IF (!(method == "li
" || method == "gn
" || method == "li_m
" ||
712 method == "gn_m
" || method == "ml
" || method == "lm
" ||
713 method == "li_cg
" || method == "gn_cg
" || method == "li_cg_m
" ||
714 method == "gn_cg_m
" || method == "lm_cg
" || method == "ml_cg
"),
715 "Valid options
for *method* are \
"nl\", \"gn\" and "
716 "\"ml\" or \"lm\".");
719 "The vector *x_norm* must have length 0 or match "
723 "All values in *x_norm* must be > 0.");
726 "The argument *max_iter* must be > 0.");
729 "The argument *stop_dx* must be > 0.");
731 if ((method ==
"ml") || (method ==
"lm") || (method ==
"lm_cg") ||
732 (method ==
"ml_cg")) {
734 "When using \"ml\", *lm_ga_setings* must be a "
735 "vector of length 6.");
737 "The vector *lm_ga_setings* can not contain any "
742 "Valid options for *clear_matrices* are 0 and 1.");
744 "Valid options for *display_progress* are 0 and 1.");
747 if (x.nelem() == 0) {
750 ws, yf, jacobian, xa, 1, 0, inversion_iterate_agenda);
752 if ((yf.nelem() == 0) || (jacobian.empty())) {
754 ws, yf, jacobian, x, 1, 0, inversion_iterate_agenda);
base min(const Array< base > &x)
Min function.
void inversion_iterate_agendaExecute(Workspace &ws, Vector &yf, Matrix &jacobian, const Vector &x, const Index jacobian_do, const Index inversion_iteration_counter, const Agenda &input_agenda)
Interface to ARTS inversion_iterate_agenda.
Vector evaluate(const Vector &xi)
Evaluate the ARTS forward model.
const Agenda * inversion_iterate_agenda_
Pointer to the inversion_iterate_agenda of the workspace.
bool reuse_jacobian_
Flag whether to reuse Jacobian from previous calculation.
AgendaWrapper(AgendaWrapper &&)=delete
const unsigned int m
Dimension of the measurement space.
AgendaWrapper(const AgendaWrapper &)=delete
Workspace * ws_
Pointer to current ARTS workspace.
AgendaWrapper & operator=(const AgendaWrapper &)=delete
ArtsVector get_measurement_vector()
Return most recently simulated measurement vector.
Vector yi_
Cached simulation result.
MatrixReference Jacobian(const Vector &xi, Vector &yi)
Evaluate forward model and compute Jacobian.
const unsigned int n
Dimension of the state space.
AgendaWrapper(Workspace *ws, unsigned int measurement_space_dimension, unsigned int state_space_dimension, ::Matrix &arts_jacobian, ::Vector &arts_y, const Agenda *inversion_iterate_agenda)
Create inversion_iterate_agendaExecute wrapper.
unsigned int iteration_counter_
AgendaWrapper & operator=(AgendaWrapper &&)=delete
MatrixReference jacobian_
Reference to the jacobian WSV.
void step(const Params &... params)
Print step to command line.
void init(Params &... params)
Initialize log output.
Vector gamma_history_
Reference to ARTS vector holding the LM gamma history.
bool linear_
Flag indicating whether forward model is linear.
bool finalized_
Flag indicating whether output has been finalized.
ArtsLog(unsigned int v, ::Vector &g, bool l=false)
Create log.
void finalize(const Params &... params)
Finalize log output.
Numeric start_cost_
Start cost (not computed by invlib)
Numeric scaling_factor_
Scaling factor for the cost.
void time(const Params &... params)
Print timing information to command line.
~ArtsLog()
Finalizes log output if necessary.
int verbosity_
Verbosity level of logger.
const TransformationMatrixType & trans_
The transformation matrix.
NormalizingSolver(const TransformationMatrixType &trans, bool apply, Params... params)
auto solve(const MatrixType &A, const VectorType &v) -> typename VectorType::ResultType
Solve linear system.
const bool apply_
Whether or not to apply the transformation.
#define ARTS_USER_ERROR_IF(condition,...)
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::MFORM > OEM_MFORM
OEM m form.
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, CG > LM_CG
Levenberg-Marquardt (LM) optimization using normed CG solver.
invlib::MatrixIdentity< Matrix > Identity
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::STANDARD, invlib::Rodgers531 > OEM_STANDARD
OEM standard form.
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, Std > LM
Levenberg-Marquardt (LM) optimization using normed ARTS QR solver.
invlib::Matrix< ArtsCovarianceMatrixWrapper > CovarianceMatrix
invlib wrapper type for ARTS the ARTS covariance class.
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
invlib::GaussNewton< Numeric, CG > GN_CG
Gauss-Newton (GN) optimization using normed CG solver.
std::vector< std::string > handle_nested_exception(const E &e, int level=0)
Handle exception encountered within invlib.
invlib::GaussNewton< Numeric, Std > GN
OEM Gauss-Newton optimization using normed ARTS QR solver.
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::NFORM > OEM_NFORM
OEM n form.
invlib::Matrix< ArtsMatrixReference<::Matrix > > MatrixReference
invlib wrapper type for ARTS matrices to be passed by reference.
void Tensor4Clip(Tensor4 &x, const Index &iq, const Numeric &limit_low, const Numeric &limit_high)
Clip Tensor4.
static std::string header()
Name to append to header line.
static std::string log(const invlib::GaussNewton< RealType, Solver > &, Vector &, size_t)
static std::string log(const invlib::LevenbergMarquardt< RealType, DampingMatrix, Solver > &g, Vector &gamma_history_, size_t i)
Returns the string to append to the log of a single step.
static std::string header()
Name to append to header line.
Log customization for different optimization methods.