Changeset 7c1091

Timestamp:

Jul 8, 2013, 2:22:03 PM (11 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:

baccf6

Parents:

c5e63a

git-author:

Frederik Heber <heber@…> (05/09/13 19:25:45)

git-committer:

Frederik Heber <heber@…> (07/08/13 14:22:03)

Message:

Extended CompoundPotential for multiple occurences of argument bunches.

• added HomologyGraph::hasGreaterEqualTimesAtomicNumber() has model might have multiple matching arguments within a single fragment.
• i.e. when there is more than one argument bunch for a given model. Hence, we check the stream whether upcoming arguments match current and next model. Note however that models are not necessarily sorted by particle_types' size. Hence, we must trick a little and need a specific map comparator.

Location:

src

Files:

• Fragmentation/Homology/HomologyGraph.cpp (modified) (1 diff)
• Fragmentation/Homology/HomologyGraph.hpp (modified) (1 diff)
• Potentials/CompoundPotential.cpp (modified) (11 diffs)
• Potentials/CompoundPotential.hpp (modified) (1 diff)
• Potentials/unittests/CompoundPotentialUnitTest.cpp (modified) (8 diffs)

src/Fragmentation/Homology/HomologyGraph.cpp

@@ -114 +114 @@
 }
 
+bool HomologyGraph::hasGreaterEqualTimesAtomicNumber(const size_t _number, const size_t _times) const
+{
+  size_t count = 0;
+  for (nodes_t::const_iterator iter = nodes.begin(); iter != nodes.end(); ++iter) {
+    if ((iter->first).getAtomicNumber() == _number)
+      count += iter->second;
+  }
+  return (count >= _times);
+}
+
 std::ostream& operator<<(std::ostream& ost, const HomologyGraph &graph)
 {

src/Fragmentation/Homology/HomologyGraph.hpp

@@ -129 +129 @@
   bool hasTimesAtomicNumber(const size_t _number, const size_t _times) const;
 
+  /** Checks whether this graph has \b greater equal \a _times nodes with \a _number
+   * atomic number.
+   *
+   * @param _number desired atomic number
+   * @param _times number this must occur
+   * @return true - graph has greater equal \a _times nodes with \a _number, false - else
+   */
+  bool hasGreaterEqualTimesAtomicNumber(const size_t _number, const size_t _times) const;
+
   /** Assignment operator for class HomologyGraph.
    *

src/Potentials/CompoundPotential.cpp

@@ -84 +84 @@
     Extractors::elementcounts_t::const_iterator countiter = counts_per_charge.begin();
     for (; countiter != counts_per_charge.end(); ++countiter)
-      if (!graph.hasTimesAtomicNumber(
+      if (!graph.hasGreaterEqualTimesAtomicNumber(
           static_cast<size_t>(countiter->first),
           static_cast<size_t>(countiter->second))
@@ -115 +115 @@
   parameters_t::const_iterator iter = _params.begin();
   BOOST_FOREACH( FunctionModel* model, models) {
-    const parameters_t &model_params = model->getParameters();
-    const size_t model_dim = model_params.size();
+    const size_t model_dim = model->getParameterDimension();
     if (dim > 0) {
       parameters_t subparams;
@@ -131 +130 @@
     }
   }
+
+#ifndef NDEBUG
+  parameters_t check_params(getParameters());
+  check_params.resize(_params.size()); // truncate to same size
+  ASSERT( check_params == _params,
+      "CompoundPotential::setParameters() - failed, mismatch in to be set "
+      +toString(_params)+" and set "+toString(check_params)+" params.");
+#endif
 }
 
@@ -140 +147 @@
   BOOST_FOREACH( const FunctionModel* model, models) {
     const CompoundPotential::parameters_t &params = model->getParameters();
-    std::copy(params.begin(), params.end(), iter);
     ASSERT( iter != parameters.end(),
         "CompoundPotential::getParameters() - iter already at end.");
-  }
+    iter = std::copy(params.begin(), params.end(), iter);
+  }
+  ASSERT( iter == parameters.end(),
+      "CompoundPotential::getParameters() - iter not at end.");
   return parameters;
 }
@@ -169 +178 @@
 }
 
-std::vector<CompoundPotential::arguments_t> CompoundPotential::splitUpArgumentsByModels(
+bool CompoundPotential::areValidArguments(
+    const SerializablePotential::ParticleTypes_t &_types,
+    const arguments_t &args) const
+{
+  // /this function does much the same as Extractors::reorderArgumentsByParticleTypes()
+  typedef std::list< argument_t > ListArguments_t;
+  ListArguments_t availableList(args.begin(), args.end());
+
+  /// basically, we have two choose any two pairs out of types but only those
+  /// where the first is less than the letter. Hence, we start the second
+  /// iterator at the current position of the first one and skip the equal case.
+  for (SerializablePotential::ParticleTypes_t::const_iterator firstiter = _types.begin();
+      firstiter != _types.end();
+      ++firstiter) {
+    for (SerializablePotential::ParticleTypes_t::const_iterator seconditer = firstiter;
+        seconditer != _types.end();
+        ++seconditer) {
+      if (seconditer == firstiter)
+        continue;
+
+      // search the right one in _args (we might allow switching places of
+      // firstiter and seconditer, as distance is symmetric).
+      // we remove the matching argument to make sure we don't pick it twice
+      ListArguments_t::iterator iter = availableList.begin();
+      for (;iter != availableList.end(); ++iter) {
+        LOG(3, "DEBUG: Current args is " << *iter << ".");
+        if ((iter->types.first == *firstiter)
+              && (iter->types.second == *seconditer)) {
+          availableList.erase(iter);
+          break;
+        }
+        else if ((iter->types.first == *seconditer)
+              && (iter->types.second == *firstiter)) {
+          availableList.erase(iter);
+          break;
+        }
+      }
+      if ( iter == availableList.end())
+        return false;
+    }
+  }
+
+  return true;
+}
+
+CompoundPotential::arguments_by_model_t CompoundPotential::splitUpArgumentsByModels(
     const arguments_t &arguments) const
 {
-  std::vector<arguments_t> partial_args(models.size());
-  std::vector<arguments_t>::iterator partialiter = partial_args.begin();
+  arguments_by_model_t partial_args;
   arguments_t::const_iterator argiter = arguments.begin();
-  BOOST_FOREACH( const SerializablePotential::ParticleTypes_t &types, particletypes_per_model) {
-    // we always expect N(N-1)/2 distances for N particle types
-    arguments_t::const_iterator enditer = argiter+(types.size()*(types.size()-1)/2);
-    ASSERT( argiter != arguments.end(),
-        "CompoundPotential::splitUpArgumentsByModels() - incorrect number of arguments.");
-    std::copy(argiter, enditer, std::back_inserter(*partialiter++));
-    argiter = enditer;
-  }
-  ASSERT( argiter == arguments.end(),
-      "CompoundPotential::splitUpArgumentsByModels() - incorrect number of arguments.");
+  particletypes_per_model_t::const_iterator typesiter = particletypes_per_model.begin();
+  models_t::const_iterator modeliter = models.begin();
+
+  // add constant model (which is always first model) with empty args if present
+  if (typesiter->empty()) {
+    partial_args.push_back(
+        std::pair<FunctionModel *, arguments_t>(*modeliter, arguments_t())
+        );
+    ++modeliter;
+    ++typesiter;
+  }
+  // then check other models
+  while (argiter != arguments.end()) {
+    if (typesiter+1 != particletypes_per_model.end()) {
+      // check whether next argument bunch is for same model or different one
+      // we extract both partial_arguments, if the latter fits, we take the latter.
+      const SerializablePotential::ParticleTypes_t &types = *typesiter;
+      const SerializablePotential::ParticleTypes_t &nexttypes = *(typesiter+1);
+
+      // we always expect N(N-1)/2 distances for N particle types
+      arguments_t::const_iterator enditer = argiter+(types.size()*(types.size()-1)/2);
+      arguments_t::const_iterator nextenditer = argiter+(nexttypes.size()*(nexttypes.size()-1)/2);
+      arguments_t args(argiter, enditer);
+      arguments_t nextargs(argiter, nextenditer);
+      if (types.size() < nexttypes.size()) {
+        if (areValidArguments(nexttypes, nextargs)) {
+          // latter matches, increment
+          ++typesiter;
+          partial_args.push_back(
+              std::make_pair(*(++modeliter), arguments_t(argiter, nextenditer))
+              );
+          argiter = nextenditer;
+        } else if (areValidArguments(types, args)) {
+          // only former matches, don't increment
+          partial_args.push_back(
+              std::make_pair(*modeliter, arguments_t(argiter, enditer))
+              );
+          argiter = enditer;
+        } else
+          ASSERT(0,
+              "CompoundPotential::splitUpArgumentsByModels() - neither type matches its argument bunch.");
+      } else {
+        if (areValidArguments(types, args)) {
+          // only former matches, don't increment
+          partial_args.push_back(
+              std::make_pair(*modeliter, arguments_t(argiter, enditer))
+              );
+          argiter = enditer;
+        } else if (areValidArguments(nexttypes, nextargs)) {
+          // latter matches, increment
+          ++typesiter;
+          partial_args.push_back(
+              std::make_pair(*(++modeliter), arguments_t(argiter, nextenditer))
+              );
+          argiter = nextenditer;
+        }
+      }
+    } else {
+      const SerializablePotential::ParticleTypes_t &types = *typesiter;
+      // we always expect N(N-1)/2 distances for N particle types
+      arguments_t::const_iterator enditer = argiter+(types.size()*(types.size()-1)/2);
+      partial_args.push_back(
+          std::make_pair(*modeliter, arguments_t(argiter, enditer))
+          );
+      argiter = enditer;
+    }
+  }
+
   return partial_args;
 }
@@ -190 +301 @@
 CompoundPotential::results_t CompoundPotential::operator()(const arguments_t &arguments) const
 {
-  // first, we have to split up the given arguments
-  std::vector<arguments_t> partial_args =
+  /// first, we have to split up the given arguments
+  arguments_by_model_t partial_args =
       splitUpArgumentsByModels(arguments);
-  // then, with each bunch of arguments, we call the specific model
+  // print split up argument list for debugging
+  if (DoLog(4)) {
+    LOG(4, "Arguments by model are: ");
+    for(arguments_by_model_t::const_iterator iter = partial_args.begin();
+        iter != partial_args.end(); ++iter) {
+      LOG(4, "\tModel with " << iter->first->getParameterDimension()
+          << " parameters " << iter->first->getParameters()
+          << " and arguments: " << iter->second);
+    }
+  }
+
+  /// then, with each bunch of arguments, we call the specific model
   results_t results(1,0.);
-  std::vector<results_t> partial_results(models.size());
-  std::transform(
-      models.begin(), models.end(),
-      partial_args.begin(),
-      partial_results.begin(),
-      boost::bind(&FunctionModel::operator(), _1, _2)
-  );
-  std::for_each(partial_results.begin(), partial_results.end(),
-      std::cout << (boost::lambda::_1)[0] << "\t");
+  std::vector<results_t> partial_results;
+  for(arguments_by_model_t::const_iterator iter = partial_args.begin();
+      iter != partial_args.end(); ++iter) {
+    partial_results.push_back(
+        (*iter->first)(iter->second)
+    );
+  }
+  // print partial results for debugging
+  if (DoLog(4)) {
+    std::stringstream output;
+    output << "Partial results are: ";
+    std::for_each(partial_results.begin(), partial_results.end(),
+        output << (boost::lambda::_1)[0] << "\t");
+    LOG(4, output.str());
+  }
+
+  /// Finally, sum up all results and return
   std::for_each(partial_results.begin(), partial_results.end(),
       results[0] += (boost::lambda::_1)[0]);
@@ -212 +342 @@
 {
   // first, we have to split up the given arguments
-  std::vector<arguments_t> partial_args =
+  arguments_by_model_t partial_args =
       splitUpArgumentsByModels(arguments);
   // then, with each bunch of arguments, we call the specific model
@@ -223 +353 @@
   std::partial_sum(dimensions.begin(), dimensions.end(), dimensions.begin());
 
-  // look for value not less than index
+  // look for first value greater than index
   std::vector<size_t>::const_iterator iter =
-      std::lower_bound(dimensions.begin(), dimensions.end(), index);
+      std::upper_bound(dimensions.begin(), dimensions.end(), index);
 
   // step forward to same model
@@ -231 +361 @@
   std::advance(modeliter,
       std::distance(const_cast<const std::vector<size_t> &>(dimensions).begin(), iter) );
-  std::vector<arguments_t>::const_iterator argiter = partial_args.begin();
-  std::advance(argiter,
-      std::distance(const_cast<const std::vector<size_t> &>(dimensions).begin(), iter) );
-
-  // evaluate with correct relative index and return
-  const size_t indexbase = (iter == dimensions.begin()) ? 0 : *(--iter);
-  CompoundPotential::results_t results =
-      (*modeliter)->parameter_derivative(*argiter, index-indexbase);
-  return results;
+
+  CompoundPotential::results_t returnresults;
+  for(arguments_by_model_t::const_iterator argiter = partial_args.begin();
+      argiter != partial_args.end(); ++argiter) {
+    const FunctionModel *model = argiter->first;
+
+    // for every matching model evaluate
+    if (model == *modeliter) {
+      // evaluate with correct relative index and return
+      const size_t indexbase = (iter == dimensions.begin()) ? 0 : *(iter-1);
+      CompoundPotential::results_t results =
+          model->parameter_derivative(argiter->second, index-indexbase);
+
+      // either set results or add
+      if (returnresults.empty())
+        returnresults = results;
+      else
+        std::transform(
+            results.begin(), results.end(),
+            returnresults.begin(),
+            returnresults.begin(),
+            std::plus<FunctionModel::result_t>());
+    }
+  }
+
+  return returnresults;
 }
 
@@ -258 +405 @@
   BOOST_FOREACH( FunctionModel* model, models) {
     const CompoundPotential::parameters_t params = model->getLowerBoxConstraints();
-    std::copy(params.begin(), params.end(), iter);
     ASSERT( iter != constraints.end(),
-        "CompoundPotential::getParameters() - iter already at end.");
-  }
+        "CompoundPotential::getLowerBoxConstraints() - iter already at end.");
+    iter = std::copy(params.begin(), params.end(), iter);
+  }
+  ASSERT( iter == constraints.end(),
+      "CompoundPotential::getLowerBoxConstraints() - iter not at end.");
   return constraints;
 }
@@ -272 +421 @@
   BOOST_FOREACH( FunctionModel* model, models) {
     const CompoundPotential::parameters_t params = model->getUpperBoxConstraints();
-    std::copy(params.begin(), params.end(), iter);
     ASSERT( iter != constraints.end(),
-        "CompoundPotential::getParameters() - iter already at end.");
-  }
+        "CompoundPotential::getUpperBoxConstraints() - iter already at end.");
+    iter = std::copy(params.begin(), params.end(), iter);
+  }
+  ASSERT( iter == constraints.end(),
+      "CompoundPotential::getUpperBoxConstraints() - iter not at end.");
   return constraints;
 }

src/Potentials/CompoundPotential.hpp

@@ -132 +132 @@
 
 private:
+  //!> typedef for split up arguments, each associated to a model
+  typedef std::vector< std::pair<FunctionModel *, arguments_t> > arguments_by_model_t;
+
   /** Helper function to split up concatenated arguments for operator() calls.
    *
    * \param arguments arguments to split up
-   * \return vector of partial arguments
+   * \return vector of partial arguments with associated model
    */
-  std::vector<arguments_t> splitUpArgumentsByModels(const arguments_t &arguments) const;
+  arguments_by_model_t splitUpArgumentsByModels(const arguments_t &arguments) const;
+
+  /** Helper function to check whether split up argument bunch matches with types.
+   *
+   * \param types types of potential to check whether args match
+   * \param args vector of argument whose types to check
+   */
+  bool areValidArguments(
+      const SerializablePotential::ParticleTypes_t &types,
+      const arguments_t &args) const;
 
 private:

src/Potentials/unittests/CompoundPotentialUnitTest.cpp

@@ -92 +92 @@
   }
 
-  // create graph
+  // create graph (i.e. this simulates a water molecule)
   {
     // add nodes
     nodes +=
-        std::make_pair(FragmentNode(0,1),1),
-        std::make_pair(FragmentNode(1,1),1);
+        std::make_pair(FragmentNode(0,2),1),
+        std::make_pair(FragmentNode(1,1),2);
 
     // add edges
     edges +=
-        std::make_pair(FragmentEdge(0,1),1);
+        std::make_pair(FragmentEdge(0,1),2);
 
     // construct graph
@@ -125 +125 @@
       4.46595;
   output +=
-      -78.5211,
-      -78.8049,
-      -78.935,
-      -78.9822,
-      -78.9867,
-      -78.971,
-      -78.9475,
-      -78.9225,
-      -78.899,
-      -78.8784;
+      2.*0.467226+d,
+      2.*0.183388+d,
+      2.*0.0532649+d,
+      2.*0.00609808+d,
+      2.*0.00160056+d,
+      2.*0.0172506+d,
+      2.*0.0407952+d,
+      2.*0.0658475+d,
+      2.*0.0893157+d,
+      2.*0.109914+d;
 
   CPPUNIT_ASSERT_EQUAL( input.size(), output.size() );
@@ -157 +157 @@
   compound.setParameters(params);
   for (size_t index = 0; index < input.size(); ++index) {
-    argument_t arg(argument_t::indices_t(0,1), argument_t::types_t(0,1), input[index]);
-    const double result = compound( FunctionModel::arguments_t(1,arg) )[0];
+    argument_t firstarg(argument_t::indices_t(0,1), argument_t::types_t(0,1), input[index]);
+    argument_t secondarg(argument_t::indices_t(0,2), argument_t::types_t(0,1), input[index]);
+    FunctionModel::arguments_t args;
+    args += firstarg,secondarg;
+    const double result = compound( args )[0];
     CPPUNIT_ASSERT(
         Helpers::isEqual(
@@ -177 +180 @@
 //  params += d,a,c,D;
 //  compound.setParameters(params);
-//  argument_t arg(argument_t::indices_t(0,1), argument_t::types_t(0,1), c);
-//  const double result = compound.derivative(FunctionModel::arguments_t(1,arg))[0]
+//  argument_t firstarg(argument_t::indices_t(0,1), argument_t::types_t(0,1), input[index]);
+//  argument_t secondarg(argument_t::indices_t(0,2), argument_t::types_t(0,1), input[index]);
+//  FunctionModel::arguments_t args;
+//  args += firstarg,secondarg;
+//  const double result = compound.derivative( args )[0]
 //  CPPUNIT_ASSERT(
 //      Helpers::isEqual(
@@ -197 +203 @@
   params += d,a,c,D;
   compound.setParameters(params);
-  argument_t arg(argument_t::indices_t(0,1), argument_t::types_t(0,1), c);
-  {
-    const double result =
-        compound.parameter_derivative(
-            FunctionModel::arguments_t(1,arg),
+  argument_t firstarg(argument_t::indices_t(0,1), argument_t::types_t(0,1), c);
+  argument_t secondarg(argument_t::indices_t(0,2), argument_t::types_t(0,1), c);
+  FunctionModel::arguments_t args;
+  args += firstarg,secondarg;
+  {
+    const double result =
+        compound.parameter_derivative(
+            args,
             0)[0];
     CPPUNIT_ASSERT(
@@ -214 +223 @@
     const double result =
         compound.parameter_derivative(
-            FunctionModel::arguments_t(1,arg),
+            args,
             1)[0];
     CPPUNIT_ASSERT(
@@ -227 +236 @@
     const double result =
         compound.parameter_derivative(
-            FunctionModel::arguments_t(1,arg),
+            args,
             2)[0];
     CPPUNIT_ASSERT(
@@ -240 +249 @@
     const double result =
         compound.parameter_derivative(
-            FunctionModel::arguments_t(1,arg),
+            args,
             3)[0];
     CPPUNIT_ASSERT(