ARTS  2.3.1285(git:92a29ea9-dirty)
m_retrieval.cc
Go to the documentation of this file.
1 #include <algorithm>
2 
3 #include "auto_md.h"
4 #include "covariance_matrix.h"
5 #include "jacobian.h"
6 
8 // Helper Functions
10 
12  const RetrievalQuantity& jq,
13  const Index rq_index,
14  const Index grid_dimensions,
15  const Sparse& covmat_block,
16  const Sparse& covmat_inv_block) {
17  Index start = covmat.nrows();
18  Index extent = covmat_block.nrows();
19  Range range(start, extent);
20  Index n_gps = 1;
21  for (Index j = 0; j < grid_dimensions; ++j) {
22  n_gps *= jq.Grids()[j].nelem();
23  }
24 
25  if (!covmat_block.empty()) {
26  if ((n_gps == extent) && (n_gps == covmat_block.ncols())) {
27  std::shared_ptr<Sparse> mat = make_shared<Sparse>(covmat_block);
28  covmat.add_correlation(
29  Block(range, range, std::make_pair(rq_index, rq_index), mat));
30  } else {
31  ostringstream os;
32  os << "The matrix in covmat_block was expected to have dimensions ["
33  << n_gps << ", " << n_gps << "] but found to have dimensions ["
34  << covmat_block.nrows() << ", " << covmat_block.ncols() << "].";
35  throw runtime_error(os.str());
36  }
37  }
38  if (!covmat_inv_block.empty()) {
39  if ((n_gps == covmat_inv_block.nrows()) &&
40  (n_gps == covmat_inv_block.ncols())) {
41  std::shared_ptr<Sparse> mat = make_shared<Sparse>(covmat_inv_block);
43  Block(range, range, std::make_pair(rq_index, rq_index), mat));
44  } else {
45  ostringstream os;
46  os << "The matrix in covmat_inv_block was expected to have dimensions ["
47  << n_gps << ", " << n_gps << "] but found to have dimensions ["
48  << covmat_block.nrows() << ", " << covmat_block.ncols() << "].";
49  throw runtime_error(os.str());
50  }
51  }
52 }
53 
55  ArrayOfRetrievalQuantity& jacobian_quantities,
56  Numeric var) {
57  Index i = jacobian_quantities.size() - 1;
58  Index start = covmat.nrows();
59  Index extent = 1;
60  Range range(start, extent);
61  std::shared_ptr<Matrix> mat = make_shared<Matrix>(1, 1);
62  mat->operator()(0, 0) = var;
63  covmat.add_correlation(Block(range, range, std::make_pair(i, i), mat));
64  mat = make_shared<Matrix>(1, 1);
65  mat->operator()(0, 0) = 1.0 / var;
67  Block(range, range, std::make_pair(i, i), mat));
68 }
69 
71 // Simple Pressure-Height Conversion
73 
74 void ZFromPSimple(Vector& z_grid, const Vector& p_grid, const Verbosity&) {
75  z_grid = Vector(p_grid.nelem());
76 
77  for (Index i = 0; i < p_grid.nelem(); ++i) {
78  if (p_grid[i] < 0.01) {
79  throw runtime_error("Pressures below 0.01 Pa are not accedpted.");
80  }
81  }
82 
83  for (Index i = 0; i < p_grid.nelem(); ++i) {
84  z_grid[i] = 16e3 * (5.0 - log10(p_grid[i]));
85  }
86 }
87 
88 void PFromZSimple(Vector& p_grid, const Vector& z_grid, const Verbosity&) {
89  p_grid = Vector(z_grid.nelem());
90 
91  for (Index i = 0; i < p_grid.nelem(); ++i) {
92  if (z_grid[i] > 120e3) {
93  throw runtime_error("Altitudes above 120 km are not accepted.");
94  }
95  }
96 
97  for (Index i = 0; i < z_grid.nelem(); ++i) {
98  p_grid[i] = pow(10.0, 5.0 - z_grid[i] / 16e3);
99  }
100 }
101 
103 // Creation of Correlation Blocks
105 
106 template <typename MatrixType>
108  Index /*m*/,
109  Index /*n*/,
110  const ArrayOfIndex& row_indices,
111  const ArrayOfIndex& column_indices,
112  const Vector& elements) {
113  matrix.insert_elements(
114  row_indices.nelem(), row_indices, column_indices, elements);
115 }
116 
117 template <>
118 void insert_elements(Matrix& matrix,
119  Index m,
120  Index n,
121  const ArrayOfIndex& row_indices,
122  const ArrayOfIndex& column_indices,
123  const Vector& elements) {
124  assert(row_indices.nelem() == column_indices.nelem());
125  assert(column_indices.nelem() == elements.nelem());
126 
127  matrix.resize(m, n);
128 
129  for (Index i = 0; i < elements.nelem(); ++i) {
130  matrix(row_indices[i], column_indices[i]) = elements[i];
131  }
132 }
133 
134 template <typename MatrixType>
136  MatrixType& block_inv,
137  const Vector& vars,
138  const Verbosity&) {
139  if (vars.empty()) {
140  throw runtime_error(
141  "Cannot pass empty vector of variances to covmat_blockSetDiagonal");
142  }
143  Index n = vars.nelem();
144  block = MatrixType(n, n);
145  block_inv = MatrixType(n, n);
146 
147  ArrayOfIndex indices(n);
148  Vector elements(n), elements_inv(n);
149  for (Index i = 0; i < n; ++i) {
150  indices[i] = i;
151  elements[i] = vars[i];
152  elements_inv[i] = 1.0 / vars[i];
153  }
154 
155  insert_elements(block, n, n, indices, indices, elements);
156  insert_elements(block_inv, n, n, indices, indices, elements_inv);
157 }
158 
159 template void covmatDiagonal(Matrix& block,
160  Matrix& block_inv,
161  const Vector& vars,
162  const Verbosity&);
163 template void covmatDiagonal(Sparse& block,
164  Sparse& block_inv,
165  const Vector& vars,
166  const Verbosity&);
167 
168 template <typename MatrixType>
170  MatrixType& block_inv,
171  const Vector& grid,
172  const Vector& sigma,
173  const Numeric& lc,
174  const Numeric& /*co*/,
175  const Verbosity&) {
176  Index n = grid.nelem();
177  ArrayOfIndex row_indices, column_indices;
178  Vector elements, elements_inv;
179 
180  if (n != sigma.nelem()) {
181  throw runtime_error("Size of grid incompatible with given variances.");
182  }
183 
184  elements = Vector(row_indices.size());
185  for (size_t i = 0; i < row_indices.size(); ++i) {
186  Numeric dz = abs(grid[row_indices[i]] - grid[column_indices[i]]);
187  Numeric e =
188  sigma[row_indices[i]] * sigma[column_indices[i]] * exp(-dz / lc);
189  elements[i] = e;
190  }
191 
192  block = MatrixType(n, n);
193  insert_elements(block, n, n, row_indices, column_indices, elements);
194  row_indices = ArrayOfIndex{};
195  column_indices = ArrayOfIndex{};
196  elements = Vector(3 * n - 2);
197 
198  Numeric dz = abs(grid[1] - grid[0]);
199  Numeric alpha = exp(-dz / lc);
200  Numeric c1 = -alpha / (1.0 - alpha * alpha);
201  Numeric c2 = 1.0 / (1.0 - alpha * alpha);
202 
203  for (Index i = 0; i < n; ++i) {
204  // Lower Tri-diagonal
205  if (i > 0) {
206  column_indices.push_back(i - 1);
207  row_indices.push_back(i);
208  elements[i * 3 - 1] = c1;
209  elements[i * 3 - 1] /= (sigma[i] * sigma[i - 1]);
210  }
211 
212  // Diagonal Elements
213  column_indices.push_back(i);
214  row_indices.push_back(i);
215  if (i == 0 || i == n - 1) {
216  elements[i * 3] = c2;
217  } else {
218  elements[i * 3] = c2 * (1.0 + alpha * alpha);
219  }
220  elements[i * 3] /= (sigma[i] * sigma[i]);
221 
222  // Upper Tri-diagonal
223  if (i < n - 1) {
224  column_indices.push_back(i + 1);
225  row_indices.push_back(i);
226  elements[i * 3 + 1] = c1;
227  elements[i * 3 + 1] /= (sigma[i] * sigma[i + 1]);
228  }
229  }
230  block_inv = MatrixType(n, n);
231  insert_elements(block_inv, n, n, row_indices, column_indices, elements);
232 }
233 
234 template void covmat1DMarkov(Matrix& block,
235  Matrix& block_inv,
236  const Vector& grid,
237  const Vector& sigma,
238  const Numeric& lc,
239  const Numeric& /*co*/,
240  const Verbosity&);
241 
242 template void covmat1DMarkov(Sparse& block,
243  Sparse& block_inv,
244  const Vector& grid,
245  const Vector& sigma,
246  const Numeric& lc,
247  const Numeric& /*co*/,
248  const Verbosity&);
249 
250 template <typename MatrixType>
251 void covmat1D(MatrixType& block,
252  const Vector& grid1,
253  const Vector& grid2,
254  const Vector& sigma1,
255  const Vector& sigma2,
256  const Vector& lc1,
257  const Vector& lc2,
258  const Numeric& co,
259  const String& fname,
260  const Verbosity&) {
261  Index m = grid1.nelem();
262  Vector sigma1_copy(sigma1), lc1_copy(lc1);
263 
264  assert((sigma1.nelem() == m) || (sigma1.nelem() == 1));
265  if (sigma1.nelem() == 1) {
266  Numeric v = sigma1[0];
267  sigma1_copy = Vector(m);
268  sigma1_copy = v;
269  }
270 
271  assert((lc1.nelem() == m) || (lc1.nelem() == 1));
272  if (lc1.nelem() == 1) {
273  Numeric v = lc1[0];
274  lc1_copy = Vector(m);
275  lc1_copy = v;
276  };
277 
278  Index n = grid2.nelem();
279  Vector sigma2_copy(sigma2), lc2_copy(lc2);
280 
281  assert((sigma2.nelem() == n) || (sigma2.nelem() == 1));
282  if (sigma2.nelem() == 1) {
283  Numeric v = sigma2[0];
284  sigma2_copy = Vector(n);
285  sigma2_copy = v;
286  }
287 
288  assert((lc2.nelem() == n) || (lc2.nelem() == 1));
289  if (lc2.nelem() == 1) {
290  Numeric v = lc2[0];
291  lc2_copy = Vector(n);
292  lc2_copy = v;
293  };
294 
295  ConstVectorView grid_view_1(grid1), sigma_view_1(sigma1_copy),
296  cls_view_1(lc1_copy);
297  ConstVectorView grid_view_2(grid2), sigma_view_2(sigma2_copy),
298  cls_view_2(lc2_copy);
299 
300  if (n == 0) {
301  n = m;
302  grid_view_2 = grid_view_1;
303  sigma_view_2 = sigma_view_1;
304  cls_view_2 = cls_view_1;
305  }
306 
307  // Correlation Functions
308  auto f_lin = [&](Index i, Index j) {
309  Numeric d = abs(grid_view_1[i] - grid_view_2[j]);
310  Numeric cl = 0.5 * (cls_view_1[i] + cls_view_2[j]);
311  return 1.0 - (1.0 - exp(-1.0)) * (d / cl);
312  };
313 
314  auto f_exp = [&](Index i, Index j) {
315  Numeric d = abs(grid_view_1[i] - grid_view_2[j]);
316  Numeric cl = 0.5 * (cls_view_1[i] + cls_view_2[j]);
317  return exp(-d / cl);
318  };
319 
320  auto f_gau = [&](Index i, Index j) {
321  Numeric d = abs(grid_view_1[i] - grid_view_2[j]);
322  Numeric cl = 0.5 * (cls_view_1[i] + cls_view_2[j]);
323  Numeric x = d / cl;
324  x *= x;
325  return exp(-x);
326  };
327 
328  ArrayOfIndex row_indices;
329  ArrayOfIndex column_indices;
330 
331  row_indices.reserve(n * m);
332  column_indices.reserve(n * m);
333 
334  std::function<Numeric(Index, Index)> f;
335 
336  if (fname == "exp") {
337  f = f_exp;
338  } else if (fname == "lin") {
339  f = f_lin;
340  } else if (fname == "gau") {
341  f = f_gau;
342  } else {
343  ostringstream os;
344  os << fname << " is not a known function name. Supported names"
345  << "are: exp, lin, gau.";
346  std::runtime_error(os.str());
347  }
348 
349  for (Index i = 0; i < m; ++i) {
350  for (Index j = 0; j < n; ++j) {
351  Numeric e = f(i, j);
352  if (e >= co) {
353  row_indices.push_back(i);
354  column_indices.push_back(j);
355  }
356  }
357  }
358 
359  Vector elements(row_indices.size());
360  for (size_t i = 0; i < row_indices.size(); ++i) {
361  Index ii = row_indices[i];
362  Index jj = column_indices[i];
363  elements[i] = sigma_view_1[ii] * sigma_view_2[jj] * f(ii, jj);
364  }
365 
366  block = MatrixType(m, n);
367  insert_elements(block, m, n, row_indices, column_indices, elements);
368 }
369 
370 template void covmat1D(Matrix& block,
371  const Vector& grid1,
372  const Vector& grid2,
373  const Vector& sigma1,
374  const Vector& sigma2,
375  const Vector& lc1,
376  const Vector& lc2,
377  const Numeric& co,
378  const String& fname,
379  const Verbosity&);
380 
381 template void covmat1D(Sparse& block,
382  const Vector& grid1,
383  const Vector& grid2,
384  const Vector& sigma1,
385  const Vector& sigma2,
386  const Vector& lc1,
387  const Vector& lc2,
388  const Numeric& co,
389  const String& fname,
390  const Verbosity&);
391 
393 // Manipulation of covmat_se and covmat_sx
395 
396 template <typename MatrixType>
398  const MatrixType& block,
399  const Index& i,
400  const Index& j,
401  const Verbosity&) {
402  Index m = block.nrows();
403  Index n = block.ncols();
404 
405  Index ii(i), jj(j);
406  if ((ii < 0) && (jj < 0)) {
407  ii = covmat_se.ndiagblocks();
408  jj = ii;
409  }
410 
411  if (j < i) {
412  throw std::runtime_error(
413  "The block must be on or above the diagonal, "
414  " i.e. *i* <= *j*.");
415  }
416 
417  Index ndiagblocks = covmat_se.ndiagblocks();
418 
419  if ((ii >= ndiagblocks)) {
420  if (ii > ndiagblocks) {
421  if (jj > ii) {
422  throw std::runtime_error(
423  "Off-diagonal block can only be added to rows that already "
424  "have a block on the diagonal.");
425  } else {
426  throw std::runtime_error(
427  "Diagonal block must be added row-by-row starting in the "
428  " upper left of the matrix.");
429  }
430  }
431  }
432 
433  if (covmat_se.has_block(ii, jj)) {
434  throw std::runtime_error("Block already present in covariance matrix.");
435  }
436 
437  if (ii == jj) {
438  if (m != n) {
439  throw std::runtime_error("Diagonal blocks must be square.");
440  }
441  Index start = covmat_se.nrows();
442  Range range(start, m);
443  std::shared_ptr<MatrixType> mat = std::make_shared<MatrixType>(block);
444  covmat_se.add_correlation(Block(range, range, std::make_pair(ii, ii), mat));
445 
446  } else {
447  const Block* b = covmat_se.get_block(ii, ii);
448  if (!b) {
449  throw std::runtime_error(
450  "Trying to add an off-diagonal block that"
451  " lacks corresponding diagonal block in the "
452  " same row.");
453  }
454  Range row_range = b->get_row_range();
455 
456  b = covmat_se.get_block(jj, jj);
457  if (!b) {
458  throw std::runtime_error(
459  "Trying to add an off-diagonal block that"
460  " lacks corresponding diagonal block in the "
461  " same column.");
462  }
463  Range column_range = b->get_column_range();
464 
465  if ((row_range.get_extent() != m) || (column_range.get_extent() != n)) {
466  throw std::runtime_error(
467  "The off-diagonal block is inconsistent "
468  "with the corresponding diagonal blocks.");
469  }
470 
471  std::shared_ptr<MatrixType> mat = std::make_shared<MatrixType>(block);
472  covmat_se.add_correlation(
473  Block(row_range, column_range, std::make_pair(ii, jj), mat));
474  }
475 }
476 
477 template void covmat_seAddBlock(CovarianceMatrix& covmat_se,
478  const Matrix& block,
479  const Index& i,
480  const Index& j,
481  const Verbosity&);
482 
483 template void covmat_seAddBlock(CovarianceMatrix& covmat_se,
484  const Sparse& block,
485  const Index& i,
486  const Index& j,
487  const Verbosity&);
488 
489 template <typename MatrixType>
491  const MatrixType& inv_block,
492  const Index& i,
493  const Index& j,
494  const Verbosity&) {
495  Index ii(i), jj(j);
496  if ((ii < 0) && (jj < 0)) {
497  ii = covmat_se.ndiagblocks() - 1;
498  jj = ii;
499  }
500 
501  Index m = inv_block.nrows();
502  Index n = inv_block.ncols();
503 
504  const Block* block = covmat_se.get_block(ii, jj);
505 
506  if (!block) {
507  throw std::runtime_error(
508  "Cannot add inverse block to the covariance "
509  " without corresponding non-inverse block.");
510  }
511 
512  if ((m != inv_block.nrows()) || (n != inv_block.ncols())) {
513  throw std::runtime_error(
514  "Dimensions of block are inconsistent with "
515  " non-inverse block.");
516  }
517 
518  Range row_range = block->get_row_range();
519  Range column_range = block->get_column_range();
520 
521  std::shared_ptr<MatrixType> mat = std::make_shared<MatrixType>(inv_block);
522  covmat_se.add_correlation_inverse(
523  Block(row_range, column_range, std::make_pair(ii, jj), mat));
524 }
525 
526 template void covmat_seAddInverseBlock(CovarianceMatrix& covmat_se,
527  const Matrix& block,
528  const Index& i,
529  const Index& j,
530  const Verbosity&);
531 
532 template void covmat_seAddInverseBlock(CovarianceMatrix& covmat_se,
533  const Sparse& block,
534  const Index& i,
535  const Index& j,
536  const Verbosity&);
537 
538 template <typename MatrixType>
539 void setCovarianceMatrix(CovarianceMatrix& covmat, const MatrixType& block) {
540  Index m = block.nrows();
541  Index n = block.ncols();
542 
543  if (n != m) {
544  throw std::runtime_error("Covariance matrix must be sqare!");
545  }
546 
547  covmat = CovarianceMatrix();
548  IndexPair indices = std::make_pair(0, 0);
549  Range range = Range(0, n);
550  std::shared_ptr<MatrixType> mat_ptr = std::make_shared<MatrixType>(block);
551  covmat.add_correlation(Block(range, range, indices, mat_ptr));
552 }
553 
554 template <>
556  const CovarianceMatrix& block) {
557  covmat = block;
558 }
559 
560 template <typename MatrixType>
562  const MatrixType& block,
563  const Verbosity& /*v*/) {
564  setCovarianceMatrix(covmat, block);
565 }
566 
567 template void covmat_seSet(CovarianceMatrix& covmat,
568  const CovarianceMatrix& block,
569  const Verbosity& /*v*/);
570 
571 template void covmat_seSet(CovarianceMatrix& covmat,
572  const Matrix& block,
573  const Verbosity& /*v*/);
574 
575 template void covmat_seSet(CovarianceMatrix& covmat,
576  const Sparse& block,
577  const Verbosity& /*v*/);
578 
579 template <typename MatrixType>
581  const MatrixType& block,
582  const Verbosity& /*v*/) {
583  setCovarianceMatrix(covmat, block);
584 }
585 
586 template void covmat_sxSet(CovarianceMatrix& covmat,
587  const CovarianceMatrix& block,
588  const Verbosity& /*v*/);
589 
590 template void covmat_sxSet(CovarianceMatrix& covmat,
591  const Matrix& block,
592  const Verbosity& /*v*/);
593 
594 template void covmat_sxSet(CovarianceMatrix& covmat,
595  const Sparse& block,
596  const Verbosity& /*v*/);
597 
598 template <typename MatrixType>
600  const ArrayOfRetrievalQuantity& jq,
601  const MatrixType& block,
602  const Index& i,
603  const Index& j,
604  const Verbosity& /*v*/) {
605  Index ii(i), jj(j);
606  if ((ii < 0) && (jj < 0)) {
607  ii = covmat_sx.ndiagblocks();
608  jj = ii;
609  if ((ii >= jq.nelem()) || (jj >= jq.nelem())) {
610  throw runtime_error(
611  "*covmat_sx* already contains more or as many diagonal"
612  " blocks as there are retrieval quantities.");
613  }
614  } else if ((ii >= jq.nelem()) || (jj >= jq.nelem())) {
615  throw runtime_error(
616  "The block indices must either be both -1 (default) or\n"
617  "non-negative and smaller than the number of retrieval \n"
618  "quantities.");
619  } else if (ii > jj) {
620  throw runtime_error(
621  "Only blocks above the diagonal can be set, hence"
622  "*i* must be less than or equal to *j*.");
623  }
624 
625  Index m = block.nrows();
626  Index n = block.ncols();
627 
628  Index jq_m = jq[ii].nelem();
629  if (jq[ii].HasAffine()) {
630  jq_m = jq[ii].TransformationMatrix().ncols();
631  }
632  Index jq_n = jq[jj].nelem();
633  if (jq[jj].HasAffine()) {
634  jq_n = jq[jj].TransformationMatrix().ncols();
635  }
636 
637  if ((m != jq_m) || (n != jq_n)) {
638  ostringstream os;
639  os << "The covariance block for retrieval quantities " << ii << " and "
640  << jj << " was expected to have dimensions [" << jq_m << ", " << jq_n
641  << "] but found to have dimensions "
642  << "[" << m << ", " << n << "].";
643 
644  throw runtime_error(os.str());
645  }
646 
648  bool any_affine;
649  jac_ranges_indices(ji, any_affine, jq);
650  Index row_start = ji[ii][0];
651  Index row_extent = ji[ii][1] - ji[ii][0] + 1;
652  Range row_range(row_start, row_extent);
653 
654  Index col_start = ji[jj][0];
655  Index col_extent = ji[jj][1] - ji[jj][0] + 1;
656  Range col_range(col_start, col_extent);
657 
658  std::shared_ptr<MatrixType> mat = make_shared<MatrixType>(block);
659  covmat_sx.add_correlation(
660  Block(row_range, col_range, std::make_pair(ii, jj), mat));
661 }
662 
663 template void covmat_sxAddBlock(CovarianceMatrix&,
665  const Matrix&,
666  const Index&,
667  const Index&,
668  const Verbosity&);
669 
670 template void covmat_sxAddBlock(CovarianceMatrix&,
672  const Sparse&,
673  const Index&,
674  const Index&,
675  const Verbosity&);
676 
677 template <typename MatrixType>
679  const ArrayOfRetrievalQuantity& jq,
680  const MatrixType& block_inv,
681  const Index& i,
682  const Index& j,
683  const Verbosity& /*v*/) {
684  Index ii(i), jj(j);
685  if ((ii < 0) && (jj < 0)) {
686  ii = jq.size() - 1;
687  jj = jq.size() - 1;
688  } else if ((ii >= jq.nelem()) || (jj >= jq.nelem())) {
689  throw runtime_error(
690  "The block indices must either be both -1 (default) or\n"
691  "non-negative and smaller than the number of retrieval \n"
692  "quantities.");
693  } else if (ii > jj) {
694  throw runtime_error(
695  "Only blocks above the diagonal can be set, hence"
696  "*i* must be less than or equal to *j*.");
697  }
698 
699  Index m = block_inv.nrows();
700  Index n = block_inv.ncols();
701 
702  Index jq_m = jq[ii].nelem();
703  if (jq[ii].HasAffine()) {
704  jq_m = jq[ii].TransformationMatrix().ncols();
705  }
706 
707  Index jq_n = jq[jj].nelem();
708  if (jq[jj].HasAffine()) {
709  jq_n = jq[jj].TransformationMatrix().ncols();
710  }
711 
712  if (!((m == jq_m) && (n == jq_n))) {
713  ostringstream os;
714  os << "The dimensions of the covariance block ( " << m;
715  os << " x " << n << " )"
716  << " with the dimensionality of ";
717  os << " retrieval quantity " << ii << " and " << jj << ", respectively.";
718 
719  throw runtime_error(os.str());
720  }
721 
722  if (!covmat_sx.has_block(ii, jj)) {
723  throw runtime_error(
724  "To add the inverse of a block the non-inverse"
725  " block must be added first.");
726  }
727 
729  bool any_affine;
730  jac_ranges_indices(ji, any_affine, jq);
731  Index row_start = ji[ii][0];
732  Index row_extent = ji[ii][1] - ji[ii][0] + 1;
733  Range row_range(row_start, row_extent);
734 
735  Index col_start = ji[jj][0];
736  Index col_extent = ji[jj][1] - ji[jj][0] + 1;
737  Range col_range(col_start, col_extent);
738 
739  std::shared_ptr<MatrixType> mat = make_shared<MatrixType>(block_inv);
740  covmat_sx.add_correlation_inverse(
741  Block(row_range, col_range, std::make_pair(ii, jj), mat));
742 }
743 
746  const Matrix&,
747  const Index&,
748  const Index&,
749  const Verbosity&);
750 
753  const Sparse&,
754  const Index&,
755  const Index&,
756  const Verbosity&);
757 
759  const CovarianceMatrix& covmat_sx,
760  const Verbosity&) {
761  x_norm = covmat_sx.diagonal();
762  for (Index i = 0; i < x_norm.nelem(); ++i) {
763  x_norm[i] = sqrt(x_norm[i]);
764  }
765 }
766 
768 // retrievalAdd Functions
770 
772  CovarianceMatrix& covmat_sx,
773  ArrayOfRetrievalQuantity& jacobian_quantities,
774  Agenda& jacobian_agenda,
775  const Index& atmosphere_dim,
776  const Sparse& covmat_block,
777  const Sparse& covmat_inv_block,
778  const Vector& p_grid,
779  const Vector& lat_grid,
780  const Vector& lon_grid,
781  const Vector& rq_p_grid,
782  const Vector& rq_lat_grid,
783  const Vector& rq_lon_grid,
784  const String& species,
785  const String& mode,
786  const Index& for_species_tag,
787  const Verbosity& verbosity) {
789  jacobian_quantities,
790  jacobian_agenda,
791  atmosphere_dim,
792  p_grid,
793  lat_grid,
794  lon_grid,
795  rq_p_grid,
796  rq_lat_grid,
797  rq_lon_grid,
798  species,
799  mode,
800  for_species_tag,
801  verbosity);
802  check_and_add_block(covmat_sx,
803  jacobian_quantities.back(),
804  jacobian_quantities.nelem() - 1,
805  atmosphere_dim,
806  covmat_block,
807  covmat_inv_block);
808 }
809 
811  CovarianceMatrix& covmat_sx,
812  ArrayOfRetrievalQuantity& jacobian_quantities,
813  Agenda& jacobian_agenda,
814  const Sparse& covmat_block,
815  const Sparse& covmat_inv_block,
816  const Vector& f_grid,
817  const Numeric& df,
818  const Verbosity& verbosity) {
820  ws, jacobian_quantities, jacobian_agenda, f_grid, df, verbosity);
821  check_and_add_block(covmat_sx,
822  jacobian_quantities.back(),
823  jacobian_quantities.nelem() - 1,
824  1,
825  covmat_block,
826  covmat_inv_block);
827 }
828 
830  CovarianceMatrix& covmat_sx,
831  ArrayOfRetrievalQuantity& jacobian_quantities,
832  Agenda& jacobian_agenda,
833  const Vector& f_grid,
834  const Sparse& covmat_block,
835  const Sparse& covmat_inv_block,
836  const Numeric& df,
837  const Verbosity& verbosity) {
839  ws, jacobian_quantities, jacobian_agenda, f_grid, df, verbosity);
840  check_and_add_block(covmat_sx,
841  jacobian_quantities.back(),
842  jacobian_quantities.nelem() - 1,
843  1,
844  covmat_block,
845  covmat_inv_block);
846 }
847 
849  CovarianceMatrix& covmat_sx,
850  ArrayOfRetrievalQuantity& jacobian_quantities,
851  Agenda& jacobian_agenda,
852  const QuantumIdentifier& catalog_identity,
853  const String& catalog_parameter,
854  const Numeric& var,
855  const Verbosity& verbosity) {
857  jacobian_quantities,
858  jacobian_agenda,
859  catalog_identity,
860  catalog_parameter,
861  verbosity);
862  add_scalar_variance(covmat_sx, jacobian_quantities, var);
863 }
864 
866  Workspace& ws,
867  CovarianceMatrix& covmat_sx,
868  ArrayOfRetrievalQuantity& jacobian_quantities,
869  Agenda& jacobian_agenda,
870  const Sparse& covmat_block,
871  const Sparse& covmat_inv_block,
872  const ArrayOfQuantumIdentifier& catalog_identities,
873  const ArrayOfString& catalog_parameters,
874  const Verbosity& verbosity) {
876  jacobian_quantities,
877  jacobian_agenda,
878  catalog_identities,
879  catalog_parameters,
880  verbosity);
881  check_and_add_block(covmat_sx,
882  jacobian_quantities.back(),
883  jacobian_quantities.nelem() - 1,
884  0,
885  covmat_block,
886  covmat_inv_block);
887 }
888 
890  CovarianceMatrix& covmat_sx,
891  ArrayOfRetrievalQuantity& jacobian_quantities,
892  Agenda& jacobian_agenda,
893  const Index& atmosphere_dim,
894  const Sparse& covmat_block,
895  const Sparse& covmat_inv_block,
896  const Vector& p_grid,
897  const Vector& lat_grid,
898  const Vector& lon_grid,
899  const Vector& rq_p_grid,
900  const Vector& rq_lat_grid,
901  const Vector& rq_lon_grid,
902  const String& component,
903  const Numeric& dB,
904  const Verbosity& verbosity) {
906  jacobian_quantities,
907  jacobian_agenda,
908  atmosphere_dim,
909  p_grid,
910  lat_grid,
911  lon_grid,
912  rq_p_grid,
913  rq_lat_grid,
914  rq_lon_grid,
915  component,
916  dB,
917  verbosity);
918  check_and_add_block(covmat_sx,
919  jacobian_quantities.back(),
920  jacobian_quantities.nelem() - 1,
921  atmosphere_dim,
922  covmat_block,
923  covmat_inv_block);
924 }
925 
927  CovarianceMatrix& covmat_sx,
928  ArrayOfRetrievalQuantity& jacobian_quantities,
929  Agenda& jacobian_agenda,
930  const Sparse& covmat_block,
931  const Sparse& covmat_inv_block,
932  const Matrix& sensor_pos,
933  const Vector& sensor_time,
934  const Index& poly_order,
935  const String& calcmode,
936  const Numeric& dza,
937  const Verbosity& verbosity) {
939  jacobian_quantities,
940  jacobian_agenda,
941  sensor_pos,
942  sensor_time,
943  poly_order,
944  calcmode,
945  dza,
946  verbosity);
947  check_and_add_block(covmat_sx,
948  jacobian_quantities.back(),
949  jacobian_quantities.nelem() - 1,
950  1,
951  covmat_block,
952  covmat_inv_block);
953 }
954 
956  CovarianceMatrix& covmat_sx,
957  ArrayOfRetrievalQuantity& jacobian_quantities,
958  Agenda& jacobian_agenda,
959  const Sparse& covmat_block,
960  const Sparse& covmat_inv_block,
961  const ArrayOfIndex& sensor_response_pol_grid,
962  const Matrix& sensor_response_dlos_grid,
963  const Matrix& sensor_pos,
964  const Index& poly_order,
965  const Index& no_pol_variation,
966  const Index& no_los_variation,
967  const Index& no_mblock_variation,
968  const Verbosity& verbosity) {
969  size_t jq_start = jacobian_quantities.size();
971  jacobian_quantities,
972  jacobian_agenda,
973  sensor_response_pol_grid,
974  sensor_response_dlos_grid,
975  sensor_pos,
976  poly_order,
977  no_pol_variation,
978  no_los_variation,
979  no_mblock_variation,
980  verbosity);
981  for (Index i = 0; i <= poly_order; ++i) {
982  check_and_add_block(covmat_sx,
983  jacobian_quantities[jq_start + i],
984  jq_start + i,
985  4,
986  covmat_block,
987  covmat_inv_block);
988  }
989 }
990 
992  CovarianceMatrix& covmat_sx,
993  ArrayOfRetrievalQuantity& jacobian_quantities,
994  Agenda& jacobian_agenda,
995  const Index& atmosphere_dim,
996  const Sparse& covmat_block,
997  const Sparse& covmat_inv_block,
998  const Vector& p_grid,
999  const Vector& lat_grid,
1000  const Vector& lon_grid,
1001  const Vector& rq_p_grid,
1002  const Vector& rq_lat_grid,
1003  const Vector& rq_lon_grid,
1004  const String& species,
1005  const String& quantity,
1006  const Verbosity& verbosity) {
1008  jacobian_quantities,
1009  jacobian_agenda,
1010  atmosphere_dim,
1011  p_grid,
1012  lat_grid,
1013  lon_grid,
1014  rq_p_grid,
1015  rq_lat_grid,
1016  rq_lon_grid,
1017  species,
1018  quantity,
1019  verbosity);
1020  check_and_add_block(covmat_sx,
1021  jacobian_quantities.back(),
1022  jacobian_quantities.nelem() - 1,
1023  atmosphere_dim,
1024  covmat_block,
1025  covmat_inv_block);
1026 }
1027 
1029  CovarianceMatrix& covmat_sx,
1030  ArrayOfRetrievalQuantity& jacobian_quantities,
1031  Agenda& jacobian_agenda,
1032  const Sparse& covmat_block,
1033  const Sparse& covmat_inv_block,
1034  const ArrayOfIndex& sensor_response_pol_grid,
1035  const Matrix& sensor_response_dlos_grid,
1036  const Matrix& sensor_pos,
1037  const Vector& period_lengths,
1038  const Index& no_pol_variation,
1039  const Index& no_los_variation,
1040  const Index& no_mblock_variation,
1041  const Verbosity& verbosity) {
1042  size_t jq_start = jacobian_quantities.size();
1043  jacobianAddSinefit(ws,
1044  jacobian_quantities,
1045  jacobian_agenda,
1046  sensor_response_pol_grid,
1047  sensor_response_dlos_grid,
1048  sensor_pos,
1049  period_lengths,
1050  no_pol_variation,
1051  no_los_variation,
1052  no_mblock_variation,
1053  verbosity);
1054  for (Index i = 0; i < period_lengths.nelem(); ++i) {
1055  check_and_add_block(covmat_sx,
1056  jacobian_quantities[jq_start + i],
1057  jq_start + i,
1058  4,
1059  covmat_block,
1060  covmat_inv_block);
1061  }
1062 }
1063 
1065  CovarianceMatrix& covmat_sx,
1066  ArrayOfRetrievalQuantity& jacobian_quantities,
1067  Agenda& jacobian_agenda,
1068  const Index& atmosphere_dim,
1069  const Sparse& covmat_block,
1070  const Sparse& covmat_inv_block,
1071  const Vector& p_grid,
1072  const Vector& lat_grid,
1073  const Vector& lon_grid,
1074  const Vector& rq_p_grid,
1075  const Vector& rq_lat_grid,
1076  const Vector& rq_lon_grid,
1077  const String& species,
1078  const Verbosity& verbosity) {
1080  jacobian_quantities,
1081  jacobian_agenda,
1082  atmosphere_dim,
1083  p_grid,
1084  lat_grid,
1085  lon_grid,
1086  rq_p_grid,
1087  rq_lat_grid,
1088  rq_lon_grid,
1089  species,
1090  verbosity);
1091  check_and_add_block(covmat_sx,
1092  jacobian_quantities.back(),
1093  jacobian_quantities.nelem() - 1,
1094  atmosphere_dim,
1095  covmat_block,
1096  covmat_inv_block);
1097 }
1098 
1100  CovarianceMatrix& covmat_sx,
1101  ArrayOfRetrievalQuantity& jacobian_quantities,
1102  Agenda& jacobian_agenda,
1103  const Index& atmosphere_dim,
1104  const Sparse& covmat_block,
1105  const Sparse& covmat_inv_block,
1106  const Vector& p_grid,
1107  const Vector& lat_grid,
1108  const Vector& lon_grid,
1109  const Vector& rq_p_grid,
1110  const Vector& rq_lat_grid,
1111  const Vector& rq_lon_grid,
1112  const String& component,
1113  const Numeric& dfrequency,
1114  const Verbosity& verbosity) {
1115  jacobianAddWind(ws,
1116  jacobian_quantities,
1117  jacobian_agenda,
1118  atmosphere_dim,
1119  p_grid,
1120  lat_grid,
1121  lon_grid,
1122  rq_p_grid,
1123  rq_lat_grid,
1124  rq_lon_grid,
1125  component,
1126  dfrequency,
1127  verbosity);
1128  check_and_add_block(covmat_sx,
1129  jacobian_quantities.back(),
1130  jacobian_quantities.nelem() - 1,
1131  atmosphere_dim,
1132  covmat_block,
1133  covmat_inv_block);
1134 }
1135 
1137  CovarianceMatrix& covmat_sx,
1138  ArrayOfRetrievalQuantity& jacobian_quantities,
1139  Agenda& jacobian_agenda,
1140  const Index& atmosphere_dim,
1141  const Sparse& covmat_block,
1142  const Sparse& covmat_inv_block,
1143  const Vector& p_grid,
1144  const Vector& lat_grid,
1145  const Vector& lon_grid,
1146  const Vector& rq_p_grid,
1147  const Vector& rq_lat_grid,
1148  const Vector& rq_lon_grid,
1149  const String& hse,
1150  const Verbosity& verbosity) {
1152  jacobian_quantities,
1153  jacobian_agenda,
1154  atmosphere_dim,
1155  p_grid,
1156  lat_grid,
1157  lon_grid,
1158  rq_p_grid,
1159  rq_lat_grid,
1160  rq_lon_grid,
1161  hse,
1162  verbosity);
1163  check_and_add_block(covmat_sx,
1164  jacobian_quantities.back(),
1165  jacobian_quantities.nelem() - 1,
1166  atmosphere_dim,
1167  covmat_block,
1168  covmat_inv_block);
1169 }
1170 
1172  CovarianceMatrix& covmat_sx,
1173  ArrayOfRetrievalQuantity& jacobian_quantities,
1174  Agenda& jacobian_agenda,
1175  const Sparse& covmat_block,
1176  const Sparse& covmat_inv_block,
1177  const Index& atmosphere_dim,
1178  const Vector& lat_grid,
1179  const Vector& lon_grid,
1180  const Vector& rq_lat_grid,
1181  const Vector& rq_lon_grid,
1182  const String& quantity,
1183  const Verbosity& verbosity) {
1185  jacobian_quantities,
1186  jacobian_agenda,
1187  atmosphere_dim,
1188  lat_grid,
1189  lon_grid,
1190  rq_lat_grid,
1191  rq_lon_grid,
1192  quantity,
1193  verbosity);
1194 
1195  check_and_add_block(covmat_sx,
1196  jacobian_quantities.back(),
1197  jacobian_quantities.nelem() - 1,
1198  atmosphere_dim - 1,
1199  covmat_block,
1200  covmat_inv_block);
1201 }
1202 
1203 /* Workspace method: Doxygen documentation will be auto-generated */
1205  Index& jacobian_do,
1206  Agenda& jacobian_agenda,
1207  Index& retrieval_checked,
1208  const CovarianceMatrix& covmat_sx,
1209  const ArrayOfRetrievalQuantity& jacobian_quantities,
1210  const Verbosity& verbosity) {
1211  jacobianClose(
1212  ws, jacobian_do, jacobian_agenda, jacobian_quantities, verbosity);
1213 
1214  ArrayOfArrayOfIndex ji_t;
1215  bool any_affine;
1216  jac_ranges_indices(ji_t, any_affine, jacobian_quantities);
1217  if (!covmat_sx.has_diagonal_blocks(ji_t)) {
1218  std::ostringstream os;
1219  os << "*covmat_sx* does not contain a diagonal block for each retrieval"
1220  " quantity in the Jacobian.\n";
1221  os << " Fails test (!covmat_sx.has_diagonal_blocks(ji_t)) for ji_t " << ji_t
1222  << "\n";
1223  throw runtime_error(os.str());
1224  }
1225  if (!covmat_sx.is_consistent(ji_t)) {
1226  throw runtime_error(
1227  "The blocks in *covmat_sx* are not consistent with the retrieval"
1228  " quantities in the Jacobian.");
1229  }
1230  retrieval_checked = true;
1231 }
1232 
1233 /* Workspace method: Doxygen documentation will be auto-generated */
1235  CovarianceMatrix& covmat_sx,
1236  Sparse& covmat_block,
1237  Sparse& covmat_inv_block,
1238  ArrayOfRetrievalQuantity& jacobian_quantities,
1239  Agenda& jacobian_agenda,
1240  const Index& initialize_jacobian,
1241  const Verbosity& verbosity) {
1242  if (initialize_jacobian == 1) {
1243  jacobianInit(jacobian_quantities, jacobian_agenda, verbosity);
1244  }
1245 
1246  covmat_block = Sparse();
1247  covmat_inv_block = Sparse();
1248  covmat_sx = CovarianceMatrix();
1249  covmat_se = CovarianceMatrix();
1250 }
1251 
1252 void retrievalErrorsExtract(Vector& retrieval_eo,
1253  Vector& retrieval_ss,
1254  const Matrix& covmat_so,
1255  const Matrix& covmat_ss,
1256  const Verbosity& /*v*/) {
1257  Index n_so = covmat_so.nrows();
1258  Index n_ss = covmat_ss.nrows();
1259 
1260  retrieval_eo.resize(n_so);
1261  for (Index i = 0; i < n_so; ++i) {
1262  retrieval_eo[i] = sqrt(covmat_so(i, i));
1263  }
1264 
1265  retrieval_ss.resize(n_ss);
1266  for (Index i = 0; i < n_ss; ++i) {
1267  retrieval_ss[i] = sqrt(covmat_ss(i, i));
1268  }
1269 }
std::pair< Index, Index > IndexPair
void retrievalAddPolyfit(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Sparse &covmat_block, const Sparse &covmat_inv_block, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Matrix &sensor_pos, const Index &poly_order, const Index &no_pol_variation, const Index &no_los_variation, const Index &no_mblock_variation, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddPolyfit.
Definition: m_retrieval.cc:955
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
void covmatDiagonal(MatrixType &block, MatrixType &block_inv, const Vector &vars, const Verbosity &)
Definition: m_retrieval.cc:135
void jacobianAddAbsSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &mode, const Index &for_species_tag, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddAbsSpecies.
Definition: m_jacobian.cc:159
void covmat_seAddBlock(CovarianceMatrix &covmat_se, const MatrixType &block, const Index &i, const Index &j, const Verbosity &)
Definition: m_retrieval.cc:397
void jacobianClose(Workspace &ws, Index &jacobian_do, Agenda &jacobian_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianClose.
Definition: m_jacobian.cc:122
void var(VectorView var, const Vector &y, const ArrayOfVector &ys, const Index start=0, const Index end=-1)
Compute the variance of the ranged ys.
Definition: raw.cc:49
const ArrayOfVector & Grids() const
Returns the grids of the retrieval.
Definition: jacobian.h:255
Range get_column_range() const
void covmat1DMarkov(MatrixType &block, MatrixType &block_inv, const Vector &grid, const Vector &sigma, const Numeric &lc, const Numeric &, const Verbosity &)
Definition: m_retrieval.cc:169
void retrievalAddMagField(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Sparse &covmat_block, const Sparse &covmat_inv_block, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &component, const Numeric &dB, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddMagField.
Definition: m_retrieval.cc:889
void covmat_seAddInverseBlock(CovarianceMatrix &covmat_se, const MatrixType &inv_block, const Index &i, const Index &j, const Verbosity &)
Definition: m_retrieval.cc:490
bool has_block(Index i, Index j)
Check if the block with indices (i, j) is contained in the covariance matrix.
void jacobianAddPointingZa(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Matrix &sensor_pos, const Vector &sensor_time, const Index &poly_order, const String &calcmode, const Numeric &dza, const Verbosity &)
WORKSPACE METHOD: jacobianAddPointingZa.
Definition: m_jacobian.cc:600
The Agenda class.
Definition: agenda_class.h:44
Index nelem() const
Number of elements.
Definition: array.h:195
invlib::Matrix< ArtsCovarianceMatrixWrapper > CovarianceMatrix
invlib wrapper type for ARTS the ARTS covariance class.
Definition: oem.h:38
void jacobianAddFreqShift(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Vector &f_grid, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: jacobianAddFreqShift.
Definition: m_jacobian.cc:271
void retrievalAddTemperature(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Sparse &covmat_block, const Sparse &covmat_inv_block, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &hse, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddTemperature.
void covmat_sxSet(CovarianceMatrix &covmat, const MatrixType &block, const Verbosity &)
Definition: m_retrieval.cc:580
void retrievalAddSpecialSpecies(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Sparse &covmat_block, const Sparse &covmat_inv_block, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddSpecialSpecies.
void retrievalDefInit(CovarianceMatrix &covmat_se, CovarianceMatrix &covmat_sx, Sparse &covmat_block, Sparse &covmat_inv_block, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Index &initialize_jacobian, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalDefInit.
Routines for setting up the jacobian.
The Vector class.
Definition: matpackI.h:860
#define abs(x)
void jacobianAddMagField(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &component, const Numeric &dB, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddMagField.
Definition: m_jacobian.cc:1557
The Sparse class.
Definition: matpackII.h:60
void jacobianAddFreqStretch(Workspace &ws, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Vector &f_grid, const Numeric &df, const Verbosity &)
WORKSPACE METHOD: jacobianAddFreqStretch.
Definition: m_jacobian.cc:425
The range class.
Definition: matpackI.h:160
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixType
Definition: matpackI.h:109
bool empty() const
Returns true if variable size is zero.
Definition: matpackI.cc:49
Index ncols() const
Returns the number of columns.
Definition: matpackII.cc:69
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
Definition: oem.h:32
void PFromZSimple(Vector &p_grid, const Vector &z_grid, const Verbosity &)
WORKSPACE METHOD: PFromZSimple.
Definition: m_retrieval.cc:88
void add_scalar_variance(CovarianceMatrix &covmat, ArrayOfRetrievalQuantity &jacobian_quantities, Numeric var)
Definition: m_retrieval.cc:54
void jacobianAddBasicCatalogParameter(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const QuantumIdentifier &catalog_identity, const String &catalog_parameter, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddBasicCatalogParameter.
Definition: m_jacobian.cc:1719
void retrievalAddSinefit(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Sparse &covmat_block, const Sparse &covmat_inv_block, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Matrix &sensor_pos, const Vector &period_lengths, const Index &no_pol_variation, const Index &no_los_variation, const Index &no_mblock_variation, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddSinefit.
bool empty() const
Returns true if variable size is zero.
Definition: matpackII.cc:61
void jacobianAddScatSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &quantity, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddScatSpecies.
Definition: m_jacobian.cc:1099
void jacobianAddTemperature(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &hse, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddTemperature.
Definition: m_jacobian.cc:1397
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
void jacobianAddWind(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &component, const Numeric &dfrequency, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddWind.
Definition: m_jacobian.cc:1476
void retrievalAddWind(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Sparse &covmat_block, const Sparse &covmat_inv_block, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &component, const Numeric &dfrequency, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddWind.
void insert_elements(MatrixType &matrix, Index, Index, const ArrayOfIndex &row_indices, const ArrayOfIndex &column_indices, const Vector &elements)
Definition: m_retrieval.cc:107
void jacobianAddSurfaceQuantity(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &quantity, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddSurfaceQuantity.
Definition: m_jacobian.cc:1338
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition: jacobian.h:120
void retrievalAddCatalogParameters(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Sparse &covmat_block, const Sparse &covmat_inv_block, const ArrayOfQuantumIdentifier &catalog_identities, const ArrayOfString &catalog_parameters, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddCatalogParameters.
Definition: m_retrieval.cc:865
bool is_consistent(const ArrayOfArrayOfIndex &jis) const
Checks that the dimensions of the blocks in the covariance matrix are consistent with ranges occupied...
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
void retrievalAddPointingZa(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Sparse &covmat_block, const Sparse &covmat_inv_block, const Matrix &sensor_pos, const Vector &sensor_time, const Index &poly_order, const String &calcmode, const Numeric &dza, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddPointingZa.
Definition: m_retrieval.cc:926
Header files of CovarianceMatrix class.
void retrievalAddFreqStretch(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Vector &f_grid, const Sparse &covmat_block, const Sparse &covmat_inv_block, const Numeric &df, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddFreqStretch.
Definition: m_retrieval.cc:829
void retrievalErrorsExtract(Vector &retrieval_eo, Vector &retrieval_ss, const Matrix &covmat_so, const Matrix &covmat_ss, const Verbosity &)
WORKSPACE METHOD: retrievalErrorsExtract.
_CS_string_type str() const
Definition: sstream.h:491
void ZFromPSimple(Vector &z_grid, const Vector &p_grid, const Verbosity &)
WORKSPACE METHOD: ZFromPSimple.
Definition: m_retrieval.cc:74
void covmat1D(MatrixType &block, const Vector &grid1, const Vector &grid2, const Vector &sigma1, const Vector &sigma2, const Vector &lc1, const Vector &lc2, const Numeric &co, const String &fname, const Verbosity &)
Definition: m_retrieval.cc:251
constexpr Index get_extent() const
Returns the extent of the range.
Definition: matpackI.h:329
void jacobianAddPolyfit(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Matrix &sensor_pos, const Index &poly_order, const Index &no_pol_variation, const Index &no_los_variation, const Index &no_mblock_variation, const Verbosity &)
WORKSPACE METHOD: jacobianAddPolyfit.
Definition: m_jacobian.cc:937
Index nrows() const
Returns the number of rows.
Definition: matpackII.cc:66
Class to identify and match lines by their quantum numbers.
Definition: quantum.h:390
void retrievalAddSurfaceQuantity(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Sparse &covmat_block, const Sparse &covmat_inv_block, const Index &atmosphere_dim, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &quantity, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddSurfaceQuantity.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
The Matrix class.
Definition: matpackI.h:1193
Range get_row_range() const
void covmat_sxAddInverseBlock(CovarianceMatrix &covmat_sx, const ArrayOfRetrievalQuantity &jq, const MatrixType &block_inv, const Index &i, const Index &j, const Verbosity &)
Definition: m_retrieval.cc:678
Vector diagonal() const
Diagonal elements as vector.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
basic_ostringstream< char, string_char_traits< char >, alloc > ostringstream
Definition: sstream.h:204
void retrievalAddFreqShift(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Sparse &covmat_block, const Sparse &covmat_inv_block, const Vector &f_grid, const Numeric &df, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddFreqShift.
Definition: m_retrieval.cc:810
void retrievalAddAbsSpecies(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Sparse &covmat_block, const Sparse &covmat_inv_block, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &mode, const Index &for_species_tag, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddAbsSpecies.
Definition: m_retrieval.cc:771
This can be used to make arrays out of anything.
Definition: array.h:40
void setCovarianceMatrix(CovarianceMatrix &covmat, const MatrixType &block)
Definition: m_retrieval.cc:539
void jacobianAddBasicCatalogParameters(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfQuantumIdentifier &catalog_identities, const ArrayOfString &catalog_parameters, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddBasicCatalogParameters.
Definition: m_jacobian.cc:1773
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
void retrievalAddCatalogParameter(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const QuantumIdentifier &catalog_identity, const String &catalog_parameter, const Numeric &var, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddCatalogParameter.
Definition: m_retrieval.cc:848
void check_and_add_block(CovarianceMatrix &covmat, const RetrievalQuantity &jq, const Index rq_index, const Index grid_dimensions, const Sparse &covmat_block, const Sparse &covmat_inv_block)
Definition: m_retrieval.cc:11
A constant view of a Vector.
Definition: matpackI.h:476
void retrievalDefClose(Workspace &ws, Index &jacobian_do, Agenda &jacobian_agenda, Index &retrieval_checked, const CovarianceMatrix &covmat_sx, const ArrayOfRetrievalQuantity &jacobian_quantities, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalDefClose.
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:77
Workspace class.
Definition: workspace_ng.h:40
const Block * get_block(Index i=-1, Index j=-1)
Return a pointer to the block with indices (i,j).
void covmat_sxExtractSqrtDiagonal(Vector &x_norm, const CovarianceMatrix &covmat_sx, const Verbosity &)
WORKSPACE METHOD: covmat_sxExtractSqrtDiagonal.
Definition: m_retrieval.cc:758
void add_correlation_inverse(Block c)
Add block inverse to covariance matrix.
Index ndiagblocks() const
! The number of diagonal blocks in the matrix excluding inverses.
void add_correlation(Block c)
Add block to covariance matrix.
void covmat_sxAddBlock(CovarianceMatrix &covmat_sx, const ArrayOfRetrievalQuantity &jq, const MatrixType &block, const Index &i, const Index &j, const Verbosity &)
Definition: m_retrieval.cc:599
void jacobianAddSinefit(Workspace &ws, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const ArrayOfIndex &sensor_response_pol_grid, const Matrix &sensor_response_dlos_grid, const Matrix &sensor_pos, const Vector &period_lengths, const Index &no_pol_variation, const Index &no_los_variation, const Index &no_mblock_variation, const Verbosity &)
WORKSPACE METHOD: jacobianAddSinefit.
Definition: m_jacobian.cc:1167
void retrievalAddScatSpecies(Workspace &ws, CovarianceMatrix &covmat_sx, ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Sparse &covmat_block, const Sparse &covmat_inv_block, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const String &quantity, const Verbosity &verbosity)
WORKSPACE METHOD: retrievalAddScatSpecies.
Definition: m_retrieval.cc:991
void covmat_seSet(CovarianceMatrix &covmat, const MatrixType &block, const Verbosity &)
Definition: m_retrieval.cc:561
bool has_diagonal_blocks(const ArrayOfArrayOfIndex &jis) const
Checks that the covariance matrix contains one diagonal block per retrieval quantity.
void jacobianAddSpecialSpecies(Workspace &, ArrayOfRetrievalQuantity &jq, Agenda &jacobian_agenda, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &rq_p_grid, const Vector &rq_lat_grid, const Vector &rq_lon_grid, const String &species, const Verbosity &verbosity)
WORKSPACE METHOD: jacobianAddSpecialSpecies.
Definition: m_jacobian.cc:1891
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
void jacobianInit(ArrayOfRetrievalQuantity &jacobian_quantities, Agenda &jacobian_agenda, const Verbosity &)
WORKSPACE METHOD: jacobianInit.
Definition: m_jacobian.cc:137
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
void resize(Index r, Index c)
Resize function.
Definition: matpackI.cc:1056