/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * 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/>.
 */

/*
 * PartialNucleiChargeFitterUnitTest.cpp
 *
 *  Created on: May 12, 2013
 *      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 "PartialNucleiChargeFitterUnitTest.hpp"

#include <boost/assign.hpp>
#include <vector>

#include "CodePatterns/Assert.hpp"
#include "LinearAlgebra/MatrixContent.hpp"
#include "LinearAlgebra/Vector.hpp"

#include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
#include "Potentials/PartialNucleiChargeFitter.hpp"

using namespace boost::assign;

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

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

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

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


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

  // create the grid with domain [0:1]^3
  const double begin[3] = { 0., 0., 0. };
  const double end[3] = { 1., 1., 1. };
  SamplingGrid::sampledvalues_t values;
  for (size_t i=0; i< DOMAINVOLUME(1)*NUMBEROFSAMPLES(2); ++i)
    values += 0.;
  SamplingGrid grid(begin, end, 2, values);

  // create position(s) in center of domain
  PartialNucleiChargeFitter::positions_t positions(1, Vector(.5,.5,.5));
  fitter = new PartialNucleiChargeFitter(grid, positions);
}


void PartialNucleiChargeFitterTest::tearDown()
{
  delete fitter;
}

/** UnitTest for constructMatrix()
 */
void PartialNucleiChargeFitterTest::constructMatrixTest()
{
  // construct potential
  MatrixContent testMatrix(DOMAINVOLUME(1)*NUMBEROFSAMPLES(2), 1);
  Vector center(.5,.5,.5);
  const double delta = 1./pow(2,2);
  Vector N;
  size_t index=0;
  for(N[0]=0.; N[0] < 1.; N[0]+=delta)
    for(N[1]=0.; N[1] < 1.; N[1]+=delta)
      for(N[2]=0.; N[2] < 1.; N[2]+=delta) {
        const double distance = N.distance(center);
        if (fabs(distance) > std::numeric_limits<double>::epsilon()*1e4 )
          testMatrix.at(index++, 0) = 1./(distance);
        else
          testMatrix.at(index++, 0) = 0.;
      }

  // check whether potential matrix is fine
  fitter->constructMatrix();
  const MatrixContent &ProblemMatrix = fitter->getMatrix();
  CPPUNIT_ASSERT_EQUAL( testMatrix, ProblemMatrix );
}

/** UnitTest for operator()
 */
void PartialNucleiChargeFitterTest::operatorTest()
{
  // solution is zero everywhere
  PartialNucleiChargeFitter::charges_t charges(1, 0.);

  // true solution is ...
  (*fitter)();
  PartialNucleiChargeFitter::charges_t return_charges = fitter->getSolutionAsCharges_t();
  CPPUNIT_ASSERT_EQUAL(charges, return_charges);
}