ARTS 2.5.4 (git: 31ce4f0e)
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_grid,
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_grid,
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:24892
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_grid, const Vector &f_backend, const Agenda &input_agenda)
Definition: auto_md.cc:25749
void water_p_eq_agendaExecute(Workspace &ws, Tensor3 &water_p_eq_field, const Tensor3 &t_field, const Agenda &input_agenda)
Definition: auto_md.cc:26034
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:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
A constant view of a Tensor3.
Definition: matpackIII.h:130
bool empty() const noexcept
Definition: matpackIII.h:152
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:140
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:143
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:146
Index ncols() const noexcept
Definition: matpackIV.h:142
bool empty() const noexcept
Definition: matpackIV.h:148
Index nrows() const noexcept
Definition: matpackIV.h:141
Index nbooks() const noexcept
Definition: matpackIV.h:139
Index npages() const noexcept
Definition: matpackIV.h:140
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
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:1261
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1010
The range class.
Definition: matpackI.h:159
The Tensor3View class.
Definition: matpackIII.h:246
The Tensor3 class.
Definition: matpackIII.h:346
The Tensor4 class.
Definition: matpackIV.h:429
The Vector class.
Definition: matpackI.h:899
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:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
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.
#define q
#define dx
#define max(a, b)
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:488
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:168
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_grid, 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:760
void flat(VectorView x, ConstMatrixView X)
flat
Definition: math_funcs.cc:706
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].
constexpr Numeric e
Elementary charge convenience name [C].
Temperature
Definition: jacobian.h:59
Particulates ENUMCLASS(Line, char, VMR, Strength, Center, ShapeG0X0, ShapeG0X1, ShapeG0X2, ShapeG0X3, ShapeD0X0, ShapeD0X1, ShapeD0X2, ShapeD0X3, ShapeG2X0, ShapeG2X1, ShapeG2X2, ShapeG2X3, ShapeD2X0, ShapeD2X1, ShapeD2X2, ShapeD2X3, ShapeFVCX0, ShapeFVCX1, ShapeFVCX2, ShapeFVCX3, ShapeETAX0, ShapeETAX1, ShapeETAX2, ShapeETAX3, ShapeYX0, ShapeYX1, ShapeYX2, ShapeYX3, ShapeGX0, ShapeGX1, ShapeGX2, ShapeGX3, ShapeDVX0, ShapeDVX1, ShapeDVX2, ShapeDVX3, ECS_SCALINGX0, ECS_SCALINGX1, ECS_SCALINGX2, ECS_SCALINGX3, ECS_BETAX0, ECS_BETAX1, ECS_BETAX2, ECS_BETAX3, ECS_LAMBDAX0, ECS_LAMBDAX1, ECS_LAMBDAX2, ECS_LAMBDAX3, ECS_DCX0, ECS_DCX1, ECS_DCX2, ECS_DCX3, NLTE) static_assert(Index(Line ScatteringString
Definition: jacobian.h:103
MagneticMagnitude
Definition: jacobian.h:61
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:101
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