ARTS 2.5.10 (git: 2f1c442c)
m_oem.cc
Go to the documentation of this file.
1/* Copyright (C) 2015
2 Patrick Eriksson <patrick.eriksson@chalmers.se>
3 Stefan Buehler <sbuehler@ltu.se>
4
5 This program is free software; you can redistribute it and/or modify it
6 under the terms of the oem::GNU General Public License as published by the
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 oem::GNU General Public License for more details.
14
15 You should have received a copy of the oem::GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18 USA. */
19
20/*===========================================================================
21 === File description
22 ===========================================================================*/
23
35/*===========================================================================
36 === External declarations
37 ===========================================================================*/
38
39#include <cmath>
40#include <iterator>
41#include <sstream>
42#include <stdexcept>
43#include <string>
44#include "array.h"
45#include "arts.h"
46#include "arts_omp.h"
47#include "auto_md.h"
48#include "jacobian.h"
49#include "math_funcs.h"
50#include "physics_funcs.h"
51#include "rte.h"
52#include "special_interp.h"
53#include "surface.h"
54#include "check_input.h"
55
56#pragma GCC diagnostic ignored "-Wconversion"
57
58#ifdef OEM_SUPPORT
59#include "oem.h"
60#endif
61
62/* Workspace method: Doxygen documentation will be auto-generated */
63void particle_bulkprop_fieldClip(Tensor4& particle_bulkprop_field,
64 const ArrayOfString& particle_bulkprop_names,
65 const String& bulkprop_name,
66 const Numeric& limit_low,
67 const Numeric& limit_high,
68 const Verbosity&) {
69 Index iq = -1;
70 if (bulkprop_name == "ALL") {
71 }
72
73 else {
74 for (Index i = 0; i < particle_bulkprop_names.nelem(); i++) {
75 if (particle_bulkprop_names[i] == bulkprop_name) {
76 iq = i;
77 break;
78 }
79 }
80 ARTS_USER_ERROR_IF (iq < 0,
81 "Could not find ", bulkprop_name,
82 " in particle_bulkprop_names.\n")
83 }
84
85 Tensor4Clip(particle_bulkprop_field, iq, limit_low, limit_high);
86}
87
88/* Workspace method: Doxygen documentation will be auto-generated */
89void vmr_fieldClip(Tensor4& vmr_field,
90 const ArrayOfArrayOfSpeciesTag& abs_species,
91 const String& species,
92 const Numeric& limit_low,
93 const Numeric& limit_high,
94 const Verbosity&) {
95 Index iq = -1;
96 if (species == "ALL") {
97 }
98
99 else {
100 for (Index i = 0; i < abs_species.nelem(); i++) {
101 if (abs_species[i].Species() == SpeciesTag(species).Spec()) {
102 iq = i;
103 break;
104 }
105 }
106 ARTS_USER_ERROR_IF (iq < 0,
107 "Could not find ", species, " in abs_species.\n")
108 }
109
110 Tensor4Clip(vmr_field, iq, limit_low, limit_high);
111}
112
113/* Workspace method: Doxygen documentation will be auto-generated */
114void xClip(Vector& x,
115 const ArrayOfRetrievalQuantity& jacobian_quantities,
116 const Index& ijq,
117 const Numeric& limit_low,
118 const Numeric& limit_high,
119 const Verbosity&) {
120 // Sizes
121 const Index nq = jacobian_quantities.nelem();
122
123 ARTS_USER_ERROR_IF (ijq < -1, "Argument *ijq* must be >= -1.");
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")
129
130 // Jacobian indices
131 ArrayOfArrayOfIndex ji;
132 {
133 bool any_affine;
134 jac_ranges_indices(ji, any_affine, jacobian_quantities);
135 }
136
137 Index ifirst = 0, ilast = x.nelem() - 1;
138 if (ijq > -1) {
139 ifirst = ji[ijq][0];
140 ilast = ji[ijq][1];
141 }
142
143 if (!std::isinf(limit_low)) {
144 for (Index i = ifirst; i <= ilast; i++) {
145 if (x[i] < limit_low) x[i] = limit_low;
146 }
147 }
148 if (!std::isinf(limit_high)) {
149 for (Index i = ifirst; i <= ilast; i++) {
150 if (x[i] > limit_high) x[i] = limit_high;
151 }
152 }
153}
154
155/* Workspace method: Doxygen documentation will be auto-generated */
156void xaStandard(Workspace& ws,
157 Vector& xa,
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,
181 const Verbosity&) {
182 // Basics
183 //
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).");
193
194 // Jacobian indices
195 ArrayOfArrayOfIndex ji;
196 {
197 bool any_affine;
198 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
199 }
200
201 // Sizes
202 const Index nq = jacobian_quantities.nelem();
203 //
204 xa.resize(ji[nq - 1][1] + 1);
205
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);
211
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,
217 gp_lat,
218 gp_lon,
219 jacobian_quantities[q],
220 atmosphere_dim,
221 p_grid,
222 lat_grid,
223 lon_grid);
224 Tensor3 t_x;
225 regrid_atmfield_by_gp(t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
226 flat(xa[ind], t_x);
227 }
228
229 // Abs species
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);
234
235 if (jacobian_quantities[q].Mode() == "rel") {
236 // This one is simple, just a vector of ones
237 xa[ind] = 1;
238 } else {
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,
242 gp_lat,
243 gp_lon,
244 jacobian_quantities[q],
245 atmosphere_dim,
246 p_grid,
247 lat_grid,
248 lon_grid);
249 Tensor3 vmr_x;
250 regrid_atmfield_by_gp(vmr_x,
251 atmosphere_dim,
252 vmr_field(isp, joker, joker, joker),
253 gp_p,
254 gp_lat,
255 gp_lon);
256
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*
261 Tensor3 t_x;
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)
265 Index i = 0;
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++) {
269 xa[ji[q][0] + i] =
270 vmr_x(i1, i2, i3) *
271 number_density(jacobian_quantities[q].Grids()[0][i1],
272 t_x(i1, i2, i3));
273 i += 1;
274 }
275 }
276 }
277 } else if (jacobian_quantities[q].Mode() == "rh") {
278 // Here we need to also interpolate *t_field*
279 Tensor3 t_x;
280 regrid_atmfield_by_gp(
281 t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
282 Tensor3 water_p_eq;
283 water_p_eq_agendaExecute(ws, water_p_eq, t_x, water_p_eq_agenda);
284 // Calculate relative humidity (vmr*p/p_sat)
285 Index i = 0;
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);
292 i += 1;
293 }
294 }
295 }
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;
300 Index i = 0;
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++) {
304 const Numeric e =
305 vmr_x(i1, i2, i3) * jacobian_quantities[q].Grids()[0][i1];
306 const Numeric r =
307 0.622 * e / (jacobian_quantities[q].Grids()[0][i1] - e);
308 xa[ji[q][0] + i] = r / (1 + r);
309 i += 1;
310 }
311 }
312 }
313 } else {
314 ARTS_ASSERT(0);
315 }
316 }
317 }
318
319 // Scattering species
320 else if (jacobian_quantities[q] == Jacobian::Special::ScatteringString) {
321 if (cloudbox_on) {
322 ARTS_USER_ERROR_IF (particle_bulkprop_field.empty(),
323 "One jacobian quantity belongs to the "
324 "scattering species category, but *particle_bulkprop_field* "
325 "is empty.");
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*.");
330
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*.")
338
339 ArrayOfGridPos gp_p, gp_lat, gp_lon;
341 gp_lat,
342 gp_lon,
343 jacobian_quantities[q],
344 atmosphere_dim,
345 p_grid,
346 lat_grid,
347 lon_grid);
348 Tensor3 pbp_x;
350 atmosphere_dim,
351 particle_bulkprop_field(isp, joker, joker, joker),
352 gp_p,
353 gp_lat,
354 gp_lon);
355 flat(xa[ind], pbp_x);
356 } else {
357 xa[ind] = 0;
358 }
359 }
360
361 // Wind
362 else if (jacobian_quantities[q].Target().isWind()) {
363 ConstTensor3View source_field(wind_u_field);
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;
368 }
369
370 // Determine grid positions for interpolation from retrieval grids back
371 // to atmospheric grids
372 ArrayOfGridPos gp_p, gp_lat, gp_lon;
374 gp_lat,
375 gp_lon,
376 jacobian_quantities[q],
377 atmosphere_dim,
378 p_grid,
379 lat_grid,
380 lon_grid);
381
382 Tensor3 wind_x;
384 wind_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
385 flat(xa[ind], wind_x);
386 }
387
388 // Magnetism
389 else if (jacobian_quantities[q].Target().isMagnetic()) {
390 if (jacobian_quantities[q] == Jacobian::Atm::MagneticMagnitude) {
391 // Determine grid positions for interpolation from retrieval grids back
392 // to atmospheric grids
393 ArrayOfGridPos gp_p, gp_lat, gp_lon;
395 gp_lat,
396 gp_lon,
397 jacobian_quantities[q],
398 atmosphere_dim,
399 p_grid,
400 lat_grid,
401 lon_grid);
402
403 //all three component's hyoptenuse is the strength
404 Tensor3 mag_u, mag_v, mag_w;
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);
411
412 Tensor3 mag_x(gp_p.nelem(), gp_lat.nelem(), gp_lon.nelem());
413 for (Index i = 0; i < gp_p.nelem(); i++)
414 for (Index j = 0; j < gp_lat.nelem(); j++)
415 for (Index k = 0; k < gp_lon.nelem(); k++)
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);
418 } else {
419 ConstTensor3View source_field(mag_u_field);
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;
424 }
425
426 // Determine grid positions for interpolation from retrieval grids back
427 // to atmospheric grids
428 ArrayOfGridPos gp_p, gp_lat, gp_lon;
430 gp_lat,
431 gp_lon,
432 jacobian_quantities[q],
433 atmosphere_dim,
434 p_grid,
435 lat_grid,
436 lon_grid);
437
438 Tensor3 mag_x;
440 mag_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
441 flat(xa[ind], mag_x);
442 }
443 }
444
445 // Surface
446 else if (jacobian_quantities[q] == Jacobian::Special::SurfaceString) {
447 surface_props_check(atmosphere_dim,
448 lat_grid,
449 lon_grid,
450 surface_props_data,
451 surface_props_names);
452 ARTS_USER_ERROR_IF (surface_props_data.empty(),
453 "One jacobian quantity belongs to the "
454 "surface category, but *surface_props_data* is empty.");
455
456 const Index isu =
457 find_first(surface_props_names, jacobian_quantities[q].Subtag());
458 ARTS_USER_ERROR_IF (isu < 0,
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*.")
463
464 ArrayOfGridPos gp_lat, gp_lon;
466 gp_lon,
467 jacobian_quantities[q],
468 atmosphere_dim,
469 lat_grid,
470 lon_grid);
471 Matrix surf_x;
473 atmosphere_dim,
474 surface_props_data(isu, joker, joker),
475 gp_lat,
476 gp_lon);
477 flat(xa[ind], surf_x);
478 }
479
480 // All variables having zero as a priori
481 // ----------------------------------------------------------------------------
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) {
486 xa[ind] = 0;
487 }
488
489 else {
491 "Found a retrieval quantity that is not yet handled by\n"
492 "internal retrievals: ", jacobian_quantities[q], '\n')
493 }
494 }
495
496 // Apply transformations
497 transform_x(xa, jacobian_quantities);
498}
499
500/* Workspace method: Doxygen documentation will be auto-generated */
502 Tensor4& vmr_field,
503 Tensor3& t_field,
504 Tensor4& particle_bulkprop_field,
505 Tensor3& wind_u_field,
506 Tensor3& wind_v_field,
507 Tensor3& wind_w_field,
508 Tensor3& mag_u_field,
509 Tensor3& mag_v_field,
510 Tensor3& mag_w_field,
511 Tensor3& surface_props_data,
512 const ArrayOfRetrievalQuantity& jacobian_quantities,
513 const Vector& x,
514 const Index& atmfields_checked,
515 const Index& atmgeom_checked,
516 const Index& atmosphere_dim,
517 const Vector& p_grid,
518 const Vector& lat_grid,
519 const Vector& lon_grid,
520 const ArrayOfArrayOfSpeciesTag& abs_species,
521 const Index& cloudbox_on,
522 const Index& cloudbox_checked,
523 const ArrayOfString& particle_bulkprop_names,
524 const ArrayOfString& surface_props_names,
525 const Agenda& water_p_eq_agenda,
526 const Verbosity&) {
527 // Basics
528 //
529 ARTS_USER_ERROR_IF (atmfields_checked != 1,
530 "The atmospheric fields must be flagged to have "
531 "passed a consistency check (atmfields_checked=1).");
532 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
533 "The atmospheric geometry must be flagged to have "
534 "passed a consistency check (atmgeom_checked=1).");
535 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
536 "The cloudbox must be flagged to have "
537 "passed a consistency check (cloudbox_checked=1).");
538
539 // Revert transformation
540 Vector x_t(x);
541 transform_x_back(x_t, jacobian_quantities);
542
543 // Main sizes
544 const Index nq = jacobian_quantities.nelem();
545
546 // Jacobian indices
548 {
549 bool any_affine;
550 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
551 }
552
553 // Check input
554 ARTS_USER_ERROR_IF (x_t.nelem() != ji[nq - 1][1] + 1,
555 "Length of *x* does not match length implied by "
556 "*jacobian_quantities*.");
557
558 // Note that when this method is called, vmr_field and other output variables
559 // have original values, i.e. matching the a priori state.
560
561 // Loop retrieval quantities
562 for (Index q = 0; q < nq; q++) {
563 // Index range of this retrieval quantity
564 const Index np = ji[q][1] - ji[q][0] + 1;
565 Range ind(ji[q][0], np);
566
567 // Atmospheric temperatures
568 // ----------------------------------------------------------------------------
569 if (jacobian_quantities[q] == Jacobian::Atm::Temperature) {
570 // Determine grid positions for interpolation from retrieval grids back
571 // to atmospheric grids
572 ArrayOfGridPos gp_p, gp_lat, gp_lon;
573 Index n_p, n_lat, n_lon;
575 gp_lat,
576 gp_lon,
577 n_p,
578 n_lat,
579 n_lon,
580 jacobian_quantities[q].Grids(),
581 atmosphere_dim,
582 p_grid,
583 lat_grid,
584 lon_grid);
585
586 // Map values in x back to t_field
587 Tensor3 t_x(n_p, n_lat, n_lon);
588 reshape(t_x, x_t[ind]);
590 t_field, atmosphere_dim, t_x, gp_p, gp_lat, gp_lon);
591 }
592
593 // Abs species
594 // ----------------------------------------------------------------------------
595 else if (jacobian_quantities[q].Target().isSpeciesVMR()) {
596 // Index position of species
597 ArrayOfSpeciesTag atag(jacobian_quantities[q].Subtag());
598 const Index isp = chk_contains("abs_species", abs_species, atag);
599
600 // Map part of x to a full atmospheric field
601 Tensor3 x_field(vmr_field.npages(), vmr_field.nrows(), vmr_field.ncols());
602 {
603 ArrayOfGridPos gp_p, gp_lat, gp_lon;
604 Index n_p, n_lat, n_lon;
606 gp_lat,
607 gp_lon,
608 n_p,
609 n_lat,
610 n_lon,
611 jacobian_quantities[q].Grids(),
612 atmosphere_dim,
613 p_grid,
614 lat_grid,
615 lon_grid);
616 //
617 Tensor3 t3_x(n_p, n_lat, n_lon);
618 reshape(t3_x, x_t[ind]);
620 x_field, atmosphere_dim, t3_x, gp_p, gp_lat, gp_lon);
621 }
622 //
623 if (jacobian_quantities[q].Mode() == "rel") {
624 // vmr = vmr0 * x
625 vmr_field(isp, joker, joker, joker) *= x_field;
626 } else if (jacobian_quantities[q].Mode() == "vmr") {
627 // vmr = x
628 vmr_field(isp, joker, joker, joker) = x_field;
629 } else if (jacobian_quantities[q].Mode() == "nd") {
630 // vmr = nd / nd_tot
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) /
636 number_density(p_grid[i1], t_field(i1, i2, i3));
637 }
638 }
639 }
640 } else if (jacobian_quantities[q].Mode() == "rh") {
641 // vmr = x * p_sat / p
642 Tensor3 water_p_eq;
643 water_p_eq_agendaExecute(ws, water_p_eq, t_field, water_p_eq_agenda);
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];
649 }
650 }
651 }
652 } else if (jacobian_quantities[q].Mode() == "q") {
653 // We have that specific humidity q, mixing ratio r and
654 // vapour pressure e, are related as
655 // q = r(1+r); r = 0.622e/(p-e); e = vmr*p;
656 // That is: vmr=e/p; e = rp/(0.622+r); r = q/(1-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];
663 }
664 }
665 }
666 } else {
667 ARTS_ASSERT(0);
668 }
669 }
670
671 // Scattering species
672 // ----------------------------------------------------------------------------
673 else if (jacobian_quantities[q] == Jacobian::Special::ScatteringString) {
674 // If no cloudbox, we assume that there is nothing to do
675 if (cloudbox_on) {
676 ARTS_USER_ERROR_IF (particle_bulkprop_field.empty(),
677 "One jacobian quantity belongs to the "
678 "scattering species category, but *particle_bulkprop_field* "
679 "is empty.");
680 ARTS_USER_ERROR_IF (particle_bulkprop_field.nbooks() !=
681 particle_bulkprop_names.nelem(),
682 "Mismatch in size between "
683 "*particle_bulkprop_field* and *particle_bulkprop_field*.");
684
685 const Index isp = find_first(particle_bulkprop_names,
686 jacobian_quantities[q].SubSubtag());
687 ARTS_USER_ERROR_IF (isp < 0,
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*.")
692
693 // Determine grid positions for interpolation from retrieval grids back
694 // to atmospheric grids
695 ArrayOfGridPos gp_p, gp_lat, gp_lon;
696 Index n_p, n_lat, n_lon;
698 gp_lat,
699 gp_lon,
700 n_p,
701 n_lat,
702 n_lon,
703 jacobian_quantities[q].Grids(),
704 atmosphere_dim,
705 p_grid,
706 lat_grid,
707 lon_grid);
708 // Map x to particle_bulkprop_field
709 Tensor3 pbfield_x(n_p, n_lat, n_lon);
710 reshape(pbfield_x, x_t[ind]);
711 Tensor3 pbfield;
713 pbfield, atmosphere_dim, pbfield_x, gp_p, gp_lat, gp_lon);
714 particle_bulkprop_field(isp, joker, joker, joker) = pbfield;
715 }
716 }
717
718 // Wind
719 // ----------------------------------------------------------------------------
720 else if (jacobian_quantities[q].Target().isWind()) {
721 // Determine grid positions for interpolation from retrieval grids back
722 // to atmospheric grids
723 ArrayOfGridPos gp_p, gp_lat, gp_lon;
724 Index n_p, n_lat, n_lon;
726 gp_lat,
727 gp_lon,
728 n_p,
729 n_lat,
730 n_lon,
731 jacobian_quantities[q].Grids(),
732 atmosphere_dim,
733 p_grid,
734 lat_grid,
735 lon_grid);
736
737 // TODO Could be done without copying.
738 Tensor3 wind_x(n_p, n_lat, n_lon);
739 reshape(wind_x, x_t[ind]);
740
741 Tensor3View target_field(wind_u_field);
742
743 Tensor3 wind_field(
744 target_field.npages(), target_field.nrows(), target_field.ncols());
746 wind_field, atmosphere_dim, wind_x, gp_p, gp_lat, gp_lon);
747
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;
754 }
755 }
756
757 // Magnetism
758 // ----------------------------------------------------------------------------
759 else if (jacobian_quantities[q].Target().isMagnetic()) {
760 // Determine grid positions for interpolation from retrieval grids back
761 // to atmospheric grids
762 ArrayOfGridPos gp_p, gp_lat, gp_lon;
763 Index n_p, n_lat, n_lon;
765 gp_lat,
766 gp_lon,
767 n_p,
768 n_lat,
769 n_lon,
770 jacobian_quantities[q].Grids(),
771 atmosphere_dim,
772 p_grid,
773 lat_grid,
774 lon_grid);
775
776 // TODO Could be done without copying.
777 Tensor3 mag_x(n_p, n_lat, n_lon);
778 reshape(mag_x, x_t[ind]);
779
780 Tensor3View target_field(mag_u_field);
781
782 Tensor3 mag_field(
783 target_field.npages(), target_field.nrows(), target_field.ncols());
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;
801 }
802 }
803 }
804 } else {
805 ARTS_USER_ERROR ( "Unsupported magnetism type");
806 }
807 }
808
809 // Surface
810 // ----------------------------------------------------------------------------
811 else if (jacobian_quantities[q] == Jacobian::Special::SurfaceString) {
812 surface_props_check(atmosphere_dim,
813 lat_grid,
814 lon_grid,
815 surface_props_data,
816 surface_props_names);
817 ARTS_USER_ERROR_IF (surface_props_data.empty(),
818 "One jacobian quantity belongs to the "
819 "surface category, but *surface_props_data* is empty.");
820
821 const Index isu =
822 find_first(surface_props_names, jacobian_quantities[q].Subtag());
823 ARTS_USER_ERROR_IF (isu < 0,
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*.")
828
829 // Determine grid positions for interpolation from retrieval grids back
830 // to atmospheric grids
831 ArrayOfGridPos gp_lat, gp_lon;
832 Index n_lat, n_lon;
834 gp_lon,
835 n_lat,
836 n_lon,
837 jacobian_quantities[q].Grids(),
838 atmosphere_dim,
839 lat_grid,
840 lon_grid);
841 // Map values in x back to surface_props_data
842 Matrix surf_x(n_lat, n_lon);
843 reshape(surf_x, x_t[ind]);
844 Matrix surf;
845 regrid_atmsurf_by_gp_oem(surf, atmosphere_dim, surf_x, gp_lat, gp_lon);
846 surface_props_data(isu, joker, joker) = surf;
847 }
848 }
849}
850
851/* Workspace method: Doxygen documentation will be auto-generated */
853 Matrix& sensor_los,
854 Vector& f_backend,
855 Vector& y_baseline,
856 Sparse& sensor_response,
857 Vector& sensor_response_f,
858 ArrayOfIndex& sensor_response_pol,
859 Matrix& sensor_response_dlos,
860 Vector& sensor_response_f_grid,
861 ArrayOfIndex& sensor_response_pol_grid,
862 Matrix& sensor_response_dlos_grid,
863 Matrix& mblock_dlos,
864 const ArrayOfRetrievalQuantity& jacobian_quantities,
865 const Vector& x,
866 const Agenda& sensor_response_agenda,
867 const Index& sensor_checked,
868 const ArrayOfTime& sensor_time,
869 const Verbosity&) {
870 // Basics
871 //
872 ARTS_USER_ERROR_IF (sensor_checked != 1,
873 "The sensor response must be flagged to have "
874 "passed a consistency check (sensor_checked=1).");
875
876 // Revert transformation
877 Vector x_t(x);
878 transform_x_back(x_t, jacobian_quantities);
879
880 // Main sizes
881 const Index nq = jacobian_quantities.nelem();
882
883 // Jacobian indices
885 {
886 bool any_affine;
887 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
888 }
889
890 // Check input
891 ARTS_USER_ERROR_IF (x_t.nelem() != ji[nq - 1][1] + 1,
892 "Length of *x* does not match length implied by "
893 "*jacobian_quantities*.");
894
895 // Flag indicating that y_baseline is not set
896 bool yb_set = false;
897
898 // Shall sensor responses be calculed?
899 bool do_sensor = false;
900
901 // Loop retrieval quantities
902 for (Index q = 0; q < nq; q++) {
903 // Index range of this retrieval quantity
904 const Index np = ji[q][1] - ji[q][0] + 1;
905
906 // Pointing off-set
907 // ----------------------------------------------------------------------------
908 if (jacobian_quantities[q].Target().isPointing()) {
909 // Handle pointing "jitter" seperately
910 if (jacobian_quantities[q].Grids()[0][0] == -1) {
911 ARTS_USER_ERROR_IF (sensor_los.nrows() != np,
912 "Mismatch between pointing jacobian and *sensor_los* found.");
913 // Simply add retrieved off-set(s) to za column of *sensor_los*
914 for (Index i = 0; i < np; i++) {
915 sensor_los(i, 0) += x_t[ji[q][0] + i];
916 }
917 }
918 // Polynomial representation
919 else {
920 ARTS_USER_ERROR_IF (sensor_los.nrows() != sensor_time.nelem(),
921 "Sizes of *sensor_los* and *sensor_time* do not match.");
922 Vector w;
923 for (Index c = 0; c < np; c++) {
924 polynomial_basis_func(w, time_vector(sensor_time), c);
925 for (Index i = 0; i < w.nelem(); i++) {
926 sensor_los(i, 0) += w[i] * x_t[ji[q][0] + c];
927 }
928 }
929 }
930 }
931
932 // Frequency shift or stretch
933 // ----------------------------------------------------------------------------
934 else if (jacobian_quantities[q].Target().isFrequency()) {
935 if (jacobian_quantities[q] == Jacobian::Sensor::FrequencyShift) {
936 ARTS_ASSERT(np == 1);
937 if (x_t[ji[q][0]] != 0) {
938 do_sensor = true;
939 f_backend += x_t[ji[q][0]];
940 }
941 } else if (jacobian_quantities[q] == Jacobian::Sensor::FrequencyStretch) {
942 ARTS_ASSERT(np == 1);
943 if (x_t[ji[q][0]] != 0) {
944 do_sensor = true;
945 Vector w;
946 polynomial_basis_func(w, f_backend, 1);
947 for (Index i = 0; i < w.nelem(); i++) {
948 f_backend[i] += w[i] * x_t[ji[q][0]];
949 }
950 }
951 } else {
952 ARTS_ASSERT(0);
953 }
954 }
955
956 // Baseline fit: polynomial or sinusoidal
957 // ----------------------------------------------------------------------------
958 else if (jacobian_quantities[q] == Jacobian::Sensor::Polyfit ||
959 jacobian_quantities[q] == Jacobian::Sensor::Sinefit) {
960 if (!yb_set) {
961 yb_set = true;
962 Index y_size = sensor_los.nrows() * sensor_response_f_grid.nelem() *
963 sensor_response_pol_grid.nelem() *
964 sensor_response_dlos_grid.nrows();
965 y_baseline.resize(y_size);
966 y_baseline = 0;
967 }
968
969 for (Index mb = 0; mb < sensor_los.nrows(); ++mb) {
970 calcBaselineFit(y_baseline,
971 x_t,
972 mb,
973 sensor_response,
974 sensor_response_pol_grid,
975 sensor_response_f_grid,
976 sensor_response_dlos_grid,
977 jacobian_quantities[q],
978 q,
979 ji);
980 }
981 }
982 }
983
984 // *y_baseline* not yet set?
985 if (!yb_set) {
986 y_baseline.resize(1);
987 y_baseline[0] = 0;
988 }
989
990 // Recalculate sensor_response?
991 if (do_sensor) {
993 sensor_response,
994 sensor_response_f,
995 sensor_response_f_grid,
996 sensor_response_pol,
997 sensor_response_pol_grid,
998 sensor_response_dlos,
999 sensor_response_dlos_grid,
1000 mblock_dlos,
1001 f_backend,
1002 sensor_response_agenda);
1003 }
1004}
1005
1006/* Workspace method: Doxygen documentation will be auto-generated */
1008 ARTS_USER_ERROR ( "Retrievals of spectroscopic variables not yet handled.");
1009}
1010
1011
1012#ifdef OEM_SUPPORT
1013/* Workspace method: Doxygen documentation will be auto-generated */
1014void OEM(Workspace& ws,
1015 Vector& x,
1016 Vector& yf,
1017 Matrix& jacobian,
1018 Matrix& dxdy,
1019 Vector& oem_diagnostics,
1020 Vector& lm_ga_history,
1021 ArrayOfString& errors,
1022 const Vector& xa,
1023 const CovarianceMatrix& covmat_sx,
1024 const Vector& y,
1025 const CovarianceMatrix& covmat_se,
1026 const ArrayOfRetrievalQuantity& jacobian_quantities,
1027 const Agenda& inversion_iterate_agenda,
1028 const String& method,
1029 const Numeric& max_start_cost,
1030 const Vector& x_norm,
1031 const Index& max_iter,
1032 const Numeric& stop_dx,
1033 const Vector& lm_ga_settings,
1034 const Index& clear_matrices,
1035 const Index& display_progress,
1036 const Verbosity&) {
1037 // Main sizes
1038 const Index n = covmat_sx.nrows();
1039 const Index m = y.nelem();
1040
1041 // Checks
1042 covmat_sx.compute_inverse();
1043 covmat_se.compute_inverse();
1044
1045 OEM_checks(ws,
1046 x,
1047 yf,
1048 jacobian,
1049 inversion_iterate_agenda,
1050 xa,
1051 covmat_sx,
1052 y,
1053 covmat_se,
1054 jacobian_quantities,
1055 method,
1056 x_norm,
1057 max_iter,
1058 stop_dx,
1059 lm_ga_settings,
1060 clear_matrices,
1061 display_progress);
1062
1063 // Size diagnostic output and init with NaNs
1064 oem_diagnostics.resize(5);
1065 oem_diagnostics = NAN;
1066 //
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;
1071 } else {
1072 lm_ga_history.resize(0);
1073 }
1074
1075 // Check for start vector and precomputed yf, jacobian
1076 if (x.nelem() != n) {
1077 x = xa;
1078 yf.resize(0);
1079 jacobian.resize(0, 0);
1080 }
1081
1082 // If no precomputed value given, we compute yf and jacobian to
1083 // compute initial cost (and use in the first OEM iteration).
1084 if (yf.nelem() == 0) {
1086 ws, yf, jacobian, xa, 1, 0, inversion_iterate_agenda);
1087 }
1088
1089 ARTS_USER_ERROR_IF (yf.nelem() not_eq y.nelem(),
1090 "Mismatch between simulated y and input y.\n"
1091 "Input y is size ", y.nelem(), " but simulated y is ",
1092 yf.nelem(), "\n"
1093 "Use your frequency grid vector and your sensor response matrix to match simulations with measurements.\n")
1094
1095 // TODO: Get this from invlib log.
1096 // Start value of cost function
1097 Numeric cost_start = NAN;
1098 if (method == "ml" || method == "lm" || display_progress ||
1099 max_start_cost > 0) {
1100 Vector dy = y;
1101 dy -= yf;
1102 Vector sdy = y;
1103 mult_inv(sdy, covmat_se, dy);
1104 Vector dx = x;
1105 dx -= xa;
1106 Vector sdx = x;
1107 mult_inv(sdx, covmat_sx, dx);
1108 cost_start = dx * sdx + dy * sdy;
1109 cost_start /= static_cast<Numeric>(m);
1110 }
1111 oem_diagnostics[1] = cost_start;
1112
1113 // Handle cases with too large start cost
1114 if (max_start_cost > 0 && cost_start > max_start_cost) {
1115 // Flag no inversion in oem_diagnostics, and let x to be undefined
1116 oem_diagnostics[0] = 99;
1117 //
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
1122 << endl;
1123 }
1124 }
1125 // Otherwise do inversion
1126 else {
1127 bool apply_norm = false;
1128 oem::Matrix T{};
1129 if (x_norm.nelem() == n) {
1130 T.resize(n, n);
1131 T *= 0.0;
1132 T.diagonal() = x_norm;
1133 for (Index i = 0; i < n; i++) {
1134 T(i, i) = x_norm[i];
1135 }
1136 apply_norm = true;
1137 }
1138
1139 oem::CovarianceMatrix Se(covmat_se), Sa(covmat_sx);
1140 oem::Vector xa_oem(xa), y_oem(y), x_oem(x);
1141 oem::AgendaWrapper aw(&ws,
1142 (unsigned int)m,
1143 (unsigned int)n,
1144 jacobian,
1145 yf,
1146 &inversion_iterate_agenda);
1147 oem::OEM_STANDARD<oem::AgendaWrapper> oem(aw, xa_oem, Sa, Se);
1148 oem::OEM_MFORM<oem::AgendaWrapper> oem_m(aw, xa_oem, Sa, Se);
1149 int oem_verbosity = static_cast<int>(display_progress);
1150
1151 int return_code = 0;
1152
1153 try {
1154 if (method == "li") {
1155 oem::Std s(T, apply_norm);
1156 oem::GN gn(stop_dx, 1, s); // Linear case, only one step.
1157 return_code = oem.compute<oem::GN, oem::ArtsLog>(
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") {
1161 oem::Std s(T, apply_norm);
1162 oem::GN gn(stop_dx, 1, s); // Linear case, only one step.
1163 return_code = oem_m.compute<oem::GN, oem::ArtsLog>(
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);
1168 oem::GN_CG gn(stop_dx, 1, cg); // Linear case, only one step.
1169 return_code = oem.compute<oem::GN_CG, oem::ArtsLog>(
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);
1174 oem::GN_CG gn(stop_dx, 1, cg); // Linear case, only one step.
1175 return_code = oem_m.compute<oem::GN_CG, oem::ArtsLog>(
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") {
1179 oem::Std s(T, apply_norm);
1180 oem::GN gn(stop_dx, (unsigned int)max_iter, s);
1181 return_code = oem.compute<oem::GN, oem::ArtsLog>(
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") {
1185 oem::Std s(T, apply_norm);
1186 oem::GN gn(stop_dx, (unsigned int)max_iter, s);
1187 return_code = oem_m.compute<oem::GN, oem::ArtsLog>(
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);
1193 return_code = oem.compute<oem::GN_CG, oem::ArtsLog>(
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);
1199 return_code = oem_m.compute<oem::GN_CG, oem::ArtsLog>(
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")) {
1203 oem::Std s(T, apply_norm);
1204 Sparse diagonal = Sparse::diagonal(covmat_sx.inverse_diagonal());
1205 CovarianceMatrix SaDiag{};
1206 SaDiag.add_correlation_inverse(Block(Range(0, n),
1207 Range(0, n),
1208 std::make_pair(0, 0),
1209 make_shared<Sparse>(diagonal)));
1211 oem::LM lm(SaInvLM, s);
1212
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]);
1221
1222 return_code = oem.compute<oem::LM&, oem::ArtsLog>(
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;
1227 }
1228 } else if ((method == "lm_cg") || (method == "ml_cg")) {
1229 oem::CG cg(T, apply_norm, 1e-10, 0);
1230
1231 Sparse diagonal = Sparse::diagonal(covmat_sx.inverse_diagonal());
1232 CovarianceMatrix SaDiag{};
1233 SaDiag.add_correlation_inverse(Block(Range(0, n),
1234 Range(0, n),
1235 std::make_pair(0, 0),
1236 make_shared<Sparse>(diagonal)));
1237 oem::LM_CG lm(SaDiag, cg);
1238
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]);
1245
1246 return_code = oem.compute<oem::LM_CG&, oem::ArtsLog>(
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;
1251 }
1252 }
1253
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);
1262 x_oem *= NAN;
1263 std::vector<std::string> sv = oem::handle_nested_exception(e);
1264 for (auto& s : sv) {
1265 std::stringstream ss{s};
1266 std::string t{};
1267 while (std::getline(ss, t)) {
1268 errors.push_back(t.c_str());
1269 }
1270 }
1271 } catch (...) {
1272 throw;
1273 }
1274
1275 x = x_oem;
1276 yf = aw.get_measurement_vector();
1277
1278 // Shall empty jacobian and dxdy be returned?
1279 if (clear_matrices) {
1280 jacobian.resize(0, 0);
1281 dxdy.resize(0, 0);
1282 } else if (oem_diagnostics[0] <= 2) {
1283 dxdy.resize(n, m);
1284 Matrix tmp1(n, m), tmp2(n, n), tmp3(n, n);
1285 mult_inv(tmp1, transpose(jacobian), covmat_se);
1286 mult(tmp2, tmp1, jacobian);
1287 add_inv(tmp2, covmat_sx);
1288 inv(tmp3, tmp2);
1289 mult(dxdy, tmp3, tmp1);
1290 }
1291 }
1292}
1293
1294/* Workspace method: Doxygen documentation will be auto-generated */
1295void covmat_soCalc(Matrix& covmat_so,
1296 const Matrix& dxdy,
1297 const CovarianceMatrix& covmat_se,
1298 const Verbosity& /*v*/) {
1299 Index n(dxdy.nrows()), m(dxdy.ncols());
1300 Matrix tmp1(m, n);
1301
1302 ARTS_USER_ERROR_IF ((m == 0) || (n == 0),
1303 "The gain matrix *dxdy* is required to compute the observation error covariance matrix.");
1304 ARTS_USER_ERROR_IF ((covmat_se.nrows() != m) || (covmat_se.ncols() != m),
1305 "The covariance matrix covmat_se has invalid dimensions.");
1306
1307 covmat_so.resize(n, n);
1308 mult(tmp1, covmat_se, transpose(dxdy));
1309 mult(covmat_so, dxdy, tmp1);
1310}
1311
1312/* Workspace method: Doxygen documentation will be auto-generated */
1313void covmat_ssCalc(Matrix& covmat_ss,
1314 const Matrix& avk,
1315 const CovarianceMatrix& covmat_sx,
1316 const Verbosity& /*v*/) {
1317 Index n(avk.ncols());
1318 Matrix tmp1(n, n), tmp2(n, n);
1319
1320 ARTS_USER_ERROR_IF (n == 0,
1321 "The averaging kernel matrix *dxdy* is required to compute the smoothing error covariance matrix.");
1322 ARTS_USER_ERROR_IF ((covmat_sx.nrows() != n) || (covmat_sx.ncols() != n),
1323 "The covariance matrix *covmat_sx* invalid dimensions.");
1324
1325 covmat_ss.resize(n, n);
1326
1327 // Sign doesn't matter since we're dealing with a quadratic form.
1328 id_mat(tmp1);
1329 tmp1 -= avk;
1330
1331 mult(tmp2, covmat_sx, transpose(tmp1));
1332 mult(covmat_ss, tmp1, tmp2);
1333}
1334
1335/* Workspace method: Doxygen documentation will be auto-generated */
1337 const CovarianceMatrix& Sc,
1338 const Verbosity& /*v*/) {
1339 S = Matrix(Sc);
1340}
1341
1342/* Workspace method: Doxygen documentation will be auto-generated */
1343void avkCalc(Matrix& avk,
1344 const Matrix& dxdy,
1345 const Matrix& jacobian,
1346 const Verbosity& /*v*/) {
1347 Index m(jacobian.nrows()), n(jacobian.ncols());
1348 ARTS_USER_ERROR_IF ((m == 0) || (n == 0),
1349 "The Jacobian matrix is empty.");
1350 ARTS_USER_ERROR_IF ((dxdy.nrows() != n) || (dxdy.ncols() != m),
1351 "Matrices have inconsistent sizes.\n"
1352 " Size of gain matrix: ", dxdy.nrows(), " x ", dxdy.ncols(),
1353 "\n"
1354 " Size of Jacobian: ", jacobian.nrows(), " x ",
1355 jacobian.ncols(), "\n")
1356
1357 avk.resize(n, n);
1358 mult(avk, dxdy, jacobian);
1359}
1360
1361#else
1362
1363void covmat_soCalc(Matrix& /* covmat_so */,
1364 const Matrix& /* dxdy */,
1365 const CovarianceMatrix& /* covmat_se*/,
1366 const Verbosity& /*v*/) {
1368 "WSM is not available because ARTS was compiled without "
1369 "OEM support.");
1370}
1371
1372void covmat_ssCalc(Matrix& /*covmat_ss*/,
1373 const Matrix& /*avk*/,
1374 const CovarianceMatrix& /*covmat_sx*/,
1375 const Verbosity& /*v*/) {
1377 "WSM is not available because ARTS was compiled without "
1378 "OEM support.");
1379}
1380
1381void avkCalc(Matrix& /* avk */,
1382 const Matrix& /* dxdy */,
1383 const Matrix& /* jacobian */,
1384 const Verbosity& /*v*/) {
1386 "WSM is not available because ARTS was compiled without "
1387 "OEM support.");
1388}
1389
1391 Vector&,
1392 Vector&,
1393 Matrix&,
1394 Matrix&,
1395 Vector&,
1396 Vector&,
1398 const Vector&,
1399 const CovarianceMatrix&,
1400 const Vector&,
1401 const CovarianceMatrix&,
1402 const Index&,
1404 const ArrayOfArrayOfIndex&,
1405 const Agenda&,
1406 const String&,
1407 const Numeric&,
1408 const Vector&,
1409 const Index&,
1410 const Numeric&,
1411 const Vector&,
1412 const Index&,
1413 const Index&,
1414 const Verbosity&) {
1416 "WSM is not available because ARTS was compiled without "
1417 "OEM support.");
1418}
1419#endif
1420
1421#if defined(OEM_SUPPORT) && 0
1422
1423#include "agenda_wrapper_mpi.h"
1424#include "oem_mpi.h"
1425
1426//
1427// Performs manipulations of workspace variables necessary for distributed
1428// retrievals with MPI:
1429//
1430// - Splits up sensor positions evenly over processes
1431// - Splits up inverse covariance matrices.
1432//
1433void MPI_Initialize(Matrix& sensor_los,
1434 Matrix& sensor_pos,
1435 Vector& sensor_time) {
1436 int initialized;
1437
1438 MPI_Initialized(&initialized);
1439 if (!initialized) {
1440 MPI_Init(nullptr, nullptr);
1441 }
1442
1443 int rank, nprocs;
1444
1445 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1446 MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
1447
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);
1452
1453 //
1454 // Split up sensor positions.
1455 //
1456
1457 if (rank < remainder) {
1458 mblock_range += 1;
1459 mblock_start += rank;
1460 } else {
1461 mblock_start += remainder;
1462 }
1463
1464 if (nmblock > 0) {
1465 Range range = Range(mblock_start, mblock_range);
1466
1467 Matrix tmp_m = sensor_los(range, joker);
1468 sensor_los = tmp_m;
1469
1470 tmp_m = sensor_pos(range, joker);
1471 sensor_pos = tmp_m;
1472
1473 Vector tmp_v = sensor_time[range];
1474 sensor_time = tmp_v;
1475 } else {
1476 sensor_los.resize(0, 0);
1477 sensor_pos.resize(0, 0);
1478 sensor_time.resize(0);
1479 }
1480}
1481
1482void OEM_MPI(Workspace& ws,
1483 Vector& x,
1484 Vector& yf,
1485 Matrix& jacobian,
1486 Matrix& dxdy,
1487 Vector& oem_diagnostics,
1488 Vector& lm_ga_history,
1489 Matrix& sensor_los,
1490 Matrix& sensor_pos,
1491 Vector& sensor_time,
1492 CovarianceMatrix& covmat_sx,
1493 CovarianceMatrix& covmat_se,
1494 const Vector& xa,
1495 const Vector& y,
1496 const ArrayOfRetrievalQuantity& jacobian_quantities,
1497 const Agenda& inversion_iterate_agenda,
1498 const String& method,
1499 const Numeric& max_start_cost,
1500 const Vector& x_norm,
1501 const Index& max_iter,
1502 const Numeric& stop_dx,
1503 const Vector& lm_ga_settings,
1504 const Index& clear_matrices,
1505 const Index& display_progress,
1506 const Verbosity& /*v*/) {
1507 // Main sizes
1508 const Index n = covmat_sx.nrows();
1509 const Index m = y.nelem();
1510
1511 // Check WSVs
1512 OEM_checks(ws,
1513 x,
1514 yf,
1515 jacobian,
1516 inversion_iterate_agenda,
1517 xa,
1518 covmat_sx,
1519 covmat_se,
1520 jacobian_quantities,
1521 method,
1522 x_norm,
1523 max_iter,
1524 stop_dx,
1525 lm_ga_settings,
1526 clear_matrices,
1527 display_progress);
1528
1529 // Calculate spectrum and Jacobian for a priori state
1530 // Jacobian is also input to the agenda, and to flag this is this first
1531 // call, this WSV must be set to be empty.
1532 jacobian.resize(0, 0);
1533
1534 // Initialize MPI environment.
1535 MPI_Initialize(sensor_los, sensor_pos, sensor_time);
1536
1537 // Setup distributed matrices.
1538 MPICovarianceMatrix SeInvMPI(covmat_se);
1539 MPICovarianceMatrix SaInvMPI(covmat_sx);
1540
1541 // Create temporary MPI vector from local results and use conversion to
1542 // standard vector to broadcast results to all processes.
1543 oem::Vector tmp;
1545 ws, tmp, jacobian, xa, 1, inversion_iterate_agenda);
1546 yf = MPIVector(tmp);
1547
1548 // Size diagnostic output and init with NaNs
1549 oem_diagnostics.resize(5);
1550 oem_diagnostics = NAN;
1551 //
1552 if (method == "ml" || method == "lm") {
1553 lm_ga_history.resize(max_iter);
1554 lm_ga_history = NAN;
1555 } else {
1556 lm_ga_history.resize(0);
1557 }
1558
1559 // Start value of cost function. Covariance matrices are already distributed
1560 // over processes, so we need to use invlib matrix algebra.
1561 Numeric cost_start = NAN;
1562 if (method == "ml" || method == "lm" || display_progress ||
1563 max_start_cost > 0) {
1564 oem::Vector dy = y;
1565 dy -= yf;
1566 cost_start = dot(dy, SeInvMPI * dy);
1567 }
1568 oem_diagnostics[1] = cost_start;
1569
1570 // Handle cases with too large start cost
1571 if (max_start_cost > 0 && cost_start > max_start_cost) {
1572 // Flag no inversion in oem_diagnostics, and let x to be undefined
1573 oem_diagnostics[0] = 99;
1574 //
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
1579 << endl;
1580 }
1581 }
1582
1583 // Otherwise do inversion
1584 else {
1585 // Size remaining output arguments
1586 x.resize(n);
1587 dxdy.resize(n, m);
1588
1589 OEMVector xa_oem(xa), y_oem(y), x_oem;
1590 oem::AgendaWrapperMPI aw(&ws, &inversion_iterate_agenda, m, n);
1591
1592 OEM_PS_PS_MPI<AgendaWrapperMPI> oem(aw, xa_oem, SaInvMPI, SeInvMPI);
1593
1594 // Call selected method
1595 int return_value = 99;
1596
1597 if (method == "li") {
1598 oem::CG cg(1e-12, 0);
1599 oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1600 return_value = oem.compute<oem::GN_CG, invlib::MPILog>(
1601 x_oem, y_oem, gn, 2 * (int)display_progress);
1602 } else if (method == "gn") {
1603 oem::CG cg(1e-12, 0);
1604 oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1605 return_value = oem.compute<oem::GN_CG, invlib::MPILog>(
1606 x_oem, y_oem, gn, 2 * (int)display_progress);
1607 } else if ((method == "lm") || (method == "ml")) {
1608 oem::CG cg(1e-12, 0);
1609 LM_CG_MPI lm(SaInvMPI, cg);
1610
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]);
1618
1619 return_value = oem.compute<oem::LM_CG_MPI, invlib::MPILog>(
1620 x_oem, y_oem, lm, 2 * (int)display_progress);
1621 }
1622
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);
1627
1628 x = x_oem;
1629 // Shall empty jacobian and dxdy be returned?
1630 if (clear_matrices && (oem_diagnostics[0])) {
1631 jacobian.resize(0, 0);
1632 dxdy.resize(0, 0);
1633 }
1634 }
1635 MPI_Finalize();
1636}
1637
1638#else
1639
1641 Vector&,
1642 Vector&,
1643 Matrix&,
1644 Matrix&,
1645 Vector&,
1646 Vector&,
1647 Matrix&,
1648 Matrix&,
1649 Vector&,
1652 const Vector&,
1653 const Vector&,
1654 const Index&,
1656 const Agenda&,
1657 const String&,
1658 const Numeric&,
1659 const Vector&,
1660 const Index&,
1661 const Numeric&,
1662 const Vector&,
1663 const Index&,
1664 const Index&,
1665 const Verbosity&) {
1667 "You have to compile ARTS with OEM support "
1668 " and enable MPI to use OEM_MPI.");
1669}
1670
1671#endif // OEM_SUPPORT && ENABLE_MPI
This file contains the definition of Array.
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:190
The global header file for ARTS.
Header file for helper functions for OpenMP.
Vector time_vector(const ArrayOfTime &times)
Converts from each Time to seconds and returns as Vector.
Definition: artstime.cc:167
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:25106
void water_p_eq_agendaExecute(Workspace &ws, Tensor3 &water_p_eq_field, const Tensor3 &t_field, const Agenda &input_agenda)
Definition: auto_md.cc:26138
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:25897
void MatrixFromCovarianceMatrix(Matrix &out, const CovarianceMatrix &in, const Verbosity &verbosity)
WORKSPACE METHOD: MatrixFromCovarianceMatrix.
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:149
The Agenda class.
Definition: agenda_class.h:69
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
A constant view of a Tensor3.
Definition: matpackIII.h:133
bool empty() const noexcept
Definition: matpackIII.h:157
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:145
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:148
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:151
Index ncols() const noexcept
Definition: matpackIV.h:146
bool empty() const noexcept
Definition: matpackIV.h:152
Index nrows() const noexcept
Definition: matpackIV.h:145
Index nbooks() const noexcept
Definition: matpackIV.h:143
Index npages() const noexcept
Definition: matpackIV.h:144
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
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.
The Matrix class.
Definition: matpackI.h:1285
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1011
The range class.
Definition: matpackI.h:160
The Tensor3View class.
Definition: matpackIII.h:251
The Tensor3 class.
Definition: matpackIII.h:352
The Tensor4 class.
Definition: matpackIV.h:435
The Vector class.
Definition: matpackI.h:910
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
Workspace class.
Definition: workspace_ng.h:53
Interface for distributed ARTS forward model.
Definition: oem_mpi.h:56
Interface to ARTS inversion_iterate_agenda.
Definition: oem.h:461
OEM log output.
Definition: oem.h:243
Normalizing solver.
Definition: oem.h:112
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,...)
Definition: debug.h:102
#define ARTS_USER_ERROR(...)
Definition: debug.h:169
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
void polynomial_basis_func(Vector &b, const Vector &x, const Index &poly_coeff)
Calculates polynomial basis functions.
Definition: jacobian.cc:848
void transform_x_back(Vector &x_t, const ArrayOfRetrievalQuantity &jqs, bool revert_functional_transforms)
Handles back-transformations of the state vector.
Definition: jacobian.cc:233
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:46
void transform_x(Vector &x, const ArrayOfRetrievalQuantity &jqs)
Handles transformations of the state vector.
Definition: jacobian.cc:156
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:872
Routines for setting up the jacobian.
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:466
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
void covmat_ssCalc(Matrix &, const Matrix &, const CovarianceMatrix &, const Verbosity &)
WORKSPACE METHOD: covmat_ssCalc.
Definition: m_oem.cc:1372
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:63
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:852
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:1640
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:1390
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:501
void x2artsSpectroscopy(const Verbosity &)
WORKSPACE METHOD: x2artsSpectroscopy.
Definition: m_oem.cc:1007
void avkCalc(Matrix &, const Matrix &, const Matrix &, const Verbosity &)
WORKSPACE METHOD: avkCalc.
Definition: m_oem.cc:1381
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:114
void covmat_soCalc(Matrix &, const Matrix &, const CovarianceMatrix &, const Verbosity &)
WORKSPACE METHOD: covmat_soCalc.
Definition: m_oem.cc:1363
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:89
void reshape(Tensor3View X, ConstVectorView x)
reshape
Definition: math_funcs.cc:782
void flat(VectorView x, ConstMatrixView X)
flat
Definition: math_funcs.cc:728
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
const Joker joker
constexpr Numeric k
Boltzmann constant convenience name [J/K].
Definition: oem.h:28
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::MFORM > OEM_MFORM
OEM m form.
Definition: oem.h:94
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, CG > LM_CG
Levenberg-Marquardt (LM) optimization using normed CG solver.
Definition: oem.h:174
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::STANDARD, invlib::Rodgers531 > OEM_STANDARD
OEM standard form.
Definition: oem.h:62
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, Std > LM
Levenberg-Marquardt (LM) optimization using normed ARTS QR solver.
Definition: oem.h:172
invlib::Matrix< ArtsCovarianceMatrixWrapper > CovarianceMatrix
invlib wrapper type for ARTS the ARTS covariance class.
Definition: oem.h:37
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:33
invlib::GaussNewton< Numeric, CG > GN_CG
Gauss-Newton (GN) optimization using normed CG solver.
Definition: oem.h:170
std::vector< std::string > handle_nested_exception(const E &e, int level=0)
Handle exception encountered within invlib.
Definition: oem.h:424
invlib::GaussNewton< Numeric, Std > GN
OEM Gauss-Newton optimization using normed ARTS QR solver.
Definition: oem.h:168
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:31
invlib::Matrix< invlib::MPIMatrix< invlib::Timer< ArtsCovarianceMatrixWrapper > > > MPICovarianceMatrix
MPI-distributed covariance matrix type.
Definition: oem_mpi.h:27
invlib::Vector< invlib::MPIVector< invlib::Timer< ArtsVector > > > MPIVector
MPI-distributed vector type.
Definition: oem_mpi.h:29
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.
Definition: oem.h:662
void Tensor4Clip(Tensor4 &x, const Index &iq, const Numeric &limit_low, const Numeric &limit_high)
Clip Tensor4.
Definition: oem.h:581
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
Definition: physics_funcs.h:70
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.
Species::Tag SpeciesTag
Definition: species_tags.h:100
The Sparse class.
Definition: matpackII.h:75
Vector diagonal() const
Diagonal elements as vector.
Definition: matpackII.cc:168
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:140
#define w
#define a
#define c