Changeset 7c1091 for src/Potentials/CompoundPotential.cpp

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.

File:

• src/Potentials/CompoundPotential.cpp (modified) (11 diffs)

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;
 }