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