/*
 * LinearInterpolationBetweenSteps.hpp
 *
 *  Created on: Feb 23, 2011
 *      Author: heber
 */

#ifndef LINEARINTERPOLATIONBETWEENSTEPS_HPP_
#define LINEARINTERPOLATIONBETWEENSTEPS_HPP_

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

#include <vector>

#include "atom.hpp"
#include "AtomSet.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Dynamics/MinimiseConstrainedPotential.hpp"
#include "molecule.hpp"
#include "parser.hpp"
#include "World.hpp"

template <class Set>
class LinearInterpolationBetweenSteps
{
public:
  LinearInterpolationBetweenSteps(AtomSetMixin<Set> &_atoms, unsigned int _MaxOuterStep) :
    MDSteps(0),
    MaxOuterStep(_MaxOuterStep),
    IsAngstroem(true),
    atoms(_atoms)
  {}
  ~LinearInterpolationBetweenSteps()
  {}

  /** Performs a linear interpolation between two desired atomic configurations with a given number of steps.
   * Note, step number is config::MaxOuterStep
   * \param *out output stream for debugging
   * \param startstep stating initial configuration in molecule::Trajectories
   * \param endstep stating final configuration in molecule::Trajectories
   * \param &prefix path and prefix
   * \param MapByIdentity if true we just use the identity to map atoms in start config to end config, if not we find mapping by \sa MinimiseConstrainedPotential()
   * \return true - success in writing step files, false - error writing files or only one step in molecule::Trajectories
   */
  bool operator()(int startstep, int endstep, std::string prefix, bool MapByIdentity)
  {
    // TODO: rewrite permutationMaps using enumeration objects
    molecule *mol = NULL;
    bool status = true;
    int MaxSteps = MaxOuterStep;
    MoleculeListClass *MoleculePerStep = new MoleculeListClass(World::getPointer());
    // Get the Permutation Map by MinimiseConstrainedPotential
    atom **PermutationMap = NULL;
    atom *Sprinter = NULL;
    if (!MapByIdentity) {
      molecule::atomSet atoms_list;
      copy(atoms.begin(), atoms.end(), atoms_list.begin());
      MinimiseConstrainedPotential Minimiser(atoms_list);
      Minimiser(PermutationMap, startstep, endstep, IsAngstroem);
    } else {
      // TODO: construction of enumeration goes here
      PermutationMap = new atom *[atoms.size()];
      for(typename AtomSetMixin<Set>::const_iterator iter = atoms.begin(); iter != atoms.end();++iter){
        PermutationMap[(*iter)->getNr()] = (*iter);
      }
    }

    // check whether we have sufficient space in Trajectories for each atom
    for(typename AtomSetMixin<Set>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
      (*iter)->ResizeTrajectory(MaxSteps);
    }
    // push endstep to last one
    for(typename AtomSetMixin<Set>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
      (*iter)->CopyStepOnStep(MaxSteps,endstep);
    }
    endstep = MaxSteps;

    // go through all steps and add the molecular configuration to the list and to the Trajectories of \a this molecule
    DoLog(1) && (Log() << Verbose(1) << "Filling intermediate " << MaxSteps << " steps with MDSteps of " << MDSteps << "." << endl);
    for (int step = 0; step <= MaxSteps; step++) {
      mol = World::getInstance().createMolecule();
      MoleculePerStep->insert(mol);
      for (typename AtomSetMixin<Set>::const_iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
        // add to molecule list
        Sprinter = mol->AddCopyAtom((*iter));
        // add to Trajectories
        Vector temp = (*iter)->getPositionAtStep(startstep) + (PermutationMap[(*iter)->getNr()]->getPositionAtStep(endstep) - (*iter)->getPositionAtStep(startstep))*((double)step/(double)MaxSteps);
        Sprinter->setPosition(temp);
        (*iter)->setAtomicVelocityAtStep(step, zeroVec);
        (*iter)->setAtomicForceAtStep(step, zeroVec);
        //Log() << Verbose(3) << step << ">=" << MDSteps-1 << endl;
      }
    }
    MDSteps = MaxSteps+1;   // otherwise new Trajectories' points aren't stored on save&exit

    // store the list to single step files
    int *SortIndex = new int[atoms.size()];
    for (int i=atoms.size(); i--; )
      SortIndex[i] = i;

    status = MoleculePerStep->OutputConfigForListOfFragments(prefix, SortIndex);
    delete[](SortIndex);

    // free and return
    delete[](PermutationMap);
    delete(MoleculePerStep);
    return status;
  };

private:
  unsigned int MDSteps;
  unsigned int MaxOuterStep;
  bool IsAngstroem;
  AtomSetMixin<Set> atoms;

};

#endif /* LINEARINTERPOLATIONBETWEENSTEPS_HPP_ */