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

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 "Potentials/EmpiricalPotential.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/PotentialFactory.hpp"
 #include "Potentials/PotentialRegistry.hpp"
 #include "Potentials/PotentialSerializer.hpp"
-#include "Potentials/PotentialTrainer.hpp"
 #include "Potentials/SerializablePotential.hpp"
-#include "World.hpp"
 
 using namespace MoleCuilder;
@@ -70 +75 @@
 /** =========== 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 =
-      PotentialTrainer::getNumbersFromElements(params.fragment.get());
+      getNumbersFromElements(params.fragment.get());
 
   // either charges and a potential is specified or a file
-  if (params.charges.get().empty()) {
-    STATUS("No charges given!");
-    return Action::failure;
+  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;
+    }
+    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 {
-    // 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);
-    }
-    {
-      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);
+    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);
+      }
     }
   }
 
   // parse homologies into container
-  const HomologyContainer &homologies = World::getInstance().getHomologies();
+  HomologyContainer &homologies = World::getInstance().getHomologies();
 
   // first we try to look into the HomologyContainer
@@ -131 +211 @@
 
   // then we ought to pick the right HomologyGraph ...
-  const HomologyGraph graph =
-      PotentialTrainer::getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers);
+  const HomologyGraph graph = getFirstGraphwithSpecifiedElements(homologies,fragmentnumbers);
   if (graph != HomologyGraph()) {
     LOG(1, "First representative graph containing fragment "
@@ -141 +220 @@
   }
 
-  // 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;
-  }
+  // 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();
+    }
+  }
+  delete model;
 
   return Action::success;