ARTS 2.5.10 (git: 2f1c442c)
jacobian.cc
Go to the documentation of this file.
1/* Copyright (C) 2004-2012 Mattias Ekstrom <ekstrom@rss.chalmers.se>
2
3 This program is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16 USA. */
17
26#include "jacobian.h"
27#include "arts.h"
28#include "arts_constants.h"
29#include "auto_md.h"
30#include "check_input.h"
31#include "global_data.h"
32#include "lin_alg.h"
33#include "physics_funcs.h"
34#include "rte.h"
35#include "special_interp.h"
36
38inline constexpr Numeric PI=Constant::pi;
39
40ostream& operator<<(ostream& os, const RetrievalQuantity& ot) {
41 return os << "\n Target = " << ot.Target()
42 << "\n Sub tag = " << ot.Subtag()
43 << "\n Mode = " << ot.Mode();
44}
45
47 bool& any_affine,
48 const ArrayOfRetrievalQuantity& jqs,
49 const bool& before_affine) {
50 jis.resize(jqs.nelem());
51
52 any_affine = false;
53
54 // Indices before affine transformation
55 if (before_affine) {
56 for (Index i = 0; i < jqs.nelem(); ++i) {
57 jis[i] = ArrayOfIndex(2);
58 if (i > 0) {
59 jis[i][0] = jis[i - 1][1] + 1;
60 } else {
61 jis[i][0] = 0;
62 }
63 const RetrievalQuantity& jq = jqs[i];
64 jis[i][1] = jis[i][0] + jq.nelem() - 1;
65 if (jq.HasAffine()) {
66 any_affine = true;
67 }
68 }
69 }
70
71 // After affine transformation
72 else {
73 for (Index i = 0; i < jqs.nelem(); ++i) {
74 jis[i] = ArrayOfIndex(2);
75 if (i > 0) {
76 jis[i][0] = jis[i - 1][1] + 1;
77 } else {
78 jis[i][0] = 0;
79 }
80 const RetrievalQuantity& jq = jqs[i];
81 if (jq.HasAffine()) {
82 jis[i][1] = jis[i][0] + jq.TransformationMatrix().ncols() - 1;
83 any_affine = true;
84 } else {
85 jis[i][1] = jis[i][0] + jq.nelem() - 1;
86 }
87 }
88 }
89}
90
92 const Vector x,
93 const ArrayOfRetrievalQuantity& jqs) {
94 // Range indices without affine trans
96 bool any_affine;
97 //
98 jac_ranges_indices(jis, any_affine, jqs, true);
99
100 Vector x_t(x);
101 transform_x_back(x_t, jqs, false);
102
103 // Apply functional transformations
104 for (Index i = 0; i < jqs.nelem(); ++i) {
105 const RetrievalQuantity& jq = jqs[i];
106 const String tfun = jq.TransformationFunc();
107 // Remember to add new functions also to transform_jacobian and transform_x_back
108 if (tfun == "") {
109 // Nothing to do
110 } else if (tfun == "log") {
111 for (Index c = jis[i][0]; c <= jis[i][1]; ++c) {
112 jacobian(joker, c) *= exp(x_t[c]);
113 }
114 } else if (tfun == "log10") {
115 for (Index c = jis[i][0]; c <= jis[i][1]; ++c) {
116 jacobian(joker, c) *= NAT_LOG_TEN * pow(10.0, x_t[c]);
117 }
118 } else if (tfun == "atanh") {
119 const Vector& pars = jq.TFuncParameters();
120 for (Index c = jis[i][0]; c <= jis[i][1]; ++c) {
121 jacobian(joker, c) *=
122 2 * (pars[1] - pars[0]) / pow(exp(-x_t[c]) + exp(x_t[c]), 2.0);
123 }
124 } else {
125 ARTS_ASSERT(0);
126 }
127 }
128
129 // Apply affine transformations
130 if (any_affine) {
132 jac_ranges_indices(jis_t, any_affine, jqs);
133
134 Matrix jacobian_t(jacobian.nrows(), jis_t.back()[1] + 1);
135
136 for (Index i = 0; i < jqs.nelem(); ++i) {
137 const RetrievalQuantity& jq = jqs[i];
138 Index col_start = jis[i][0];
139 Index col_extent = jis[i][1] - jis[i][0] + 1;
140 Range col_range(col_start, col_extent);
141 Index col_start_t = jis_t[i][0];
142 Index col_extent_t = jis_t[i][1] - jis_t[i][0] + 1;
143 Range col_range_t(col_start_t, col_extent_t);
144 if (jq.HasAffine()) {
145 mult(jacobian_t(joker, col_range_t),
146 jacobian(joker, col_range),
148 } else {
149 jacobian_t(joker, col_range_t) = jacobian(joker, col_range);
150 }
151 }
152 swap(jacobian_t, jacobian);
153 }
154}
155
157 // Range indices without affine trans
159 bool any_affine;
160 //
161 jac_ranges_indices(jis, any_affine, jqs, true);
162
163 // Apply functional transformations
164 for (Index i = 0; i < jqs.nelem(); ++i) {
165 const RetrievalQuantity& jq = jqs[i];
166 const String tfun = jq.TransformationFunc();
167 // Remember to add new functions also to transform_jacobian and transform_x_back
168 if (tfun == "") {
169 // Nothing to do
170 } else if (tfun == "log") {
171 const Vector& pars = jq.TFuncParameters();
172 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
173 ARTS_USER_ERROR_IF (x[r] <= pars[0],
174 "log-transformation selected for retrieval quantity with\n"
175 "index ", i, " (0-based), but at least one value <= z_min\n"
176 "found for this quantity. This is not allowed.")
177 x[r] = log(x[r] - pars[0]);
178 }
179 } else if (tfun == "log10") {
180 const Vector& pars = jq.TFuncParameters();
181 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
182 ARTS_USER_ERROR_IF (x[r] <= 0,
183 "log10-transformation selected for retrieval quantity with\n"
184 "index ", i, " (0-based), but at least one value <= z_min\n"
185 "found for this quantity. This is not allowed.")
186 x[r] = log10(x[r] - pars[0]);
187 }
188 } else if (tfun == "atanh") {
189 const Vector& pars = jq.TFuncParameters();
190 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
191 ARTS_USER_ERROR_IF (x[r] <= pars[0],
192 "atanh-transformation selected for retrieval quantity with\n"
193 "index ", i, " (0-based), but at least one value <= z_min\n"
194 "found for this quantity. This is not allowed.")
195 ARTS_USER_ERROR_IF (x[r] >= pars[1],
196 "atanh-transformation selected for retrieval quantity with\n"
197 "index ", i, " (0-based), but at least one value is\n"
198 ">= z_max. This is not allowed.")
199 x[r] = atanh(2 * (x[r] - pars[0]) / (pars[1] - pars[0]) - 1);
200 }
201 } else {
202 ARTS_ASSERT(0);
203 }
204 }
205
206 // Apply affine transformations
207 if (any_affine) {
209 jac_ranges_indices(jis_t, any_affine, jqs);
210
211 Vector x_t(jis_t.back()[1] + 1);
212
213 for (Index i = 0; i < jqs.nelem(); ++i) {
214 const RetrievalQuantity& jq = jqs[i];
215 Index col_start = jis[i][0];
216 Index col_extent = jis[i][1] - jis[i][0] + 1;
217 Range col_range(col_start, col_extent);
218 Index col_start_t = jis_t[i][0];
219 Index col_extent_t = jis_t[i][1] - jis_t[i][0] + 1;
220 Range col_range_t(col_start_t, col_extent_t);
221 if (jq.HasAffine()) {
222 Vector t(x[col_range]);
223 t -= jq.OffsetVector();
224 mult(x_t[col_range_t], transpose(jq.TransformationMatrix()), t);
225 } else {
226 x_t[col_range_t] = x[col_range];
227 }
228 }
229 swap(x, x_t);
230 }
231}
232
234 const ArrayOfRetrievalQuantity& jqs,
235 bool revert_functional_transforms) {
236 // Range indices without affine trans
238 bool any_affine;
239 //
240 jac_ranges_indices(jis, any_affine, jqs, true);
241
242 // Revert affine transformations
243 // Apply affine transformations
244 if (any_affine) {
246 jac_ranges_indices(jis_t, any_affine, jqs);
247
248 Vector x(jis.back()[1] + 1);
249
250 for (Index i = 0; i < jqs.nelem(); ++i) {
251 const RetrievalQuantity& jq = jqs[i];
252 Index col_start = jis[i][0];
253 Index col_extent = jis[i][1] - jis[i][0] + 1;
254 Range col_range(col_start, col_extent);
255 Index col_start_t = jis_t[i][0];
256 Index col_extent_t = jis_t[i][1] - jis_t[i][0] + 1;
257 Range col_range_t(col_start_t, col_extent_t);
258 if (jq.HasAffine()) {
259 mult(x[col_range], jq.TransformationMatrix(), x_t[col_range_t]);
260 x[col_range] += jq.OffsetVector();
261 } else {
262 x[col_range] = x_t[col_range_t];
263 }
264 }
265 swap(x_t, x);
266 }
267
268 if (revert_functional_transforms) {
269 // Revert functional transformations
270 for (Index i = 0; i < jqs.nelem(); ++i) {
271 const RetrievalQuantity& jq = jqs[i];
272 const String tfun = jq.TransformationFunc();
273 // Remember to add new functions also to transform_jacobian and transform_x_back
274 if (tfun == "") {
275 // Nothing to do
276 } else if (tfun == "log") {
277 const Vector& pars = jq.TFuncParameters();
278 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
279 x_t[r] = pars[0] + exp(x_t[r]);
280 }
281 } else if (tfun == "log10") {
282 const Vector& pars = jq.TFuncParameters();
283 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
284 x_t[r] = pars[0] + pow(10.0, x_t[r]);
285 }
286 } else if (tfun == "atanh") {
287 const Vector& pars = jq.TFuncParameters();
288 for (Index r = jis[i][0]; r <= jis[i][1]; ++r) {
289 x_t[r] = pars[0] + ((pars[1] - pars[0]) / 2) * (1 + tanh(x_t[r]));
290 }
291 } else {
292 ARTS_ASSERT(0);
293 }
294 }
295 }
296}
297
298/*===========================================================================
299 === Help sub-functions to handle analytical jacobians (in alphabetical order)
300 ===========================================================================*/
301
302// Small help function, to make the code below cleaner
304 ConstMatrixView diy_dq,
305 const Numeric& w) {
306 for (Index irow = 0; irow < diy_dx.nrows(); irow++) {
307 for (Index icol = 0; icol < diy_dx.ncols(); icol++) {
308 diy_dx(irow, icol) += w * diy_dq(irow, icol);
309 }
310 }
311}
312
314 const RetrievalQuantity& jacobian_quantity,
315 ConstTensor3View diy_dpath,
316 const Index& atmosphere_dim,
317 const Ppath& ppath,
318 ConstVectorView ppath_p) {
319 ARTS_ASSERT(jacobian_quantity.Grids().nelem() == atmosphere_dim or jacobian_quantity.Grids().empty());
320
321 // We want here an extrapolation to infinity ->
322 // extremly high extrapolation factor
323 const Numeric extpolfac = 1.0e99;
324
325 if (ppath.np > 1) // Otherwise nothing to do here
326 {
327 // Pressure
328 Index nr1 = jacobian_quantity.Grids().empty() ? 0 : jacobian_quantity.Grids()[0].nelem();
329 ArrayOfGridPos gp_p(ppath.np);
330 if (nr1 > 1) {
331 p2gridpos(gp_p, jacobian_quantity.Grids()[0], ppath_p, extpolfac);
333 } else {
334 gp4length1grid(gp_p);
335 }
336
337 // Latitude
338 Index nr2 = 1;
339 ArrayOfGridPos gp_lat;
340 if (atmosphere_dim > 1) {
341 gp_lat.resize(ppath.np);
342 nr2 = jacobian_quantity.Grids().empty() ? 0 : jacobian_quantity.Grids()[1].nelem();
343 if (nr2 > 1) {
344 gridpos(gp_lat,
345 jacobian_quantity.Grids()[1],
346 ppath.pos(joker, 1),
347 extpolfac);
349 } else {
350 gp4length1grid(gp_lat);
351 }
352 }
353
354 // Longitude
355 ArrayOfGridPos gp_lon;
356 if (atmosphere_dim > 2) {
357 Index nr3 = jacobian_quantity.Grids().empty() ? 0 : jacobian_quantity.Grids()[2].nelem();
358 gp_lon.resize(ppath.np);
359 if (nr3 > 1) {
360 gridpos(gp_lon,
361 jacobian_quantity.Grids()[2],
362 ppath.pos(joker, 2),
363 extpolfac);
365 } else {
366 gp4length1grid(gp_lon);
367 }
368 }
369
370 //- 1D
371 if (atmosphere_dim == 1) {
372 for (Index ip = 0; ip < ppath.np; ip++) {
373 if (gp_p[ip].fd[1] > 0) {
374 from_dpath_to_dx(diy_dx(gp_p[ip].idx, joker, joker),
375 diy_dpath(ip, joker, joker),
376 gp_p[ip].fd[1]);
377 }
378 if (gp_p[ip].fd[0] > 0) {
379 from_dpath_to_dx(diy_dx(gp_p[ip].idx + 1, joker, joker),
380 diy_dpath(ip, joker, joker),
381 gp_p[ip].fd[0]);
382 }
383 }
384 }
385
386 //- 2D
387 else if (atmosphere_dim == 2) {
388 for (Index ip = 0; ip < ppath.np; ip++) {
389 Index ix = nr1 * gp_lat[ip].idx + gp_p[ip].idx;
390 // Low lat, low p
391 if (gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[1] > 0)
392 from_dpath_to_dx(diy_dx(ix, joker, joker),
393 diy_dpath(ip, joker, joker),
394 gp_lat[ip].fd[1] * gp_p[ip].fd[1]);
395 // Low lat, high p
396 if (gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[0] > 0)
397 from_dpath_to_dx(diy_dx(ix + 1, joker, joker),
398 diy_dpath(ip, joker, joker),
399 gp_lat[ip].fd[1] * gp_p[ip].fd[0]);
400 // High lat, low p
401 if (gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[1] > 0)
402 from_dpath_to_dx(diy_dx(ix + nr1, joker, joker),
403 diy_dpath(ip, joker, joker),
404 gp_lat[ip].fd[0] * gp_p[ip].fd[1]);
405 // High lat, high p
406 if (gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[0] > 0)
407 from_dpath_to_dx(diy_dx(ix + nr1 + 1, joker, joker),
408 diy_dpath(ip, joker, joker),
409 gp_lat[ip].fd[0] * gp_p[ip].fd[0]);
410 }
411 }
412
413 //- 3D
414 else if (atmosphere_dim == 3) {
415 for (Index ip = 0; ip < ppath.np; ip++) {
416 Index ix =
417 nr2 * nr1 * gp_lon[ip].idx + nr1 * gp_lat[ip].idx + gp_p[ip].idx;
418 // Low lon, low lat, low p
419 if (gp_lon[ip].fd[1] > 0 && gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[1] > 0)
421 diy_dx(ix, joker, joker),
422 diy_dpath(ip, joker, joker),
423 gp_lon[ip].fd[1] * gp_lat[ip].fd[1] * gp_p[ip].fd[1]);
424 // Low lon, low lat, high p
425 if (gp_lon[ip].fd[1] > 0 && gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[0] > 0)
427 diy_dx(ix + 1, joker, joker),
428 diy_dpath(ip, joker, joker),
429 gp_lon[ip].fd[1] * gp_lat[ip].fd[1] * gp_p[ip].fd[0]);
430 // Low lon, high lat, low p
431 if (gp_lon[ip].fd[1] > 0 && gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[1] > 0)
433 diy_dx(ix + nr1, joker, joker),
434 diy_dpath(ip, joker, joker),
435 gp_lon[ip].fd[1] * gp_lat[ip].fd[0] * gp_p[ip].fd[1]);
436 // Low lon, high lat, high p
437 if (gp_lon[ip].fd[1] > 0 && gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[0] > 0)
439 diy_dx(ix + nr1 + 1, joker, joker),
440 diy_dpath(ip, joker, joker),
441 gp_lon[ip].fd[1] * gp_lat[ip].fd[0] * gp_p[ip].fd[0]);
442
443 // Increase *ix* (to be valid for high lon level)
444 ix += nr2 * nr1;
445
446 // High lon, low lat, low p
447 if (gp_lon[ip].fd[0] > 0 && gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[1] > 0)
449 diy_dx(ix, joker, joker),
450 diy_dpath(ip, joker, joker),
451 gp_lon[ip].fd[0] * gp_lat[ip].fd[1] * gp_p[ip].fd[1]);
452 // High lon, low lat, high p
453 if (gp_lon[ip].fd[0] > 0 && gp_lat[ip].fd[1] > 0 && gp_p[ip].fd[0] > 0)
455 diy_dx(ix + 1, joker, joker),
456 diy_dpath(ip, joker, joker),
457 gp_lon[ip].fd[0] * gp_lat[ip].fd[1] * gp_p[ip].fd[0]);
458 // High lon, high lat, low p
459 if (gp_lon[ip].fd[0] > 0 && gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[1] > 0)
461 diy_dx(ix + nr1, joker, joker),
462 diy_dpath(ip, joker, joker),
463 gp_lon[ip].fd[0] * gp_lat[ip].fd[0] * gp_p[ip].fd[1]);
464 // High lon, high lat, high p
465 if (gp_lon[ip].fd[0] > 0 && gp_lat[ip].fd[0] > 0 && gp_p[ip].fd[0] > 0)
467 diy_dx(ix + nr1 + 1, joker, joker),
468 diy_dpath(ip, joker, joker),
469 gp_lon[ip].fd[0] * gp_lat[ip].fd[0] * gp_p[ip].fd[0]);
470 }
471 }
472 }
473}
474
476 const RetrievalQuantity& jacobian_quantity,
477 ConstMatrixView diy_dpos,
478 const Index& atmosphere_dim,
479 ConstVectorView rtp_pos) {
480 ARTS_ASSERT(jacobian_quantity.Grids().nelem() ==
481 max(atmosphere_dim - 1, Index(1)));
482 ARTS_ASSERT(rtp_pos.nelem() == atmosphere_dim);
483
484 // We want here an extrapolation to infinity ->
485 // extremly high extrapolation factor
486 const Numeric extpolfac = 1.0e99;
487
488 // Handle 1D separately
489 if (atmosphere_dim == 1) {
490 diy_dx(0, joker, joker) = diy_dpos;
491 return;
492 }
493
494 // Latitude
495 Index nr1 = 1;
496 ArrayOfGridPos gp_lat;
497 {
498 gp_lat.resize(1);
499 nr1 = jacobian_quantity.Grids()[0].nelem();
500 if (nr1 > 1) {
501 gridpos(gp_lat,
502 jacobian_quantity.Grids()[0],
503 Vector(1, rtp_pos[1]),
504 extpolfac);
506 } else {
507 gp4length1grid(gp_lat);
508 }
509 }
510
511 // Longitude
512 ArrayOfGridPos gp_lon;
513 if (atmosphere_dim > 2) {
514 gp_lon.resize(1);
515 if (jacobian_quantity.Grids()[1].nelem() > 1) {
516 gridpos(gp_lon,
517 jacobian_quantity.Grids()[1],
518 Vector(1, rtp_pos[2]),
519 extpolfac);
521 } else {
522 gp4length1grid(gp_lon);
523 }
524 }
525
526 //- 2D
527 if (atmosphere_dim == 2) {
528 if (gp_lat[0].fd[1] > 0) {
529 from_dpath_to_dx(diy_dx(gp_lat[0].idx, joker, joker),
530 diy_dpos(joker, joker),
531 gp_lat[0].fd[1]);
532 }
533 if (gp_lat[0].fd[0] > 0) {
534 from_dpath_to_dx(diy_dx(gp_lat[0].idx + 1, joker, joker),
535 diy_dpos(joker, joker),
536 gp_lat[0].fd[0]);
537 }
538 }
539 //- 3D
540 else {
541 Index ix = nr1 * gp_lon[0].idx + gp_lat[0].idx;
542 // Low lon, low lat
543 if (gp_lon[0].fd[1] > 0 && gp_lat[0].fd[1] > 0)
544 from_dpath_to_dx(diy_dx(ix, joker, joker),
545 diy_dpos(joker, joker),
546 gp_lon[0].fd[1] * gp_lat[0].fd[1]);
547 // Low lon, high lat
548 if (gp_lon[0].fd[1] > 0 && gp_lat[0].fd[0] > 0)
549 from_dpath_to_dx(diy_dx(ix + 1, joker, joker),
550 diy_dpos(joker, joker),
551 gp_lon[0].fd[1] * gp_lat[0].fd[0]);
552 // High lon, low lat
553 if (gp_lon[0].fd[0] > 0 && gp_lat[0].fd[1] > 0)
554 from_dpath_to_dx(diy_dx(ix + nr1, joker, joker),
555 diy_dpos(joker, joker),
556 gp_lon[0].fd[0] * gp_lat[0].fd[1]);
557 // High lon, high lat
558 if (gp_lon[0].fd[0] > 0 && gp_lat[0].fd[0] > 0)
559 from_dpath_to_dx(diy_dx(ix + nr1 + 1, joker, joker),
560 diy_dpos(joker, joker),
561 gp_lon[0].fd[0] * gp_lat[0].fd[0]);
562 }
563}
564
566 const ArrayOfArrayOfSpeciesTag& abs_species) {
567 ArrayOfIndex aoi(jacobian_quantities.nelem(), -1);
568
570 if (jacobian_quantities[iq] == Jacobian::Line::VMR) {
571 auto p = std::find_if(abs_species.cbegin(), abs_species.cend(),
572 [qid=jacobian_quantities[iq].QuantumIdentity()](auto& specs){
573 return std::any_of(specs.cbegin(), specs.cend(),
574 [qid](auto& spec){return qid.Isotopologue() == spec.Isotopologue();});
575 });
576 if (p not_eq abs_species.cend()) {
577 aoi[iq] = Index(abs_species.cend() - p);
578 } else {
580 "Could not find ",
581 jacobian_quantities[iq].Subtag(),
582 " in species of abs_species.\n")
583 }
584 } else if (jacobian_quantities[iq] == Jacobian::Special::ArrayOfSpeciesTagVMR) {
585 ArrayOfSpeciesTag atag(jacobian_quantities[iq].Subtag());
586 aoi[iq] = chk_contains("abs_species", abs_species, atag);
587 } else if (jacobian_quantities[iq] == Jacobian::Atm::Particulates) {
588 aoi[iq] = -9999;
589 } else if (jacobian_quantities[iq] == Jacobian::Atm::Electrons) {
590 aoi[iq] = -9999;
591 }
592 )
593
594 return aoi;
595}
596
597ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity& jacobian_quantities, Index np, Index nf, Index ns, bool active) {
598 ArrayOfTensor3 diy_dpath(jacobian_quantities.nelem());
599
600 const Index nn = active ? np * nf : nf;
601 FOR_ANALYTICAL_JACOBIANS_DO(diy_dpath[iq] = Tensor3(np, nn, ns, 0.0);)
602
603 return diy_dpath;
604}
605
606ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity& jacobian_quantities, Index np, Index nf, Index ns, bool active) {
607 ArrayOfTensor3 diy_dx(jacobian_quantities.nelem());
608
609 bool any_affine;
610 ArrayOfArrayOfIndex jacobian_indices;
611 jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
612
613 const Index nn = active ? np * nf : nf;
614 FOR_ANALYTICAL_JACOBIANS_DO2(diy_dx[iq] = Tensor3(jacobian_indices[iq][1] - jacobian_indices[iq][0] + 1, nn, ns, 0.0);)
615
616 return diy_dx;
617}
618
620 const ArrayOfString& scat_species,
621 const bool cloudbox_on) {
622 ArrayOfIndex aoi(jacobian_quantities.nelem(), -1);
623
625 if (cloudbox_on and jacobian_quantities[iq] == Jacobian::Special::ScatteringString) {
626 aoi[iq] = find_first(scat_species, jacobian_quantities[iq].Subtag());
627 ARTS_USER_ERROR_IF (aoi[iq] < 0,
628 "Jacobian quantity with index ", iq, " refers to\n"
629 " ", jacobian_quantities[iq].Subtag(),
630 "\nbut this species could not be found in *scat_species*.")
631 }
632 )
633
634 return aoi;
635}
636
637/*===========================================================================
638 === Other functions, in alphabetical order
639 ===========================================================================*/
640
642 ostringstream& os,
643 const Vector& p_grid,
644 const Vector& lat_grid,
645 const Vector& lon_grid,
646 const Vector& p_retr,
647 const Vector& lat_retr,
648 const Vector& lon_retr,
649 const String& p_retr_name,
650 const String& lat_retr_name,
651 const String& lon_retr_name,
652 const Index& dim) {
653 ARTS_ASSERT(grids.nelem() == dim);
654
655 if (p_retr.nelem() == 0) {
656 os << "The grid vector *" << p_retr_name << "* is empty,"
657 << " at least one pressure level\n"
658 << "should be specified.";
659 return false;
660 } else if (!is_decreasing(p_retr)) {
661 os << "The pressure grid vector *" << p_retr_name << "* is not a\n"
662 << "strictly decreasing vector, which is required.";
663 return false;
664 } else if (p_grid.nelem() == 1 and p_grid.nelem() == p_retr.nelem()) {
665 if (p_grid[0] not_eq p_retr[0]) {
666 os << "Mismatching 1-long grids for " << p_retr_name;
667 return false;
668 }
669
670 // Necessary repeat but grids are OK
671 grids[0] = p_retr;
672 } else if (log(p_retr[0]) > 1.5 * log(p_grid[0]) - 0.5 * log(p_grid[1]) ||
673 log(p_retr[p_retr.nelem() - 1]) <
674 1.5 * log(p_grid[p_grid.nelem() - 1]) -
675 0.5 * log(p_grid[p_grid.nelem() - 2])) {
676 os << "The grid vector *" << p_retr_name << "* is not covered by the\n"
677 << "corresponding atmospheric grid.";
678 return false;
679 } else {
680 // Pressure grid ok, add it to grids
681 grids[0] = p_retr;
682 }
683
684 if (dim >= 2) {
685 // If 2D and 3D atmosphere, check latitude grid
686 if (lat_retr.nelem() == 0) {
687 os << "The grid vector *" << lat_retr_name << "* is empty,"
688 << " at least one latitude\n"
689 << "should be specified for a 2D/3D atmosphere.";
690 return false;
691 } else if (!is_increasing(lat_retr)) {
692 os << "The latitude grid vector *" << lat_retr_name << "* is not a\n"
693 << "strictly increasing vector, which is required.";
694 return false;
695 } else if (lat_grid.nelem() == 1 and lat_grid.nelem() == lat_retr.nelem()) {
696 if (lat_grid[0] not_eq lat_retr[0]) {
697 os << "Mismatching 1-long grids for " << lat_retr_name;
698 return false;
699 }
700
701 // Necessary repeat but grids are OK
702 grids[1] = lat_retr;
703 } else if (lat_retr[0] < 1.5 * lat_grid[0] - 0.5 * lat_grid[1] ||
704 lat_retr[lat_retr.nelem() - 1] >
705 1.5 * lat_grid[lat_grid.nelem() - 1] -
706 0.5 * lat_grid[lat_grid.nelem() - 2]) {
707 os << "The grid vector *" << lat_retr_name << "* is not covered by the\n"
708 << "corresponding atmospheric grid.";
709 return false;
710 } else {
711 // Latitude grid ok, add it to grids
712 grids[1] = lat_retr;
713 }
714 if (dim == 3) {
715 // For 3D atmosphere check longitude grid
716 if (lon_retr.nelem() == 0) {
717 os << "The grid vector *" << lon_retr_name << "* is empty,"
718 << " at least one longitude\n"
719 << "should be specified for a 3D atmosphere.";
720 return false;
721 } else if (!is_increasing(lon_retr)) {
722 os << "The longitude grid vector *" << lon_retr_name << "* is not a\n"
723 << "strictly increasing vector, which is required.";
724 return false;
725 } else if (lon_grid.nelem() == 1 and
726 lon_grid.nelem() == lon_retr.nelem()) {
727 if (lon_grid[0] not_eq lon_retr[0]) {
728 os << "Mismatching 1-long grids for " << lon_retr_name;
729 return false;
730 }
731
732 // Necessary repeat but grids are OK
733 grids[2] = lon_retr;
734 } else if (lon_retr[0] < 1.5 * lon_grid[0] - 0.5 * lon_grid[1] ||
735 lon_retr[lon_retr.nelem() - 1] >
736 1.5 * lon_grid[lon_grid.nelem() - 1] -
737 0.5 * lon_grid[lon_grid.nelem() - 2]) {
738 os << "The grid vector *" << lon_retr_name
739 << "* is not covered by the\n"
740 << "corresponding atmospheric grid.";
741 return false;
742 } else {
743 // Longitude grid ok, add it to grids
744 grids[2] = lon_retr;
745 }
746 }
747 }
748 return true;
749}
750
752 ostringstream& os,
753 const Vector& lat_grid,
754 const Vector& lon_grid,
755 const Vector& lat_retr,
756 const Vector& lon_retr,
757 const String& lat_retr_name,
758 const String& lon_retr_name,
759 const Index& dim) {
760 ARTS_ASSERT(grids.nelem() == max(dim - 1, Index(1)));
761
762 if (dim == 1) {
763 // Here we only need to create a length 1 dummy grid
764 grids[0].resize(1);
765 grids[0][0] = 0;
766 }
767
768 if (dim >= 2) {
769 // If 2D and 3D atmosphere, check latitude grid
770 if (lat_retr.nelem() == 0) {
771 os << "The grid vector *" << lat_retr_name << "* is empty,"
772 << " at least one latitude\n"
773 << "should be specified for a 2D/3D atmosphere.";
774 return false;
775 } else if (!is_increasing(lat_retr)) {
776 os << "The latitude grid vector *" << lat_retr_name << "* is not a\n"
777 << "strictly increasing vector, which is required.";
778 return false;
779 } else if (lat_grid.nelem() == 1 and lat_grid.nelem() == lat_retr.nelem()) {
780 if (lat_grid[0] not_eq lat_retr[0]) {
781 os << "Mismatching 1-long grids for " << lat_retr_name;
782 return false;
783 }
784
785 // Necessary repeat but grids are OK
786 grids[0] = lat_retr;
787 } else if (lat_retr[0] < 1.5 * lat_grid[0] - 0.5 * lat_grid[1] ||
788 lat_retr[lat_retr.nelem() - 1] >
789 1.5 * lat_grid[lat_grid.nelem() - 1] -
790 0.5 * lat_grid[lat_grid.nelem() - 2]) {
791 os << "The grid vector *" << lat_retr_name << "* is not covered by the\n"
792 << "corresponding atmospheric grid.";
793 return false;
794 } else {
795 // Latitude grid ok, add it to grids
796 grids[0] = lat_retr;
797 }
798
799 if (dim == 3) {
800 // For 3D atmosphere check longitude grid
801 if (lon_retr.nelem() == 0) {
802 os << "The grid vector *" << lon_retr_name << "* is empty,"
803 << " at least one longitude\n"
804 << "should be specified for a 3D atmosphere.";
805 return false;
806 } else if (!is_increasing(lon_retr)) {
807 os << "The longitude grid vector *" << lon_retr_name << "* is not a\n"
808 << "strictly increasing vector, which is required.";
809 return false;
810 } else if (lon_grid.nelem() == 1 and
811 lon_grid.nelem() == lon_retr.nelem()) {
812 if (lon_grid[0] not_eq lon_retr[0]) {
813 os << "Mismatching 1-long grids for " << lon_retr_name;
814 return false;
815 }
816
817 // Necessary repeat but grids are OK
818 grids[1] = lon_retr;
819 } else if (lon_retr[0] < 1.5 * lon_grid[0] - 0.5 * lon_grid[1] ||
820 lon_retr[lon_retr.nelem() - 1] >
821 1.5 * lon_grid[lon_grid.nelem() - 1] -
822 0.5 * lon_grid[lon_grid.nelem() - 2]) {
823 os << "The grid vector *" << lon_retr_name
824 << "* is not covered by the\n"
825 << "corresponding atmospheric grid.";
826 return false;
827 } else {
828 // Longitude grid ok, add it to grids
829 grids[1] = lon_retr;
830 }
831 }
832 }
833 return true;
834}
835
837 for (Index i = 0; i < gp.nelem(); i++) {
838 if (gp[i].fd[0] < 0) {
839 gp[i].fd[0] = 0;
840 gp[i].fd[1] = 1;
841 } else if (gp[i].fd[0] > 1) {
842 gp[i].fd[0] = 1;
843 gp[i].fd[1] = 0;
844 }
845 }
846}
847
849 const Vector& x,
850 const Index& poly_coeff) {
851 const Index l = x.nelem();
852
853 ARTS_ASSERT(l > poly_coeff);
854
855 if (b.nelem() != l) b.resize(l);
856
857 if (poly_coeff == 0) {
858 b = 1.0;
859 } else {
860 const Numeric xmin = min(x);
861 const Numeric dx = 0.5 * (max(x) - xmin);
862 //
863 for (Index i = 0; i < l; i++) {
864 b[i] = (x[i] - xmin) / dx - 1.0;
865 b[i] = pow(b[i], int(poly_coeff));
866 }
867 //
868 b -= mean(b);
869 }
870}
871
872void calcBaselineFit(Vector& y_baseline,
873 const Vector& x,
874 const Index& mblock_index,
875 const Sparse& sensor_response,
876 const ArrayOfIndex& sensor_response_pol_grid,
877 const Vector& sensor_response_f_grid,
878 const Matrix& sensor_response_dlos_grid,
879 const RetrievalQuantity& rq,
880 const Index rq_index,
881 const ArrayOfArrayOfIndex& jacobian_indices) {
882 bool is_sine_fit = false;
883 if (rq == Jacobian::Sensor::Polyfit) {
884 is_sine_fit = false;
885 } else if (rq == Jacobian::Sensor::Sinefit) {
886 is_sine_fit = true;
887 } else {
889 "Retrieval quantity is neither a polynomial or a sine "
890 " baseline fit.");
891 }
892
893 // Size and check of sensor_response
894 //
895 const Index nf = sensor_response_f_grid.nelem();
896 const Index npol = sensor_response_pol_grid.nelem();
897 const Index nlos = sensor_response_dlos_grid.nrows();
898
899 // Evaluate basis functions for fits.
900 Vector w, s, c;
901 if (is_sine_fit) {
902 s.resize(nf);
903 c.resize(nf);
904 Numeric period = rq.Grids()[0][0];
905 for (Index f = 0; f < nf; f++) {
906 Numeric a = (sensor_response_f_grid[f] - sensor_response_f_grid[0]) * 2 *
907 PI / period;
908 s[f] = sin(a);
909 c[f] = cos(a);
910 }
911 } else {
912 Numeric poly_coeff = rq.Grids()[0][0];
914 w, sensor_response_f_grid, static_cast<Index>(poly_coeff));
915 }
916
917 // Compute baseline
918 ArrayOfVector jg = rq.Grids();
919 const Index n1 = jg[1].nelem();
920 const Index n2 = jg[2].nelem();
921 const Index n3 = jg[3].nelem();
922 const Range rowind = get_rowindex_for_mblock(sensor_response, mblock_index);
923 const Index row4 = rowind.get_start();
924 Index col4 = jacobian_indices[rq_index][0];
925
926 if (n3 > 1) {
927 col4 += mblock_index * n2 * n1;
928 }
929
930 for (Index l = 0; l < nlos; l++) {
931 const Index row3 = row4 + l * nf * npol;
932 const Index col3 = col4 + l * n1 * (is_sine_fit ? 2 : 1);
933
934 for (Index f = 0; f < nf; f++) {
935 const Index row2 = row3 + f * npol;
936
937 for (Index p = 0; p < npol; p++) {
938 Index col1 = col3;
939 if (n1 > 1) {
940 col1 += p;
941 }
942 if (is_sine_fit) {
943 y_baseline[row2 + p] += x[col1] * s[f] + x[col1 + 1] * c[f];
944 } else {
945 y_baseline[row2 + p] += w[f] * x[col1];
946 }
947 }
948 }
949 }
950}
951
953 const String& unit,
954 const Numeric& vmr,
955 const Numeric& p,
956 const Numeric& t) {
957 if (unit == "rel" || unit == "logrel") {
958 x = 1;
959 } else if (unit == "vmr") {
960 if (vmr == 0) {
961 x = 0;
962 return;
963 }
964 x = 1 / vmr;
965 } else if (unit == "nd") {
966 if (vmr == 0) {
967 x = 0;
968 return;
969 }
970 x = 1 / (vmr * number_density(p, t));
971 } else {
973 "Allowed options for gas species jacobians are "
974 "\"rel\", \"vmr\" and \"nd\".\nYou have selected: ",
975 unit, '\n')
976 }
977}
978
980 const String& unit,
981 const Numeric& vmr,
982 const Numeric& p,
983 const Numeric& t) {
984 if (unit == "rel" || unit == "logrel") {
985 x = vmr;
986 } else if (unit == "vmr") {
987 x = 1;
988 } else if (unit == "nd") {
989 x = 1 / number_density(p, t);
990 } else {
992 "Allowed options for gas species jacobians are "
993 "\"rel\", \"vmr\" and \"nd\".\nYou have selected: ",
994 unit, '\n')
995 }
996}
997
998//======================================================================
999// Propmat partials descriptions
1000//======================================================================
1001
1002bool is_wind_parameter(const RetrievalQuantity& t) noexcept {
1003 return t.Target().isWind();
1004}
1005
1007 return t.Target().isWind() or t.Target().isFrequency();
1008}
1009
1011 return t == Jacobian::Atm::MagneticMagnitude;
1012}
1013
1014bool is_nlte_parameter(const RetrievalQuantity& t) noexcept {
1015 return t == Jacobian::Line::NLTE;
1016}
1017
1018#define ISLINESHAPETYPE(X) \
1019 bool is_pressure_broadening_##X(const RetrievalQuantity& t) noexcept { \
1020 return t == Jacobian::Line::Shape##X##X0 or \
1021 t == Jacobian::Line::Shape##X##X1 or \
1022 t == Jacobian::Line::Shape##X##X2 or \
1023 t == Jacobian::Line::Shape##X##X3; \
1024 }
1029ISLINESHAPETYPE(FVC)
1030ISLINESHAPETYPE(ETA)
1034#undef ISLINESHAPETYPE
1035
1036#define VARISLINESHAPEPARAM(X, Y) (t == Jacobian::Line::Shape##X##Y)
1038 return VARISLINESHAPEPARAM(G0, X0) or VARISLINESHAPEPARAM(D0, X0) or
1039 VARISLINESHAPEPARAM(G2, X0) or VARISLINESHAPEPARAM(D2, X0) or
1040 VARISLINESHAPEPARAM(FVC, X0) or VARISLINESHAPEPARAM(ETA, X0) or
1041 VARISLINESHAPEPARAM(Y, X0) or VARISLINESHAPEPARAM(G, X0) or
1042 VARISLINESHAPEPARAM(DV, X0);
1043}
1044
1046 return VARISLINESHAPEPARAM(G0, X1) or VARISLINESHAPEPARAM(D0, X1) or
1047 VARISLINESHAPEPARAM(G2, X1) or VARISLINESHAPEPARAM(D2, X1) or
1048 VARISLINESHAPEPARAM(FVC, X1) or VARISLINESHAPEPARAM(ETA, X1) or
1049 VARISLINESHAPEPARAM(Y, X1) or VARISLINESHAPEPARAM(G, X1) or
1050 VARISLINESHAPEPARAM(DV, X1);
1051}
1052
1054 return VARISLINESHAPEPARAM(G0, X2) or VARISLINESHAPEPARAM(D0, X2) or
1055 VARISLINESHAPEPARAM(G2, X2) or VARISLINESHAPEPARAM(D2, X2) or
1056 VARISLINESHAPEPARAM(FVC, X2) or VARISLINESHAPEPARAM(ETA, X2) or
1057 VARISLINESHAPEPARAM(Y, X2) or VARISLINESHAPEPARAM(G, X2) or
1058 VARISLINESHAPEPARAM(DV, X2);
1059}
1060
1062 return VARISLINESHAPEPARAM(G0, X3) or VARISLINESHAPEPARAM(D0, X3) or
1063 VARISLINESHAPEPARAM(G2, X3) or VARISLINESHAPEPARAM(D2, X3) or
1064 VARISLINESHAPEPARAM(FVC, X3) or VARISLINESHAPEPARAM(ETA, X3) or
1065 VARISLINESHAPEPARAM(Y, X3) or VARISLINESHAPEPARAM(G, X3) or
1066 VARISLINESHAPEPARAM(DV, X3);
1067}
1068#undef VARISLINESHAPEPARAM
1069
1071 const RetrievalQuantity& t) noexcept {
1075}
1076
1083}
1084
1085bool is_line_parameter(const RetrievalQuantity& t) noexcept {
1086 return t == Jacobian::Type::Line;
1087}
1088
1090 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return (j == Jacobian::Atm::Temperature or j == Jacobian::Special::ArrayOfSpeciesTagVMR or j == Jacobian::Line::VMR or is_frequency_parameter(j));});
1091}
1092
1094 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return (j == Jacobian::Atm::Temperature or j == Jacobian::Line::VMR or j == Jacobian::Special::ArrayOfSpeciesTagVMR or is_frequency_parameter(j));});
1095}
1096
1098 ARTS_USER_ERROR_IF (std::any_of(js.cbegin(), js.cend(), [](auto& j){return is_line_parameter(j);}),
1099 "Line specific parameters are not supported while using continuum tags.\nWe do not track what lines are in the continuum.\n")
1100 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return (j == Jacobian::Atm::Temperature or j == Jacobian::Special::ArrayOfSpeciesTagVMR or is_frequency_parameter(j));});
1101}
1102
1104 ARTS_USER_ERROR_IF (std::any_of(js.cbegin(), js.cend(), [](auto& j){return is_line_parameter(j);}),
1105 "Line specific parameters are not supported while\n using the relaxation matrix line mixing routine.\n We do not yet track individual lines in the relaxation matrix calculations.\n")
1106 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return (j == Jacobian::Atm::Temperature or is_frequency_parameter(j));});
1107}
1108
1110 ARTS_USER_ERROR_IF (std::any_of(js.cbegin(), js.cend(), [](auto& j){return is_line_parameter(j);}),
1111 "Line specific parameters are not supported while using Lookup table.\nWe do not track lines in the Lookup.\n")
1112 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return (j == Jacobian::Atm::Temperature or j == Jacobian::Special::ArrayOfSpeciesTagVMR or is_frequency_parameter(j));});
1113}
1114
1116 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return not (j == Jacobian::Type::Sensor);});
1117}
1118
1120 if (rq == Jacobian::Line::VMR) { // Single tag
1121 for (auto& st : ast) {
1122 if ((rq.QuantumIdentity().isotopologue_index == st.spec_ind) or
1123 (rq.QuantumIdentity().Species() == st.Spec() and
1125 st.Isotopologue().isotname == Species::Joker)))
1126 return true;
1127 }
1128 } else {
1129 return rq.Target().species_array_id == ast;
1130 }
1131
1132 return false;
1133}
1134
1135bool species_match(const RetrievalQuantity& rq, const Species::Species species) {
1136 if (rq == Jacobian::Line::VMR and rq.QuantumIdentity().Species() == species)
1137 return true;
1138 return false;
1139}
1140
1142 const Species::IsotopeRecord& ir) {
1143 auto ir2 = rq.QuantumIdentity().Isotopologue();
1144 if (ir.spec == ir2.spec and (ir.isotname == Species::Joker or ir.isotname == ir2.isotname))
1145 return true;
1146 else
1147 return false;
1148}
1149
1151 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return j == Jacobian::Atm::Temperature;});
1152}
1153
1155 const QuantumIdentifier& line_qid) noexcept {
1156 auto p = std::find_if(js.cbegin(), js.cend(), [&line_qid](auto& j) {
1157 return j == Jacobian::Line::VMR
1158 and j.QuantumIdentity().Species() == line_qid.Species()
1159 and j.QuantumIdentity().Isotopologue() == line_qid.Isotopologue();}
1160 );
1161 if (p not_eq js.cend())
1162 return {true, p -> QuantumIdentity()};
1163 else
1164 return {false, line_qid};
1165}
1166
1168 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return j == Jacobian::Line::Center;});
1169}
1170
1172 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return is_wind_parameter(j);});
1173}
1174
1176 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return is_frequency_parameter(j);});
1177}
1178
1180 return std::any_of(js.cbegin(), js.cend(), [](auto& j){return j.is_mag();});
1181}
1182
1184 auto p = std::find_if(js.cbegin(), js.cend(), [](auto& j){return j == Jacobian::Atm::Temperature;});
1185 if (p not_eq js.cend())
1186 return p -> Target().perturbation;
1187 else
1188 return std::numeric_limits<Numeric>::quiet_NaN();
1189}
1190
1192 auto p = std::find_if(js.cbegin(), js.cend(), [](auto& j){return is_frequency_parameter(j);});
1193 if (p not_eq js.cend())
1194 return p -> Target().perturbation;
1195 else
1196 return std::numeric_limits<Numeric>::quiet_NaN();
1197}
1198
1200 auto p = std::find_if(js.cbegin(), js.cend(), [](auto& j){return j.is_mag();});
1201 if (p not_eq js.cend())
1202 return p -> Target().perturbation;
1203 else
1204 return std::numeric_limits<Numeric>::quiet_NaN();
1205}
1206
1207namespace Jacobian {
1208std::ostream& operator<<(std::ostream& os, const Target& x) {
1209 os << x.TargetType() << ' ';
1210 switch (toType(x.TargetType())) {
1211 case Type::Atm:
1212 os << x.atm;
1213 break;
1214 case Type::Line:
1215 os << x.line;
1216 break;
1217 case Type::Sensor:
1218 os << x.sensor;
1219 break;
1220 case Type::Special:
1221 os << x.special;
1222 break;
1223 case Type::FINAL:
1224 os << "BAD STATE";
1225 break;
1226 }
1227 if (x.needQuantumIdentity()) os << ' ' << x.qid;
1228 if (x.needArrayOfSpeciesTag()) os << ' ' << x.species_array_id;
1229 if (x.needString()) os << ' ' << x.string_id;
1230 if (std::isnormal(x.perturbation)) os << ' ' << x.perturbation;
1231
1232 return os;
1233}
1234} // namespace Jacobian
Array< Index > ArrayOfIndex
An array of Index.
Definition: array.h:136
Index find_first(const Array< base > &x, const base &w)
Find first occurance.
Definition: array.h:190
base max(const Array< base > &x)
Max function.
Definition: array.h:145
base min(const Array< base > &x)
Min function.
Definition: array.h:161
The global header file for ARTS.
Constants of physical expressions as constexpr.
TimeStep mean(const ArrayOfTimeStep &dt)
Returns the mean time step.
Definition: artstime.cc:190
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
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Matrix.
Definition: matpackI.h:1065
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
A constant view of a Tensor3.
Definition: matpackIII.h:133
A constant view of a Vector.
Definition: matpackI.h:521
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The MatrixView class.
Definition: matpackI.h:1188
The Matrix class.
Definition: matpackI.h:1285
The range class.
Definition: matpackI.h:160
constexpr Index get_start() const noexcept
Returns the start index of the range.
Definition: matpackI.h:341
Deals with internal derivatives, Jacobian definition, and OEM calculations.
Definition: jacobian.h:325
const Matrix & TransformationMatrix() const
Definition: jacobian.h:499
bool HasAffine() const
Definition: jacobian.h:496
const String & Mode() const
Returns the mode.
Definition: jacobian.h:394
const String & TransformationFunc() const
Definition: jacobian.h:497
Index nelem() const
Number of elements in the grids.
Definition: jacobian.h:422
const Vector & TFuncParameters() const
Definition: jacobian.h:498
const Vector & OffsetVector() const
Definition: jacobian.h:500
Jacobian::Target & Target()
Get the Jacobian Target.
Definition: jacobian.h:431
const QuantumIdentifier & QuantumIdentity() const
Returns the identity of this Jacobian.
Definition: jacobian.h:445
const ArrayOfVector & Grids() const
Returns the grids of the retrieval.
Definition: jacobian.h:408
const String & Subtag() const
Returns the sub-tag.
Definition: jacobian.h:365
The Tensor3View class.
Definition: matpackIII.h:251
The Tensor3 class.
Definition: matpackIII.h:352
The Vector class.
Definition: matpackI.h:910
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
#define ARTS_USER_ERROR(...)
Definition: debug.h:169
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void gp4length1grid(ArrayOfGridPos &gp)
Grid position matching a grid of length 1.
void dxdvmrscf(Numeric &x, const String &unit, const Numeric &vmr, const Numeric &p, const Numeric &t)
Scale factor for conversion of derivatives with respect to VMR.
Definition: jacobian.cc:979
bool species_iso_match(const RetrievalQuantity &rq, const Species::IsotopeRecord &ir)
Returns if the Retrieval quantity is VMR derivative for all the species in the species tags.
Definition: jacobian.cc:1141
ArrayOfIndex get_pointers_for_analytical_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:565
void vmrunitscf(Numeric &x, const String &unit, const Numeric &vmr, const Numeric &p, const Numeric &t)
Scale factor for conversion between gas species units.
Definition: jacobian.cc:952
#define VARISLINESHAPEPARAM(X, Y)
Definition: jacobian.cc:1036
bool supports_propmat_clearsky(const ArrayOfRetrievalQuantity &js)
Returns if the array supports propagation matrix derivatives.
Definition: jacobian.cc:1115
void diy_from_pos_to_rgrids(Tensor3View diy_dx, const RetrievalQuantity &jacobian_quantity, ConstMatrixView diy_dpos, const Index &atmosphere_dim, ConstVectorView rtp_pos)
diy_from_pos_to_rgrids
Definition: jacobian.cc:475
bool do_line_center_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a line center derivative.
Definition: jacobian.cc:1167
ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:606
void polynomial_basis_func(Vector &b, const Vector &x, const Index &poly_coeff)
Calculates polynomial basis functions.
Definition: jacobian.cc:848
bool supports_CIA(const ArrayOfRetrievalQuantity &js)
Returns if the array supports CIA derivatives.
Definition: jacobian.cc:1089
void transform_jacobian(Matrix &jacobian, const Vector x, const ArrayOfRetrievalQuantity &jqs)
Applies both functional and affine transformations.
Definition: jacobian.cc:91
bool is_wind_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a wind parameter in propagation matrix calculations.
Definition: jacobian.cc:1002
bool is_lineshape_parameter_X3(const RetrievalQuantity &t) noexcept
Definition: jacobian.cc:1061
bool is_lineshape_parameter_X1(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a X1 derivative.
Definition: jacobian.cc:1045
bool is_lineshape_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a G0, D0, G2, D2, ETA, FVC, Y, G, DV derivative.
Definition: jacobian.cc:1077
bool supports_relaxation_matrix(const ArrayOfRetrievalQuantity &js)
Returns if the array supports relaxation matrix derivatives.
Definition: jacobian.cc:1103
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a frequency derivative.
Definition: jacobian.cc:1175
bool is_derived_magnetic_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a derived magnetic parameter.
Definition: jacobian.cc:1010
bool supports_continuum(const ArrayOfRetrievalQuantity &js)
Returns if the array supports continuum derivatives.
Definition: jacobian.cc:1097
bool species_match(const RetrievalQuantity &rq, const ArrayOfSpeciesTag &ast)
Returns if the Retrieval quantity is VMR derivative for all the species in the species tags.
Definition: jacobian.cc:1119
bool do_magnetic_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a magnetic derivative.
Definition: jacobian.cc:1179
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
Definition: jacobian.cc:1150
#define ISLINESHAPETYPE(X)
Definition: jacobian.cc:1018
ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:597
bool supports_hitran_xsec(const ArrayOfRetrievalQuantity &js)
Returns if the array supports HITRAN cross-section derivatives.
Definition: jacobian.cc:1093
bool supports_lookup(const ArrayOfRetrievalQuantity &js)
Returns if the array supports lookup table derivatives.
Definition: jacobian.cc:1109
ArrayOfIndex get_pointers_for_scat_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfString &scat_species, const bool cloudbox_on)
Help function for analytical jacobian calculations.
Definition: jacobian.cc:619
bool is_lineshape_parameter_X2(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a X2 derivative.
Definition: jacobian.cc:1053
void transform_x_back(Vector &x_t, const ArrayOfRetrievalQuantity &jqs, bool revert_functional_transforms)
Handles back-transformations of the state vector.
Definition: jacobian.cc:233
Numeric magnetic_field_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the magnetic field perturbation if it exists.
Definition: jacobian.cc:1199
constexpr Numeric NAT_LOG_TEN
Definition: jacobian.cc:37
void diy_from_path_to_rgrids(Tensor3View diy_dx, const RetrievalQuantity &jacobian_quantity, ConstTensor3View diy_dpath, const Index &atmosphere_dim, const Ppath &ppath, ConstVectorView ppath_p)
Maps jacobian data for points along the propagation path, to jacobian retrieval grid data.
Definition: jacobian.cc:313
bool is_nlte_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a NLTE parameter.
Definition: jacobian.cc:1014
bool do_wind_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants a wind-based frequency derivative derivative.
Definition: jacobian.cc:1171
bool is_frequency_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a frequency parameter in propagation matrix calculations.
Definition: jacobian.cc:1006
constexpr Numeric PI
Definition: jacobian.cc:38
bool check_retrieval_grids(ArrayOfVector &grids, ostringstream &os, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &p_retr, const Vector &lat_retr, const Vector &lon_retr, const String &p_retr_name, const String &lat_retr_name, const String &lon_retr_name, const Index &dim)
Check that the retrieval grids are defined for each atmosphere dim.
Definition: jacobian.cc:641
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the temperature perturbation if it exists.
Definition: jacobian.cc:1183
jacobianVMRcheck do_vmr_jacobian(const ArrayOfRetrievalQuantity &js, const QuantumIdentifier &line_qid) noexcept
Returns the required info for VMR Jacobian.
Definition: jacobian.cc:1154
void jac_ranges_indices(ArrayOfArrayOfIndex &jis, bool &any_affine, const ArrayOfRetrievalQuantity &jqs, const bool &before_affine)
Determines the index range inside x and the Jacobian for each retrieval quantity.
Definition: jacobian.cc:46
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity &js) noexcept
Returns the frequency perturbation if it exists.
Definition: jacobian.cc:1191
void from_dpath_to_dx(MatrixView diy_dx, ConstMatrixView diy_dq, const Numeric &w)
Definition: jacobian.cc:303
void jacobian_type_extrapol(ArrayOfGridPos &gp)
Adopts grid positions to extrapolation used for jacobians.
Definition: jacobian.cc:836
bool is_lineshape_parameter_X0(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a X0 derivative.
Definition: jacobian.cc:1037
void transform_x(Vector &x, const ArrayOfRetrievalQuantity &jqs)
Handles transformations of the state vector.
Definition: jacobian.cc:156
ostream & operator<<(ostream &os, const RetrievalQuantity &ot)
Definition: jacobian.cc:40
void calcBaselineFit(Vector &y_baseline, const Vector &x, const Index &mblock_index, const Sparse &sensor_response, const ArrayOfIndex &sensor_response_pol_grid, const Vector &sensor_response_f_grid, const Matrix &sensor_response_dlos_grid, const RetrievalQuantity &rq, const Index rq_index, const ArrayOfArrayOfIndex &jacobian_indices)
Calculate baseline fit.
Definition: jacobian.cc:872
bool is_lineshape_parameter_bar_linemixing(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a G0, D0, G2, D2, ETA, FVC derivative.
Definition: jacobian.cc:1070
bool is_line_parameter(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is related to the absorption line.
Definition: jacobian.cc:1085
Routines for setting up the jacobian.
bool is_pressure_broadening_D0(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a D0 derivative.
bool is_pressure_broadening_Y(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a Y derivative.
bool is_pressure_broadening_DV(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a DV derivative.
bool is_pressure_broadening_G0(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a G0 derivative.
bool is_pressure_broadening_FVC(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a FVC derivative.
bool is_pressure_broadening_D2(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a D2 derivative.
bool is_pressure_broadening_G(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a G derivative.
bool is_pressure_broadening_G2(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a G0 derivative.
bool is_pressure_broadening_ETA(const RetrievalQuantity &t) noexcept
Returns if the Retrieval quantity is a ETA derivative.
#define FOR_ANALYTICAL_JACOBIANS_DO2(what_to_do)
Definition: jacobian.h:555
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:547
Linear algebra functions.
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:218
bool is_decreasing(ConstVectorView x)
Checks if a vector is sorted in reversed order and is strictly decreasing.
Definition: logic.cc:255
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
void swap(ComplexVector &v1, ComplexVector &v2) noexcept
Swaps two objects.
const Joker joker
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code:
constexpr Numeric ln_10
Natural logarithm of 10.
std::ostream & operator<<(std::ostream &os, const Target &x)
Definition: jacobian.cc:1208
constexpr std::string_view Joker
Definition: isotopologues.h:13
This file contains declerations of functions of physical character.
constexpr Numeric number_density(Numeric p, Numeric t) noexcept
number_density
Definition: physics_funcs.h:70
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
Range get_rowindex_for_mblock(const Sparse &sensor_response, const Index &mblock_index)
Returns the "range" of y corresponding to a measurement block.
Definition: rte.cc:1114
Declaration of functions in rte.cc.
void p2gridpos(ArrayOfGridPos &gp, ConstVectorView old_pgrid, ConstVectorView new_pgrid, const Numeric &extpolfac)
Calculates grid positions for pressure values.
Header file for special_interp.cc.
Holds all information required for individual partial derivatives.
Definition: jacobian.h:108
Numeric perturbation
Perturbations for methods where theoretical computations are impossible or plain slow.
Definition: jacobian.h:125
String string_id
ID for some of the Special types of partial derivatives.
Definition: jacobian.h:134
bool needArrayOfSpeciesTag() const noexcept
Does this type need the ArrayOfSpeciesTag?
Definition: jacobian.h:307
Line line
Type of line quantity.
Definition: jacobian.h:116
ArrayOfSpeciesTag species_array_id
ID for some of the Special types of partial derivatives.
Definition: jacobian.h:131
std::string_view TargetType() const noexcept
Return type as string.
Definition: jacobian.h:199
Sensor sensor
Type of sensor quantity.
Definition: jacobian.h:119
Atm atm
Type of atm quantity.
Definition: jacobian.h:113
bool needString() const noexcept
Does this type need the String?
Definition: jacobian.h:312
QuantumIdentifier qid
ID for the Line types of partial derivatives.
Definition: jacobian.h:128
bool needQuantumIdentity() const noexcept
Does this type need the QuantumIdentifier?
Definition: jacobian.h:302
Special special
Type of special quantity.
Definition: jacobian.h:122
The structure to describe a propagation path and releated quantities.
Definition: ppath_struct.h:17
Index np
Number of points describing the ppath.
Definition: ppath_struct.h:21
Matrix pos
The distance between start pos and the last position in pos.
Definition: ppath_struct.h:33
A logical struct for global quantum numbers with species identifiers.
Species::Species Species() const noexcept
Species::IsotopeRecord Isotopologue() const noexcept
The Sparse class.
Definition: matpackII.h:75
Struct containing all information needed about one isotope.
Definition: isotopologues.h:16
Species spec
Species type as defined in species.h.
Definition: isotopologues.h:18
std::string_view isotname
A custom name that is unique for this Species type.
Definition: isotopologues.h:21
Deals with whether or not we should do a VMR derivative.
Definition: jacobian.h:1277
#define w
#define a
#define c
#define b