ARTS  2.4.0(git:4fb77825)
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 extern const String ABSSPECIES_MAINTAG;
60 extern const String TEMPERATURE_MAINTAG;
61 extern const String POINTING_MAINTAG;
62 extern const String POINTING_SUBTAG_A;
63 extern const String FREQUENCY_MAINTAG;
64 extern const String FREQUENCY_SUBTAG_0;
65 extern const String FREQUENCY_SUBTAG_1;
66 extern const String POLYFIT_MAINTAG;
67 extern const String SCATSPECIES_MAINTAG;
68 extern const String SINEFIT_MAINTAG;
69 extern const String SURFACE_MAINTAG;
70 extern const String WIND_MAINTAG;
71 extern const String MAGFIELD_MAINTAG;
72 
73 /* Workspace method: Doxygen documentation will be auto-generated */
76  const String& bulkprop_name,
77  const Numeric& limit_low,
78  const Numeric& limit_high,
79  const Verbosity&) {
80  Index iq = -1;
81  if (bulkprop_name == "ALL") {
82  }
83 
84  else {
85  for (Index i = 0; i < particle_bulkprop_names.nelem(); i++) {
86  if (particle_bulkprop_names[i] == bulkprop_name) {
87  iq = i;
88  break;
89  }
90  }
91  if (iq < 0) {
92  ostringstream os;
93  os << "Could not find " << bulkprop_name
94  << " in particle_bulkprop_names.\n";
95  throw std::runtime_error(os.str());
96  }
97  }
98 
99  Tensor4Clip(particle_bulkprop_field, iq, limit_low, limit_high);
100 }
101 
102 /* Workspace method: Doxygen documentation will be auto-generated */
105  const String& species,
106  const Numeric& limit_low,
107  const Numeric& limit_high,
108  const Verbosity&) {
109  Index iq = -1;
110  if (species == "ALL") {
111  }
112 
113  else {
114  for (Index i = 0; i < abs_species.nelem(); i++) {
115  if (abs_species[i][0].Species() == SpeciesTag(species).Species()) {
116  iq = i;
117  break;
118  }
119  }
120  if (iq < 0) {
121  ostringstream os;
122  os << "Could not find " << species << " in abs_species.\n";
123  throw std::runtime_error(os.str());
124  }
125  }
126 
127  Tensor4Clip(vmr_field, iq, limit_low, limit_high);
128 }
129 
130 /* Workspace method: Doxygen documentation will be auto-generated */
131 void xClip(Vector& x,
133  const Index& ijq,
134  const Numeric& limit_low,
135  const Numeric& limit_high,
136  const Verbosity&) {
137  // Sizes
138  const Index nq = jacobian_quantities.nelem();
139 
140  if (ijq < -1) throw runtime_error("Argument *ijq* must be >= -1.");
141  if (ijq >= nq) {
142  ostringstream os;
143  os << "Argument *ijq* is too high.\n"
144  << "You have selected index: " << ijq << "\n"
145  << "but the number of quantities is only: " << nq << "\n"
146  << "(Note that zero-based indexing is used)\n";
147  throw runtime_error(os.str());
148  }
149 
150  // Jacobian indices
152  {
153  bool any_affine;
154  jac_ranges_indices(ji, any_affine, jacobian_quantities);
155  }
156 
157  Index ifirst = 0, ilast = x.nelem() - 1;
158  if (ijq > -1) {
159  ifirst = ji[ijq][0];
160  ilast = ji[ijq][1];
161  }
162 
163  if (!std::isinf(limit_low)) {
164  for (Index i = ifirst; i <= ilast; i++) {
165  if (x[i] < limit_low) x[i] = limit_low;
166  }
167  }
168  if (!std::isinf(limit_high)) {
169  for (Index i = ifirst; i <= ilast; i++) {
170  if (x[i] > limit_high) x[i] = limit_high;
171  }
172  }
173 }
174 
175 /* Workspace method: Doxygen documentation will be auto-generated */
177  Vector& xa,
179  const Index& atmfields_checked,
180  const Index& atmgeom_checked,
181  const Index& atmosphere_dim,
182  const Vector& p_grid,
183  const Vector& lat_grid,
184  const Vector& lon_grid,
185  const Tensor3& t_field,
186  const Tensor4& vmr_field,
188  const Index& cloudbox_on,
189  const Index& cloudbox_checked,
192  const Tensor3& wind_u_field,
193  const Tensor3& wind_v_field,
194  const Tensor3& wind_w_field,
195  const Tensor3& mag_u_field,
196  const Tensor3& mag_v_field,
197  const Tensor3& mag_w_field,
200  const Agenda& water_p_eq_agenda,
201  const Verbosity&) {
202  // Basics
203  //
204  if (atmfields_checked != 1)
205  throw runtime_error(
206  "The atmospheric fields must be flagged to have "
207  "passed a consistency check (atmfields_checked=1).");
208  if (atmgeom_checked != 1)
209  throw runtime_error(
210  "The atmospheric geometry must be flagged to have "
211  "passed a consistency check (atmgeom_checked=1).");
212  if (cloudbox_checked != 1)
213  throw runtime_error(
214  "The cloudbox must be flagged to have "
215  "passed a consistency check (cloudbox_checked=1).");
216 
217  // Jacobian indices
219  {
220  bool any_affine;
221  jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
222  }
223 
224  // Sizes
225  const Index nq = jacobian_quantities.nelem();
226  //
227  xa.resize(ji[nq - 1][1] + 1);
228 
229  // Loop retrieval quantities and fill *xa*
230  for (Index q = 0; q < jacobian_quantities.nelem(); q++) {
231  // Index range of this retrieval quantity
232  const Index np = ji[q][1] - ji[q][0] + 1;
233  Range ind(ji[q][0], np);
234 
235  // Atmospheric temperatures
236  if (jacobian_quantities[q].MainTag() == TEMPERATURE_MAINTAG) {
237  // Here we need to interpolate *t_field*
238  ArrayOfGridPos gp_p, gp_lat, gp_lon;
240  gp_lat,
241  gp_lon,
244  p_grid,
245  lat_grid,
246  lon_grid);
247  Tensor3 t_x;
248  regrid_atmfield_by_gp(t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
249  flat(xa[ind], t_x);
250  }
251 
252  // Abs species
253  else if (jacobian_quantities[q].MainTag() == ABSSPECIES_MAINTAG) {
254  // Index position of species
255  ArrayOfSpeciesTag atag;
257  const Index isp = chk_contains("abs_species", abs_species, atag);
258 
259  if (jacobian_quantities[q].Mode() == "rel") {
260  // This one is simple, just a vector of ones
261  xa[ind] = 1;
262  } else {
263  // For all remaining options we need to interpolate *vmr_field*
264  ArrayOfGridPos gp_p, gp_lat, gp_lon;
266  gp_lat,
267  gp_lon,
270  p_grid,
271  lat_grid,
272  lon_grid);
273  Tensor3 vmr_x;
274  regrid_atmfield_by_gp(vmr_x,
276  vmr_field(isp, joker, joker, joker),
277  gp_p,
278  gp_lat,
279  gp_lon);
280 
281  if (jacobian_quantities[q].Mode() == "vmr") {
282  flat(xa[ind], vmr_x);
283  } else if (jacobian_quantities[q].Mode() == "nd") {
284  // Here we need to also interpolate *t_field*
285  Tensor3 t_x;
287  t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
288  // Calculate number density for species (vmr*nd_tot)
289  Index i = 0;
290  for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
291  for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
292  for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
293  xa[ji[q][0] + i] =
294  vmr_x(i1, i2, i3) *
295  number_density(jacobian_quantities[q].Grids()[0][i1],
296  t_x(i1, i2, i3));
297  i += 1;
298  }
299  }
300  }
301  } else if (jacobian_quantities[q].Mode() == "rh") {
302  // Here we need to also interpolate *t_field*
303  Tensor3 t_x;
305  t_x, atmosphere_dim, t_field, gp_p, gp_lat, gp_lon);
306  Tensor3 water_p_eq;
307  water_p_eq_agendaExecute(ws, water_p_eq, t_x, water_p_eq_agenda);
308  // Calculate relative humidity (vmr*p/p_sat)
309  Index i = 0;
310  for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
311  for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
312  for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
313  xa[ji[q][0] + i] = vmr_x(i1, i2, i3) *
314  jacobian_quantities[q].Grids()[0][i1] /
315  water_p_eq(i1, i2, i3);
316  i += 1;
317  }
318  }
319  }
320  } else if (jacobian_quantities[q].Mode() == "q") {
321  // Calculate specific humidity q, from mixing ratio r and
322  // vapour pressure e, as
323  // q = r(1+r); r = 0.622e/(p-e); e = vmr*p;
324  Index i = 0;
325  for (Index i3 = 0; i3 < vmr_x.ncols(); i3++) {
326  for (Index i2 = 0; i2 < vmr_x.nrows(); i2++) {
327  for (Index i1 = 0; i1 < vmr_x.npages(); i1++) {
328  const Numeric e =
329  vmr_x(i1, i2, i3) * jacobian_quantities[q].Grids()[0][i1];
330  const Numeric r =
331  0.622 * e / (jacobian_quantities[q].Grids()[0][i1] - e);
332  xa[ji[q][0] + i] = r / (1 + r);
333  i += 1;
334  }
335  }
336  }
337  } else {
338  assert(0);
339  }
340  }
341  }
342 
343  // Scattering species
344  else if (jacobian_quantities[q].MainTag() == SCATSPECIES_MAINTAG) {
345  if (cloudbox_on) {
346  if (particle_bulkprop_field.empty()) {
347  throw runtime_error(
348  "One jacobian quantity belongs to the "
349  "scattering species category, but *particle_bulkprop_field* "
350  "is empty.");
351  }
352  if (particle_bulkprop_field.nbooks() !=
353  particle_bulkprop_names.nelem()) {
354  throw runtime_error(
355  "Mismatch in size between "
356  "*particle_bulkprop_field* and *particle_bulkprop_names*.");
357  }
358 
360  jacobian_quantities[q].SubSubtag());
361  if (isp < 0) {
362  ostringstream os;
363  os << "Jacobian quantity with index " << q << " covers a "
364  << "scattering species, and the field quantity is set to \""
365  << jacobian_quantities[q].SubSubtag() << "\", but this quantity "
366  << "could not found in *particle_bulkprop_names*.";
367  throw runtime_error(os.str());
368  }
369 
370  ArrayOfGridPos gp_p, gp_lat, gp_lon;
372  gp_lat,
373  gp_lon,
376  p_grid,
377  lat_grid,
378  lon_grid);
379  Tensor3 pbp_x;
380  regrid_atmfield_by_gp(pbp_x,
383  gp_p,
384  gp_lat,
385  gp_lon);
386  flat(xa[ind], pbp_x);
387  } else {
388  xa[ind] = 0;
389  }
390  }
391 
392  // Wind
393  else if (jacobian_quantities[q].MainTag() == WIND_MAINTAG) {
394  ConstTensor3View source_field(wind_u_field);
395  if (jacobian_quantities[q].Subtag() == "v") {
396  source_field = wind_v_field;
397  } else if (jacobian_quantities[q].Subtag() == "w") {
398  source_field = wind_w_field;
399  }
400 
401  // Determine grid positions for interpolation from retrieval grids back
402  // to atmospheric grids
403  ArrayOfGridPos gp_p, gp_lat, gp_lon;
405  gp_lat,
406  gp_lon,
409  p_grid,
410  lat_grid,
411  lon_grid);
412 
413  Tensor3 wind_x;
415  wind_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
416  flat(xa[ind], wind_x);
417  }
418 
419  // Magnetism
420  else if (jacobian_quantities[q].MainTag() == MAGFIELD_MAINTAG) {
421  if (jacobian_quantities[q].Subtag() == "strength") {
422  // Determine grid positions for interpolation from retrieval grids back
423  // to atmospheric grids
424  ArrayOfGridPos gp_p, gp_lat, gp_lon;
426  gp_lat,
427  gp_lon,
430  p_grid,
431  lat_grid,
432  lon_grid);
433 
434  //all three component's hyoptenuse is the strength
435  Tensor3 mag_u, mag_v, mag_w;
437  mag_u, atmosphere_dim, mag_u_field, gp_p, gp_lat, gp_lon);
439  mag_v, atmosphere_dim, mag_v_field, gp_p, gp_lat, gp_lon);
441  mag_w, atmosphere_dim, mag_w_field, gp_p, gp_lat, gp_lon);
442 
443  Tensor3 mag_x(gp_p.nelem(), gp_lat.nelem(), gp_lon.nelem());
444  for (Index i = 0; i < gp_p.nelem(); i++)
445  for (Index j = 0; j < gp_lat.nelem(); j++)
446  for (Index k = 0; k < gp_lon.nelem(); k++)
447  mag_x(i, j, k) = std::hypot(
448  std::hypot(mag_u(i, j, k), mag_u(i, j, k)),
449  mag_w(i, j, k)); //nb, should remove one hypot for c++17
450  flat(xa[ind], mag_x);
451  } else {
452  ConstTensor3View source_field(mag_u_field);
453  if (jacobian_quantities[q].Subtag() == "v") {
454  source_field = mag_v_field;
455  } else if (jacobian_quantities[q].Subtag() == "w") {
456  source_field = mag_w_field;
457  } else if (jacobian_quantities[q].Subtag() == "u") {
458  } else
459  throw runtime_error("Unsupported magnetism type");
460 
461  // Determine grid positions for interpolation from retrieval grids back
462  // to atmospheric grids
463  ArrayOfGridPos gp_p, gp_lat, gp_lon;
465  gp_lat,
466  gp_lon,
469  p_grid,
470  lat_grid,
471  lon_grid);
472 
473  Tensor3 mag_x;
475  mag_x, atmosphere_dim, source_field, gp_p, gp_lat, gp_lon);
476  flat(xa[ind], mag_x);
477  }
478  }
479 
480  // Surface
481  else if (jacobian_quantities[q].MainTag() == SURFACE_MAINTAG) {
483  lat_grid,
484  lon_grid,
487  if (surface_props_data.empty()) {
488  throw runtime_error(
489  "One jacobian quantity belongs to the "
490  "surface category, but *surface_props_data* is empty.");
491  }
492 
493  const Index isu =
495  if (isu < 0) {
496  ostringstream os;
497  os << "Jacobian quantity with index " << q << " covers a "
498  << "surface property, and the field Subtag is set to \""
499  << jacobian_quantities[q].Subtag() << "\", but this quantity "
500  << "could not found in *surface_props_names*.";
501  throw runtime_error(os.str());
502  }
503 
504  ArrayOfGridPos gp_lat, gp_lon;
505  get_gp_atmsurf_to_rq(gp_lat,
506  gp_lon,
509  lat_grid,
510  lon_grid);
511  Matrix surf_x;
515  gp_lat,
516  gp_lon);
517  flat(xa[ind], surf_x);
518  }
519 
520  // All variables having zero as a priori
521  // ----------------------------------------------------------------------------
522  else if (jacobian_quantities[q].MainTag() == POINTING_MAINTAG ||
523  jacobian_quantities[q].MainTag() == FREQUENCY_MAINTAG ||
524  jacobian_quantities[q].MainTag() == POLYFIT_MAINTAG ||
525  jacobian_quantities[q].MainTag() == SINEFIT_MAINTAG) {
526  xa[ind] = 0;
527  }
528 
529  else {
530  ostringstream os;
531  os << "Found a retrieval quantity that is not yet handled by\n"
532  << "internal retrievals: " << jacobian_quantities[q].MainTag() << endl;
533  throw runtime_error(os.str());
534  }
535  }
536 
537  // Apply transformations
539 }
540 
541 /* Workspace method: Doxygen documentation will be auto-generated */
544  Tensor3& t_field,
554  const Vector& x,
555  const Index& atmfields_checked,
556  const Index& atmgeom_checked,
557  const Index& atmosphere_dim,
558  const Vector& p_grid,
559  const Vector& lat_grid,
560  const Vector& lon_grid,
562  const Index& cloudbox_on,
563  const Index& cloudbox_checked,
566  const Agenda& water_p_eq_agenda,
567  const Verbosity&) {
568  // Basics
569  //
570  if (atmfields_checked != 1)
571  throw runtime_error(
572  "The atmospheric fields must be flagged to have "
573  "passed a consistency check (atmfields_checked=1).");
574  if (atmgeom_checked != 1)
575  throw runtime_error(
576  "The atmospheric geometry must be flagged to have "
577  "passed a consistency check (atmgeom_checked=1).");
578  if (cloudbox_checked != 1)
579  throw runtime_error(
580  "The cloudbox must be flagged to have "
581  "passed a consistency check (cloudbox_checked=1).");
582 
583  // Revert transformation
584  Vector x_t(x);
586 
587  // Main sizes
588  const Index nq = jacobian_quantities.nelem();
589 
590  // Jacobian indices
592  {
593  bool any_affine;
594  jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
595  }
596 
597  // Check input
598  if (x_t.nelem() != ji[nq - 1][1] + 1)
599  throw runtime_error(
600  "Length of *x* does not match length implied by "
601  "*jacobian_quantities*.");
602 
603  // Note that when this method is called, vmr_field and other output variables
604  // have original values, i.e. matching the a priori state.
605 
606  // Loop retrieval quantities
607  for (Index q = 0; q < nq; q++) {
608  // Index range of this retrieval quantity
609  const Index np = ji[q][1] - ji[q][0] + 1;
610  Range ind(ji[q][0], np);
611 
612  // Atmospheric temperatures
613  // ----------------------------------------------------------------------------
614  if (jacobian_quantities[q].MainTag() == TEMPERATURE_MAINTAG) {
615  // Determine grid positions for interpolation from retrieval grids back
616  // to atmospheric grids
617  ArrayOfGridPos gp_p, gp_lat, gp_lon;
618  Index n_p, n_lat, n_lon;
620  gp_lat,
621  gp_lon,
622  n_p,
623  n_lat,
624  n_lon,
625  jacobian_quantities[q].Grids(),
627  p_grid,
628  lat_grid,
629  lon_grid);
630 
631  // Map values in x back to t_field
632  Tensor3 t_x(n_p, n_lat, n_lon);
633  reshape(t_x, x_t[ind]);
635  t_field, atmosphere_dim, t_x, gp_p, gp_lat, gp_lon);
636  }
637 
638  // Abs species
639  // ----------------------------------------------------------------------------
640  else if (jacobian_quantities[q].MainTag() == ABSSPECIES_MAINTAG) {
641  // Index position of species
642  ArrayOfSpeciesTag atag;
644  const Index isp = chk_contains("abs_species", abs_species, atag);
645 
646  // Map part of x to a full atmospheric field
647  Tensor3 x_field(vmr_field.npages(), vmr_field.nrows(), vmr_field.ncols());
648  {
649  ArrayOfGridPos gp_p, gp_lat, gp_lon;
650  Index n_p, n_lat, n_lon;
652  gp_lat,
653  gp_lon,
654  n_p,
655  n_lat,
656  n_lon,
657  jacobian_quantities[q].Grids(),
659  p_grid,
660  lat_grid,
661  lon_grid);
662  //
663  Tensor3 t3_x(n_p, n_lat, n_lon);
664  reshape(t3_x, x_t[ind]);
666  x_field, atmosphere_dim, t3_x, gp_p, gp_lat, gp_lon);
667  }
668  //
669  if (jacobian_quantities[q].Mode() == "rel") {
670  // vmr = vmr0 * x
671  vmr_field(isp, joker, joker, joker) *= x_field;
672  } else if (jacobian_quantities[q].Mode() == "vmr") {
673  // vmr = x
674  vmr_field(isp, joker, joker, joker) = x_field;
675  } else if (jacobian_quantities[q].Mode() == "nd") {
676  // vmr = nd / nd_tot
677  for (Index i3 = 0; i3 < vmr_field.ncols(); i3++) {
678  for (Index i2 = 0; i2 < vmr_field.nrows(); i2++) {
679  for (Index i1 = 0; i1 < vmr_field.npages(); i1++) {
680  vmr_field(isp, i1, i2, i3) =
681  x_field(i1, i2, i3) /
682  number_density(p_grid[i1], t_field(i1, i2, i3));
683  }
684  }
685  }
686  } else if (jacobian_quantities[q].Mode() == "rh") {
687  // vmr = x * p_sat / p
688  Tensor3 water_p_eq;
690  for (Index i3 = 0; i3 < vmr_field.ncols(); i3++) {
691  for (Index i2 = 0; i2 < vmr_field.nrows(); i2++) {
692  for (Index i1 = 0; i1 < vmr_field.npages(); i1++) {
693  vmr_field(isp, i1, i2, i3) =
694  x_field(i1, i2, i3) * water_p_eq(i1, i2, i3) / p_grid[i1];
695  }
696  }
697  }
698  } else if (jacobian_quantities[q].Mode() == "q") {
699  // We have that specific humidity q, mixing ratio r and
700  // vapour pressure e, are related as
701  // q = r(1+r); r = 0.622e/(p-e); e = vmr*p;
702  // That is: vmr=e/p; e = rp/(0.622+r); r = q/(1-q)
703  for (Index i3 = 0; i3 < vmr_field.ncols(); i3++) {
704  for (Index i2 = 0; i2 < vmr_field.nrows(); i2++) {
705  for (Index i1 = 0; i1 < vmr_field.npages(); i1++) {
706  const Numeric r = x_field(i1, i2, i3) / (1 - x_field(i1, i2, i3));
707  const Numeric e = r * p_grid[i1] / (0.622 + r);
708  vmr_field(isp, i1, i2, i3) = e / p_grid[i1];
709  }
710  }
711  }
712  } else {
713  assert(0);
714  }
715  }
716 
717  // Scattering species
718  // ----------------------------------------------------------------------------
719  else if (jacobian_quantities[q].MainTag() == SCATSPECIES_MAINTAG) {
720  // If no cloudbox, we assume that there is nothing to do
721  if (cloudbox_on) {
722  if (particle_bulkprop_field.empty()) {
723  throw runtime_error(
724  "One jacobian quantity belongs to the "
725  "scattering species category, but *particle_bulkprop_field* "
726  "is empty.");
727  }
728  if (particle_bulkprop_field.nbooks() !=
729  particle_bulkprop_names.nelem()) {
730  throw runtime_error(
731  "Mismatch in size between "
732  "*particle_bulkprop_field* and *particle_bulkprop_field*.");
733  }
734 
736  jacobian_quantities[q].SubSubtag());
737  if (isp < 0) {
738  ostringstream os;
739  os << "Jacobian quantity with index " << q << " covers a "
740  << "scattering species, and the field quantity is set to \""
741  << jacobian_quantities[q].SubSubtag() << "\", but this quantity "
742  << "could not found in *particle_bulkprop_names*.";
743  throw runtime_error(os.str());
744  }
745 
746  // Determine grid positions for interpolation from retrieval grids back
747  // to atmospheric grids
748  ArrayOfGridPos gp_p, gp_lat, gp_lon;
749  Index n_p, n_lat, n_lon;
751  gp_lat,
752  gp_lon,
753  n_p,
754  n_lat,
755  n_lon,
756  jacobian_quantities[q].Grids(),
758  p_grid,
759  lat_grid,
760  lon_grid);
761  // Map x to particle_bulkprop_field
762  Tensor3 pbfield_x(n_p, n_lat, n_lon);
763  reshape(pbfield_x, x_t[ind]);
764  Tensor3 pbfield;
766  pbfield, atmosphere_dim, pbfield_x, gp_p, gp_lat, gp_lon);
767  particle_bulkprop_field(isp, joker, joker, joker) = pbfield;
768  }
769  }
770 
771  // Wind
772  // ----------------------------------------------------------------------------
773  else if (jacobian_quantities[q].MainTag() == WIND_MAINTAG) {
774  // Determine grid positions for interpolation from retrieval grids back
775  // to atmospheric grids
776  ArrayOfGridPos gp_p, gp_lat, gp_lon;
777  Index n_p, n_lat, n_lon;
779  gp_lat,
780  gp_lon,
781  n_p,
782  n_lat,
783  n_lon,
784  jacobian_quantities[q].Grids(),
786  p_grid,
787  lat_grid,
788  lon_grid);
789 
790  // TODO Could be done without copying.
791  Tensor3 wind_x(n_p, n_lat, n_lon);
792  reshape(wind_x, x_t[ind]);
793 
794  Tensor3View target_field(wind_u_field);
795 
796  Tensor3 wind_field(
797  target_field.npages(), target_field.nrows(), target_field.ncols());
799  wind_field, atmosphere_dim, wind_x, gp_p, gp_lat, gp_lon);
800 
801  if (jacobian_quantities[q].Subtag() == "u") {
802  wind_u_field = wind_field;
803  } else if (jacobian_quantities[q].Subtag() == "v") {
804  wind_v_field = wind_field;
805  } else if (jacobian_quantities[q].Subtag() == "w") {
806  wind_w_field = wind_field;
807  }
808  }
809 
810  // Magnetism
811  // ----------------------------------------------------------------------------
812  else if (jacobian_quantities[q].MainTag() == MAGFIELD_MAINTAG) {
813  // Determine grid positions for interpolation from retrieval grids back
814  // to atmospheric grids
815  ArrayOfGridPos gp_p, gp_lat, gp_lon;
816  Index n_p, n_lat, n_lon;
818  gp_lat,
819  gp_lon,
820  n_p,
821  n_lat,
822  n_lon,
823  jacobian_quantities[q].Grids(),
825  p_grid,
826  lat_grid,
827  lon_grid);
828 
829  // TODO Could be done without copying.
830  Tensor3 mag_x(n_p, n_lat, n_lon);
831  reshape(mag_x, x_t[ind]);
832 
833  Tensor3View target_field(mag_u_field);
834 
835  Tensor3 mag_field(
836  target_field.npages(), target_field.nrows(), target_field.ncols());
838  mag_field, atmosphere_dim, mag_x, gp_p, gp_lat, gp_lon);
839  if (jacobian_quantities[q].Subtag() == "u") {
840  mag_u_field = mag_field;
841  } else if (jacobian_quantities[q].Subtag() == "v") {
842  mag_v_field = mag_field;
843  } else if (jacobian_quantities[q].Subtag() == "w") {
844  mag_w_field = mag_field;
845  } else if (jacobian_quantities[q].Subtag() == "strength") {
846  for (Index i = 0; i < n_p; i++) {
847  for (Index j = 0; j < n_lat; j++) {
848  for (Index k = 0; k < n_lon; k++) {
849  Numeric scale =
850  mag_x(i, j, k) /
851  std::hypot(
852  std::hypot(mag_u_field(i, j, k), mag_v_field(i, j, k)),
853  mag_w_field(i, j, k)); // nb,remove one hypot for c++17
854  mag_u_field(i, j, k) *= scale;
855  mag_v_field(i, j, k) *= scale;
856  mag_w_field(i, j, k) *= scale;
857  }
858  }
859  }
860  } else
861  throw runtime_error("Unsupported magnetism type");
862  }
863 
864  // Surface
865  // ----------------------------------------------------------------------------
866  else if (jacobian_quantities[q].MainTag() == SURFACE_MAINTAG) {
868  lat_grid,
869  lon_grid,
872  if (surface_props_data.empty()) {
873  throw runtime_error(
874  "One jacobian quantity belongs to the "
875  "surface category, but *surface_props_data* is empty.");
876  }
877 
878  const Index isu =
880  if (isu < 0) {
881  ostringstream os;
882  os << "Jacobian quantity with index " << q << " covers a "
883  << "surface property, and the field Subtag is set to \""
884  << jacobian_quantities[q].Subtag() << "\", but this quantity "
885  << "could not found in *surface_props_names*.";
886  throw runtime_error(os.str());
887  }
888 
889  // Determine grid positions for interpolation from retrieval grids back
890  // to atmospheric grids
891  ArrayOfGridPos gp_lat, gp_lon;
892  Index n_lat, n_lon;
893  get_gp_rq_to_atmgrids(gp_lat,
894  gp_lon,
895  n_lat,
896  n_lon,
897  jacobian_quantities[q].Grids(),
899  lat_grid,
900  lon_grid);
901  // Map values in x back to surface_props_data
902  Matrix surf_x(n_lat, n_lon);
903  reshape(surf_x, x_t[ind]);
904  Matrix surf;
905  regrid_atmsurf_by_gp_oem(surf, atmosphere_dim, surf_x, gp_lat, gp_lon);
906  surface_props_data(isu, joker, joker) = surf;
907  }
908  }
909 }
910 
911 /* Workspace method: Doxygen documentation will be auto-generated */
914  Vector& f_backend,
925  const Vector& x,
927  const Index& sensor_checked,
928  const Vector& sensor_time,
929  const Verbosity&) {
930  // Basics
931  //
932  if (sensor_checked != 1)
933  throw runtime_error(
934  "The sensor response must be flagged to have "
935  "passed a consistency check (sensor_checked=1).");
936 
937  // Revert transformation
938  Vector x_t(x);
940 
941  // Main sizes
942  const Index nq = jacobian_quantities.nelem();
943 
944  // Jacobian indices
946  {
947  bool any_affine;
948  jac_ranges_indices(ji, any_affine, jacobian_quantities, true);
949  }
950 
951  // Check input
952  if (x_t.nelem() != ji[nq - 1][1] + 1)
953  throw runtime_error(
954  "Length of *x* does not match length implied by "
955  "*jacobian_quantities*.");
956 
957  // Flag indicating that y_baseline is not set
958  bool yb_set = false;
959 
960  // Shall sensor responses be calculed?
961  bool do_sensor = false;
962 
963  // Loop retrieval quantities
964  for (Index q = 0; q < nq; q++) {
965  // Index range of this retrieval quantity
966  const Index np = ji[q][1] - ji[q][0] + 1;
967 
968  // Pointing off-set
969  // ----------------------------------------------------------------------------
970  if (jacobian_quantities[q].MainTag() == POINTING_MAINTAG) {
971  if (jacobian_quantities[q].Subtag() != POINTING_SUBTAG_A) {
972  ostringstream os;
973  os << "Only pointing off-sets treated by *jacobianAddPointingZa* "
974  << "are so far handled.";
975  throw runtime_error(os.str());
976  }
977  // Handle pointing "jitter" seperately
978  if (jacobian_quantities[q].Grids()[0][0] == -1) {
979  if (sensor_los.nrows() != np)
980  throw runtime_error(
981  "Mismatch between pointing jacobian and *sensor_los* found.");
982  // Simply add retrieved off-set(s) to za column of *sensor_los*
983  for (Index i = 0; i < np; i++) {
984  sensor_los(i, 0) += x_t[ji[q][0] + i];
985  }
986  }
987  // Polynomial representation
988  else {
989  if (sensor_los.nrows() != sensor_time.nelem())
990  throw runtime_error(
991  "Sizes of *sensor_los* and *sensor_time* do not match.");
992  Vector w;
993  for (Index c = 0; c < np; c++) {
995  for (Index i = 0; i < w.nelem(); i++) {
996  sensor_los(i, 0) += w[i] * x_t[ji[q][0] + c];
997  }
998  }
999  }
1000  }
1001 
1002  // Frequncy shift or stretch
1003  // ----------------------------------------------------------------------------
1004  else if (jacobian_quantities[q].MainTag() == FREQUENCY_MAINTAG) {
1005  if (jacobian_quantities[q].Subtag() == FREQUENCY_SUBTAG_0) {
1006  assert(np == 1);
1007  if (x_t[ji[q][0]] != 0) {
1008  do_sensor = true;
1009  f_backend += x_t[ji[q][0]];
1010  }
1011  } else if (jacobian_quantities[q].Subtag() == FREQUENCY_SUBTAG_1) {
1012  assert(np == 1);
1013  if (x_t[ji[q][0]] != 0) {
1014  do_sensor = true;
1015  Vector w;
1017  for (Index i = 0; i < w.nelem(); i++) {
1018  f_backend[i] += w[i] * x_t[ji[q][0]];
1019  }
1020  }
1021  } else {
1022  assert(0);
1023  }
1024  }
1025 
1026  // Baseline fit: polynomial or sinusoidal
1027  // ----------------------------------------------------------------------------
1028  else if (jacobian_quantities[q].MainTag() == POLYFIT_MAINTAG ||
1029  jacobian_quantities[q].MainTag() == SINEFIT_MAINTAG) {
1030  if (!yb_set) {
1031  yb_set = true;
1032  Index y_size = sensor_los.nrows() * sensor_response_f_grid.nelem() *
1033  sensor_response_pol_grid.nelem() *
1034  sensor_response_dlos_grid.nrows();
1035  y_baseline.resize(y_size);
1036  y_baseline = 0;
1037  }
1038 
1039  for (Index mb = 0; mb < sensor_los.nrows(); ++mb) {
1041  x_t,
1042  mb,
1048  q,
1049  ji);
1050  }
1051  }
1052  }
1053 
1054  // *y_baseline* not yet set?
1055  if (!yb_set) {
1056  y_baseline.resize(1);
1057  y_baseline[0] = 0;
1058  }
1059 
1060  // Recalculate sensor_response?
1061  if (do_sensor) {
1071  f_backend,
1073  }
1074 }
1075 
1076 /* Workspace method: Doxygen documentation will be auto-generated */
1078  throw runtime_error("Retrievals of spectroscopic variables not yet handled.");
1079 }
1080 
1081 
1082 #ifdef OEM_SUPPORT
1083 /* Workspace method: Doxygen documentation will be auto-generated */
1084 void OEM(Workspace& ws,
1085  Vector& x,
1086  Vector& yf,
1087  Matrix& jacobian,
1088  Matrix& dxdy,
1091  ArrayOfString& errors,
1092  const Vector& xa,
1093  const CovarianceMatrix& covmat_sx,
1094  const Vector& y,
1095  const CovarianceMatrix& covmat_se,
1098  const String& method,
1099  const Numeric& max_start_cost,
1100  const Vector& x_norm,
1101  const Index& max_iter,
1102  const Numeric& stop_dx,
1103  const Vector& lm_ga_settings,
1104  const Index& clear_matrices,
1105  const Index& display_progress,
1106  const Verbosity&) {
1107  // Main sizes
1108  const Index n = covmat_sx.nrows();
1109  const Index m = y.nelem();
1110 
1111  // Checks
1112  covmat_sx.compute_inverse();
1113  covmat_se.compute_inverse();
1114 
1115  OEM_checks(ws,
1116  x,
1117  yf,
1118  jacobian,
1120  xa,
1121  covmat_sx,
1122  y,
1123  covmat_se,
1125  method,
1126  x_norm,
1127  max_iter,
1128  stop_dx,
1129  lm_ga_settings,
1130  clear_matrices,
1131  display_progress);
1132 
1133  // Size diagnostic output and init with NaNs
1134  oem_diagnostics.resize(5);
1135  oem_diagnostics = NAN;
1136  //
1137  if (method == "ml" || method == "lm" || method == "ml_cg" ||
1138  method == "lm_cg") {
1139  lm_ga_history.resize(max_iter + 1);
1140  lm_ga_history = NAN;
1141  } else {
1142  lm_ga_history.resize(0);
1143  }
1144 
1145  // Check for start vector and precomputed yf, jacobian
1146  if (x.nelem() != n) {
1147  x = xa;
1148  yf.resize(0);
1149  jacobian.resize(0, 0);
1150  }
1151 
1152  // If no precomputed value given, we compute yf and jacobian to
1153  // compute initial cost (and use in the first OEM iteration).
1154  if (yf.nelem() == 0) {
1156  ws, yf, jacobian, xa, 1, 0, inversion_iterate_agenda);
1157  }
1158 
1159  if (yf.nelem() not_eq y.nelem()) {
1160  std::ostringstream os;
1161  os << "Mismatch between simulated y and input y.\n";
1162  os << "Input y is size " << y.nelem() << " but simulated y is "
1163  << yf.nelem() << "\n";
1164  os << "Use your frequency grid vector and your sensor response matrix to match simulations with measurements.\n";
1165  throw std::runtime_error(os.str());
1166  }
1167 
1168  // TODO: Get this from invlib log.
1169  // Start value of cost function
1170  Numeric cost_start = NAN;
1171  if (method == "ml" || method == "lm" || display_progress ||
1172  max_start_cost > 0) {
1173  Vector dy = y;
1174  dy -= yf;
1175  Vector sdy = y;
1176  mult(sdy, covmat_se, dy);
1177  Vector dx = x;
1178  dx -= xa;
1179  Vector sdx = x;
1180  mult(sdx, covmat_sx, dx);
1181  cost_start = dx * sdx + dy * sdy;
1182  cost_start /= static_cast<Numeric>(m);
1183  }
1184  oem_diagnostics[1] = cost_start;
1185 
1186  // Handle cases with too large start cost
1187  if (max_start_cost > 0 && cost_start > max_start_cost) {
1188  // Flag no inversion in oem_diagnostics, and let x to be undefined
1189  oem_diagnostics[0] = 99;
1190  //
1191  if (display_progress) {
1192  cout << "\n No OEM inversion, too high start cost:\n"
1193  << " Set limit : " << max_start_cost << endl
1194  << " Found value : " << cost_start << endl
1195  << endl;
1196  }
1197  }
1198  // Otherwise do inversion
1199  else {
1200  bool apply_norm = false;
1201  oem::Matrix T{};
1202  if (x_norm.nelem() == n) {
1203  T.resize(n, n);
1204  T *= 0.0;
1205  T.diagonal() = x_norm;
1206  for (Index i = 0; i < n; i++) {
1207  T(i, i) = x_norm[i];
1208  }
1209  apply_norm = true;
1210  }
1211 
1213  oem::Vector xa_oem(xa), y_oem(y), x_oem(x);
1214  oem::AgendaWrapper aw(&ws,
1215  (unsigned int)m,
1216  (unsigned int)n,
1217  jacobian,
1218  yf,
1220  oem::OEM_STANDARD<oem::AgendaWrapper> oem(aw, xa_oem, Sa, Se);
1221  oem::OEM_MFORM<oem::AgendaWrapper> oem_m(aw, xa_oem, Sa, Se);
1222  int oem_verbosity = static_cast<int>(display_progress);
1223 
1224  int return_code = 0;
1225 
1226  try {
1227  if (method == "li") {
1228  oem::Std s(T, apply_norm);
1229  oem::GN gn(stop_dx, 1, s); // Linear case, only one step.
1230  return_code = oem.compute<oem::GN, oem::ArtsLog>(
1231  x_oem, y_oem, gn, oem_verbosity, lm_ga_history, true);
1232  oem_diagnostics[0] = static_cast<Index>(return_code);
1233  } else if (method == "li_m") {
1234  oem::Std s(T, apply_norm);
1235  oem::GN gn(stop_dx, 1, s); // Linear case, only one step.
1236  return_code = oem_m.compute<oem::GN, oem::ArtsLog>(
1237  x_oem, y_oem, gn, oem_verbosity, lm_ga_history, true);
1238  oem_diagnostics[0] = static_cast<Index>(return_code);
1239  } else if (method == "li_cg") {
1240  oem::CG cg(T, apply_norm, 1e-10, 0);
1241  oem::GN_CG gn(stop_dx, 1, cg); // Linear case, only one step.
1242  return_code = oem.compute<oem::GN_CG, oem::ArtsLog>(
1243  x_oem, y_oem, gn, oem_verbosity, lm_ga_history, true);
1244  oem_diagnostics[0] = static_cast<Index>(return_code);
1245  } else if (method == "li_cg_m") {
1246  oem::CG cg(T, apply_norm, 1e-10, 0);
1247  oem::GN_CG gn(stop_dx, 1, cg); // Linear case, only one step.
1248  return_code = oem_m.compute<oem::GN_CG, oem::ArtsLog>(
1249  x_oem, y_oem, gn, oem_verbosity, lm_ga_history, true);
1250  oem_diagnostics[0] = static_cast<Index>(return_code);
1251  } else if (method == "gn") {
1252  oem::Std s(T, apply_norm);
1253  oem::GN gn(stop_dx, (unsigned int)max_iter, s);
1254  return_code = oem.compute<oem::GN, oem::ArtsLog>(
1255  x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1256  oem_diagnostics[0] = static_cast<Index>(return_code);
1257  } else if (method == "gn_m") {
1258  oem::Std s(T, apply_norm);
1259  oem::GN gn(stop_dx, (unsigned int)max_iter, s);
1260  return_code = oem_m.compute<oem::GN, oem::ArtsLog>(
1261  x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1262  oem_diagnostics[0] = static_cast<Index>(return_code);
1263  } else if (method == "gn_cg") {
1264  oem::CG cg(T, apply_norm, 1e-10, 0);
1265  oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1266  return_code = oem.compute<oem::GN_CG, oem::ArtsLog>(
1267  x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1268  oem_diagnostics[0] = static_cast<Index>(return_code);
1269  } else if (method == "gn_cg_m") {
1270  oem::CG cg(T, apply_norm, 1e-10, 0);
1271  oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1272  return_code = oem_m.compute<oem::GN_CG, oem::ArtsLog>(
1273  x_oem, y_oem, gn, oem_verbosity, lm_ga_history);
1274  oem_diagnostics[0] = static_cast<Index>(return_code);
1275  } else if ((method == "lm") || (method == "ml")) {
1276  oem::Std s(T, apply_norm);
1277  Sparse diagonal = Sparse::diagonal(covmat_sx.inverse_diagonal());
1278  CovarianceMatrix SaDiag{};
1279  SaDiag.add_correlation_inverse(Block(Range(0, n),
1280  Range(0, n),
1281  std::make_pair(0, 0),
1282  make_shared<Sparse>(diagonal)));
1283  oem::CovarianceMatrix SaInvLM = inv(oem::CovarianceMatrix(SaDiag));
1284  oem::LM lm(SaInvLM, s);
1285 
1286  lm.set_tolerance(stop_dx);
1287  lm.set_maximum_iterations((unsigned int)max_iter);
1288  lm.set_lambda(lm_ga_settings[0]);
1289  lm.set_lambda_decrease(lm_ga_settings[1]);
1290  lm.set_lambda_increase(lm_ga_settings[2]);
1291  lm.set_lambda_maximum(lm_ga_settings[3]);
1292  lm.set_lambda_threshold(lm_ga_settings[4]);
1293  lm.set_lambda_constraint(lm_ga_settings[5]);
1294 
1295  return_code = oem.compute<oem::LM&, oem::ArtsLog>(
1296  x_oem, y_oem, lm, oem_verbosity, lm_ga_history);
1297  oem_diagnostics[0] = static_cast<Index>(return_code);
1298  if (lm.get_lambda() > lm.get_lambda_maximum()) {
1299  oem_diagnostics[0] = 2;
1300  }
1301  } else if ((method == "lm_cg") || (method == "ml_cg")) {
1302  oem::CG cg(T, apply_norm, 1e-10, 0);
1303 
1304  Sparse diagonal = Sparse::diagonal(covmat_sx.inverse_diagonal());
1305  CovarianceMatrix SaDiag{};
1306  SaDiag.add_correlation_inverse(Block(Range(0, n),
1307  Range(0, n),
1308  std::make_pair(0, 0),
1309  make_shared<Sparse>(diagonal)));
1310  oem::LM_CG lm(SaDiag, cg);
1311 
1312  lm.set_maximum_iterations((unsigned int)max_iter);
1313  lm.set_lambda(lm_ga_settings[0]);
1314  lm.set_lambda_decrease(lm_ga_settings[1]);
1315  lm.set_lambda_increase(lm_ga_settings[2]);
1316  lm.set_lambda_threshold(lm_ga_settings[3]);
1317  lm.set_lambda_maximum(lm_ga_settings[4]);
1318 
1319  return_code = oem.compute<oem::LM_CG&, oem::ArtsLog>(
1320  x_oem, y_oem, lm, oem_verbosity, lm_ga_history);
1321  oem_diagnostics[0] = static_cast<Index>(return_code);
1322  if (lm.get_lambda() > lm.get_lambda_maximum()) {
1323  oem_diagnostics[0] = 2;
1324  }
1325  }
1326 
1327  oem_diagnostics[2] = oem.cost / static_cast<Numeric>(m);
1328  oem_diagnostics[3] = oem.cost_y / static_cast<Numeric>(m);
1329  oem_diagnostics[4] = static_cast<Numeric>(oem.iterations);
1330  } catch (const std::exception& e) {
1331  oem_diagnostics[0] = 9;
1332  oem_diagnostics[2] = oem.cost;
1333  oem_diagnostics[3] = oem.cost_y;
1334  oem_diagnostics[4] = static_cast<Numeric>(oem.iterations);
1335  x_oem *= NAN;
1336  std::vector<std::string> sv = oem::handle_nested_exception(e);
1337  for (auto& s : sv) {
1338  std::stringstream ss{s};
1339  std::string t{};
1340  while (std::getline(ss, t)) {
1341  errors.push_back(t.c_str());
1342  }
1343  }
1344  } catch (...) {
1345  throw;
1346  }
1347 
1348  x = x_oem;
1349  yf = aw.get_measurement_vector();
1350 
1351  // Shall empty jacobian and dxdy be returned?
1352  if (clear_matrices) {
1353  jacobian.resize(0, 0);
1354  dxdy.resize(0, 0);
1355  } else if (oem_diagnostics[0] <= 2) {
1356  dxdy.resize(n, m);
1357  Matrix tmp1(n, m), tmp2(n, n), tmp3(n, n);
1359  mult(tmp2, tmp1, jacobian);
1360  add_inv(tmp2, covmat_sx);
1361  inv(tmp3, tmp2);
1362  mult(dxdy, tmp3, tmp1);
1363  }
1364  }
1365 }
1366 
1367 /* Workspace method: Doxygen documentation will be auto-generated */
1369  const Matrix& dxdy,
1370  const CovarianceMatrix& covmat_se,
1371  const Verbosity& /*v*/) {
1372  Index n(dxdy.nrows()), m(dxdy.ncols());
1373  Matrix tmp1(m, n);
1374 
1375  if ((m == 0) || (n == 0)) {
1376  throw runtime_error(
1377  "The gain matrix *dxdy* is required to compute the observation error covariance matrix.");
1378  }
1379  if ((covmat_se.nrows() != m) || (covmat_se.ncols() != m)) {
1380  throw runtime_error(
1381  "The covariance matrix covmat_se has invalid dimensions.");
1382  }
1383 
1384  covmat_so.resize(n, n);
1385  mult(tmp1, covmat_se, transpose(dxdy));
1386  mult(covmat_so, dxdy, tmp1);
1387 }
1388 
1389 /* Workspace method: Doxygen documentation will be auto-generated */
1391  const Matrix& avk,
1392  const CovarianceMatrix& covmat_sx,
1393  const Verbosity& /*v*/) {
1394  Index n(avk.ncols());
1395  Matrix tmp1(n, n), tmp2(n, n);
1396 
1397  if (n == 0) {
1398  throw runtime_error(
1399  "The averaging kernel matrix *dxdy* is required to compute the smoothing error covariance matrix.");
1400  }
1401  if ((covmat_sx.nrows() != n) || (covmat_sx.ncols() != n)) {
1402  throw runtime_error(
1403  "The covariance matrix *covmat_sx* invalid dimensions.");
1404  }
1405 
1406  covmat_ss.resize(n, n);
1407 
1408  // Sign doesn't matter since we're dealing with a quadratic form.
1409  id_mat(tmp1);
1410  tmp1 -= avk;
1411 
1412  mult(tmp2, covmat_sx, transpose(tmp1));
1413  mult(covmat_ss, tmp1, tmp2);
1414 }
1415 
1416 /* Workspace method: Doxygen documentation will be auto-generated */
1418  const CovarianceMatrix& Sc,
1419  const Verbosity& /*v*/) {
1420  S = Matrix(Sc);
1421 }
1422 
1423 /* Workspace method: Doxygen documentation will be auto-generated */
1424 void avkCalc(Matrix& avk,
1425  const Matrix& dxdy,
1426  const Matrix& jacobian,
1427  const Verbosity& /*v*/) {
1428  Index m(jacobian.nrows()), n(jacobian.ncols());
1429  if ((m == 0) || (n == 0))
1430  throw runtime_error("The Jacobian matrix is empty.");
1431  if ((dxdy.nrows() != n) || (dxdy.ncols() != m)) {
1432  ostringstream os;
1433  os << "Matrices have inconsistent sizes.\n"
1434  << " Size of gain matrix: " << dxdy.nrows() << " x " << dxdy.ncols()
1435  << "\n"
1436  << " Size of Jacobian: " << jacobian.nrows() << " x "
1437  << jacobian.ncols() << "\n";
1438  throw runtime_error(os.str());
1439  }
1440 
1441  avk.resize(n, n);
1442  mult(avk, dxdy, jacobian);
1443 }
1444 
1445 #else
1446 
1447 void covmat_soCalc(Matrix& /* covmat_so */,
1448  const Matrix& /* dxdy */,
1449  const CovarianceMatrix& /* covmat_se*/,
1450  const Verbosity& /*v*/) {
1451  throw runtime_error(
1452  "WSM is not available because ARTS was compiled without "
1453  "OEM support.");
1454 }
1455 
1456 void covmat_ssCalc(Matrix& /*covmat_ss*/,
1457  const Matrix& /*avk*/,
1458  const CovarianceMatrix& /*covmat_sx*/,
1459  const Verbosity& /*v*/) {
1460  throw runtime_error(
1461  "WSM is not available because ARTS was compiled without "
1462  "OEM support.");
1463 }
1464 
1465 void avkCalc(Matrix& /* avk */,
1466  const Matrix& /* dxdy */,
1467  const Matrix& /* jacobian */,
1468  const Verbosity& /*v*/) {
1469  throw runtime_error(
1470  "WSM is not available because ARTS was compiled without "
1471  "OEM support.");
1472 }
1473 
1475  Vector&,
1476  Vector&,
1477  Matrix&,
1478  Matrix&,
1479  Vector&,
1480  Vector&,
1481  ArrayOfString&,
1482  const Vector&,
1483  const CovarianceMatrix&,
1484  const Vector&,
1485  const CovarianceMatrix&,
1486  const Index&,
1487  const ArrayOfRetrievalQuantity&,
1488  const ArrayOfArrayOfIndex&,
1489  const Agenda&,
1490  const String&,
1491  const Numeric&,
1492  const Vector&,
1493  const Index&,
1494  const Numeric&,
1495  const Vector&,
1496  const Index&,
1497  const Index&,
1498  const Verbosity&) {
1499  throw runtime_error(
1500  "WSM is not available because ARTS was compiled without "
1501  "OEM support.");
1502 }
1503 #endif
1504 
1505 #if defined(OEM_SUPPORT) && 0
1506 
1507 #include "agenda_wrapper_mpi.h"
1508 #include "oem_mpi.h"
1509 
1510 //
1511 // Performs manipulations of workspace variables necessary for distributed
1512 // retrievals with MPI:
1513 //
1514 // - Splits up sensor positions evenly over processes
1515 // - Splits up inverse covariance matrices.
1516 //
1517 void MPI_Initialize(Matrix& sensor_los,
1518  Matrix& sensor_pos,
1519  Vector& sensor_time) {
1520  int initialized;
1521 
1522  MPI_Initialized(&initialized);
1523  if (!initialized) {
1524  MPI_Init(nullptr, nullptr);
1525  }
1526 
1527  int rank, nprocs;
1528 
1529  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
1530  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
1531 
1532  int nmblock = (int)sensor_pos.nrows();
1533  int mblock_range = nmblock / nprocs;
1534  int mblock_start = mblock_range * rank;
1535  int remainder = nmblock % std::max(mblock_range, nprocs);
1536 
1537  //
1538  // Split up sensor positions.
1539  //
1540 
1541  if (rank < remainder) {
1542  mblock_range += 1;
1543  mblock_start += rank;
1544  } else {
1545  mblock_start += remainder;
1546  }
1547 
1548  if (nmblock > 0) {
1549  Range range = Range(mblock_start, mblock_range);
1550 
1551  Matrix tmp_m = sensor_los(range, joker);
1552  sensor_los = tmp_m;
1553 
1554  tmp_m = sensor_pos(range, joker);
1555  sensor_pos = tmp_m;
1556 
1557  Vector tmp_v = sensor_time[range];
1558  sensor_time = tmp_v;
1559  } else {
1560  sensor_los.resize(0, 0);
1561  sensor_pos.resize(0, 0);
1562  sensor_time.resize(0);
1563  }
1564 }
1565 
1566 void OEM_MPI(Workspace& ws,
1567  Vector& x,
1568  Vector& yf,
1569  Matrix& jacobian,
1570  Matrix& dxdy,
1573  Matrix& sensor_los,
1574  Matrix& sensor_pos,
1578  const Vector& xa,
1579  const Vector& y,
1582  const String& method,
1583  const Numeric& max_start_cost,
1584  const Vector& x_norm,
1585  const Index& max_iter,
1586  const Numeric& stop_dx,
1587  const Vector& lm_ga_settings,
1588  const Index& clear_matrices,
1589  const Index& display_progress,
1590  const Verbosity& /*v*/) {
1591  // Main sizes
1592  const Index n = covmat_sx.nrows();
1593  const Index m = y.nelem();
1594 
1595  // Check WSVs
1596  OEM_checks(ws,
1597  x,
1598  yf,
1599  jacobian,
1601  xa,
1602  covmat_sx,
1603  covmat_se,
1605  method,
1606  x_norm,
1607  max_iter,
1608  stop_dx,
1609  lm_ga_settings,
1610  clear_matrices,
1611  display_progress);
1612 
1613  // Calculate spectrum and Jacobian for a priori state
1614  // Jacobian is also input to the agenda, and to flag this is this first
1615  // call, this WSV must be set to be empty.
1616  jacobian.resize(0, 0);
1617 
1618  // Initialize MPI environment.
1619  MPI_Initialize(sensor_los, sensor_pos, sensor_time);
1620 
1621  // Setup distributed matrices.
1622  MPICovarianceMatrix SeInvMPI(covmat_se);
1623  MPICovarianceMatrix SaInvMPI(covmat_sx);
1624 
1625  // Create temporary MPI vector from local results and use conversion to
1626  // standard vector to broadcast results to all processes.
1627  oem::Vector tmp;
1629  ws, tmp, jacobian, xa, 1, inversion_iterate_agenda);
1630  yf = MPIVector(tmp);
1631 
1632  // Size diagnostic output and init with NaNs
1633  oem_diagnostics.resize(5);
1634  oem_diagnostics = NAN;
1635  //
1636  if (method == "ml" || method == "lm") {
1637  lm_ga_history.resize(max_iter);
1638  lm_ga_history = NAN;
1639  } else {
1640  lm_ga_history.resize(0);
1641  }
1642 
1643  // Start value of cost function. Covariance matrices are already distributed
1644  // over processes, so we need to use invlib matrix algebra.
1645  Numeric cost_start = NAN;
1646  if (method == "ml" || method == "lm" || display_progress ||
1647  max_start_cost > 0) {
1648  oem::Vector dy = y;
1649  dy -= yf;
1650  cost_start = dot(dy, SeInvMPI * dy);
1651  }
1652  oem_diagnostics[1] = cost_start;
1653 
1654  // Handle cases with too large start cost
1655  if (max_start_cost > 0 && cost_start > max_start_cost) {
1656  // Flag no inversion in oem_diagnostics, and let x to be undefined
1657  oem_diagnostics[0] = 99;
1658  //
1659  if (display_progress) {
1660  cout << "\n No OEM inversion, too high start cost:\n"
1661  << " Set limit : " << max_start_cost << endl
1662  << " Found value : " << cost_start << endl
1663  << endl;
1664  }
1665  }
1666 
1667  // Otherwise do inversion
1668  else {
1669  // Size remaining output arguments
1670  x.resize(n);
1671  dxdy.resize(n, m);
1672 
1673  OEMVector xa_oem(xa), y_oem(y), x_oem;
1675 
1676  OEM_PS_PS_MPI<AgendaWrapperMPI> oem(aw, xa_oem, SaInvMPI, SeInvMPI);
1677 
1678  // Call selected method
1679  int return_value = 99;
1680 
1681  if (method == "li") {
1682  oem::CG cg(1e-12, 0);
1683  oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1684  return_value = oem.compute<oem::GN_CG, invlib::MPILog>(
1685  x_oem, y_oem, gn, 2 * (int)display_progress);
1686  } else if (method == "gn") {
1687  oem::CG cg(1e-12, 0);
1688  oem::GN_CG gn(stop_dx, (unsigned int)max_iter, cg);
1689  return_value = oem.compute<oem::GN_CG, invlib::MPILog>(
1690  x_oem, y_oem, gn, 2 * (int)display_progress);
1691  } else if ((method == "lm") || (method == "ml")) {
1692  oem::CG cg(1e-12, 0);
1693  LM_CG_MPI lm(SaInvMPI, cg);
1694 
1695  lm.set_tolerance(stop_dx);
1696  lm.set_maximum_iterations((unsigned int)max_iter);
1697  lm.set_lambda(lm_ga_settings[0]);
1698  lm.set_lambda_decrease(lm_ga_settings[1]);
1699  lm.set_lambda_increase(lm_ga_settings[2]);
1700  lm.set_lambda_threshold(lm_ga_settings[3]);
1701  lm.set_lambda_maximum(lm_ga_settings[4]);
1702 
1703  return_value = oem.compute<oem::LM_CG_MPI, invlib::MPILog>(
1704  x_oem, y_oem, lm, 2 * (int)display_progress);
1705  }
1706 
1707  oem_diagnostics[0] = return_value;
1708  oem_diagnostics[2] = oem.cost;
1709  oem_diagnostics[3] = oem.cost_y;
1710  oem_diagnostics[4] = static_cast<Numeric>(oem.iterations);
1711 
1712  x = x_oem;
1713  // Shall empty jacobian and dxdy be returned?
1714  if (clear_matrices && (oem_diagnostics[0])) {
1715  jacobian.resize(0, 0);
1716  dxdy.resize(0, 0);
1717  }
1718  }
1719  MPI_Finalize();
1720 }
1721 
1722 #else
1723 
1725  Vector&,
1726  Vector&,
1727  Matrix&,
1728  Matrix&,
1729  Vector&,
1730  Vector&,
1731  Matrix&,
1732  Matrix&,
1733  Vector&,
1736  const Vector&,
1737  const Vector&,
1738  const Index&,
1739  const ArrayOfRetrievalQuantity&,
1740  const Agenda&,
1741  const String&,
1742  const Numeric&,
1743  const Vector&,
1744  const Index&,
1745  const Numeric&,
1746  const Vector&,
1747  const Index&,
1748  const Index&,
1749  const Verbosity&) {
1750  throw runtime_error(
1751  "You have to compile ARTS with OEM support "
1752  " and enable MPI to use OEM_MPI.");
1753 }
1754 
1755 #endif // OEM_SUPPORT && ENABLE_MPI
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
ARTS::Var::wind_v_field
Tensor3 wind_v_field(Workspace &ws) noexcept
Definition: autoarts.h:7246
SCATSPECIES_MAINTAG
const String SCATSPECIES_MAINTAG
ARTS::Var::covmat_ss
Matrix covmat_ss(Workspace &ws) noexcept
Definition: autoarts.h:2919
ARTS::Var::sensor_response_agenda
Agenda sensor_response_agenda(Workspace &ws) noexcept
Definition: autoarts.h:6322
ARTS::Var::sensor_response_dlos_grid
Matrix sensor_response_dlos_grid(Workspace &ws) noexcept
Definition: autoarts.h:6356
auto_md.h
vmr_fieldClip
void vmr_fieldClip(Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const String &species, const Numeric &limit_low, const Numeric &limit_high, const Verbosity &)
WORKSPACE METHOD: vmr_fieldClip.
Definition: m_oem.cc:103
ARTS::Var::atmgeom_checked
Index atmgeom_checked(Workspace &ws) noexcept
Definition: autoarts.h:2489
ARTS::Var::jacobian_quantities
ArrayOfRetrievalQuantity jacobian_quantities(Workspace &ws) noexcept
Definition: autoarts.h:3924
array_species_tag_from_string
void array_species_tag_from_string(ArrayOfSpeciesTag &tags, const String &names)
Converts a String to ArrayOfSpeciesTag.
Definition: abs_species_tags.cc:586
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
surface.h
id_mat
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2240
ARTS::Var::sensor_response
Sparse sensor_response(Workspace &ws) noexcept
Definition: autoarts.h:6300
ARTS::Var::dxdy
Matrix dxdy(Workspace &ws) noexcept
Definition: autoarts.h:3368
SINEFIT_MAINTAG
const String SINEFIT_MAINTAG
oem::AgendaWrapper
Interface to ARTS inversion_iterate_agenda.
Definition: oem.h:462
regrid_atmfield_by_gp
void regrid_atmfield_by_gp(Tensor3 &field_new, const Index &atmosphere_dim, ConstTensor3View field_old, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regrids an atmospheric field, for precalculated grid positions.
Definition: special_interp.cc:310
ARTS::Var::lat_grid
Vector lat_grid(Workspace &ws) noexcept
Definition: autoarts.h:3962
joker
const Joker joker
surface_props_check
void surface_props_check(const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names)
Peforms basic checks of surface_props_data and surface_props_names
Definition: surface.cc:139
inv
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
oem::Matrix
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:34
ARTS::Var::y_baseline
Vector y_baseline(Workspace &ws) noexcept
Definition: autoarts.h:7447
avkCalc
void avkCalc(Matrix &, const Matrix &, const Matrix &, const Verbosity &)
WORKSPACE METHOD: avkCalc.
Definition: m_oem.cc:1465
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
oem::MPICovarianceMatrix
invlib::Matrix< invlib::MPIMatrix< invlib::Timer< ArtsCovarianceMatrixWrapper > >> MPICovarianceMatrix
MPI-distributed covariance matrix type.
Definition: oem_mpi.h:27
ARTS::Var::xa
Vector xa(Workspace &ws) noexcept
Definition: autoarts.h:7359
Sparse
The Sparse class.
Definition: matpackII.h:60
CovarianceMatrix
Definition: covariance_matrix.h:226
ARTS::Var::wind_u_field
Tensor3 wind_u_field(Workspace &ws) noexcept
Definition: autoarts.h:7211
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
polynomial_basis_func
void polynomial_basis_func(Vector &b, const Vector &x, const Index &poly_coeff)
Calculates polynomial basis functions.
Definition: jacobian.cc:897
OEM_checks
void OEM_checks(Workspace &ws, Vector &x, Vector &yf, Matrix &jacobian, const Agenda &inversion_iterate_agenda, const Vector &xa, const CovarianceMatrix &covmat_sx, const Vector &y, const CovarianceMatrix &covmat_se, const ArrayOfRetrievalQuantity &jacobian_quantities, const String &method, const Vector &x_norm, const Index &max_iter, const Numeric &stop_dx, const Vector &lm_ga_settings, const Index &clear_matrices, const Index &display_progress)
Error checking for OEM method.
Definition: oem.h:666
regrid_atmsurf_by_gp_oem
void regrid_atmsurf_by_gp_oem(Matrix &field_new, const Index &atmosphere_dim, ConstMatrixView field_old, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regridding of surface field OEM-type.
Definition: special_interp.cc:658
ARTS::Var::mag_w_field
Tensor3 mag_w_field(Workspace &ws) noexcept
Definition: autoarts.h:4197
Species
QuantumIdentifier::QType Index LowerQuantumNumbers Species
Definition: arts_api_classes.cc:255
covmat_ssCalc
void covmat_ssCalc(Matrix &, const Matrix &, const CovarianceMatrix &, const Verbosity &)
WORKSPACE METHOD: covmat_ssCalc.
Definition: m_oem.cc:1456
FREQUENCY_SUBTAG_0
const String FREQUENCY_SUBTAG_0
oem::CovarianceMatrix
invlib::Matrix< ArtsCovarianceMatrixWrapper > CovarianceMatrix
invlib wrapper type for ARTS the ARTS covariance class.
Definition: oem.h:38
ARTS::Var::inversion_iterate_agenda
Agenda inversion_iterate_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3635
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
array.h
This file contains the definition of Array.
FREQUENCY_SUBTAG_1
const String FREQUENCY_SUBTAG_1
Agenda
The Agenda class.
Definition: agenda_class.h:44
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
transform_x_back
void transform_x_back(Vector &x_t, const ArrayOfRetrievalQuantity &jqs, bool revert_functional_transforms)
Handles back-transformations of the state vector.
Definition: jacobian.cc:257
ARTS::Var::particle_bulkprop_names
ArrayOfString particle_bulkprop_names(Workspace &ws) noexcept
Definition: autoarts.h:4823
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:239
x2artsAtmAndSurf
void x2artsAtmAndSurf(Workspace &ws, Tensor4 &vmr_field, Tensor3 &t_field, Tensor4 &particle_bulkprop_field, Tensor3 &wind_u_field, Tensor3 &wind_v_field, Tensor3 &wind_w_field, Tensor3 &mag_u_field, Tensor3 &mag_v_field, Tensor3 &mag_w_field, Tensor3 &surface_props_data, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &x, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const ArrayOfArrayOfSpeciesTag &abs_species, const Index &cloudbox_on, const Index &cloudbox_checked, const ArrayOfString &particle_bulkprop_names, const ArrayOfString &surface_props_names, const Agenda &water_p_eq_agenda, const Verbosity &)
WORKSPACE METHOD: x2artsAtmAndSurf.
Definition: m_oem.cc:542
get_gp_atmgrids_to_rq
void get_gp_atmgrids_to_rq(ArrayOfGridPos &gp_p, ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, const RetrievalQuantity &rq, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric fields to retrieval grids.
Definition: special_interp.cc:366
Array< String >
oem::GN
invlib::GaussNewton< Numeric, Std > GN
OEM Gauss-Newton optimization using normed ARTS QR solver.
Definition: oem.h:169
SpeciesTag
A tag group can consist of the sum of several of these.
Definition: abs_species_tags.h:44
xaStandard
void xaStandard(Workspace &ws, Vector &xa, const ArrayOfRetrievalQuantity &jacobian_quantities, 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 Tensor3 &t_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Index &cloudbox_on, const Index &cloudbox_checked, const Tensor4 &particle_bulkprop_field, const ArrayOfString &particle_bulkprop_names, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Tensor3 &surface_props_data, const ArrayOfString &surface_props_names, const Agenda &water_p_eq_agenda, const Verbosity &)
WORKSPACE METHOD: xaStandard.
Definition: m_oem.cc:176
flat
void flat(VectorView x, ConstMatrixView X)
flat
Definition: math_funcs.cc:707
number_density
Numeric number_density(const Numeric &p, const Numeric &t)
number_density
Definition: physics_funcs.cc:219
ARTS::Var::covmat_so
Matrix covmat_so(Workspace &ws) noexcept
Definition: autoarts.h:2905
POINTING_SUBTAG_A
const String POINTING_SUBTAG_A
ARTS::Var::sensor_response_f
Vector sensor_response_f(Workspace &ws) noexcept
Definition: autoarts.h:6376
x2artsSpectroscopy
void x2artsSpectroscopy(const Verbosity &)
WORKSPACE METHOD: x2artsSpectroscopy.
Definition: m_oem.cc:1077
oem::GN_CG
invlib::GaussNewton< Numeric, CG > GN_CG
Gauss-Newton (GN) optimization using normed CG solver.
Definition: oem.h:171
inversion_iterate_agendaExecute
void inversion_iterate_agendaExecute(Workspace &ws, Vector &yf, Matrix &jacobian, const Vector &x, const Index jacobian_do, const Index inversion_iteration_counter, const Agenda &input_agenda)
Definition: auto_md.cc:23951
transform_x
void transform_x(Vector &x, const ArrayOfRetrievalQuantity &jqs)
Handles transformations of the state vector.
Definition: jacobian.cc:168
MAGFIELD_MAINTAG
const String MAGFIELD_MAINTAG
WIND_MAINTAG
const String WIND_MAINTAG
reshape
void reshape(Tensor3View X, ConstVectorView x)
reshape
Definition: math_funcs.cc:761
ARTS::Var::sensor_time
Vector sensor_time(Workspace &ws) noexcept
Definition: autoarts.h:6444
my_basic_string< char >
ARTS::Var::lm_ga_history
Vector lm_ga_history(Workspace &ws) noexcept
Definition: autoarts.h:4024
oem::ArtsLog
OEM log output.
Definition: oem.h:244
oem::OEM_MFORM
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::MFORM > OEM_MFORM
OEM m form.
Definition: oem.h:95
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
CovarianceMatrix::add_correlation_inverse
void add_correlation_inverse(Block c)
Add block inverse to covariance matrix.
Definition: covariance_matrix.cc:512
ARTS::Var::vmr_field
Tensor4 vmr_field(Workspace &ws) noexcept
Definition: autoarts.h:7130
physics_funcs.h
This file contains declerations of functions of physical character.
ARTS::Var::abs_species
ArrayOfArrayOfSpeciesTag abs_species(Workspace &ws) noexcept
Definition: autoarts.h:2157
POINTING_MAINTAG
const String POINTING_MAINTAG
ARTS::Var::atmfields_checked
Index atmfields_checked(Workspace &ws) noexcept
Definition: autoarts.h:2473
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
ARTS::Var::surface_props_data
Tensor3 surface_props_data(Workspace &ws) noexcept
Definition: autoarts.h:6745
jacobian.h
Routines for setting up the jacobian.
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:49
xClip
void xClip(Vector &x, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &ijq, const Numeric &limit_low, const Numeric &limit_high, const Verbosity &)
WORKSPACE METHOD: xClip.
Definition: m_oem.cc:131
oem::handle_nested_exception
std::vector< std::string > handle_nested_exception(const E &e, int level=0)
Handle exception encountered within invlib.
Definition: oem.h:425
ARTS::Var::yf
Vector yf(Workspace &ws) noexcept
Definition: autoarts.h:7658
oem.h
Defines the ARTS interface to the invlib library.
ARTS::Var::covmat_sx
CovarianceMatrix covmat_sx(Workspace &ws) noexcept
Definition: autoarts.h:2945
math_funcs.h
ARTS::Var::f_backend
Vector f_backend(Workspace &ws) noexcept
Definition: autoarts.h:3418
ARTS::Var::covmat_se
CovarianceMatrix covmat_se(Workspace &ws) noexcept
Definition: autoarts.h:2879
ARTS::Var::sensor_pos
Matrix sensor_pos(Workspace &ws) noexcept
Definition: autoarts.h:6279
max
#define max(a, b)
Definition: legacy_continua.cc:20629
ARTS::Var::sensor_response_f_grid
Vector sensor_response_f_grid(Workspace &ws) noexcept
Definition: autoarts.h:6392
calcBaselineFit
void calcBaselineFit(Vector &y_baseline, const Vector &x, const Index &mblock_index, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const RetrievalQuantity &rq, const Index rq_index, const ArrayOfArrayOfIndex &jacobian_indices)
Calculate baseline fit.
Definition: jacobian.cc:921
POLYFIT_MAINTAG
const String POLYFIT_MAINTAG
sensor_response_agendaExecute
void sensor_response_agendaExecute(Workspace &ws, Sparse &sensor_response, Vector &sensor_response_f, Vector &sensor_response_f_grid, ArrayOfIndex &sensor_response_pol, ArrayOfIndex &sensor_response_pol_grid, Matrix &sensor_response_dlos, Matrix &sensor_response_dlos_grid, Matrix &mblock_dlos_grid, const Vector &f_backend, const Agenda &input_agenda)
Definition: auto_md.cc:24922
mult_inv
void mult_inv(MatrixView C, ConstMatrixView A, const CovarianceMatrix &B)
Definition: covariance_matrix.cc:574
ConstTensor3View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
ARTS::Var::lon_grid
Vector lon_grid(Workspace &ws) noexcept
Definition: autoarts.h:4090
Tensor4Clip
void Tensor4Clip(Tensor4 &x, const Index &iq, const Numeric &limit_low, const Numeric &limit_high)
Clip Tensor4.
Definition: oem.h:582
Sparse::diagonal
Vector diagonal() const
Diagonal elements as vector.
Definition: matpackII.cc:172
water_p_eq_agendaExecute
void water_p_eq_agendaExecute(Workspace &ws, Tensor3 &water_p_eq_field, const Tensor3 &t_field, const Agenda &input_agenda)
Definition: auto_md.cc:25260
oem::Vector
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
ARTS::Var::mblock_dlos_grid
Matrix mblock_dlos_grid(Workspace &ws) noexcept
Definition: autoarts.h:4247
Range
The range class.
Definition: matpackI.h:160
ARTS::Var::sensor_response_pol
ArrayOfIndex sensor_response_pol(Workspace &ws) noexcept
Definition: autoarts.h:6411
oem
Definition: oem.h:29
ARTS::Var::sensor_response_pol_grid
ArrayOfIndex sensor_response_pol_grid(Workspace &ws) noexcept
Definition: autoarts.h:6428
q
#define q
Definition: legacy_continua.cc:21712
ARTS::Var::oem_diagnostics
Vector oem_diagnostics(Workspace &ws) noexcept
Definition: autoarts.h:4718
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
ARTS::Var::particle_bulkprop_field
Tensor4 particle_bulkprop_field(Workspace &ws) noexcept
Definition: autoarts.h:4810
OEM
void OEM(Workspace &, Vector &, Vector &, Matrix &, Matrix &, Vector &, Vector &, ArrayOfString &, const Vector &, const CovarianceMatrix &, const Vector &, const CovarianceMatrix &, const Index &, const ArrayOfRetrievalQuantity &, const ArrayOfArrayOfIndex &, const Agenda &, const String &, const Numeric &, const Vector &, const Index &, const Numeric &, const Vector &, const Index &, const Index &, const Verbosity &)
Definition: m_oem.cc:1474
rte.h
Declaration of functions in rte.cc.
ARTS::Var::mag_v_field
Tensor3 mag_v_field(Workspace &ws) noexcept
Definition: autoarts.h:4164
ARTS::Var::t_field
Tensor3 t_field(Workspace &ws) noexcept
Definition: autoarts.h:6947
Workspace
Workspace class.
Definition: workspace_ng.h:40
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
oem::OEM_STANDARD
invlib::MAP< ForwardModel, Matrix, CovarianceMatrix, CovarianceMatrix, Vector, Formulation::STANDARD, invlib::Rodgers531 > OEM_STANDARD
OEM standard form.
Definition: oem.h:63
add_inv
void add_inv(MatrixView A, const CovarianceMatrix &B)
Definition: covariance_matrix.cc:611
dx
#define dx
Definition: legacy_continua.cc:22180
oem::LM_CG
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, CG > LM_CG
Levenberg-Marquardt (LM) optimization using normed CG solver.
Definition: oem.h:175
special_interp.h
Header file for special_interp.cc.
oem::AgendaWrapperMPI
Interface for distributed ARTS forward model.
Definition: oem_mpi.h:56
TEMPERATURE_MAINTAG
const String TEMPERATURE_MAINTAG
x2artsSensor
void x2artsSensor(Workspace &ws, Matrix &sensor_los, Vector &f_backend, Vector &y_baseline, Sparse &sensor_response, Vector &sensor_response_f, ArrayOfIndex &sensor_response_pol, Matrix &sensor_response_dlos, Vector &sensor_response_f_grid, ArrayOfIndex &sensor_response_pol_grid, Matrix &sensor_response_dlos_grid, Matrix &mblock_dlos_grid, const ArrayOfRetrievalQuantity &jacobian_quantities, const Vector &x, const Agenda &sensor_response_agenda, const Index &sensor_checked, const Vector &sensor_time, const Verbosity &)
WORKSPACE METHOD: x2artsSensor.
Definition: m_oem.cc:912
oem::NormalizingSolver
Normalizing solver.
Definition: oem.h:113
ARTS::Var::wind_w_field
Tensor3 wind_w_field(Workspace &ws) noexcept
Definition: autoarts.h:7279
Block
Definition: covariance_matrix.h:62
chk_contains
Index chk_contains(const String &x_name, const Array< T > &x, const T &what)
Check if an array contains a value.
Definition: check_input.h:149
particle_bulkprop_fieldClip
void particle_bulkprop_fieldClip(Tensor4 &particle_bulkprop_field, const ArrayOfString &particle_bulkprop_names, const String &bulkprop_name, const Numeric &limit_low, const Numeric &limit_high, const Verbosity &)
WORKSPACE METHOD: particle_bulkprop_fieldClip.
Definition: m_oem.cc:74
ARTS::Var::avk
Matrix avk(Workspace &ws) noexcept
Definition: autoarts.h:2522
OEM_MPI
void OEM_MPI(Workspace &, Vector &, Vector &, Matrix &, Matrix &, Vector &, Vector &, Matrix &, Matrix &, Vector &, CovarianceMatrix &, CovarianceMatrix &, const Vector &, const Vector &, const Index &, const ArrayOfRetrievalQuantity &, const Agenda &, const String &, const Numeric &, const Vector &, const Index &, const Numeric &, const Vector &, const Index &, const Index &, const Verbosity &)
Definition: m_oem.cc:1724
ARTS::Var::sensor_checked
Index sensor_checked(Workspace &ws) noexcept
Definition: autoarts.h:6144
MatrixFromCovarianceMatrix
void MatrixFromCovarianceMatrix(Matrix &out, const CovarianceMatrix &in, const Verbosity &verbosity)
WORKSPACE METHOD: MatrixFromCovarianceMatrix.
ARTS::Var::jacobian
Matrix jacobian(Workspace &ws) noexcept
Definition: autoarts.h:3892
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
covmat_soCalc
void covmat_soCalc(Matrix &, const Matrix &, const CovarianceMatrix &, const Verbosity &)
WORKSPACE METHOD: covmat_soCalc.
Definition: m_oem.cc:1447
ARTS::Var::sensor_response_dlos
Matrix sensor_response_dlos(Workspace &ws) noexcept
Definition: autoarts.h:6340
transpose
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
Definition: complex.cc:1509
find_first
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:292
get_gp_rq_to_atmgrids
void get_gp_rq_to_atmgrids(ArrayOfGridPos &gp_p, ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, Index &n_p, Index &n_lat, Index &n_lon, const ArrayOfVector &ret_grids, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric fields to retrieval grids.
Definition: special_interp.cc:413
ARTS::Var::mag_u_field
Tensor3 mag_u_field(Workspace &ws) noexcept
Definition: autoarts.h:4130
oem::LM
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, Std > LM
Levenberg-Marquardt (LM) optimization using normed ARTS QR solver.
Definition: oem.h:173
Vector
The Vector class.
Definition: matpackI.h:860
arts_omp.h
Header file for helper functions for OpenMP.
SURFACE_MAINTAG
const String SURFACE_MAINTAG
get_gp_atmsurf_to_rq
void get_gp_atmsurf_to_rq(ArrayOfGridPos &gp_lat, ArrayOfGridPos &gp_lon, const RetrievalQuantity &rq, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid)
Determines grid positions for regridding of atmospheric surfaces to retrieval grids.
Definition: special_interp.cc:392
jac_ranges_indices
void jac_ranges_indices(ArrayOfArrayOfIndex &jis, bool &any_affine, const ArrayOfRetrievalQuantity &jqs, const bool &before_affine)
Determines the index range inside x and the Jacobian for each retrieval quantity.
Definition: jacobian.cc:58
ARTS::Var::water_p_eq_agenda
Agenda water_p_eq_agenda(Workspace &ws) noexcept
Definition: autoarts.h:7165
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
ARTS::Var::sensor_los
Matrix sensor_los(Workspace &ws) noexcept
Definition: autoarts.h:6201
ARTS::Var::cloudbox_checked
Index cloudbox_checked(Workspace &ws) noexcept
Definition: autoarts.h:2650
regrid_atmfield_by_gp_oem
void regrid_atmfield_by_gp_oem(Tensor3 &field_new, const Index &atmosphere_dim, ConstTensor3View field_old, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Regridding of atmospheric field OEM-type.
Definition: special_interp.cc:507
FREQUENCY_MAINTAG
const String FREQUENCY_MAINTAG
oem::MPIVector
invlib::Vector< invlib::MPIVector< invlib::Timer< ArtsVector > >> MPIVector
MPI-distributed vector type.
Definition: oem_mpi.h:29
ABSSPECIES_MAINTAG
const String ABSSPECIES_MAINTAG
arts.h
The global header file for ARTS.
ARTS::Var::cloudbox_on
Index cloudbox_on(Workspace &ws) noexcept
Definition: autoarts.h:2782
oem_mpi.h
Optimal estimation method for MPI-distributed retrieval.
ARTS::Var::surface_props_names
ArrayOfString surface_props_names(Workspace &ws) noexcept
Definition: autoarts.h:6759