Changeset 98d166 for src

Timestamp:

Oct 19, 2014, 5:08:36 PM (10 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

f1eabd

Parents:

9bce96

git-author:

Frederik Heber <heber@…> (09/11/14 20:38:05)

git-committer:

Frederik Heber <heber@…> (10/19/14 17:08:36)

Message:

Split FitPotentialAction into two actions.

• FitPotentialAction fits just a single potential.
• new FitCompoundPotentialAction fits via a potential file.
• PotentialTrainer contains all helper functions and common code and is only compiled when CONDLEVMAR.

Location:

src

Files:

• Actions/GlobalListOfActions.hpp (modified) (1 diff)
• Actions/Makefile.am (modified) (1 diff)
• Actions/PotentialAction/FitCompoundPotentialAction.cpp (added)
• Actions/PotentialAction/FitCompoundPotentialAction.def (added)
• Actions/PotentialAction/FitCompoundPotentialAction.hpp (added)
• Actions/PotentialAction/FitPotentialAction.cpp (modified) (4 diffs)
• Actions/PotentialAction/FitPotentialAction.def (modified) (2 diffs)
• Potentials/Makefile.am (modified) (1 diff)
• Potentials/PotentialTrainer.cpp (added)
• Potentials/PotentialTrainer.hpp (added)

src/Actions/GlobalListOfActions.hpp ¶

@@ -148 +148 @@
 #define GLOBALLISTOFACTIONS_LEVMAR \
     BOOST_PP_SEQ_PUSH_BACK( \
-        GLOBALLISTOFACTIONS_initial, \
-        PotentialFitPotential \
+                BOOST_PP_SEQ_PUSH_BACK( \
+                                GLOBALLISTOFACTIONS_initial, \
+                                PotentialFitPotential \
+                ), \
+        PotentialFitCompoundPotential \
     )
 #else

src/Actions/Makefile.am ¶

@@ -363 +363 @@
 if CONDLEVMAR
 POTENTIALACTIONSOURCE += \
+  Actions/PotentialAction/FitCompoundPotentialAction.cpp \
   Actions/PotentialAction/FitPotentialAction.cpp
 POTENTIALACTIONHEADER += \
+  Actions/PotentialAction/FitCompoundPotentialAction.hpp \
   Actions/PotentialAction/FitPotentialAction.hpp
 POTENTIALACTIONDEFS += \
+  Actions/PotentialAction/FitCompoundPotentialAction.def \
   Actions/PotentialAction/FitPotentialAction.def
 endif

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
+  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;
+  }
 
-  /******************** 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();
-    }
+  // store the potentials
+  // print fitted potentials
+  std::stringstream potentials;
+  PotentialSerializer serialize(potentials);
+  serialize();
+  LOG(1, "STATUS: Resulting parameters are " << std::endl << potentials.str());
+  if (!params.potential_file.get().string().empty()) {
+    PotentialTrainer::storePotentials(params.potential_file.get());
   }
-  delete model;
 
   return Action::success;

src/Actions/PotentialAction/FitPotentialAction.def ¶

@@ -20 +20 @@
 #include "Parameters/Validators/Specific/FileSuffixValidator.hpp"
 #include "Parameters/Validators/Specific/FilePresentValidator.hpp"
+#include "Parameters/Validators/Specific/FileSuffixValidator.hpp"
 #include "Parameters/Validators/Specific/PotentialTypeValidator.hpp"
 
@@ -25 +26 @@
 // ValueStorage by the token "Z" -> first column: int, Z, "Z"
 // "undefine" if no parameters are required, use (NOPARAM_DEFAULT) for each (undefined) default value
-#define paramtypes (boost::filesystem::path)(std::string)(boost::filesystem::path)(std::vector<const element *>)(std::vector<const element *>)(unsigned int)(double)
-#define paramtokens ("training-file")("potential-type")("potential-file")("potential-charges")("fragment-charges")("take-best-of")("set-threshold")
-#define paramdescriptions ("optional file to write training data to")("potential type to fit")("optional potential file specifying multiple potentials to fit")("charges specifying the potential")("charges specifying the fragment")("take the best among this many approximations")("Require L2 error to be smaller than threshold, overrides number of attempts")
-#define paramdefaults (PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(PARAM_DEFAULT(3))(PARAM_DEFAULT(1.))
-#define paramreferences (training_file)(potentialtype)(potential_file)(charges)(fragment)(best_of_howmany)(threshold)
+#define paramtypes (boost::filesystem::path)(boost::filesystem::path)(std::string)(std::vector<const element *>)(std::vector<const element *>)(unsigned int)(double)
+#define paramtokens ("training-file")("potential-file")("potential-type")("potential-charges")("fragment-charges")("take-best-of")("set-threshold")
+#define paramdescriptions ("optional file to write training data to")("optional potential file to wrote resulting potentials to")("potential type to fit")("charges specifying the potential")("charges specifying the fragment")("take the best among this many approximations")("Require L2 error to be smaller than threshold, overrides number of attempts")
+#define paramdefaults (PARAM_DEFAULT(""))(PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(PARAM_DEFAULT(3))(PARAM_DEFAULT(1.))
+#define paramreferences (training_file)(potential_file)(potentialtype)(charges)(fragment)(best_of_howmany)(threshold)
 #define paramvalids \
 (DummyValidator<boost::filesystem::path>()) \
-(EmptyStringValidator() || PotentialTypeValidator()) \
 (!FilePresentValidator() || FileSuffixValidator("potentials")) \
+(PotentialTypeValidator()) \
 (STLVectorValidator< std::vector<const element *> >(0,99, ElementValidator())) \
 (STLVectorValidator< std::vector<const element *> >(1,99, ElementValidator())) \

src/Potentials/Makefile.am ¶

@@ -53 +53 @@
   Potentials/Specifics/PotentialTypes.def \
   Potentials/Specifics/PotentialTypes.undef
+
+if CONDLEVMAR
+POTENTIALSHEADER += \
+  Potentials/PotentialTrainer.hpp
+POTENTIALSSOURCE = \
+  Potentials/PotentialTrainer.cpp
+endif
 
 # add here headers for real potentials (i.e. the PotentialFactory should instantiate)