Changeset 825d33

Timestamp:

Apr 10, 2018, 6:43:30 AM (7 years ago)

Author:

Frederik Heber <frederik.heber@…>

Branches:

AutomationFragmentation_failures, Candidate_v1.6.1, ChemicalSpaceEvaluator, Exclude_Hydrogens_annealWithBondGraph, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_contraction-expansion, Gui_displays_atomic_force_velocity, PythonUI_with_named_parameters, StoppableMakroAction, TremoloParser_IncreasedPrecision

Children:

c1c6cb

Parents:

9861d0

git-author:

Frederik Heber <frederik.heber@…> (06/29/17 14:53:38)

git-committer:

Frederik Heber <frederik.heber@…> (04/10/18 06:43:30)

Message:

Extracted calculation of weights per atom into BondVectors.

Location:

src/Dynamics

Files:

• BondVectors.cpp (modified) (3 diffs)
• BondVectors.hpp (modified) (1 diff)
• ForceAnnealing.hpp (modified) (2 diffs)

src/Dynamics/BondVectors.cpp

@@ -39 +39 @@
 
 #include <algorithm>
+#include <functional>
 #include <iterator>
+#include <numeric>
 
 #include "CodePatterns/Assert.hpp"
@@ -46 +48 @@
 #include "Atom/atom.hpp"
 #include "Bond/bond.hpp"
+#include "Helpers/defs.hpp"
 
 void BondVectors::recalculateBondVectorsAtStep(
@@ -107 +110 @@
   return BondVectors;
 }
+
+BondVectors::weights_t BondVectors::getWeightsForAtomAtStep(
+    const atom &_walker,
+    const size_t &_step) const
+{
+  const std::vector<Vector> BondVectors =
+      getAtomsBondVectorsAtStep(_walker, _step);
+
+  weights_t weights;
+  for (std::vector<Vector>::const_iterator iter = BondVectors.begin();
+      iter != BondVectors.end(); ++iter) {
+    std::vector<double> scps;
+    scps.reserve(BondVectors.size());
+    std::transform(
+        BondVectors.begin(), BondVectors.end(),
+        std::back_inserter(scps),
+        boost::bind(static_cast< double (*)(double) >(&fabs),
+            boost::bind(&Vector::ScalarProduct, boost::cref(*iter), _1))
+        );
+    const double scp_sum = std::accumulate(scps.begin(), scps.end(), 0.);
+    ASSERT( (scp_sum-1.) > -MYEPSILON,
+        "ForceAnnealing() - sum of weights must be equal or larger one but is "
+        +toString(scp_sum));
+    weights.push_back( 1./scp_sum );
+  }
+  LOG(4, "DEBUG: Weights for atom #" << _walker.getId() << ": " << weights);
+
+  // 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) {
+    std::vector<double> scps;
+    scps.reserve(BondVectors.size());
+    std::transform(
+        BondVectors.begin(), BondVectors.end(),
+        weights.begin(),
+        std::back_inserter(scps),
+        boost::bind(static_cast< double (*)(double) >(&fabs),
+            boost::bind(std::multiplies<double>(),
+                boost::bind(&Vector::ScalarProduct, boost::cref(*iter), _1),
+                _2))
+        );
+    const double scp_sum = std::accumulate(scps.begin(), scps.end(), 0.);
+    ASSERT( fabs(scp_sum - 1.) < MYEPSILON,
+        "ForceAnnealing::operator() - for BondVector "+toString(*iter)
+        +" we have weighted scalar product of "+toString(scp_sum)+" != 1.");
+  }
+#endif
+  return weights;
+}

src/Dynamics/BondVectors.hpp

@@ -83 +83 @@
       const size_t &_step) const;
 
+  //!> typedef for the weights for the Bondvectors of a single atom
+  typedef std::deque<double> weights_t;
+
+  /** Calculates the weights for a frame where each Bondvector of the
+   * given atom is a vector of the frame.
+   *
+   * The idea is that we can represent any vector by appropriate weights such
+   * that is still sums up to one.
+   *
+   * \param _walker atom to get BondVectors for
+   * \param _step time step for which the bond vector is request
+   */
+  weights_t getWeightsForAtomAtStep(
+      const atom &_walker,
+      const size_t &_step) const;
+
 private:
   /** Calculates the bond vector for each bond in the internal container.

src/Dynamics/ForceAnnealing.hpp

@@ -292 +292 @@
 
     // for each atom we need to gather weights and then project the gradient
-    typedef std::deque<double> weights_t;
-    typedef std::map<atomId_t, weights_t > weights_per_atom_t;
+    typedef std::map<atomId_t, BondVectors::weights_t > weights_per_atom_t;
     std::vector<weights_per_atom_t> weights_per_atom(2);
     for (size_t timestep = 0; timestep <= 1; ++timestep) {
@@ -313 +312 @@
 
           // gather subset of BondVectors for the current atom
-          std::vector<Vector> BondVectors = bv.getAtomsBondVectorsAtStep(walker, CurrentStep);
+          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
-          std::pair<weights_per_atom_t::iterator, bool> inserter =
+          const std::pair<weights_per_atom_t::iterator, bool> inserter =
               weights_per_atom[timestep].insert(
-                  std::make_pair(walker.getId(), weights_t()) );
+                  std::make_pair(walker.getId(),
+                  bv.getWeightsForAtomAtStep(walker, CurrentStep)) );
           ASSERT( inserter.second,
               "ForceAnnealing::operator() - weight map for atom "+toString(walker)
               +" and time step "+toString(timestep)+" 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;
-            scps.reserve(BondVectors.size());
-            std::transform(
-                BondVectors.begin(), BondVectors.end(),
-                std::back_inserter(scps),
-                boost::bind(static_cast< double (*)(double) >(&fabs),
-                    boost::bind(&Vector::ScalarProduct, boost::cref(*iter), _1))
-                );
-            const double scp_sum = std::accumulate(scps.begin(), scps.end(), 0.);
-            ASSERT( (scp_sum-1.) > -MYEPSILON,
-                "ForceAnnealing() - sum of weights must be equal or larger one but is "
-                +toString(scp_sum));
-            weights.push_back( 1./scp_sum );
-          }
-          LOG(4, "DEBUG: Weights for atom #" << walker.getId() << ": " << weights);
-
-          // 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) {
-            std::vector<double> scps;
-            scps.reserve(BondVectors.size());
-            std::transform(
-                BondVectors.begin(), BondVectors.end(),
-                weights.begin(),
-                std::back_inserter(scps),
-                boost::bind(static_cast< double (*)(double) >(&fabs),
-                    boost::bind(std::multiplies<double>(),
-                        boost::bind(&Vector::ScalarProduct, boost::cref(*iter), _1),
-                        _2))
-                );
-            const double scp_sum = std::accumulate(scps.begin(), scps.end(), 0.);
-            ASSERT( fabs(scp_sum - 1.) < MYEPSILON,
-                "ForceAnnealing::operator() - for BondVector "+toString(*iter)
-                +" we have weighted scalar product of "+toString(scp_sum)+" != 1.");
-          }
-#endif
+          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
           {
-            weights_t::const_iterator weightiter = weights.begin();
+            BondVectors::weights_t::const_iterator weightiter = weights.begin();
             std::vector<Vector>::const_iterator vectoriter = BondVectors.begin();
             Vector forcesum_check;