source: src/Dynamics/ForceAnnealing.hpp@ f433ec

/*
 * ForceAnnealing.hpp
 *
 *  Created on: Aug 02, 2014
 *      Author: heber
 */

#ifndef FORCEANNEALING_HPP_
#define FORCEANNEALING_HPP_

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

#include <algorithm>
#include <functional>
#include <iterator>

#include <boost/bind.hpp>

#include "Atom/atom.hpp"
#include "Atom/AtomSet.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "Dynamics/AtomicForceManipulator.hpp"
#include "Dynamics/BondVectors.hpp"
#include "Fragmentation/ForceMatrix.hpp"
#include "Graph/BoostGraphCreator.hpp"
#include "Graph/BoostGraphHelpers.hpp"
#include "Graph/BreadthFirstSearchGatherer.hpp"
#include "Helpers/helpers.hpp"
#include "Helpers/defs.hpp"
#include "LinearAlgebra/LinearSystemOfEquations.hpp"
#include "LinearAlgebra/MatrixContent.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/VectorContent.hpp"
#include "Thermostats/ThermoStatContainer.hpp"
#include "Thermostats/Thermostat.hpp"
#include "World.hpp"

/** This class is the essential build block for performing structural optimization.
 *
 * Sadly, we have to use some static instances as so far values cannot be passed
 * between actions. Hence, we need to store the current step and the adaptive-
 * step width (we cannot perform a line search, as we have no control over the
 * calculation of the forces).
 *
 * However, we do use the bond graph, i.e. if a single atom needs to be shifted
 * to the left, then the whole molecule left of it is shifted, too. This is
 * controlled by the \a max_distance parameter.
 */
template <class T>
class ForceAnnealing : public AtomicForceManipulator<T>
{
public:
  /** Constructor of class ForceAnnealing.
   *
   * \note We use a fixed delta t of 1.
   *
   * \param _atoms set of atoms to integrate
   * \param _Deltat time step width in atomic units
   * \param _IsAngstroem whether length units are in angstroem or bohr radii
   * \param _maxSteps number of optimization steps to perform
   * \param _max_distance up to this bond order is bond graph taken into account.
   */
  ForceAnnealing(
      AtomSetMixin<T> &_atoms,
      const double _Deltat,
      bool _IsAngstroem,
      const size_t _maxSteps,
      const int _max_distance,
      const double _damping_factor) :
    AtomicForceManipulator<T>(_atoms, _Deltat, _IsAngstroem),
    maxSteps(_maxSteps),
    max_distance(_max_distance),
    damping_factor(_damping_factor)
  {}

  /** Destructor of class ForceAnnealing.
   *
   */
  ~ForceAnnealing()
  {}

  /** Performs Gradient optimization.
   *
   * We assume that forces have just been calculated.
   *
   *
   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
   * \param offset offset in matrix file to the first force component
   * \todo This is not yet checked if it is correctly working with DoConstrainedMD set >0.
   */
  void operator()(
      const int _CurrentTimeStep,
      const size_t _offset,
      const bool _UseBondgraph)
  {
    // make sum of forces equal zero
    AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(
        _offset,
        _CurrentTimeStep-1>=0 ? _CurrentTimeStep - 1 : 0);

    // are we in initial step? Then set static entities
    Vector maxComponents(zeroVec);
    if (currentStep == 0) {
      currentDeltat = AtomicForceManipulator<T>::Deltat;
      currentStep = 1;
      LOG(2, "DEBUG: Initial step, setting values, current step is #" << currentStep);

      // always use atomic annealing on first step
      anneal(_CurrentTimeStep, _offset, maxComponents);
    } else {
      ++currentStep;
      LOG(2, "DEBUG: current step is #" << currentStep);

      // bond graph annealing is always followed by a normal annealing
      if (_UseBondgraph)
        annealWithBondGraph(_CurrentTimeStep, _offset, maxComponents);
      anneal(_CurrentTimeStep, _offset, maxComponents);
    }


    LOG(1, "STATUS: Largest remaining force components at step #"
        << currentStep << " are " << maxComponents);

    // are we in final step? Remember to reset static entities
    if (currentStep == maxSteps) {
      LOG(2, "DEBUG: Final step, resetting values");
      reset();
    }
  }

  /** Helper function to calculate the Barzilai-Borwein stepwidth.
   *
   * \param _PositionDifference difference in position between current and last step
   * \param _GradientDifference difference in gradient between current and last step
   * \return step width according to Barzilai-Borwein
   */
  double getBarzilaiBorweinStepwidth(const Vector &_PositionDifference, const Vector &_GradientDifference)
  {
    double stepwidth = 0.;
    if (_GradientDifference.NormSquared() > MYEPSILON)
      stepwidth = fabs(_PositionDifference.ScalarProduct(_GradientDifference))/
          _GradientDifference.NormSquared();
    if (fabs(stepwidth) < 1e-10) {
      // dont' warn in first step, deltat usage normal
      if (currentStep != 1)
        ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
      stepwidth = currentDeltat;
    }
    return stepwidth;
  }

  /** Performs Gradient optimization on the atoms.
   *
   * We assume that forces have just been calculated.
   *
   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
   * \param offset offset in matrix file to the first force component
   * \param maxComponents to be filled with maximum force component over all atoms
   */
  void anneal(
      const int CurrentTimeStep,
      const size_t offset,
      Vector &maxComponents)
  {
    bool deltat_decreased = false;
    for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
      // atom's force vector gives steepest descent direction
      const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
      const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
      const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1 >= 0 ? CurrentTimeStep - 1 : 0);
      const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
      LOG(4, "DEBUG: oldPosition for atom " << **iter << " is " << oldPosition);
      LOG(4, "DEBUG: currentPosition for atom " << **iter << " is " << currentPosition);
      LOG(4, "DEBUG: oldGradient for atom " << **iter << " is " << oldGradient);
      LOG(4, "DEBUG: currentGradient for atom " << **iter << " is " << currentGradient);
//      LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);

      // we use Barzilai-Borwein update with position reversed to get descent
      const double stepwidth = getBarzilaiBorweinStepwidth(
          currentPosition - oldPosition, currentGradient - oldGradient);
      Vector PositionUpdate = stepwidth * currentGradient;
      LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);

      // extract largest components for showing progress of annealing
      for(size_t i=0;i<NDIM;++i)
        maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));

      // steps may go back and forth again (updates are of same magnitude but
      // have different sign: Check whether this is the case and one step with
      // deltat to interrupt this sequence
      const Vector PositionDifference = currentPosition - oldPosition;
      if ((currentStep > 1) && (!PositionDifference.IsZero()))
        if ((PositionUpdate.ScalarProduct(PositionDifference) < 0)
            && (fabs(PositionUpdate.NormSquared()-PositionDifference.NormSquared()) < 1e-3)) {
          // for convergence we want a null sequence here, too
          if (!deltat_decreased) {
            deltat_decreased = true;
            currentDeltat = .5*currentDeltat;
          }
          LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate
              << " > " << PositionDifference
              << ", using deltat: " << currentDeltat);
          PositionUpdate = currentDeltat * currentGradient;
      }

      // finally set new values
      (*iter)->setPosition(currentPosition + PositionUpdate);
    }
  }

  // knowing the number of bonds in total, we can setup the storage for the
  // projected forces
  enum whichatom_t {
    leftside=0,
    rightside=1,
    MAX_sides
  };

  /** Helper function to put bond force into a container.
   *
   * \param _walker atom
   * \param _current_bond current bond of \a _walker
   * \param _timestep time step
   * \param _force calculated bond force
   * \param _bv bondvectors for obtaining the correct index
   * \param _projected_forces container
   */
  static void ForceStore(
      const atom &_walker,
      const bond::ptr &_current_bond,
      const size_t &_timestep,
      const double _force,
      const BondVectors &_bv,
      std::vector< // time step
            std::vector< // which bond side
              std::vector<double> > // over all bonds
                > &_projected_forces)
  {
    std::vector<double> &forcelist = (&_walker == _current_bond->leftatom) ?
        _projected_forces[_timestep][leftside] : _projected_forces[_timestep][rightside];
    const size_t index = _bv.getIndexForBond(_current_bond);
    ASSERT( index != (size_t)-1,
        "ForceAnnealing() - could not find bond "+toString(*_current_bond)
        +" in bondvectors");
    forcelist[index] = _force;
  }

  /** Performs Gradient optimization on the bonds.
   *
   * We assume that forces have just been calculated. These forces are projected
   * onto the bonds and these are annealed subsequently by moving atoms in the
   * bond neighborhood on either side conjunctively.
   *
   *
   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
   * \param offset offset in matrix file to the first force component
   * \param maxComponents to be filled with maximum force component over all atoms
   */
  void annealWithBondGraph(
      const int CurrentTimeStep,
      const size_t offset,
      Vector &maxComponents)
  {
    // get nodes on either side of selected bond via BFS discovery
    BoostGraphCreator BGcreator;
    BGcreator.createFromRange(
        AtomicForceManipulator<T>::atoms.begin(),
        AtomicForceManipulator<T>::atoms.end(),
        AtomicForceManipulator<T>::atoms.size(),
        BreadthFirstSearchGatherer::AlwaysTruePredicate);
    BreadthFirstSearchGatherer NodeGatherer(BGcreator);

    /// We assume that a force is local, i.e. a bond is too short yet and hence
    /// the atom needs to be moved. However, all the adjacent (bound) atoms might
    /// already be at the perfect distance. If we just move the atom alone, we ruin
    /// all the other bonds. Hence, it would be sensible to move every atom found
    /// through the bond graph in the direction of the force as well by the same
    /// PositionUpdate. This is almost what we are going to do.

    /// One issue is: If we need to shorten bond, then we use the PositionUpdate
    /// also on the the other bond partner already. This is because it is in the
    /// direction of the bond. Therefore, the update is actually performed twice on
    /// each bond partner, i.e. the step size is twice as large as it should be.
    /// This problem only occurs when bonds need to be shortened, not when they
    /// need to be made longer (then the force vector is facing the other
    /// direction than the bond vector).
    /// As a remedy we need to average the force on either end of the bond and
    /// check whether each gradient points inwards out or outwards with respect
    /// to the bond and then shift accordingly.

    /// One more issue is that the projection onto the bond directions does not
    /// recover the gradient but may be larger as the bond directions are a
    /// generating system and not a basis (e.g. 3 bonds on a plane where 2 would
    /// suffice to span the plane). To this end, we need to account for the
    /// overestimation and obtain a weighting for each bond.

    // initialize helper class for bond vectors using bonds from range of atoms
    BondVectors bv;
    bv.setFromAtomRange< T >(
        AtomicForceManipulator<T>::atoms.begin(),
        AtomicForceManipulator<T>::atoms.end(),
        CurrentTimeStep);
    const BondVectors::container_t &sorted_bonds = bv.getSorted();

    std::vector< // time step
      std::vector< // which bond side
        std::vector<double> > // over all bonds
          > projected_forces(2); // one for leftatoms, one for rightatoms (and for both time steps)
    for (size_t i=0;i<2;++i) {
      projected_forces[i].resize(MAX_sides);
      for (size_t j=0;j<MAX_sides;++j)
        projected_forces[i][j].resize(sorted_bonds.size(), 0.);
    }

    // for each atom we need to gather weights and then project the gradient
    typedef std::map<atomId_t, BondVectors::weights_t > weights_per_atom_t;
    std::vector<weights_per_atom_t> weights_per_atom(2);
    typedef std::map<atomId_t, Vector> RemnantGradient_per_atom_t;
    RemnantGradient_per_atom_t RemnantGradient_per_atom;
    for (size_t timestep = 0; timestep <= 1; ++timestep) {
      const size_t CurrentStep = CurrentTimeStep-timestep-1 >= 0 ? CurrentTimeStep-timestep-1 : 0;
      LOG(2, "DEBUG: CurrentTimeStep is " << CurrentTimeStep
          << ", timestep is " << timestep
          << ", and CurrentStep is " << CurrentStep);

      for(typename AtomSetMixin<T>::const_iterator iter = AtomicForceManipulator<T>::atoms.begin();
          iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
        const atom &walker = *(*iter);
        const Vector &walkerGradient = walker.getAtomicForceAtStep(CurrentStep);
        LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
            << walker << " is " << walkerGradient << " with magnitude of "
            << walkerGradient.Norm());

        const BondList& ListOfBonds = walker.getListOfBonds();
        if (walkerGradient.Norm() > MYEPSILON) {

          // gather subset of BondVectors for the current atom
          const std::vector<Vector> BondVectors =
              bv.getAtomsBondVectorsAtStep(walker, CurrentStep);

          // go through all its bonds and calculate what magnitude is represented
          // by the others i.e. sum of scalar products against other bonds
          const std::pair<weights_per_atom_t::iterator, bool> inserter =
              weights_per_atom[timestep].insert(
                  std::make_pair(walker.getId(),
                  bv.getWeightsForAtomAtStep(walker, BondVectors, CurrentStep)) );
          ASSERT( inserter.second,
              "ForceAnnealing::operator() - weight map for atom "+toString(walker)
              +" and time step "+toString(timestep)+" already filled?");
          BondVectors::weights_t &weights = inserter.first->second;
          ASSERT( weights.size() == ListOfBonds.size(),
              "ForceAnnealing::operator() - number of weights "
              +toString(weights.size())+" does not match number of bonds "
              +toString(ListOfBonds.size())+", error in calculation?");

          // projected gradient over all bonds and place in one of projected_forces
          // using the obtained weights
          BondVectors::forcestore_t forcestoring =
              boost::bind(&ForceAnnealing::ForceStore, _1, _2, _3, _4,
                  boost::cref(bv), boost::ref(projected_forces));
          const Vector RemnantGradient = bv.getRemnantGradientForAtomAtStep(
              walker, BondVectors, weights, timestep, forcestoring
          );
          RemnantGradient_per_atom.insert( std::make_pair(walker.getId(), RemnantGradient) );
        } else {
          LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
              << MYEPSILON << " for atom " << walker);
          // note that projected_forces is initialized to full length and filled
          // with zeros. Hence, nothing to do here
        }
      }
    }

    // step through each bond and shift the atoms
    std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end

    LOG(3, "DEBUG: current step is " << currentStep << ", given time step is " << CurrentTimeStep);
    const BondVectors::mapped_t bondvectors = bv.getBondVectorsAtStep(CurrentTimeStep);

    for (BondVectors::container_t::const_iterator bondsiter = sorted_bonds.begin();
        bondsiter != sorted_bonds.end(); ++bondsiter) {
      const bond::ptr &current_bond = *bondsiter;
      const size_t index = bv.getIndexForBond(current_bond);
      const atom* bondatom[MAX_sides] = {
          current_bond->leftatom,
          current_bond->rightatom
      };

      // remove the edge
#ifndef NDEBUG
      const bool status =
#endif
          BGcreator.removeEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());
      ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");

      // gather nodes for either atom
      BoostGraphHelpers::Nodeset_t bondside_set[MAX_sides];
      BreadthFirstSearchGatherer::distance_map_t distance_map[MAX_sides];
      for (size_t side=leftside;side<MAX_sides;++side) {
        bondside_set[side] = NodeGatherer(bondatom[side]->getId(), max_distance);
        distance_map[side] = NodeGatherer.getDistances();
        std::sort(bondside_set[side].begin(), bondside_set[side].end());
      }

      // re-add edge
      BGcreator.addEdge(bondatom[leftside]->getId(), bondatom[rightside]->getId());

      // do for both leftatom and rightatom of bond
      for (size_t side = leftside; side < MAX_sides; ++side) {
        const double &bondforce = projected_forces[0][side][index];
        const double &oldbondforce = projected_forces[1][side][index];
        const double bondforcedifference = fabs(bondforce - oldbondforce);
        LOG(4, "DEBUG: bondforce for " << (side == leftside ? "left" : "right")
            << " side of bond is " << bondforce);
        LOG(4, "DEBUG: oldbondforce for " << (side == leftside ? "left" : "right")
            << " side of bond is " << oldbondforce);
        // if difference or bondforce itself is zero, do nothing
        if ((fabs(bondforce) < MYEPSILON) || (fabs(bondforcedifference) < MYEPSILON))
          continue;

        // get BondVector to bond
        const BondVectors::mapped_t::const_iterator bviter =
            bondvectors.find(current_bond);
        ASSERT( bviter != bondvectors.end(),
            "ForceAnnealing() - cannot find current_bond ?");
        ASSERT( fabs(bviter->second.Norm() -1.) < MYEPSILON,
            "ForceAnnealing() - norm of BondVector is not one");
        const Vector &BondVector = bviter->second;

        // calculate gradient and position differences for stepwidth
        const Vector currentGradient = bondforce * BondVector;
        LOG(4, "DEBUG: current projected gradient for "
            << (side == leftside ? "left" : "right") << " side of bond is " << currentGradient);
        const Vector &oldPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
        const Vector &currentPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0);
        const Vector PositionDifference = currentPosition - oldPosition;
        LOG(4, "DEBUG: old position is " << oldPosition);
        LOG(4, "DEBUG: current position is " << currentPosition);
        LOG(4, "DEBUG: difference in position is " << PositionDifference);
        LOG(4, "DEBUG: bondvector is " << BondVector);
        const double projected_PositionDifference = PositionDifference.ScalarProduct(BondVector);
        LOG(4, "DEBUG: difference in position projected onto bondvector is "
            << projected_PositionDifference);
        LOG(4, "DEBUG: abs. difference in forces is " << bondforcedifference);

        // calculate step width
        double stepwidth =
            fabs(projected_PositionDifference)/bondforcedifference;
        if (fabs(stepwidth) < 1e-10) {
          // dont' warn in first step, deltat usage normal
          if (currentStep != 1)
            ELOG(1, "INFO: Barzilai-Borwein stepwidth is zero, using deltat " << currentDeltat << " instead.");
          stepwidth = currentDeltat;
        }
        Vector PositionUpdate = stepwidth * currentGradient;
        LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);

        // add PositionUpdate for all nodes in the bondside_set
        for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set[side].begin();
            setiter != bondside_set[side].end(); ++setiter) {
          const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
            = distance_map[side].find(*setiter);
          ASSERT( diter != distance_map[side].end(),
              "ForceAnnealing() - could not find distance to an atom.");
          const double factor = pow(damping_factor, diter->second+1);
          LOG(3, "DEBUG: Update for atom #" << *setiter << " will be "
              << factor << "*" << PositionUpdate);
          if (GatheredUpdates.count((*setiter))) {
            GatheredUpdates[(*setiter)] += factor*PositionUpdate;
          } else {
            GatheredUpdates.insert(
                std::make_pair(
                    (*setiter),
                    factor*PositionUpdate) );
          }
        }
      }
    }

    for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
      atom &walker = *(*iter);
      // extract largest components for showing progress of annealing
      const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep-1 : 0);
      for(size_t i=0;i<NDIM;++i)
        maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));

      // reset force vector for next step except on final one
      if (currentStep != maxSteps)
        walker.setAtomicForce(zeroVec);
    }

    // apply the gathered updates and set remnant gradients for atomic annealing
    for (std::map<atomId_t, Vector>::const_iterator iter = GatheredUpdates.begin();
        iter != GatheredUpdates.end(); ++iter) {
      const atomId_t &atomid = iter->first;
      const Vector &update = iter->second;
      atom* const walker = World::getInstance().getAtom(AtomById(atomid));
      ASSERT( walker != NULL,
          "ForceAnnealing() - walker with id "+toString(atomid)+" has suddenly disappeared.");
      LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
          << ", namely " << *walker);
      walker->setPosition(
          walker->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0)
          + update);
      walker->setAtomicVelocity(update);
      walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
    }
  }

  /** Reset function to unset static entities and artificial velocities.
   *
   */
  void reset()
  {
    currentDeltat = 0.;
    currentStep = 0;
  }

private:
  //!> contains the current step in relation to maxsteps
  static size_t currentStep;
  //!> contains the maximum number of steps, determines initial and final step with currentStep
  size_t maxSteps;
  static double currentDeltat;
  //!> minimum deltat for internal while loop (adaptive step width)
  static double MinimumDeltat;
  //!> contains the maximum bond graph distance up to which shifts of a single atom are spread
  const int max_distance;
  //!> the shifted is dampened by this factor with the power of the bond graph distance to the shift causing atom
  const double damping_factor;
};

template <class T>
double ForceAnnealing<T>::currentDeltat = 0.;
template <class T>
size_t ForceAnnealing<T>::currentStep = 0;
template <class T>
double ForceAnnealing<T>::MinimumDeltat = 1e-8;

#endif /* FORCEANNEALING_HPP_ */