ARTS  2.4.0(git:4fb77825)
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 #include "array.h"
21 #include "interpolation.h"
22 #include "interpolation_poly.h"
23 #include "math_funcs.h"
24 #include "matpackVII.h"
25 #include "xml_io.h"
26 
27 void test01() {
28  cout << "Simple interpolation cases\n"
29  << "--------------------------\n";
30  // Vector og(5,5,-1); // 5,4,3,2,1
31  Vector og(1, 5, +1); // 1, 2, 3, 4, 5
32  Vector ng(2, 5, 0.25); // 2.0, 2,25, 2.5, 2.75, 3.0
33 
34  cout << "og:\n" << og << "\n";
35  cout << "ng:\n" << ng << "\n";
36 
37  // To store the grid positions:
38  ArrayOfGridPos gp(ng.nelem());
39 
40  gridpos(gp, og, ng);
41  cout << "gp:\n" << gp << "\n";
42 
43  cout << "1D:\n"
44  << "---\n";
45  {
46  // To store interpolation weights:
47  Matrix itw(gp.nelem(), 2);
48  interpweights(itw, gp);
49 
50  cout << "itw:\n" << itw << "\n";
51 
52  // Original field:
53  Vector of(og.nelem(), 0);
54  of[2] = 10; // 0, 0, 10, 0, 0
55 
56  cout << "of:\n" << of << "\n";
57 
58  // Interpolated field:
59  Vector nf(ng.nelem());
60 
61  interp(nf, itw, of, gp);
62 
63  cout << "nf:\n" << nf << "\n";
64  }
65 
66  cout << "Blue 2D:\n"
67  << "--------\n";
68  {
69  // To store interpolation weights:
70  Matrix itw(gp.nelem(), 4);
71  interpweights(itw, gp, gp);
72 
73  cout << "itw:\n" << itw << "\n";
74 
75  // Original field:
76  Matrix of(og.nelem(), og.nelem(), 0);
77  of(2, 2) = 10; // 0 Matrix with 10 in the middle
78 
79  cout << "of:\n" << of << "\n";
80 
81  // Interpolated field:
82  Vector nf(ng.nelem());
83 
84  interp(nf, itw, of, gp, gp);
85 
86  cout << "nf:\n" << nf << "\n";
87  }
88 
89  cout << "Blue 6D:\n"
90  << "--------\n";
91  {
92  // To store interpolation weights:
93  Matrix itw(gp.nelem(), 64);
94  interpweights(itw, gp, gp, gp, gp, gp, gp);
95 
96  // cout << "itw:\n" << itw << "\n";
97 
98  // Original field:
99  Index n = og.nelem();
100  Tensor6 of(n, n, n, n, n, n, 0);
101  of(2, 2, 2, 2, 2, 2) = 10; // 0 Tensor with 10 in the middle
102 
103  // cout << "of:\n" << of << "\n";
104 
105  // Interpolated field:
106  Vector nf(ng.nelem());
107 
108  interp(nf, itw, of, gp, gp, gp, gp, gp, gp);
109 
110  cout << "nf:\n" << nf << "\n";
111  }
112 
113  cout << "Green 2D:\n"
114  << "---------\n";
115  {
116  // To store interpolation weights:
117  Tensor3 itw(gp.nelem(), gp.nelem(), 4);
118  interpweights(itw, gp, gp);
119 
120  for (Index i = 0; i < itw.ncols(); ++i)
121  cout << "itw " << i << ":\n"
122  << itw(Range(joker), Range(joker), i) << "\n";
123 
124  // Original field:
125  Matrix of(og.nelem(), og.nelem(), 0);
126  of(2, 2) = 10; // 0 Matrix with 10 in the middle
127 
128  cout << "of:\n" << of << "\n";
129 
130  // Interpolated field:
131  Matrix nf(ng.nelem(), ng.nelem());
132 
133  interp(nf, itw, of, gp, gp);
134 
135  cout << "nf:\n" << nf << "\n";
136  }
137 
138  cout << "Green 6D:\n"
139  << "---------\n";
140  {
141  // To store interpolation weights:
142  Tensor7 itw(gp.nelem(),
143  gp.nelem(),
144  gp.nelem(),
145  gp.nelem(),
146  gp.nelem(),
147  gp.nelem(),
148  64);
149  interpweights(itw, gp, gp, gp, gp, gp, gp);
150 
151  // Original field:
152  Tensor6 of(og.nelem(),
153  og.nelem(),
154  og.nelem(),
155  og.nelem(),
156  og.nelem(),
157  og.nelem(),
158  0);
159  of(2, 2, 2, 2, 2, 2) = 10; // 0 Tensor with 10 in the middle
160 
161  cout << "Middle slice of of:\n"
162  << of(2, 2, 2, 2, Range(joker), Range(joker)) << "\n";
163 
164  // Interpolated field:
165  Tensor6 nf(
166  ng.nelem(), ng.nelem(), ng.nelem(), ng.nelem(), ng.nelem(), ng.nelem());
167 
168  interp(nf, itw, of, gp, gp, gp, gp, gp, gp);
169 
170  cout << "Last slice of nf:\n"
171  << nf(4, 4, 4, 4, Range(joker), Range(joker)) << "\n";
172  }
173 }
174 
175 void test02(Index n) {
176  cout << "Test whether for loop or iterator are quicker\n"
177  << "a) for loop\n"
178  << "---------------------------------------------\n";
179 
180  Vector a(n);
181  for (Index i = 0; i < a.nelem(); ++i) a[i] = (Numeric)i;
182 }
183 
184 void test03(Index n) {
185  cout << "Test whether for loop or iterator are quicker\n"
186  << "b) iterator\n"
187  << "---------------------------------------------\n";
188 
189  Vector a(n);
190  Iterator1D ai = a.begin();
191  const Iterator1D ae = a.end();
192  Index i = 0;
193  for (; ai != ae; ++ai, ++i) *ai = (Numeric)i;
194 }
195 
196 // Result: Both are almost equally fast, with a slight advantage of
197 // the for loop if compiler optimization is enabled.
198 
199 void test04() {
200  cout << "Green type interpolation of all pages of a Tensor3\n";
201 
202  // The original Tensor is called a, the new one n.
203 
204  // 10 pages, 20 rows, 30 columns, all grids are: 1,2,3
205  Vector a_pgrid(1, 3, 1), a_rgrid(1, 3, 1), a_cgrid(1, 3, 1);
206  Tensor3 a(a_pgrid.nelem(), a_rgrid.nelem(), a_cgrid.nelem());
207 
208  a = 0;
209  // Put some simple numbers in the middle of each page:
210  a(0, 1, 1) = 10;
211  a(1, 1, 1) = 20;
212  a(2, 1, 1) = 30;
213 
214  // New row and column grids:
215  // 1, 1.5, 2, 2.5, 3
216  Vector n_rgrid(1, 5, .5), n_cgrid(1, 5, .5);
217  Tensor3 n(a_pgrid.nelem(), n_rgrid.nelem(), n_cgrid.nelem());
218 
219  // So, n has the same number of pages as a, but more rows and columns.
220 
221  // Get the grid position arrays:
222  ArrayOfGridPos n_rgp(n_rgrid.nelem()); // For row grid positions.
223  ArrayOfGridPos n_cgp(n_cgrid.nelem()); // For column grid positions.
224 
225  gridpos(n_rgp, a_rgrid, n_rgrid);
226  gridpos(n_cgp, a_cgrid, n_cgrid);
227 
228  // Get the interpolation weights:
229  Tensor3 itw(n_rgrid.nelem(), n_cgrid.nelem(), 4);
230  interpweights(itw, n_rgp, n_cgp);
231 
232  // Do a "green" interpolation for all pages of the Tensor a:
233 
234  for (Index i = 0; i < a.npages(); ++i) {
235  // Select the current page of both a and n:
236  ConstMatrixView ap = a(i, Range(joker), Range(joker));
237  MatrixView np = n(i, Range(joker), Range(joker));
238 
239  // Do the interpolation:
240  interp(np, itw, ap, n_rgp, n_cgp);
241 
242  // Note that this is efficient, because interpolation weights and
243  // grid positions are re-used.
244  }
245 
246  cout << "Original field:\n";
247  for (Index i = 0; i < a.npages(); ++i)
248  cout << "page " << i << ":\n" << a(i, Range(joker), Range(joker)) << "\n";
249 
250  cout << "Interpolated field:\n";
251  for (Index i = 0; i < n.npages(); ++i)
252  cout << "page " << i << ":\n" << n(i, Range(joker), Range(joker)) << "\n";
253 }
254 
255 void test05() {
256  cout << "Very simple interpolation case\n";
257 
258  Vector og(1, 5, +1); // 1, 2, 3, 4, 5
259  Vector ng(2, 5, 0.25); // 2.0, 2,25, 2.5, 2.75, 3.0
260 
261  cout << "Original grid:\n" << og << "\n";
262  cout << "New grid:\n" << ng << "\n";
263 
264  // To store the grid positions:
265  ArrayOfGridPos gp(ng.nelem());
266 
267  gridpos(gp, og, ng);
268  cout << "Grid positions:\n" << gp;
269 
270  // To store interpolation weights:
271  Matrix itw(gp.nelem(), 2);
272  interpweights(itw, gp);
273 
274  cout << "Interpolation weights:\n" << itw << "\n";
275 
276  // Original field:
277  Vector of(og.nelem(), 0);
278  of[2] = 10; // 0, 0, 10, 0, 0
279 
280  cout << "Original field:\n" << of << "\n";
281 
282  // Interpolated field:
283  Vector nf(ng.nelem());
284 
285  interp(nf, itw, of, gp);
286 
287  cout << "New field:\n" << nf << "\n";
288 }
289 
290 void test06() {
291  cout << "Simple extrapolation cases\n"
292  << "--------------------------\n";
293  // Vector og(5,5,-1); // 5,4,3,2,1
294  Vector og(1, 5, +1); // 1, 2, 3, 4, 5
295  Vector ng{0.9, 1.5, 3, 4.5, 5.1}; // 0.9, 1.5, 3, 4.5, 5.1
296 
297  cout << "og:\n" << og << "\n";
298  cout << "ng:\n" << ng << "\n";
299 
300  // To store the grid positions:
301  ArrayOfGridPos gp(ng.nelem());
302 
303  gridpos(gp, og, ng);
304  cout << "gp:\n" << gp << "\n";
305 
306  cout << "1D:\n"
307  << "---\n";
308  {
309  // To store interpolation weights:
310  Matrix itw(gp.nelem(), 2);
311  interpweights(itw, gp);
312 
313  cout << "itw:\n" << itw << "\n";
314 
315  // Original field:
316  Vector of(og.nelem(), 0);
317  for (Index i = 0; i < og.nelem(); ++i)
318  of[i] = (Numeric)(10 * (i + 1)); // 10, 20, 30, 40, 50
319 
320  cout << "of:\n" << of << "\n";
321 
322  // Interpolated field:
323  Vector nf(ng.nelem());
324 
325  interp(nf, itw, of, gp);
326 
327  cout << "nf:\n" << nf << "\n";
328  }
329 
330  cout << "Blue 2D:\n"
331  << "--------\n";
332  {
333  // To store interpolation weights:
334  Matrix itw(gp.nelem(), 4);
335  interpweights(itw, gp, gp);
336 
337  cout << "itw:\n" << itw << "\n";
338 
339  // Original field:
340  Matrix of(og.nelem(), og.nelem(), 0);
341  of(2, 2) = 10; // 0 Matrix with 10 in the middle
342 
343  cout << "of:\n" << of << "\n";
344 
345  // Interpolated field:
346  Vector nf(ng.nelem());
347 
348  interp(nf, itw, of, gp, gp);
349 
350  cout << "nf:\n" << nf << "\n";
351  }
352 
353  cout << "Blue 6D:\n"
354  << "--------\n";
355  {
356  // To store interpolation weights:
357  Matrix itw(gp.nelem(), 64);
358  interpweights(itw, gp, gp, gp, gp, gp, gp);
359 
360  // cout << "itw:\n" << itw << "\n";
361 
362  // Original field:
363  Index n = og.nelem();
364  Tensor6 of(n, n, n, n, n, n, 0);
365  of(2, 2, 2, 2, 2, 2) = 10; // 0 Tensor with 10 in the middle
366 
367  // cout << "of:\n" << of << "\n";
368 
369  // Interpolated field:
370  Vector nf(ng.nelem());
371 
372  interp(nf, itw, of, gp, gp, gp, gp, gp, gp);
373 
374  cout << "nf:\n" << nf << "\n";
375  }
376 
377  cout << "Green 2D:\n"
378  << "---------\n";
379  {
380  // To store interpolation weights:
381  Tensor3 itw(gp.nelem(), gp.nelem(), 4);
382  interpweights(itw, gp, gp);
383 
384  for (Index i = 0; i < itw.ncols(); ++i)
385  cout << "itw " << i << ":\n"
386  << itw(Range(joker), Range(joker), i) << "\n";
387 
388  // Original field:
389  Matrix of(og.nelem(), og.nelem(), 0);
390  of(2, 2) = 10; // 0 Matrix with 10 in the middle
391 
392  cout << "of:\n" << of << "\n";
393 
394  // Interpolated field:
395  Matrix nf(ng.nelem(), ng.nelem());
396 
397  interp(nf, itw, of, gp, gp);
398 
399  cout << "nf:\n" << nf << "\n";
400  }
401 
402  cout << "Green 6D:\n"
403  << "---------\n";
404  {
405  // To store interpolation weights:
406  Tensor7 itw(gp.nelem(),
407  gp.nelem(),
408  gp.nelem(),
409  gp.nelem(),
410  gp.nelem(),
411  gp.nelem(),
412  64);
413  interpweights(itw, gp, gp, gp, gp, gp, gp);
414 
415  // Original field:
416  Tensor6 of(og.nelem(),
417  og.nelem(),
418  og.nelem(),
419  og.nelem(),
420  og.nelem(),
421  og.nelem(),
422  0);
423  of(2, 2, 2, 2, 2, 2) = 10; // 0 Tensor with 10 in the middle
424 
425  cout << "Middle slice of of:\n"
426  << of(2, 2, 2, 2, Range(joker), Range(joker)) << "\n";
427 
428  // Interpolated field:
429  Tensor6 nf(
430  ng.nelem(), ng.nelem(), ng.nelem(), ng.nelem(), ng.nelem(), ng.nelem());
431 
432  interp(nf, itw, of, gp, gp, gp, gp, gp, gp);
433 
434  cout << "Last slice of nf:\n"
435  << nf(4, 4, 4, 4, Range(joker), Range(joker)) << "\n";
436  }
437 }
438 
439 // Test polynomial interpolation (included by CE)
440 void test07() {
441  // FileType ftype = FILE_TYPE_ASCII;
442  Vector new_x(0, 21, +0.25);
443  Vector x(0, 10, +1);
444 
445  ArrayOfGridPos gp(new_x.nelem());
446  gridpos(gp, x, new_x);
447 
448  Vector y1(x.nelem());
449  Vector y2(x.nelem());
450  Vector y3(x.nelem());
451 
452  for (Index i = 0; i < x.nelem(); i++) {
453  // linear function
454  y1[i] = 3 * x[i];
455  // cubic function
456  y2[i] = pow(x[i], 3) + 2;
457  // trigonometric function
458  y3[i] = sin(x[i]);
459  }
460 
461  // Linear interpolation:
462  Matrix itw(gp.nelem(), 2);
463  interpweights(itw, gp);
464 
465  Vector y1_lin(new_x.nelem());
466  Vector y2_lin(new_x.nelem());
467  Vector y3_lin(new_x.nelem());
468 
469  interp(y1_lin, itw, y1, gp);
470  interp(y2_lin, itw, y2, gp);
471  interp(y3_lin, itw, y3, gp);
472 
473  cout << "y1_lin = [" << y1_lin << "];\n";
474  cout << "y2_lin = [" << y2_lin << "];\n";
475  cout << "y3_lin = [" << y3_lin << "];\n";
476 
477  // Cubic interpolation:
478  Vector y1_cub(new_x.nelem());
479  Vector y2_cub(new_x.nelem());
480  Vector y3_cub(new_x.nelem());
481 
482  for (Index i = 0; i < new_x.nelem(); i++) {
483  y1_cub[i] = interp_poly(x, y1, new_x[i], gp[i]);
484  y2_cub[i] = interp_poly(x, y2, new_x[i], gp[i]);
485  y3_cub[i] = interp_poly(x, y3, new_x[i], gp[i]);
486  }
487 
488  cout << "y1_cub = [" << y1_cub << "];\n";
489  cout << "y2_cub = [" << y2_cub << "];\n";
490  cout << "y3_cub = [" << y3_cub << "];\n";
491 
492  // Stefan's new polynomial interpolation routines:
493  Index order = 2;
494 
495  Vector y1_new(new_x.nelem());
496  Vector y2_new(new_x.nelem());
497  Vector y3_new(new_x.nelem());
498 
499  ArrayOfGridPosPoly gpp(new_x.nelem());
500  gridpos_poly(gpp, x, new_x, order);
501  Matrix itwp(new_x.nelem(), order + 1);
502  interpweights(itwp, gpp);
503 
504  interp(y1_new, itwp, y1, gpp);
505  interp(y2_new, itwp, y2, gpp);
506  interp(y3_new, itwp, y3, gpp);
507 
508  cout << "y1_new = [" << y1_new << "];\n";
509  cout << "y2_new = [" << y2_new << "];\n";
510  cout << "y3_new = [" << y3_new << "];\n";
511 }
512 
513 void test08() {
514  cout << "Very simple interpolation case for the "
515  << "new higher order polynomials.\n";
516 
517  Vector og(1, 5, +1); // 1, 2, 3, 4, 5
518  Vector ng(2, 9, 0.25); // 2.0, 2,25, 2.5, 2.75, 3.0 ... 4.0
519 
520  cout << "Original grid:\n" << og << "\n";
521  cout << "New grid:\n" << ng << "\n";
522 
523  // To store the grid positions:
524  ArrayOfGridPosPoly gp(ng.nelem());
525 
526  Index order = 0; // Interpolation order.
527 
528  gridpos_poly(gp, og, ng, order);
529  cout << "Grid positions:\n" << gp;
530 
531  // To store interpolation weights:
532  Matrix itw(gp.nelem(), order + 1);
533  interpweights(itw, gp);
534 
535  cout << "Interpolation weights:\n" << itw << "\n";
536 
537  // Original field:
538  Vector of(og.nelem(), 0);
539  of[2] = 10; // 0, 0, 10, 0, 0
540 
541  cout << "Original field:\n" << of << "\n";
542 
543  // Interpolated field:
544  Vector nf(ng.nelem());
545 
546  interp(nf, itw, of, gp);
547 
548  cout << "New field (order=" << order << "):\n" << nf << "\n";
549 
550  cout << "All orders systematically:\n";
551  for (order = 0; order < 5; ++order) {
552  gridpos_poly(gp, og, ng, order);
553  itw.resize(gp.nelem(), order + 1);
554  interpweights(itw, gp);
555  interp(nf, itw, of, gp);
556 
557  cout << "order " << order << ": ";
558  for (Index i = 0; i < nf.nelem(); ++i) cout << setw(8) << nf[i] << " ";
559  cout << "\n";
560  }
561 }
562 
563 int main() { test08(); }
Matrix
The Matrix class.
Definition: matpackI.h:1193
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
interp_poly
Numeric interp_poly(ConstVectorView x, ConstVectorView y, const Numeric &x_i, const GridPos &gp)
Polynomial interpolation.
Definition: interpolation.cc:2723
gridpos
void gridpos(ArrayOfGridPos &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Numeric &extpolfac)
Set up a grid position Array.
Definition: interpolation.cc:156
test02
void test02(Index n)
Definition: test_interpolation.cc:175
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
test08
void test08()
Definition: test_interpolation.cc:513
joker
const Joker joker
interpolation.h
Header file for interpolation.cc.
interp
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
Definition: interpolation.cc:970
test01
void test01()
Definition: test_interpolation.cc:27
array.h
This file contains the definition of Array.
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
Array
This can be used to make arrays out of anything.
Definition: array.h:108
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
test07
void test07()
Definition: test_interpolation.cc:440
test04
void test04()
Definition: test_interpolation.cc:199
Iterator1D
The iterator class for sub vectors.
Definition: matpackI.h:360
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
math_funcs.h
test06
void test06()
Definition: test_interpolation.cc:290
interpolation_poly.h
Header file for interpolation_poly.cc.
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
test05
void test05()
Definition: test_interpolation.cc:255
Range
The range class.
Definition: matpackI.h:160
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
VectorView::end
Iterator1D end()
Return iterator behind last element.
Definition: matpackI.cc:148
main
int main()
Definition: test_interpolation.cc:563
VectorView::begin
Iterator1D begin()
Return iterator to first element.
Definition: matpackI.cc:144
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Tensor6
The Tensor6 class.
Definition: matpackVI.h:1088
gridpos_poly
void gridpos_poly(ArrayOfGridPosPoly &gp, ConstVectorView old_grid, ConstVectorView new_grid, const Index order, const Numeric &extpolfac)
The maximum difference from 1 that we allow for a sum check.
Definition: interpolation_poly.cc:120
Vector
The Vector class.
Definition: matpackI.h:860
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
test03
void test03(Index n)
Definition: test_interpolation.cc:184
interpweights
void interpweights(VectorView itw, const GridPos &tc)
Red 1D interpolation weights.
Definition: interpolation.cc:741
Tensor7
The Tensor7 class.
Definition: matpackVII.h:2382
matpackVII.h
xml_io.h
This file contains basic functions to handle XML data files.