56#pragma GCC diagnostic ignored "-Wconversion"
65 const String& bulkprop_name,
70 if (bulkprop_name ==
"ALL") {
74 for (
Index i = 0; i < particle_bulkprop_names.
nelem(); i++) {
75 if (particle_bulkprop_names[i] == bulkprop_name) {
81 "Could not find ", bulkprop_name,
82 " in particle_bulkprop_names.\n")
85 Tensor4Clip(particle_bulkprop_field, iq, limit_low, limit_high);
96 if (species ==
"ALL") {
100 for (
Index i = 0; i < abs_species.
nelem(); i++) {
107 "Could not find ", species,
" in abs_species.\n")
124 ARTS_USER_ERROR_IF (ijq >= nq,
125 "Argument *ijq* is too high.\n
"
126 "You have selected index:
", ijq, "\n
"
127 "but the number of quantities is only:
", nq, "\n
"
128 "(Note that zero-based indexing is used)\n
")
131 ArrayOfArrayOfIndex ji;
134 jac_ranges_indices(ji, any_affine, jacobian_quantities);
137 Index ifirst = 0, ilast = x.nelem() - 1;
143 if (!std::isinf(limit_low)) {
144 for (Index i = ifirst; i <= ilast; i++) {
145 if (x[i] < limit_low) x[i] = limit_low;
148 if (!std::isinf(limit_high)) {
149 for (Index i = ifirst; i <= ilast; i++) {
150 if (x[i] > limit_high) x[i] = limit_high;
155/* Workspace method: Doxygen documentation will be auto-generated */
156void xaStandard(Workspace& ws,
158 const ArrayOfRetrievalQuantity& jacobian_quantities,
159 const Index& atmfields_checked,
160 const Index& atmgeom_checked,
161 const Index& atmosphere_dim,
162 const Vector& p_grid,
163 const Vector& lat_grid,
164 const Vector& lon_grid,
165 const Tensor3& t_field,
166 const Tensor4& vmr_field,
167 const ArrayOfArrayOfSpeciesTag& abs_species,
168 const Index& cloudbox_on,
169 const Index& cloudbox_checked,
170 const Tensor4& particle_bulkprop_field,
171 const ArrayOfString& particle_bulkprop_names,
172 const Tensor3& wind_u_field,
173 const Tensor3& wind_v_field,
174 const Tensor3& wind_w_field,
175 const Tensor3& mag_u_field,
176 const Tensor3& mag_v_field,
177 const Tensor3& mag_w_field,
178 const Tensor3& surface_props_data,
179 const ArrayOfString& surface_props_names,
180 const Agenda& water_p_eq_agenda,
184 ARTS_USER_ERROR_IF (atmfields_checked != 1,
185 "The atmospheric fields must be flagged to have
"
186 "passed
a consistency check (atmfields_checked=1).
");
187 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
188 "The atmospheric geometry must be flagged to have
"
189 "passed
a consistency check (atmgeom_checked=1).
");
190 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
191 "The cloudbox must be flagged to have
"
192 "passed
a consistency check (cloudbox_checked=1).
");
195 ArrayOfArrayOfIndex ji;
198 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
202 const Index nq = jacobian_quantities.nelem();
204 xa.resize(ji[nq - 1][1] + 1);
206 // Loop retrieval quantities and fill *xa*
207 for (Index q = 0; q < jacobian_quantities.nelem(); q++) {
208 // Index range of this retrieval quantity
209 const Index np = ji[q][1] - ji[q][0] + 1;
210 Range ind(ji[q][0], np);
212 // Atmospheric temperatures
213 if (jacobian_quantities[q] == Jacobian::Atm::Temperature) {
214 // Here we need to interpolate *t_field*
215 ArrayOfGridPos gp_p, gp_lat, gp_lon;
216 get_gp_atmgrids_to_rq(gp_p,
219 jacobian_quantities[q],
225 regrid_atmfield_by_gp(t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
230 else if (jacobian_quantities[q].Target().isSpeciesVMR()) {
231 // Index position of species
232 ArrayOfSpeciesTag atag(jacobian_quantities[q].Subtag());
233 const Index isp = chk_contains("abs_species
", abs_species, atag);
235 if (jacobian_quantities[q].Mode() == "rel
") {
236 // This one is simple, just a vector of ones
239 // For all remaining options we need to interpolate *vmr_field*
240 ArrayOfGridPos gp_p, gp_lat, gp_lon;
241 get_gp_atmgrids_to_rq(gp_p,
244 jacobian_quantities[q],
250 regrid_atmfield_by_gp(vmr_x,
252 vmr_field(isp, joker, joker, joker),
257 if (jacobian_quantities[q].Mode() == "vmr
") {
258 flat(xa[ind], vmr_x);
259 } else if (jacobian_quantities[q].Mode() == "nd
") {
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);
264 // Calculate number density for species (vmr*nd_tot)
266 for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
267 for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
268 for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
271 number_density(jacobian_quantities[q].Grids()[0][i1],
277 } else if (jacobian_quantities[q].Mode() == "rh
") {
278 // Here we need to also interpolate *t_field*
280 regrid_atmfield_by_gp(
281 t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
283 water_p_eq_agendaExecute(ws, water_p_eq, t_x, water_p_eq_agenda);
284 // Calculate relative humidity (vmr*p/p_sat)
286 for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
287 for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
288 for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
289 xa[ji[q][0] + i] = vmr_x(i1, i2, i3) *
290 jacobian_quantities[q].Grids()[0][i1] /
291 water_p_eq(i1, i2, i3);
296 } else if (jacobian_quantities[q].Mode() == "q
") {
297 // Calculate specific humidity q, from mixing ratio r and
298 // vapour pressure e, as
299 // q = r(1+r); r = 0.622e/(p-e); e = vmr*p;
301 for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
302 for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
303 for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
305 vmr_x(i1, i2, i3) * jacobian_quantities[q].Grids()[0][i1];
307 0.622 * e / (jacobian_quantities[q].Grids()[0][i1] - e);
308 xa[ji[q][0] + i] = r / (1 + r);
319 // Scattering species
320 else if (jacobian_quantities[q] == Jacobian::Special::ScatteringString) {
322 ARTS_USER_ERROR_IF (particle_bulkprop_field.empty(),
323 "One jacobian quantity belongs to the
"
324 "scattering species category, but *particle_bulkprop_field*
"
326 ARTS_USER_ERROR_IF (particle_bulkprop_field.nbooks() !=
327 particle_bulkprop_names.nelem(),
328 "Mismatch in size between
"
329 "*particle_bulkprop_field* and *particle_bulkprop_names*.
");
331 const Index isp = find_first(particle_bulkprop_names,
332 jacobian_quantities[q].SubSubtag());
333 ARTS_USER_ERROR_IF (isp < 0,
334 "Jacobian quantity with index
", q, " covers
a "
335 "scattering species, and the field quantity is set to \
"",
336 jacobian_quantities[q].SubSubtag(),
"\", but this quantity "
337 "could not found in *particle_bulkprop_names*.")
343 jacobian_quantities[q],
355 flat(xa[ind], pbp_x);
362 else if (jacobian_quantities[q].Target().isWind()) {
364 if (jacobian_quantities[q] == Jacobian::Atm::WindV) {
365 source_field = wind_v_field;
366 }
else if (jacobian_quantities[q] == Jacobian::Atm::WindW) {
367 source_field = wind_w_field;
376 jacobian_quantities[q],
384 wind_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
385 flat(xa[ind], wind_x);
389 else if (jacobian_quantities[q].Target().isMagnetic()) {
390 if (jacobian_quantities[q] == Jacobian::Atm::MagneticMagnitude) {
397 jacobian_quantities[q],
406 mag_u, atmosphere_dim, mag_u_field, gp_p, gp_lat, gp_lon);
408 mag_v, atmosphere_dim, mag_v_field, gp_p, gp_lat, gp_lon);
410 mag_w, atmosphere_dim, mag_w_field, gp_p, gp_lat, gp_lon);
416 mag_x(i, j, k) = std::hypot(mag_u(i, j, k), mag_v(i, j, k), mag_w(i, j, k));
417 flat(xa[ind], mag_x);
420 if (jacobian_quantities[q] == Jacobian::Atm::MagneticV) {
421 source_field = mag_v_field;
422 }
else if (jacobian_quantities[q] == Jacobian::Atm::MagneticW) {
423 source_field = mag_w_field;
432 jacobian_quantities[q],
440 mag_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
441 flat(xa[ind], mag_x);
446 else if (jacobian_quantities[q] == Jacobian::Special::SurfaceString) {
451 surface_props_names);
453 "One jacobian quantity belongs to the "
454 "surface category, but *surface_props_data* is empty.");
457 find_first(surface_props_names, jacobian_quantities[q].Subtag());
459 "Jacobian quantity with index ", q,
" covers a "
460 "surface property, and the field Subtag is set to \"",
461 jacobian_quantities[q].Subtag(),
"\", but this quantity "
462 "could not found in *surface_props_names*.")
467 jacobian_quantities[q],
477 flat(xa[ind], surf_x);
482 else if (jacobian_quantities[q].Target().isPointing() ||
483 jacobian_quantities[q].Target().isFrequency() ||
484 jacobian_quantities[q] ==
Jacobian::Sensor::Polyfit ||
485 jacobian_quantities[q] ==
Jacobian::Sensor::Sinefit) {
491 "Found a retrieval quantity that is not yet handled by\n"
492 "internal retrievals: ", jacobian_quantities[q],
'\n')
504 Tensor4& particle_bulkprop_field,
514 const Index& atmfields_checked,
515 const Index& atmgeom_checked,
516 const Index& atmosphere_dim,
521 const Index& cloudbox_on,
522 const Index& cloudbox_checked,
525 const Agenda& water_p_eq_agenda,
530 "The atmospheric fields must be flagged to have "
531 "passed a consistency check (atmfields_checked=1).");
533 "The atmospheric geometry must be flagged to have "
534 "passed a consistency check (atmgeom_checked=1).");
536 "The cloudbox must be flagged to have "
537 "passed a consistency check (cloudbox_checked=1).");
555 "Length of *x* does not match length implied by "
556 "*jacobian_quantities*.");
562 for (
Index q = 0; q < nq; q++) {
564 const Index np = ji[q][1] - ji[q][0] + 1;
565 Range ind(ji[q][0], np);
569 if (jacobian_quantities[q] == Jacobian::Atm::Temperature) {
573 Index n_p, n_lat, n_lon;
580 jacobian_quantities[q].Grids(),
587 Tensor3 t_x(n_p, n_lat, n_lon);
590 t_field, atmosphere_dim, t_x, gp_p, gp_lat, gp_lon);
595 else if (jacobian_quantities[q].Target().isSpeciesVMR()) {
604 Index n_p, n_lat, n_lon;
611 jacobian_quantities[q].Grids(),
617 Tensor3 t3_x(n_p, n_lat, n_lon);
620 x_field, atmosphere_dim, t3_x, gp_p, gp_lat, gp_lon);
623 if (jacobian_quantities[q].Mode() ==
"rel") {
626 }
else if (jacobian_quantities[q].Mode() ==
"vmr") {
629 }
else if (jacobian_quantities[q].Mode() ==
"nd") {
631 for (
Index i3 = 0; i3 < vmr_field.
ncols(); i3++) {
632 for (
Index i2 = 0; i2 < vmr_field.
nrows(); i2++) {
633 for (
Index i1 = 0; i1 < vmr_field.
npages(); i1++) {
634 vmr_field(isp, i1, i2, i3) =
635 x_field(i1, i2, i3) /
640 }
else if (jacobian_quantities[q].Mode() ==
"rh") {
644 for (
Index i3 = 0; i3 < vmr_field.
ncols(); i3++) {
645 for (
Index i2 = 0; i2 < vmr_field.
nrows(); i2++) {
646 for (
Index i1 = 0; i1 < vmr_field.
npages(); i1++) {
647 vmr_field(isp, i1, i2, i3) =
648 x_field(i1, i2, i3) * water_p_eq(i1, i2, i3) / p_grid[i1];
652 }
else if (jacobian_quantities[q].Mode() ==
"q") {
657 for (
Index i3 = 0; i3 < vmr_field.
ncols(); i3++) {
658 for (
Index i2 = 0; i2 < vmr_field.
nrows(); i2++) {
659 for (
Index i1 = 0; i1 < vmr_field.
npages(); i1++) {
660 const Numeric r = x_field(i1, i2, i3) / (1 - x_field(i1, i2, i3));
661 const Numeric e = r * p_grid[i1] / (0.622 + r);
662 vmr_field(isp, i1, i2, i3) = e / p_grid[i1];
673 else if (jacobian_quantities[q] == Jacobian::Special::ScatteringString) {
677 "One jacobian quantity belongs to the "
678 "scattering species category, but *particle_bulkprop_field* "
681 particle_bulkprop_names.
nelem(),
682 "Mismatch in size between "
683 "*particle_bulkprop_field* and *particle_bulkprop_field*.");
686 jacobian_quantities[q].SubSubtag());
688 "Jacobian quantity with index ", q,
" covers a "
689 "scattering species, and the field quantity is set to \"",
690 jacobian_quantities[q].SubSubtag(),
"\", but this quantity "
691 "could not found in *particle_bulkprop_names*.")
696 Index n_p, n_lat, n_lon;
703 jacobian_quantities[q].Grids(),
709 Tensor3 pbfield_x(n_p, n_lat, n_lon);
713 pbfield, atmosphere_dim, pbfield_x, gp_p, gp_lat, gp_lon);
720 else if (jacobian_quantities[q].Target().isWind()) {
724 Index n_p, n_lat, n_lon;
731 jacobian_quantities[q].Grids(),
738 Tensor3 wind_x(n_p, n_lat, n_lon);
746 wind_field, atmosphere_dim, wind_x, gp_p, gp_lat, gp_lon);
748 if (jacobian_quantities[q] == Jacobian::Atm::WindU) {
749 wind_u_field = wind_field;
750 }
else if (jacobian_quantities[q] == Jacobian::Atm::WindV) {
751 wind_v_field = wind_field;
752 }
else if (jacobian_quantities[q] == Jacobian::Atm::WindW) {
753 wind_w_field = wind_field;
759 else if (jacobian_quantities[q].Target().isMagnetic()) {
763 Index n_p, n_lat, n_lon;
770 jacobian_quantities[q].Grids(),
777 Tensor3 mag_x(n_p, n_lat, n_lon);
785 mag_field, atmosphere_dim, mag_x, gp_p, gp_lat, gp_lon);
786 if (jacobian_quantities[q] == Jacobian::Atm::MagneticU) {
787 mag_u_field = mag_field;
788 }
else if (jacobian_quantities[q] == Jacobian::Atm::MagneticV) {
789 mag_v_field = mag_field;
790 }
else if (jacobian_quantities[q] == Jacobian::Atm::MagneticW) {
791 mag_w_field = mag_field;
792 }
else if (jacobian_quantities[q] == Jacobian::Atm::MagneticMagnitude) {
793 for (
Index i = 0; i < n_p; i++) {
794 for (
Index j = 0; j < n_lat; j++) {
795 for (
Index k = 0; k < n_lon; k++) {
796 const Numeric scale = mag_x(i, j, k) /
797 std::hypot(mag_u_field(i, j, k), mag_v_field(i, j, k), mag_w_field(i, j, k));
798 mag_u_field(i, j, k) *= scale;
799 mag_v_field(i, j, k) *= scale;
800 mag_w_field(i, j, k) *= scale;
811 else if (jacobian_quantities[q] == Jacobian::Special::SurfaceString) {
816 surface_props_names);
818 "One jacobian quantity belongs to the "
819 "surface category, but *surface_props_data* is empty.");
822 find_first(surface_props_names, jacobian_quantities[q].Subtag());
824 "Jacobian quantity with index ", q,
" covers a "
825 "surface property, and the field Subtag is set to \"",
826 jacobian_quantities[q].Subtag(),
"\", but this quantity "
827 "could not found in *surface_props_names*.")
837 jacobian_quantities[q].Grids(),
842 Matrix surf_x(n_lat, n_lon);
857 Vector& sensor_response_f,
859 Matrix& sensor_response_dlos,
860 Vector& sensor_response_f_grid,
862 Matrix& sensor_response_dlos_grid,
866 const Agenda& sensor_response_agenda,
867 const Index& sensor_checked,
873 "The sensor response must be flagged to have "
874 "passed a consistency check (sensor_checked=1).");
892 "Length of *x* does not match length implied by "
893 "*jacobian_quantities*.");
899 bool do_sensor =
false;
902 for (
Index q = 0; q < nq; q++) {
904 const Index np = ji[q][1] - ji[q][0] + 1;
908 if (jacobian_quantities[q].Target().isPointing()) {
910 if (jacobian_quantities[q].Grids()[0][0] == -1) {
912 "Mismatch between pointing jacobian and *sensor_los* found.");
914 for (
Index i = 0; i < np; i++) {
915 sensor_los(i, 0) += x_t[ji[q][0] + i];
921 "Sizes of *sensor_los* and *sensor_time* do not match.");
925 for (
Index i = 0; i <
w.nelem(); i++) {
926 sensor_los(i, 0) +=
w[i] * x_t[ji[q][0] +
c];
934 else if (jacobian_quantities[q].Target().isFrequency()) {
935 if (jacobian_quantities[q] == Jacobian::Sensor::FrequencyShift) {
937 if (x_t[ji[q][0]] != 0) {
939 f_backend += x_t[ji[q][0]];
941 }
else if (jacobian_quantities[q] == Jacobian::Sensor::FrequencyStretch) {
943 if (x_t[ji[q][0]] != 0) {
947 for (
Index i = 0; i <
w.nelem(); i++) {
948 f_backend[i] +=
w[i] * x_t[ji[q][0]];
958 else if (jacobian_quantities[q] == Jacobian::Sensor::Polyfit ||
959 jacobian_quantities[q] == Jacobian::Sensor::Sinefit) {
963 sensor_response_pol_grid.
nelem() *
964 sensor_response_dlos_grid.
nrows();
965 y_baseline.
resize(y_size);
969 for (
Index mb = 0; mb < sensor_los.
nrows(); ++mb) {
974 sensor_response_pol_grid,
975 sensor_response_f_grid,
976 sensor_response_dlos_grid,
977 jacobian_quantities[q],
995 sensor_response_f_grid,
997 sensor_response_pol_grid,
998 sensor_response_dlos,
999 sensor_response_dlos_grid,
1002 sensor_response_agenda);
1008 ARTS_USER_ERROR (
"Retrievals of spectroscopic variables not yet handled.");
1027 const Agenda& inversion_iterate_agenda,
1029 const Numeric& max_start_cost,
1031 const Index& max_iter,
1033 const Vector& lm_ga_settings,
1034 const Index& clear_matrices,
1035 const Index& display_progress,
1049 inversion_iterate_agenda,
1054 jacobian_quantities,
1064 oem_diagnostics.
resize(5);
1065 oem_diagnostics = NAN;
1067 if (method ==
"ml" || method ==
"lm" || method ==
"ml_cg" ||
1068 method ==
"lm_cg") {
1069 lm_ga_history.
resize(max_iter + 1);
1070 lm_ga_history = NAN;
1076 if (x.
nelem() != n) {
1084 if (yf.
nelem() == 0) {
1086 ws, yf, jacobian, xa, 1, 0, inversion_iterate_agenda);
1090 "Mismatch between simulated y and input y.\n"
1091 "Input y is size ", y.
nelem(),
" but simulated y is ",
1093 "Use your frequency grid vector and your sensor response matrix to match simulations with measurements.\n")
1098 if (method == "ml" || method == "lm" || display_progress ||
1099 max_start_cost > 0) {
1108 cost_start = dx * sdx + dy * sdy;
1109 cost_start /=
static_cast<Numeric>(m);
1111 oem_diagnostics[1] = cost_start;
1114 if (max_start_cost > 0 && cost_start > max_start_cost) {
1116 oem_diagnostics[0] = 99;
1118 if (display_progress) {
1119 cout <<
"\n No OEM inversion, too high start cost:\n"
1120 <<
" Set limit : " << max_start_cost << endl
1121 <<
" Found value : " << cost_start << endl
1127 bool apply_norm =
false;
1129 if (x_norm.
nelem() == n) {
1132 T.diagonal() = x_norm;
1133 for (
Index i = 0; i < n; i++) {
1134 T(i, i) = x_norm[i];
1146 &inversion_iterate_agenda);
1149 int oem_verbosity =
static_cast<int>(display_progress);
1151 int return_code = 0;
1154 if (method ==
"li") {
1158 x_oem, y_oem, gn, oem_verbosity, lm_ga_history,
true);
1159 oem_diagnostics[0] =
static_cast<Index>(return_code);
1160 }
else if (method ==
"li_m") {
1164 x_oem, y_oem, gn, oem_verbosity, lm_ga_history,
true);
1165 oem_diagnostics[0] =
static_cast<Index>(return_code);
1166 }
else if (method ==
"li_cg") {
1167 oem::CG cg(T, apply_norm, 1e-10, 0);
1170 x_oem, y_oem, gn, oem_verbosity, lm_ga_history,
true);
1171 oem_diagnostics[0] =
static_cast<Index>(return_code);
1172 }
else if (method ==
"li_cg_m") {
1173 oem::CG cg(T, apply_norm, 1e-10, 0);
1176 x_oem, y_oem, gn, oem_verbosity, lm_ga_history,
true);
1177 oem_diagnostics[0] =
static_cast<Index>(return_code);
1178 }
else if (method ==
"gn") {
1180 oem::GN gn(stop_dx, (
unsigned int)max_iter, s);
1182 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1183 oem_diagnostics[0] =
static_cast<Index>(return_code);
1184 }
else if (method ==
"gn_m") {
1186 oem::GN gn(stop_dx, (
unsigned int)max_iter, s);
1188 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1189 oem_diagnostics[0] =
static_cast<Index>(return_code);
1190 }
else if (method ==
"gn_cg") {
1191 oem::CG cg(T, apply_norm, 1e-10, 0);
1192 oem::GN_CG gn(stop_dx, (
unsigned int)max_iter, cg);
1194 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1195 oem_diagnostics[0] =
static_cast<Index>(return_code);
1196 }
else if (method ==
"gn_cg_m") {
1197 oem::CG cg(T, apply_norm, 1e-10, 0);
1198 oem::GN_CG gn(stop_dx, (
unsigned int)max_iter, cg);
1200 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1201 oem_diagnostics[0] =
static_cast<Index>(return_code);
1202 }
else if ((method ==
"lm") || (method ==
"ml")) {
1208 std::make_pair(0, 0),
1209 make_shared<Sparse>(diagonal)));
1213 lm.set_tolerance(stop_dx);
1214 lm.set_maximum_iterations((
unsigned int)max_iter);
1215 lm.set_lambda(lm_ga_settings[0]);
1216 lm.set_lambda_decrease(lm_ga_settings[1]);
1217 lm.set_lambda_increase(lm_ga_settings[2]);
1218 lm.set_lambda_maximum(lm_ga_settings[3]);
1219 lm.set_lambda_threshold(lm_ga_settings[4]);
1220 lm.set_lambda_constraint(lm_ga_settings[5]);
1223 x_oem, y_oem, lm, oem_verbosity, lm_ga_history);
1224 oem_diagnostics[0] =
static_cast<Index>(return_code);
1225 if (lm.get_lambda() > lm.get_lambda_maximum()) {
1226 oem_diagnostics[0] = 2;
1228 }
else if ((method ==
"lm_cg") || (method ==
"ml_cg")) {
1229 oem::CG cg(T, apply_norm, 1e-10, 0);
1235 std::make_pair(0, 0),
1236 make_shared<Sparse>(diagonal)));
1239 lm.set_maximum_iterations((
unsigned int)max_iter);
1240 lm.set_lambda(lm_ga_settings[0]);
1241 lm.set_lambda_decrease(lm_ga_settings[1]);
1242 lm.set_lambda_increase(lm_ga_settings[2]);
1243 lm.set_lambda_threshold(lm_ga_settings[3]);
1244 lm.set_lambda_maximum(lm_ga_settings[4]);
1247 x_oem, y_oem, lm, oem_verbosity, lm_ga_history);
1248 oem_diagnostics[0] =
static_cast<Index>(return_code);
1249 if (lm.get_lambda() > lm.get_lambda_maximum()) {
1250 oem_diagnostics[0] = 2;
1254 oem_diagnostics[2] =
oem.cost /
static_cast<Numeric>(m);
1255 oem_diagnostics[3] =
oem.cost_y /
static_cast<Numeric>(m);
1256 oem_diagnostics[4] =
static_cast<Numeric>(
oem.iterations);
1257 }
catch (
const std::exception& e) {
1258 oem_diagnostics[0] = 9;
1259 oem_diagnostics[2] =
oem.cost;
1260 oem_diagnostics[3] =
oem.cost_y;
1261 oem_diagnostics[4] =
static_cast<Numeric>(
oem.iterations);
1264 for (
auto& s : sv) {
1265 std::stringstream ss{s};
1267 while (std::getline(ss, t)) {
1268 errors.push_back(t.c_str());
1276 yf = aw.get_measurement_vector();
1279 if (clear_matrices) {
1282 }
else if (oem_diagnostics[0] <= 2) {
1284 Matrix tmp1(n, m), tmp2(n, n), tmp3(n, n);
1286 mult(tmp2, tmp1, jacobian);
1289 mult(dxdy, tmp3, tmp1);
1303 "The gain matrix *dxdy* is required to compute the observation error covariance matrix.");
1305 "The covariance matrix covmat_se has invalid dimensions.");
1309 mult(covmat_so, dxdy, tmp1);
1318 Matrix tmp1(n, n), tmp2(n, n);
1321 "The averaging kernel matrix *dxdy* is required to compute the smoothing error covariance matrix.");
1323 "The covariance matrix *covmat_sx* invalid dimensions.");
1332 mult(covmat_ss, tmp1, tmp2);
1349 "The Jacobian matrix is empty.");
1351 "Matrices have inconsistent sizes.\n"
1352 " Size of gain matrix: ", dxdy.
nrows(),
" x ", dxdy.
ncols(),
1354 " Size of Jacobian: ", jacobian.
nrows(),
" x ",
1355 jacobian.
ncols(),
"\n")
1358 mult(avk, dxdy, jacobian);
1368 "WSM is not available because ARTS was compiled without "
1377 "WSM is not available because ARTS was compiled without "
1386 "WSM is not available because ARTS was compiled without "
1416 "WSM is not available because ARTS was compiled without "
1421#if defined(OEM_SUPPORT) && 0
1423#include "agenda_wrapper_mpi.h"
1433void MPI_Initialize(
Matrix& sensor_los,
1438 MPI_Initialized(&initialized);
1440 MPI_Init(
nullptr,
nullptr);
1445 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1446 MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
1448 int nmblock = (int)sensor_pos.
nrows();
1449 int mblock_range = nmblock / nprocs;
1450 int mblock_start = mblock_range * rank;
1451 int remainder = nmblock % std::max(mblock_range, nprocs);
1457 if (rank < remainder) {
1459 mblock_start += rank;
1461 mblock_start += remainder;
1465 Range range =
Range(mblock_start, mblock_range);
1470 tmp_m = sensor_pos(range,
joker);
1473 Vector tmp_v = sensor_time[range];
1474 sensor_time = tmp_v;
1497 const Agenda& inversion_iterate_agenda,
1499 const Numeric& max_start_cost,
1501 const Index& max_iter,
1503 const Vector& lm_ga_settings,
1504 const Index& clear_matrices,
1505 const Index& display_progress,
1516 inversion_iterate_agenda,
1520 jacobian_quantities,
1535 MPI_Initialize(sensor_los, sensor_pos, sensor_time);
1545 ws, tmp, jacobian, xa, 1, inversion_iterate_agenda);
1549 oem_diagnostics.
resize(5);
1550 oem_diagnostics = NAN;
1552 if (method ==
"ml" || method ==
"lm") {
1553 lm_ga_history.
resize(max_iter);
1554 lm_ga_history = NAN;
1562 if (method ==
"ml" || method ==
"lm" || display_progress ||
1563 max_start_cost > 0) {
1566 cost_start = dot(dy, SeInvMPI * dy);
1568 oem_diagnostics[1] = cost_start;
1571 if (max_start_cost > 0 && cost_start > max_start_cost) {
1573 oem_diagnostics[0] = 99;
1575 if (display_progress) {
1576 cout <<
"\n No OEM inversion, too high start cost:\n"
1577 <<
" Set limit : " << max_start_cost << endl
1578 <<
" Found value : " << cost_start << endl
1589 OEMVector xa_oem(xa), y_oem(y), x_oem;
1592 OEM_PS_PS_MPI<AgendaWrapperMPI>
oem(aw, xa_oem, SaInvMPI, SeInvMPI);
1595 int return_value = 99;
1597 if (method ==
"li") {
1599 oem::GN_CG gn(stop_dx, (
unsigned int)max_iter, cg);
1601 x_oem, y_oem, gn, 2 * (int)display_progress);
1602 }
else if (method ==
"gn") {
1604 oem::GN_CG gn(stop_dx, (
unsigned int)max_iter, cg);
1606 x_oem, y_oem, gn, 2 * (int)display_progress);
1607 }
else if ((method ==
"lm") || (method ==
"ml")) {
1609 LM_CG_MPI lm(SaInvMPI, cg);
1611 lm.set_tolerance(stop_dx);
1612 lm.set_maximum_iterations((
unsigned int)max_iter);
1613 lm.set_lambda(lm_ga_settings[0]);
1614 lm.set_lambda_decrease(lm_ga_settings[1]);
1615 lm.set_lambda_increase(lm_ga_settings[2]);
1616 lm.set_lambda_threshold(lm_ga_settings[3]);
1617 lm.set_lambda_maximum(lm_ga_settings[4]);
1619 return_value =
oem.compute<oem::LM_CG_MPI, invlib::MPILog>(
1620 x_oem, y_oem, lm, 2 * (int)display_progress);
1623 oem_diagnostics[0] = return_value;
1624 oem_diagnostics[2] =
oem.cost;
1625 oem_diagnostics[3] =
oem.cost_y;
1626 oem_diagnostics[4] =
static_cast<Numeric>(
oem.iterations);
1630 if (clear_matrices && (oem_diagnostics[0])) {
1667 "You have to compile ARTS with OEM support "
1668 " 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 &out, const CovarianceMatrix &in, const Verbosity &verbosity)
WORKSPACE METHOD: MatrixFromCovarianceMatrix.
Index nelem() const ARTS_NOEXCEPT
Index nrows() const noexcept
Index ncols() const noexcept
A constant view of a Tensor3.
bool empty() const noexcept
Index npages() const
Returns the number of pages.
Index nrows() const
Returns the number of rows.
Index ncols() const
Returns the number of columns.
Index ncols() const noexcept
bool empty() const noexcept
Index nrows() const noexcept
Index nbooks() const noexcept
Index npages() const noexcept
Index nelem() const noexcept
Returns the number of elements.
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.
void resize(Index r, Index c)
Resize function.
void resize(Index n)
Resize function.
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)
#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 id_mat(MatrixView I)
Identity Matrix.
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
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
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
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.
Vector diagonal() const
Diagonal elements as vector.
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