Changeset 90e540 for src/Dynamics/ForceAnnealing.hpp

Timestamp:

Jul 20, 2017, 9:38:39 AM (8 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

ForceAnnealing_with_BondGraph_continued

Children:

ccbdf4

Parents:

6d7c73

git-author:

Frederik Heber <frederik.heber@…> (07/18/17 22:24:12)

git-committer:

Frederik Heber <frederik.heber@…> (07/20/17 09:38:39)

Message:

We now obtain weights via levmar minimization.

• this should yield the best possible weights within the interval of [1/n,1.].
• note that we cannot always get an exact solution because of this constraint.

File:

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

src/Dynamics/ForceAnnealing.hpp

@@ -117 +117 @@
       LOG(2, "DEBUG: current step is #" << currentStep);
 
+      // bond graph annealing is always followed by a normal annealing
       if (_UseBondgraph)
         annealWithBondGraph(_CurrentTimeStep, _offset, maxComponents);
-      else
-        anneal(_CurrentTimeStep, _offset, maxComponents);
+      anneal(_CurrentTimeStep, _offset, maxComponents);
     }
 
@@ -203 +203 @@
   }
 
+  // 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.
    *
@@ -260 +297 @@
     const BondVectors::container_t &sorted_bonds = bv.getSorted();
 
-    // 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
-    };
     std::vector< // time step
       std::vector< // which bond side
@@ -280 +310 @@
     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) {
       // TODO: We have this extra step in between because of DoOutput copying
@@ -304 +336 @@
               weights_per_atom[timestep].insert(
                   std::make_pair(walker.getId(),
-                  bv.getWeightsForAtomAtStep(walker, CurrentStep)) );
+                  bv.getWeightsForAtomAtStep(walker, BondVectors, CurrentStep)) );
           ASSERT( inserter.second,
               "ForceAnnealing::operator() - weight map for atom "+toString(walker)
@@ -316 +348 @@
           // projected gradient over all bonds and place in one of projected_forces
           // using the obtained weights
-          {
-            BondVectors::weights_t::const_iterator weightiter = weights.begin();
-            std::vector<Vector>::const_iterator vectoriter = BondVectors.begin();
-            Vector forcesum_check;
-            for(BondList::const_iterator bonditer = ListOfBonds.begin();
-                bonditer != ListOfBonds.end(); ++bonditer, ++weightiter, ++vectoriter) {
-              const bond::ptr &current_bond = *bonditer;
-              const Vector &BondVector = *vectoriter;
-
-              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] = (*weightiter)*walkerGradient.ScalarProduct(BondVector);
-              LOG(4, "DEBUG: BondVector " << BondVector << " receives projected force of "
-                  << forcelist[index]);
-              forcesum_check += forcelist[index] * BondVector;
-            }
-            ASSERT( weightiter == weights.end(),
-                "ForceAnnealing - weightiter is not at end when it should be.");
-            ASSERT( vectoriter == BondVectors.end(),
-                "ForceAnnealing - vectoriter is not at end when it should be.");
-            LOG(3, "DEBUG: sum of projected forces is " << forcesum_check);
-          }
-
+          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 "
@@ -465 +477 @@
     }
 
-    // apply the gathered updates
+    // 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) {
@@ -476 +488 @@
           << ", namely " << *walker);
       walker->setPosition( walker->getPosition() + update );
+      walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
     }
   }