source: src/Dynamics/ForceAnnealing.hpp@ 8f50b4

/*
 * 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 <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 "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);

    // 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);

      if (_UseBondgraph)
        annealWithBondGraph(_CurrentTimeStep, _offset, maxComponents);
      else
        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)
  {
    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]));

      // are we in initial step? Then don't check against velocity
      if ((currentStep > 1) && (!(*iter)->getAtomicVelocity().IsZero()))
        // update with currentDelta tells us how the current gradient relates to
        // the last one: If it has become larger, reduce currentDelta
        if ((PositionUpdate.ScalarProduct((*iter)->getAtomicVelocity()) < 0)
            && (currentDeltat > MinimumDeltat)) {
          currentDeltat = .5*currentDeltat;
          LOG(2, "DEBUG: Upgrade in other direction: " << PositionUpdate.NormSquared()
              << " > " << (*iter)->getAtomicVelocity().NormSquared()
              << ", decreasing deltat: " << currentDeltat);
          PositionUpdate = currentDeltat * currentGradient;
      }
      // finally set new values
      (*iter)->setPosition(currentPosition + PositionUpdate);
      (*iter)->setAtomicVelocity(PositionUpdate);
      //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
//        (*iter)->VelocityVerletUpdateU((*iter)->getId(), CurrentTimeStep-1, Deltat, IsAngstroem);
    }
  }

  /** 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. t where \f$ t + \Delta t \f$ is 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
//    std::vector<atomId_t> atomids;
//    for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
//        iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
//      atomids.push_back((*iter)->getId());
//    }
//    ASSERT( atomids.size() == AtomicForceManipulator<T>::atoms.size(),
//        "ForceAnnealing() - could not gather all atomic ids?");
    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 more 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.

    // gather weights
    typedef std::deque<double> weights_t;
    typedef std::map<atomId_t, weights_t > weights_per_atom_t;
    std::vector<weights_per_atom_t> weights_per_atom(2);
    for (size_t timestep = 0; timestep <= 1; ++timestep) {
      const size_t CurrentStep = CurrentTimeStep-2*timestep >= 0 ? CurrentTimeStep - 2*timestep : 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);

        if (walkerGradient.Norm() > MYEPSILON) {

          // gather BondVector and projected gradient over all bonds
          const BondList& ListOfBonds = walker.getListOfBondsAtStep(CurrentStep);
          std::vector<double> projected_forces;
          std::vector<Vector> BondVectors;
          projected_forces.reserve(ListOfBonds.size());
          for(BondList::const_iterator bonditer = ListOfBonds.begin();
              bonditer != ListOfBonds.end(); ++bonditer) {
            const bond::ptr &current_bond = *bonditer;
            BondVectors.push_back(
                current_bond->leftatom->getPositionAtStep(CurrentStep)
                    - current_bond->rightatom->getPositionAtStep(CurrentStep));
            Vector &BondVector = BondVectors.back();
            BondVector.Normalize();
            projected_forces.push_back( walkerGradient.ScalarProduct(BondVector) );
          }

          // go through all bonds and check what magnitude is represented by the others
          // i.e. sum of scalar products against other bonds
          std::pair<weights_per_atom_t::iterator, bool> inserter =
              weights_per_atom[timestep-1].insert(
                  std::make_pair(walker.getId(), weights_t()) );
          ASSERT( inserter.second,
              "ForceAnnealing::operator() - weight map for atom "+toString(walker)
              +" and time step "+toString(timestep-1)+" already filled?");
          weights_t &weights = inserter.first->second;
          for (std::vector<Vector>::const_iterator iter = BondVectors.begin();
              iter != BondVectors.end(); ++iter) {
            std::vector<double> scps;
            std::transform(
                BondVectors.begin(), BondVectors.end(),
                std::back_inserter(scps),
                boost::bind(&Vector::ScalarProduct, boost::cref(*iter), _1)
                );
            const double scp_sum = std::accumulate(scps.begin(), scps.end(), 0.);
            weights.push_back( 1./scp_sum );
          }
          // for testing we check whether all weighted scalar products now come out as 1.
#ifndef NDEBUG
          for (std::vector<Vector>::const_iterator iter = BondVectors.begin();
              iter != BondVectors.end(); ++iter) {
            double scp_sum = 0.;
            weights_t::const_iterator weightiter = weights.begin();
            for (std::vector<Vector>::const_iterator otheriter = BondVectors.begin();
                otheriter != BondVectors.end(); ++otheriter, ++weightiter) {
              scp_sum += (*weightiter)*(*iter).ScalarProduct(*otheriter);
            }
            ASSERT( fabs(scp_sum - 1.) < MYEPSILON,
                "ForceAnnealing::operator() - for BondVector "+toString(*iter)
                +" we have weighted scalar product of "+toString(scp_sum)+" != 1.");
          }
#endif
        } else {
          LOG(2, "DEBUG: Gradient is " << walkerGradient << " less than "
              << MYEPSILON << " for atom " << walker);
        }
      }
    }

    // step through each bond and shift the atoms
    std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
//    for (size_t i = 0;i<bondvector.size();++i) {
//      const atom* bondatom[2] = {
//           bondvector[i]->leftatom,
//           bondvector[i]->rightatom};
//      const double &bondforce = bondforces[i];
//      const double &oldbondforce = oldbondforces[i];
//      const double bondforcedifference = (bondforce - oldbondforce);
//      Vector BondVector = (bondatom[0]->getPosition() - bondatom[1]->getPosition());
//      BondVector.Normalize();
//      double stepwidth = 0.;
//      for (size_t n=0;n<2;++n) {
//        const Vector oldPosition = bondatom[n]->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
//        const Vector currentPosition = bondatom[n]->getPosition();
//        stepwidth += fabs((currentPosition - oldPosition).ScalarProduct(BondVector))/bondforcedifference;
//      }
//      stepwidth = stepwidth/2;
//      Vector PositionUpdate = stepwidth * BondVector;
//      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.");
//        PositionUpdate = currentDeltat * BondVector;
//      }
//      LOG(3, "DEBUG: Update would be " << PositionUpdate);
//
//      // remove the edge
//#ifndef NDEBUG
//      const bool status =
//#endif
//          BGcreator.removeEdge(bondatom[0]->getId(), bondatom[1]->getId());
//      ASSERT( status, "ForceAnnealing() - edge to found bond is not present?");
//
//      // gather nodes for either atom
//      BoostGraphHelpers::Nodeset_t bondside_set[2];
//      BreadthFirstSearchGatherer::distance_map_t distance_map[2];
//      for (size_t n=0;n<2;++n) {
//        bondside_set[n] = NodeGatherer(bondatom[n]->getId(), max_distance);
//        distance_map[n] = NodeGatherer.getDistances();
//        std::sort(bondside_set[n].begin(), bondside_set[n].end());
//      }
//
//      // re-add edge
//      BGcreator.addEdge(bondatom[0]->getId(), bondatom[1]->getId());
//
//      // add PositionUpdate for all nodes in the bondside_set
//      for (size_t n=0;n<2;++n) {
//        for (BoostGraphHelpers::Nodeset_t::const_iterator setiter = bondside_set[n].begin();
//            setiter != bondside_set[n].end(); ++setiter) {
//          const BreadthFirstSearchGatherer::distance_map_t::const_iterator diter
//            = distance_map[n].find(*setiter);
//          ASSERT( diter != distance_map[n].end(),
//              "ForceAnnealing() - could not find distance to an atom.");
//          const double factor = pow(damping_factor, diter->second);
//          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) );
//          }
//        }
//        // invert for other atom
//        PositionUpdate *= -1;
//      }
//    }
//    delete[] bondforces;
//    delete[] oldbondforces;

    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.getAtomicForce();
      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
    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->getPosition() + update );
    }
  }

  /** 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_ */