Changeset 0588e9

Timestamp:

Oct 17, 2013, 7:24:06 AM (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:

f58e56

Parents:

675cd6

git-author:

Frederik Heber <heber@…> (09/26/13 16:19:05)

git-committer:

Frederik Heber <heber@…> (10/17/13 07:24:06)

Message:

AnalyseFragmentationResultsAction now also calculates and prints cycle correction.

• we calculate the energy with and without determined cycles per level and print each energy.
• we make use of the ValueMask for this.
• new parameter fragment-prefix allows specifying a file containg the cycle keysets.
• added regression test for FragmentMolecule's DoCyclesFull.

Files:

• src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.cpp (modified) (5 diffs)
• src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.def (modified) (2 diffs)
• tests/Python/AllActions/options.dat (modified) (1 diff)
• tests/regression/Fragmentation/FragmentMolecule/pre/benzene.pdb (added)
• tests/regression/Fragmentation/FragmentMolecule/testsuite-fragmentation-fragment-molecule.at (modified) (1 diff)

src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.cpp

@@ -42 +42 @@
 
 #include <boost/foreach.hpp>
+#include <boost/lambda/lambda.hpp>
 #include <boost/mpl/remove.hpp>
 
+#include <algorithm>
 #include <fstream>
 #include <iostream>
@@ -100 +102 @@
 }
 
+/** Print cycle correction from received results.
+ *
+ * @param results summed up results container
+ */
+void printReceivedCycleResults(
+    const FragmentationShortRangeResults &results)
+{
+  typedef boost::mpl::remove<
+        boost::mpl::remove<MPQCDataEnergyVector_t, MPQCDataFused::energy_eigenvalues>::type,
+      MPQCDataFused::energy_eigenhistogram>::type
+    MPQCDataEnergyVector_noeigenvalues_t;
+  const std::string energyresult =
+      writeTable<MPQCDataEnergyMap_t, MPQCDataEnergyVector_noeigenvalues_t >()(
+          results.Result_Energy_fused, results.getMaxLevel());
+  LOG(2, "DEBUG: Energy table is \n" << energyresult);
+  std::string filename;
+  filename += FRAGMENTPREFIX + std::string("_CycleEnergy.dat");
+  writeToFile(filename, energyresult);
+}
+
 /** Print (short range) energy, forces, and timings from received results.
  *
@@ -324 +346 @@
 }
 
+// this it taken from
+// http://stackoverflow.com/questions/2291802/is-there-a-c-iterator-that-can-iterate-over-a-file-line-by-line
+namespace detail
+{
+  /** Extend the string class by a friend function.
+   *
+   */
+  class Line : public std::string
+  {
+    friend std::istream & operator>>(std::istream & is, Line & line)
+    {
+      return std::getline(is, line);
+    }
+  };
+}
+
+/** Parse the given stream line-by-line, passing each to \a dest.
+ *
+ * \param is stream to parse line-wise
+ * \param dest output iterator
+ */
+template<class OutIt>
+void read_lines(std::istream& is, OutIt dest)
+{
+    typedef std::istream_iterator<detail::Line> InIt;
+    std::copy(InIt(is), InIt(), dest);
+}
+
+
 ActionState::ptr FragmentationAnalyseFragmentationResultsAction::performCall() {
 
-  // if file is given, parse from file into resultscontainer
+  /// if file is given, parse from file into ResultsContainer
   FragmentationResultContainer& container = FragmentationResultContainer::getInstance();
   if (!params.resultsfile.get().empty()) {
@@ -342 +393 @@
   }
 
-  // get data and keysets from ResultsContainer
+  /// get data and keysets from ResultsContainer
   const std::map<JobId_t, MPQCData> &shortrangedata = container.getShortRangeResults();
   const KeySetsContainer &keysets = container.getKeySets();
@@ -354 +405 @@
   }
 
+  /// calculate normal contributions with (if present) cycles coming at their
+  /// respective bond order.
   std::vector<bool> ValueMask(shortrangedata.size(), true);
   FragmentationShortRangeResults shortrangeresults(shortrangedata, keysets, forcekeysets, ValueMask);
   shortrangeresults(shortrangedata);
   printReceivedShortResults(shortrangeresults);
+
+  /// now do we need to calculate the cycle contribution
+  // check whether there are cycles in container or else in file
+  KeySetsContainer cycles = container.getCycles();
+  if (cycles.KeySets.empty()) {
+    // parse from file if cycles is empty
+    boost::filesystem::path filename(
+        params.prefix.get() + std::string(CYCLEKEYSETFILE));
+    if (boost::filesystem::exists(filename)) {
+      LOG(1, "INFO: Parsing cycles file " << filename.string() << ".");
+      // parse file line by line
+      std::ifstream File;
+      File.open(filename.string().c_str());
+      typedef std::istream_iterator<detail::Line> InIt;
+      for (InIt iter = InIt(File); iter != InIt(); ++iter) {
+        KeySetsContainer::IntVector cycle;
+        std::stringstream line(*iter);
+        while (line.good()) {
+          int id;
+          line >> id >> ws;
+          cycle.push_back(id);
+        }
+        if (!cycle.empty()) {
+          LOG(2, "DEBUG: Adding cycle " << cycle << ".");
+          cycles.insert( cycle, cycle.size());
+        }
+      }
+      File.close();
+    } else {
+      LOG(1, "INFO: Cycles file not found at " << filename.string() << ".");
+    }
+  }
+
+  // copy the shortrangeresults such that private MaxLevel is set in
+  // FragmentationShortRangeResults
+  FragmentationShortRangeResults cycleresults(shortrangeresults);
+  // calculate energy if cycles are calculated fully at each level already
+  if (!cycles.KeySets.empty()) {
+    // gather cycle sizes
+    std::vector<size_t> cyclesizes(cycles.KeySets.size());
+    std::transform(
+        cycles.KeySets.begin(), cycles.KeySets.end(),
+        cyclesizes.begin(),
+        boost::bind(&KeySetsContainer::IntVector::size, boost::lambda::_1)
+    );
+    // get maximum
+    std::vector<size_t>::const_iterator maximum_size =
+        std::max_element(cyclesizes.begin(), cyclesizes.end());
+
+    // loop from level 1 to max ring size and gather corrections
+    for (size_t level = 1; level <= *maximum_size; ++level) {
+      // create ValueMask for this level by stepping through each keyset and checking size
+      std::vector<bool> localValueMask(shortrangedata.size(), false);
+      size_t index=0;
+      KeySetsContainer::ArrayOfIntVectors::const_iterator keysetsiter = keysets.KeySets.begin();
+      KeySetsContainer::ArrayOfIntVectors::const_iterator cyclesiter = cycles.KeySets.begin();
+      for (; (keysetsiter != keysets.KeySets.end()) && (cyclesiter != cycles.KeySets.end());) {
+        if (cyclesiter->size() > keysetsiter->size()) {
+          // add if not greater than level in size
+          if ((*keysetsiter).size() <= level)
+            localValueMask[index] = true;
+          ++keysetsiter;
+          ++index;
+        } else if (cyclesiter->size() < keysetsiter->size()) {
+          ++cyclesiter;
+        } else { // both sets have same size
+          if (*cyclesiter > *keysetsiter) {
+            // add if not greater than level in size
+            if ((*keysetsiter).size() <= level)
+              localValueMask[index] = true;
+            ++keysetsiter;
+            ++index;
+          } else if (*cyclesiter < *keysetsiter) {
+            ++cyclesiter;
+          } else {
+            // also always add all cycles
+            localValueMask[index] = true;
+            ++cyclesiter;
+            ++keysetsiter;
+            ++index;
+          }
+        }
+      }
+      // activate rest if desired by level
+      for (; keysetsiter != keysets.KeySets.end(); ++keysetsiter) {
+        if ((*keysetsiter).size() <= level)
+          localValueMask[index] = true;
+        ++index;
+      }
+      LOG(2, "DEBUG: ValueMask for cycle correction at level " << level << " is "
+          << localValueMask << ".");
+      // create FragmentationShortRangeResults
+      FragmentationShortRangeResults localresults(shortrangedata, keysets, forcekeysets, localValueMask);
+      // and perform summation
+      localresults(shortrangedata);
+      // finally, extract the corrections
+      cycleresults.Result_Energy_fused[level-1] =
+          localresults.Result_Energy_fused.back();
+      cycleresults.Result_Time_fused[level-1] =
+          localresults.Result_Time_fused.back();
+      cycleresults.Result_Force_fused[level-1] =
+          localresults.Result_Force_fused.back();
+    }
+    printReceivedCycleResults(cycleresults);
+  }
 
   // adding obtained forces

src/Actions/FragmentationAction/AnalyseFragmentationResultsAction.def

@@ -8 +8 @@
 // all includes and forward declarations necessary for non-integral types below
 #include <boost/filesystem/path.hpp>
+#include <string>
 
 #include "Parameters/Validators/DummyValidator.hpp"
@@ -14 +15 @@
 // 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)(bool)
-#define paramtokens ("fragment-resultfile")("store-grids")
-#define paramdescriptions ("parse fragment results from the given file")("whether to store sampled grids in homology files")
-#define paramdefaults (PARAM_DEFAULT(""))(PARAM_DEFAULT("0"))
-#define paramreferences (resultsfile)(DoStoreGrids)
+#define paramtypes (boost::filesystem::path)(bool)(std::string)
+#define paramtokens ("fragment-resultfile")("store-grids")("fragment-prefix")
+#define paramdescriptions ("parse fragment results from the given file")("whether to store sampled grids in homology files")("prefix of each fragment file")
+#define paramdefaults (PARAM_DEFAULT(""))(PARAM_DEFAULT("0"))(NOPARAM_DEFAULT)
+#define paramreferences (resultsfile)(DoStoreGrids)(prefix)
 #define paramvalids \
 (DummyValidator< boost::filesystem::path >()) \
-(DummyValidator< bool >())
+(DummyValidator< bool >()) \
+(DummyValidator< std::string >())
 
 #undef statetypes

tests/Python/AllActions/options.dat

@@ -82 +82 @@
 fragment-molecule       "./"
 fragment-path   "test/"
+fragment-prefix "BondFragment"
 fragment-resultfile     "results.dat"
 grid-level      "5"

tests/regression/Fragmentation/FragmentMolecule/testsuite-fragmentation-fragment-molecule.at

@@ -34 +34 @@
 
 AT_CLEANUP
+
+AT_SETUP([Fragmentation - Fragmentation with cycles])
+AT_KEYWORDS([fragmentation fragment-molecule cycle])
+
+file=benzene.pdb
+AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Fragmentation/FragmentMolecule/pre/$file $file], 0)
+AT_CHECK([chmod u+w $file], 0, [ignore], [ignore])
+AT_CHECK([../../molecuilder -i $file --update-molecules --select-all-atoms --fragment-molecule ./BondFragment --DoCyclesFull 1 --distance 2 --order 2 --output-types "xyz"], 0, [ignore], [ignore])
+AT_CHECK([ls -l BondFragment??.xyz | wc -l], 0, [13
+], [ignore])
+AT_CHECK([cat <BondFragment12.xyz | wc -l], 0, [14
+], [ignore])
+AT_CHECK([cat <BondFragmentCycles.dat | wc -w], 0, [6
+], [ignore])
+
+AT_CLEANUP