testSplineInterpolator.cc

Go to the documentation of this file.

#include <cmath>
#include <iostream>
#include <vector>
#include <fstream>

#include <utl/SplineInterpolator.h>

#include <tst/Verify.h>
#include <cppunit/extensions/HelperMacros.h>

using namespace std;
using namespace utl;
using namespace tst;
using boost::multi_array;
using namespace utl::Spline;

class SplineInterpolatorTest : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(SplineInterpolatorTest);
  CPPUNIT_TEST(test1d);
  CPPUNIT_TEST(test1dmulti);
  CPPUNIT_TEST(test2d);
  CPPUNIT_TEST(test2dmulti);
  CPPUNIT_TEST(test3d);
  CPPUNIT_TEST(test3dmulti);
  CPPUNIT_TEST_SUITE_END();

public:
  void setUp() {
    fNX = 17;
    fNY = 13;
    fNZ = 7;
    fOrder = 3;
    fOver = 1;

    fX0 = -1.2;
    fX1 =  5.8;

    fY0 = -7.7;
    fY1 =  1.3;

    fZ0 = -3.9;
    fZ1 =  4.2;
  }

  void tearDown() { }

  void
  test1d()
  {
    vector<double> vs(fNX);
    for (size_t i=0;i<fNX;++i)
    {
      const double d = double(i)/(fNX-1);
      const double x = fX1*d + fX0*(1-d);
      vs[i] = ExpectedValue(x);
    }

    Interpolator1D ip(fX0,fX1,vs,BoundaryCondition(eNatural));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(), fX0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(), fX1, 1e-7));

    for (size_t i = 0; i < fNX*fOver; ++i) {
      const double d = double(i)/(fNX*fOver-1);
      const double x = fX1*d + fX0*(1-d);
      const double vt = ExpectedValue(x);
      CPPUNIT_ASSERT(Verify<CloseTo>(ip(x), vt, 1e-7));
    }
  }

  void
  test1dmulti()
  {
    multi_array<double,2> vs(boost::extents[2][fNX]);
    for (size_t i=0;i<fNX;++i)
    {
      const double d = double(i)/(fNX-1);
      const double x = fX1*d + fX0*(1-d);
      const double v = ExpectedValue(x);
      vs[0][i] = v;
      vs[1][i] = 2*v;
    }

    VectorInterpolator1D ip(fX0,fX1,vs);
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(), fX0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(), fX1, 1e-7));

    vector<double> out(2);
    for (size_t i = 0; i < fNX*fOver; ++i) {
      const double d = double(i)/(fNX*fOver-1);
      const double x = fX1*d + fX0*(1-d);
      ip(out,x);
      const double vt = ExpectedValue(x);
      CPPUNIT_ASSERT(Verify<CloseTo>(out[0], vt  , 1e-7));
      CPPUNIT_ASSERT(Verify<CloseTo>(out[1], 2*vt, 1e-7));
    }
  }

  void
  test2d()
  {
    multi_array<double,2> vs(boost::extents[fNX][fNY]);
    for (size_t i = 0; i < fNX; ++i)
      for (size_t j = 0; j < fNY; ++j)
    {
      const double di = double(i)/(fNX-1);
      const double dj = double(j)/(fNY-1);
      const double x = fX1*di + fX0*(1-di);
      const double y = fY1*dj + fY0*(1-dj);
      vs[i][j] = ExpectedValue(x) + ExpectedValue(y);
    }

    Interpolator2D ip(fX0,fX1,fY0,fY1,vs);
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(0), fX0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(0), fX1, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(1), fY0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(1), fY1, 1e-7));

    for (size_t i = 0; i < fNX*fOver; ++i)
      for (size_t j = 0; j < fNY*fOver; ++j)
    {
      const double di = double(i)/(fNX*fOver-1);
      const double dj = double(j)/(fNY*fOver-1);
      const double x = fX1*di + fX0*(1-di);
      const double y = fY1*dj + fY0*(1-dj);
      const double vt = ExpectedValue(x) + ExpectedValue(y);
      CPPUNIT_ASSERT(Verify<CloseTo>(ip(x,y), vt, 1e-7));
    }
  }

  void
  test2dmulti()
  {
    multi_array<double,3> vs(boost::extents[2][fNX][fNY]);
    for (size_t i = 0; i < fNX; ++i)
      for (size_t j = 0; j < fNY; ++j)
    {
      const double di = double(i)/(fNX-1);
      const double dj = double(j)/(fNY-1);
      const double x = fX1*di + fX0*(1-di);
      const double y = fY1*dj + fY0*(1-dj);
      const double vt = ExpectedValue(x) + ExpectedValue(y);
      vs[0][i][j] = vt;
      vs[1][i][j] = 2*vt;
    }

    VectorInterpolator2D ip(fX0,fX1,fY0,fY1,vs);
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(0), fX0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(0), fX1, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(1), fY0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(1), fY1, 1e-7));

    vector<double> out(2);
    for (size_t i = 0; i < fNX*fOver; ++i)
      for (size_t j = 0; j < fNY*fOver; ++j)
    {
      const double di = double(i)/(fNX*fOver-1);
      const double dj = double(j)/(fNY*fOver-1);
      const double x = fX1*di + fX0*(1-di);
      const double y = fY1*dj + fY0*(1-dj);
      ip(out,x,y);
      const double vt = ExpectedValue(x) + ExpectedValue(y);
      CPPUNIT_ASSERT(Verify<CloseTo>(out[0], vt  , 1e-7));
      CPPUNIT_ASSERT(Verify<CloseTo>(out[1], 2*vt, 1e-7));
    }
  }

  void
  test3d()
  {
    multi_array<double,3> vs(boost::extents[fNX][fNY][fNZ]);
    for (size_t i = 0; i < fNX; ++i)
      for (size_t j = 0; j < fNY; ++j)
        for (size_t k = 0; k < fNZ; ++k)
    {
      const double di = double(i)/(fNX-1);
      const double dj = double(j)/(fNY-1);
      const double dk = double(k)/(fNZ-1);
      const double x = fX1*di + fX0*(1-di);
      const double y = fY1*dj + fY0*(1-dj);
      const double z = fZ1*dk + fZ0*(1-dk);
      vs[i][j][k] = ExpectedValue(x)
                  + ExpectedValue(y)
                  + ExpectedValue(z);
    }

    Interpolator3D ip(fX0,fX1,fY0,fY1,fZ0,fZ1,vs);
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(0), fX0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(0), fX1, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(1), fY0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(1), fY1, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(2), fZ0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(2), fZ1, 1e-7));

    for (size_t i = 0; i < fNX*fOver; ++i)
      for (size_t j = 0; j < fNY*fOver; ++j)
        for (size_t k = 0; k < fNZ*fOver; ++k)
    {
      const double di = double(i)/(fNX*fOver-1);
      const double dj = double(j)/(fNY*fOver-1);
      const double dk = double(k)/(fNZ*fOver-1);
      const double x = fX1*di + fX0*(1-di);
      const double y = fY1*dj + fY0*(1-dj);
      const double z = fZ1*dk + fZ0*(1-dk);
      const double vt = ExpectedValue(x)
                      + ExpectedValue(y)
                      + ExpectedValue(z);
      CPPUNIT_ASSERT(Verify<CloseTo>(ip(x,y,z), vt, 1e-7));
    }
  }

  void
  test3dmulti()
  {
    multi_array<double,4> vs(boost::extents[2][fNX][fNY][fNZ]);
    for (size_t i = 0; i < fNX; ++i)
      for (size_t j = 0; j < fNY; ++j)
        for (size_t k = 0; k < fNZ; ++k)
    {
      const double di = double(i)/(fNX-1);
      const double dj = double(j)/(fNY-1);
      const double dk = double(k)/(fNZ-1);
      const double x = fX1*di + fX0*(1-di);
      const double y = fY1*dj + fY0*(1-dj);
      const double z = fZ1*dk + fZ0*(1-dk);
      const double vt = ExpectedValue(x)
                      + ExpectedValue(y)
                      + ExpectedValue(z);
      vs[0][i][j][k] = vt;
      vs[1][i][j][k] = 2*vt;
    }

    VectorInterpolator3D ip(fX0,fX1,fY0,fY1,fZ0,fZ1,vs);
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(0), fX0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(0), fX1, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(1), fY0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(1), fY1, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStart(2), fZ0, 1e-7));
    CPPUNIT_ASSERT(Verify<CloseTo>(ip.GetStop(2), fZ1, 1e-7));

    vector<double> out(2);
    for (size_t i = 0; i < fNX*fOver; ++i)
      for (size_t j = 0; j < fNY*fOver; ++j)
        for (size_t k = 0; k < fNZ*fOver; ++k)
    {
      const double di = double(i)/(fNX*fOver-1);
      const double dj = double(j)/(fNY*fOver-1);
      const double dk = double(k)/(fNZ*fOver-1);
      const double x = fX1*di + fX0*(1-di);
      const double y = fY1*dj + fY0*(1-dj);
      const double z = fZ1*dk + fZ0*(1-dk);
      ip(out,x,y,z);
      const double vt = ExpectedValue(x)
                      + ExpectedValue(y)
                      + ExpectedValue(z);
      CPPUNIT_ASSERT(Verify<CloseTo>(out[0], vt, 1e-7));
      CPPUNIT_ASSERT(Verify<CloseTo>(out[1], 2*vt, 1e-7));
    }
  }

private:

  size_t fNX, fNY, fNZ, fOrder, fOver;
  double fX0, fX1, fY0, fY1, fZ0, fZ1;

  double
  ExpectedValue(const double x)
    const
  {
    const double factors[] = { 2.9, 6.1, 9.5 };
    double result = 0;
    for (size_t o=0; o<fOrder; ++o)
      result += factors[o]*pow(x, int(o));
    return result;
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(SplineInterpolatorTest);


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: