/*
 * 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 "CodePatterns/Log.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);
}

#ifdef HAVE_INLINE
inline
#endif
static double getBoundaryCaseFactor(const int N[NDIM], const int &_length)
{
  static double third_two=::pow(4., 1./3.);
  double returnthreshold = 1.;
  for (size_t i=0;i<NDIM;++i)
    if ((N[i] == 0) || (N[i] == _length-1))
      returnthreshold *= third_two;
  return returnthreshold;
}

/** UnitTest for padWithZerosForEvenNumberedSamples()
 *
 */
void SamplingGridTest::padWithZerosForEvenNumberedSamplesTest()
{
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 2., 2., 2. };

  // we check for window of odd numbers
  {
  //    LOG(2, "Checking with difference " << grid_difference);
    const int gridpoints = 1*2*2;
    const int extended_gridpoints = 2*4*4;
    SamplingGrid::sampledvalues_t values(gridpoints, 1.); // length in one axis is odd
    SamplingGrid::sampledvalues_t checkvalues(extended_gridpoints, 0.);
    {
      int N[NDIM];
      int index = 0;
      for (N[0]=2;N[0]<3;++N[0]) {
        for (N[1]=0;N[1]<1;++N[1]) {
          for (N[2]=0;N[2]<4;++N[2]) {
            checkvalues[index++] = 0.;
          }
        }
        for (N[1]=1;N[1]<3;++N[1]) {
          checkvalues[index++] = 0.;
          for (N[2]=1;N[2]<3;++N[2]) {
            checkvalues[index++] = 1.;
          }
          checkvalues[index++] = 0.;
        }
        for (N[1]=3;N[1]<4;++N[1]) {
          for (N[2]=0;N[2]<4;++N[2]) {
            checkvalues[index++] = 0.;
          }
        }
      }
      for (N[0]=3;N[0]<4;++N[0]) {
        for (N[1]=0;N[1]<4;++N[1]) {
          for (N[2]=0;N[2]<4;++N[2]) {
            checkvalues[index++] = 0.;
          }
        }
      }
    }

    /** Here, we must initialize the larger grid with zero window first and
     * then add the values into the smaller window.
     */
    SamplingGrid grid(begin, end, 2);
    const double window_begin[NDIM] = {
        grid.getNearestHigherGridPoint(1.0, 0), // index 2
        grid.getNearestHigherGridPoint(0.5, 1), // index 1 -> 0
        grid.getNearestHigherGridPoint(0.5, 2) }; // index 1 -> 0
    const double window_end[NDIM] = {
        grid.getNearestLowerGridPoint(1.5, 0), // index 3 -> 4
        grid.getNearestLowerGridPoint(1.5, 1), // index 3 -> 4
        grid.getNearestLowerGridPoint(1.5, 2) }; // index 3 -> 4
    grid.setWindow(window_begin, window_end);
    grid.addOntoWindow(window_begin, window_end, values, +1.);
    grid.padWithZerosForEvenNumberedSamples();

    SamplingGrid checkgrid(begin, end, 2);
    const double check_window_begin[NDIM] = {
        checkgrid.getNearestHigherGridPoint(1., 0),
        checkgrid.getNearestHigherGridPoint(0., 1),
        checkgrid.getNearestHigherGridPoint(0., 2) };
    const double check_window_end[NDIM] = {
        checkgrid.getNearestHigherGridPoint(2., 0),
        checkgrid.getNearestHigherGridPoint(2., 1),
        checkgrid.getNearestHigherGridPoint(2., 2) };
    checkgrid.setWindow(check_window_begin, check_window_end);
    checkgrid.addOntoWindow(check_window_begin, check_window_end, checkvalues, +1.);

    std::cout << " grid value " << grid.sampled_grid << std::endl;
    std::cout << " checkgrid value " << checkgrid.sampled_grid << std::endl;
    CPPUNIT_ASSERT_EQUAL( checkgrid.sampled_grid.size(), grid.sampled_grid.size());
    CPPUNIT_ASSERT_EQUAL( checkgrid, grid );
  }
}

/** UnitTest for downsample()
 */
void SamplingGridTest::downsample_gridTest()
{
  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, ever larger largegrids
  const int base_level = 2;
  for (int grid_difference = 1; grid_difference <= 3; ++grid_difference) {
//    LOG(2, "Checking with difference " << grid_difference);
    const int length_c = ::pow(2, base_level);
    const int length_f = ::pow(2, base_level+grid_difference);
    SamplingGrid::sampledvalues_t checkvalues((int)::pow(length_c, (int)NDIM), 0.);
    SamplingGrid::sampledvalues_t othervalues((int)::pow(length_f, (int)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, base_level+grid_difference, othervalues);
//    std::cout << " largegrid value " << largegrid.sampled_grid << std::endl;
    SamplingGrid smallgrid(begin, end, base_level);
    SamplingGrid::downsample(smallgrid, largegrid, base_level);
    SamplingGrid checkgrid(begin, end, base_level, checkvalues);
//    std::cout << " smallgrid value " << smallgrid.sampled_grid << std::endl;
//    std::cout << " checkgrid value " << checkgrid.sampled_grid << std::endl;
    const double threshold = 3.2e-2;
    // we skip the first value as checkvalues[0] is always zero, hard to get from downsampling
    for (N[0]=1; 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 =
              threshold*(1.+grid_difference/4.)*getBoundaryCaseFactor(N, length_c);
          const int index = calculateIndex(N, length_c);
//          std::cout << "Comparing "
//              << fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index]))
//              << " < " << check_threshold << std::endl;
          CPPUNIT_ASSERT(
              fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index])) < check_threshold);
//          if (fabs(checkgrid.sampled_grid[index]) > 1e-10) {
//            std::cout << "Comparing "
//                << fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index])/checkgrid.sampled_grid[index])
//                << " < " << check_threshold << std::endl;
//            CPPUNIT_ASSERT(
//                fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index])/checkgrid.sampled_grid[index]) < check_threshold);
//          }
        }
  }
}

/** UnitTest for downsample()
 */
void SamplingGridTest::downsample_smallerwindowTest()
{
  const double begin[NDIM] = { 0., 0., 0. };
  const double end[NDIM] = { 2., 2., 2. };

  // same grid, window is half of grid, different level by one
  const int base_level = 3;
  for (int grid_difference = 1; grid_difference <= 3; ++grid_difference) {
//    LOG(2, "Checking with difference " << grid_difference);
    const int window_length_c = ::pow(2, base_level-1);
    const int window_length_f = ::pow(2, base_level+grid_difference-1);
    SamplingGrid::sampledvalues_t checkvalues((int)::pow(window_length_c, (int)NDIM), 0.);
    SamplingGrid::sampledvalues_t othervalues((int)::pow(window_length_f, (int)NDIM), 0.);
    int N[NDIM];
    for (N[0]=0; N[0]< window_length_f; ++N[0]) {
      for (N[1]=0; N[1]< window_length_f; ++N[1]) {
        for (N[2]=0; N[2]< window_length_f; ++N[2]) {
          const int index = calculateIndex(N, window_length_f);
          const double dist = calculateDistanceSquared(N, window_length_f);
          othervalues[index] = cos(M_PI*dist/1.5);
        }
      }
    }
    for (N[0]=0; N[0]< window_length_c; ++N[0]) {
      for (N[1]=0; N[1]< window_length_c; ++N[1]) {
        for (N[2]=0; N[2]< window_length_c; ++N[2]) {
          const int index = calculateIndex(N, window_length_c);
          const double dist = calculateDistanceSquared(N, window_length_c);
          checkvalues[index] = cos(M_PI*dist/1.5);
        }
      }
    }

    /** Here, we must initialize the larger grid with zero window first and
     * then add the values into the smaller window.
     */
    SamplingGrid largegrid(begin, end, base_level+grid_difference);
    const double window_begin[NDIM] = {
        largegrid.getNearestHigherGridPoint(0.5, 0),
        largegrid.getNearestHigherGridPoint(0.5, 1),
        largegrid.getNearestHigherGridPoint(0.5, 2) };
    const double window_end[NDIM] = {
        largegrid.getNearestLowerGridPoint(1.5, 0),
        largegrid.getNearestLowerGridPoint(1.5, 1),
        largegrid.getNearestLowerGridPoint(1.5, 2) };
    largegrid.setWindow(window_begin, window_end);
    largegrid.addOntoWindow(window_begin, window_end, othervalues, +1.);
//    std::cout << " largegrid value " << largegrid.sampled_grid << std::endl;

    // smallgrid is downsample from full large grid
    SamplingGrid smallgrid(begin, end, base_level);
    SamplingGrid::downsample(smallgrid, largegrid, base_level);
//    std::cout << " smallgrid value " << smallgrid.sampled_grid << std::endl;

    // checkgrid is created in the same way as the large grid (only with smaller level)
    SamplingGrid checkgrid(begin, end, base_level);
    checkgrid.setWindow(window_begin, window_end);
    checkgrid.addOntoWindow(window_begin, window_end, checkvalues, +1.);
//    std::cout << " checkgrid value " << checkgrid.sampled_grid << std::endl;

    // then we compare over the full length
    /** Note that the threshold does not get better with increasing grid_difference,
     * on the contrary! For a grid_difference of 1, the finer grid is exact as it
     * directly sampled from the function. For a larger grid_difference, the finer
     * grid (being one level away from the desired coarse grid) is just an
     * approximation although obtained from a finer sampling. Hence, the larger the
     * grid_difference, the worse is the smallgrid with respect to checkgrid.
     */
    const double threshold = 3.2e-2;
    // we skip the first value as checkvalues[0] is always zero, hard to get from downsampling
    for (N[0]=1; N[0]< window_length_c; ++N[0])
      for (N[1]=0; N[1]< window_length_c; ++N[1])
        for (N[2]=0; N[2]< window_length_c; ++N[2]) {
          const double check_threshold =
              threshold*(1.+grid_difference/4.)*getBoundaryCaseFactor(N, window_length_c);
          const int index = calculateIndex(N, window_length_c);
          const double difference = fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index]));
//          std::cout << "Comparing |"
//              << scientific
//              << smallgrid.sampled_grid[index] << " - "
//              << checkgrid.sampled_grid[index] << "| = "
//              << fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index]))
//              << " < " << check_threshold << " => "
//              << fixed << setprecision(2)
//              << 100*difference/threshold
//              << " %, compare to " << 100*(1.+grid_difference/4.)*getBoundaryCaseFactor(N, window_length_c)
//              << " %" << std::endl;
          CPPUNIT_ASSERT( difference < check_threshold );
        }
  }
}