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