Changes in src/Actions/PotentialAction/FitPotentialAction.cpp [707a2b:0ea063]

File:

• src/Actions/PotentialAction/FitPotentialAction.cpp (modified) (4 diffs)

src/Actions/PotentialAction/FitPotentialAction.cpp

@@ -55 +55 @@
 #include "Fragmentation/Homology/HomologyGraph.hpp"
 #include "Fragmentation/Summation/SetValues/Fragment.hpp"
-#include "FunctionApproximation/Extractors.hpp"
-#include "FunctionApproximation/FunctionApproximation.hpp"
-#include "FunctionApproximation/FunctionModel.hpp"
-#include "FunctionApproximation/TrainingData.hpp"
-#include "FunctionApproximation/writeDistanceEnergyTable.hpp"
-#include "Potentials/CompoundPotential.hpp"
-#include "Potentials/Exceptions.hpp"
-#include "Potentials/PotentialDeserializer.hpp"
+#include "Potentials/EmpiricalPotential.hpp"
 #include "Potentials/PotentialFactory.hpp"
 #include "Potentials/PotentialRegistry.hpp"
 #include "Potentials/PotentialSerializer.hpp"
+#include "Potentials/PotentialTrainer.hpp"
 #include "Potentials/SerializablePotential.hpp"
+#include "World.hpp"
 
 using namespace MoleCuilder;
@@ -75 +70 @@
 /** =========== define the function ====================== */
 
-HomologyGraph getFirstGraphwithSpecifiedElements(
-    const HomologyContainer &homologies,
-    const SerializablePotential::ParticleTypes_t &types)
-{
-  ASSERT( !types.empty(),
-      "getFirstGraphwithSpecifiedElements() - charges is empty?");
-  // create charges
-  Fragment::charges_t charges;
-  charges.resize(types.size());
-  std::transform(types.begin(), types.end(),
-      charges.begin(), boost::lambda::_1);
-  // convert into count map
-  Extractors::elementcounts_t counts_per_charge =
-      Extractors::_detail::getElementCounts(charges);
-  ASSERT( !counts_per_charge.empty(),
-      "getFirstGraphwithSpecifiedElements() - charge counts are empty?");
-  LOG(2, "DEBUG: counts_per_charge is " << counts_per_charge << ".");
-  // we want to check each (unique) key only once
-  for (HomologyContainer::const_key_iterator iter = homologies.key_begin();
-      iter != homologies.key_end(); iter = homologies.getNextKey(iter)) {
-    // check if every element has the right number of counts
-    Extractors::elementcounts_t::const_iterator countiter = counts_per_charge.begin();
-    for (; countiter != counts_per_charge.end(); ++countiter)
-      if (!(*iter).hasTimesAtomicNumber(
-          static_cast<size_t>(countiter->first),
-          static_cast<size_t>(countiter->second))
-          )
-        break;
-    if( countiter == counts_per_charge.end())
-      return *iter;
-  }
-  return HomologyGraph();
-}
-
-SerializablePotential::ParticleTypes_t getNumbersFromElements(
-    const std::vector<const element *> &fragment)
-{
-  SerializablePotential::ParticleTypes_t fragmentnumbers;
-  std::transform(fragment.begin(), fragment.end(), std::back_inserter(fragmentnumbers),
-      boost::bind(&element::getAtomicNumber, _1));
-  return fragmentnumbers;
-}
-
-
 ActionState::ptr PotentialFitPotentialAction::performCall() {
   // fragment specifies the homology fragment to use
   SerializablePotential::ParticleTypes_t fragmentnumbers =
-      getNumbersFromElements(params.fragment.get());
+      PotentialTrainer::getNumbersFromElements(params.fragment.get());
 
   // either charges and a potential is specified or a file
-  if (boost::filesystem::exists(params.potential_file.get())) {
-    std::ifstream returnstream(params.potential_file.get().string().c_str());
-    if (returnstream.good()) {
-      try {
-        PotentialDeserializer deserialize(returnstream);
-        deserialize();
-      } catch (SerializablePotentialMissingValueException &e) {
-        if (const std::string *key = boost::get_error_info<SerializablePotentialKey>(e))
-          STATUS("Missing value when parsing information for potential "+*key+".");
-        else
-          STATUS("Missing value parsing information for potential with unknown key.");
-        return Action::failure;
-      } catch (SerializablePotentialIllegalKeyException &e) {
-        if (const std::string *key = boost::get_error_info<SerializablePotentialKey>(e))
-          STATUS("Illegal key parsing information for potential "+*key+".");
-        else
-          STATUS("Illegal key parsing information for potential with unknown key.");
-        return Action::failure;
-      }
-    } else {
-      STATUS("Failed to parse from "+params.potential_file.get().string()+".");
-      return Action::failure;
+  if (params.charges.get().empty()) {
+    STATUS("No charges given!");
+    return Action::failure;
+  } else {
+    // charges specify the potential type
+    SerializablePotential::ParticleTypes_t chargenumbers =
+        PotentialTrainer::getNumbersFromElements(params.charges.get());
+
+    LOG(0, "STATUS: I'm training now a " << params.potentialtype.get()
+        << " potential on charges " << chargenumbers << " on data from World's homologies.");
+
+    // register desired potential and an additional constant one
+    {
+      EmpiricalPotential *potential =
+          PotentialFactory::getInstance().createInstance(
+              params.potentialtype.get(),
+              chargenumbers);
+      // check whether such a potential already exists
+      const std::string potential_name = potential->getName();
+      if (PotentialRegistry::getInstance().isPresentByName(potential_name)) {
+        delete potential;
+        potential = PotentialRegistry::getInstance().getByName(potential_name);
+      } else
+        PotentialRegistry::getInstance().registerInstance(potential);
     }
-    returnstream.close();
-
-    LOG(0, "STATUS: I'm training now a set of potentials parsed from "
-        << params.potential_file.get().string() << " on a fragment "
-        << fragmentnumbers << " on data from World's homologies.");
-
-  } else {
-    if (params.charges.get().empty()) {
-      STATUS("Neither charges nor potential file given!");
-      return Action::failure;
-    } else {
-      // charges specify the potential type
-      SerializablePotential::ParticleTypes_t chargenumbers =
-          getNumbersFromElements(params.charges.get());
-
-      LOG(0, "STATUS: I'm training now a " << params.potentialtype.get()
-          << " potential on charges " << chargenumbers << " on data from World's homologies.");
-
-      // register desired potential and an additional constant one
-      {
-        EmpiricalPotential *potential =
-            PotentialFactory::getInstance().createInstance(
-                params.potentialtype.get(),
-                chargenumbers);
-        // check whether such a potential already exists
-        const std::string potential_name = potential->getName();
-        if (PotentialRegistry::getInstance().isPresentByName(potential_name)) {
-          delete potential;
-          potential = PotentialRegistry::getInstance().getByName(potential_name);
-        } else
-          PotentialRegistry::getInstance().registerInstance(potential);
-      }
-      {
-        EmpiricalPotential *constant =
-            PotentialFactory::getInstance().createInstance(
-                std::string("constant"),
-                SerializablePotential::ParticleTypes_t());
-        // check whether such a potential already exists
-        const std::string constant_name = constant->getName();
-        if (PotentialRegistry::getInstance().isPresentByName(constant_name)) {
-          delete constant;
-          constant = PotentialRegistry::getInstance().getByName(constant_name);
-        } else
-          PotentialRegistry::getInstance().registerInstance(constant);
-      }
+    {
+      EmpiricalPotential *constant =
+          PotentialFactory::getInstance().createInstance(
+              std::string("constant"),
+              SerializablePotential::ParticleTypes_t());
+      // check whether such a potential already exists
+      const std::string constant_name = constant->getName();
+      if (PotentialRegistry::getInstance().isPresentByName(constant_name)) {
+        delete constant;
+        constant = PotentialRegistry::getInstance().getByName(constant_name);
+      } else
+        PotentialRegistry::getInstance().registerInstance(constant);
     }
   }
 
   // parse homologies into container
-  HomologyContainer &homologies = World::getInstance().getHomologies();
+  const HomologyContainer &homologies = World::getInstance().getHomologies();
 
   // first we try to look into the HomologyContainer
@@ -211 +131 @@
 
   // then we ought to pick the right HomologyGraph ...
-  const HomologyGraph graph = getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers);
+  const HomologyGraph graph =
+      PotentialTrainer::getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers);
   if (graph != HomologyGraph()) {
     LOG(1, "First representative graph containing fragment "
@@ -220 +141 @@
   }
 
-  // fit potential
-  FunctionModel *model = new CompoundPotential(graph);
-  ASSERT( model != NULL,
-      "PotentialFitPotentialAction::performCall() - model is NULL.");
-
-  /******************** TRAINING ********************/
-  // fit potential
-  FunctionModel::parameters_t bestparams(model->getParameterDimension(), 0.);
-  {
-    // Afterwards we go through all of this type and gather the distance and the energy value
-    TrainingData data(model->getSpecificFilter());
-    data(homologies.getHomologousGraphs(graph));
-
-    // print distances and energies if desired for debugging
-    if (!data.getTrainingInputs().empty()) {
-      // print which distance is which
-      size_t counter=1;
-      if (DoLog(3)) {
-        const FunctionModel::arguments_t &inputs = data.getAllArguments()[0];
-        for (FunctionModel::arguments_t::const_iterator iter = inputs.begin();
-            iter != inputs.end(); ++iter) {
-          const argument_t &arg = *iter;
-          LOG(3, "DEBUG: distance " << counter++ << " is between (#"
-              << arg.indices.first << "c" << arg.types.first << ","
-              << arg.indices.second << "c" << arg.types.second << ").");
-        }
-      }
-
-      // print table
-      if (params.training_file.get().string().empty()) {
-        LOG(3, "DEBUG: I gathered the following training data:\n" <<
-            _detail::writeDistanceEnergyTable(data.getDistanceEnergyTable()));
-      } else {
-        std::ofstream trainingstream(params.training_file.get().string().c_str());
-        if (trainingstream.good()) {
-          LOG(3, "DEBUG: Writing training data to file " <<
-              params.training_file.get().string() << ".");
-          trainingstream << _detail::writeDistanceEnergyTable(data.getDistanceEnergyTable());
-        }
-        trainingstream.close();
-      }
-    }
-
-    if ((params.threshold.get() < 1) && (params.best_of_howmany.isSet()))
-      ELOG(2, "threshold parameter always overrules max_runs, both are specified.");
-    // now perform the function approximation by optimizing the model function
-    FunctionApproximation approximator(data, *model);
-    if (model->isBoxConstraint() && approximator.checkParameterDerivatives()) {
-      double l2error = std::numeric_limits<double>::max();
-      // seed with current time
-      srand((unsigned)time(0));
-      unsigned int runs=0;
-      // threshold overrules max_runs
-      const double threshold = params.threshold.get();
-      const unsigned int max_runs = (threshold >= 1.) ?
-          (params.best_of_howmany.isSet() ? params.best_of_howmany.get() : 1) : 0;
-      LOG(1, "INFO: Maximum runs is " << max_runs << " and threshold set to " << threshold << ".");
-      do {
-        // generate new random initial parameter values
-        model->setParametersToRandomInitialValues(data);
-        LOG(1, "INFO: Initial parameters of run " << runs << " are "
-            << model->getParameters() << ".");
-        approximator(FunctionApproximation::ParameterDerivative);
-        LOG(1, "INFO: Final parameters of run " << runs << " are "
-            << model->getParameters() << ".");
-        const double new_l2error = data.getL2Error(*model);
-        if (new_l2error < l2error) {
-          // store currently best parameters
-          l2error = new_l2error;
-          bestparams = model->getParameters();
-          LOG(1, "STATUS: New fit from run " << runs
-              << " has better error of " << l2error << ".");
-        }
-      } while (( ++runs < max_runs) || (l2error > threshold));
-      // reset parameters from best fit
-      model->setParameters(bestparams);
-      LOG(1, "INFO: Best parameters with L2 error of "
-          << l2error << " are " << model->getParameters() << ".");
-    } else {
-      STATUS("No required parameter derivatives for a box constraint minimization known.");
-      return Action::failure;
-    }
-
-    // create a map of each fragment with error.
-    HomologyContainer::range_t fragmentrange = homologies.getHomologousGraphs(graph);
-    TrainingData::L2ErrorConfigurationIndexMap_t WorseFragmentMap =
-        data.getWorstFragmentMap(*model, fragmentrange);
-    LOG(0, "RESULT: WorstFragmentMap " << WorseFragmentMap << ".");
-
-    // print fitted potentials
-    std::stringstream potentials;
-    PotentialSerializer serialize(potentials);
-    serialize();
-    LOG(1, "STATUS: Resulting parameters are " << std::endl << potentials.str());
-    std::ofstream returnstream(params.potential_file.get().string().c_str());
-    if (returnstream.good()) {
-      returnstream << potentials.str();
-    }
+  // training
+  PotentialTrainer trainer;
+  const bool status = trainer(
+      homologies,
+      graph,
+      params.training_file.get(),
+      params.threshold.get(),
+      params.best_of_howmany.get());
+  if (!status) {
+    STATUS("No required parameter derivatives for a box constraint minimization known.");
+    return Action::failure;
   }
-  delete model;
 
   return Action::success;