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