ARTS 2.5.0 (git: 9ee3ac6c)
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
55#ifdef OEM_SUPPORT
56#include "oem.h"
57#endif
58
59/* Workspace method: Doxygen documentation will be auto-generated */
60void particle_bulkprop_fieldClip(Tensor4& particle_bulkprop_field,
61 const ArrayOfString& particle_bulkprop_names,
62 const String& bulkprop_name,
63 const Numeric& limit_low,
64 const Numeric& limit_high,
65 const Verbosity&) {
66 Index iq = -1;
67 if (bulkprop_name == "ALL") {
68 }
69
70 else {
71 for (Index i = 0; i < particle_bulkprop_names.nelem(); i++) {
72 if (particle_bulkprop_names[i] == bulkprop_name) {
73 iq = i;
74 break;
75 }
76 }
77 ARTS_USER_ERROR_IF (iq < 0,
78 "Could not find ", bulkprop_name,
79 " in particle_bulkprop_names.\n")
80 }
81
82 Tensor4Clip(particle_bulkprop_field, iq, limit_low, limit_high);
83}
84
85/* Workspace method: Doxygen documentation will be auto-generated */
86void vmr_fieldClip(Tensor4& vmr_field,
87 const ArrayOfArrayOfSpeciesTag& abs_species,
88 const String& species,
89 const Numeric& limit_low,
90 const Numeric& limit_high,
91 const Verbosity&) {
92 Index iq = -1;
93 if (species == "ALL") {
94 }
95
96 else {
97 for (Index i = 0; i < abs_species.nelem(); i++) {
98 if (abs_species[i].Species() == SpeciesTag(species).Spec()) {
99 iq = i;
100 break;
101 }
102 }
103 ARTS_USER_ERROR_IF (iq < 0,
104 "Could not find ", species, " in abs_species.\n")
105 }
106
107 Tensor4Clip(vmr_field, iq, limit_low, limit_high);
108}
109
110/* Workspace method: Doxygen documentation will be auto-generated */
111void xClip(Vector& x,
112 const ArrayOfRetrievalQuantity& jacobian_quantities,
113 const Index& ijq,
114 const Numeric& limit_low,
115 const Numeric& limit_high,
116 const Verbosity&) {
117 // Sizes
118 const Index nq = jacobian_quantities.nelem();
119
120 ARTS_USER_ERROR_IF (ijq < -1, "Argument *ijq* must be >= -1.");
121 ARTS_USER_ERROR_IF (ijq >= nq,
122 "Argument *ijq* is too high.\n"
123 "You have selected index: ", ijq, "\n"
124 "but the number of quantities is only: ", nq, "\n"
125 "(Note that zero-based indexing is used)\n")
126
127 // Jacobian indices
128 ArrayOfArrayOfIndex ji;
129 {
130 bool any_affine;
131 jac_ranges_indices(ji, any_affine, jacobian_quantities);
132 }
133
134 Index ifirst = 0, ilast = x.nelem() - 1;
135 if (ijq > -1) {
136 ifirst = ji[ijq][0];
137 ilast = ji[ijq][1];
138 }
139
140 if (!std::isinf(limit_low)) {
141 for (Index i = ifirst; i <= ilast; i++) {
142 if (x[i] < limit_low) x[i] = limit_low;
143 }
144 }
145 if (!std::isinf(limit_high)) {
146 for (Index i = ifirst; i <= ilast; i++) {
147 if (x[i] > limit_high) x[i] = limit_high;
148 }
149 }
150}
151
152/* Workspace method: Doxygen documentation will be auto-generated */
153void xaStandard(Workspace& ws,
154 Vector& xa,
155 const ArrayOfRetrievalQuantity& jacobian_quantities,
156 const Index& atmfields_checked,
157 const Index& atmgeom_checked,
158 const Index& atmosphere_dim,
159 const Vector& p_grid,
160 const Vector& lat_grid,
161 const Vector& lon_grid,
162 const Tensor3& t_field,
163 const Tensor4& vmr_field,
164 const ArrayOfArrayOfSpeciesTag& abs_species,
165 const Index& cloudbox_on,
166 const Index& cloudbox_checked,
167 const Tensor4& particle_bulkprop_field,
168 const ArrayOfString& particle_bulkprop_names,
169 const Tensor3& wind_u_field,
170 const Tensor3& wind_v_field,
171 const Tensor3& wind_w_field,
172 const Tensor3& mag_u_field,
173 const Tensor3& mag_v_field,
174 const Tensor3& mag_w_field,
175 const Tensor3& surface_props_data,
176 const ArrayOfString& surface_props_names,
177 const Agenda& water_p_eq_agenda,
178 const Verbosity&) {
179 // Basics
180 //
181 ARTS_USER_ERROR_IF (atmfields_checked != 1,
182 "The atmospheric fields must be flagged to have "
183 "passed a consistency check (atmfields_checked=1).");
184 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
185 "The atmospheric geometry must be flagged to have "
186 "passed a consistency check (atmgeom_checked=1).");
187 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
188 "The cloudbox must be flagged to have "
189 "passed a consistency check (cloudbox_checked=1).");
190
191 // Jacobian indices
192 ArrayOfArrayOfIndex ji;
193 {
194 bool any_affine;
195 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
196 }
197
198 // Sizes
199 const Index nq = jacobian_quantities.nelem();
200 //
201 xa.resize(ji[nq - 1][1] + 1);
202
203 // Loop retrieval quantities and fill *xa*
204 for (Index q = 0; q < jacobian_quantities.nelem(); q++) {
205 // Index range of this retrieval quantity
206 const Index np = ji[q][1] - ji[q][0] + 1;
207 Range ind(ji[q][0], np);
208
209 // Atmospheric temperatures
210 if (jacobian_quantities[q] == Jacobian::Atm::Temperature) {
211 // Here we need to interpolate *t_field*
212 ArrayOfGridPos gp_p, gp_lat, gp_lon;
213 get_gp_atmgrids_to_rq(gp_p,
214 gp_lat,
215 gp_lon,
216 jacobian_quantities[q],
217 atmosphere_dim,
218 p_grid,
219 lat_grid,
220 lon_grid);
221 Tensor3 t_x;
222 regrid_atmfield_by_gp(t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
223 flat(xa[ind], t_x);
224 }
225
226 // Abs species
227 else if (jacobian_quantities[q].Target().isSpeciesVMR()) {
228 // Index position of species
229 ArrayOfSpeciesTag atag(jacobian_quantities[q].Subtag());
230 const Index isp = chk_contains("abs_species", abs_species, atag);
231
232 if (jacobian_quantities[q].Mode() == "rel") {
233 // This one is simple, just a vector of ones
234 xa[ind] = 1;
235 } else {
236 // For all remaining options we need to interpolate *vmr_field*
237 ArrayOfGridPos gp_p, gp_lat, gp_lon;
238 get_gp_atmgrids_to_rq(gp_p,
239 gp_lat,
240 gp_lon,
241 jacobian_quantities[q],
242 atmosphere_dim,
243 p_grid,
244 lat_grid,
245 lon_grid);
246 Tensor3 vmr_x;
247 regrid_atmfield_by_gp(vmr_x,
248 atmosphere_dim,
249 vmr_field(isp, joker, joker, joker),
250 gp_p,
251 gp_lat,
252 gp_lon);
253
254 if (jacobian_quantities[q].Mode() == "vmr") {
255 flat(xa[ind], vmr_x);
256 } else if (jacobian_quantities[q].Mode() == "nd") {
257 // Here we need to also interpolate *t_field*
258 Tensor3 t_x;
259 regrid_atmfield_by_gp(
260 t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
261 // Calculate number density for species (vmr*nd_tot)
262 Index i = 0;
263 for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
264 for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
265 for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
266 xa[ji[q][0] + i] =
267 vmr_x(i1, i2, i3) *
268 number_density(jacobian_quantities[q].Grids()[0][i1],
269 t_x(i1, i2, i3));
270 i += 1;
271 }
272 }
273 }
274 } else if (jacobian_quantities[q].Mode() == "rh") {
275 // Here we need to also interpolate *t_field*
276 Tensor3 t_x;
277 regrid_atmfield_by_gp(
278 t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
279 Tensor3 water_p_eq;
280 water_p_eq_agendaExecute(ws, water_p_eq, t_x, water_p_eq_agenda);
281 // Calculate relative humidity (vmr*p/p_sat)
282 Index i = 0;
283 for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
284 for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
285 for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
286 xa[ji[q][0] + i] = vmr_x(i1, i2, i3) *
287 jacobian_quantities[q].Grids()[0][i1] /
288 water_p_eq(i1, i2, i3);
289 i += 1;
290 }
291 }
292 }
293 } else if (jacobian_quantities[q].Mode() == "q") {
294 // Calculate specific humidity q, from mixing ratio r and
295 // vapour pressure e, as
296 // q = r(1+r); r = 0.622e/(p-e); e = vmr*p;
297 Index i = 0;
298 for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
299 for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
300 for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
301 const Numeric e =
302 vmr_x(i1, i2, i3) * jacobian_quantities[q].Grids()[0][i1];
303 const Numeric r =
304 0.622 * e / (jacobian_quantities[q].Grids()[0][i1] - e);
305 xa[ji[q][0] + i] = r / (1 + r);
306 i += 1;
307 }
308 }
309 }
310 } else {
311 ARTS_ASSERT(0);
312 }
313 }
314 }
315
316 // Scattering species
317 else if (jacobian_quantities[q] == Jacobian::Special::ScatteringString) {
318 if (cloudbox_on) {
319 ARTS_USER_ERROR_IF (particle_bulkprop_field.empty(),
320 "One jacobian quantity belongs to the "
321 "scattering species category, but *particle_bulkprop_field* "
322 "is empty.");
323 ARTS_USER_ERROR_IF (particle_bulkprop_field.nbooks() !=
324 particle_bulkprop_names.nelem(),
325 "Mismatch in size between "
326 "*particle_bulkprop_field* and *particle_bulkprop_names*.");
327
328 const Index isp = find_first(particle_bulkprop_names,
329 jacobian_quantities[q].SubSubtag());
330 ARTS_USER_ERROR_IF (isp < 0,
331 "Jacobian quantity with index ", q, " covers a "
332 "scattering species, and the field quantity is set to \"",
333 jacobian_quantities[q].SubSubtag(), "\", but this quantity "
334 "could not found in *particle_bulkprop_names*.")
335
336 ArrayOfGridPos gp_p, gp_lat, gp_lon;
338 gp_lat,
339 gp_lon,
340 jacobian_quantities[q],
341 atmosphere_dim,
342 p_grid,
343 lat_grid,
344 lon_grid);
345 Tensor3 pbp_x;
347 atmosphere_dim,
348 particle_bulkprop_field(isp, joker, joker, joker),
349 gp_p,
350 gp_lat,
351 gp_lon);
352 flat(xa[ind], pbp_x);
353 } else {
354 xa[ind] = 0;
355 }
356 }
357
358 // Wind
359 else if (jacobian_quantities[q].Target().isWind()) {
360 ConstTensor3View source_field(wind_u_field);
361 if (jacobian_quantities[q] == Jacobian::Atm::WindV) {
362 source_field = wind_v_field;
363 } else if (jacobian_quantities[q] == Jacobian::Atm::WindW) {
364 source_field = wind_w_field;
365 }
366
367 // Determine grid positions for interpolation from retrieval grids back
368 // to atmospheric grids
369 ArrayOfGridPos gp_p, gp_lat, gp_lon;
371 gp_lat,
372 gp_lon,
373 jacobian_quantities[q],
374 atmosphere_dim,
375 p_grid,
376 lat_grid,
377 lon_grid);
378
379 Tensor3 wind_x;
381 wind_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
382 flat(xa[ind], wind_x);
383 }
384
385 // Magnetism
386 else if (jacobian_quantities[q].Target().isMagnetic()) {
387 if (jacobian_quantities[q] == Jacobian::Atm::MagneticMagnitude) {
388 // Determine grid positions for interpolation from retrieval grids back
389 // to atmospheric grids
390 ArrayOfGridPos gp_p, gp_lat, gp_lon;
392 gp_lat,
393 gp_lon,
394 jacobian_quantities[q],
395 atmosphere_dim,
396 p_grid,
397 lat_grid,
398 lon_grid);
399
400 //all three component's hyoptenuse is the strength
401 Tensor3 mag_u, mag_v, mag_w;
403 mag_u, atmosphere_dim, mag_u_field, gp_p, gp_lat, gp_lon);
405 mag_v, atmosphere_dim, mag_v_field, gp_p, gp_lat, gp_lon);
407 mag_w, atmosphere_dim, mag_w_field, gp_p, gp_lat, gp_lon);
408
409 Tensor3 mag_x(gp_p.nelem(), gp_lat.nelem(), gp_lon.nelem());
410 for (Index i = 0; i < gp_p.nelem(); i++)
411 for (Index j = 0; j < gp_lat.nelem(); j++)
412 for (Index k = 0; k < gp_lon.nelem(); k++)
413 mag_x(i, j, k) = std::hypot(mag_u(i, j, k), mag_v(i, j, k), mag_w(i, j, k));
414 flat(xa[ind], mag_x);
415 } else {
416 ConstTensor3View source_field(mag_u_field);
417 if (jacobian_quantities[q] == Jacobian::Atm::MagneticV) {
418 source_field = mag_v_field;
419 } else if (jacobian_quantities[q] == Jacobian::Atm::MagneticW) {
420 source_field = mag_w_field;
421 }
422
423 // Determine grid positions for interpolation from retrieval grids back
424 // to atmospheric grids
425 ArrayOfGridPos gp_p, gp_lat, gp_lon;
427 gp_lat,
428 gp_lon,
429 jacobian_quantities[q],
430 atmosphere_dim,
431 p_grid,
432 lat_grid,
433 lon_grid);
434
435 Tensor3 mag_x;
437 mag_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
438 flat(xa[ind], mag_x);
439 }
440 }
441
442 // Surface
443 else if (jacobian_quantities[q] == Jacobian::Special::SurfaceString) {
444 surface_props_check(atmosphere_dim,
445 lat_grid,
446 lon_grid,
447 surface_props_data,
448 surface_props_names);
449 ARTS_USER_ERROR_IF (surface_props_data.empty(),
450 "One jacobian quantity belongs to the "
451 "surface category, but *surface_props_data* is empty.");
452
453 const Index isu =
454 find_first(surface_props_names, jacobian_quantities[q].Subtag());
455 ARTS_USER_ERROR_IF (isu < 0,
456 "Jacobian quantity with index ", q, " covers a "
457 "surface property, and the field Subtag is set to \"",
458 jacobian_quantities[q].Subtag(), "\", but this quantity "
459 "could not found in *surface_props_names*.")
460
461 ArrayOfGridPos gp_lat, gp_lon;
463 gp_lon,
464 jacobian_quantities[q],
465 atmosphere_dim,
466 lat_grid,
467 lon_grid);
468 Matrix surf_x;
470 atmosphere_dim,
471 surface_props_data(isu, joker, joker),
472 gp_lat,
473 gp_lon);
474 flat(xa[ind], surf_x);
475 }
476
477 // All variables having zero as a priori
478 // ----------------------------------------------------------------------------
479 else if (jacobian_quantities[q].Target().isPointing() ||
480 jacobian_quantities[q].Target().isFrequency() ||
481 jacobian_quantities[q] == Jacobian::Sensor::Polyfit ||
482 jacobian_quantities[q] == Jacobian::Sensor::Sinefit) {
483 xa[ind] = 0;
484 }
485
486 else {
488 "Found a retrieval quantity that is not yet handled by\n"
489 "internal retrievals: ", jacobian_quantities[q], '\n')
490 }
491 }
492
493 // Apply transformations
494 transform_x(xa, jacobian_quantities);
495}
496
497/* Workspace method: Doxygen documentation will be auto-generated */
499 Tensor4& vmr_field,
500 Tensor3& t_field,
501 Tensor4& particle_bulkprop_field,
502 Tensor3& wind_u_field,
503 Tensor3& wind_v_field,
504 Tensor3& wind_w_field,
505 Tensor3& mag_u_field,
506 Tensor3& mag_v_field,
507 Tensor3& mag_w_field,
508 Tensor3& surface_props_data,
509 const ArrayOfRetrievalQuantity& jacobian_quantities,
510 const Vector& x,
511 const Index& atmfields_checked,
512 const Index& atmgeom_checked,
513 const Index& atmosphere_dim,
514 const Vector& p_grid,
515 const Vector& lat_grid,
516 const Vector& lon_grid,
517 const ArrayOfArrayOfSpeciesTag& abs_species,
518 const Index& cloudbox_on,
519 const Index& cloudbox_checked,
520 const ArrayOfString& particle_bulkprop_names,
521 const ArrayOfString& surface_props_names,
522 const Agenda& water_p_eq_agenda,
523 const Verbosity&) {
524 // Basics
525 //
526 ARTS_USER_ERROR_IF (atmfields_checked != 1,
527 "The atmospheric fields must be flagged to have "
528 "passed a consistency check (atmfields_checked=1).");
529 ARTS_USER_ERROR_IF (atmgeom_checked != 1,
530 "The atmospheric geometry must be flagged to have "
531 "passed a consistency check (atmgeom_checked=1).");
532 ARTS_USER_ERROR_IF (cloudbox_checked != 1,
533 "The cloudbox must be flagged to have "
534 "passed a consistency check (cloudbox_checked=1).");
535
536 // Revert transformation
537 Vector x_t(x);
538 transform_x_back(x_t, jacobian_quantities);
539
540 // Main sizes
541 const Index nq = jacobian_quantities.nelem();
542
543 // Jacobian indices
545 {
546 bool any_affine;
547 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
548 }
549
550 // Check input
551 ARTS_USER_ERROR_IF (x_t.nelem() != ji[nq - 1][1] + 1,
552 "Length of *x* does not match length implied by "
553 "*jacobian_quantities*.");
554
555 // Note that when this method is called, vmr_field and other output variables
556 // have original values, i.e. matching the a priori state.
557
558 // Loop retrieval quantities
559 for (Index q = 0; q < nq; q++) {
560 // Index range of this retrieval quantity
561 const Index np = ji[q][1] - ji[q][0] + 1;
562 Range ind(ji[q][0], np);
563
564 // Atmospheric temperatures
565 // ----------------------------------------------------------------------------
566 if (jacobian_quantities[q] == Jacobian::Atm::Temperature) {
567 // Determine grid positions for interpolation from retrieval grids back
568 // to atmospheric grids
569 ArrayOfGridPos gp_p, gp_lat, gp_lon;
570 Index n_p, n_lat, n_lon;
572 gp_lat,
573 gp_lon,
574 n_p,
575 n_lat,
576 n_lon,
577 jacobian_quantities[q].Grids(),
578 atmosphere_dim,
579 p_grid,
580 lat_grid,
581 lon_grid);
582
583 // Map values in x back to t_field
584 Tensor3 t_x(n_p, n_lat, n_lon);
585 reshape(t_x, x_t[ind]);
587 t_field, atmosphere_dim, t_x, gp_p, gp_lat, gp_lon);
588 }
589
590 // Abs species
591 // ----------------------------------------------------------------------------
592 else if (jacobian_quantities[q].Target().isSpeciesVMR()) {
593 // Index position of species
594 ArrayOfSpeciesTag atag(jacobian_quantities[q].Subtag());
595 const Index isp = chk_contains("abs_species", abs_species, atag);
596
597 // Map part of x to a full atmospheric field
598 Tensor3 x_field(vmr_field.npages(), vmr_field.nrows(), vmr_field.ncols());
599 {
600 ArrayOfGridPos gp_p, gp_lat, gp_lon;
601 Index n_p, n_lat, n_lon;
603 gp_lat,
604 gp_lon,
605 n_p,
606 n_lat,
607 n_lon,
608 jacobian_quantities[q].Grids(),
609 atmosphere_dim,
610 p_grid,
611 lat_grid,
612 lon_grid);
613 //
614 Tensor3 t3_x(n_p, n_lat, n_lon);
615 reshape(t3_x, x_t[ind]);
617 x_field, atmosphere_dim, t3_x, gp_p, gp_lat, gp_lon);
618 }
619 //
620 if (jacobian_quantities[q].Mode() == "rel") {
621 // vmr = vmr0 * x
622 vmr_field(isp, joker, joker, joker) *= x_field;
623 } else if (jacobian_quantities[q].Mode() == "vmr") {
624 // vmr = x
625 vmr_field(isp, joker, joker, joker) = x_field;
626 } else if (jacobian_quantities[q].Mode() == "nd") {
627 // vmr = nd / nd_tot
628 for (Index i3 = 0; i3 < vmr_field.ncols(); i3++) {
629 for (Index i2 = 0; i2 < vmr_field.nrows(); i2++) {
630 for (Index i1 = 0; i1 < vmr_field.npages(); i1++) {
631 vmr_field(isp, i1, i2, i3) =
632 x_field(i1, i2, i3) /
633 number_density(p_grid[i1], t_field(i1, i2, i3));
634 }
635 }
636 }
637 } else if (jacobian_quantities[q].Mode() == "rh") {
638 // vmr = x * p_sat / p
639 Tensor3 water_p_eq;
640 water_p_eq_agendaExecute(ws, water_p_eq, t_field, water_p_eq_agenda);
641 for (Index i3 = 0; i3 < vmr_field.ncols(); i3++) {
642 for (Index i2 = 0; i2 < vmr_field.nrows(); i2++) {
643 for (Index i1 = 0; i1 < vmr_field.npages(); i1++) {
644 vmr_field(isp, i1, i2, i3) =
645 x_field(i1, i2, i3) * water_p_eq(i1, i2, i3) / p_grid[i1];
646 }
647 }
648 }
649 } else if (jacobian_quantities[q].Mode() == "q") {
650 // We have that specific humidity q, mixing ratio r and
651 // vapour pressure e, are related as
652 // q = r(1+r); r = 0.622e/(p-e); e = vmr*p;
653 // That is: vmr=e/p; e = rp/(0.622+r); r = q/(1-q)
654 for (Index i3 = 0; i3 < vmr_field.ncols(); i3++) {
655 for (Index i2 = 0; i2 < vmr_field.nrows(); i2++) {
656 for (Index i1 = 0; i1 < vmr_field.npages(); i1++) {
657 const Numeric r = x_field(i1, i2, i3) / (1 - x_field(i1, i2, i3));
658 const Numeric e = r * p_grid[i1] / (0.622 + r);
659 vmr_field(isp, i1, i2, i3) = e / p_grid[i1];
660 }
661 }
662 }
663 } else {
664 ARTS_ASSERT(0);
665 }
666 }
667
668 // Scattering species
669 // ----------------------------------------------------------------------------
670 else if (jacobian_quantities[q] == Jacobian::Special::ScatteringString) {
671 // If no cloudbox, we assume that there is nothing to do
672 if (cloudbox_on) {
673 ARTS_USER_ERROR_IF (particle_bulkprop_field.empty(),
674 "One jacobian quantity belongs to the "
675 "scattering species category, but *particle_bulkprop_field* "
676 "is empty.");
677 ARTS_USER_ERROR_IF (particle_bulkprop_field.nbooks() !=
678 particle_bulkprop_names.nelem(),
679 "Mismatch in size between "
680 "*particle_bulkprop_field* and *particle_bulkprop_field*.");
681
682 const Index isp = find_first(particle_bulkprop_names,
683 jacobian_quantities[q].SubSubtag());
684 ARTS_USER_ERROR_IF (isp < 0,
685 "Jacobian quantity with index ", q, " covers a "
686 "scattering species, and the field quantity is set to \"",
687 jacobian_quantities[q].SubSubtag(), "\", but this quantity "
688 "could not found in *particle_bulkprop_names*.")
689
690 // Determine grid positions for interpolation from retrieval grids back
691 // to atmospheric grids
692 ArrayOfGridPos gp_p, gp_lat, gp_lon;
693 Index n_p, n_lat, n_lon;
695 gp_lat,
696 gp_lon,
697 n_p,
698 n_lat,
699 n_lon,
700 jacobian_quantities[q].Grids(),
701 atmosphere_dim,
702 p_grid,
703 lat_grid,
704 lon_grid);
705 // Map x to particle_bulkprop_field
706 Tensor3 pbfield_x(n_p, n_lat, n_lon);
707 reshape(pbfield_x, x_t[ind]);
708 Tensor3 pbfield;
710 pbfield, atmosphere_dim, pbfield_x, gp_p, gp_lat, gp_lon);
711 particle_bulkprop_field(isp, joker, joker, joker) = pbfield;
712 }
713 }
714
715 // Wind
716 // ----------------------------------------------------------------------------
717 else if (jacobian_quantities[q].Target().isWind()) {
718 // Determine grid positions for interpolation from retrieval grids back
719 // to atmospheric grids
720 ArrayOfGridPos gp_p, gp_lat, gp_lon;
721 Index n_p, n_lat, n_lon;
723 gp_lat,
724 gp_lon,
725 n_p,
726 n_lat,
727 n_lon,
728 jacobian_quantities[q].Grids(),
729 atmosphere_dim,
730 p_grid,
731 lat_grid,
732 lon_grid);
733
734 // TODO Could be done without copying.
735 Tensor3 wind_x(n_p, n_lat, n_lon);
736 reshape(wind_x, x_t[ind]);
737
738 Tensor3View target_field(wind_u_field);
739
740 Tensor3 wind_field(
741 target_field.npages(), target_field.nrows(), target_field.ncols());
743 wind_field, atmosphere_dim, wind_x, gp_p, gp_lat, gp_lon);
744
745 if (jacobian_quantities[q] == Jacobian::Atm::WindU) {
746 wind_u_field = wind_field;
747 } else if (jacobian_quantities[q] == Jacobian::Atm::WindV) {
748 wind_v_field = wind_field;
749 } else if (jacobian_quantities[q] == Jacobian::Atm::WindW) {
750 wind_w_field = wind_field;
751 }
752 }
753
754 // Magnetism
755 // ----------------------------------------------------------------------------
756 else if (jacobian_quantities[q].Target().isMagnetic()) {
757 // Determine grid positions for interpolation from retrieval grids back
758 // to atmospheric grids
759 ArrayOfGridPos gp_p, gp_lat, gp_lon;
760 Index n_p, n_lat, n_lon;
762 gp_lat,
763 gp_lon,
764 n_p,
765 n_lat,
766 n_lon,
767 jacobian_quantities[q].Grids(),
768 atmosphere_dim,
769 p_grid,
770 lat_grid,
771 lon_grid);
772
773 // TODO Could be done without copying.
774 Tensor3 mag_x(n_p, n_lat, n_lon);
775 reshape(mag_x, x_t[ind]);
776
777 Tensor3View target_field(mag_u_field);
778
779 Tensor3 mag_field(
780 target_field.npages(), target_field.nrows(), target_field.ncols());
782 mag_field, atmosphere_dim, mag_x, gp_p, gp_lat, gp_lon);
783 if (jacobian_quantities[q] == Jacobian::Atm::MagneticU) {
784 mag_u_field = mag_field;
785 } else if (jacobian_quantities[q] == Jacobian::Atm::MagneticV) {
786 mag_v_field = mag_field;
787 } else if (jacobian_quantities[q] == Jacobian::Atm::MagneticW) {
788 mag_w_field = mag_field;
789 } else if (jacobian_quantities[q] == Jacobian::Atm::MagneticMagnitude) {
790 for (Index i = 0; i < n_p; i++) {
791 for (Index j = 0; j < n_lat; j++) {
792 for (Index k = 0; k < n_lon; k++) {
793 const Numeric scale = mag_x(i, j, k) /
794 std::hypot(mag_u_field(i, j, k), mag_v_field(i, j, k), mag_w_field(i, j, k));
795 mag_u_field(i, j, k) *= scale;
796 mag_v_field(i, j, k) *= scale;
797 mag_w_field(i, j, k) *= scale;
798 }
799 }
800 }
801 } else {
802 ARTS_USER_ERROR ( "Unsupported magnetism type");
803 }
804 }
805
806 // Surface
807 // ----------------------------------------------------------------------------
808 else if (jacobian_quantities[q] == Jacobian::Special::SurfaceString) {
809 surface_props_check(atmosphere_dim,
810 lat_grid,
811 lon_grid,
812 surface_props_data,
813 surface_props_names);
814 ARTS_USER_ERROR_IF (surface_props_data.empty(),
815 "One jacobian quantity belongs to the "
816 "surface category, but *surface_props_data* is empty.");
817
818 const Index isu =
819 find_first(surface_props_names, jacobian_quantities[q].Subtag());
820 ARTS_USER_ERROR_IF (isu < 0,
821 "Jacobian quantity with index ", q, " covers a "
822 "surface property, and the field Subtag is set to \"",
823 jacobian_quantities[q].Subtag(), "\", but this quantity "
824 "could not found in *surface_props_names*.")
825
826 // Determine grid positions for interpolation from retrieval grids back
827 // to atmospheric grids
828 ArrayOfGridPos gp_lat, gp_lon;
829 Index n_lat, n_lon;
831 gp_lon,
832 n_lat,
833 n_lon,
834 jacobian_quantities[q].Grids(),
835 atmosphere_dim,
836 lat_grid,
837 lon_grid);
838 // Map values in x back to surface_props_data
839 Matrix surf_x(n_lat, n_lon);
840 reshape(surf_x, x_t[ind]);
841 Matrix surf;
842 regrid_atmsurf_by_gp_oem(surf, atmosphere_dim, surf_x, gp_lat, gp_lon);
843 surface_props_data(isu, joker, joker) = surf;
844 }
845 }
846}
847
848/* Workspace method: Doxygen documentation will be auto-generated */
850 Matrix& sensor_los,
851 Vector& f_backend,
852 Vector& y_baseline,
853 Sparse& sensor_response,
854 Vector& sensor_response_f,
855 ArrayOfIndex& sensor_response_pol,
856 Matrix& sensor_response_dlos,
857 Vector& sensor_response_f_grid,
858 ArrayOfIndex& sensor_response_pol_grid,
859 Matrix& sensor_response_dlos_grid,
860 Matrix& mblock_dlos_grid,
861 const ArrayOfRetrievalQuantity& jacobian_quantities,
862 const Vector& x,
863 const Agenda& sensor_response_agenda,
864 const Index& sensor_checked,
865 const ArrayOfTime& sensor_time,
866 const Verbosity&) {
867 // Basics
868 //
869 ARTS_USER_ERROR_IF (sensor_checked != 1,
870 "The sensor response must be flagged to have "
871 "passed a consistency check (sensor_checked=1).");
872
873 // Revert transformation
874 Vector x_t(x);
875 transform_x_back(x_t, jacobian_quantities);
876
877 // Main sizes
878 const Index nq = jacobian_quantities.nelem();
879
880 // Jacobian indices
882 {
883 bool any_affine;
884 jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
885 }
886
887 // Check input
888 ARTS_USER_ERROR_IF (x_t.nelem() != ji[nq - 1][1] + 1,
889 "Length of *x* does not match length implied by "
890 "*jacobian_quantities*.");
891
892 // Flag indicating that y_baseline is not set
893 bool yb_set = false;
894
895 // Shall sensor responses be calculed?
896 bool do_sensor = false;
897
898 // Loop retrieval quantities
899 for (Index q = 0; q < nq; q++) {
900 // Index range of this retrieval quantity
901 const Index np = ji[q][1] - ji[q][0] + 1;
902
903 // Pointing off-set
904 // ----------------------------------------------------------------------------
905 if (jacobian_quantities[q].Target().isPointing()) {
906 // Handle pointing "jitter" seperately
907 if (jacobian_quantities[q].Grids()[0][0] == -1) {
908 ARTS_USER_ERROR_IF (sensor_los.nrows() != np,
909 "Mismatch between pointing jacobian and *sensor_los* found.");
910 // Simply add retrieved off-set(s) to za column of *sensor_los*
911 for (Index i = 0; i < np; i++) {
912 sensor_los(i, 0) += x_t[ji[q][0] + i];
913 }
914 }
915 // Polynomial representation
916 else {
917 ARTS_USER_ERROR_IF (sensor_los.nrows() != sensor_time.nelem(),
918 "Sizes of *sensor_los* and *sensor_time* do not match.");
919 Vector w;
920 for (Index c = 0; c < np; c++) {
921 polynomial_basis_func(w, time_vector(sensor_time), c);
922 for (Index i = 0; i < w.nelem(); i++) {
923 sensor_los(i, 0) += w[i] * x_t[ji[q][0] + c];
924 }
925 }
926 }
927 }
928
929 // Frequency shift or stretch
930 // ----------------------------------------------------------------------------
931 else if (jacobian_quantities[q].Target().isFrequency()) {
932 if (jacobian_quantities[q] == Jacobian::Sensor::FrequencyShift) {
933 ARTS_ASSERT(np == 1);
934 if (x_t[ji[q][0]] != 0) {
935 do_sensor = true;
936 f_backend += x_t[ji[q][0]];
937 }
938 } else if (jacobian_quantities[q] == Jacobian::Sensor::FrequencyStretch) {
939 ARTS_ASSERT(np == 1);
940 if (x_t[ji[q][0]] != 0) {
941 do_sensor = true;
942 Vector w;
943 polynomial_basis_func(w, f_backend, 1);
944 for (Index i = 0; i < w.nelem(); i++) {
945 f_backend[i] += w[i] * x_t[ji[q][0]];
946 }
947 }
948 } else {
949 ARTS_ASSERT(0);
950 }
951 }
952
953 // Baseline fit: polynomial or sinusoidal
954 // ----------------------------------------------------------------------------
955 else if (jacobian_quantities[q] == Jacobian::Sensor::Polyfit ||
956 jacobian_quantities[q] == Jacobian::Sensor::Sinefit) {
957 if (!yb_set) {
958 yb_set = true;
959 Index y_size = sensor_los.nrows() * sensor_response_f_grid.nelem() *
960 sensor_response_pol_grid.nelem() *
961 sensor_response_dlos_grid.nrows();
962 y_baseline.resize(y_size);
963 y_baseline = 0;
964 }
965
966 for (Index mb = 0; mb < sensor_los.nrows(); ++mb) {
967 calcBaselineFit(y_baseline,
968 x_t,
969 mb,
970 sensor_response,
971 sensor_response_pol_grid,
972 sensor_response_f_grid,
973 sensor_response_dlos_grid,
974 jacobian_quantities[q],
975 q,
976 ji);
977 }
978 }
979 }
980
981 // *y_baseline* not yet set?
982 if (!yb_set) {
983 y_baseline.resize(1);
984 y_baseline[0] = 0;
985 }
986
987 // Recalculate sensor_response?
988 if (do_sensor) {
990 sensor_response,
991 sensor_response_f,
992 sensor_response_f_grid,
993 sensor_response_pol,
994 sensor_response_pol_grid,
995 sensor_response_dlos,
996 sensor_response_dlos_grid,
997 mblock_dlos_grid,
998 f_backend,
999 sensor_response_agenda);
1000 }
1001}
1002
1003/* Workspace method: Doxygen documentation will be auto-generated */
1005 ARTS_USER_ERROR ( "Retrievals of spectroscopic variables not yet handled.");
1006}
1007
1008
1009#ifdef OEM_SUPPORT
1010/* Workspace method: Doxygen documentation will be auto-generated */
1011void OEM(Workspace& ws,
1012 Vector& x,
1013 Vector& yf,
1014 Matrix& jacobian,
1015 Matrix& dxdy,
1016 Vector& oem_diagnostics,
1017 Vector& lm_ga_history,
1018 ArrayOfString& errors,
1019 const Vector& xa,
1020 const CovarianceMatrix& covmat_sx,
1021 const Vector& y,
1022 const CovarianceMatrix& covmat_se,
1023 const ArrayOfRetrievalQuantity& jacobian_quantities,
1024 const Agenda& inversion_iterate_agenda,
1025 const String& method,
1026 const Numeric& max_start_cost,
1027 const Vector& x_norm,
1028 const Index& max_iter,
1029 const Numeric& stop_dx,
1030 const Vector& lm_ga_settings,
1031 const Index& clear_matrices,
1032 const Index& display_progress,
1033 const Verbosity&) {
1034 // Main sizes
1035 const Index n = covmat_sx.nrows();
1036 const Index m = y.nelem();
1037
1038 // Checks
1039 covmat_sx.compute_inverse();
1040 covmat_se.compute_inverse();
1041
1042 OEM_checks(ws,
1043 x,
1044 yf,
1045 jacobian,
1046 inversion_iterate_agenda,
1047 xa,
1048 covmat_sx,
1049 y,
1050 covmat_se,
1051 jacobian_quantities,
1052 method,
1053 x_norm,
1054 max_iter,
1055 stop_dx,
1056 lm_ga_settings,
1057 clear_matrices,
1058 display_progress);
1059
1060 // Size diagnostic output and init with NaNs
1061 oem_diagnostics.resize(5);
1062 oem_diagnostics = NAN;
1063 //
1064 if (method == "ml" || method == "lm" || method == "ml_cg" ||
1065 method == "lm_cg") {
1066 lm_ga_history.resize(max_iter + 1);
1067 lm_ga_history = NAN;
1068 } else {
1069 lm_ga_history.resize(0);
1070 }
1071
1072 // Check for start vector and precomputed yf, jacobian
1073 if (x.nelem() != n) {
1074 x = xa;
1075 yf.resize(0);
1076 jacobian.resize(0, 0);
1077 }
1078
1079 // If no precomputed value given, we compute yf and jacobian to
1080 // compute initial cost (and use in the first OEM iteration).
1081 if (yf.nelem() == 0) {
1083 ws, yf, jacobian, xa, 1, 0, inversion_iterate_agenda);
1084 }
1085
1086 ARTS_USER_ERROR_IF (yf.nelem() not_eq y.nelem(),
1087 "Mismatch between simulated y and input y.\n"
1088 "Input y is size ", y.nelem(), " but simulated y is ",
1089 yf.nelem(), "\n"
1090 "Use your frequency grid vector and your sensor response matrix to match simulations with measurements.\n")
1091
1092 // TODO: Get this from invlib log.
1093 // Start value of cost function
1094 Numeric cost_start = NAN;
1095 if (method == "ml" || method == "lm" || display_progress ||
1096 max_start_cost > 0) {
1097 Vector dy = y;
1098 dy -= yf;
1099 Vector sdy = y;
1100 mult_inv(sdy, covmat_se, dy);
1101 Vector dx = x;
1102 dx -= xa;
1103 Vector sdx = x;
1104 mult_inv(sdx, covmat_sx, dx);
1105 cost_start = dx * sdx + dy * sdy;
1106 cost_start /= static_cast<Numeric>(m);
1107 }
1108 oem_diagnostics[1] = cost_start;
1109
1110 // Handle cases with too large start cost
1111 if (max_start_cost > 0 && cost_start > max_start_cost) {
1112 // Flag no inversion in oem_diagnostics, and let x to be undefined
1113 oem_diagnostics[0] = 99;
1114 //
1115 if (display_progress) {
1116 cout << "\n No OEM inversion, too high start cost:\n"
1117 << " Set limit : " << max_start_cost << endl
1118 << " Found value : " << cost_start << endl
1119 << endl;
1120 }
1121 }
1122 // Otherwise do inversion
1123 else {
1124 bool apply_norm = false;
1125 oem::Matrix T{};
1126 if (x_norm.nelem() == n) {
1127 T.resize(n, n);
1128 T *= 0.0;
1129 T.diagonal() = x_norm;
1130 for (Index i = 0; i < n; i++) {
1131 T(i, i) = x_norm[i];
1132 }
1133 apply_norm = true;
1134 }
1135
1136 oem::CovarianceMatrix Se(covmat_se), Sa(covmat_sx);
1137 oem::Vector xa_oem(xa), y_oem(y), x_oem(x);
1138 oem::AgendaWrapper aw(&ws,
1139 (unsigned int)m,
1140 (unsigned int)n,
1141 jacobian,
1142 yf,
1143 &inversion_iterate_agenda);
1144 oem::OEM_STANDARD<oem::AgendaWrapper> oem(aw, xa_oem, Sa, Se);
1145 oem::OEM_MFORM<oem::AgendaWrapper> oem_m(aw, xa_oem, Sa, Se);
1146 int oem_verbosity = static_cast<int>(display_progress);
1147
1148 int return_code = 0;
1149
1150 try {
1151 if (method == "li") {
1152 oem::Std s(T, apply_norm);
1153 oem::GN gn(stop_dx, 1, s); // Linear case, only one step.
1154 return_code = oem.compute<oem::GN, oem::ArtsLog>(
1155 x_oem, y_oem, gn, oem_verbosity, lm_ga_history, true);
1156 oem_diagnostics[0] = static_cast<Index>(return_code);
1157 } else if (method == "li_m") {
1158 oem::Std s(T, apply_norm);
1159 oem::GN gn(stop_dx, 1, s); // Linear case, only one step.
1160 return_code = oem_m.compute<oem::GN, oem::ArtsLog>(
1161 x_oem, y_oem, gn, oem_verbosity, lm_ga_history, true);
1162 oem_diagnostics[0] = static_cast<Index>(return_code);
1163 } else if (method == "li_cg") {
1164 oem::CG cg(T, apply_norm, 1e-10, 0);
1165 oem::GN_CG gn(stop_dx, 1, cg); // Linear case, only one step.
1166 return_code = oem.compute<oem::GN_CG, oem::ArtsLog>(
1167 x_oem, y_oem, gn, oem_verbosity, lm_ga_history, true);
1168 oem_diagnostics[0] = static_cast<Index>(return_code);
1169 } else if (method == "li_cg_m") {
1170 oem::CG cg(T, apply_norm, 1e-10, 0);
1171 oem::GN_CG gn(stop_dx, 1, cg); // Linear case, only one step.
1172 return_code = oem_m.compute<oem::GN_CG, oem::ArtsLog>(
1173 x_oem, y_oem, gn, oem_verbosity, lm_ga_history, true);
1174 oem_diagnostics[0] = static_cast<Index>(return_code);
1175 } else if (method == "gn") {
1176 oem::Std s(T, apply_norm);
1177 oem::GN gn(stop_dx, (unsigned int)max_iter, s);
1178 return_code = oem.compute<oem::GN, oem::ArtsLog>(
1179 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1180 oem_diagnostics[0] = static_cast<Index>(return_code);
1181 } else if (method == "gn_m") {
1182 oem::Std s(T, apply_norm);
1183 oem::GN gn(stop_dx, (unsigned int)max_iter, s);
1184 return_code = oem_m.compute<oem::GN, oem::ArtsLog>(
1185 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1186 oem_diagnostics[0] = static_cast<Index>(return_code);
1187 } else if (method == "gn_cg") {
1188 oem::CG cg(T, apply_norm, 1e-10, 0);
1189 oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1190 return_code = oem.compute<oem::GN_CG, oem::ArtsLog>(
1191 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1192 oem_diagnostics[0] = static_cast<Index>(return_code);
1193 } else if (method == "gn_cg_m") {
1194 oem::CG cg(T, apply_norm, 1e-10, 0);
1195 oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1196 return_code = oem_m.compute<oem::GN_CG, oem::ArtsLog>(
1197 x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1198 oem_diagnostics[0] = static_cast<Index>(return_code);
1199 } else if ((method == "lm") || (method == "ml")) {
1200 oem::Std s(T, apply_norm);
1201 Sparse diagonal = Sparse::diagonal(covmat_sx.inverse_diagonal());
1202 CovarianceMatrix SaDiag{};
1203 SaDiag.add_correlation_inverse(Block(Range(0, n),
1204 Range(0, n),
1205 std::make_pair(0, 0),
1206 make_shared<Sparse>(diagonal)));
1208 oem::LM lm(SaInvLM, s);
1209
1210 lm.set_tolerance(stop_dx);
1211 lm.set_maximum_iterations((unsigned int)max_iter);
1212 lm.set_lambda(lm_ga_settings[0]);
1213 lm.set_lambda_decrease(lm_ga_settings[1]);
1214 lm.set_lambda_increase(lm_ga_settings[2]);
1215 lm.set_lambda_maximum(lm_ga_settings[3]);
1216 lm.set_lambda_threshold(lm_ga_settings[4]);
1217 lm.set_lambda_constraint(lm_ga_settings[5]);
1218
1219 return_code = oem.compute<oem::LM&, oem::ArtsLog>(
1220 x_oem, y_oem, lm, oem_verbosity, lm_ga_history);
1221 oem_diagnostics[0] = static_cast<Index>(return_code);
1222 if (lm.get_lambda() > lm.get_lambda_maximum()) {
1223 oem_diagnostics[0] = 2;
1224 }
1225 } else if ((method == "lm_cg") || (method == "ml_cg")) {
1226 oem::CG cg(T, apply_norm, 1e-10, 0);
1227
1228 Sparse diagonal = Sparse::diagonal(covmat_sx.inverse_diagonal());
1229 CovarianceMatrix SaDiag{};
1230 SaDiag.add_correlation_inverse(Block(Range(0, n),
1231 Range(0, n),
1232 std::make_pair(0, 0),
1233 make_shared<Sparse>(diagonal)));
1234 oem::LM_CG lm(SaDiag, cg);
1235
1236 lm.set_maximum_iterations((unsigned int)max_iter);
1237 lm.set_lambda(lm_ga_settings[0]);
1238 lm.set_lambda_decrease(lm_ga_settings[1]);
1239 lm.set_lambda_increase(lm_ga_settings[2]);
1240 lm.set_lambda_threshold(lm_ga_settings[3]);
1241 lm.set_lambda_maximum(lm_ga_settings[4]);
1242
1243 return_code = oem.compute<oem::LM_CG&, oem::ArtsLog>(
1244 x_oem, y_oem, lm, oem_verbosity, lm_ga_history);
1245 oem_diagnostics[0] = static_cast<Index>(return_code);
1246 if (lm.get_lambda() > lm.get_lambda_maximum()) {
1247 oem_diagnostics[0] = 2;
1248 }
1249 }
1250
1251 oem_diagnostics[2] = oem.cost / static_cast<Numeric>(m);
1252 oem_diagnostics[3] = oem.cost_y / static_cast<Numeric>(m);
1253 oem_diagnostics[4] = static_cast<Numeric>(oem.iterations);
1254 } catch (const std::exception& e) {
1255 oem_diagnostics[0] = 9;
1256 oem_diagnostics[2] = oem.cost;
1257 oem_diagnostics[3] = oem.cost_y;
1258 oem_diagnostics[4] = static_cast<Numeric>(oem.iterations);
1259 x_oem *= NAN;
1260 std::vector<std::string> sv = oem::handle_nested_exception(e);
1261 for (auto& s : sv) {
1262 std::stringstream ss{s};
1263 std::string t{};
1264 while (std::getline(ss, t)) {
1265 errors.push_back(t.c_str());
1266 }
1267 }
1268 } catch (...) {
1269 throw;
1270 }
1271
1272 x = x_oem;
1273 yf = aw.get_measurement_vector();
1274
1275 // Shall empty jacobian and dxdy be returned?
1276 if (clear_matrices) {
1277 jacobian.resize(0, 0);
1278 dxdy.resize(0, 0);
1279 } else if (oem_diagnostics[0] <= 2) {
1280 dxdy.resize(n, m);
1281 Matrix tmp1(n, m), tmp2(n, n), tmp3(n, n);
1282 mult_inv(tmp1, transpose(jacobian), covmat_se);
1283 mult(tmp2, tmp1, jacobian);
1284 add_inv(tmp2, covmat_sx);
1285 inv(tmp3, tmp2);
1286 mult(dxdy, tmp3, tmp1);
1287 }
1288 }
1289}
1290
1291/* Workspace method: Doxygen documentation will be auto-generated */
1292void covmat_soCalc(Matrix& covmat_so,
1293 const Matrix& dxdy,
1294 const CovarianceMatrix& covmat_se,
1295 const Verbosity& /*v*/) {
1296 Index n(dxdy.nrows()), m(dxdy.ncols());
1297 Matrix tmp1(m, n);
1298
1299 ARTS_USER_ERROR_IF ((m == 0) || (n == 0),
1300 "The gain matrix *dxdy* is required to compute the observation error covariance matrix.");
1301 ARTS_USER_ERROR_IF ((covmat_se.nrows() != m) || (covmat_se.ncols() != m),
1302 "The covariance matrix covmat_se has invalid dimensions.");
1303
1304 covmat_so.resize(n, n);
1305 mult(tmp1, covmat_se, transpose(dxdy));
1306 mult(covmat_so, dxdy, tmp1);
1307}
1308
1309/* Workspace method: Doxygen documentation will be auto-generated */
1310void covmat_ssCalc(Matrix& covmat_ss,
1311 const Matrix& avk,
1312 const CovarianceMatrix& covmat_sx,
1313 const Verbosity& /*v*/) {
1314 Index n(avk.ncols());
1315 Matrix tmp1(n, n), tmp2(n, n);
1316
1317 ARTS_USER_ERROR_IF (n == 0,
1318 "The averaging kernel matrix *dxdy* is required to compute the smoothing error covariance matrix.");
1319 ARTS_USER_ERROR_IF ((covmat_sx.nrows() != n) || (covmat_sx.ncols() != n),
1320 "The covariance matrix *covmat_sx* invalid dimensions.");
1321
1322 covmat_ss.resize(n, n);
1323
1324 // Sign doesn't matter since we're dealing with a quadratic form.
1325 id_mat(tmp1);
1326 tmp1 -= avk;
1327
1328 mult(tmp2, covmat_sx, transpose(tmp1));
1329 mult(covmat_ss, tmp1, tmp2);
1330}
1331
1332/* Workspace method: Doxygen documentation will be auto-generated */
1334 const CovarianceMatrix& Sc,
1335 const Verbosity& /*v*/) {
1336 S = Matrix(Sc);
1337}
1338
1339/* Workspace method: Doxygen documentation will be auto-generated */
1340void avkCalc(Matrix& avk,
1341 const Matrix& dxdy,
1342 const Matrix& jacobian,
1343 const Verbosity& /*v*/) {
1344 Index m(jacobian.nrows()), n(jacobian.ncols());
1345 ARTS_USER_ERROR_IF ((m == 0) || (n == 0),
1346 "The Jacobian matrix is empty.");
1347 ARTS_USER_ERROR_IF ((dxdy.nrows() != n) || (dxdy.ncols() != m),
1348 "Matrices have inconsistent sizes.\n"
1349 " Size of gain matrix: ", dxdy.nrows(), " x ", dxdy.ncols(),
1350 "\n"
1351 " Size of Jacobian: ", jacobian.nrows(), " x ",
1352 jacobian.ncols(), "\n")
1353
1354 avk.resize(n, n);
1355 mult(avk, dxdy, jacobian);
1356}
1357
1358#else
1359
1360void covmat_soCalc(Matrix& /* covmat_so */,
1361 const Matrix& /* dxdy */,
1362 const CovarianceMatrix& /* covmat_se*/,
1363 const Verbosity& /*v*/) {
1365 "WSM is not available because ARTS was compiled without "
1366 "OEM support.");
1367}
1368
1369void covmat_ssCalc(Matrix& /*covmat_ss*/,
1370 const Matrix& /*avk*/,
1371 const CovarianceMatrix& /*covmat_sx*/,
1372 const Verbosity& /*v*/) {
1374 "WSM is not available because ARTS was compiled without "
1375 "OEM support.");
1376}
1377
1378void avkCalc(Matrix& /* avk */,
1379 const Matrix& /* dxdy */,
1380 const Matrix& /* jacobian */,
1381 const Verbosity& /*v*/) {
1383 "WSM is not available because ARTS was compiled without "
1384 "OEM support.");
1385}
1386
1388 Vector&,
1389 Vector&,
1390 Matrix&,
1391 Matrix&,
1392 Vector&,
1393 Vector&,
1395 const Vector&,
1396 const CovarianceMatrix&,
1397 const Vector&,
1398 const CovarianceMatrix&,
1399 const Index&,
1401 const ArrayOfArrayOfIndex&,
1402 const Agenda&,
1403 const String&,
1404 const Numeric&,
1405 const Vector&,
1406 const Index&,
1407 const Numeric&,
1408 const Vector&,
1409 const Index&,
1410 const Index&,
1411 const Verbosity&) {
1413 "WSM is not available because ARTS was compiled without "
1414 "OEM support.");
1415}
1416#endif
1417
1418#if defined(OEM_SUPPORT) && 0
1419
1420#include "agenda_wrapper_mpi.h"
1421#include "oem_mpi.h"
1422
1423//
1424// Performs manipulations of workspace variables necessary for distributed
1425// retrievals with MPI:
1426//
1427// - Splits up sensor positions evenly over processes
1428// - Splits up inverse covariance matrices.
1429//
1430void MPI_Initialize(Matrix& sensor_los,
1431 Matrix& sensor_pos,
1432 Vector& sensor_time) {
1433 int initialized;
1434
1435 MPI_Initialized(&initialized);
1436 if (!initialized) {
1437 MPI_Init(nullptr, nullptr);
1438 }
1439
1440 int rank, nprocs;
1441
1442 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1443 MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
1444
1445 int nmblock = (int)sensor_pos.nrows();
1446 int mblock_range = nmblock / nprocs;
1447 int mblock_start = mblock_range * rank;
1448 int remainder = nmblock % std::max(mblock_range, nprocs);
1449
1450 //
1451 // Split up sensor positions.
1452 //
1453
1454 if (rank < remainder) {
1455 mblock_range += 1;
1456 mblock_start += rank;
1457 } else {
1458 mblock_start += remainder;
1459 }
1460
1461 if (nmblock > 0) {
1462 Range range = Range(mblock_start, mblock_range);
1463
1464 Matrix tmp_m = sensor_los(range, joker);
1465 sensor_los = tmp_m;
1466
1467 tmp_m = sensor_pos(range, joker);
1468 sensor_pos = tmp_m;
1469
1470 Vector tmp_v = sensor_time[range];
1471 sensor_time = tmp_v;
1472 } else {
1473 sensor_los.resize(0, 0);
1474 sensor_pos.resize(0, 0);
1475 sensor_time.resize(0);
1476 }
1477}
1478
1479void OEM_MPI(Workspace& ws,
1480 Vector& x,
1481 Vector& yf,
1482 Matrix& jacobian,
1483 Matrix& dxdy,
1484 Vector& oem_diagnostics,
1485 Vector& lm_ga_history,
1486 Matrix& sensor_los,
1487 Matrix& sensor_pos,
1488 Vector& sensor_time,
1489 CovarianceMatrix& covmat_sx,
1490 CovarianceMatrix& covmat_se,
1491 const Vector& xa,
1492 const Vector& y,
1493 const ArrayOfRetrievalQuantity& jacobian_quantities,
1494 const Agenda& inversion_iterate_agenda,
1495 const String& method,
1496 const Numeric& max_start_cost,
1497 const Vector& x_norm,
1498 const Index& max_iter,
1499 const Numeric& stop_dx,
1500 const Vector& lm_ga_settings,
1501 const Index& clear_matrices,
1502 const Index& display_progress,
1503 const Verbosity& /*v*/) {
1504 // Main sizes
1505 const Index n = covmat_sx.nrows();
1506 const Index m = y.nelem();
1507
1508 // Check WSVs
1509 OEM_checks(ws,
1510 x,
1511 yf,
1512 jacobian,
1513 inversion_iterate_agenda,
1514 xa,
1515 covmat_sx,
1516 covmat_se,
1517 jacobian_quantities,
1518 method,
1519 x_norm,
1520 max_iter,
1521 stop_dx,
1522 lm_ga_settings,
1523 clear_matrices,
1524 display_progress);
1525
1526 // Calculate spectrum and Jacobian for a priori state
1527 // Jacobian is also input to the agenda, and to flag this is this first
1528 // call, this WSV must be set to be empty.
1529 jacobian.resize(0, 0);
1530
1531 // Initialize MPI environment.
1532 MPI_Initialize(sensor_los, sensor_pos, sensor_time);
1533
1534 // Setup distributed matrices.
1535 MPICovarianceMatrix SeInvMPI(covmat_se);
1536 MPICovarianceMatrix SaInvMPI(covmat_sx);
1537
1538 // Create temporary MPI vector from local results and use conversion to
1539 // standard vector to broadcast results to all processes.
1540 oem::Vector tmp;
1542 ws, tmp, jacobian, xa, 1, inversion_iterate_agenda);
1543 yf = MPIVector(tmp);
1544
1545 // Size diagnostic output and init with NaNs
1546 oem_diagnostics.resize(5);
1547 oem_diagnostics = NAN;
1548 //
1549 if (method == "ml" || method == "lm") {
1550 lm_ga_history.resize(max_iter);
1551 lm_ga_history = NAN;
1552 } else {
1553 lm_ga_history.resize(0);
1554 }
1555
1556 // Start value of cost function. Covariance matrices are already distributed
1557 // over processes, so we need to use invlib matrix algebra.
1558 Numeric cost_start = NAN;
1559 if (method == "ml" || method == "lm" || display_progress ||
1560 max_start_cost > 0) {
1561 oem::Vector dy = y;
1562 dy -= yf;
1563 cost_start = dot(dy, SeInvMPI * dy);
1564 }
1565 oem_diagnostics[1] = cost_start;
1566
1567 // Handle cases with too large start cost
1568 if (max_start_cost > 0 && cost_start > max_start_cost) {
1569 // Flag no inversion in oem_diagnostics, and let x to be undefined
1570 oem_diagnostics[0] = 99;
1571 //
1572 if (display_progress) {
1573 cout << "\n No OEM inversion, too high start cost:\n"
1574 << " Set limit : " << max_start_cost << endl
1575 << " Found value : " << cost_start << endl
1576 << endl;
1577 }
1578 }
1579
1580 // Otherwise do inversion
1581 else {
1582 // Size remaining output arguments
1583 x.resize(n);
1584 dxdy.resize(n, m);
1585
1586 OEMVector xa_oem(xa), y_oem(y), x_oem;
1587 oem::AgendaWrapperMPI aw(&ws, &inversion_iterate_agenda, m, n);
1588
1589 OEM_PS_PS_MPI<AgendaWrapperMPI> oem(aw, xa_oem, SaInvMPI, SeInvMPI);
1590
1591 // Call selected method
1592 int return_value = 99;
1593
1594 if (method == "li") {
1595 oem::CG cg(1e-12, 0);
1596 oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1597 return_value = oem.compute<oem::GN_CG, invlib::MPILog>(
1598 x_oem, y_oem, gn, 2 * (int)display_progress);
1599 } else if (method == "gn") {
1600 oem::CG cg(1e-12, 0);
1601 oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1602 return_value = oem.compute<oem::GN_CG, invlib::MPILog>(
1603 x_oem, y_oem, gn, 2 * (int)display_progress);
1604 } else if ((method == "lm") || (method == "ml")) {
1605 oem::CG cg(1e-12, 0);
1606 LM_CG_MPI lm(SaInvMPI, cg);
1607
1608 lm.set_tolerance(stop_dx);
1609 lm.set_maximum_iterations((unsigned int)max_iter);
1610 lm.set_lambda(lm_ga_settings[0]);
1611 lm.set_lambda_decrease(lm_ga_settings[1]);
1612 lm.set_lambda_increase(lm_ga_settings[2]);
1613 lm.set_lambda_threshold(lm_ga_settings[3]);
1614 lm.set_lambda_maximum(lm_ga_settings[4]);
1615
1616 return_value = oem.compute<oem::LM_CG_MPI, invlib::MPILog>(
1617 x_oem, y_oem, lm, 2 * (int)display_progress);
1618 }
1619
1620 oem_diagnostics[0] = return_value;
1621 oem_diagnostics[2] = oem.cost;
1622 oem_diagnostics[3] = oem.cost_y;
1623 oem_diagnostics[4] = static_cast<Numeric>(oem.iterations);
1624
1625 x = x_oem;
1626 // Shall empty jacobian and dxdy be returned?
1627 if (clear_matrices && (oem_diagnostics[0])) {
1628 jacobian.resize(0, 0);
1629 dxdy.resize(0, 0);
1630 }
1631 }
1632 MPI_Finalize();
1633}
1634
1635#else
1636
1638 Vector&,
1639 Vector&,
1640 Matrix&,
1641 Matrix&,
1642 Vector&,
1643 Vector&,
1644 Matrix&,
1645 Matrix&,
1646 Vector&,
1649 const Vector&,
1650 const Vector&,
1651 const Index&,
1653 const Agenda&,
1654 const String&,
1655 const Numeric&,
1656 const Vector&,
1657 const Index&,
1658 const Numeric&,
1659 const Vector&,
1660 const Index&,
1661 const Index&,
1662 const Verbosity&) {
1664 "You have to compile ARTS with OEM support "
1665 " and enable MPI to use OEM_MPI.");
1666}
1667
1668#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:292
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:23250
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:24281
void water_p_eq_agendaExecute(Workspace &ws, Tensor3 &water_p_eq_field, const Tensor3 &t_field, const Agenda &input_agenda)
Definition: auto_md.cc:24619
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:44
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
A constant view of a Tensor3.
Definition: matpackIII.h:132
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
bool empty() const
Check if variable is empty.
Definition: matpackIII.cc:36
Index npages() const
Returns the number of pages.
Definition: matpackIV.cc:58
Index nbooks() const
Returns the number of books.
Definition: matpackIV.cc:55
Index ncols() const
Returns the number of columns.
Definition: matpackIV.cc:64
bool empty() const
Check if variable is empty.
Definition: matpackIV.cc:50
Index nrows() const
Returns the number of rows.
Definition: matpackIV.cc:61
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
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:1225
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056
The range class.
Definition: matpackI.h:165
The Sparse class.
Definition: matpackII.h:67
Vector diagonal() const
Diagonal elements as vector.
Definition: matpackII.cc:172
The Tensor3View class.
Definition: matpackIII.h:239
The Tensor3 class.
Definition: matpackIII.h:339
The Tensor4 class.
Definition: matpackIV.h:421
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
Workspace class.
Definition: workspace_ng.h:40
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:847
void transform_x_back(Vector &x_t, const ArrayOfRetrievalQuantity &jqs, bool revert_functional_transforms)
Handles back-transformations of the state vector.
Definition: jacobian.cc:232
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:45
void transform_x(Vector &x, const ArrayOfRetrievalQuantity &jqs)
Handles transformations of the state vector.
Definition: jacobian.cc:155
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:871
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:2235
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:1369
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:60
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:849
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:1637
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:1387
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:498
void x2artsSpectroscopy(const Verbosity &)
WORKSPACE METHOD: x2artsSpectroscopy.
Definition: m_oem.cc:1004
void avkCalc(Matrix &, const Matrix &, const Matrix &, const Verbosity &)
WORKSPACE METHOD: avkCalc.
Definition: m_oem.cc:1378
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:111
void covmat_soCalc(Matrix &, const Matrix &, const CovarianceMatrix &, const Verbosity &)
WORKSPACE METHOD: covmat_soCalc.
Definition: m_oem.cc:1360
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:86
void reshape(Tensor3View X, ConstVectorView x)
reshape
Definition: math_funcs.cc:735
void flat(VectorView x, ConstMatrixView X)
flat
Definition: math_funcs.cc:681
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
const Joker joker
Particulates Polyfit
Definition: jacobian.h:86
Particulates FrequencyStretch
Definition: jacobian.h:85
Temperature
Definition: jacobian.h:57
Particulates FrequencyShift
Definition: jacobian.h:84
Particulates Sinefit
Definition: jacobian.h:87
MagneticMagnitude
Definition: jacobian.h:59
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
#define w
#define a
#define c
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:99
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