Context Navigation

-              r94f567
+              r9340ee
 #include "CodePatterns/MemDebug.hpp"
+#include <boost/assign.hpp>
 #include <boost/filesystem.hpp>
 #include <boost/program_options.hpp>
 …
 #include "Helpers/defs.hpp"
 #include "Potentials/Specifics/PairPotential_Morse.hpp"
+#include "Potentials/Specifics/PairPotential_Angle.hpp"
 #include "Potentials/Specifics/SaturationPotential.hpp"
 namespace po = boost::program_options;
+using namespace boost::assign;
+HomologyGraph getFirstGraphWithThreeCarbons(const HomologyContainer &homologies)
+{
+  FragmentNode SaturatedCarbon(6,4); // carbon has atomic number 6 and should have 4 bonds for C3H8
+  FragmentNode DanglingCarbon(6,3); // carbon has atomic number 6 and should have 3 pure bonds for C3H8
+  for (HomologyContainer::container_t::const_iterator iter =
+      homologies.begin(); iter != homologies.end(); ++iter) {
+    if ((iter->first.hasNode(SaturatedCarbon,2)) && (iter->first.hasNode(DanglingCarbon,1)))
+      return iter->first;
+  }
+  return HomologyGraph();
+}
 HomologyGraph getFirstGraphWithTwoCarbons(const HomologyContainer &homologies)
 …
+}
+double
+function_angle(
+    const double &r_ij,
+    const double &r_ik,
+    const double &r_jk
+  )
+{
+//  Info info(__func__);
+  const double angle = pow(r_ij,2.) + pow(r_ik,2.) - pow(r_jk,2.);
+  const double divisor = 2.* r_ij * r_ik;
+//  LOG(2, "DEBUG: cos(theta)= " << angle/divisor);
+  if (divisor == 0.)
+    return 0.;
+  else
+    return angle/divisor;
+}
 int main(int argc, char **argv)
 …
     LOG(1, "INFO: graph " << iter->first << " has Fragment "
         << iter->second.first << " and associated energy " << iter->second.second << ".");
+  }
+  /******************** Angle TRAINING ********************/
+  {
+    // then we ought to pick the right HomologyGraph ...
+    const HomologyGraph graph = getFirstGraphWithThreeCarbons(homologies);
+    LOG(1, "First representative graph containing three saturated carbons is " << graph << ".");
+    // Afterwards we go through all of this type and gather the distance and the energy value
+    typedef std::pair<
+        FunctionApproximation::inputs_t,
+        FunctionApproximation::outputs_t> InputOutputVector_t;
+    InputOutputVector_t DistanceEnergyVector;
+    std::pair<HomologyContainer::const_iterator, HomologyContainer::const_iterator> range =
+        homologies.getHomologousGraphs(graph);
+    for (HomologyContainer::const_iterator fragiter = range.first; fragiter != range.second; ++fragiter) {
+      // get distance out of Fragment
+      const double &energy = fragiter->second.second;
+      const Fragment &fragment = fragiter->second.first;
+      const Fragment::charges_t charges = fragment.getCharges();
+      const Fragment::positions_t positions = fragment.getPositions();
+      std::vector< std::pair<Vector, size_t> > DistanceVectors;
+      for (Fragment::charges_t::const_iterator chargeiter = charges.begin();
+          chargeiter != charges.end(); ++chargeiter) {
+        if (*chargeiter == 6) {
+          Fragment::positions_t::const_iterator positer = positions.begin();
+          const size_t steps = std::distance(charges.begin(), chargeiter);
+          std::advance(positer, steps);
+          DistanceVectors.push_back(
+              std::make_pair(Vector((*positer)[0], (*positer)[1], (*positer)[2]),
+                  steps));
+        }
+      }
+      if (DistanceVectors.size() == (size_t)3) {
+        FunctionModel::arguments_t args(3);
+        // we require specific ordering of the carbons: ij, ik, jk
+        typedef std::vector< std::pair<size_t, size_t> > indices_t;
+        indices_t indices;
+        indices += std::make_pair(0,1), std::make_pair(0,2), std::make_pair(1,2);
+        // create the three arguments
+        for (indices_t::const_iterator iter = indices.begin(); iter != indices.end(); ++iter) {
+          const size_t &firstindex = iter->first;
+          const size_t &secondindex = iter->second;
+          argument_t &arg = args[(size_t)std::distance(const_cast<const indices_t&>(indices).begin(), iter)];
+          arg.indices.first = DistanceVectors[firstindex].second;
+          arg.indices.second = DistanceVectors[secondindex].second;
+          arg.distance = DistanceVectors[firstindex].first.distance(DistanceVectors[secondindex].first);
+          arg.globalid = DistanceEnergyVector.first.size();
+        }
+        // make largest distance last to create correct angle
+        // (this would normally depend on the order of the nodes in the subgraph)
+        std::list<argument_t> sorted_args;
+        double greatestdistance = 0.;
+        for(FunctionModel::arguments_t::const_iterator iter = args.begin(); iter != args.end(); ++iter)
+          greatestdistance = std::max(greatestdistance, iter->distance);
+        for(FunctionModel::arguments_t::const_iterator iter = args.begin(); iter != args.end(); ++iter)
+          if (iter->distance == greatestdistance)
+            sorted_args.push_back(*iter);
+          else
+            sorted_args.push_front(*iter);
+        // and add the training pair
+        DistanceEnergyVector.first.push_back( FunctionModel::arguments_t(sorted_args.begin(), sorted_args.end()) );
+        DistanceEnergyVector.second.push_back( FunctionModel::results_t(1,energy) );
+      } else {
+        ELOG(2, "main() - found not exactly three carbon atoms in fragment "
+            << fragment << ".");
+      }
+    }
+    // print training data for debugging
+    {
+      LOG(1, "INFO: I gathered the following (" << DistanceEnergyVector.first.size()
+          << "," << DistanceEnergyVector.second.size() << ") data pairs: ");
+      FunctionApproximation::inputs_t::const_iterator initer = DistanceEnergyVector.first.begin();
+      FunctionApproximation::outputs_t::const_iterator outiter = DistanceEnergyVector.second.begin();
+      for (; initer != DistanceEnergyVector.first.end(); ++initer, ++outiter) {
+        std::stringstream stream;
+        const double cos_angle = function_angle((*initer)[0].distance,(*initer)[1].distance,(*initer)[2].distance);
+        for (size_t index = 0; index < (*initer).size(); ++index)
+           stream << " (" << (*initer)[index].indices.first << "," << (*initer)[index].indices.second
+              << ") " << (*initer)[index].distance;
+        stream << " with energy " << *outiter << " and cos(angle) " << cos_angle;
+        LOG(1, "INFO:" << stream.str());
+      }
+    }
+    // NOTICE that distance are in bohrradi as they come from MPQC!
+    // now perform the function approximation by optimizing the model function
+    FunctionModel::parameters_t params(PairPotential_Angle::MAXPARAMS, 0.);
+    params[PairPotential_Angle::energy_offset] =  -1.;
+    params[PairPotential_Angle::spring_constant] =  1.;
+    params[PairPotential_Angle::equilibrium_distance] =  0.2;
+    PairPotential_Angle angle;
+    LOG(0, "INFO: Initial parameters are " << params << ".");
+    angle.setParameters(params);
+    FunctionModel &model = angle;
+    FunctionApproximation approximator(
+        DistanceEnergyVector.first.begin()->size(),
+        DistanceEnergyVector.second.begin()->size(),
+        model);
+    approximator.setTrainingData(DistanceEnergyVector.first,DistanceEnergyVector.second);
+    if (model.isBoxConstraint() && approximator.checkParameterDerivatives())
+      approximator(FunctionApproximation::ParameterDerivative);
+    else
+      ELOG(0, "We require parameter derivatives for a box constraint minimization.");
+    params = model.getParameters();
+    LOG(0, "RESULT: Best parameters are " << params << ".");
+  }

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 9340ee for src/LevMartester.cpp

Legend:

src/LevMartester.cpp

Download in other formats: