ARTS 2.5.11 (git: 725533f0)
m_oem.cc
Go to the documentation of this file.
1/*===========================================================================
2 === File description
3 ===========================================================================*/
4
16/*===========================================================================
17 === External declarations
18 ===========================================================================*/
19
20#include <cmath>
21#include <iterator>
22#include <sstream>
23#include <stdexcept>
24#include <string>
25#include "array.h"
26#include "arts.h"
27#include "arts_omp.h"
28#include "auto_md.h"
29#include "debug.h"
30#include "jacobian.h"
31#include "math_funcs.h"
32#include "physics_funcs.h"
33#include "rte.h"
34#include "special_interp.h"
35#include "surface.h"
36#include "check_input.h"
37
38#pragma GCC diagnostic ignored "-Wconversion"
39
40#ifdef OEM_SUPPORT
41#include "oem.h"
42#endif
43
44/* Workspace method: Doxygen documentation will be auto-generated */
45void particle_bulkprop_fieldClip(Tensor4& particle_bulkprop_field,
46 const ArrayOfString& particle_bulkprop_names,
47 const String& bulkprop_name,
48 const Numeric& limit_low,
49 const Numeric& limit_high,
50 const Verbosity&) {
51 Index iq = -1;
52 if (bulkprop_name == "ALL") {
53 }
54
55 else {
56 for (Index i = 0; i < particle_bulkprop_names.nelem(); i++) {
57 if (particle_bulkprop_names[i] == bulkprop_name) {
58 iq = i;
59 break;
60 }
61 }
62 ARTS_USER_ERROR_IF (iq < 0,
63 "Could not find ", bulkprop_name,
64 " in particle_bulkprop_names.\n")
65 }
66
67 Tensor4Clip(particle_bulkprop_field, iq, limit_low, limit_high);
68}
69
70/* Workspace method: Doxygen documentation will be auto-generated */
71void vmr_fieldClip(Tensor4& vmr_field,
72 const ArrayOfArrayOfSpeciesTag& abs_species,
73 const String& species,
74 const Numeric& limit_low,
75 const Numeric& limit_high,
76 const Verbosity&) {
77 Index iq = -1;
78 if (species == "ALL") {
79 }
80
81 else {
82 for (Index i = 0; i < abs_species.nelem(); i++) {
83 if (abs_species[i].Species() == SpeciesTag(species).Spec()) {
84 iq = i;
85 break;
86 }
87 }
88 ARTS_USER_ERROR_IF (iq < 0,
89 "Could not find ", species, " in abs_species.\n")
90 }
91
92 Tensor4Clip(vmr_field, iq, limit_low, limit_high);
93}
94
95/* Workspace method: Doxygen documentation will be auto-generated */
96void xClip(Vector& x,
97 const ArrayOfRetrievalQuantity& jacobian_quantities,
98 const Index& ijq,
99 const Numeric& limit_low,
100 const Numeric& limit_high,
101 const Verbosity&) {
102 // Sizes
103 const Index nq = jacobian_quantities.nelem();
104
105 ARTS_USER_ERROR_IF (ijq < -1, "Argument *ijq* must be >= -1.");
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")
111
112 // Jacobian indices
113 ArrayOfArrayOfIndex ji;
114 {
115 bool any_affine;
116 jac_ranges_indices(ji, any_affine, jacobian_quantities);
117 }
118
119 Index ifirst = 0, ilast = x.nelem() - 1;
120 if (ijq > -1) {
121 ifirst = ji[ijq][0];
122 ilast = ji[ijq][1];
123 }
124
125 if (!std::isinf(limit_low)) {
126 for (Index i = ifirst; i <= ilast; i++) {
127 if (x[i] < limit_low) x[i] = limit_low;
128 }
129 }
130 if (!std::isinf(limit_high)) {
131 for (Index i = ifirst; i <= ilast; i++) {
132 if (x[i] > limit_high) x[i] = limit_high;
133 }
134 }
135}
136
137/* Workspace method: Doxygen documentation will be auto-generated */
138void xaStandard(Workspace& ws,
139 Vector& xa,
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,
163 const Verbosity&) {
164 // Basics
165 //
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).");
175
176 // Jacobian indices
177 ArrayOfArrayOfIndex ji;
178 {
179 bool any_affine;
180 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
181 }
182
183 // Sizes
184 const Index nq = jacobian_quantities.nelem();
185 //
186 xa.resize(ji[nq - 1][1] + 1);
187
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);
193
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,
199 gp_lat,
200 gp_lon,
201 jacobian_quantities[q],
202 atmosphere_dim,
203 p_grid,
204 lat_grid,
205 lon_grid);
206 Tensor3 t_x;
207 regrid_atmfield_by_gp(t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
208 flat(xa[ind], t_x);
209 }
210
211 // Abs species
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);
216
217 if (jacobian_quantities[q].Mode() == "rel") {
218 // This one is simple, just a vector of ones
219 xa[ind] = 1;
220 } else {
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,
224 gp_lat,
225 gp_lon,
226 jacobian_quantities[q],
227 atmosphere_dim,
228 p_grid,
229 lat_grid,
230 lon_grid);
231 Tensor3 vmr_x;
232 regrid_atmfield_by_gp(vmr_x,
233 atmosphere_dim,
234 vmr_field(isp, joker, joker, joker),
235 gp_p,
236 gp_lat,
237 gp_lon);
238
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*
243 Tensor3 t_x;
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)
247 Index i = 0;
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++) {
251 xa[ji[q][0] + i] =
252 vmr_x(i1, i2, i3) *
253 number_density(jacobian_quantities[q].Grids()[0][i1],
254 t_x(i1, i2, i3));
255 i += 1;
256 }
257 }
258 }
259 } else if (jacobian_quantities[q].Mode() == "rh") {
260 // Here we need to also interpolate *t_field*
261 Tensor3 t_x;
262 regrid_atmfield_by_gp(
263 t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
264 Tensor3 water_p_eq;
265 water_p_eq_agendaExecute(ws, water_p_eq, t_x, water_p_eq_agenda);
266 // Calculate relative humidity (vmr*p/p_sat)
267 Index i = 0;
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);
274 i += 1;
275 }
276 }
277 }
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;
282 Index i = 0;
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++) {
286 const Numeric e =
287 vmr_x(i1, i2, i3) * jacobian_quantities[q].Grids()[0][i1];
288 const Numeric r =
289 0.622 * e / (jacobian_quantities[q].Grids()[0][i1] - e);
290 xa[ji[q][0] + i] = r / (1 + r);
291 i += 1;
292 }
293 }
294 }
295 } else {
296 ARTS_ASSERT(0);
297 }
298 }
299 }
300
301 // Scattering species
302 else if (jacobian_quantities[q] == Jacobian::Special::ScatteringString) {
303 if (cloudbox_on) {
304 ARTS_USER_ERROR_IF (particle_bulkprop_field.empty(),
305 "One jacobian quantity belongs to the "
306 "scattering species category, but *particle_bulkprop_field* "
307 "is empty.");
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*.");
312
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*.")
320
321 ArrayOfGridPos gp_p, gp_lat, gp_lon;
322 get_gp_atmgrids_to_rq(gp_p,
323 gp_lat,
324 gp_lon,
325 jacobian_quantities[q],
326 atmosphere_dim,
327 p_grid,
328 lat_grid,
329 lon_grid);
330 Tensor3 pbp_x;
331 regrid_atmfield_by_gp(pbp_x,
332 atmosphere_dim,
333 particle_bulkprop_field(isp, joker, joker, joker),
334 gp_p,
335 gp_lat,
336 gp_lon);
337 flat(xa[ind], pbp_x);
338 } else {
339 xa[ind] = 0;
340 }
341 }
342
343 // Wind
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);
350 }
351
352 // Determine grid positions for interpolation from retrieval grids back
353 // to atmospheric grids
354 ArrayOfGridPos gp_p, gp_lat, gp_lon;
355 get_gp_atmgrids_to_rq(gp_p,
356 gp_lat,
357 gp_lon,
358 jacobian_quantities[q],
359 atmosphere_dim,
360 p_grid,
361 lat_grid,
362 lon_grid);
363
364 Tensor3 wind_x;
365 regrid_atmfield_by_gp(
366 wind_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
367 flat(xa[ind], wind_x);
368 }
369
370 // Magnetism
371 else if (jacobian_quantities[q].Target().isMagnetic()) {
372 if (jacobian_quantities[q] == Jacobian::Atm::MagneticMagnitude) {
373 // Determine grid positions for interpolation from retrieval grids back
374 // to atmospheric grids
375 ArrayOfGridPos gp_p, gp_lat, gp_lon;
376 get_gp_atmgrids_to_rq(gp_p,
377 gp_lat,
378 gp_lon,
379 jacobian_quantities[q],
380 atmosphere_dim,
381 p_grid,
382 lat_grid,
383 lon_grid);
384
385 //all three component's hyoptenuse is the strength
386 Tensor3 mag_u, mag_v, mag_w;
387 regrid_atmfield_by_gp(
388 mag_u, atmosphere_dim, mag_u_field, gp_p, gp_lat, gp_lon);
389 regrid_atmfield_by_gp(
390 mag_v, atmosphere_dim, mag_v_field, gp_p, gp_lat, gp_lon);
391 regrid_atmfield_by_gp(
392 mag_w, atmosphere_dim, mag_w_field, gp_p, gp_lat, gp_lon);
393
394 Tensor3 mag_x(gp_p.nelem(), gp_lat.nelem(), gp_lon.nelem());
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);
400 } else {
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);
406 }
407
408 // Determine grid positions for interpolation from retrieval grids back
409 // to atmospheric grids
410 ArrayOfGridPos gp_p, gp_lat, gp_lon;
411 get_gp_atmgrids_to_rq(gp_p,
412 gp_lat,
413 gp_lon,
414 jacobian_quantities[q],
415 atmosphere_dim,
416 p_grid,
417 lat_grid,
418 lon_grid);
419
420 Tensor3 mag_x;
421 regrid_atmfield_by_gp(
422 mag_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
423 flat(xa[ind], mag_x);
424 }
425 }
426
427 // Surface
428 else if (jacobian_quantities[q] == Jacobian::Special::SurfaceString) {
429 surface_props_check(atmosphere_dim,
430 lat_grid,
431 lon_grid,
432 surface_props_data,
433 surface_props_names);
434 ARTS_USER_ERROR_IF (surface_props_data.empty(),
435 "One jacobian quantity belongs to the "
436 "surface category, but *surface_props_data* is empty.");
437
438 const Index isu =
439 find_first(surface_props_names, jacobian_quantities[q].Subtag());
440 ARTS_USER_ERROR_IF (isu < 0,
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*.")
445
446 ArrayOfGridPos gp_lat, gp_lon;
447 get_gp_atmsurf_to_rq(gp_lat,
448 gp_lon,
449 jacobian_quantities[q],
450 atmosphere_dim,
451 lat_grid,
452 lon_grid);
453 Matrix surf_x;
454 regrid_atmsurf_by_gp_oem(surf_x,
455 atmosphere_dim,
456 surface_props_data(isu, joker, joker),
457 gp_lat,
458 gp_lon);
459 flat(xa[ind], surf_x);
460 }
461
462 // All variables having zero as a priori
463 // ----------------------------------------------------------------------------
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) {
468 xa[ind] = 0;
469 }
470
471 else {
472 ARTS_USER_ERROR (
473 "Found a retrieval quantity that is not yet handled by\n"
474 "internal retrievals: ", jacobian_quantities[q], '\n')
475 }
476 }
477
478 // Apply transformations
479 transform_x(xa, jacobian_quantities);
480}
481
482/* Workspace method: Doxygen documentation will be auto-generated */
483void x2artsAtmAndSurf(Workspace& ws,
484 Tensor4& vmr_field,
485 Tensor3& t_field,
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,
494 const ArrayOfRetrievalQuantity& jacobian_quantities,
495 const Vector& x,
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,
502 const ArrayOfArrayOfSpeciesTag& abs_species,
503 const Index& cloudbox_on,
504 const Index& cloudbox_checked,
505 const ArrayOfString& particle_bulkprop_names,
506 const ArrayOfString& surface_props_names,
507 const Agenda& water_p_eq_agenda,
508 const Verbosity&) {
509 // Basics
510 //
511 ARTS_USER_ERROR_IF (atmfields_checked != 1,
512 "The atmospheric fields must be flagged to have "
513 "passed a consistency check (atmfields_checked=1).");
514 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
515 "The atmospheric geometry must be flagged to have "
516 "passed a consistency check (atmgeom_checked=1).");
517 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
518 "The cloudbox must be flagged to have "
519 "passed a consistency check (cloudbox_checked=1).");
520
521 // Revert transformation
522 Vector x_t(x);
523 transform_x_back(x_t, jacobian_quantities);
524
525 // Main sizes
526 const Index nq = jacobian_quantities.nelem();
527
528 // Jacobian indices
529 ArrayOfArrayOfIndex ji;
530 {
531 bool any_affine;
532 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
533 }
534
535 // Check input
536 ARTS_USER_ERROR_IF (x_t.nelem() != ji[nq - 1][1] + 1,
537 "Length of *x* does not match length implied by "
538 "*jacobian_quantities*.");
539
540 // Note that when this method is called, vmr_field and other output variables
541 // have original values, i.e. matching the a priori state.
542
543 // Loop retrieval quantities
544 for (Index q = 0; q < nq; q++) {
545 // Index range of this retrieval quantity
546 const Index np = ji[q][1] - ji[q][0] + 1;
547 Range ind(ji[q][0], np);
548
549 // Atmospheric temperatures
550 // ----------------------------------------------------------------------------
551 if (jacobian_quantities[q] == Jacobian::Atm::Temperature) {
552 // Determine grid positions for interpolation from retrieval grids back
553 // to atmospheric grids
554 ArrayOfGridPos gp_p, gp_lat, gp_lon;
555 Index n_p, n_lat, n_lon;
556 get_gp_rq_to_atmgrids(gp_p,
557 gp_lat,
558 gp_lon,
559 n_p,
560 n_lat,
561 n_lon,
562 jacobian_quantities[q].Grids(),
563 atmosphere_dim,
564 p_grid,
565 lat_grid,
566 lon_grid);
567
568 // Map values in x back to t_field
569 Tensor3 t_x(n_p, n_lat, n_lon);
570 reshape(t_x, x_t[ind]);
571 regrid_atmfield_by_gp_oem(
572 t_field, atmosphere_dim, t_x, gp_p, gp_lat, gp_lon);
573 }
574
575 // Abs species
576 // ----------------------------------------------------------------------------
577 else if (jacobian_quantities[q].Target().isSpeciesVMR()) {
578 // Index position of species
579 ArrayOfSpeciesTag atag(jacobian_quantities[q].Subtag());
580 const Index isp = chk_contains("abs_species", abs_species, atag);
581
582 // Map part of x to a full atmospheric field
583 Tensor3 x_field(vmr_field.npages(), vmr_field.nrows(), vmr_field.ncols());
584 {
585 ArrayOfGridPos gp_p, gp_lat, gp_lon;
586 Index n_p, n_lat, n_lon;
587 get_gp_rq_to_atmgrids(gp_p,
588 gp_lat,
589 gp_lon,
590 n_p,
591 n_lat,
592 n_lon,
593 jacobian_quantities[q].Grids(),
594 atmosphere_dim,
595 p_grid,
596 lat_grid,
597 lon_grid);
598 //
599 Tensor3 t3_x(n_p, n_lat, n_lon);
600 reshape(t3_x, x_t[ind]);
601 regrid_atmfield_by_gp_oem(
602 x_field, atmosphere_dim, t3_x, gp_p, gp_lat, gp_lon);
603 }
604 //
605 if (jacobian_quantities[q].Mode() == "rel") {
606 // vmr = vmr0 * x
607 vmr_field(isp, joker, joker, joker) *= x_field;
608 } else if (jacobian_quantities[q].Mode() == "vmr") {
609 // vmr = x
610 vmr_field(isp, joker, joker, joker) = x_field;
611 } else if (jacobian_quantities[q].Mode() == "nd") {
612 // vmr = nd / nd_tot
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) /
618 number_density(p_grid[i1], t_field(i1, i2, i3));
619 }
620 }
621 }
622 } else if (jacobian_quantities[q].Mode() == "rh") {
623 // vmr = x * p_sat / p
624 Tensor3 water_p_eq;
625 water_p_eq_agendaExecute(ws, water_p_eq, t_field, water_p_eq_agenda);
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];
631 }
632 }
633 }
634 } else if (jacobian_quantities[q].Mode() == "q") {
635 // We have that specific humidity q, mixing ratio r and
636 // vapour pressure e, are related as
637 // q = r(1+r); r = 0.622e/(p-e); e = vmr*p;
638 // That is: vmr=e/p; e = rp/(0.622+r); r = q/(1-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];
645 }
646 }
647 }
648 } else {
649 ARTS_ASSERT(0);
650 }
651 }
652
653 // Scattering species
654 // ----------------------------------------------------------------------------
655 else if (jacobian_quantities[q] == Jacobian::Special::ScatteringString) {
656 // If no cloudbox, we assume that there is nothing to do
657 if (cloudbox_on) {
658 ARTS_USER_ERROR_IF (particle_bulkprop_field.empty(),
659 "One jacobian quantity belongs to the "
660 "scattering species category, but *particle_bulkprop_field* "
661 "is empty.");
662 ARTS_USER_ERROR_IF (particle_bulkprop_field.nbooks() !=
663 particle_bulkprop_names.nelem(),
664 "Mismatch in size between "
665 "*particle_bulkprop_field* and *particle_bulkprop_field*.");
666
667 const Index isp = find_first(particle_bulkprop_names,
668 jacobian_quantities[q].SubSubtag());
669 ARTS_USER_ERROR_IF (isp < 0,
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*.")
674
675 // Determine grid positions for interpolation from retrieval grids back
676 // to atmospheric grids
677 ArrayOfGridPos gp_p, gp_lat, gp_lon;
678 Index n_p, n_lat, n_lon;
679 get_gp_rq_to_atmgrids(gp_p,
680 gp_lat,
681 gp_lon,
682 n_p,
683 n_lat,
684 n_lon,
685 jacobian_quantities[q].Grids(),
686 atmosphere_dim,
687 p_grid,
688 lat_grid,
689 lon_grid);
690 // Map x to particle_bulkprop_field
691 Tensor3 pbfield_x(n_p, n_lat, n_lon);
692 reshape(pbfield_x, x_t[ind]);
693 Tensor3 pbfield;
694 regrid_atmfield_by_gp_oem(
695 pbfield, atmosphere_dim, pbfield_x, gp_p, gp_lat, gp_lon);
696 particle_bulkprop_field(isp, joker, joker, joker) = pbfield;
697 }
698 }
699
700 // Wind
701 // ----------------------------------------------------------------------------
702 else if (jacobian_quantities[q].Target().isWind()) {
703 // Determine grid positions for interpolation from retrieval grids back
704 // to atmospheric grids
705 ArrayOfGridPos gp_p, gp_lat, gp_lon;
706 Index n_p, n_lat, n_lon;
707 get_gp_rq_to_atmgrids(gp_p,
708 gp_lat,
709 gp_lon,
710 n_p,
711 n_lat,
712 n_lon,
713 jacobian_quantities[q].Grids(),
714 atmosphere_dim,
715 p_grid,
716 lat_grid,
717 lon_grid);
718
719 // TODO Could be done without copying.
720 Tensor3 wind_x(n_p, n_lat, n_lon);
721 reshape(wind_x, x_t[ind]);
722
723 Tensor3View target_field(wind_u_field);
724
725 Tensor3 wind_field(
726 target_field.npages(), target_field.nrows(), target_field.ncols());
727 regrid_atmfield_by_gp_oem(
728 wind_field, atmosphere_dim, wind_x, gp_p, gp_lat, gp_lon);
729
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;
736 }
737 }
738
739 // Magnetism
740 // ----------------------------------------------------------------------------
741 else if (jacobian_quantities[q].Target().isMagnetic()) {
742 // Determine grid positions for interpolation from retrieval grids back
743 // to atmospheric grids
744 ArrayOfGridPos gp_p, gp_lat, gp_lon;
745 Index n_p, n_lat, n_lon;
746 get_gp_rq_to_atmgrids(gp_p,
747 gp_lat,
748 gp_lon,
749 n_p,
750 n_lat,
751 n_lon,
752 jacobian_quantities[q].Grids(),
753 atmosphere_dim,
754 p_grid,
755 lat_grid,
756 lon_grid);
757
758 // TODO Could be done without copying.
759 Tensor3 mag_x(n_p, n_lat, n_lon);
760 reshape(mag_x, x_t[ind]);
761
762 Tensor3View target_field(mag_u_field);
763
764 Tensor3 mag_field(
765 target_field.npages(), target_field.nrows(), target_field.ncols());
766 regrid_atmfield_by_gp_oem(
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;
783 }
784 }
785 }
786 } else {
787 ARTS_USER_ERROR ( "Unsupported magnetism type");
788 }
789 }
790
791 // Surface
792 // ----------------------------------------------------------------------------
793 else if (jacobian_quantities[q] == Jacobian::Special::SurfaceString) {
794 surface_props_check(atmosphere_dim,
795 lat_grid,
796 lon_grid,
797 surface_props_data,
798 surface_props_names);
799 ARTS_USER_ERROR_IF (surface_props_data.empty(),
800 "One jacobian quantity belongs to the "
801 "surface category, but *surface_props_data* is empty.");
802
803 const Index isu =
804 find_first(surface_props_names, jacobian_quantities[q].Subtag());
805 ARTS_USER_ERROR_IF (isu < 0,
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*.")
810
811 // Determine grid positions for interpolation from retrieval grids back
812 // to atmospheric grids
813 ArrayOfGridPos gp_lat, gp_lon;
814 Index n_lat, n_lon;
815 get_gp_rq_to_atmgrids(gp_lat,
816 gp_lon,
817 n_lat,
818 n_lon,
819 jacobian_quantities[q].Grids(),
820 atmosphere_dim,
821 lat_grid,
822 lon_grid);
823 // Map values in x back to surface_props_data
824 Matrix surf_x(n_lat, n_lon);
825 reshape(surf_x, x_t[ind]);
826 Matrix surf;
827 regrid_atmsurf_by_gp_oem(surf, atmosphere_dim, surf_x, gp_lat, gp_lon);
828 surface_props_data(isu, joker, joker) = surf;
829 }
830 }
831}
832
833/* Workspace method: Doxygen documentation will be auto-generated */
834void x2artsSensor(Workspace& ws,
835 Matrix& sensor_los,
836 Vector& f_backend,
837 Vector& y_baseline,
838 Sparse& sensor_response,
839 Vector& sensor_response_f,
840 ArrayOfIndex& sensor_response_pol,
841 Matrix& sensor_response_dlos,
842 Vector& sensor_response_f_grid,
843 ArrayOfIndex& sensor_response_pol_grid,
844 Matrix& sensor_response_dlos_grid,
845 Matrix& mblock_dlos,
846 const ArrayOfRetrievalQuantity& jacobian_quantities,
847 const Vector& x,
848 const Agenda& sensor_response_agenda,
849 const Index& sensor_checked,
850 const ArrayOfTime& sensor_time,
851 const Verbosity&) {
852 // Basics
853 //
854 ARTS_USER_ERROR_IF (sensor_checked != 1,
855 "The sensor response must be flagged to have "
856 "passed a consistency check (sensor_checked=1).");
857
858 // Revert transformation
859 Vector x_t(x);
860 transform_x_back(x_t, jacobian_quantities);
861
862 // Main sizes
863 const Index nq = jacobian_quantities.nelem();
864
865 // Jacobian indices
866 ArrayOfArrayOfIndex ji;
867 {
868 bool any_affine;
869 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
870 }
871
872 // Check input
873 ARTS_USER_ERROR_IF (x_t.nelem() != ji[nq - 1][1] + 1,
874 "Length of *x* does not match length implied by "
875 "*jacobian_quantities*.");
876
877 // Flag indicating that y_baseline is not set
878 bool yb_set = false;
879
880 // Shall sensor responses be calculed?
881 bool do_sensor = false;
882
883 // Loop retrieval quantities
884 for (Index q = 0; q < nq; q++) {
885 // Index range of this retrieval quantity
886 const Index np = ji[q][1] - ji[q][0] + 1;
887
888 // Pointing off-set
889 // ----------------------------------------------------------------------------
890 if (jacobian_quantities[q].Target().isPointing()) {
891 // Handle pointing "jitter" seperately
892 if (jacobian_quantities[q].Grids()[0][0] == -1) {
893 ARTS_USER_ERROR_IF (sensor_los.nrows() != np,
894 "Mismatch between pointing jacobian and *sensor_los* found.");
895 // Simply add retrieved off-set(s) to za column of *sensor_los*
896 for (Index i = 0; i < np; i++) {
897 sensor_los(i, 0) += x_t[ji[q][0] + i];
898 }
899 }
900 // Polynomial representation
901 else {
902 ARTS_USER_ERROR_IF (sensor_los.nrows() != sensor_time.nelem(),
903 "Sizes of *sensor_los* and *sensor_time* do not match.");
904 Vector w;
905 for (Index c = 0; c < np; c++) {
906 polynomial_basis_func(w, time_vector(sensor_time), c);
907 for (Index i = 0; i < w.nelem(); i++) {
908 sensor_los(i, 0) += w[i] * x_t[ji[q][0] + c];
909 }
910 }
911 }
912 }
913
914 // Frequency shift or stretch
915 // ----------------------------------------------------------------------------
916 else if (jacobian_quantities[q].Target().isFrequency()) {
917 if (jacobian_quantities[q] == Jacobian::Sensor::FrequencyShift) {
918 ARTS_ASSERT(np == 1);
919 if (x_t[ji[q][0]] != 0) {
920 do_sensor = true;
921 f_backend += x_t[ji[q][0]];
922 }
923 } else if (jacobian_quantities[q] == Jacobian::Sensor::FrequencyStretch) {
924 ARTS_ASSERT(np == 1);
925 if (x_t[ji[q][0]] != 0) {
926 do_sensor = true;
927 Vector w;
928 polynomial_basis_func(w, f_backend, 1);
929 for (Index i = 0; i < w.nelem(); i++) {
930 f_backend[i] += w[i] * x_t[ji[q][0]];
931 }
932 }
933 } else {
934 ARTS_ASSERT(0);
935 }
936 }
937
938 // Baseline fit: polynomial or sinusoidal
939 // ----------------------------------------------------------------------------
940 else if (jacobian_quantities[q] == Jacobian::Sensor::Polyfit ||
941 jacobian_quantities[q] == Jacobian::Sensor::Sinefit) {
942 if (!yb_set) {
943 yb_set = true;
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);
948 y_baseline = 0;
949 }
950
951 for (Index mb = 0; mb < sensor_los.nrows(); ++mb) {
952 calcBaselineFit(y_baseline,
953 x_t,
954 mb,
955 sensor_response,
956 sensor_response_pol_grid,
957 sensor_response_f_grid,
958 sensor_response_dlos_grid,
959 jacobian_quantities[q],
960 q,
961 ji);
962 }
963 }
964 }
965
966 // *y_baseline* not yet set?
967 if (!yb_set) {
968 y_baseline.resize(1);
969 y_baseline[0] = 0;
970 }
971
972 // Recalculate sensor_response?
973 if (do_sensor) {
974 sensor_response_agendaExecute(ws,
975 sensor_response,
976 sensor_response_f,
977 sensor_response_f_grid,
978 sensor_response_pol,
979 sensor_response_pol_grid,
980 sensor_response_dlos,
981 sensor_response_dlos_grid,
982 mblock_dlos,
983 f_backend,
984 sensor_response_agenda);
985 }
986}
987
988/* Workspace method: Doxygen documentation will be auto-generated */
989void x2artsSpectroscopy(const Verbosity&) {
990 ARTS_USER_ERROR ( "Retrievals of spectroscopic variables not yet handled.");
991}
992
993
994#ifdef OEM_SUPPORT
995/* Workspace method: Doxygen documentation will be auto-generated */
996void OEM(Workspace& ws,
997 Vector& x,
998 Vector& yf,
999 Matrix& jacobian,
1000 Matrix& dxdy,
1001 Vector& oem_diagnostics,
1002 Vector& lm_ga_history,
1003 ArrayOfString& errors,
1004 const Vector& xa,
1005 const CovarianceMatrix& covmat_sx,
1006 const Vector& y,
1007 const CovarianceMatrix& covmat_se,
1008 const ArrayOfRetrievalQuantity& jacobian_quantities,
1009 const Agenda& inversion_iterate_agenda,
1010 const String& method,
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,
1018 const Verbosity&) {
1019 // Main sizes
1020 const Index n = covmat_sx.nrows();
1021 const Index m = y.nelem();
1022
1023 // Checks
1024 covmat_sx.compute_inverse();
1025 covmat_se.compute_inverse();
1026
1027 OEM_checks(ws,
1028 x,
1029 yf,
1030 jacobian,
1031 inversion_iterate_agenda,
1032 xa,
1033 covmat_sx,
1034 y,
1035 covmat_se,
1036 jacobian_quantities,
1037 method,
1038 x_norm,
1039 max_iter,
1040 stop_dx,
1041 lm_ga_settings,
1042 clear_matrices,
1043 display_progress);
1044
1045 // Size diagnostic output and init with NaNs
1046 oem_diagnostics.resize(5);
1047 oem_diagnostics = NAN;
1048 //
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;
1053 } else {
1054 lm_ga_history.resize(0);
1055 }
1056
1057 // Check for start vector and precomputed yf, jacobian
1058 if (x.nelem() != n) {
1059 x = xa;
1060 yf.resize(0);
1061 jacobian.resize(0, 0);
1062 }
1063
1064 // If no precomputed value given, we compute yf and jacobian to
1065 // compute initial cost (and use in the first OEM iteration).
1066 if (yf.nelem() == 0) {
1067 inversion_iterate_agendaExecute(
1068 ws, yf, jacobian, xa, 1, 0, inversion_iterate_agenda);
1069 }
1070
1071 ARTS_USER_ERROR_IF (yf.nelem() not_eq y.nelem(),
1072 "Mismatch between simulated y and input y.\n"
1073 "Input y is size ", y.nelem(), " but simulated y is ",
1074 yf.nelem(), "\n"
1075 "Use your frequency grid vector and your sensor response matrix to match simulations with measurements.\n")
1076
1077 // TODO: Get this from invlib log.
1078 // Start value of cost function
1079 Numeric cost_start = NAN;
1080 if (method == "ml" || method == "lm" || display_progress ||
1081 max_start_cost > 0) {
1082 Vector dy = y;
1083 dy -= yf;
1084 Vector sdy = y;
1085 mult_inv(ExhaustiveMatrixView{sdy}, covmat_se, ExhaustiveMatrixView{dy});
1086 Vector dx = x;
1087 dx -= xa;
1088 Vector sdx = x;
1089 mult_inv(ExhaustiveMatrixView{sdx}, covmat_sx, ExhaustiveMatrixView{dx});
1090 cost_start = dx * sdx + dy * sdy;
1091 cost_start /= static_cast<Numeric>(m);
1092 }
1093 oem_diagnostics[1] = cost_start;
1094
1095 // Handle cases with too large start cost
1096 if (max_start_cost > 0 && cost_start > max_start_cost) {
1097 // Flag no inversion in oem_diagnostics, and let x to be undefined
1098 oem_diagnostics[0] = 99;
1099 //
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
1104 << endl;
1105 }
1106 }
1107 // Otherwise do inversion
1108 else {
1109 bool apply_norm = false;
1110 oem::Matrix T{};
1111 if (x_norm.nelem() == n) {
1112 T.resize(n, n);
1113 T *= 0.0;
1114 T.diagonal() = x_norm;
1115 for (Index i = 0; i < n; i++) {
1116 T(i, i) = x_norm[i];
1117 }
1118 apply_norm = true;
1119 }
1120
1121 oem::CovarianceMatrix Se(covmat_se), Sa(covmat_sx);
1122 oem::Vector xa_oem(xa), y_oem(y), x_oem(x);
1123 oem::AgendaWrapper aw(&ws,
1124 (unsigned int)m,
1125 (unsigned int)n,
1126 jacobian,
1127 yf,
1128 &inversion_iterate_agenda);
1129 oem::OEM_STANDARD<oem::AgendaWrapper> oem(aw, xa_oem, Sa, Se);
1130 oem::OEM_MFORM<oem::AgendaWrapper> oem_m(aw, xa_oem, Sa, Se);
1131 int oem_verbosity = static_cast<int>(display_progress);
1132
1133 int return_code = 0;
1134
1135 try {
1136 if (method == "li") {
1137 oem::Std s(T, apply_norm);
1138 oem::GN gn(stop_dx, 1, s); // Linear case, only one step.
1139 return_code = oem.compute<oem::GN, oem::ArtsLog>(
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") {
1143 ARTS_USER_ERROR(method, " is not supported")
1144 oem::Std s(T, apply_norm);
1145 oem::GN gn(stop_dx, 1, s); // Linear case, only one step.
1146// return_code = oem_m.compute<oem::GN, oem::ArtsLog>(
1147// x_oem, y_oem, gn, oem_verbosity, lm_ga_history, true);
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);
1151 oem::GN_CG gn(stop_dx, 1, cg); // Linear case, only one step.
1152 return_code = oem.compute<oem::GN_CG, oem::ArtsLog>(
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);
1157 oem::GN_CG gn(stop_dx, 1, cg); // Linear case, only one step.
1158 return_code = oem_m.compute<oem::GN_CG, oem::ArtsLog>(
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") {
1162 oem::Std s(T, apply_norm);
1163 oem::GN gn(stop_dx, (unsigned int)max_iter, s);
1164 return_code = oem.compute<oem::GN, oem::ArtsLog>(
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") {
1168 ARTS_USER_ERROR(method, " is not supported")
1169 oem::Std s(T, apply_norm);
1170 oem::GN gn(stop_dx, (unsigned int)max_iter, s);
1171// return_code = oem_m.compute<oem::GN, oem::ArtsLog>(
1172// x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
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);
1177 return_code = oem.compute<oem::GN_CG, oem::ArtsLog>(
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);
1183 return_code = oem_m.compute<oem::GN_CG, oem::ArtsLog>(
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")) {
1187 oem::Std s(T, apply_norm);
1188 Sparse diagonal = Sparse::diagonal(covmat_sx.inverse_diagonal());
1189 CovarianceMatrix SaDiag{};
1190 SaDiag.add_correlation_inverse(Block(Range(0, n),
1191 Range(0, n),
1192 std::make_pair(0, 0),
1193 make_shared<Sparse>(diagonal)));
1194 oem::CovarianceMatrix SaInvLM = inv(oem::CovarianceMatrix(SaDiag));
1195 oem::LM lm(SaInvLM, s);
1196
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]);
1205
1206 return_code = oem.compute<oem::LM&, oem::ArtsLog>(
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;
1211 }
1212 } else if ((method == "lm_cg") || (method == "ml_cg")) {
1213 oem::CG cg(T, apply_norm, 1e-10, 0);
1214
1215 Sparse diagonal = Sparse::diagonal(covmat_sx.inverse_diagonal());
1216 CovarianceMatrix SaDiag{};
1217 SaDiag.add_correlation_inverse(Block(Range(0, n),
1218 Range(0, n),
1219 std::make_pair(0, 0),
1220 make_shared<Sparse>(diagonal)));
1221 oem::LM_CG lm(SaDiag, cg);
1222
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]);
1229
1230 return_code = oem.compute<oem::LM_CG&, oem::ArtsLog>(
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;
1235 }
1236 }
1237
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);
1246 x_oem *= NAN;
1247 std::vector<std::string> sv = oem::handle_nested_exception(e);
1248 for (auto& s : sv) {
1249 std::stringstream ss{s};
1250 std::string t{};
1251 while (std::getline(ss, t)) {
1252 errors.push_back(t.c_str());
1253 }
1254 }
1255 } catch (...) {
1256 throw;
1257 }
1258
1259 x = x_oem;
1260 yf = aw.get_measurement_vector();
1261
1262 // Shall empty jacobian and dxdy be returned?
1263 if (clear_matrices) {
1264 jacobian.resize(0, 0);
1265 dxdy.resize(0, 0);
1266 } else if (oem_diagnostics[0] <= 2) {
1267 dxdy.resize(n, m);
1268 Matrix tmp1(n, m), tmp2(n, n), tmp3(n, n);
1269 mult_inv(tmp1, transpose(jacobian), covmat_se);
1270 mult(tmp2, tmp1, jacobian);
1271 add_inv(tmp2, covmat_sx);
1272 inv(tmp3, tmp2);
1273 mult(dxdy, tmp3, tmp1);
1274 }
1275 }
1276}
1277
1278/* Workspace method: Doxygen documentation will be auto-generated */
1279void covmat_soCalc(Matrix& covmat_so,
1280 const Matrix& dxdy,
1281 const CovarianceMatrix& covmat_se,
1282 const Verbosity& /*v*/) {
1283 Index n(dxdy.nrows()), m(dxdy.ncols());
1284 Matrix tmp1(m, n);
1285
1286 ARTS_USER_ERROR_IF ((m == 0) || (n == 0),
1287 "The gain matrix *dxdy* is required to compute the observation error covariance matrix.");
1288 ARTS_USER_ERROR_IF ((covmat_se.nrows() != m) || (covmat_se.ncols() != m),
1289 "The covariance matrix covmat_se has invalid dimensions.");
1290
1291 covmat_so.resize(n, n);
1292 mult(tmp1, covmat_se, transpose(dxdy));
1293 mult(covmat_so, dxdy, tmp1);
1294}
1295
1296/* Workspace method: Doxygen documentation will be auto-generated */
1297void covmat_ssCalc(Matrix& covmat_ss,
1298 const Matrix& avk,
1299 const CovarianceMatrix& covmat_sx,
1300 const Verbosity& /*v*/) {
1301 Index n(avk.ncols());
1302 Matrix tmp1(n, n), tmp2(n, n);
1303
1304 ARTS_USER_ERROR_IF (n == 0,
1305 "The averaging kernel matrix *dxdy* is required to compute the smoothing error covariance matrix.");
1306 ARTS_USER_ERROR_IF ((covmat_sx.nrows() != n) || (covmat_sx.ncols() != n),
1307 "The covariance matrix *covmat_sx* invalid dimensions.");
1308
1309 covmat_ss.resize(n, n);
1310
1311 // Sign doesn't matter since we're dealing with a quadratic form.
1312 id_mat(tmp1);
1313 tmp1 -= avk;
1314
1315 mult(tmp2, covmat_sx, transpose(tmp1));
1316 mult(covmat_ss, tmp1, tmp2);
1317}
1318
1319/* Workspace method: Doxygen documentation will be auto-generated */
1320void MatrixFromCovarianceMatrix(Matrix& S,
1321 const CovarianceMatrix& Sc,
1322 const Verbosity& /*v*/) {
1323 S = Matrix(Sc);
1324}
1325
1326/* Workspace method: Doxygen documentation will be auto-generated */
1327void avkCalc(Matrix& avk,
1328 const Matrix& dxdy,
1329 const Matrix& jacobian,
1330 const Verbosity& /*v*/) {
1331 Index m(jacobian.nrows()), n(jacobian.ncols());
1332 ARTS_USER_ERROR_IF ((m == 0) || (n == 0),
1333 "The Jacobian matrix is empty.");
1334 ARTS_USER_ERROR_IF ((dxdy.nrows() != n) || (dxdy.ncols() != m),
1335 "Matrices have inconsistent sizes.\n"
1336 " Size of gain matrix: ", dxdy.nrows(), " x ", dxdy.ncols(),
1337 "\n"
1338 " Size of Jacobian: ", jacobian.nrows(), " x ",
1339 jacobian.ncols(), "\n")
1340
1341 avk.resize(n, n);
1342 mult(avk, dxdy, jacobian);
1343}
1344
1345#else
1346
1347void covmat_soCalc(Matrix& /* covmat_so */,
1348 const Matrix& /* dxdy */,
1349 const CovarianceMatrix& /* covmat_se*/,
1350 const Verbosity& /*v*/) {
1351 ARTS_USER_ERROR (
1352 "WSM is not available because ARTS was compiled without "
1353 "OEM support.");
1354}
1355
1356void covmat_ssCalc(Matrix& /*covmat_ss*/,
1357 const Matrix& /*avk*/,
1358 const CovarianceMatrix& /*covmat_sx*/,
1359 const Verbosity& /*v*/) {
1360 ARTS_USER_ERROR (
1361 "WSM is not available because ARTS was compiled without "
1362 "OEM support.");
1363}
1364
1365void avkCalc(Matrix& /* avk */,
1366 const Matrix& /* dxdy */,
1367 const Matrix& /* jacobian */,
1368 const Verbosity& /*v*/) {
1369 ARTS_USER_ERROR (
1370 "WSM is not available because ARTS was compiled without "
1371 "OEM support.");
1372}
1373
1374void OEM(Workspace&,
1375 Vector&,
1376 Vector&,
1377 Matrix&,
1378 Matrix&,
1379 Vector&,
1380 Vector&,
1381 ArrayOfString&,
1382 const Vector&,
1383 const CovarianceMatrix&,
1384 const Vector&,
1385 const CovarianceMatrix&,
1386 const Index&,
1387 const ArrayOfRetrievalQuantity&,
1388 const ArrayOfArrayOfIndex&,
1389 const Agenda&,
1390 const String&,
1391 const Numeric&,
1392 const Vector&,
1393 const Index&,
1394 const Numeric&,
1395 const Vector&,
1396 const Index&,
1397 const Index&,
1398 const Verbosity&) {
1399 ARTS_USER_ERROR (
1400 "WSM is not available because ARTS was compiled without "
1401 "OEM support.");
1402}
1403#endif
1404
1405#if defined(OEM_SUPPORT) && 0
1406
1407#include "agenda_wrapper_mpi.h"
1408#include "oem_mpi.h"
1409
1410//
1411// Performs manipulations of workspace variables necessary for distributed
1412// retrievals with MPI:
1413//
1414// - Splits up sensor positions evenly over processes
1415// - Splits up inverse covariance matrices.
1416//
1417void MPI_Initialize(Matrix& sensor_los,
1418 Matrix& sensor_pos,
1419 Vector& sensor_time) {
1420 int initialized;
1421
1422 MPI_Initialized(&initialized);
1423 if (!initialized) {
1424 MPI_Init(nullptr, nullptr);
1425 }
1426
1427 int rank, nprocs;
1428
1429 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1430 MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
1431
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);
1436
1437 //
1438 // Split up sensor positions.
1439 //
1440
1441 if (rank < remainder) {
1442 mblock_range += 1;
1443 mblock_start += rank;
1444 } else {
1445 mblock_start += remainder;
1446 }
1447
1448 if (nmblock > 0) {
1449 Range range = Range(mblock_start, mblock_range);
1450
1451 Matrix tmp_m = sensor_los(range, joker);
1452 sensor_los = tmp_m;
1453
1454 tmp_m = sensor_pos(range, joker);
1455 sensor_pos = tmp_m;
1456
1457 Vector tmp_v = sensor_time[range];
1458 sensor_time = tmp_v;
1459 } else {
1460 sensor_los.resize(0, 0);
1461 sensor_pos.resize(0, 0);
1462 sensor_time.resize(0);
1463 }
1464}
1465
1466void OEM_MPI(Workspace& ws,
1467 Vector& x,
1468 Vector& yf,
1469 Matrix& jacobian,
1470 Matrix& dxdy,
1471 Vector& oem_diagnostics,
1472 Vector& lm_ga_history,
1473 Matrix& sensor_los,
1474 Matrix& sensor_pos,
1475 Vector& sensor_time,
1476 CovarianceMatrix& covmat_sx,
1477 CovarianceMatrix& covmat_se,
1478 const Vector& xa,
1479 const Vector& y,
1480 const ArrayOfRetrievalQuantity& jacobian_quantities,
1481 const Agenda& inversion_iterate_agenda,
1482 const String& method,
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,
1490 const Verbosity& /*v*/) {
1491 // Main sizes
1492 const Index n = covmat_sx.nrows();
1493 const Index m = y.nelem();
1494
1495 // Check WSVs
1496 OEM_checks(ws,
1497 x,
1498 yf,
1499 jacobian,
1500 inversion_iterate_agenda,
1501 xa,
1502 covmat_sx,
1503 covmat_se,
1504 jacobian_quantities,
1505 method,
1506 x_norm,
1507 max_iter,
1508 stop_dx,
1509 lm_ga_settings,
1510 clear_matrices,
1511 display_progress);
1512
1513 // Calculate spectrum and Jacobian for a priori state
1514 // Jacobian is also input to the agenda, and to flag this is this first
1515 // call, this WSV must be set to be empty.
1516 jacobian.resize(0, 0);
1517
1518 // Initialize MPI environment.
1519 MPI_Initialize(sensor_los, sensor_pos, sensor_time);
1520
1521 // Setup distributed matrices.
1522 MPICovarianceMatrix SeInvMPI(covmat_se);
1523 MPICovarianceMatrix SaInvMPI(covmat_sx);
1524
1525 // Create temporary MPI vector from local results and use conversion to
1526 // standard vector to broadcast results to all processes.
1527 oem::Vector tmp;
1528 inversion_iterate_agendaExecute(
1529 ws, tmp, jacobian, xa, 1, inversion_iterate_agenda);
1530 yf = MPIVector(tmp);
1531
1532 // Size diagnostic output and init with NaNs
1533 oem_diagnostics.resize(5);
1534 oem_diagnostics = NAN;
1535 //
1536 if (method == "ml" || method == "lm") {
1537 lm_ga_history.resize(max_iter);
1538 lm_ga_history = NAN;
1539 } else {
1540 lm_ga_history.resize(0);
1541 }
1542
1543 // Start value of cost function. Covariance matrices are already distributed
1544 // over processes, so we need to use invlib matrix algebra.
1545 Numeric cost_start = NAN;
1546 if (method == "ml" || method == "lm" || display_progress ||
1547 max_start_cost > 0) {
1548 oem::Vector dy = y;
1549 dy -= yf;
1550 cost_start = dot(dy, SeInvMPI * dy);
1551 }
1552 oem_diagnostics[1] = cost_start;
1553
1554 // Handle cases with too large start cost
1555 if (max_start_cost > 0 && cost_start > max_start_cost) {
1556 // Flag no inversion in oem_diagnostics, and let x to be undefined
1557 oem_diagnostics[0] = 99;
1558 //
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
1563 << endl;
1564 }
1565 }
1566
1567 // Otherwise do inversion
1568 else {
1569 // Size remaining output arguments
1570 x.resize(n);
1571 dxdy.resize(n, m);
1572
1573 OEMVector xa_oem(xa), y_oem(y), x_oem;
1574 oem::AgendaWrapperMPI aw(&ws, &inversion_iterate_agenda, m, n);
1575
1576 OEM_PS_PS_MPI<AgendaWrapperMPI> oem(aw, xa_oem, SaInvMPI, SeInvMPI);
1577
1578 // Call selected method
1579 int return_value = 99;
1580
1581 if (method == "li") {
1582 oem::CG cg(1e-12, 0);
1583 oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1584 return_value = oem.compute<oem::GN_CG, invlib::MPILog>(
1585 x_oem, y_oem, gn, 2 * (int)display_progress);
1586 } else if (method == "gn") {
1587 oem::CG cg(1e-12, 0);
1588 oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1589 return_value = oem.compute<oem::GN_CG, invlib::MPILog>(
1590 x_oem, y_oem, gn, 2 * (int)display_progress);
1591 } else if ((method == "lm") || (method == "ml")) {
1592 oem::CG cg(1e-12, 0);
1593 LM_CG_MPI lm(SaInvMPI, cg);
1594
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]);
1602
1603 return_value = oem.compute<oem::LM_CG_MPI, invlib::MPILog>(
1604 x_oem, y_oem, lm, 2 * (int)display_progress);
1605 }
1606
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);
1611
1612 x = x_oem;
1613 // Shall empty jacobian and dxdy be returned?
1614 if (clear_matrices && (oem_diagnostics[0])) {
1615 jacobian.resize(0, 0);
1616 dxdy.resize(0, 0);
1617 }
1618 }
1619 MPI_Finalize();
1620}
1621
1622#else
1623
1624void OEM_MPI(Workspace&,
1625 Vector&,
1626 Vector&,
1627 Matrix&,
1628 Matrix&,
1629 Vector&,
1630 Vector&,
1631 Matrix&,
1632 Matrix&,
1633 Vector&,
1634 CovarianceMatrix&,
1635 CovarianceMatrix&,
1636 const Vector&,
1637 const Vector&,
1638 const Index&,
1639 const ArrayOfRetrievalQuantity&,
1640 const Agenda&,
1641 const String&,
1642 const Numeric&,
1643 const Vector&,
1644 const Index&,
1645 const Numeric&,
1646 const Vector&,
1647 const Index&,
1648 const Index&,
1649 const Verbosity&) {
1650 ARTS_USER_ERROR (
1651 "You have to compile ARTS with OEM support "
1652 " and enable MPI to use OEM_MPI.");
1653}
1654
1655#endif // OEM_SUPPORT && ENABLE_MPI
array.h
This file contains the definition of Array.
find_first
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition array.h:173
arts.h
The global header file for ARTS.
arts_omp.h
Header file for helper functions for OpenMP.
time_vector
Vector time_vector(const ArrayOfTime &times)
Converts from each Time to seconds and returns as Vector.
Definition artstime.cc:149
inversion_iterate_agendaExecute
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:26882
water_p_eq_agendaExecute
void water_p_eq_agendaExecute(Workspace &ws, Tensor3 &water_p_eq_field, const Tensor3 &t_field, const Agenda &input_agenda)
Definition auto_md.cc:27914
sensor_response_agendaExecute
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)
Definition auto_md.cc:27673
auto_md.h
MatrixFromCovarianceMatrix
void MatrixFromCovarianceMatrix(Matrix &output, const CovarianceMatrix &input, const Verbosity &verbosity)
WORKSPACE METHOD: MatrixFromCovarianceMatrix.
check_input.h
chk_contains
Index chk_contains(const String &x_name, const Array< T > &x, const T &what)
Check if an array contains a value.
Definition check_input.h:131
Agenda
The Agenda class.
Definition agenda_class.h:52
Array::nelem
Index nelem() const ARTS_NOEXCEPT
Definition array.h:75
CovarianceMatrix::compute_inverse
void compute_inverse() const
Compute the inverse of this correlation matrix.
Definition covariance_matrix.cc:388
CovarianceMatrix::add_correlation_inverse
void add_correlation_inverse(Block c)
Add block inverse to covariance matrix.
Definition covariance_matrix.cc:509
CovarianceMatrix::inverse_diagonal
Vector inverse_diagonal() const
Diagonal of the inverse of the covariance matrix as vector.
Definition covariance_matrix.cc:526
Workspace
Workspace class.
Definition workspace_ng.h:36
my_basic_string< char >
oem::AgendaWrapperMPI
Interface for distributed ARTS forward model.
Definition oem_mpi.h:56
oem::AgendaWrapper
Interface to ARTS inversion_iterate_agenda.
Definition oem.h:461
oem::ArtsLog
OEM log output.
Definition oem.h:243
oem::NormalizingSolver
Normalizing solver.
Definition oem.h:112
mult
void mult(MatrixView C, ConstMatrixView A, const Block &B)
Definition covariance_matrix.cc:22
add_inv
void add_inv(MatrixView A, const CovarianceMatrix &B)
Definition covariance_matrix.cc:608
mult_inv
void mult_inv(MatrixView C, ConstMatrixView A, const CovarianceMatrix &B)
Definition covariance_matrix.cc:571
debug.h
Helper macros for debugging.
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition debug.h:86
ARTS_USER_ERROR
#define ARTS_USER_ERROR(...)
Definition debug.h:153
ARTS_USER_ERROR_IF
#define ARTS_USER_ERROR_IF(condition,...)
Definition debug.h:137
polynomial_basis_func
void polynomial_basis_func(Vector &b, const Vector &x, const Index &poly_coeff)
Calculates polynomial basis functions.
Definition jacobian.cc:831
transform_x_back
void transform_x_back(Vector &x_t, const ArrayOfRetrievalQuantity &jqs, bool revert_functional_transforms)
Handles back-transformations of the state vector.
Definition jacobian.cc:216
jac_ranges_indices
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.
Definition jacobian.cc:29
transform_x
void transform_x(Vector &x, const ArrayOfRetrievalQuantity &jqs)
Handles transformations of the state vector.
Definition jacobian.cc:139
calcBaselineFit
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.
Definition jacobian.cc:855
jacobian.h
Routines for setting up the jacobian.
covmat_ssCalc
void covmat_ssCalc(Matrix &, const Matrix &, const CovarianceMatrix &, const Verbosity &)
WORKSPACE METHOD: covmat_ssCalc.
Definition m_oem.cc:1356
particle_bulkprop_fieldClip
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.
Definition m_oem.cc:45
x2artsSensor
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.
Definition m_oem.cc:834
OEM_MPI
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 &)
Definition m_oem.cc:1624
OEM
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 &)
Definition m_oem.cc:1374
x2artsAtmAndSurf
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.
Definition m_oem.cc:483
x2artsSpectroscopy
void x2artsSpectroscopy(const Verbosity &)
WORKSPACE METHOD: x2artsSpectroscopy.
Definition m_oem.cc:989
avkCalc
void avkCalc(Matrix &, const Matrix &, const Matrix &, const Verbosity &)
WORKSPACE METHOD: avkCalc.
Definition m_oem.cc:1365
xClip
void xClip(Vector &x, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &ijq, const Numeric &limit_low, const Numeric &limit_high, const Verbosity &)
WORKSPACE METHOD: xClip.
Definition m_oem.cc:96
covmat_soCalc
void covmat_soCalc(Matrix &, const Matrix &, const CovarianceMatrix &, const Verbosity &)
WORKSPACE METHOD: covmat_soCalc.
Definition m_oem.cc:1347
vmr_fieldClip
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.
Definition m_oem.cc:71
reshape
void reshape(Tensor3View X, ConstVectorView x)
reshape
Definition math_funcs.cc:763
flat
void flat(VectorView x, ConstMatrixView X)
flat
Definition math_funcs.cc:709
math_funcs.h
Constant::k
constexpr Numeric k
Boltzmann constant convenience name [J/K].
Definition arts_constants.h:152
oem
Definition oem.h:28
oem::OEM_MFORM
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::MFORM > OEM_MFORM
OEM m form.
Definition oem.h:94
oem::LM_CG
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, CG > LM_CG
Levenberg-Marquardt (LM) optimization using normed CG solver.
Definition oem.h:174
oem::OEM_STANDARD
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::STANDARD, invlib::Rodgers531 > OEM_STANDARD
OEM standard form.
Definition oem.h:62
oem::LM
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, Std > LM
Levenberg-Marquardt (LM) optimization using normed ARTS QR solver.
Definition oem.h:172
oem::CovarianceMatrix
invlib::Matrix< ArtsCovarianceMatrixWrapper > CovarianceMatrix
invlib wrapper type for ARTS the ARTS covariance class.
Definition oem.h:37
oem::Matrix
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition oem.h:33
oem::GN_CG
invlib::GaussNewton< Numeric, CG > GN_CG
Gauss-Newton (GN) optimization using normed CG solver.
Definition oem.h:170
oem::handle_nested_exception
std::vector< std::string > handle_nested_exception(const E &e, int level=0)
Handle exception encountered within invlib.
Definition oem.h:424
oem::GN
invlib::GaussNewton< Numeric, Std > GN
OEM Gauss-Newton optimization using normed ARTS QR solver.
Definition oem.h:168
oem::Vector
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition oem.h:31
oem::MPICovarianceMatrix
invlib::Matrix< invlib::MPIMatrix< invlib::Timer< ArtsCovarianceMatrixWrapper > > > MPICovarianceMatrix
MPI-distributed covariance matrix type.
Definition oem_mpi.h:27
oem::MPIVector
invlib::Vector< invlib::MPIVector< invlib::Timer< ArtsVector > > > MPIVector
MPI-distributed vector type.
Definition oem_mpi.h:29
oem.h
Defines the ARTS interface to the invlib library.
OEM_checks
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.
Definition oem.h:662
Tensor4Clip
void Tensor4Clip(Tensor4 &x, const Index &iq, const Numeric &limit_low, const Numeric &limit_high)
Clip Tensor4.
Definition oem.h:581
oem_mpi.h
Optimal estimation method for MPI-distributed retrieval.
physics_funcs.h
This file contains declerations of functions of physical character.
number_density
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition physics_funcs.h:51
rte.h
Declaration of functions in rte.cc.
regrid_atmfield_by_gp
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.
Definition special_interp.cc:293
regrid_atmfield_by_gp_oem
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.
Definition special_interp.cc:490
get_gp_atmsurf_to_rq
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.
Definition special_interp.cc:375
get_gp_rq_to_atmgrids
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.
Definition special_interp.cc:396
regrid_atmsurf_by_gp_oem
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.
Definition special_interp.cc:641
get_gp_atmgrids_to_rq
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.
Definition special_interp.cc:349
special_interp.h
Header file for special_interp.cc.
SpeciesTag
Species::Tag SpeciesTag
Definition species_tags.h:94
surface_props_check
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
Definition surface.cc:122
w
#define w
a
#define a
c
#define c