Changeset f433ec for src/Dynamics/BondVectors.cpp

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.

File:

• src/Dynamics/BondVectors.cpp (modified) (4 diffs)

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;
+}