Changeset b3aaf4 for src/Dynamics

Timestamp:

Jun 20, 2018, 8:21:13 AM (7 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

Candidate_v1.6.1, ChemicalSpaceEvaluator, Exclude_Hydrogens_annealWithBondGraph, ForceAnnealing_with_BondGraph_contraction-expansion

Children:

31baf0

Parents:

883d98

git-author:

Frederik Heber <frederik.heber@…> (12/05/17 09:10:11)

git-committer:

Frederik Heber <frederik.heber@…> (06/20/18 08:21:13)

Message:

Preventing cancelling updates for optimization with bond graph.

• we check norm of single updates against gathered update. If the latter is very small, we perform a regular optimization step.
• extracted method annealAtom_BB() to this end.
• Moreover, we cap the maximum BB update at 0.5 as it tends to overshoot.
• TESTS: This fixes the issue with 2-body Ising model Python regression test where this issue occurs.

File:

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

src/Dynamics/ForceAnnealing.hpp

@@ -17 +17 @@
 #include <functional>
 #include <iterator>
+#include <math.h>
 
 #include <boost/bind.hpp>
@@ -167 +168 @@
   {
     double stepwidth = 0.;
-    if (_GradientDifference.NormSquared() > MYEPSILON)
+    if (_GradientDifference.Norm() > MYEPSILON)
       stepwidth = fabs(_PositionDifference.ScalarProduct(_GradientDifference))/
           _GradientDifference.NormSquared();
@@ -200 +201 @@
         iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
       // atom's force vector gives steepest descent direction
-      const Vector currentPosition = (*iter)->getPositionAtStep(CurrentTimeStep);
-      const Vector currentGradient = (*iter)->getAtomicForceAtStep(CurrentTimeStep);
+      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);
@@ -248 +249 @@
   }
 
+  /** Performs Gradient optimization on a single atom using BarzilaiBorwein step width.
+   *
+   * \param _atom atom to anneal
+   * \param OldTimeStep old time step
+   * \param CurrentTimeStep current time step whose gradient we've just calculated
+   * \param TimeStepToSet time step to update (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
+   */
+  void annealAtom_BarzilaiBorwein(
+      atom * const _atom,
+      const int &OldTimeStep,
+      const int &CurrentTimeStep,
+      const int &TimeStepToSet
+      )
+  {
+    // atom's force vector gives steepest descent direction
+    const Vector &oldPosition = _atom->getPositionAtStep(OldTimeStep);
+    const Vector &currentPosition = _atom->getPositionAtStep(CurrentTimeStep);
+    const Vector &oldGradient = _atom->getAtomicForceAtStep(OldTimeStep);
+    const Vector &currentGradient = _atom->getAtomicForceAtStep(CurrentTimeStep);
+    LOG(4, "DEBUG: oldPosition for atom #" << _atom->getId() << " is " << oldPosition);
+    LOG(4, "DEBUG: currentPosition for atom #" << _atom->getId() << " is " << currentPosition);
+    LOG(4, "DEBUG: oldGradient for atom #" << _atom->getId() << " is " << oldGradient);
+    LOG(4, "DEBUG: currentGradient for atom #" << _atom->getId() << " is " << currentGradient);
+//      LOG(4, "DEBUG: Force for atom #" << _atom->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 = getBarzilaiBorweinStepwidth(PositionDifference, GradientDifference);
+    Vector PositionUpdate = stepwidth * currentGradient;
+    LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
+
+    // finally set new values
+    _atom->setPositionAtStep(TimeStepToSet, currentPosition + PositionUpdate);
+  }
+
   /** 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().
+   * time steps present, i.e. this must be preceded by a simple anneal().
    *
    * We assume that forces have just been calculated.
@@ -274 +311 @@
     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);
+
+      annealAtom_BarzilaiBorwein(*iter, OldTimeStep, CurrentTimeStep, _TimeStep);
+
+      // extract largest components for showing progress of annealing
       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 = getBarzilaiBorweinStepwidth(PositionDifference, GradientDifference);
-      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]));
-
-      // finally set new values
-      (*iter)->setPositionAtStep(_TimeStep, currentPosition + PositionUpdate);
     }
 
@@ -399 +419 @@
 
     std::map<atomId_t, Vector> GatheredUpdates; //!< gathers all updates which are applied at the end
+    std::map<atomId_t, double> LargestUpdate_per_Atom; //!< check whether updates cancelled each other
     for(typename AtomSetMixin<T>::iterator iter = AtomicForceManipulator<T>::atoms.begin();
         iter != AtomicForceManipulator<T>::atoms.end(); ++iter) {
@@ -428 +449 @@
       }
       Vector PositionUpdate = stepwidth * currentGradient;
+      // cap updates (if non-zero) at 0.2 angstroem. BB tends to overshoot.
+      for (size_t i=0;i<NDIM;++i)
+        if (fabs(PositionUpdate[i]) > MYEPSILON)
+          PositionUpdate[i] = std::min(0.5, fabs(PositionUpdate[i]))*PositionUpdate[i]/fabs(PositionUpdate[i]);
       LOG(3, "DEBUG: Update would be " << stepwidth << "*" << currentGradient << " = " << PositionUpdate);
 
@@ -558 +583 @@
         // hence, the decision is only based on sides once we have picked a side
         // depending on all bonds associated with have same good type.
-        const double sign = sides[first_per_side[sideno]];
+
+        // away from gradient (minus) and contracting
+        // or towards gradient (plus) and expanding
+        // gather all on same side and remove
+        const double sign =
+            (sides[first_per_side[sideno]] == types[first_per_side[sideno]])
+            ? sides[first_per_side[sideno]] : -1.*sides[first_per_side[sideno]];
+
         LOG(4, "DEBUG: Removing edges from side with sign " << sign);
         BondList::const_iterator bonditer = ListOfBonds.begin();
@@ -603 +635 @@
           if (GatheredUpdates.count((*setiter))) {
             GatheredUpdates[(*setiter)] += factor*PositionUpdate;
+            LargestUpdate_per_Atom[(*setiter)] =
+                std::max(LargestUpdate_per_Atom[(*setiter)], factor*PositionUpdate.Norm());
           } else {
             GatheredUpdates.insert(
@@ -608 +642 @@
                     (*setiter),
                     factor*PositionUpdate) );
+            LargestUpdate_per_Atom.insert(
+                std::make_pair(
+                    (*setiter),
+                    factor*PositionUpdate.Norm()) );
           }
         }
@@ -617 +655 @@
                 walker.getId(),
                 PositionUpdate) );
+        LargestUpdate_per_Atom.insert(
+            std::make_pair(
+                walker.getId(),
+                PositionUpdate.Norm()) );
       }
     }
@@ -624 +666 @@
       atom &walker = *(*iter);
       // extract largest components for showing progress of annealing
-      const Vector currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep);
+      const Vector &currentGradient = walker.getAtomicForceAtStep(CurrentTimeStep);
       for(size_t i=0;i<NDIM;++i)
         maxComponents[i] = std::max(maxComponents[i], fabs(currentGradient[i]));
@@ -653 +695 @@
       for (size_t i=0;i<NDIM;++i)
         LargestUpdate[i] = std::max(LargestUpdate[i], fabs(update[i]));
-      walker->setPositionAtStep(_TimeStep,
-          walker->getPositionAtStep(CurrentTimeStep) + update); // - CommonTranslation);
+
+      std::map<atomId_t, double>::const_iterator largestiter = LargestUpdate_per_Atom.find(atomid);
+      ASSERT( largestiter != LargestUpdate_per_Atom.end(),
+          "ForceAnnealing() - walker with id "+toString(atomid)+" not in LargestUpdates.");
+      // if we had large updates but their sum is very small
+      if (update.Norm()/largestiter->second > MYEPSILON) {
+        walker->setPositionAtStep(_TimeStep,
+            walker->getPositionAtStep(CurrentTimeStep) + update); // - CommonTranslation);
+      } else {
+        // then recalc update with simple anneal
+        LOG(2, "WARNING: Updates on atom " << *iter << " cancel themselves, performing simple anneal step.");
+        annealAtom_BarzilaiBorwein(walker, OldTimeStep, CurrentTimeStep, _TimeStep);
+      }
     }
     LOG(1, "STATUS: Largest absolute update components are " << LargestUpdate);