ARTS 2.5.0 (git: 9ee3ac6c)
oem.h
Go to the documentation of this file.
1
12#ifndef _ARTS_OEM_H_
13#define _ARTS_OEM_H_
14
15#include <type_traits>
16
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"
23
25// Type Aliases
27
28namespace oem {
29
31using Vector = invlib::Vector<ArtsVector>;
33using Matrix = invlib::Matrix<ArtsMatrix>;
35using MatrixReference = invlib::Matrix<ArtsMatrixReference<::Matrix>>;
37using CovarianceMatrix = invlib::Matrix<ArtsCovarianceMatrixWrapper>;
38using Identity = invlib::MatrixIdentity<Matrix>;
39
41// OEM Formulations
43
45using invlib::Formulation;
46
55template <typename ForwardModel>
56using OEM_STANDARD = invlib::MAP<ForwardModel,
57 Matrix,
60 Vector,
61 Formulation::STANDARD,
62 invlib::Rodgers531>;
63
72template <typename ForwardModel>
73using OEM_NFORM = invlib::MAP<ForwardModel,
74 Matrix,
77 Vector,
78 Formulation::NFORM>;
79
88template <typename ForwardModel>
89using OEM_MFORM = invlib::MAP<ForwardModel,
90 Matrix,
93 Vector,
94 Formulation::MFORM>;
95
97// Solvers
99
110template <typename TransformationMatrixType,
111 typename SolverType = invlib::Standard>
112class NormalizingSolver : SolverType {
113 public:
114 template <typename... Params>
115 NormalizingSolver(const TransformationMatrixType &trans,
116 bool apply,
117 Params... params)
118 : SolverType(params...), apply_(apply), trans_(trans) {}
119
130 template <typename MatrixType, typename VectorType>
131 auto solve(const MatrixType &A, const VectorType &v) ->
132 typename VectorType::ResultType {
133 typename VectorType::ResultType w;
134 if (apply_) {
135 typename VectorType::ResultType vv = trans_ * v;
136 auto &&ww = SolverType::solve(trans_ * A * trans_, vv);
137 w = trans_ * ww;
138 } else {
139 w = SolverType::solve(A, v);
140 VectorType u = v - A * w;
141 }
142 return w;
143 }
144
145 private:
147 const bool apply_ = false;
149 const TransformationMatrixType &trans_;
150};
151
158
166
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>;
175
177// Custom Log Class
179
186template <typename T>
188
194template <typename RealType, typename DampingMatrix, typename Solver>
196 invlib::LevenbergMarquardt<RealType, DampingMatrix, Solver>> {
198 static constexpr auto name = "Levenberg-Marquardt";
199
202 std::string out = "Gamma Factor";
203 return out;
204 }
205
207 static std::string log(
208 const invlib::LevenbergMarquardt<RealType, DampingMatrix, Solver> &g,
209 Vector &gamma_history_,
210 size_t i) {
211 std::string lambda = std::to_string(g.get_lambda());
212 std::string out(15 - std::min<size_t>(lambda.size(), 15), ' ');
213 out += lambda;
214 gamma_history_[i] = g.get_lambda();
215 return out;
216 }
217};
218
220template <typename RealType, typename Solver>
221struct OptimizerLog<invlib::GaussNewton<RealType, Solver>> {
223 static constexpr auto name = "Gauss-Newton";
224
226 static std::string header() { return ""; }
227
228 static std::string log(const invlib::GaussNewton<RealType, Solver> &,
229 Vector &,
230 size_t) {
231 return "";
232 }
233};
234
242template <invlib::LogType type>
243class ArtsLog {
244 public:
252 ArtsLog(unsigned int v, ::Vector &g, bool l = false)
253 : verbosity_(v), gamma_history_(g), linear_(l), finalized_(false) {}
254
257 if ((verbosity_ >= 1) && (!finalized_)) {
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;
263 }
264 }
265
271 template <typename... Params>
272 void init(Params &... params) {
273 if (verbosity_ >= 1) {
274 std::tuple<Params &...> tuple(params...);
275
276 auto &y = std::get<4>(tuple);
277 scaling_factor_ = 1.0 / static_cast<Numeric>(y.nelem());
278 std::cout << std::endl;
279 std::cout << invlib::center("MAP Computation") << std::endl;
280
281 // Print formulation.
282 int formulation = static_cast<int>(std::get<6>(tuple));
283 switch (formulation) {
284 case 0:
285 std::cout << "Formulation: Standard" << std::endl;
286 break;
287 case 1:
288 std::cout << "Formulation: N-Form" << std::endl;
289 break;
290
291 case 2:
292 std::cout << "Formulation: M-Form" << std::endl;
293 break;
294 }
295
296 // Print optimization method.
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;
302
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.";
307 std::cout << std::setw(15) << OptimizerLog<OptimizationType>::header();
308 std::cout << std::endl << invlib::separator() << std::endl;
309 }
310 }
311
317 template <typename... Params>
318 void step(const Params &... params) {
319 if (verbosity_ >= 1) {
320 std::tuple<const Params &...> tuple(params...);
321 using OptimizationType = typename std::decay<
322 typename std::tuple_element<5, decltype(tuple)>::type>::type;
323
324 auto step_number = std::get<0>(tuple);
325 std::cout << std::setw(5) << step_number;
326 if (step_number == 0) {
327 start_cost_ = std::get<1>(tuple);
328 }
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);
332
333 if (std::isnan(std::get<4>(tuple))) {
334 std::cout << std::setw(15) << " ";
335 } else {
336 std::cout << std::setw(15) << std::get<4>(tuple);
337 }
338 std::cout << OptimizerLog<OptimizationType>::log(
339 std::get<5>(tuple), gamma_history_, std::get<0>(tuple));
340 std::cout << std::endl;
341 }
342 }
343
349 template <typename... Params>
350 void finalize(const Params &... params) {
351 if (verbosity_ >= 1) {
352 std::cout << invlib::separator() << std::endl;
353
354 std::tuple<const Params &...> tuple(params...);
355 std::cout << std::endl;
356
357 std::cout << "Total number of steps: ";
358 std::cout << std::get<1>(tuple) << std::endl;
359 std::cout << "Final scaled cost function value: ";
360 std::cout << std::get<2>(tuple) * scaling_factor_ << std::endl;
361
362 bool converged = std::get<0>(tuple);
363 if (converged) {
364 std::cout << "OEM computation converged." << std::endl;
365 } else if (linear_) {
366 std::cout << "Linear OEM computation finished." << std::endl;
367 } else {
368 std::cout << "OEM computation DID NOT converge!" << std::endl;
369 }
370 }
371
372 finalized_ = true;
373 }
374
376 template <typename... Params>
377 void time(const Params &... params) {
378 if (verbosity_ >= 1) {
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;
387
388 std::cout << std::endl;
389 std::cout << invlib::center("----") << std::endl;
390 std::cout << std::endl;
391 }
392 }
393
394 private:
404 bool linear_ = false;
406 bool finalized_ = false;
407};
408
410// Exception Handling
412
423template <typename E>
424std::vector<std::string> handle_nested_exception(const E &e, int level = 0) {
425 const std::exception *re;
426 std::vector<std::string> errors{};
427
428 re = dynamic_cast<const std::exception *>(&e);
429 if (re) {
430 std::string s{};
431
432 // If invlib level, extend error description.
433 if (level == 0) {
434 s = "Run-time error in oem computation: ";
435 }
436
437 s += re->what();
438 errors.push_back(s);
439 }
440
441 try {
442 std::rethrow_if_nested(e);
443 } catch (const std::exception &ne) {
444 std::vector<std::string> sv(handle_nested_exception(ne, level + 1));
445 errors.insert(errors.end(), sv.begin(), sv.end());
446 } catch (...) {
447 }
448 return errors;
449}
450
452// Forward model interface
454
462 public:
464 const unsigned int m = 0;
466 const unsigned int n = 0;
467
483 unsigned int measurement_space_dimension,
484 unsigned int state_space_dimension,
485 ::Matrix &arts_jacobian,
486 ::Vector &arts_y,
487 const Agenda *inversion_iterate_agenda)
488 : m(measurement_space_dimension),
489 n(state_space_dimension),
490 inversion_iterate_agenda_(inversion_iterate_agenda),
492 jacobian_(arts_jacobian),
493 reuse_jacobian_((arts_jacobian.nrows() != 0) &&
494 (arts_jacobian.ncols() != 0) && (arts_y.nelem() != 0)),
495 ws_(ws),
496 yi_(arts_y) {}
497
502 ArtsVector get_measurement_vector() { return yi_; }
503
504 AgendaWrapper(const AgendaWrapper &) = delete;
508
521 if (!reuse_jacobian_) {
524 yi = yi_;
526 } else {
527 reuse_jacobian_ = false;
528 yi = yi_;
529 }
530 return jacobian_;
531 }
532
542 Vector evaluate(const Vector &xi) {
543 if (!reuse_jacobian_) {
544 Matrix dummy;
546 yi_,
547 dummy,
548 xi,
549 0,
552 } else {
553 reuse_jacobian_ = false;
554 }
555 return yi_;
556 }
557
558 private:
561 unsigned int iteration_counter_;
570};
571} // namespace oem
572
573
582 const Index& iq,
583 const Numeric& limit_low,
584 const Numeric& limit_high) {
585 // Sizes
586 const Index nq = x.nbooks();
587
588 ARTS_USER_ERROR_IF (iq < -1, "Argument *iq* must be >= -1.");
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")
594
595 Index ifirst = 0, ilast = nq - 1;
596 if (iq > -1) {
597 ifirst = iq;
598 ilast = iq;
599 }
600
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;
607 }
608 }
609 }
610 }
611 }
612
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;
619 }
620 }
621 }
622 }
623 }
624}
625
627// OEM error checking
629
662void OEM_checks(Workspace& ws,
663 Vector& x,
664 Vector& yf,
665 Matrix& jacobian,
666 const Agenda& inversion_iterate_agenda,
667 const Vector& xa,
668 const CovarianceMatrix& covmat_sx,
669 const Vector& y,
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();
682
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.");
699
700 ArrayOfArrayOfIndex jacobian_indices;
701 bool any_affine;
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*).");
709
710 // Check GINs
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\".");
717
718 ARTS_USER_ERROR_IF (!(x_norm.nelem() == 0 || x_norm.nelem() == n),
719 "The vector *x_norm* must have length 0 or match "
720 "*covmat_sx*.");
721
722 ARTS_USER_ERROR_IF (x_norm.nelem() > 0 && min(x_norm) <= 0,
723 "All values in *x_norm* must be > 0.");
724
725 ARTS_USER_ERROR_IF (max_iter <= 0,
726 "The argument *max_iter* must be > 0.");
727
728 ARTS_USER_ERROR_IF (stop_dx <= 0,
729 "The argument *stop_dx* must be > 0.");
730
731 if ((method == "ml") || (method == "lm") || (method == "lm_cg") ||
732 (method == "ml_cg")) {
733 ARTS_USER_ERROR_IF (lm_ga_settings.nelem() != 6,
734 "When using \"ml\", *lm_ga_setings* must be a "
735 "vector of length 6.");
736 ARTS_USER_ERROR_IF (min(lm_ga_settings) < 0,
737 "The vector *lm_ga_setings* can not contain any "
738 "negative value.");
739 }
740
741 ARTS_USER_ERROR_IF (clear_matrices < 0 || clear_matrices > 1,
742 "Valid options for *clear_matrices* are 0 and 1.");
743 ARTS_USER_ERROR_IF (display_progress < 0 || display_progress > 1,
744 "Valid options for *display_progress* are 0 and 1.");
745
746 // If necessary compute yf and jacobian.
747 if (x.nelem() == 0) {
748 x = xa;
750 ws, yf, jacobian, xa, 1, 0, inversion_iterate_agenda);
751 }
752 if ((yf.nelem() == 0) || (jacobian.empty())) {
754 ws, yf, jacobian, x, 1, 0, inversion_iterate_agenda);
755 }
756}
757
758#endif // _ARTS_OEM_H_
Index nrows
Index ncols
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)
Definition: auto_md.cc:23250
The Agenda class.
Definition: agenda_class.h:44
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:55
The Matrix class.
Definition: matpackI.h:1225
The Tensor4 class.
Definition: matpackIV.h:421
The Vector class.
Definition: matpackI.h:876
Workspace class.
Definition: workspace_ng.h:40
Interface to ARTS inversion_iterate_agenda.
Definition: oem.h:461
Vector evaluate(const Vector &xi)
Evaluate the ARTS forward model.
Definition: oem.h:542
const Agenda * inversion_iterate_agenda_
Pointer to the inversion_iterate_agenda of the workspace.
Definition: oem.h:560
bool reuse_jacobian_
Flag whether to reuse Jacobian from previous calculation.
Definition: oem.h:565
AgendaWrapper(AgendaWrapper &&)=delete
const unsigned int m
Dimension of the measurement space.
Definition: oem.h:464
AgendaWrapper(const AgendaWrapper &)=delete
Workspace * ws_
Pointer to current ARTS workspace.
Definition: oem.h:567
AgendaWrapper & operator=(const AgendaWrapper &)=delete
ArtsVector get_measurement_vector()
Return most recently simulated measurement vector.
Definition: oem.h:502
Vector yi_
Cached simulation result.
Definition: oem.h:569
MatrixReference Jacobian(const Vector &xi, Vector &yi)
Evaluate forward model and compute Jacobian.
Definition: oem.h:520
const unsigned int n
Dimension of the state space.
Definition: oem.h:466
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.
Definition: oem.h:482
unsigned int iteration_counter_
Definition: oem.h:561
AgendaWrapper & operator=(AgendaWrapper &&)=delete
MatrixReference jacobian_
Reference to the jacobian WSV.
Definition: oem.h:563
OEM log output.
Definition: oem.h:243
void step(const Params &... params)
Print step to command line.
Definition: oem.h:318
void init(Params &... params)
Initialize log output.
Definition: oem.h:272
Vector gamma_history_
Reference to ARTS vector holding the LM gamma history.
Definition: oem.h:398
bool linear_
Flag indicating whether forward model is linear.
Definition: oem.h:404
bool finalized_
Flag indicating whether output has been finalized.
Definition: oem.h:406
ArtsLog(unsigned int v, ::Vector &g, bool l=false)
Create log.
Definition: oem.h:252
void finalize(const Params &... params)
Finalize log output.
Definition: oem.h:350
Numeric start_cost_
Start cost (not computed by invlib)
Definition: oem.h:402
Numeric scaling_factor_
Scaling factor for the cost.
Definition: oem.h:400
void time(const Params &... params)
Print timing information to command line.
Definition: oem.h:377
~ArtsLog()
Finalizes log output if necessary.
Definition: oem.h:256
int verbosity_
Verbosity level of logger.
Definition: oem.h:396
Normalizing solver.
Definition: oem.h:112
const TransformationMatrixType & trans_
The transformation matrix.
Definition: oem.h:149
NormalizingSolver(const TransformationMatrixType &trans, bool apply, Params... params)
Definition: oem.h:115
auto solve(const MatrixType &A, const VectorType &v) -> typename VectorType::ResultType
Solve linear system.
Definition: oem.h:131
const bool apply_
Whether or not to apply the transformation.
Definition: oem.h:147
void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v)
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
#define min(a, b)
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixType
Definition: matpackI.h:114
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Index nelem(const Lines &l)
Number of lines.
char Type type
Matrix matrix(const std::array< Numeric, 196 > data)
Create the square matrix from the static data.
Definition: igrf13.cc:206
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
constexpr bool isnan(double d) noexcept
Definition: nonstd.h:39
Definition: oem.h:28
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::MFORM > OEM_MFORM
OEM m form.
Definition: oem.h:94
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, CG > LM_CG
Levenberg-Marquardt (LM) optimization using normed CG solver.
Definition: oem.h:174
invlib::MatrixIdentity< Matrix > Identity
Definition: oem.h:38
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::STANDARD, invlib::Rodgers531 > OEM_STANDARD
OEM standard form.
Definition: oem.h:62
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, Std > LM
Levenberg-Marquardt (LM) optimization using normed ARTS QR solver.
Definition: oem.h:172
invlib::Matrix< ArtsCovarianceMatrixWrapper > CovarianceMatrix
invlib wrapper type for ARTS the ARTS covariance class.
Definition: oem.h:37
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:33
invlib::GaussNewton< Numeric, CG > GN_CG
Gauss-Newton (GN) optimization using normed CG solver.
Definition: oem.h:170
std::vector< std::string > handle_nested_exception(const E &e, int level=0)
Handle exception encountered within invlib.
Definition: oem.h:424
invlib::GaussNewton< Numeric, Std > GN
OEM Gauss-Newton optimization using normed ARTS QR solver.
Definition: oem.h:168
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::NFORM > OEM_NFORM
OEM n form.
Definition: oem.h:78
invlib::Matrix< ArtsMatrixReference<::Matrix > > MatrixReference
invlib wrapper type for ARTS matrices to be passed by reference.
Definition: oem.h:35
void Tensor4Clip(Tensor4 &x, const Index &iq, const Numeric &limit_low, const Numeric &limit_high)
Clip Tensor4.
Definition: oem.h:581
#define u
#define v
#define w
#define a