ARTS 2.5.4 (git: 31ce4f0e)
test_interpolation.cc
Go to the documentation of this file.
1/* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.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
18#include <cmath>
19#include <iostream>
20
21#include "array.h"
22#include "arts_conversions.h"
23#include "auto_md.h"
24#include "interpolation.h"
26#include "math_funcs.h"
27#include "matpackVII.h"
28#include "xml_io.h"
29
30void test01() {
31 cout << "Simple interpolation cases\n"
32 << "--------------------------\n";
33 // Vector og(5,5,-1); // 5,4,3,2,1
34 Vector og(1, 5, +1); // 1, 2, 3, 4, 5
35 Vector ng(2, 5, 0.25); // 2.0, 2,25, 2.5, 2.75, 3.0
36
37 cout << "og:\n" << og << "\n";
38 cout << "ng:\n" << ng << "\n";
39
40 // To store the grid positions:
41 ArrayOfGridPos gp(ng.nelem());
42
43 gridpos(gp, og, ng);
44 cout << "gp:\n" << gp << "\n";
45
46 cout << "1D:\n"
47 << "---\n";
48 {
49 // To store interpolation weights:
50 Matrix itw(gp.nelem(), 2);
51 interpweights(itw, gp);
52
53 cout << "itw:\n" << itw << "\n";
54
55 // Original field:
56 Vector of(og.nelem(), 0);
57 of[2] = 10; // 0, 0, 10, 0, 0
58
59 cout << "of:\n" << of << "\n";
60
61 // Interpolated field:
62 Vector nf(ng.nelem());
63
64 interp(nf, itw, of, gp);
65
66 cout << "nf:\n" << nf << "\n";
67 }
68
69 cout << "Blue 2D:\n"
70 << "--------\n";
71 {
72 // To store interpolation weights:
73 Matrix itw(gp.nelem(), 4);
74 interpweights(itw, gp, gp);
75
76 cout << "itw:\n" << itw << "\n";
77
78 // Original field:
79 Matrix of(og.nelem(), og.nelem(), 0);
80 of(2, 2) = 10; // 0 Matrix with 10 in the middle
81
82 cout << "of:\n" << of << "\n";
83
84 // Interpolated field:
85 Vector nf(ng.nelem());
86
87 interp(nf, itw, of, gp, gp);
88
89 cout << "nf:\n" << nf << "\n";
90 }
91
92 cout << "Blue 6D:\n"
93 << "--------\n";
94 {
95 // To store interpolation weights:
96 Matrix itw(gp.nelem(), 64);
97 interpweights(itw, gp, gp, gp, gp, gp, gp);
98
99 // cout << "itw:\n" << itw << "\n";
100
101 // Original field:
102 Index n = og.nelem();
103 Tensor6 of(n, n, n, n, n, n, 0);
104 of(2, 2, 2, 2, 2, 2) = 10; // 0 Tensor with 10 in the middle
105
106 // cout << "of:\n" << of << "\n";
107
108 // Interpolated field:
109 Vector nf(ng.nelem());
110
111 interp(nf, itw, of, gp, gp, gp, gp, gp, gp);
112
113 cout << "nf:\n" << nf << "\n";
114 }
115
116 cout << "Green 2D:\n"
117 << "---------\n";
118 {
119 // To store interpolation weights:
120 Tensor3 itw(gp.nelem(), gp.nelem(), 4);
121 interpweights(itw, gp, gp);
122
123 for (Index i = 0; i < itw.ncols(); ++i)
124 cout << "itw " << i << ":\n"
125 << itw(Range(joker), Range(joker), i) << "\n";
126
127 // Original field:
128 Matrix of(og.nelem(), og.nelem(), 0);
129 of(2, 2) = 10; // 0 Matrix with 10 in the middle
130
131 cout << "of:\n" << of << "\n";
132
133 // Interpolated field:
134 Matrix nf(ng.nelem(), ng.nelem());
135
136 interp(nf, itw, of, gp, gp);
137
138 cout << "nf:\n" << nf << "\n";
139 }
140
141 cout << "Green 6D:\n"
142 << "---------\n";
143 {
144 // To store interpolation weights:
145 Tensor7 itw(gp.nelem(), gp.nelem(), gp.nelem(), gp.nelem(), gp.nelem(),
146 gp.nelem(), 64);
147 interpweights(itw, gp, gp, gp, gp, gp, gp);
148
149 // Original field:
150 Tensor6 of(og.nelem(), og.nelem(), og.nelem(), og.nelem(), og.nelem(),
151 og.nelem(), 0);
152 of(2, 2, 2, 2, 2, 2) = 10; // 0 Tensor with 10 in the middle
153
154 cout << "Middle slice of of:\n"
155 << of(2, 2, 2, 2, Range(joker), Range(joker)) << "\n";
156
157 // Interpolated field:
158 Tensor6 nf(ng.nelem(), ng.nelem(), ng.nelem(), ng.nelem(), ng.nelem(),
159 ng.nelem());
160
161 interp(nf, itw, of, gp, gp, gp, gp, gp, gp);
162
163 cout << "Last slice of nf:\n"
164 << nf(4, 4, 4, 4, Range(joker), Range(joker)) << "\n";
165 }
166}
167
168void test02(Index n) {
169 cout << "Test whether for loop or iterator are quicker\n"
170 << "a) for loop\n"
171 << "---------------------------------------------\n";
172
173 Vector a(n);
174 for (Index i = 0; i < a.nelem(); ++i) a[i] = (Numeric)i;
175}
176
177void test03(Index n) {
178 cout << "Test whether for loop or iterator are quicker\n"
179 << "b) iterator\n"
180 << "---------------------------------------------\n";
181
182 Vector a(n);
183 Iterator1D ai = a.begin();
184 const Iterator1D ae = a.end();
185 Index i = 0;
186 for (; ai != ae; ++ai, ++i) *ai = (Numeric)i;
187}
188
189// Result: Both are almost equally fast, with a slight advantage of
190// the for loop if compiler optimization is enabled.
191
192void test04() {
193 cout << "Green type interpolation of all pages of a Tensor3\n";
194
195 // The original Tensor is called a, the new one n.
196
197 // 10 pages, 20 rows, 30 columns, all grids are: 1,2,3
198 Vector a_pgrid(1, 3, 1), a_rgrid(1, 3, 1), a_cgrid(1, 3, 1);
199 Tensor3 a(a_pgrid.nelem(), a_rgrid.nelem(), a_cgrid.nelem());
200
201 a = 0;
202 // Put some simple numbers in the middle of each page:
203 a(0, 1, 1) = 10;
204 a(1, 1, 1) = 20;
205 a(2, 1, 1) = 30;
206
207 // New row and column grids:
208 // 1, 1.5, 2, 2.5, 3
209 Vector n_rgrid(1, 5, .5), n_cgrid(1, 5, .5);
210 Tensor3 n(a_pgrid.nelem(), n_rgrid.nelem(), n_cgrid.nelem());
211
212 // So, n has the same number of pages as a, but more rows and columns.
213
214 // Get the grid position arrays:
215 ArrayOfGridPos n_rgp(n_rgrid.nelem()); // For row grid positions.
216 ArrayOfGridPos n_cgp(n_cgrid.nelem()); // For column grid positions.
217
218 gridpos(n_rgp, a_rgrid, n_rgrid);
219 gridpos(n_cgp, a_cgrid, n_cgrid);
220
221 // Get the interpolation weights:
222 Tensor3 itw(n_rgrid.nelem(), n_cgrid.nelem(), 4);
223 interpweights(itw, n_rgp, n_cgp);
224
225 // Do a "green" interpolation for all pages of the Tensor a:
226
227 for (Index i = 0; i < a.npages(); ++i) {
228 // Select the current page of both a and n:
230 MatrixView np = n(i, Range(joker), Range(joker));
231
232 // Do the interpolation:
233 interp(np, itw, ap, n_rgp, n_cgp);
234
235 // Note that this is efficient, because interpolation weights and
236 // grid positions are re-used.
237 }
238
239 cout << "Original field:\n";
240 for (Index i = 0; i < a.npages(); ++i)
241 cout << "page " << i << ":\n" << a(i, Range(joker), Range(joker)) << "\n";
242
243 cout << "Interpolated field:\n";
244 for (Index i = 0; i < n.npages(); ++i)
245 cout << "page " << i << ":\n" << n(i, Range(joker), Range(joker)) << "\n";
246}
247
248void test05() {
249 cout << "Very simple interpolation case\n";
250
251 Vector og(1, 5, +1); // 1, 2, 3, 4, 5
252 Vector ng(2, 5, 0.25); // 2.0, 2,25, 2.5, 2.75, 3.0
253
254 cout << "Original grid:\n" << og << "\n";
255 cout << "New grid:\n" << ng << "\n";
256
257 // To store the grid positions:
258 ArrayOfGridPos gp(ng.nelem());
259
260 gridpos(gp, og, ng);
261 cout << "Grid positions:\n" << gp;
262
263 // To store interpolation weights:
264 Matrix itw(gp.nelem(), 2);
265 interpweights(itw, gp);
266
267 cout << "Interpolation weights:\n" << itw << "\n";
268
269 // Original field:
270 Vector of(og.nelem(), 0);
271 of[2] = 10; // 0, 0, 10, 0, 0
272
273 cout << "Original field:\n" << of << "\n";
274
275 // Interpolated field:
276 Vector nf(ng.nelem());
277
278 interp(nf, itw, of, gp);
279
280 cout << "New field:\n" << nf << "\n";
281}
282
283void test06() {
284 cout << "Simple extrapolation cases\n"
285 << "--------------------------\n";
286 // Vector og(5,5,-1); // 5,4,3,2,1
287 Vector og(1, 5, +1); // 1, 2, 3, 4, 5
288 Vector ng{0.9, 1.5, 3, 4.5, 5.1}; // 0.9, 1.5, 3, 4.5, 5.1
289
290 cout << "og:\n" << og << "\n";
291 cout << "ng:\n" << ng << "\n";
292
293 // To store the grid positions:
294 ArrayOfGridPos gp(ng.nelem());
295
296 gridpos(gp, og, ng);
297 cout << "gp:\n" << gp << "\n";
298
299 cout << "1D:\n"
300 << "---\n";
301 {
302 // To store interpolation weights:
303 Matrix itw(gp.nelem(), 2);
304 interpweights(itw, gp);
305
306 cout << "itw:\n" << itw << "\n";
307
308 // Original field:
309 Vector of(og.nelem(), 0);
310 for (Index i = 0; i < og.nelem(); ++i)
311 of[i] = (Numeric)(10 * (i + 1)); // 10, 20, 30, 40, 50
312
313 cout << "of:\n" << of << "\n";
314
315 // Interpolated field:
316 Vector nf(ng.nelem());
317
318 interp(nf, itw, of, gp);
319
320 cout << "nf:\n" << nf << "\n";
321 }
322
323 cout << "Blue 2D:\n"
324 << "--------\n";
325 {
326 // To store interpolation weights:
327 Matrix itw(gp.nelem(), 4);
328 interpweights(itw, gp, gp);
329
330 cout << "itw:\n" << itw << "\n";
331
332 // Original field:
333 Matrix of(og.nelem(), og.nelem(), 0);
334 of(2, 2) = 10; // 0 Matrix with 10 in the middle
335
336 cout << "of:\n" << of << "\n";
337
338 // Interpolated field:
339 Vector nf(ng.nelem());
340
341 interp(nf, itw, of, gp, gp);
342
343 cout << "nf:\n" << nf << "\n";
344 }
345
346 cout << "Blue 6D:\n"
347 << "--------\n";
348 {
349 // To store interpolation weights:
350 Matrix itw(gp.nelem(), 64);
351 interpweights(itw, gp, gp, gp, gp, gp, gp);
352
353 // cout << "itw:\n" << itw << "\n";
354
355 // Original field:
356 Index n = og.nelem();
357 Tensor6 of(n, n, n, n, n, n, 0);
358 of(2, 2, 2, 2, 2, 2) = 10; // 0 Tensor with 10 in the middle
359
360 // cout << "of:\n" << of << "\n";
361
362 // Interpolated field:
363 Vector nf(ng.nelem());
364
365 interp(nf, itw, of, gp, gp, gp, gp, gp, gp);
366
367 cout << "nf:\n" << nf << "\n";
368 }
369
370 cout << "Green 2D:\n"
371 << "---------\n";
372 {
373 // To store interpolation weights:
374 Tensor3 itw(gp.nelem(), gp.nelem(), 4);
375 interpweights(itw, gp, gp);
376
377 for (Index i = 0; i < itw.ncols(); ++i)
378 cout << "itw " << i << ":\n"
379 << itw(Range(joker), Range(joker), i) << "\n";
380
381 // Original field:
382 Matrix of(og.nelem(), og.nelem(), 0);
383 of(2, 2) = 10; // 0 Matrix with 10 in the middle
384
385 cout << "of:\n" << of << "\n";
386
387 // Interpolated field:
388 Matrix nf(ng.nelem(), ng.nelem());
389
390 interp(nf, itw, of, gp, gp);
391
392 cout << "nf:\n" << nf << "\n";
393 }
394
395 cout << "Green 6D:\n"
396 << "---------\n";
397 {
398 // To store interpolation weights:
399 Tensor7 itw(gp.nelem(), gp.nelem(), gp.nelem(), gp.nelem(), gp.nelem(),
400 gp.nelem(), 64);
401 interpweights(itw, gp, gp, gp, gp, gp, gp);
402
403 // Original field:
404 Tensor6 of(og.nelem(), og.nelem(), og.nelem(), og.nelem(), og.nelem(),
405 og.nelem(), 0);
406 of(2, 2, 2, 2, 2, 2) = 10; // 0 Tensor with 10 in the middle
407
408 cout << "Middle slice of of:\n"
409 << of(2, 2, 2, 2, Range(joker), Range(joker)) << "\n";
410
411 // Interpolated field:
412 Tensor6 nf(ng.nelem(), ng.nelem(), ng.nelem(), ng.nelem(), ng.nelem(),
413 ng.nelem());
414
415 interp(nf, itw, of, gp, gp, gp, gp, gp, gp);
416
417 cout << "Last slice of nf:\n"
418 << nf(4, 4, 4, 4, Range(joker), Range(joker)) << "\n";
419 }
420}
421
422void test09() {
423 // Original and new grid
424 Vector og(1, 5, +1); // 1, 2, 3, 4, 5
425 Vector ng(2, 9, 0.25); // 2.0, 2,25, 2.5, 2.75, 3.0 ... 4.0
426
427 // Original data
428 Vector yi{5, -2, 50, 2, 1};
429
430 std::cout << og << '\n' << yi << '\n' << '\n';
431
432 // Do the fit point-by-point:
433 for (auto x : ng) {
434 // Simple linear interpolation:
435 GridPos gp;
436 gridpos(gp, og, x);
437
438 // General Lagrange interpolation, special case interpolation order 1:
439 const LagrangeInterpolation lag(0, x, og);
440
441 std::cout << "gp point: " << gp << "lag point:\n" << lag;
442
443 // Linear interpolation weights:
444 Vector iwgp(2);
445 interpweights(iwgp, gp);
446
447 // General Lagrange interpolation weights:
448 const Vector iwlag = interpweights(lag);
449
450 std::cout << "gp iw: " << iwgp << "\nlag iw: " << iwlag << '\n';
451
452 std::cout << "gp res: " << interp(iwgp, yi, gp)
453 << "\nlag res: " << interp(yi, iwlag, lag) << '\n'
454 << '\n';
455 }
456
457 // Linear interpolation of all points
458 ArrayOfGridPos gp(ng.nelem());
459 gridpos(gp, og, ng);
460 Matrix gp_iw(gp.nelem(), 2);
461 interpweights(gp_iw, gp);
462 Vector gp_y(ng.nelem());
463 interp(gp_y, gp_iw, yi, gp);
464 std::cout << "gp: " << gp_y << '\n';
465
466 // General Lagrange interpolation, special case interpolation order 1:
467 const auto lag = Interpolation::LagrangeVector(ng, og, 1);
468 std::cout << "lag: " << reinterp(yi, interpweights(lag), lag) << '\n';
469}
470
471void test11() {
472 constexpr int N = 20;
473 Vector x0(1, N, +1); // 1, 2, 3, 4, 5 ... 10
474 Vector x1(1, N, +2); // 1, 2, 3, 4, 5 ... 10
475 Vector x0n(100);
476 Vector x1n(100);
477 nlinspace(x0n, x0[0], x0[N - 1], 100);
478 nlinspace(x1n, x1[0], x1[N - 1], 100);
479 Matrix yi(N, N);
480 for (Index i = 0; i < N; i++)
481 for (Index j = 0; j < N; j++)
482 yi(i, j) = std::exp(-x0[i] / 3.14 + x1[j] / 3.14);
483
484 std::cerr << x0 << '\n';
485 std::cerr << x1 << '\n';
486 std::cerr << yi << '\n';
487
488 constexpr Index order = 5;
489 {
490 const auto lag0 = Interpolation::LagrangeVector(
491 x0n, x0, order, 0.5, false, Interpolation::GridType::Standard);
492 const auto lag1 = Interpolation::LagrangeVector(
493 x1n, x1, order, 0.5, false, Interpolation::GridType::Standard);
494 const auto iwlag = interpweights(lag0, lag1);
495 std::cout << x0n << '\n'
496 << x1n << '\n'
497 << reinterp(yi, iwlag, lag0, lag1) << '\n';
498 }
499}
500
501constexpr bool is_around(Numeric x, Numeric x0, Numeric e = 1e-12) {
502 return x - x0 < e and x - x0 > -e;
503}
504
505void test12() {
506 constexpr int N = 10;
507 constexpr std::array<Numeric, N> y{1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
508 constexpr std::array<Numeric, N> y2{1, 2 * 2, 3 * 3, 4 * 4, 5 * 5,
509 6 * 6, 7 * 7, 8 * 8, 9 * 9, 10 * 10};
510 constexpr std::array<Numeric, N> y3{
511 1, 2 * 2 * 2, 3 * 3 * 3, 4 * 4 * 4, 5 * 5 * 5,
512 6 * 6 * 6, 7 * 7 * 7, 8 * 8 * 8, 9 * 9 * 9, 10 * 10 * 10};
513 constexpr Numeric x = 1.5;
514
515 // Set up the interpolation Lagranges
516 constexpr Interpolation::FixedLagrange<1> lin(0, x,
517 std::array<Numeric, 2>{1, 2});
519 0, x, std::array<Numeric, 3>{1, 2, 3});
521 0, x, std::array<Numeric, 4>{1, 2, 3, 4});
522
523 // Set up the interpolation weights
524 constexpr auto lin_iw = interpweights(lin);
525 constexpr auto sqr_iw = interpweights(sqr);
526 constexpr auto cub_iw = interpweights(cub);
527
528 // Get interpolation value
529 constexpr auto lin_lin = interp(y, lin_iw, lin);
530 constexpr auto lin_sqr = interp(y2, lin_iw, lin);
531 constexpr auto lin_cub = interp(y3, lin_iw, lin);
532 constexpr auto sqr_lin = interp(y, sqr_iw, sqr);
533 constexpr auto sqr_sqr = interp(y2, sqr_iw, sqr);
534 constexpr auto sqr_cub = interp(y3, sqr_iw, sqr);
535 constexpr auto cub_lin = interp(y, cub_iw, cub);
536 constexpr auto cub_sqr = interp(y2, cub_iw, cub);
537 constexpr auto cub_cub = interp(y3, cub_iw, cub);
538
539 // Compile-time check that these values are good
540 static_assert(is_around(lin_lin, x));
541 static_assert(is_around(sqr_sqr, x * x));
542 static_assert(is_around(cub_cub, x * x * x));
543
544 // Should output 1.5 2.5 4.5 1.5 2.25 3 1.5 2.25 3.375
545 std::cout << lin_lin << ' ' << lin_sqr << ' ' << lin_cub << ' ' << sqr_lin
546 << ' ' << sqr_sqr << ' ' << sqr_cub << ' ' << cub_lin << ' '
547 << cub_sqr << ' ' << cub_cub << '\n';
548}
549
550void test13() {
551 constexpr int N = 10;
552 constexpr std::array<Numeric, N> y{1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
553 constexpr std::array<Numeric, N> y2{1, 2 * 2, 3 * 3, 4 * 4, 5 * 5,
554 6 * 6, 7 * 7, 8 * 8, 9 * 9, 10 * 10};
555 constexpr std::array<Numeric, N> y3{
556 1, 2 * 2 * 2, 3 * 3 * 3, 4 * 4 * 4, 5 * 5 * 5,
557 6 * 6 * 6, 7 * 7 * 7, 8 * 8 * 8, 9 * 9 * 9, 10 * 10 * 10};
558 constexpr Numeric x = 2;
559 constexpr Numeric dx = 1e-2;
560
561 // Set up the interpolation Lagranges
562 constexpr FixedLagrangeInterpolation<1> lin(
563 0, x, std::array<Numeric, 2>{1, 2}, true);
564 constexpr FixedLagrangeInterpolation<2> sqr(
565 0, x, std::array<Numeric, 3>{1, 2, 3}, true);
566 constexpr FixedLagrangeInterpolation<3> cub(
567 0, x, std::array<Numeric, 4>{1, 2, 3, 4}, true);
568 constexpr FixedLagrangeInterpolation<1> dlin(
569 0, x + dx, std::array<Numeric, 2>{1, 2}, true);
570 constexpr FixedLagrangeInterpolation<2> dsqr(
571 0, x + dx, std::array<Numeric, 3>{1, 2, 3}, true);
572 constexpr FixedLagrangeInterpolation<3> dcub(
573 0, x + dx, std::array<Numeric, 4>{1, 2, 3, 4}, true);
574
575 // Set up the interpolation weights
576 constexpr auto lin_iw = interpweights(lin);
577 constexpr auto sqr_iw = interpweights(sqr);
578 constexpr auto cub_iw = interpweights(cub);
579 constexpr auto dlin_iw = dinterpweights(lin);
580 constexpr auto dsqr_iw = dinterpweights(sqr);
581 constexpr auto dcub_iw = dinterpweights(cub);
582 constexpr auto alt_lin_iw = interpweights(dlin);
583 constexpr auto alt_sqr_iw = interpweights(dsqr);
584 constexpr auto alt_cub_iw = interpweights(dcub);
585
586 // Get interpolation value
587 constexpr auto lin_lin = interp(y, lin_iw, lin);
588 constexpr auto sqr_sqr = interp(y2, sqr_iw, sqr);
589 constexpr auto cub_cub = interp(y3, cub_iw, cub);
590 constexpr auto dlin_lin = interp(y, dlin_iw, dlin);
591 constexpr auto dsqr_sqr = interp(y2, dsqr_iw, dsqr);
592 constexpr auto dcub_cub = interp(y3, dcub_iw, dcub);
593 constexpr auto alt_lin_lin = interp(y, alt_lin_iw, dlin);
594 constexpr auto alt_sqr_sqr = interp(y2, alt_sqr_iw, dsqr);
595 constexpr auto alt_cub_cub = interp(y3, alt_cub_iw, dcub);
596
597 // Compile-time check that these values are good
598 static_assert(is_around(lin_lin, x));
599 static_assert(is_around(sqr_sqr, x * x));
600 static_assert(is_around(cub_cub, x * x * x));
601 static_assert(is_around(dlin_lin, 1));
602 static_assert(is_around(dsqr_sqr, 2 * x));
603 static_assert(is_around(dcub_cub, 3 * x * x));
604
605 std::cout << lin_lin << ' ' << ' ' << sqr_sqr << ' ' << ' ' << ' ' << cub_cub
606 << '\n'
607 << dlin_lin << ' ' << ' ' << ' ' << ' ' << dsqr_sqr << ' ' << ' '
608 << ' ' << ' ' << ' ' << ' ' << dcub_cub << '\n'
609 << alt_lin_lin << ' ' << alt_sqr_sqr << ' ' << alt_cub_cub << '\n'
610 << (alt_lin_lin - lin_lin) / dx << ' ' << ' ' << ' ' << ' '
611 << (alt_sqr_sqr - sqr_sqr) / dx << ' ' << ' ' << ' '
612 << (alt_cub_cub - cub_cub) / dx << '\n';
613}
614
615void test14() {
616 Vector y{1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
617 Vector y2{1, 2 * 2, 3 * 3, 4 * 4, 5 * 5, 6 * 6, 7 * 7, 8 * 8, 9 * 9, 10 * 10};
618 Vector y3{1, 2 * 2 * 2, 3 * 3 * 3, 4 * 4 * 4, 5 * 5 * 5,
619 6 * 6 * 6, 7 * 7 * 7, 8 * 8 * 8, 9 * 9 * 9, 10 * 10 * 10};
620
621 Vector x{0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5};
622
623 // Set up the interpolation Lagranges
624 auto lin = Interpolation::FixedLagrangeVector<1>(
625 x, y, true, Interpolation::GridType::Standard);
626 auto sqr = Interpolation::FixedLagrangeVector<2>(
627 x, y, true, Interpolation::GridType::Standard);
628 auto cub = Interpolation::FixedLagrangeVector<3>(
629 x, y, true, Interpolation::GridType::Standard);
630
631 // Set up the interpolation weights
632 auto lin_iw = interpweights(lin);
633 auto sqr_iw = interpweights(sqr);
634 auto cub_iw = interpweights(cub);
635 auto dlin_iw = dinterpweights(lin);
636 auto dsqr_iw = dinterpweights(sqr);
637 auto dcub_iw = dinterpweights(cub);
638
639 // Get interpolation value
640 auto lin_lin = reinterp(y, lin_iw, lin);
641 auto sqr_sqr = reinterp(y2, sqr_iw, sqr);
642 auto cub_cub = reinterp(y3, cub_iw, cub);
643 auto dlin_lin = reinterp(y, dlin_iw, lin);
644 auto dsqr_sqr = reinterp(y2, dsqr_iw, sqr);
645 auto dcub_cub = reinterp(y3, dcub_iw, cub);
646
647 // Print the values and derivatives
648 std::cout << lin_lin << '\n' << sqr_sqr << '\n' << cub_cub << '\n';
649 std::cout << dlin_lin << '\n' << dsqr_sqr << '\n' << dcub_cub << '\n';
650}
651
652void test16() {
653 const Vector xasc{2, 3, 4};
654 const Vector xdes{4, 3, 2};
655 Vector x{2.25, 3.25, 2.35};
656
657 auto asc = Interpolation::LagrangeVector(x, xasc, 1, 1.0, true,
658 Interpolation::GridType::Standard);
659 auto des = Interpolation::LagrangeVector(x, xdes, 1, 1.0, true,
660 Interpolation::GridType::Standard);
661 auto fasc = Interpolation::FixedLagrangeVector<1>(
662 x, xasc, true, Interpolation::GridType::Standard);
663 auto fdes = Interpolation::FixedLagrangeVector<1>(
664 x, xdes, true, Interpolation::GridType::Standard);
665 for (Index i = 0; i < x.size(); i++) {
666 std::cout << x[i] << ": " << asc[i] << ' ' << '-' << ' ' << des[i] << ' '
667 << " xstart [asc des]: " << xasc[asc[i].pos] << ' '
668 << xdes[des[i].pos] << '\n';
669 }
670 std::cout << '\n';
671 for (Index i = 0; i < x.size(); i++) {
672 std::cout << x[i] << ": " << fasc[i] << ' ' << '-' << ' ' << fdes[i] << ' '
673 << " xstart [asc des]: " << xasc[fasc[i].pos] << ' '
674 << xdes[fdes[i].pos] << '\n';
675 }
676}
677
678void test17() {
679 const Index N = 500;
680 Vector x(N);
681 Verbosity verbosity;
682 VectorNLinSpace(x, N, 0, Constant::two_pi, verbosity);
683 Vector y = x;
684 for (auto& f : y) f = std::sin(f);
685 for (Numeric n = -Constant::two_pi; n <= 2 * Constant::two_pi; n += 0.1) {
686 auto lag = Interpolation::Lagrange(0, n, x, 1, true,
687 Interpolation::GridType::Cyclic,
688 {0, Constant::two_pi});
690 0, n, x, true, Interpolation::GridType::Cyclic,
691 {0, Constant::two_pi});
692 auto lag_iw = interpweights(lag);
693 auto flag_iw = interpweights(flag);
694 auto dlag_iw = dinterpweights(lag);
695 auto dflag_iw = dinterpweights(flag);
696 std::cout << n << ' ' << interp(y, lag_iw, lag) << ' '
697 << interp(y, flag_iw, flag) << ' ' << std::sin(n) << ' '
698 << interp(y, dlag_iw, lag) << ' ' << interp(y, dflag_iw, flag)
699 << ' ' << std::cos(n) << '\n';
700 }
701}
702
703void test18() {
704 const Index N = 500;
705 Vector x(N);
706 Verbosity verbosity;
707 VectorNLinSpace(x, N, -180, 180, verbosity);
708 Vector y = x;
709 for (auto& f : y) f = Conversion::sind(f);
710 for (Numeric n = -3 * 180; n <= 3 * 180; n += 0.1) {
711 auto lag = Interpolation::Lagrange(0, n, x, 1, true,
712 Interpolation::GridType::Cyclic,
713 {-180, 180});
715 0, n, x, true, Interpolation::GridType::Cyclic, {-180, 180});
716 auto lag_iw = interpweights(lag);
717 auto flag_iw = interpweights(flag);
718 auto dlag_iw = dinterpweights(lag);
719 auto dflag_iw = dinterpweights(flag);
720 std::cout << n << ' ' << interp(y, lag_iw, lag) << ' '
721 << interp(y, flag_iw, flag) << ' ' << Conversion::sind(n) << ' '
722 << interp(y, dlag_iw, lag) << ' ' << interp(y, dflag_iw, flag)
723 << ' ' << Conversion::cosd(n) << '\n';
724 }
725}
726
727void test19() {
728 const Index N = 500;
729 Vector x(N);
730 Verbosity verbosity;
731 VectorNLinSpace(x, N, 0, 0.5, verbosity);
732 Vector y = x;
733 for (auto& f : y) f = Conversion::sind(720 * f);
734 for (Numeric n = -0.5; n <= 1.5; n += 0.01) {
735 auto lag =
736 Interpolation::Lagrange(0, n, x, 1, true,
737 Interpolation::GridType::Cyclic, {0, 0.5});
739 0, n, x, true, Interpolation::GridType::Cyclic, {0, 0.5});
740 auto lag_iw = interpweights(lag);
741 auto flag_iw = interpweights(flag);
742 auto dlag_iw = dinterpweights(lag);
743 auto dflag_iw = dinterpweights(flag);
744 std::cout << n << ' ' << interp(y, lag_iw, lag) << ' '
745 << interp(y, flag_iw, flag) << ' ' << Conversion::sind(720 * n)
746 << ' ' << interp(y, dlag_iw, lag) << ' '
747 << interp(y, dflag_iw, flag) << ' ' << Conversion::cosd(720 * n)
748 << '\n';
749 }
750}
751
752void test20() {
753 const Index N = 500;
754 Vector x(N);
755 Verbosity verbosity;
756 VectorNLinSpace(x, N, -0.123, 0.456, verbosity);
757 Vector y = x;
758 for (auto& f : y) f = Conversion::sind(360 / (0.456 + 0.123) * f);
759 for (Numeric n = -0.5; n <= 1.5; n += 0.01) {
760 auto lag = Interpolation::Lagrange(0, n, x, 1, true,
761 Interpolation::GridType::Cyclic,
762 {-0.123, 0.456});
764 0, n, x, true, Interpolation::GridType::Cyclic, {-0.123, 0.456});
765 auto lag_iw = interpweights(lag);
766 auto flag_iw = interpweights(flag);
767 auto dlag_iw = dinterpweights(lag);
768 auto dflag_iw = dinterpweights(flag);
769 std::cout << n << ' ' << interp(y, lag_iw, lag) << ' '
770 << interp(y, flag_iw, flag) << ' '
771 << Conversion::sind(360 / (0.456 + 0.123) * n) << ' '
772 << interp(y, dlag_iw, lag) << ' ' << interp(y, dflag_iw, flag)
773 << ' ' << Conversion::cosd(360 / (0.456 + 0.123) * n) << '\n';
774 }
775}
776
777void test21() {
778 const Index N = 500;
779 Vector x(N);
780 Verbosity verbosity;
781 VectorNLinSpace(x, N, 0.05, 0.45, verbosity);
782 Vector y = x;
783 for (auto& f : y) f = Conversion::sind(720 * f);
784 for (Numeric n = -0.5; n <= 1.5; n += 0.01) {
785 auto lag =
786 Interpolation::Lagrange(0, n, x, 1, true,
787 Interpolation::GridType::Cyclic, {0, 0.5});
789 0, n, x, true, Interpolation::GridType::Cyclic, {0, 0.5});
790 auto lag_iw = interpweights(lag);
791 auto flag_iw = interpweights(flag);
792 auto dlag_iw = dinterpweights(lag);
793 auto dflag_iw = dinterpweights(flag);
794 std::cout << n << ' ' << interp(y, lag_iw, lag) << ' '
795 << interp(y, flag_iw, flag) << ' ' << Conversion::sind(720 * n)
796 << ' ' << interp(y, dlag_iw, lag) << ' '
797 << interp(y, dflag_iw, flag) << ' ' << Conversion::cosd(720 * n)
798 << '\n';
799 }
800}
801
802void test22() {
803 const Index N = 500;
804 Vector x(N);
805 Verbosity verbosity;
806 VectorNLinSpace(x, N, Constant::two_pi, 0, verbosity);
807 Vector y = x;
808 for (auto& f : y) f = std::sin(f);
809 for (Numeric n = -Constant::two_pi; n <= 2 * Constant::two_pi; n += 0.1) {
810 auto lag = Interpolation::Lagrange(0, n, x, 1, true,
811 Interpolation::GridType::Cyclic,
812 {0, Constant::two_pi});
814 0, n, x, true, Interpolation::GridType::Cyclic,
815 {0, Constant::two_pi});
816 auto lag_iw = interpweights(lag);
817 auto flag_iw = interpweights(flag);
818 auto dlag_iw = dinterpweights(lag);
819 auto dflag_iw = dinterpweights(flag);
820 std::cout << n << ' ' << interp(y, lag_iw, lag) << ' '
821 << interp(y, flag_iw, flag) << ' ' << std::sin(n) << ' '
822 << interp(y, dlag_iw, lag) << ' ' << interp(y, dflag_iw, flag)
823 << ' ' << std::cos(n) << '\n';
824 }
825}
826
827void test23() {
828 const Index N = 500;
829 Vector x(N);
830 Verbosity verbosity;
831 VectorNLinSpace(x, N, 0.45, 0.05, verbosity);
832 Vector y = x;
833 for (auto& f : y) f = Conversion::sind(720 * f);
834 for (Numeric n = -0.5; n <= 1.5; n += 0.01) {
835 auto lag =
836 Interpolation::Lagrange(0, n, x, 1, true,
837 Interpolation::GridType::Cyclic, {0, 0.5});
839 0, n, x, true, Interpolation::GridType::Cyclic, {0, 0.5});
840 auto lag_iw = interpweights(lag);
841 auto flag_iw = interpweights(flag);
842 auto dlag_iw = dinterpweights(lag);
843 auto dflag_iw = dinterpweights(flag);
844 std::cout << n << ' ' << interp(y, lag_iw, lag) << ' '
845 << interp(y, flag_iw, flag) << ' ' << Conversion::sind(720 * n)
846 << ' ' << interp(y, dlag_iw, lag) << ' '
847 << interp(y, dflag_iw, flag) << ' ' << Conversion::cosd(720 * n)
848 << '\n';
849 }
850}
851
852void test25() {
853 const Index N = 500;
854 Vector x(N);
855 Verbosity verbosity;
856 VectorNLinSpace(x, N, 30, 150, verbosity);
857 Vector y = x;
858 for (auto& f : y) f = 15*f*f + f*f*f;
859 for (Numeric n = 0; n <= 180; n += 0.01) {
860 const auto lag =
861 Interpolation::Lagrange(0, n, x, 5, true, Interpolation::GridType::CosDeg);
862 std::cout << n << ' ' << interp(y, interpweights(lag), lag) << ' ' << interp(y, dinterpweights(lag), lag) << '\n';
863 }
864}
865
866void test26() {
867 auto f = [](Numeric x){return Math::pow2(Interpolation::cyclic_clamp(x, {-2, 2})) - 4;};
868 auto df = [](Numeric x){return 2*Interpolation::cyclic_clamp(x, {-2, 2});};
869
870 constexpr Index N = 9;
871 constexpr std::array<Numeric, N> xi{-2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2};
872 constexpr std::array<Numeric, N> yi{f(-2), f(-1.5), f(-1), f(-0.5), f(0), f(0.5), f(1), f(1.5), f(2)};
873 constexpr Index O1 = 3;
874
875 // Test for a few values of interpolation
876 {
877 constexpr Numeric x = -1.75;
878 constexpr Interpolation::FixedLagrange<O1> cyc(0, x, xi, true, Interpolation::GridType::Cyclic, {-2, 2});
879 static_assert(f(x) == interp(yi, interpweights(cyc), cyc));
880 static_assert(df(x) == interp(yi, dinterpweights(cyc), cyc));
881 constexpr Interpolation::FixedLagrange<O1> lin(0, x, xi, true, Interpolation::GridType::Standard);
882 static_assert(f(x) == interp(yi, interpweights(lin), lin));
883 static_assert(df(x) == interp(yi, dinterpweights(lin), lin));
884 }
885 {
886 constexpr Numeric x = -1.25;
887 constexpr Interpolation::FixedLagrange<O1> cyc(0, x, xi, true, Interpolation::GridType::Cyclic, {-2, 2});
888 static_assert(f(x) == interp(yi, interpweights(cyc), cyc));
889 static_assert(df(x) == interp(yi, dinterpweights(cyc), cyc));
890 constexpr Interpolation::FixedLagrange<O1> lin(0, x, xi, true, Interpolation::GridType::Standard);
891 static_assert(f(x) == interp(yi, interpweights(lin), lin));
892 static_assert(df(x) == interp(yi, dinterpweights(lin), lin));
893 }
894 {
895 constexpr Numeric x = -0.25;
896 constexpr Interpolation::FixedLagrange<O1> cyc(0, x, xi, true, Interpolation::GridType::Cyclic, {-2, 2});
897 static_assert(f(x) == interp(yi, interpweights(cyc), cyc));
898 static_assert(df(x) == interp(yi, dinterpweights(cyc), cyc));
899 constexpr Interpolation::FixedLagrange<O1> lin(0, x, xi, true, Interpolation::GridType::Standard);
900 static_assert(f(x) == interp(yi, interpweights(lin), lin));
901 static_assert(df(x) == interp(yi, dinterpweights(lin), lin));
902 }
903 {
904 constexpr Numeric x = 1;
905 constexpr Interpolation::FixedLagrange<O1> cyc(0, x, xi, true, Interpolation::GridType::Cyclic, {-2, 2});
906 static_assert(f(x) == interp(yi, interpweights(cyc), cyc));
907 static_assert(df(x) == interp(yi, dinterpweights(cyc), cyc));
908 constexpr Interpolation::FixedLagrange<O1> lin(0, x, xi, true, Interpolation::GridType::Standard);
909 static_assert(f(x) == interp(yi, interpweights(lin), lin));
910 static_assert(df(x) == interp(yi, dinterpweights(lin), lin));
911 }
912 {
913 constexpr Numeric x = -2;
914 constexpr Interpolation::FixedLagrange<O1> cyc(0, x, xi, true, Interpolation::GridType::Cyclic, {-2, 2});
915 static_assert(f(x) == interp(yi, interpweights(cyc), cyc));
916 static_assert(df(x) == interp(yi, dinterpweights(cyc), cyc));
917 constexpr Interpolation::FixedLagrange<O1> lin(0, x, xi, true, Interpolation::GridType::Standard);
918 static_assert(f(x) == interp(yi, interpweights(lin), lin));
919 static_assert(df(x) == interp(yi, dinterpweights(lin), lin));
920 }
921 {
922 constexpr Numeric x = 0;
923 constexpr Interpolation::FixedLagrange<O1> cyc(0, x, xi, true, Interpolation::GridType::Cyclic, {-2, 2});
924 static_assert(f(x) == interp(yi, interpweights(cyc), cyc));
925 static_assert(df(x) == interp(yi, dinterpweights(cyc), cyc));
926 constexpr Interpolation::FixedLagrange<O1> lin(0, x, xi, true, Interpolation::GridType::Standard);
927 static_assert(f(x) == interp(yi, interpweights(lin), lin));
928 static_assert(df(x) == interp(yi, dinterpweights(lin), lin));
929 }
930 {
931 constexpr Numeric x = -4;
932 constexpr Interpolation::FixedLagrange<O1> cyc(0, x, xi, true, Interpolation::GridType::Cyclic, {-2, 2});
933 static_assert(f(x) == interp(yi, interpweights(cyc), cyc));
934 static_assert(df(x) == interp(yi, dinterpweights(cyc), cyc));
935 }
936 {
937 constexpr Numeric x = 4;
938 constexpr Interpolation::FixedLagrange<O1> cyc(0, x, xi, true, Interpolation::GridType::Cyclic, {-2, 2});
939 static_assert(f(x) == interp(yi, interpweights(cyc), cyc));
940 static_assert(df(x) == interp(yi, dinterpweights(cyc), cyc));
941 }
942
943 constexpr Index O2 = 3;
944 std::cout << "x f(x) interp(x) df(x) dinterp(x)\n";
945 for (Numeric X=-3; X<3; X+=0.025) {
946 const Interpolation::FixedLagrange<O2> cyc(0, X, xi, true, Interpolation::GridType::Cyclic, {-2, 2});
947 const Interpolation::FixedLagrange<O2> lin(0, X, xi, true, Interpolation::GridType::Standard);
948 std::cout << X << ' ' << f(X) << ' ' << interp(yi, interpweights(cyc), cyc) << ' ' << df(X) << ' ' << interp(yi, dinterpweights(cyc), cyc)
949 << ' ' << interp(yi, interpweights(lin), lin) << ' ' << interp(yi, dinterpweights(lin), lin) << '\n';
950 }
951}
952
953void test27() {
954 static_assert(Interpolation::toGridType("Cyclic") == Interpolation::GridType::Cyclic);
955 static_assert(Interpolation::toString(Interpolation::GridType::Cyclic) == "Cyclic");
956 for (auto a : Interpolation::enumstrs::GridTypeNames)
957 std::cout << a.size() << '\n';
958 for (auto a : Interpolation::enumstrs::GridTypeNames)
959 std::cout << a << '\n';
960 for (auto a : Interpolation::enumstrs::GridTypeNames)
961 std::cout << a << '\n';
962}
963
964void test28() {
965 using Math::pow2;
966 using Math::pow3;
967 using Conversion::sind;
968 using Conversion::cosd;
969
970 // Old Grid of pressure, latitude, and longitude
971 Vector pre; VectorNLogSpace(pre, 10, 1e5, 1e-1, Verbosity());
972 Vector lat; VectorNLinSpace(lat, 5, -80, 80, Verbosity());
973 Vector lon; VectorNLinSpace(lon, 4, -170, 170, Verbosity());
974
975 // New Grids (pressure will be reduced)
976 Vector newpre(1, pre[pre.nelem()/2]);
977 Vector newlat; VectorNLinSpace(newlat, 4, - 90, 90, Verbosity());
978 Vector newlon; VectorNLinSpace(newlon, 3, -180, 180, Verbosity());
979
980 // Old Data given some values
981 Tensor3 data(pre.nelem(), lat.nelem(), lon.nelem());
982 for (Index i=0; i<pre.nelem(); i++) {
983 for (Index j=0; j<lat.nelem(); j++) {
984 for (Index k=0; k<lon.nelem(); k++) {
985 data(i, j, k) =
986 std::log(pre[i]) * sind(lat[j]) * pow3(cosd(lon[k]));
987 }
988 }
989 }
990
991 // Create Lagrange interpolation values
992 const ArrayOfLagrangeInterpolation lag_pre =
993 Interpolation::LagrangeVector(newpre, pre, 2, 1e9, false,
994 Interpolation::GridType::Log);
995 const ArrayOfLagrangeInterpolation lag_lat =
996 Interpolation::LagrangeVector(newlat, lat, 3, 1e9, false,
997 Interpolation::GridType::SinDeg);
998 const ArrayOfLagrangeInterpolation lag_lon =
999 Interpolation::LagrangeVector(newlon, lon, 1, 1e9, false,
1000 Interpolation::GridType::Cyclic, {-180, 180});
1001
1002 // Create the interpolation weights
1003 const auto lag_iw = interpweights(lag_pre, lag_lat, lag_lon);
1004
1005 // Create and reduce the new data to dim 1 and 2
1006 const Matrix newdata = reinterp(data,
1007 lag_iw, lag_pre, lag_lat, lag_lon).reduce_rank<1, 2>();
1008
1009 // Print the data
1010 std::cout << data(pre.nelem()/2, joker, joker)
1011 << '\n' << '\n' << newdata << '\n';
1012}
1013
1014int main() { test09(); }
This file contains the definition of Array.
Common ARTS conversions.
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Matrix.
Definition: matpackI.h:1043
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:140
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:146
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The iterator class for sub vectors.
Definition: matpackI.h:390
The MatrixView class.
Definition: matpackI.h:1164
The Matrix class.
Definition: matpackI.h:1261
The range class.
Definition: matpackI.h:159
The Tensor3 class.
Definition: matpackIII.h:346
The Tensor6 class.
Definition: matpackVI.h:1099
The Tensor7 class.
Definition: matpackVII.h:2399
The Vector class.
Definition: matpackI.h:899
constexpr std::string_view toString(EnergyLevelMapType x) noexcept
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Header file for interpolation.cc.
long int flag
#define dx
void VectorNLogSpace(Vector &x, const Index &n, const Numeric &start, const Numeric &stop, const Verbosity &verbosity)
WORKSPACE METHOD: VectorNLogSpace.
void VectorNLinSpace(Vector &x, const Index &n, const Numeric &start, const Numeric &stop, const Verbosity &verbosity)
WORKSPACE METHOD: VectorNLinSpace.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:227
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
const Joker joker
constexpr Numeric two_pi
Two times pi.
constexpr Numeric k
Boltzmann constant convenience name [J/K].
constexpr Numeric e
Elementary charge convenience name [C].
auto sind(auto x) noexcept
Returns the sine of the deg2rad of the input.
auto cosd(auto x) noexcept
Returns the cosine of the deg2rad of the input.
constexpr Numeric cyclic_clamp(const Numeric x, const std::pair< Numeric, Numeric > xlim) noexcept
void reinterp(VectorView out, const ConstVectorView &iy, const Grid< Vector, 1 > &iw, const Array< Lagrange > &dim0)
Vector dinterpweights(const Lagrange &dim0)
Array< Lagrange > LagrangeVector(const ConstVectorView &xs, const ConstVectorView &xi, const Index polyorder, const Numeric extrapol, const bool do_derivs, const GridType type, const std::pair< Numeric, Numeric > cycle)
constexpr auto pow3(auto x) noexcept
power of three
constexpr auto pow2(auto x) noexcept
power of two
#define N
Definition: rng.cc:164
Structure to store a grid position.
Definition: interpolation.h:73
A Lagrange interpolation computer.
void test26()
void test05()
void test18()
void test04()
void test20()
void test27()
void test03(Index n)
void test25()
void test02(Index n)
void test23()
void test06()
void test16()
void test12()
void test19()
void test17()
constexpr bool is_around(Numeric x, Numeric x0, Numeric e=1e-12)
void test13()
void test09()
void test14()
void test01()
void test11()
int main()
void test28()
void test22()
void test21()
#define a
This file contains basic functions to handle XML data files.