Changeset f433ec for src/Dynamics

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:

07d4b1

Parents:

7b4e67

git-author:

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

git-committer:

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

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.

Location:

src/Dynamics

Files:

• BondVectors.cpp (modified) (4 diffs)
• BondVectors.hpp (modified) (4 diffs)
• ForceAnnealing.hpp (modified) (8 diffs)
• unittests/BondVectorsUnitTest.cpp (modified) (2 diffs)
• unittests/Makefile.am (modified) (1 diff)

src/Dynamics/BondVectors.cpp

@@ -45 +45 @@
 #include "CodePatterns/Assert.hpp"
 #include "CodePatterns/Log.hpp"
+
+#include <levmar.h>
 
 #include "Atom/atom.hpp"
@@ -111 +113 @@
 }
 
+struct WeightsMinimization {
+
+  static void evaluate(double *p, double *x, int m, int n, void *data)
+  {
+    // current weights in p, output to x
+    double *matrix = static_cast<double*>(data);
+    for(size_t i=0;i<(size_t)n;++i) {
+      x[i] = 0.;
+      for(size_t j=0;j<(size_t)n;++j)
+        x[i] += p[j]*matrix[i*n+j];
+//      x[i] = .5*x[i]*x[i];
+    }
+  }
+
+  static void evaluate_derivative(double *p, double *jac, int m, int n, void *data)
+  {
+    // current weights in p, output to x
+    double *matrix = static_cast<double*>(data);
+//    // tmp = (Bx - 1)
+//    double *tmp = new double[n];
+//    for(size_t i=0;i<(size_t)n;++i) {
+//      tmp[i] = -1.;
+//      for(size_t j=0;j<(size_t)n;++j)
+//        tmp[i] += p[j]*matrix[i*n+j];
+//    }
+//    // tmp(i) * B_(ij)
+//    for(size_t i=0;i<(size_t)n;++i)
+//      for(size_t j=0;j<(size_t)n;++j)
+//        jac[i*n+j] = tmp[i]*matrix[i*n+j];
+//    delete[] tmp ;
+    for(size_t i=0;i<(size_t)n;++i)
+      for(size_t j=0;j<(size_t)n;++j)
+        jac[i*n+j] = matrix[i*n+j];
+  }
+
+  static double* getMatrixFromBondVectors(
+      const std::vector<Vector> &_bondvectors
+      )
+  {
+    const size_t n = _bondvectors.size();
+    double *matrix = new double[n*n];
+    size_t i=0;
+    for (std::vector<Vector>::const_iterator iter = _bondvectors.begin();
+        iter != _bondvectors.end(); ++iter, ++i) {
+      size_t j=0;
+      for (std::vector<Vector>::const_iterator otheriter = _bondvectors.begin();
+          otheriter != _bondvectors.end(); ++otheriter, ++j) {
+        // only magnitude is important not the sign
+        matrix[i*n+j] = fabs((*iter).ScalarProduct(*otheriter));
+      }
+    }
+
+    return matrix;
+  }
+};
+
+
+BondVectors::weights_t BondVectors::getWeightsForAtomAtStep(
+    const atom &_walker,
+    const std::vector<Vector> &_bondvectors,
+    const size_t &_step) const
+{
+  // let levmar optimize
+  register int i, j;
+  int ret;
+  double *p;
+  double *x;
+  int n=_bondvectors.size();
+  double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
+
+  double *matrix = WeightsMinimization::getMatrixFromBondVectors(_bondvectors);
+
+  weights_t weights(n, 0.);
+
+  // minim. options [\tau, \epsilon1, \epsilon2, \epsilon3]. Respectively the scale factor for initial \mu,
+  // * stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2.
+  opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
+  opts[4]= LM_DIFF_DELTA; // relevant only if the Jacobian is approximated using finite differences; specifies forward differencing
+  //opts[4]=-LM_DIFF_DELTA; // specifies central differencing to approximate Jacobian; more accurate but more expensive to compute!
+
+  // prepare initial values for weights
+  p = new double[n];
+  for (i=0;i<n;++i)
+    p[i] = 1.;
+  // prepare output value, i.e. the row sums
+  x = new double[n];
+  for (i=0;i<n;++i)
+    x[i] = 1.;
+
+  {
+    double *work, *covar;
+    work=(double *)malloc((LM_DIF_WORKSZ(n, n)+n*n)*sizeof(double));
+    if(!work){
+      ELOG(0, "BondVectors::getWeightsForAtomAtStep() - memory allocation request failed.");
+      return weights;
+    }
+    covar=work+LM_DIF_WORKSZ(n, n);
+
+    // give this pointer as additional data to construct function pointer in
+    // LevMarCallback and call
+    double *lb = new double[n];
+    double *ub = new double[n];
+    for (i=0;i<n;++i) {
+      lb[i] = 1./(double)n;
+      ub[i] = 1.;
+    }
+    ret=dlevmar_bc_der(
+        &WeightsMinimization::evaluate,
+        &WeightsMinimization::evaluate_derivative,
+        p, x, n, n, lb, ub, NULL, 1000, opts, info, work, covar, matrix); // no Jacobian, caller allocates work memory, covariance estimated
+    delete[] lb;
+    delete[] ub;
+
+    if (0)
+    {
+      std::stringstream covar_msg;
+      covar_msg << "Covariance of the fit:\n";
+      for(i=0; i<n; ++i){
+        for(j=0; j<n; ++j)
+          covar_msg << covar[i*n+j] << " ";
+        covar_msg << std::endl;
+      }
+      covar_msg << std::endl;
+      LOG(1, "INFO: " << covar_msg.str());
+    }
+
+    free(work);
+  }
+
+  if (DoLog(4)) {
+   std::stringstream result_msg;
+   result_msg << "Levenberg-Marquardt returned " << ret << " in " << info[5] << " iter, reason " << info[6] << "\nSolution: ";
+   for(i=0; i<n; ++i)
+     result_msg << p[i] << " ";
+   result_msg << "\n\nMinimization info:\n";
+   std::vector<std::string> infonames(LM_INFO_SZ);
+   infonames[0] = std::string("||e||_2 at initial p");
+   infonames[1] = std::string("||e||_2");
+   infonames[2] = std::string("||J^T e||_inf");
+   infonames[3] = std::string("||Dp||_2");
+   infonames[4] = std::string("mu/max[J^T J]_ii");
+   infonames[5] = std::string("# iterations");
+   infonames[6] = std::string("reason for termination");
+   infonames[7] = std::string(" # function evaluations");
+   infonames[8] = std::string(" # Jacobian evaluations");
+   infonames[9] = std::string(" # linear systems solved");
+   for(i=0; i<LM_INFO_SZ; ++i)
+      result_msg << infonames[i] << ": " << info[i] << " ";
+    result_msg << std::endl;
+    LOG(4, "DEBUG: " << result_msg.str());
+  }
+
+  std::copy(p, p+n, weights.begin());
+  LOG(4, "DEBUG: Weights for atom #" << _walker.getId() << ": " << weights);
+
+  delete[] p;
+  delete[] x;
+  delete[] matrix;
+
+  return weights;
+}
+
 BondVectors::weights_t BondVectors::getWeightsForAtomAtStep(
     const atom &_walker,
@@ -118 +282 @@
       getAtomsBondVectorsAtStep(_walker, _step);
 
-  weights_t weights;
-  for (std::vector<Vector>::const_iterator iter = BondVectors.begin();
-      iter != BondVectors.end(); ++iter) {
+  return getWeightsForAtomAtStep(_walker, BondVectors, _step);
+}
+
+Vector BondVectors::getRemnantGradientForAtomAtStep(
+    const atom &_walker,
+    const std::vector<Vector> _BondVectors,
+    const BondVectors::weights_t &_weights,
+    const size_t &_step,
+    forcestore_t _forcestore) const
+{
+  const Vector &walkerGradient = _walker.getAtomicForceAtStep(_step);
+  BondVectors::weights_t::const_iterator weightiter = _weights.begin();
+  std::vector<Vector>::const_iterator vectoriter = _BondVectors.begin();
+  const BondList& ListOfBonds = _walker.getListOfBonds();
+
+  Vector forcesum;
+  for(BondList::const_iterator bonditer = ListOfBonds.begin();
+      bonditer != ListOfBonds.end(); ++bonditer, ++weightiter, ++vectoriter) {
+    const bond::ptr &current_bond = *bonditer;
+    const Vector &BondVector = *vectoriter;
+
+    const double temp = (*weightiter)*walkerGradient.ScalarProduct(BondVector);
+    _forcestore(_walker, current_bond, _step, temp);
+    LOG(4, "DEBUG: BondVector " << BondVector << " receives projected force of "
+        << temp);
+    forcesum += temp * BondVector;
+  }
+  ASSERT( weightiter == _weights.end(),
+      "BondVectors::getRemnantGradientForAtomAtStep() - weightiter is not at end when it should be.");
+  ASSERT( vectoriter == _BondVectors.end(),
+      "BondVectors::getRemnantGradientForAtomAtStep() - vectoriter is not at end when it should be.");
+
+  return walkerGradient-forcesum;
+}
+
+bool BondVectors::getCheckWeightSumForAtomAtStep(
+    const atom &_walker,
+    const std::vector<Vector> _BondVectors,
+    const BondVectors::weights_t &_weights,
+    const size_t &_step) const
+{
+  bool status = true;
+  for (std::vector<Vector>::const_iterator iter = _BondVectors.begin();
+      iter != _BondVectors.end(); ++iter) {
     std::vector<double> scps;
-    scps.reserve(BondVectors.size());
+    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(),
+        _BondVectors.begin(), _BondVectors.end(),
+        _weights.begin(),
         std::back_inserter(scps),
         boost::bind(static_cast< double (*)(double) >(&fabs),
@@ -153 +338 @@
         );
     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;
-}
+    if (fabs(scp_sum -1.) > MYEPSILON)
+      status = false;
+  }
+
+  return status;
+}

src/Dynamics/BondVectors.hpp

@@ -17 +17 @@
 #include <map>
 #include <vector>
+
+#include <boost/function.hpp>
 
 #include "CodePatterns/Assert.hpp"
@@ -93 +95 @@
    *
    * \param _walker atom to get BondVectors for
+   * \param _bondvectors precalculated bond vectors for given \a _walker
+   * \param _step time step for which the bond vector is request
+   */
+  weights_t getWeightsForAtomAtStep(
+      const atom &_walker,
+      const std::vector<Vector> &_bondvectors,
+      const size_t &_step) const;
+
+  /** 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
    */
@@ -98 +115 @@
       const atom &_walker,
       const size_t &_step) const;
+
+  /** Function typedef to store the bond gradient into a specific container
+   * depending on the atom, its current bond and the time step.
+   */
+  typedef boost::function<void (
+      const atom &,
+      const bond::ptr &,
+      const size_t &,
+      const double)> forcestore_t;
+
+  /** Function calculates the remaining part of the atomic gradient that is
+   * not captured by the sum of the force along the Bond Vectors.
+   *
+   * \param _walker atom to get BondVectors for
+   * \param _BondVectors precalculated bond vectors for given \a _walker
+   * \param _weights weight per bond vector (as it is a frame, not a basis)
+   * \param _step time step for which the bond vector is request
+   * \param _forcestore additional function which may be used to store each
+   *        calculated bond force in a bound container
+   */
+  Vector getRemnantGradientForAtomAtStep(
+      const atom &_walker,
+      const std::vector<Vector> _BondVectors,
+      const BondVectors::weights_t &_weights,
+      const size_t &_step,
+      forcestore_t _forcestore) const;
 
 private:
@@ -105 +148 @@
    */
   void recalculateBondVectorsAtStep(const size_t &_step) const;
+
+  /** Helper function to check whether weights sum up to one for each
+   * Bond Vector.
+   *
+   * \param _walker atom to get BondVectors for
+   * \param _BondVectors precalculated bond vectors for given \a _walker
+   * \param _weights weight per bond vector (as it is a frame, not a basis)
+   * \param _step time step for which the bond vector is request
+   */
+  bool getCheckWeightSumForAtomAtStep(
+      const atom &_walker,
+      const std::vector<Vector> _BondVectors,
+      const BondVectors::weights_t &_weights,
+      const size_t &_step) const;
 
 private:

src/Dynamics/ForceAnnealing.hpp

@@ -119 +119 @@
       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);
     }
 
@@ -217 +217 @@
   }
 
+  // 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.
    *
@@ -274 +311 @@
     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
@@ -294 +324 @@
     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) {
       const size_t CurrentStep = CurrentTimeStep-timestep-1 >= 0 ? CurrentTimeStep-timestep-1 : 0;
@@ -320 +352 @@
               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)
@@ -332 +364 @@
           // 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 "
@@ -487 +499 @@
     }
 
-    // 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) {
@@ -501 +513 @@
           + update);
       walker->setAtomicVelocity(update);
-//      walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
+      walker->setAtomicForce( RemnantGradient_per_atom[walker->getId()] );
     }
   }

src/Dynamics/unittests/BondVectorsUnitTest.cpp

@@ -197 +197 @@
   CPPUNIT_ASSERT( fabs(weight_sum - 1.) < MYEPSILON );
   // check weight
-  CPPUNIT_ASSERT_EQUAL( weights[0], 1. );
+  CPPUNIT_ASSERT( fabs(weight_sum - 1.) < 1e-10 );
 }
 
@@ -289 +289 @@
   // check number of weights
   CPPUNIT_ASSERT_EQUAL( weights.size(), (size_t)5 );
-  // check sum of weights: one linear independent, two dependent vectors = 1 + 2*0.5
-  const double weight_sum = std::accumulate(weights.begin(), weights.end(), 0.);
-  CPPUNIT_ASSERT( fabs(weight_sum - 2.) < 1e-10 );
-}
+  // check sum of weights
+  const double weight_sum = std::accumulate(weights.begin(), weights.end(), 0.);
+  CPPUNIT_ASSERT( fabs(weights[0] - .372244) < 1e-6 );
+  CPPUNIT_ASSERT( fabs(weights[1] - .529694) < 1e-6 );
+  CPPUNIT_ASSERT( fabs(weights[2] - .2) < 1e-6 );
+  CPPUNIT_ASSERT( fabs(weights[3] - .248464) < 1e-6 );
+  CPPUNIT_ASSERT( fabs(weights[4] - .248464) < 1e-6 );
+}

src/Dynamics/unittests/Makefile.am

@@ -11 +11 @@
   BondVectorsUnitTest
 
-XFAIL_TESTS += BondVectorsUnitTest
-  
 TESTS += $(DYNAMICSTESTS)
 check_PROGRAMS += $(DYNAMICSTESTS)