test_interpolation.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2012 Stefan Buehler <sbuehler@ltu.se>

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */

#include <iostream>
#include <cmath>
#include "array.h"
#include "matpackVII.h"
#include "interpolation.h"
#include "interpolation_poly.h"
#include "make_vector.h"
#include "xml_io.h"
#include "math_funcs.h"

void test01()
{
  cout << "Simple interpolation cases\n"
       << "--------------------------\n";
  //  Vector og(5,5,-1);                // 5,4,3,2,1
  Vector og(1,5,+1);            // 1, 2, 3, 4, 5
  Vector ng(2,5,0.25);          // 2.0, 2,25, 2.5, 2.75, 3.0

  cout << "og:\n" << og << "\n";
  cout << "ng:\n" << ng << "\n";

  // To store the grid positions:
  ArrayOfGridPos gp(ng.nelem());

  gridpos(gp,og,ng);
  cout << "gp:\n" << gp << "\n";

  cout << "1D:\n"
       << "---\n";
  {
    // To store interpolation weights:
    Matrix itw(gp.nelem(),2);
    interpweights(itw,gp);
    
    cout << "itw:\n" << itw << "\n";

    // Original field:
    Vector of(og.nelem(),0);
    of[2] = 10;                 // 0, 0, 10, 0, 0

    cout << "of:\n" << of << "\n";

    // Interpolated field:
    Vector nf(ng.nelem());

    interp(nf, itw, of, gp);

    cout << "nf:\n" << nf << "\n";
  }

  cout << "Blue 2D:\n"
       << "--------\n";
  {
    // To store interpolation weights:
    Matrix itw(gp.nelem(),4);
    interpweights(itw,gp,gp);
    
    cout << "itw:\n" << itw << "\n";

    // Original field:
    Matrix of(og.nelem(),og.nelem(),0);
    of(2,2) = 10;                       // 0 Matrix with 10 in the middle

    cout << "of:\n" << of << "\n";

    // Interpolated field:
    Vector nf(ng.nelem());

    interp(nf, itw, of, gp, gp);

    cout << "nf:\n" << nf << "\n";
  }

  cout << "Blue 6D:\n"
       << "--------\n";
  {
    // To store interpolation weights:
    Matrix itw(gp.nelem(),64);
    interpweights(itw,gp,gp,gp,gp,gp,gp);
    
    //    cout << "itw:\n" << itw << "\n";

    // Original field:
    Index n = og.nelem();
    Tensor6 of(n,n,n,n,n,n,0);
    of(2,2,2,2,2,2) = 10;                       // 0 Tensor with 10 in the middle

    //    cout << "of:\n" << of << "\n";

    // Interpolated field:
    Vector nf(ng.nelem());

    interp(nf, itw, of, gp, gp, gp, gp, gp, gp);

    cout << "nf:\n" << nf << "\n";
  }

  cout << "Green 2D:\n"
       << "---------\n";
  {
    // To store interpolation weights:
    Tensor3 itw(gp.nelem(),gp.nelem(),4);
    interpweights(itw,gp,gp);
    
    for ( Index i=0; i<itw.ncols(); ++i )
      cout << "itw " << i << ":\n" << itw(Range(joker),Range(joker),i) << "\n";

    // Original field:
    Matrix of(og.nelem(),og.nelem(),0);
    of(2,2) = 10;                       // 0 Matrix with 10 in the middle

    cout << "of:\n" << of << "\n";

    // Interpolated field:
    Matrix nf(ng.nelem(),ng.nelem());

    interp(nf, itw, of, gp, gp);

    cout << "nf:\n" << nf << "\n";
  }

  cout << "Green 6D:\n"
       << "---------\n";
  {
    // To store interpolation weights:
    Tensor7 itw(gp.nelem(),
                gp.nelem(),
                gp.nelem(),
                gp.nelem(),
                gp.nelem(),
                gp.nelem(),
                64);
    interpweights(itw,gp,gp,gp,gp,gp,gp);
    
    // Original field:
    Tensor6 of(og.nelem(),og.nelem(),og.nelem(),og.nelem(),og.nelem(),og.nelem(),0);
    of(2,2,2,2,2,2) = 10;                       // 0 Tensor with 10 in the middle

    cout << "Middle slice of of:\n" << of(2,2,2,2,Range(joker),Range(joker)) << "\n";

    // Interpolated field:
    Tensor6 nf(ng.nelem(),ng.nelem(),ng.nelem(),ng.nelem(),ng.nelem(),ng.nelem());

    interp(nf, itw, of, gp, gp, gp, gp, gp, gp);

    cout << "Last slice of nf:\n" << nf(4,4,4,4,Range(joker),Range(joker)) << "\n";
  }
}

void test02(Index n)
{
  cout << "Test whether for loop or iterator are quicker\n"
       << "a) for loop\n"
       << "---------------------------------------------\n";

  Vector a(n);
  for (Index i=0; i<a.nelem(); ++i)
    a[i] = (Numeric)i;
}

void test03(Index n)
{
  cout << "Test whether for loop or iterator are quicker\n"
       << "b) iterator\n"
       << "---------------------------------------------\n";

  Vector a(n);
  Iterator1D ai=a.begin();
  const Iterator1D ae=a.end();
  Index i=0;
  for ( ; ai!=ae; ++ai, ++i )
    *ai = (Numeric)i;
}

// Result: Both are almost equally fast, with a slight advantage of
// the for loop if compiler optimization is enabled.

void test04()
{
  cout << "Green type interpolation of all pages of a Tensor3\n";

  // The original Tensor is called a, the new one n. 

  // 10 pages, 20 rows, 30 columns, all grids are: 1,2,3
  Vector  a_pgrid(1,3,1), a_rgrid(1,3,1), a_cgrid(1,3,1); 
  Tensor3 a( a_pgrid.nelem(), a_rgrid.nelem(), a_cgrid.nelem() ); 

  a = 0;
  // Put some simple numbers in the middle of each page:
  a(0,1,1) = 10;
  a(1,1,1) = 20;
  a(2,1,1) = 30;

  // New row and column grids:
  // 1, 1.5, 2, 2.5, 3
  Vector  n_rgrid(1,5,.5), n_cgrid(1,5,.5); 
  Tensor3 n( a_pgrid.nelem(), n_rgrid.nelem(), n_cgrid.nelem() ); 

  // So, n has the same number of pages as a, but more rows and columns.

  // Get the grid position arrays:
  ArrayOfGridPos n_rgp(n_rgrid.nelem()); // For row grid positions.
  ArrayOfGridPos n_cgp(n_cgrid.nelem()); // For column grid positions.

  gridpos( n_rgp, a_rgrid, n_rgrid );
  gridpos( n_cgp, a_cgrid, n_cgrid );

  // Get the interpolation weights:
  Tensor3 itw( n_rgrid.nelem(), n_cgrid.nelem(), 4 );
  interpweights( itw, n_rgp, n_cgp );

  // Do a "green" interpolation for all pages of the Tensor a:

  for ( Index i=0; i<a.npages(); ++i )
    {
      // Select the current page of both a and n:
      ConstMatrixView ap = a( i, Range(joker), Range(joker) );
      MatrixView      np = n( i, Range(joker), Range(joker) );

      // Do the interpolation:
      interp( np, itw, ap, n_rgp, n_cgp );

      // Note that this is efficient, because interpolation weights and
      // grid positions are re-used.
    }

  cout << "Original field:\n";
  for ( Index i=0; i<a.npages(); ++i )
      cout << "page " << i << ":\n" << a(i,Range(joker),Range(joker)) << "\n";

  cout << "Interpolated field:\n";
  for ( Index i=0; i<n.npages(); ++i )
      cout << "page " << i << ":\n" << n(i,Range(joker),Range(joker)) << "\n";
}

void test05()
{
  cout << "Very simple interpolation case\n";

  Vector og(1,5,+1);            // 1, 2, 3, 4, 5
  Vector ng(2,5,0.25);          // 2.0, 2,25, 2.5, 2.75, 3.0

  cout << "Original grid:\n" << og << "\n";
  cout << "New grid:\n" << ng << "\n";

  // To store the grid positions:
  ArrayOfGridPos gp(ng.nelem());

  gridpos(gp,og,ng);
  cout << "Grid positions:\n" << gp;

  // To store interpolation weights:
  Matrix itw(gp.nelem(),2);
  interpweights(itw,gp);
    
  cout << "Interpolation weights:\n" << itw << "\n";

  // Original field:
  Vector of(og.nelem(),0);
  of[2] = 10;                   // 0, 0, 10, 0, 0

  cout << "Original field:\n" << of << "\n";

  // Interpolated field:
  Vector nf(ng.nelem());

  interp(nf, itw, of, gp);

  cout << "New field:\n" << nf << "\n";
}

void test06()
{
  cout << "Simple extrapolation cases\n"
       << "--------------------------\n";
  //  Vector og(5,5,-1);                // 5,4,3,2,1
  Vector og(1,5,+1);                // 1, 2, 3, 4, 5
  MakeVector ng(0.9,1.5,3,4.5,5.1); // 0.9, 1.5, 3, 4.5, 5.1

  cout << "og:\n" << og << "\n";
  cout << "ng:\n" << ng << "\n";

  // To store the grid positions:
  ArrayOfGridPos gp(ng.nelem());

  gridpos(gp,og,ng);
  cout << "gp:\n" << gp << "\n";

  cout << "1D:\n"
       << "---\n";
  {
    // To store interpolation weights:
    Matrix itw(gp.nelem(),2);
    interpweights(itw,gp);
    
    cout << "itw:\n" << itw << "\n";

    // Original field:
    Vector of(og.nelem(),0);
    for ( Index i=0; i<og.nelem(); ++i )
      of[i] = (Numeric)(10*(i+1));                 // 10, 20, 30, 40, 50

    cout << "of:\n" << of << "\n";

    // Interpolated field:
    Vector nf(ng.nelem());

    interp(nf, itw, of, gp);

    cout << "nf:\n" << nf << "\n";
  }

  cout << "Blue 2D:\n"
       << "--------\n";
  {
    // To store interpolation weights:
    Matrix itw(gp.nelem(),4);
    interpweights(itw,gp,gp);
    
    cout << "itw:\n" << itw << "\n";

    // Original field:
    Matrix of(og.nelem(),og.nelem(),0);
    of(2,2) = 10;                       // 0 Matrix with 10 in the middle

    cout << "of:\n" << of << "\n";

    // Interpolated field:
    Vector nf(ng.nelem());

    interp(nf, itw, of, gp, gp);

    cout << "nf:\n" << nf << "\n";
  }

  cout << "Blue 6D:\n"
       << "--------\n";
  {
    // To store interpolation weights:
    Matrix itw(gp.nelem(),64);
    interpweights(itw,gp,gp,gp,gp,gp,gp);
    
    //    cout << "itw:\n" << itw << "\n";

    // Original field:
    Index n = og.nelem();
    Tensor6 of(n,n,n,n,n,n,0);
    of(2,2,2,2,2,2) = 10;                       // 0 Tensor with 10 in the middle

    //    cout << "of:\n" << of << "\n";

    // Interpolated field:
    Vector nf(ng.nelem());

    interp(nf, itw, of, gp, gp, gp, gp, gp, gp);

    cout << "nf:\n" << nf << "\n";
  }

  cout << "Green 2D:\n"
       << "---------\n";
  {
    // To store interpolation weights:
    Tensor3 itw(gp.nelem(),gp.nelem(),4);
    interpweights(itw,gp,gp);
    
    for ( Index i=0; i<itw.ncols(); ++i )
      cout << "itw " << i << ":\n" << itw(Range(joker),Range(joker),i) << "\n";

    // Original field:
    Matrix of(og.nelem(),og.nelem(),0);
    of(2,2) = 10;                       // 0 Matrix with 10 in the middle

    cout << "of:\n" << of << "\n";

    // Interpolated field:
    Matrix nf(ng.nelem(),ng.nelem());

    interp(nf, itw, of, gp, gp);

    cout << "nf:\n" << nf << "\n";
  }

  cout << "Green 6D:\n"
       << "---------\n";
  {
    // To store interpolation weights:
    Tensor7 itw(gp.nelem(),
                gp.nelem(),
                gp.nelem(),
                gp.nelem(),
                gp.nelem(),
                gp.nelem(),
                64);
    interpweights(itw,gp,gp,gp,gp,gp,gp);
    
    // Original field:
    Tensor6 of(og.nelem(),og.nelem(),og.nelem(),og.nelem(),og.nelem(),og.nelem(),0);
    of(2,2,2,2,2,2) = 10;                       // 0 Tensor with 10 in the middle

    cout << "Middle slice of of:\n" << of(2,2,2,2,Range(joker),Range(joker)) << "\n";

    // Interpolated field:
    Tensor6 nf(ng.nelem(),ng.nelem(),ng.nelem(),ng.nelem(),ng.nelem(),ng.nelem());

    interp(nf, itw, of, gp, gp, gp, gp, gp, gp);

    cout << "Last slice of nf:\n" << nf(4,4,4,4,Range(joker),Range(joker)) << "\n";
  }
}

// Test polynomial interpolation (included by CE)
void test07()
{
  // FileType ftype = FILE_TYPE_ASCII;
  Vector new_x(0, 21, +0.25);
  Vector x(0, 10, +1);

  ArrayOfGridPos gp(new_x.nelem());
  gridpos(gp, x, new_x);
  
  Vector y1(x.nelem());
  Vector y2(x.nelem());
  Vector y3(x.nelem());
  
  for(Index i = 0; i < x.nelem(); i++)
    { 
      // linear function
      y1[i] = 3*x[i];
      // cubic function
      y2[i] = pow(x[i],3) + 2;
      // trigonometric function
      y3[i] = sin(x[i]);
    }
  
  // Linear interpolation:
   Matrix itw(gp.nelem(),2);
   interpweights(itw, gp);

   Vector y1_lin(new_x.nelem());
   Vector y2_lin(new_x.nelem());
   Vector y3_lin(new_x.nelem());

   interp(y1_lin, itw, y1, gp);
   interp(y2_lin, itw, y2, gp);
   interp(y3_lin, itw, y3, gp);

   cout << "y1_lin = ["<< y1_lin << "];\n";
   cout << "y2_lin = ["<< y2_lin << "];\n";
   cout << "y3_lin = ["<< y3_lin << "];\n";

   // Cubic interpolation:
   Vector y1_cub(new_x.nelem());
   Vector y2_cub(new_x.nelem());
   Vector y3_cub(new_x.nelem());

   for(Index i = 0; i < new_x.nelem(); i++)
     {
       y1_cub[i] = interp_poly(x, y1, new_x[i], gp[i]);
       y2_cub[i] = interp_poly(x, y2, new_x[i], gp[i]);
       y3_cub[i] = interp_poly(x, y3, new_x[i], gp[i]);
     }
   
   cout << "y1_cub = ["<< y1_cub << "];\n";
   cout << "y2_cub = ["<< y2_cub << "];\n";
   cout << "y3_cub = ["<< y3_cub << "];\n";

   // Stefan's new polynomial interpolation routines:
   Index order = 2;

   Vector y1_new(new_x.nelem());
   Vector y2_new(new_x.nelem());
   Vector y3_new(new_x.nelem());

   ArrayOfGridPosPoly gpp(new_x.nelem());
   gridpos_poly(gpp, x, new_x, order);   
   Matrix itwp(new_x.nelem(), order+1);
   interpweights(itwp, gpp);

   interp(y1_new, itwp, y1, gpp);
   interp(y2_new, itwp, y2, gpp);
   interp(y3_new, itwp, y3, gpp);

   cout << "y1_new = ["<< y1_new << "];\n";
   cout << "y2_new = ["<< y2_new << "];\n";
   cout << "y3_new = ["<< y3_new << "];\n";
}
   
void test08()
{
  cout << "Very simple interpolation case for the "
       << "new higher order polynomials.\n";

  Vector og(1,5,+1);            // 1, 2, 3, 4, 5
  Vector ng(2,9,0.25);          // 2.0, 2,25, 2.5, 2.75, 3.0 ... 4.0

  cout << "Original grid:\n" << og << "\n";
  cout << "New grid:\n" << ng << "\n";

  // To store the grid positions:
  ArrayOfGridPosPoly gp(ng.nelem());

  Index order=0;                // Interpolation order.

  gridpos_poly(gp,og,ng,order);
  cout << "Grid positions:\n" << gp;

  // To store interpolation weights:
  Matrix itw(gp.nelem(),order+1);
  interpweights(itw,gp);
    
  cout << "Interpolation weights:\n" << itw << "\n";

  // Original field:
  Vector of(og.nelem(),0);
  of[2] = 10;                   // 0, 0, 10, 0, 0

  cout << "Original field:\n" << of << "\n";

  // Interpolated field:
  Vector nf(ng.nelem());

  interp(nf, itw, of, gp);

  cout << "New field (order=" << order << "):\n" << nf << "\n";  

  cout << "All orders systematically:\n";
  for (order=0; order<5; ++order)
    {
      gridpos_poly(gp,og,ng,order);
      itw.resize(gp.nelem(),order+1);
      interpweights(itw,gp);
      interp(nf, itw, of, gp);

      cout << "order " << order << ": ";
      for (Index i=0; i<nf.nelem(); ++i)
        cout << setw(8) << nf[i] << " ";
      cout << "\n";
    }
}

int main()
{
  test08();
}