Changeset fcc860 for src/Dynamics/ForceAnnealing.hpp

Timestamp:

Aug 3, 2017, 10:44:57 AM (8 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

ForceAnnealing_with_BondGraph_continued

Children:

4919f0

Parents:

93fd2a6

git-author:

Frederik Heber <frederik.heber@…> (08/03/17 09:24:07)

git-committer:

Frederik Heber <frederik.heber@…> (08/03/17 10:44:57)

Message:

ForceAnnealing functions now have better readable time step variables.

• _TimeStep is the parameter, Old.. and CurrentTimeStep designate the two reference time steps for calculating step width.
• Split functions into simple step width and using BB step width.

File:

• src/Dynamics/ForceAnnealing.hpp (modified) (16 diffs)

src/Dynamics/ForceAnnealing.hpp

@@ -92 +92 @@
    *
    *
-   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \param _TimeStep 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 int _TimeStep,
       const size_t _offset,
       const bool _UseBondgraph)
   {
+    const int CurrentTimeStep = _TimeStep-1;
+    ASSERT( CurrentTimeStep >= 0,
+        "ForceAnnealing::operator() - a new time step (upon which we work) must already have been copied.");
+
     // make sum of forces equal zero
     AtomicForceManipulator<T>::correctForceMatrixForFixedCenterOfMass(
         _offset,
-        _CurrentTimeStep-1>=0 ? _CurrentTimeStep - 1 : 0);
+        CurrentTimeStep);
 
     // are we in initial step? Then set static entities
@@ -114 +118 @@
 
       // always use atomic annealing on first step
-      maxComponents = anneal(_CurrentTimeStep);
+      maxComponents = anneal(_TimeStep);
     } else {
       ++currentStep;
@@ -121 +125 @@
       // bond graph annealing is always followed by a normal annealing
       if (_UseBondgraph)
-        maxComponents = annealWithBondGraph(_CurrentTimeStep);
+        maxComponents = annealWithBondGraph_BarzilaiBorwein(_TimeStep);
       // cannot store RemnantGradient in Atom's Force as it ruins BB stepwidth calculation
       else
-        maxComponents = anneal(_CurrentTimeStep);
+        maxComponents = anneal_BarzilaiBorwein(_TimeStep);
     }
 
@@ -141 +145 @@
    * 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 _TimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
    * \return to be filled with maximum force component over all atoms
    */
   Vector anneal(
-      const int CurrentTimeStep)
+      const int _TimeStep)
   {
+    const int CurrentTimeStep = _TimeStep-1;
+    ASSERT( CurrentTimeStep >= 0,
+        "ForceAnnealing::anneal() - a new time step (upon which we work) must already have been copied.");
+
+    LOG(1, "STATUS: performing simple anneal with default stepwidth " << currentDeltat << " at step #" << currentStep);
+
     Vector maxComponents;
     bool deltat_decreased = false;
@@ -152 +162 @@
                     iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
                         // atom's force vector gives steepest descent direction
-                        const Vector oldPosition = (*iter)->getPositionAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
-                        const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0);
-                        const Vector oldGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-2 >= 0 ? CurrentTimeStep - 2 : 0);
-                        const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0);
-      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);
+                        const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
+                        const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
+      LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition);
+      LOG(4, "DEBUG: currentGradient for atom #" << (*iter)->getId() << " is " << currentGradient);
 //                      LOG(4, "DEBUG: Force for atom " << **iter << " is " << currentGradient);
 
       // we use Barzilai-Borwein update with position reversed to get descent
-      const Vector PositionDifference = currentPosition - oldPosition;
-                        const Vector GradientDifference = (currentGradient - oldGradient);
-                        double stepwidth = 0.;
-                        if (GradientDifference.Norm() > 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;
-                        }
+                        double stepwidth = currentDeltat;
       Vector PositionUpdate = stepwidth * currentGradient;
                         LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
@@ -186 +181 @@
                         // have different sign: Check whether this is the case and one step with
                         // deltat to interrupt this sequence
-                        if ((currentStep > 1) && (!PositionDifference.IsZero()))
+                        if (currentStep > 1) {
+                    const int OldTimeStep = _TimeStep-2;
+                    ASSERT( OldTimeStep >= 0,
+                        "ForceAnnealing::anneal() - if currentStep is "+toString(currentStep)
+                        +", then there should be at least three time steps.");
+              const Vector oldPosition = (*iter)->getPositionAtStep(OldTimeStep);
+              const Vector PositionDifference = currentPosition - oldPosition;
+        LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
+        LOG(4, "DEBUG: PositionDifference for atom #" << (*iter)->getId() << " is " << PositionDifference);
                                 if ((PositionUpdate.ScalarProduct(PositionDifference) < 0)
                                     && (fabs(PositionUpdate.NormSquared()-PositionDifference.NormSquared()) < 1e-3)) {
@@ -198 +201 @@
                                       << ", using deltat: " << currentDeltat);
                                   PositionUpdate = currentDeltat * currentGradient;
+                                }
                         }
 
                         // finally set new values
-                        (*iter)->setPosition(currentPosition + PositionUpdate);
+                        (*iter)->setPositionAtStep(_TimeStep, currentPosition + PositionUpdate);
                 }
+
+    return maxComponents;
+  }
+
+  /** Performs Gradient optimization on the atoms using BarzilaiBorwein step width.
+   *
+   * \note this can only be called when there are at least two optimization
+   * time steps present, i.e. this must be preceeded by a simple anneal().
+   *
+   * We assume that forces have just been calculated.
+   *
+   * \param _TimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \return to be filled with maximum force component over all atoms
+   */
+  Vector anneal_BarzilaiBorwein(
+      const int _TimeStep)
+  {
+    const int OldTimeStep = _TimeStep-2;
+    const int CurrentTimeStep = _TimeStep-1;
+    ASSERT( OldTimeStep >= 0,
+        "ForceAnnealing::anneal_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
+    ASSERT(currentStep > 1,
+        "ForceAnnealing::anneal_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
+
+    LOG(1, "STATUS: performing BarzilaiBorwein anneal at step #" << currentStep);
+
+    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(OldTimeStep);
+      const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
+      const Vector oldGradient = (*iter)->getAtomicForceAtStep(OldTimeStep);
+      const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
+      LOG(4, "DEBUG: oldPosition for atom #" << (*iter)->getId() << " is " << oldPosition);
+      LOG(4, "DEBUG: currentPosition for atom #" << (*iter)->getId() << " is " << currentPosition);
+      LOG(4, "DEBUG: oldGradient for atom #" << (*iter)->getId() << " is " << oldGradient);
+      LOG(4, "DEBUG: currentGradient for atom #" << (*iter)->getId() << " is " << currentGradient);
+//      LOG(4, "DEBUG: Force for atom #" << (*iter)->getId() << " is " << currentGradient);
+
+      // we use Barzilai-Borwein update with position reversed to get descent
+      const Vector PositionDifference = currentPosition - oldPosition;
+      const Vector GradientDifference = (currentGradient - oldGradient);
+      double stepwidth = 0.;
+      if (GradientDifference.Norm() > 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;
+      }
+      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)
+        if (currentGradient[i] > maxComponents[i])
+          maxComponents[i] = 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
+      if (!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)->setPositionAtStep(_TimeStep, currentPosition + PositionUpdate);
+    }
 
     return maxComponents;
@@ -244 +330 @@
   }
 
-  /** Performs Gradient optimization on the bonds.
+  /** Performs Gradient optimization on the bonds with BarzilaiBorwein stepwdith.
+   *
+   * \note this can only be called when there are at least two optimization
+   * time steps present, i.e. this must be preceeded by a simple anneal().
    *
    * We assume that forces have just been calculated. These forces are projected
@@ -251 +340 @@
    *
    *
-   * \param CurrentTimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   * \param _TimeStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
    * \param maxComponents to be filled with maximum force component over all atoms
    */
-  Vector annealWithBondGraph(
-      const int CurrentTimeStep)
+  Vector annealWithBondGraph_BarzilaiBorwein(
+      const int _TimeStep)
   {
+    const int OldTimeStep = _TimeStep-2;
+    const int CurrentTimeStep = _TimeStep-1;
+    ASSERT(currentStep > 1,
+        "annealWithBondGraph_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth.");
+    ASSERT(OldTimeStep >= 0,
+        "annealWithBondGraph_BarzilaiBorwein() - we need two present optimization steps to compute stepwidth, and the new one to update on already present.");
+
+    LOG(1, "STATUS: performing BarzilaiBorwein anneal on bonds at step #" << currentStep);
+
     Vector maxComponents;
 
@@ -297 +395 @@
         AtomicForceManipulator<T>::atoms.begin(),
         AtomicForceManipulator<T>::atoms.end(),
-        CurrentTimeStep);
+        _TimeStep);
     const BondVectors::container_t &sorted_bonds = bv.getSorted();
 
@@ -316 +414 @@
     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
+      const size_t ReferenceTimeStep = _TimeStep-timestep-1;
+      LOG(2, "DEBUG: given time step is " << _TimeStep
           << ", timestep is " << timestep
-          << ", and CurrentStep is " << CurrentStep);
+          << ", and ReferenceTimeStep is " << ReferenceTimeStep);
 
       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);
+        const Vector &walkerGradient = walker.getAtomicForceAtStep(ReferenceTimeStep);
         LOG(3, "DEBUG: Gradient of atom #" << walker.getId() << ", namely "
             << walker << " is " << walkerGradient);
@@ -333 +431 @@
           // gather subset of BondVectors for the current atom
           const std::vector<Vector> BondVectors =
-              bv.getAtomsBondVectorsAtStep(walker, CurrentStep);
+              bv.getAtomsBondVectorsAtStep(walker, ReferenceTimeStep);
 
           // go through all its bonds and calculate what magnitude is represented
@@ -340 +438 @@
               weights_per_atom[timestep].insert(
                   std::make_pair(walker.getId(),
-                  bv.getWeightsForAtomAtStep(walker, BondVectors, CurrentStep)) );
+                  bv.getWeightsForAtomAtStep(walker, BondVectors, ReferenceTimeStep)) );
           ASSERT( inserter.second,
               "ForceAnnealing::operator() - weight map for atom "+toString(walker)
@@ -428 +526 @@
         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 &oldPosition = bondatom[side]->getPositionAtStep(OldTimeStep);
+        const Vector &currentPosition = bondatom[side]->getPositionAtStep(CurrentTimeStep);
         const Vector PositionDifference = currentPosition - oldPosition;
         LOG(4, "DEBUG: old position is " << oldPosition);
@@ -478 +576 @@
       atom &walker = *(*iter);
       // extract largest components for showing progress of annealing
-      const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep-1 : 0);
+      const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep);
       for(size_t i=0;i<NDIM;++i)
         if (currentGradient[i] > maxComponents[i])
@@ -505 +603 @@
       LOG(3, "DEBUG: Applying update " << update << " to atom #" << atomid
           << ", namely " << *walker);
-      walker->setPosition(
-          walker->getPositionAtStep(CurrentTimeStep-1>=0 ? CurrentTimeStep - 1 : 0)
-          + update - CommonTranslation);
+      walker->setPositionAtStep(_TimeStep,
+          walker->getPositionAtStep(CurrentTimeStep) + update - CommonTranslation);
 //      walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
     }