38#pragma GCC diagnostic ignored "-Wconversion"
47 const String& bulkprop_name,
48 const Numeric& limit_low,
49 const Numeric& limit_high,
52 if (bulkprop_name ==
"ALL") {
56 for (Index i = 0; i < particle_bulkprop_names.
nelem(); i++) {
57 if (particle_bulkprop_names[i] == bulkprop_name) {
63 "Could not find ", bulkprop_name,
64 " in particle_bulkprop_names.\n")
67 Tensor4Clip(particle_bulkprop_field, iq, limit_low, limit_high);
74 const Numeric& limit_low,
75 const Numeric& limit_high,
78 if (species ==
"ALL") {
82 for (Index i = 0; i < abs_species.
nelem(); i++) {
89 "Could not find ", species,
" in abs_species.\n")
99 const Numeric& limit_low,
100 const Numeric& limit_high,
103 const Index nq = jacobian_quantities.
nelem();
106 ARTS_USER_ERROR_IF (ijq >= nq,
107 "Argument *ijq* is too high.\n
"
108 "You have selected index:
", ijq, "\n
"
109 "but the number of quantities is only:
", nq, "\n
"
110 "(Note that zero-based indexing is used)\n
")
113 ArrayOfArrayOfIndex ji;
116 jac_ranges_indices(ji, any_affine, jacobian_quantities);
119 Index ifirst = 0, ilast = x.nelem() - 1;
125 if (!std::isinf(limit_low)) {
126 for (Index i = ifirst; i <= ilast; i++) {
127 if (x[i] < limit_low) x[i] = limit_low;
130 if (!std::isinf(limit_high)) {
131 for (Index i = ifirst; i <= ilast; i++) {
132 if (x[i] > limit_high) x[i] = limit_high;
137/* Workspace method: Doxygen documentation will be auto-generated */
138void xaStandard(Workspace& ws,
140 const ArrayOfRetrievalQuantity& jacobian_quantities,
141 const Index& atmfields_checked,
142 const Index& atmgeom_checked,
143 const Index& atmosphere_dim,
144 const Vector& p_grid,
145 const Vector& lat_grid,
146 const Vector& lon_grid,
147 const Tensor3& t_field,
148 const Tensor4& vmr_field,
149 const ArrayOfArrayOfSpeciesTag& abs_species,
150 const Index& cloudbox_on,
151 const Index& cloudbox_checked,
152 const Tensor4& particle_bulkprop_field,
153 const ArrayOfString& particle_bulkprop_names,
154 const Tensor3& wind_u_field,
155 const Tensor3& wind_v_field,
156 const Tensor3& wind_w_field,
157 const Tensor3& mag_u_field,
158 const Tensor3& mag_v_field,
159 const Tensor3& mag_w_field,
160 const Tensor3& surface_props_data,
161 const ArrayOfString& surface_props_names,
162 const Agenda& water_p_eq_agenda,
166 ARTS_USER_ERROR_IF (atmfields_checked != 1,
167 "The atmospheric fields must be flagged to have
"
168 "passed
a consistency check (atmfields_checked=1).
");
169 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
170 "The atmospheric geometry must be flagged to have
"
171 "passed
a consistency check (atmgeom_checked=1).
");
172 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
173 "The cloudbox must be flagged to have
"
174 "passed
a consistency check (cloudbox_checked=1).
");
177 ArrayOfArrayOfIndex ji;
180 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
184 const Index nq = jacobian_quantities.nelem();
186 xa.resize(ji[nq - 1][1] + 1);
188 // Loop retrieval quantities and fill *xa*
189 for (Index q = 0; q < jacobian_quantities.nelem(); q++) {
190 // Index range of this retrieval quantity
191 const Index np = ji[q][1] - ji[q][0] + 1;
192 Range ind(ji[q][0], np);
194 // Atmospheric temperatures
195 if (jacobian_quantities[q] == Jacobian::Atm::Temperature) {
196 // Here we need to interpolate *t_field*
197 ArrayOfGridPos gp_p, gp_lat, gp_lon;
198 get_gp_atmgrids_to_rq(gp_p,
201 jacobian_quantities[q],
207 regrid_atmfield_by_gp(t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
212 else if (jacobian_quantities[q].Target().isSpeciesVMR()) {
213 // Index position of species
214 ArrayOfSpeciesTag atag(jacobian_quantities[q].Subtag());
215 const Index isp = chk_contains("abs_species
", abs_species, atag);
217 if (jacobian_quantities[q].Mode() == "rel
") {
218 // This one is simple, just a vector of ones
221 // For all remaining options we need to interpolate *vmr_field*
222 ArrayOfGridPos gp_p, gp_lat, gp_lon;
223 get_gp_atmgrids_to_rq(gp_p,
226 jacobian_quantities[q],
232 regrid_atmfield_by_gp(vmr_x,
234 vmr_field(isp, joker, joker, joker),
239 if (jacobian_quantities[q].Mode() == "vmr
") {
240 flat(xa[ind], vmr_x);
241 } else if (jacobian_quantities[q].Mode() == "nd
") {
242 // Here we need to also interpolate *t_field*
244 regrid_atmfield_by_gp(
245 t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
246 // Calculate number density for species (vmr*nd_tot)
248 for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
249 for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
250 for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
253 number_density(jacobian_quantities[q].Grids()[0][i1],
259 } else if (jacobian_quantities[q].Mode() == "rh
") {
260 // Here we need to also interpolate *t_field*
262 regrid_atmfield_by_gp(
263 t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
265 water_p_eq_agendaExecute(ws, water_p_eq, t_x, water_p_eq_agenda);
266 // Calculate relative humidity (vmr*p/p_sat)
268 for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
269 for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
270 for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
271 xa[ji[q][0] + i] = vmr_x(i1, i2, i3) *
272 jacobian_quantities[q].Grids()[0][i1] /
273 water_p_eq(i1, i2, i3);
278 } else if (jacobian_quantities[q].Mode() == "q
") {
279 // Calculate specific humidity q, from mixing ratio r and
280 // vapour pressure e, as
281 // q = r(1+r); r = 0.622e/(p-e); e = vmr*p;
283 for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
284 for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
285 for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
287 vmr_x(i1, i2, i3) * jacobian_quantities[q].Grids()[0][i1];
289 0.622 * e / (jacobian_quantities[q].Grids()[0][i1] - e);
290 xa[ji[q][0] + i] = r / (1 + r);
301 // Scattering species
302 else if (jacobian_quantities[q] == Jacobian::Special::ScatteringString) {
304 ARTS_USER_ERROR_IF (particle_bulkprop_field.empty(),
305 "One jacobian quantity belongs to the
"
306 "scattering species category, but *particle_bulkprop_field*
"
308 ARTS_USER_ERROR_IF (particle_bulkprop_field.nbooks() !=
309 particle_bulkprop_names.nelem(),
310 "Mismatch in size between
"
311 "*particle_bulkprop_field* and *particle_bulkprop_names*.
");
313 const Index isp = find_first(particle_bulkprop_names,
314 jacobian_quantities[q].SubSubtag());
315 ARTS_USER_ERROR_IF (isp < 0,
316 "Jacobian quantity with index
", q, " covers
a "
317 "scattering species, and the field quantity is set to \
"",
318 jacobian_quantities[q].SubSubtag(),
"\", but this quantity "
319 "could not found in *particle_bulkprop_names*.")
325 jacobian_quantities[q],
333 particle_bulkprop_field(isp, joker, joker, joker),
337 flat(xa[ind], pbp_x);
344 else if (jacobian_quantities[q].Target().isWind()) {
345 ConstTensor3View source_field(wind_u_field);
346 if (jacobian_quantities[q] == Jacobian::Atm::WindV) {
347 source_field.set(wind_v_field);
348 }
else if (jacobian_quantities[q] == Jacobian::Atm::WindW) {
349 source_field.set(wind_w_field);
358 jacobian_quantities[q],
366 wind_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
367 flat(xa[ind], wind_x);
371 else if (jacobian_quantities[q].Target().isMagnetic()) {
372 if (jacobian_quantities[q] == Jacobian::Atm::MagneticMagnitude) {
379 jacobian_quantities[q],
386 Tensor3 mag_u, mag_v, mag_w;
388 mag_u, atmosphere_dim, mag_u_field, gp_p, gp_lat, gp_lon);
390 mag_v, atmosphere_dim, mag_v_field, gp_p, gp_lat, gp_lon);
392 mag_w, atmosphere_dim, mag_w_field, gp_p, gp_lat, gp_lon);
395 for (Index i = 0; i < gp_p.
nelem(); i++)
396 for (Index j = 0; j < gp_lat.
nelem(); j++)
397 for (Index k = 0;
k < gp_lon.
nelem();
k++)
398 mag_x(i, j, k) = std::hypot(mag_u(i, j, k), mag_v(i, j, k), mag_w(i, j, k));
399 flat(xa[ind], mag_x);
401 ConstTensor3View source_field(mag_u_field);
402 if (jacobian_quantities[q] == Jacobian::Atm::MagneticV) {
403 source_field.set(mag_v_field);
404 }
else if (jacobian_quantities[q] == Jacobian::Atm::MagneticW) {
405 source_field.set(mag_w_field);
414 jacobian_quantities[q],
422 mag_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
423 flat(xa[ind], mag_x);
428 else if (jacobian_quantities[q] == Jacobian::Special::SurfaceString) {
433 surface_props_names);
435 "One jacobian quantity belongs to the "
436 "surface category, but *surface_props_data* is empty.");
439 find_first(surface_props_names, jacobian_quantities[q].Subtag());
441 "Jacobian quantity with index ", q,
" covers a "
442 "surface property, and the field Subtag is set to \"",
443 jacobian_quantities[q].Subtag(),
"\", but this quantity "
444 "could not found in *surface_props_names*.")
449 jacobian_quantities[q],
456 surface_props_data(isu, joker, joker),
459 flat(xa[ind], surf_x);
464 else if (jacobian_quantities[q].Target().isPointing() ||
465 jacobian_quantities[q].Target().isFrequency() ||
466 jacobian_quantities[q] ==
Jacobian::Sensor::Polyfit ||
467 jacobian_quantities[q] ==
Jacobian::Sensor::Sinefit) {
473 "Found a retrieval quantity that is not yet handled by\n"
474 "internal retrievals: ", jacobian_quantities[q],
'\n')
486 Tensor4& particle_bulkprop_field,
487 Tensor3& wind_u_field,
488 Tensor3& wind_v_field,
489 Tensor3& wind_w_field,
490 Tensor3& mag_u_field,
491 Tensor3& mag_v_field,
492 Tensor3& mag_w_field,
493 Tensor3& surface_props_data,
496 const Index& atmfields_checked,
497 const Index& atmgeom_checked,
498 const Index& atmosphere_dim,
499 const Vector& p_grid,
500 const Vector& lat_grid,
501 const Vector& lon_grid,
503 const Index& cloudbox_on,
504 const Index& cloudbox_checked,
507 const Agenda& water_p_eq_agenda,
512 "The atmospheric fields must be flagged to have "
513 "passed a consistency check (atmfields_checked=1).");
515 "The atmospheric geometry must be flagged to have "
516 "passed a consistency check (atmgeom_checked=1).");
518 "The cloudbox must be flagged to have "
519 "passed a consistency check (cloudbox_checked=1).");
526 const Index nq = jacobian_quantities.
nelem();
537 "Length of *x* does not match length implied by "
538 "*jacobian_quantities*.");
544 for (Index q = 0; q < nq; q++) {
546 const Index np = ji[q][1] - ji[q][0] + 1;
547 Range ind(ji[q][0], np);
551 if (jacobian_quantities[q] == Jacobian::Atm::Temperature) {
555 Index n_p, n_lat, n_lon;
562 jacobian_quantities[q].Grids(),
569 Tensor3 t_x(n_p, n_lat, n_lon);
572 t_field, atmosphere_dim, t_x, gp_p, gp_lat, gp_lon);
577 else if (jacobian_quantities[q].Target().isSpeciesVMR()) {
580 const Index isp =
chk_contains(
"abs_species", abs_species, atag);
583 Tensor3 x_field(vmr_field.npages(), vmr_field.nrows(), vmr_field.ncols());
586 Index n_p, n_lat, n_lon;
593 jacobian_quantities[q].Grids(),
599 Tensor3 t3_x(n_p, n_lat, n_lon);
602 x_field, atmosphere_dim, t3_x, gp_p, gp_lat, gp_lon);
605 if (jacobian_quantities[q].Mode() ==
"rel") {
607 vmr_field(isp, joker, joker, joker) *= x_field;
608 }
else if (jacobian_quantities[q].Mode() ==
"vmr") {
610 vmr_field(isp, joker, joker, joker) = x_field;
611 }
else if (jacobian_quantities[q].Mode() ==
"nd") {
613 for (Index i3 = 0; i3 < vmr_field.ncols(); i3++) {
614 for (Index i2 = 0; i2 < vmr_field.nrows(); i2++) {
615 for (Index i1 = 0; i1 < vmr_field.npages(); i1++) {
616 vmr_field(isp, i1, i2, i3) =
617 x_field(i1, i2, i3) /
622 }
else if (jacobian_quantities[q].Mode() ==
"rh") {
626 for (Index i3 = 0; i3 < vmr_field.ncols(); i3++) {
627 for (Index i2 = 0; i2 < vmr_field.nrows(); i2++) {
628 for (Index i1 = 0; i1 < vmr_field.npages(); i1++) {
629 vmr_field(isp, i1, i2, i3) =
630 x_field(i1, i2, i3) * water_p_eq(i1, i2, i3) / p_grid[i1];
634 }
else if (jacobian_quantities[q].Mode() ==
"q") {
639 for (Index i3 = 0; i3 < vmr_field.ncols(); i3++) {
640 for (Index i2 = 0; i2 < vmr_field.nrows(); i2++) {
641 for (Index i1 = 0; i1 < vmr_field.npages(); i1++) {
642 const Numeric r = x_field(i1, i2, i3) / (1 - x_field(i1, i2, i3));
643 const Numeric e = r * p_grid[i1] / (0.622 + r);
644 vmr_field(isp, i1, i2, i3) = e / p_grid[i1];
655 else if (jacobian_quantities[q] == Jacobian::Special::ScatteringString) {
659 "One jacobian quantity belongs to the "
660 "scattering species category, but *particle_bulkprop_field* "
663 particle_bulkprop_names.
nelem(),
664 "Mismatch in size between "
665 "*particle_bulkprop_field* and *particle_bulkprop_field*.");
667 const Index isp =
find_first(particle_bulkprop_names,
668 jacobian_quantities[q].SubSubtag());
670 "Jacobian quantity with index ", q,
" covers a "
671 "scattering species, and the field quantity is set to \"",
672 jacobian_quantities[q].SubSubtag(),
"\", but this quantity "
673 "could not found in *particle_bulkprop_names*.")
678 Index n_p, n_lat, n_lon;
685 jacobian_quantities[q].Grids(),
691 Tensor3 pbfield_x(n_p, n_lat, n_lon);
695 pbfield, atmosphere_dim, pbfield_x, gp_p, gp_lat, gp_lon);
696 particle_bulkprop_field(isp, joker, joker, joker) = pbfield;
702 else if (jacobian_quantities[q].Target().isWind()) {
706 Index n_p, n_lat, n_lon;
713 jacobian_quantities[q].Grids(),
720 Tensor3 wind_x(n_p, n_lat, n_lon);
723 Tensor3View target_field(wind_u_field);
726 target_field.npages(), target_field.nrows(), target_field.ncols());
728 wind_field, atmosphere_dim, wind_x, gp_p, gp_lat, gp_lon);
730 if (jacobian_quantities[q] == Jacobian::Atm::WindU) {
731 wind_u_field = wind_field;
732 }
else if (jacobian_quantities[q] == Jacobian::Atm::WindV) {
733 wind_v_field = wind_field;
734 }
else if (jacobian_quantities[q] == Jacobian::Atm::WindW) {
735 wind_w_field = wind_field;
741 else if (jacobian_quantities[q].Target().isMagnetic()) {
745 Index n_p, n_lat, n_lon;
752 jacobian_quantities[q].Grids(),
759 Tensor3 mag_x(n_p, n_lat, n_lon);
762 Tensor3View target_field(mag_u_field);
765 target_field.npages(), target_field.nrows(), target_field.ncols());
767 mag_field, atmosphere_dim, mag_x, gp_p, gp_lat, gp_lon);
768 if (jacobian_quantities[q] == Jacobian::Atm::MagneticU) {
769 mag_u_field = mag_field;
770 }
else if (jacobian_quantities[q] == Jacobian::Atm::MagneticV) {
771 mag_v_field = mag_field;
772 }
else if (jacobian_quantities[q] == Jacobian::Atm::MagneticW) {
773 mag_w_field = mag_field;
774 }
else if (jacobian_quantities[q] == Jacobian::Atm::MagneticMagnitude) {
775 for (Index i = 0; i < n_p; i++) {
776 for (Index j = 0; j < n_lat; j++) {
777 for (Index k = 0; k < n_lon; k++) {
778 const Numeric scale = mag_x(i, j, k) /
779 std::hypot(mag_u_field(i, j, k), mag_v_field(i, j, k), mag_w_field(i, j, k));
780 mag_u_field(i, j, k) *= scale;
781 mag_v_field(i, j, k) *= scale;
782 mag_w_field(i, j, k) *= scale;
793 else if (jacobian_quantities[q] == Jacobian::Special::SurfaceString) {
798 surface_props_names);
800 "One jacobian quantity belongs to the "
801 "surface category, but *surface_props_data* is empty.");
804 find_first(surface_props_names, jacobian_quantities[q].Subtag());
806 "Jacobian quantity with index ", q,
" covers a "
807 "surface property, and the field Subtag is set to \"",
808 jacobian_quantities[q].Subtag(),
"\", but this quantity "
809 "could not found in *surface_props_names*.")
819 jacobian_quantities[q].Grids(),
824 Matrix surf_x(n_lat, n_lon);
828 surface_props_data(isu, joker, joker) = surf;
838 Sparse& sensor_response,
839 Vector& sensor_response_f,
841 Matrix& sensor_response_dlos,
842 Vector& sensor_response_f_grid,
844 Matrix& sensor_response_dlos_grid,
848 const Agenda& sensor_response_agenda,
849 const Index& sensor_checked,
855 "The sensor response must be flagged to have "
856 "passed a consistency check (sensor_checked=1).");
863 const Index nq = jacobian_quantities.
nelem();
874 "Length of *x* does not match length implied by "
875 "*jacobian_quantities*.");
881 bool do_sensor =
false;
884 for (Index q = 0; q < nq; q++) {
886 const Index np = ji[q][1] - ji[q][0] + 1;
890 if (jacobian_quantities[q].Target().isPointing()) {
892 if (jacobian_quantities[q].Grids()[0][0] == -1) {
894 "Mismatch between pointing jacobian and *sensor_los* found.");
896 for (Index i = 0; i < np; i++) {
897 sensor_los(i, 0) += x_t[ji[q][0] + i];
903 "Sizes of *sensor_los* and *sensor_time* do not match.");
905 for (Index
c = 0;
c < np;
c++) {
907 for (Index i = 0; i <
w.nelem(); i++) {
908 sensor_los(i, 0) +=
w[i] * x_t[ji[q][0] +
c];
916 else if (jacobian_quantities[q].Target().isFrequency()) {
917 if (jacobian_quantities[q] == Jacobian::Sensor::FrequencyShift) {
919 if (x_t[ji[q][0]] != 0) {
921 f_backend += x_t[ji[q][0]];
923 }
else if (jacobian_quantities[q] == Jacobian::Sensor::FrequencyStretch) {
925 if (x_t[ji[q][0]] != 0) {
929 for (Index i = 0; i <
w.nelem(); i++) {
930 f_backend[i] +=
w[i] * x_t[ji[q][0]];
940 else if (jacobian_quantities[q] == Jacobian::Sensor::Polyfit ||
941 jacobian_quantities[q] == Jacobian::Sensor::Sinefit) {
944 Index y_size = sensor_los.nrows() * sensor_response_f_grid.nelem() *
945 sensor_response_pol_grid.
nelem() *
946 sensor_response_dlos_grid.nrows();
947 y_baseline.resize(y_size);
951 for (Index mb = 0; mb < sensor_los.nrows(); ++mb) {
956 sensor_response_pol_grid,
957 sensor_response_f_grid,
958 sensor_response_dlos_grid,
959 jacobian_quantities[q],
968 y_baseline.resize(1);
977 sensor_response_f_grid,
979 sensor_response_pol_grid,
980 sensor_response_dlos,
981 sensor_response_dlos_grid,
984 sensor_response_agenda);
990 ARTS_USER_ERROR (
"Retrievals of spectroscopic variables not yet handled.");
1001 Vector& oem_diagnostics,
1002 Vector& lm_ga_history,
1009 const Agenda& inversion_iterate_agenda,
1011 const Numeric& max_start_cost,
1012 const Vector& x_norm,
1013 const Index& max_iter,
1014 const Numeric& stop_dx,
1015 const Vector& lm_ga_settings,
1016 const Index& clear_matrices,
1017 const Index& display_progress,
1020 const Index n = covmat_sx.
nrows();
1021 const Index m = y.nelem();
1031 inversion_iterate_agenda,
1036 jacobian_quantities,
1046 oem_diagnostics.resize(5);
1047 oem_diagnostics = NAN;
1049 if (method ==
"ml" || method ==
"lm" || method ==
"ml_cg" ||
1050 method ==
"lm_cg") {
1051 lm_ga_history.resize(max_iter + 1);
1052 lm_ga_history = NAN;
1054 lm_ga_history.resize(0);
1058 if (x.nelem() != n) {
1061 jacobian.resize(0, 0);
1066 if (yf.nelem() == 0) {
1068 ws, yf, jacobian, xa, 1, 0, inversion_iterate_agenda);
1072 "Mismatch between simulated y and input y.\n"
1073 "Input y is size ", y.nelem(),
" but simulated y is ",
1075 "Use your frequency grid vector and your sensor response matrix to match simulations with measurements.\n")
1079 Numeric cost_start = NAN;
1080 if (method == "ml" || method == "lm" || display_progress ||
1081 max_start_cost > 0) {
1085 mult_inv(ExhaustiveMatrixView{sdy}, covmat_se, ExhaustiveMatrixView{dy});
1089 mult_inv(ExhaustiveMatrixView{sdx}, covmat_sx, ExhaustiveMatrixView{dx});
1090 cost_start = dx * sdx + dy * sdy;
1091 cost_start /=
static_cast<Numeric
>(m);
1093 oem_diagnostics[1] = cost_start;
1096 if (max_start_cost > 0 && cost_start > max_start_cost) {
1098 oem_diagnostics[0] = 99;
1100 if (display_progress) {
1101 cout <<
"\n No OEM inversion, too high start cost:\n"
1102 <<
" Set limit : " << max_start_cost << endl
1103 <<
" Found value : " << cost_start << endl
1109 bool apply_norm =
false;
1111 if (x_norm.nelem() == n) {
1114 T.diagonal() = x_norm;
1115 for (Index i = 0; i < n; i++) {
1116 T(i, i) = x_norm[i];
1128 &inversion_iterate_agenda);
1131 int oem_verbosity =
static_cast<int>(display_progress);
1133 int return_code = 0;
1136 if (method ==
"li") {
1140 x_oem, y_oem, gn, oem_verbosity, lm_ga_history,
true);
1141 oem_diagnostics[0] =
static_cast<Index
>(return_code);
1142 }
else if (method ==
"li_m") {
1144 oem::Std s(T, apply_norm);
1145 oem::GN gn(stop_dx, 1, s);
1148 oem_diagnostics[0] = static_cast<Index>(return_code);
1149 } else if (method == "li_cg") {
1150 oem::CG cg(T, apply_norm, 1e-10, 0);
1153 x_oem, y_oem, gn, oem_verbosity, lm_ga_history,
true);
1154 oem_diagnostics[0] =
static_cast<Index
>(return_code);
1155 }
else if (method ==
"li_cg_m") {
1156 oem::CG cg(T, apply_norm, 1e-10, 0);
1159 x_oem, y_oem, gn, oem_verbosity, lm_ga_history,
true);
1160 oem_diagnostics[0] =
static_cast<Index
>(return_code);
1161 }
else if (method ==
"gn") {
1163 oem::GN gn(stop_dx, (
unsigned int)max_iter, s);
1165 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1166 oem_diagnostics[0] =
static_cast<Index
>(return_code);
1167 }
else if (method ==
"gn_m") {
1169 oem::Std s(T, apply_norm);
1170 oem::GN gn(stop_dx, (
unsigned int)max_iter, s);
1173 oem_diagnostics[0] = static_cast<Index>(return_code);
1174 } else if (method == "gn_cg") {
1175 oem::CG cg(T, apply_norm, 1e-10, 0);
1176 oem::GN_CG gn(stop_dx, (
unsigned int)max_iter, cg);
1178 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1179 oem_diagnostics[0] =
static_cast<Index
>(return_code);
1180 }
else if (method ==
"gn_cg_m") {
1181 oem::CG cg(T, apply_norm, 1e-10, 0);
1182 oem::GN_CG gn(stop_dx, (
unsigned int)max_iter, cg);
1184 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1185 oem_diagnostics[0] =
static_cast<Index
>(return_code);
1186 }
else if ((method ==
"lm") || (method ==
"ml")) {
1192 std::make_pair(0, 0),
1193 make_shared<Sparse>(diagonal)));
1197 lm.set_tolerance(stop_dx);
1198 lm.set_maximum_iterations((
unsigned int)max_iter);
1199 lm.set_lambda(lm_ga_settings[0]);
1200 lm.set_lambda_decrease(lm_ga_settings[1]);
1201 lm.set_lambda_increase(lm_ga_settings[2]);
1202 lm.set_lambda_maximum(lm_ga_settings[3]);
1203 lm.set_lambda_threshold(lm_ga_settings[4]);
1204 lm.set_lambda_constraint(lm_ga_settings[5]);
1207 x_oem, y_oem, lm, oem_verbosity, lm_ga_history);
1208 oem_diagnostics[0] =
static_cast<Index
>(return_code);
1209 if (lm.get_lambda() > lm.get_lambda_maximum()) {
1210 oem_diagnostics[0] = 2;
1212 }
else if ((method ==
"lm_cg") || (method ==
"ml_cg")) {
1213 oem::CG cg(T, apply_norm, 1e-10, 0);
1219 std::make_pair(0, 0),
1220 make_shared<Sparse>(diagonal)));
1223 lm.set_maximum_iterations((
unsigned int)max_iter);
1224 lm.set_lambda(lm_ga_settings[0]);
1225 lm.set_lambda_decrease(lm_ga_settings[1]);
1226 lm.set_lambda_increase(lm_ga_settings[2]);
1227 lm.set_lambda_threshold(lm_ga_settings[3]);
1228 lm.set_lambda_maximum(lm_ga_settings[4]);
1231 x_oem, y_oem, lm, oem_verbosity, lm_ga_history);
1232 oem_diagnostics[0] =
static_cast<Index
>(return_code);
1233 if (lm.get_lambda() > lm.get_lambda_maximum()) {
1234 oem_diagnostics[0] = 2;
1238 oem_diagnostics[2] =
oem.cost /
static_cast<Numeric
>(m);
1239 oem_diagnostics[3] =
oem.cost_y /
static_cast<Numeric
>(m);
1240 oem_diagnostics[4] =
static_cast<Numeric
>(
oem.iterations);
1241 }
catch (
const std::exception& e) {
1242 oem_diagnostics[0] = 9;
1243 oem_diagnostics[2] =
oem.cost;
1244 oem_diagnostics[3] =
oem.cost_y;
1245 oem_diagnostics[4] =
static_cast<Numeric
>(
oem.iterations);
1248 for (
auto& s : sv) {
1249 std::stringstream ss{s};
1251 while (std::getline(ss, t)) {
1252 errors.push_back(t.c_str());
1260 yf = aw.get_measurement_vector();
1263 if (clear_matrices) {
1264 jacobian.resize(0, 0);
1266 }
else if (oem_diagnostics[0] <= 2) {
1268 Matrix tmp1(n, m), tmp2(n, n), tmp3(n, n);
1269 mult_inv(tmp1, transpose(jacobian), covmat_se);
1270 mult(tmp2, tmp1, jacobian);
1273 mult(dxdy, tmp3, tmp1);
1283 Index n(dxdy.nrows()), m(dxdy.ncols());
1287 "The gain matrix *dxdy* is required to compute the observation error covariance matrix.");
1289 "The covariance matrix covmat_se has invalid dimensions.");
1291 covmat_so.resize(n, n);
1292 mult(tmp1, covmat_se, transpose(dxdy));
1293 mult(covmat_so, dxdy, tmp1);
1301 Index n(avk.ncols());
1302 Matrix tmp1(n, n), tmp2(n, n);
1305 "The averaging kernel matrix *dxdy* is required to compute the smoothing error covariance matrix.");
1307 "The covariance matrix *covmat_sx* invalid dimensions.");
1309 covmat_ss.resize(n, n);
1315 mult(tmp2, covmat_sx, transpose(tmp1));
1316 mult(covmat_ss, tmp1, tmp2);
1329 const Matrix& jacobian,
1331 Index m(jacobian.nrows()), n(jacobian.ncols());
1333 "The Jacobian matrix is empty.");
1335 "Matrices have inconsistent sizes.\n"
1336 " Size of gain matrix: ", dxdy.nrows(),
" x ", dxdy.ncols(),
1338 " Size of Jacobian: ", jacobian.nrows(),
" x ",
1339 jacobian.ncols(),
"\n")
1342 mult(avk, dxdy, jacobian);
1352 "WSM is not available because ARTS was compiled without "
1361 "WSM is not available because ARTS was compiled without "
1370 "WSM is not available because ARTS was compiled without "
1400 "WSM is not available because ARTS was compiled without "
1405#if defined(OEM_SUPPORT) && 0
1407#include "agenda_wrapper_mpi.h"
1417void MPI_Initialize(Matrix& sensor_los,
1419 Vector& sensor_time) {
1422 MPI_Initialized(&initialized);
1424 MPI_Init(
nullptr,
nullptr);
1429 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1430 MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
1432 int nmblock = (int)sensor_pos.nrows();
1433 int mblock_range = nmblock / nprocs;
1434 int mblock_start = mblock_range * rank;
1435 int remainder = nmblock % std::max(mblock_range, nprocs);
1441 if (rank < remainder) {
1443 mblock_start += rank;
1445 mblock_start += remainder;
1449 Range range = Range(mblock_start, mblock_range);
1451 Matrix tmp_m = sensor_los(range, joker);
1454 tmp_m = sensor_pos(range, joker);
1457 Vector tmp_v = sensor_time[range];
1458 sensor_time = tmp_v;
1460 sensor_los.resize(0, 0);
1461 sensor_pos.resize(0, 0);
1462 sensor_time.resize(0);
1471 Vector& oem_diagnostics,
1472 Vector& lm_ga_history,
1475 Vector& sensor_time,
1481 const Agenda& inversion_iterate_agenda,
1483 const Numeric& max_start_cost,
1484 const Vector& x_norm,
1485 const Index& max_iter,
1486 const Numeric& stop_dx,
1487 const Vector& lm_ga_settings,
1488 const Index& clear_matrices,
1489 const Index& display_progress,
1492 const Index n = covmat_sx.
nrows();
1493 const Index m = y.nelem();
1500 inversion_iterate_agenda,
1504 jacobian_quantities,
1516 jacobian.resize(0, 0);
1519 MPI_Initialize(sensor_los, sensor_pos, sensor_time);
1529 ws, tmp, jacobian, xa, 1, inversion_iterate_agenda);
1533 oem_diagnostics.resize(5);
1534 oem_diagnostics = NAN;
1536 if (method ==
"ml" || method ==
"lm") {
1537 lm_ga_history.resize(max_iter);
1538 lm_ga_history = NAN;
1540 lm_ga_history.resize(0);
1545 Numeric cost_start = NAN;
1546 if (method ==
"ml" || method ==
"lm" || display_progress ||
1547 max_start_cost > 0) {
1550 cost_start = dot(dy, SeInvMPI * dy);
1552 oem_diagnostics[1] = cost_start;
1555 if (max_start_cost > 0 && cost_start > max_start_cost) {
1557 oem_diagnostics[0] = 99;
1559 if (display_progress) {
1560 cout <<
"\n No OEM inversion, too high start cost:\n"
1561 <<
" Set limit : " << max_start_cost << endl
1562 <<
" Found value : " << cost_start << endl
1573 OEMVector xa_oem(xa), y_oem(y), x_oem;
1576 OEM_PS_PS_MPI<AgendaWrapperMPI>
oem(aw, xa_oem, SaInvMPI, SeInvMPI);
1579 int return_value = 99;
1581 if (method ==
"li") {
1583 oem::GN_CG gn(stop_dx, (
unsigned int)max_iter, cg);
1585 x_oem, y_oem, gn, 2 * (int)display_progress);
1586 }
else if (method ==
"gn") {
1588 oem::GN_CG gn(stop_dx, (
unsigned int)max_iter, cg);
1590 x_oem, y_oem, gn, 2 * (int)display_progress);
1591 }
else if ((method ==
"lm") || (method ==
"ml")) {
1593 LM_CG_MPI lm(SaInvMPI, cg);
1595 lm.set_tolerance(stop_dx);
1596 lm.set_maximum_iterations((
unsigned int)max_iter);
1597 lm.set_lambda(lm_ga_settings[0]);
1598 lm.set_lambda_decrease(lm_ga_settings[1]);
1599 lm.set_lambda_increase(lm_ga_settings[2]);
1600 lm.set_lambda_threshold(lm_ga_settings[3]);
1601 lm.set_lambda_maximum(lm_ga_settings[4]);
1603 return_value =
oem.compute<oem::LM_CG_MPI, invlib::MPILog>(
1604 x_oem, y_oem, lm, 2 * (int)display_progress);
1607 oem_diagnostics[0] = return_value;
1608 oem_diagnostics[2] =
oem.cost;
1609 oem_diagnostics[3] =
oem.cost_y;
1610 oem_diagnostics[4] =
static_cast<Numeric
>(
oem.iterations);
1614 if (clear_matrices && (oem_diagnostics[0])) {
1615 jacobian.resize(0, 0);
1651 "You have to compile ARTS with OEM support "
1652 " and enable MPI to use OEM_MPI.");
This file contains the definition of Array.
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
The global header file for ARTS.
Header file for helper functions for OpenMP.
Vector time_vector(const ArrayOfTime ×)
Converts from each Time to seconds and returns as Vector.
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)
void water_p_eq_agendaExecute(Workspace &ws, Tensor3 &water_p_eq_field, const Tensor3 &t_field, const Agenda &input_agenda)
void sensor_response_agendaExecute(Workspace &ws, Sparse &sensor_response, Vector &sensor_response_f, Vector &sensor_response_f_grid, ArrayOfIndex &sensor_response_pol, ArrayOfIndex &sensor_response_pol_grid, Matrix &sensor_response_dlos, Matrix &sensor_response_dlos_grid, Matrix &mblock_dlos, const Vector &f_backend, const Agenda &input_agenda)
void MatrixFromCovarianceMatrix(Matrix &output, const CovarianceMatrix &input, const Verbosity &verbosity)
WORKSPACE METHOD: MatrixFromCovarianceMatrix.
Index nelem() const ARTS_NOEXCEPT
void compute_inverse() const
Compute the inverse of this correlation matrix.
void add_correlation_inverse(Block c)
Add block inverse to covariance matrix.
Vector inverse_diagonal() const
Diagonal of the inverse of the covariance matrix as vector.
Interface for distributed ARTS forward model.
Interface to ARTS inversion_iterate_agenda.
void mult(MatrixView C, ConstMatrixView A, const Block &B)
void add_inv(MatrixView A, const CovarianceMatrix &B)
void mult_inv(MatrixView C, ConstMatrixView A, const CovarianceMatrix &B)
Helper macros for debugging.
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
void polynomial_basis_func(Vector &b, const Vector &x, const Index &poly_coeff)
Calculates polynomial basis functions.
void transform_x_back(Vector &x_t, const ArrayOfRetrievalQuantity &jqs, bool revert_functional_transforms)
Handles back-transformations of the state vector.
void jac_ranges_indices(ArrayOfArrayOfIndex &jis, bool &any_affine, const ArrayOfRetrievalQuantity &jqs, const bool &before_affine)
Determines the index range inside x and the Jacobian for each retrieval quantity.
void transform_x(Vector &x, const ArrayOfRetrievalQuantity &jqs)
Handles transformations of the state vector.
void calcBaselineFit(Vector &y_baseline, const Vector &x, const Index &mblock_index, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const RetrievalQuantity &rq, const Index rq_index, const ArrayOfArrayOfIndex &jacobian_indices)
Calculate baseline fit.
Routines for setting up the jacobian.
void covmat_ssCalc(Matrix &, const Matrix &, const CovarianceMatrix &, const Verbosity &)
WORKSPACE METHOD: covmat_ssCalc.
void particle_bulkprop_fieldClip(Tensor4 &particle_bulkprop_field, const ArrayOfString &particle_bulkprop_names, const String &bulkprop_name, const Numeric &limit_low, const Numeric &limit_high, const Verbosity &)
WORKSPACE METHOD: particle_bulkprop_fieldClip.
void x2artsSensor(Workspace &ws, Matrix &sensor_los, Vector &f_backend, Vector &y_baseline, Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, ArrayOfIndex &sensor_response_pol_grid, Matrix &sensor_response_dlos_grid, Matrix &mblock_dlos, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &x, const Agenda &sensor_response_agenda, const Index &sensor_checked, const ArrayOfTime &sensor_time, const Verbosity &)
WORKSPACE METHOD: x2artsSensor.
void OEM_MPI(Workspace &, Vector &, Vector &, Matrix &, Matrix &, Vector &, Vector &, Matrix &, Matrix &, Vector &, CovarianceMatrix &, CovarianceMatrix &, const Vector &, const Vector &, const Index &, const ArrayOfRetrievalQuantity &, const Agenda &, const String &, const Numeric &, const Vector &, const Index &, const Numeric &, const Vector &, const Index &, const Index &, const Verbosity &)
void OEM(Workspace &, Vector &, Vector &, Matrix &, Matrix &, Vector &, Vector &, ArrayOfString &, const Vector &, const CovarianceMatrix &, const Vector &, const CovarianceMatrix &, const Index &, const ArrayOfRetrievalQuantity &, const ArrayOfArrayOfIndex &, const Agenda &, const String &, const Numeric &, const Vector &, const Index &, const Numeric &, const Vector &, const Index &, const Index &, const Verbosity &)
void x2artsAtmAndSurf(Workspace &ws, Tensor4 &vmr_field, Tensor3 &t_field, Tensor4 &particle_bulkprop_field, Tensor3 &wind_u_field, Tensor3 &wind_v_field, Tensor3 &wind_w_field, Tensor3 &mag_u_field, Tensor3 &mag_v_field, Tensor3 &mag_w_field, Tensor3 &surface_props_data, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &x, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const Index &cloudbox_on, const Index &cloudbox_checked, const ArrayOfString &particle_bulkprop_names, const ArrayOfString &surface_props_names, const Agenda &water_p_eq_agenda, const Verbosity &)
WORKSPACE METHOD: x2artsAtmAndSurf.
void x2artsSpectroscopy(const Verbosity &)
WORKSPACE METHOD: x2artsSpectroscopy.
void avkCalc(Matrix &, const Matrix &, const Matrix &, const Verbosity &)
WORKSPACE METHOD: avkCalc.
void xClip(Vector &x, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &ijq, const Numeric &limit_low, const Numeric &limit_high, const Verbosity &)
WORKSPACE METHOD: xClip.
void covmat_soCalc(Matrix &, const Matrix &, const CovarianceMatrix &, const Verbosity &)
WORKSPACE METHOD: covmat_soCalc.
void vmr_fieldClip(Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const String &species, const Numeric &limit_low, const Numeric &limit_high, const Verbosity &)
WORKSPACE METHOD: vmr_fieldClip.
void reshape(Tensor3View X, ConstVectorView x)
reshape
void flat(VectorView x, ConstMatrixView X)
flat
constexpr Numeric k
Boltzmann constant convenience name [J/K].
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::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::Matrix< invlib::MPIMatrix< invlib::Timer< ArtsCovarianceMatrixWrapper > > > MPICovarianceMatrix
MPI-distributed covariance matrix type.
invlib::Vector< invlib::MPIVector< invlib::Timer< ArtsVector > > > MPIVector
MPI-distributed vector type.
Defines the ARTS interface to the invlib library.
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)
Error checking for OEM method.
void Tensor4Clip(Tensor4 &x, const Index &iq, const Numeric &limit_low, const Numeric &limit_high)
Clip Tensor4.
Optimal estimation method for MPI-distributed retrieval.
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Declaration of functions in rte.cc.
void regrid_atmfield_by_gp(Tensor3 &field_new, const Index &atmosphere_dim, ConstTensor3View field_old, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regrids an atmospheric field, for precalculated grid positions.
void regrid_atmfield_by_gp_oem(Tensor3 &field_new, const Index &atmosphere_dim, ConstTensor3View field_old, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regridding of atmospheric field OEM-type.
void get_gp_atmsurf_to_rq(ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, const RetrievalQuantity &rq, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric surfaces to retrieval grids.
void get_gp_rq_to_atmgrids(ArrayOfGridPos &gp_p, ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, Index &n_p, Index &n_lat, Index &n_lon, const ArrayOfVector &ret_grids, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric fields to retrieval grids.
void regrid_atmsurf_by_gp_oem(Matrix &field_new, const Index &atmosphere_dim, ConstMatrixView field_old, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regridding of surface field OEM-type.
void get_gp_atmgrids_to_rq(ArrayOfGridPos &gp_p, ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, const RetrievalQuantity &rq, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric fields to retrieval grids.
Header file for special_interp.cc.
void surface_props_check(const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names)
Peforms basic checks of surface_props_data and surface_props_names