Context Navigation

-              rac868ee
+              r52e7da
+}
-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)
+{
 …
     // 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());
+      }
+    }
+    TrainingData AngleData(
+        boost::bind(&Extractors::reorderArgumentsByIncreasingDistance,
+            boost::bind(&Extractors::gatherAllSymmetricDistanceArguments,
+                boost::bind(&Extractors::gatherPositionOfTuples,
+                    _1, Fragment::charges_t(3,6.)
+                ), _2 // gather carbon triples
+            )
+          )
+        );
+    AngleData(homologies.getHomologousGraphs(graph));
+    LOG(1, "INFO: I gathered the following training data: " << AngleData);
     // NOTICE that distance are in bohrradi as they come from MPQC!
 …
     angle.setParameters(params);
     FunctionModel &model = angle;
+    FunctionApproximation approximator(
+        DistanceEnergyVector.first.begin()->size(),
+        DistanceEnergyVector.second.begin()->size(),
+        model);
+    approximator.setTrainingData(DistanceEnergyVector.first,DistanceEnergyVector.second);
+    FunctionApproximation approximator(AngleData, model);
     if (model.isBoxConstraint() && approximator.checkParameterDerivatives())
       approximator(FunctionApproximation::ParameterDerivative);
 …
     TrainingData MorseData(
         boost::bind(&Extractors::gatherAllSymmetricDistanceArguments,
             boost::bind(&Extractors::gatherDistanceOfTuples,
+            boost::bind(&Extractors::gatherPositionOfTuples,
                 _1, Fragment::charges_t(2,6.)
             ), _2 // gather first carbon pair

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

Context Navigation

Changeset 52e7da for src/LevMartester.cpp

Legend:

src/LevMartester.cpp

Download in other formats: