source: src/Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.cpp@ 06653a

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 Frederik Heber. All rights reserved.
 * 
 *
 *   This file is part of MoleCuilder.
 *
 *    MoleCuilder 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 of the License, or
 *    (at your option) any later version.
 *
 *    MoleCuilder 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 MoleCuilder.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * SamplingGridUnitTest.cpp
 *
 *  Created on: Jul 29, 2012
 *      Author: heber
 */

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

using namespace std;

#include <cppunit/CompilerOutputter.h>
#include <cppunit/extensions/TestFactoryRegistry.h>
#include <cppunit/ui/text/TestRunner.h>

// include headers that implement a archive in simple text format
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

#include "SamplingGridUnitTest.hpp"

#include "CodePatterns/Assert.hpp"

#include <boost/assign.hpp>

#include <cmath>
#include <numeric>

#ifdef HAVE_TESTRUNNER
#include "UnitTestMain.hpp"
#endif /*HAVE_TESTRUNNER*/

using namespace boost::assign;

/********************************************** Test classes **************************************/

const double grid_value=1.;

#define NUMBEROFSAMPLES(n) (size_t)(pow(pow(2,n),(int)NDIM))
#define DOMAINVOLUME(l) (size_t)pow(l,(int)NDIM)

// Registers the fixture into the 'registry'
CPPUNIT_TEST_SUITE_REGISTRATION( SamplingGridTest );


void SamplingGridTest::setUp()
{
  // failing asserts should be thrown
  ASSERT_DO(Assert::Throw);

  // create the grid
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 1., 1., 1. };
  for (size_t i=0; i< DOMAINVOLUME(1)*NUMBEROFSAMPLES(2); ++i)
    values += grid_value;
  grid = new SamplingGrid(begin, end, 2, values);
  CPPUNIT_ASSERT_EQUAL( grid_value, *(grid->sampled_grid.begin()) );
}


void SamplingGridTest::tearDown()
{
  delete grid;
}

/** UnitTest on equivalent combination of props and values
 */
void SamplingGridTest::equivalentGrids_Test()
{
  // check illegal grid
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 1., 1., 1. };
  SamplingGridProperties illegal_props(begin, end, 5);
  SamplingGridProperties legal_props(begin, end, 2);
  CPPUNIT_ASSERT( !grid->isEquivalent(illegal_props) );
  CPPUNIT_ASSERT( grid->isEquivalent(legal_props) );
  SamplingGrid::sampledvalues_t illegal_values;
  for (size_t i=0; i< NUMBEROFSAMPLES(1); ++i)
    illegal_values += 1.5;
  SamplingGrid::sampledvalues_t legal_values;
  for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
    legal_values += 1.5;
#ifndef NDEBUG
  // throws because props and size of illegal_values don't match
  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
  CPPUNIT_ASSERT_THROW( SamplingGrid illegal_grid(illegal_props, illegal_values), Assert::AssertionFailure );
  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
  CPPUNIT_ASSERT_THROW( SamplingGrid illegal_grid(legal_props, illegal_values), Assert::AssertionFailure );
  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
  CPPUNIT_ASSERT_THROW( SamplingGrid illegal_grid(illegal_props, legal_values), Assert::AssertionFailure );
#endif
  // check that grid is still the same
  for (SamplingGrid::sampledvalues_t::const_iterator iter = grid->sampled_grid.begin();
      iter != grid->sampled_grid.end(); ++iter)
    CPPUNIT_ASSERT_EQUAL( grid_value, *iter );
}

/** UnitTest on compatible combination of props and values
 */
void SamplingGridTest::compatibleGrids_Test()
{
  // check illegal grid
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 1., 1., 1. };
  const double otherend[NDIM] = { 2.5, 2.5, 2.5 };
  SamplingGridProperties illegal_props(begin, otherend, 5);
  SamplingGridProperties legal_props(begin, end, 3);
  CPPUNIT_ASSERT( !grid->isCompatible(illegal_props) );
  CPPUNIT_ASSERT( grid->isCompatible(legal_props) );
  SamplingGrid::sampledvalues_t illegal_values;
  for (size_t i=0; i< NUMBEROFSAMPLES(1); ++i)
    illegal_values += 1.5;
  SamplingGrid::sampledvalues_t legal_values;
  for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
    legal_values += 1.5;
#ifndef NDEBUG
  // throws because props and size of illegal_values don't match
  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
  CPPUNIT_ASSERT_THROW( SamplingGrid illegal_grid(illegal_props, illegal_values), Assert::AssertionFailure );
  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
  CPPUNIT_ASSERT_THROW( SamplingGrid illegal_grid(legal_props, illegal_values), Assert::AssertionFailure );
  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
  CPPUNIT_ASSERT_THROW( SamplingGrid illegal_grid(illegal_props, legal_values), Assert::AssertionFailure );
#endif
  // check that grid is still the same
  for (SamplingGrid::sampledvalues_t::const_iterator iter = grid->sampled_grid.begin();
      iter != grid->sampled_grid.end(); ++iter)
    CPPUNIT_ASSERT_EQUAL( grid_value, *iter );
}

/** UnitTest for isCongruent()
 */
void SamplingGridTest::isCongruent_Test()
{
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 2., 2., 2. };
  const double otherbegin[NDIM] = { 0.1, 0.1, 0.1 };
  const double otherend[NDIM] = { 1., 1., 1. };
  SamplingGridProperties illegal_begin_props(otherbegin, end, 5);
  SamplingGridProperties illegal_end_props(begin, otherend, 5);
  SamplingGridProperties illegal_level_props(begin, end, 5);
  SamplingGridProperties legal_props(begin, end, 2);

  // differing windows
//  const double begin_window[NDIM] = { 0.5, 0.5, 0.5 };
//  const double end_window[NDIM] = { 1., 1., 1. };
  const double otherbegin_window[NDIM] = { 0.45, 0.45, 0.45 };
  const double otherend_window[NDIM] = { 1.05, 1.05, 1.05 };

  // check that incompatible grid are also incongruent
  SamplingGrid default_grid(legal_props);
  // note that we always construct a temporary SamplingGrid from given props
  CPPUNIT_ASSERT( default_grid.isEquivalent(illegal_begin_props) == default_grid.isCongruent(illegal_begin_props));
  CPPUNIT_ASSERT( default_grid.isEquivalent(illegal_end_props) == default_grid.isCongruent(illegal_end_props));
  CPPUNIT_ASSERT( default_grid.isEquivalent(illegal_level_props) == default_grid.isCongruent(illegal_level_props));
  CPPUNIT_ASSERT( default_grid.isEquivalent(legal_props) == default_grid.isCongruent(legal_props) );

  default_grid.setWindowSize(begin, end);
  // same window
  {
    SamplingGrid illegal_begin_grid(illegal_begin_props);
    SamplingGrid illegal_end_grid(illegal_end_props);
    SamplingGrid illegal_level_grid(illegal_level_props);
    SamplingGrid legal_grid(legal_props);
    illegal_begin_grid.setWindowSize(begin, end);
    illegal_end_grid.setWindowSize(begin, end);
    illegal_level_grid.setWindowSize(begin, end);
    legal_grid.setWindowSize(begin, end);
    CPPUNIT_ASSERT( !illegal_begin_grid.isCongruent(default_grid) );
    CPPUNIT_ASSERT( !illegal_end_grid.isCongruent(default_grid) );
    CPPUNIT_ASSERT( !illegal_level_grid.isCongruent(default_grid) );
    CPPUNIT_ASSERT( legal_grid.isCongruent(default_grid) );
  }

  // different begin
  {
    SamplingGrid illegal_begin_grid(illegal_begin_props);
    SamplingGrid illegal_end_grid(illegal_end_props);
    SamplingGrid illegal_level_grid(illegal_level_props);
    SamplingGrid legal_grid(legal_props);
    illegal_begin_grid.setWindowSize(otherbegin_window, end);
    illegal_end_grid.setWindowSize(otherbegin_window, end);
    illegal_level_grid.setWindowSize(otherbegin_window, end);
    legal_grid.setWindowSize(begin, end);
    CPPUNIT_ASSERT( !illegal_begin_grid.isCongruent(legal_grid) );
    CPPUNIT_ASSERT( !illegal_end_grid.isCongruent(legal_grid) );
    CPPUNIT_ASSERT( !illegal_level_grid.isCongruent(legal_grid) );
    CPPUNIT_ASSERT( legal_grid.isCongruent(default_grid) );
  }

  // different end
  {
    SamplingGrid illegal_begin_grid(illegal_begin_props);
    SamplingGrid illegal_end_grid(illegal_end_props);
    SamplingGrid illegal_level_grid(illegal_level_props);
    SamplingGrid legal_grid(legal_props);
    illegal_begin_grid.setWindowSize(begin, otherend_window);
    illegal_end_grid.setWindowSize(begin, otherend_window);
    illegal_level_grid.setWindowSize(begin, otherend_window);
    legal_grid.setWindowSize(begin, end);
    CPPUNIT_ASSERT( !illegal_begin_grid.isCongruent(legal_grid) );
    CPPUNIT_ASSERT( !illegal_end_grid.isCongruent(legal_grid) );
    CPPUNIT_ASSERT( !illegal_level_grid.isCongruent(legal_grid) );
    CPPUNIT_ASSERT( legal_grid.isCongruent(default_grid) );
  }

  // different begin and end
  {
    SamplingGrid illegal_begin_grid(illegal_begin_props);
    SamplingGrid illegal_end_grid(illegal_end_props);
    SamplingGrid illegal_level_grid(illegal_level_props);
    SamplingGrid legal_grid(legal_props);
    illegal_begin_grid.setWindowSize(otherbegin_window, otherend_window);
    illegal_end_grid.setWindowSize(otherbegin_window, otherend_window);
    illegal_level_grid.setWindowSize(otherbegin_window, otherend_window);
    legal_grid.setWindowSize(begin, end);
    CPPUNIT_ASSERT( !illegal_begin_grid.isCongruent(legal_grid) );
    CPPUNIT_ASSERT( !illegal_end_grid.isCongruent(legal_grid) );
    CPPUNIT_ASSERT( !illegal_level_grid.isCongruent(legal_grid) );
    CPPUNIT_ASSERT( legal_grid.isCongruent(default_grid) );
  }
}

/** UnitTest for integral()
 */
void SamplingGridTest::integral_Test()
{
  double sum = 0.;
  sum = std::accumulate( grid->sampled_grid.begin(), grid->sampled_grid.end(), sum );
  CPPUNIT_ASSERT_EQUAL( sum*grid->getVolume()/grid->getWindowGridPoints(), grid->integral() );
}

/** UnitTest for getVolume()
 */
void SamplingGridTest::getVolume_Test()
{
  CPPUNIT_ASSERT_EQUAL( 1., grid->getVolume() );
}

/** UnitTest for getWindowSize()
 */
void SamplingGridTest::getWindowSize_Test()
{
  // check size of default grid
  CPPUNIT_ASSERT_EQUAL( (size_t)(NUMBEROFSAMPLES(grid->level)), grid->getWindowGridPoints() );

  // create another one and check its size, too
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 1., 1., 1. };
  for (size_t level = 3; level<=6; ++level) {
    values.clear();
    for (size_t i=0; i< NUMBEROFSAMPLES(level); ++i)
      values += grid_value;
    delete grid;
    // use other pointer in case something fails
    SamplingGrid *tmpgrid = new SamplingGrid(begin, end, level, values);
    grid = tmpgrid;
    CPPUNIT_ASSERT_EQUAL( (size_t)NUMBEROFSAMPLES(level), grid->getWindowGridPoints() );
  }
}

/** UnitTest for extendWindow()
 */
void SamplingGridTest::extendWindow_Test()
{
  // we have a grid with size of one, extend to twice the size and check
  const double begin[NDIM] = { 0., 0., 0. };
  const double size = 2.;
  const double end[NDIM] = { size, size, size };
  double offset[NDIM];
  for (offset[0] = 0.; offset[0] <= 1.; offset[0] += .5)
    for (offset[1] = 0.; offset[1] <= 1.; offset[1] += .5)
      for (offset[2] = 0.; offset[2] <= 1.; offset[2] += .5) {
        const double window_begin[NDIM] = { 0.+offset[0], 0.+offset[1], 0.+offset[2]};
        const double window_end[NDIM] = { 1.+offset[0], 1.+offset[1], 1.+offset[2]};
        SamplingGrid newgrid(begin, end, 2);
        newgrid.setWindowSize(window_begin, window_end);
        // resize values by hand to new window size. Otherwise they get zero'd.
        newgrid.sampled_grid = values;
        newgrid.sampled_grid.resize(NUMBEROFSAMPLES(1));
        newgrid.extendWindow(begin, end);

        // check integral
        CPPUNIT_ASSERT_EQUAL(
            grid_value/(NUMBEROFSAMPLES(grid->level)/NUMBEROFSAMPLES(grid->level)),
            grid->integral()
            );

        // check number of points
        CPPUNIT_ASSERT_EQUAL(
            (size_t)NUMBEROFSAMPLES(grid->level),
            grid->getWindowGridPoints()
            );
      }
}

/** UnitTest for extendWindow() with asymmetric values
 */
void SamplingGridTest::extendWindow_asymmetric_Test()
{
  std::cout << "SamplingGridTest::extendWindow_asymmetric_Test()" << std::endl;
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 2., 2., 2. };
  double offset[NDIM];
  for (offset[0] = 0.; offset[0] <= 1.; offset[0] += .5)
    for (offset[1] = 0.; offset[1] <= 1.; offset[1] += .5)
      for (offset[2] = 0.; offset[2] <= 1.; offset[2] += .5) {
        const double window_begin[NDIM] = { 0.+offset[0], 0.+offset[1], 0.+offset[2]};
        const double window_end[NDIM] = { 1.+offset[0], 1.+offset[1], 1.+offset[2]};
        SamplingGrid newgrid(begin, end, 2);
        CPPUNIT_ASSERT_EQUAL( (size_t)0, newgrid.getWindowGridPoints() );
        newgrid.setWindowSize(window_begin, window_end);
        // window size is only half of domain size
        const size_t max_samples = NUMBEROFSAMPLES(newgrid.level)*pow(0.5,(int)NDIM);
        for (size_t i=0; i< max_samples; ++i)
          newgrid.sampled_grid += grid_value*i;
        const size_t sum_weight = (max_samples)*(max_samples-1)/2;
        const double integral = newgrid.integral();
        newgrid.extendWindow(begin, end);

        // check that integral has remained the same
        CPPUNIT_ASSERT_EQUAL( integral, newgrid.integral() );
        CPPUNIT_ASSERT_EQUAL( grid_value*sum_weight/DOMAINVOLUME(2), newgrid.integral() );
      }
}

/** UnitTest for addOntoWindow()
 */
void SamplingGridTest::addOntoWindow_Test()
{
  // first window is from (0,0,0) to (1,1,1)
  CPPUNIT_ASSERT_EQUAL( 1.*grid_value, grid->integral() );

  // create values for half-sized window
  values.clear();
  for (size_t i=0; i< (size_t)pow(.5*pow(2,2),(int)NDIM); ++i)
    values += grid_value;

  // check that too large a window throws
#ifndef NDEBUG
  const double begin[NDIM] = { .5, .5, .5 };
  const double wrongend[NDIM] = { 1.5, 1.5, 1.5 };
  std::cout << "The following assertion is intended and does not indicate a failure of the test." << std::endl;
  CPPUNIT_ASSERT_THROW( grid->addOntoWindow(begin, wrongend, values, +1.), Assert::AssertionFailure );
#endif

  // create another window from (.5,.5,.5) to (1., 1., 1.)
  double offset[NDIM];
  for (offset[0] = 0.; offset[0] <= .5; offset[0] += .5)
    for (offset[1] = 0.; offset[1] <= .5; offset[1] += .5)
      for (offset[2] = 0.; offset[2] <= .5; offset[2] += .5) {
        const double window_begin[NDIM] = { 0.+offset[0], 0.+offset[1], 0.+offset[2]};
        const double window_end[NDIM] = { .5+offset[0], .5+offset[1], .5+offset[2]};

        SamplingGrid newgrid(*grid);
        // now perform working operation
        newgrid.addOntoWindow(window_begin, window_end, values, +1.);

        // check integral to be one and one eighth times the old value
        CPPUNIT_ASSERT_EQUAL( (1.+pow(.5,(int)NDIM))*grid_value, newgrid.integral() );
      }
}

/** UnitTest for addOntoWindow() with asymmetric values
 */
void SamplingGridTest::addOntoWindow_asymmetric_Test()
{
  const size_t size = grid->end[0]-grid->begin[0];
  // check with asymmetric values
  grid->sampled_grid.clear();
  grid->sampled_grid.resize(DOMAINVOLUME(size)*NUMBEROFSAMPLES(grid->level), 0.);

  for (size_t i=0;i<grid->level*(grid->end[0]-grid->begin[0]);++i)
    grid->sampled_grid[(i*2+0)*2+0] += .5*grid_value;
  for (size_t i=0;i<grid->level*(grid->end[1]-grid->begin[1]);++i)
    grid->sampled_grid[(0*2+i)*2+0] += 1.*grid_value;
  for (size_t i=0;i<grid->level*(grid->end[2]-grid->begin[2]);++i)
    grid->sampled_grid[(0*2+0)*2+i] += 1.5*grid_value;

  const double integral = grid->integral();

  // now perform working operation
  grid->addOntoWindow(grid->begin, grid->end, values, +1.);
  // values is equal to integral of 1.
  CPPUNIT_ASSERT_EQUAL( 1.+integral, grid->integral() );
}

/** UnitTest for operator+=()
 */
void SamplingGridTest::operatorPlusEqual_Test()
{
  // create other grid
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 1., 1., 1. };
  SamplingGrid::sampledvalues_t othervalues;
  const double othergrid_value = 1.5;
  for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
    othervalues += othergrid_value;
  SamplingGrid othergrid(begin, end, 2, othervalues);
  CPPUNIT_ASSERT_EQUAL( othergrid_value, *(othergrid.sampled_grid.begin()) );

  // perform operation
  CPPUNIT_ASSERT_NO_THROW( *grid += othergrid );

  // check the contents of the grid
  const double sum = grid_value+othergrid_value;
  for (SamplingGrid::sampledvalues_t::const_iterator iter = grid->sampled_grid.begin();
      iter != grid->sampled_grid.end(); ++iter)
    CPPUNIT_ASSERT_EQUAL( sum, *iter );
}

/** UnitTest for operator-=()
 */
void SamplingGridTest::operatorMinusEqual_Test()
{
  // create other grid
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 1., 1., 1. };
  SamplingGrid::sampledvalues_t othervalues;
  const double othergrid_value = 1.5;
  for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
    othervalues += othergrid_value;
  SamplingGrid othergrid(begin, end, 2, othervalues);
  CPPUNIT_ASSERT_EQUAL( othergrid_value, *(othergrid.sampled_grid.begin()) );

  // perform operation
  CPPUNIT_ASSERT_NO_THROW( *grid -= othergrid );

  // check the contents of the grid
  const double difference = grid_value-othergrid_value;
  for (SamplingGrid::sampledvalues_t::const_iterator iter = grid->sampled_grid.begin();
      iter != grid->sampled_grid.end(); ++iter)
    CPPUNIT_ASSERT_EQUAL( difference, *iter );
}


/** UnitTest for operator==()
 */
void SamplingGridTest::equality_Test()
{
  const double begin[NDIM] = { 0., 0., 0. };
  const double otherbegin[NDIM] = { .5, 0.1, -0.5 };
  const double end[NDIM] = { 1., 1., 1. };
  const double otherend[NDIM] = { 2., 2., 2. };
  const double begin_window[NDIM] = { 0., 0., 0. };
  const double otherbegin_window[NDIM] = { .75, 0.1, 0. };
  const double end_window[NDIM] = { 1., 1., 1. };
  const double otherend_window[NDIM] = { .9, .9, .9 };

  // create other grid
  SamplingGrid samegrid(begin, end, 2, values);
  SamplingGrid differentbegin(otherbegin, end, 2, values);
  SamplingGrid differentend(begin, otherend, 2, values);
  SamplingGrid::sampledvalues_t morevalues;
  {
    const double othergrid_value = 1.5;
    for (size_t i=0; i< NUMBEROFSAMPLES(4); ++i)
      morevalues += othergrid_value;
  }
  SamplingGrid differentlevel(begin, end, 4, morevalues);
  SamplingGrid::sampledvalues_t othervalues;
  {
    const double othergrid_value = 1.5;
    for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
      othervalues += othergrid_value;
  }
  SamplingGrid differentvalues(begin, end, 2, othervalues);

  // check for same level, begin, and end
  CPPUNIT_ASSERT( *grid == *grid );
  CPPUNIT_ASSERT( *grid == samegrid );
  CPPUNIT_ASSERT( samegrid == *grid);
  CPPUNIT_ASSERT( *grid != differentbegin);
  CPPUNIT_ASSERT( differentbegin != *grid);
  CPPUNIT_ASSERT( *grid != differentend );
  CPPUNIT_ASSERT( differentend != *grid );
  CPPUNIT_ASSERT( *grid != differentlevel );
  CPPUNIT_ASSERT( differentlevel != *grid );

  // check for all but differing values
  CPPUNIT_ASSERT( *grid != differentvalues );
  CPPUNIT_ASSERT( differentvalues != *grid );

  // check for different windows
  SamplingGrid differentwindowbegin(begin, end, 2);
  differentwindowbegin.setWindow(otherbegin_window, end_window);
  SamplingGrid differentwindowend(begin, end, 2);
  differentwindowend.setWindow(begin_window, otherend_window);
  CPPUNIT_ASSERT( *grid != differentwindowbegin );
  CPPUNIT_ASSERT( differentwindowbegin != *grid );
  CPPUNIT_ASSERT( *grid != differentwindowend );
  CPPUNIT_ASSERT( differentwindowend != *grid );

  // check against ZeroInstance
  CPPUNIT_ASSERT( *grid != ZeroInstance<SamplingGrid>() );
  CPPUNIT_ASSERT( samegrid != ZeroInstance<SamplingGrid>() );
  CPPUNIT_ASSERT( differentbegin != ZeroInstance<SamplingGrid>() );
  CPPUNIT_ASSERT( differentend != ZeroInstance<SamplingGrid>() );
  CPPUNIT_ASSERT( differentlevel != ZeroInstance<SamplingGrid>() );
  CPPUNIT_ASSERT( differentvalues != ZeroInstance<SamplingGrid>() );
  CPPUNIT_ASSERT( differentwindowbegin != ZeroInstance<SamplingGrid>() );
  CPPUNIT_ASSERT( differentwindowend != ZeroInstance<SamplingGrid>() );
}

/** UnitTest for serialization
 */
void SamplingGridTest::serializeTest()
{
  // serialize
  std::stringstream outputstream;
  boost::archive::text_oarchive oa(outputstream);
  oa << grid;

  // deserialize
  SamplingGrid *samegrid = NULL;
  std::stringstream returnstream(outputstream.str());
  boost::archive::text_iarchive ia(returnstream);
  ia >> samegrid;

  CPPUNIT_ASSERT( samegrid != NULL );
  CPPUNIT_ASSERT( *grid == *samegrid );

  delete samegrid;
}

#ifdef HAVE_INLINE
inline
#endif
static int calculateIndex(const int N[NDIM], const int &_length)
{
  return N[2] + N[1]*_length + N[0]*_length*_length;
}

#ifdef HAVE_INLINE
inline
#endif
static double calculateDistanceSquared(const int N[NDIM], const int &_length)
{
  return
      ::pow(N[0]/(double)_length-.5,2)
      + ::pow(N[1]/(double)_length-.5,2)
      + ::pow(N[2]/(double)_length-.5,2);
}
/** UnitTest for downsample()
 */
void SamplingGridTest::downsampleTest()
{
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 1., 1., 1. };

  // simple test, one level difference, same value everywhere
  {
    SamplingGrid::sampledvalues_t checkvalues;
    SamplingGrid::sampledvalues_t othervalues;
    const double othergrid_value = 1.5;
    for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
      checkvalues += othergrid_value;
    for (size_t i=0; i< NUMBEROFSAMPLES(3); ++i)
      othervalues += othergrid_value;

    SamplingGrid largegrid(begin, end, 3, othervalues);
//    std::cout << " largegrid value " << largegrid.sampled_grid << std::endl;
    SamplingGrid smallgrid(begin, end, 2);
    SamplingGrid::downsample(smallgrid, largegrid, 2);
    SamplingGrid checkgrid(begin, end, 2, checkvalues);
//    std::cout << " smallgrid value " << smallgrid.sampled_grid << std::endl;
//    std::cout << " checkgrid value " << checkgrid.sampled_grid << std::endl;
    CPPUNIT_ASSERT_EQUAL( smallgrid, checkgrid );
  }

  // simple test, two level difference, same value everywhere
  {
    SamplingGrid::sampledvalues_t checkvalues;
    SamplingGrid::sampledvalues_t othervalues;
    const double othergrid_value = 1.5;
    for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
      checkvalues += othergrid_value;
    for (size_t i=0; i< NUMBEROFSAMPLES(4); ++i)
      othervalues += othergrid_value;

    SamplingGrid largegrid(begin, end, 4, othervalues);
//    std::cout << " largegrid value " << largegrid.sampled_grid << std::endl;
    SamplingGrid smallgrid(begin, end, 2);
    SamplingGrid::downsample(smallgrid, largegrid, 2);
    SamplingGrid checkgrid(begin, end, 2, checkvalues);
//    std::cout << " smallgrid value " << smallgrid.sampled_grid << std::endl;
//    std::cout << " checkgrid value " << checkgrid.sampled_grid << std::endl;
    CPPUNIT_ASSERT_EQUAL( smallgrid, checkgrid );
  }

  // same grid, window equals grids, different values
  {
    const int length_c = ::pow(2,2);
    const int length_f = ::pow(2,3);
    SamplingGrid::sampledvalues_t checkvalues((int)::pow(length_c, NDIM), 0.);
    SamplingGrid::sampledvalues_t othervalues((int)::pow(length_f, NDIM), 0.);
    int N[NDIM];
    for (N[0]=0; N[0]< length_f; ++N[0]) {
      for (N[1]=0; N[1]< length_f; ++N[1]) {
        for (N[2]=0; N[2]< length_f; ++N[2]) {
          const int index = calculateIndex(N, length_f);
          const double dist = calculateDistanceSquared(N, length_f);
          othervalues[index] = cos(M_PI*dist/1.5);
        }
      }
    }
    for (N[0]=0; N[0]< length_c; ++N[0]) {
      for (N[1]=0; N[1]< length_c; ++N[1]) {
        for (N[2]=0; N[2]< length_c; ++N[2]) {
          const int index = calculateIndex(N, length_c);
          const double dist = calculateDistanceSquared(N, length_c);
          checkvalues[index] = cos(M_PI*dist/1.5);
        }
      }
    }

    SamplingGrid largegrid(begin, end, 3, othervalues);
//    std::cout << " largegrid value " << largegrid.sampled_grid << std::endl;
    SamplingGrid smallgrid(begin, end, 2);
    SamplingGrid::downsample(smallgrid, largegrid, 2);
    SamplingGrid checkgrid(begin, end, 2, checkvalues);
//    std::cout << " smallgrid value " << smallgrid.sampled_grid << std::endl;
//    std::cout << " checkgrid value " << checkgrid.sampled_grid << std::endl;
    static const double threshold = 1.2e-1;
    for (N[0]=0; N[0]< length_c; ++N[0])
      for (N[1]=0; N[1]< length_c; ++N[1])
        for (N[2]=0; N[2]< length_c; ++N[2]) {
          const double check_threshold =
                    (((N[0] == 0) || (N[0] == length_c-1))
                  && ((N[1] == 0) || (N[1] == length_c-1))
                  && ((N[2] == 0) || (N[2] == length_c-1)))
              ? 2.*threshold : threshold;
          const int index = calculateIndex(N, length_c);
//          std::cout << "Comparing "
//              << fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index]))
//              << " < " << threshold << std::endl;
          CPPUNIT_ASSERT(
              fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index])) < check_threshold);
        }
  }
}