Changeset bb8b0a for src/Actions/PotentialAction/FitPartialChargesAction.cpp

Timestamp:

Mar 20, 2017, 3:03:17 PM (8 years ago)

Author:

Frederik Heber <heber@…>

Branches:

FitPartialCharges_GlobalError, PartialCharges_OrthogonalSummation

Children:

3607c3

Parents:

5b3781

git-author:

Frederik Heber <heber@…> (06/12/16 19:02:32)

git-committer:

Frederik Heber <heber@…> (03/20/17 15:03:17)

Message:

First version of orthogonal summation of partial charges.

• implemented directly in FitPartialChargesAction::performCall().

File:

• src/Actions/PotentialAction/FitPartialChargesAction.cpp (modified) (5 diffs)

src/Actions/PotentialAction/FitPartialChargesAction.cpp

@@ -49 +49 @@
 #include <boost/filesystem.hpp>
 #include <boost/foreach.hpp>
+#include <boost/mpl/for_each.hpp>
+
 #include <algorithm>
 #include <functional>
@@ -65 +67 @@
 #include "Fragmentation/Homology/HomologyContainer.hpp"
 #include "Fragmentation/Homology/HomologyGraph.hpp"
+#include "Fragmentation/Summation/Containers/PartialChargesFused.hpp"
+#include "Fragmentation/Summation/Containers/PartialChargesMap.hpp"
 #include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
+#include "Fragmentation/Summation/ZeroInstanceInitializer.hpp"
 #include "FunctionApproximation/Extractors.hpp"
 #include "Potentials/PartialNucleiChargeFitter.hpp"
@@ -400 +405 @@
 }
 
+struct KeySetSizeComp {
+  bool operator() (const KeySet &a, const KeySet &b) { return a.size()<b.size(); }
+} keyset_comparator;
+
+#include <boost/fusion/sequence.hpp>
+#include <boost/mpl/for_each.hpp>
+
+#include "Fragmentation/Summation/AllLevelOrthogonalSummator.hpp"
+#include "Fragmentation/Summation/IndexSet.hpp"
+#include "Fragmentation/Summation/IndexSetContainer.hpp"
+#include "Fragmentation/Summation/SubsetMap.hpp"
+#include "Fragmentation/Summation/Containers/PartialChargesFused.hpp"
+#include "Fragmentation/Summation/Containers/PartialChargesMap.hpp"
+#include "Fragmentation/Summation/SetValues/IndexedPartialCharges.hpp"
+
 ActionState::ptr PotentialFitPartialChargesAction::performCall()
 {
@@ -409 +429 @@
     return Action::failure;
   }
+  LOG(3, "There are  " << selected_atoms.size() << " selected atoms.");
 
   /// obtain possible fragments to each selected atom
@@ -442 +463 @@
   /// obtain average charge for each atom the fitted charges over all its fragments
   detail::fitted_charges_t fitted_charges;
-  for (std::vector<const atom *>::const_iterator atomiter = selected_atoms.begin();
-      atomiter != selected_atoms.end(); ++atomiter) {
-    const atomId_t walkerid = (*atomiter)->getId();
-    const double average_charge = fitAverageChargeToAtom(
-        *atomiter, atomfragments, keyset_graphs, fittedcharges_per_fragment);
-
-    if (average_charge != 0.) {
-      LOG(2, "DEBUG: For atom " << **atomiter << " we have an average charge of "
-          << average_charge);
-
-      fitted_charges.insert( std::make_pair(walkerid, average_charge) );
-    }
+//  for (std::vector<const atom *>::const_iterator atomiter = selected_atoms.begin();
+//      atomiter != selected_atoms.end(); ++atomiter) {
+//    const atomId_t walkerid = (*atomiter)->getId();
+//    const double average_charge = fitAverageChargeToAtom(
+//        *atomiter, atomfragments, keyset_graphs, fittedcharges_per_fragment);
+//
+//    if (average_charge != 0.) {
+//      LOG(2, "DEBUG: For atom " << **atomiter << " we have an average charge of "
+//          << average_charge);
+//
+//      fitted_charges.insert( std::make_pair(walkerid, average_charge) );
+//    }
+//  }
+  {
+    // TODO:: convert the data into boost::fusion::map format
+    std::map<JobId_t, PartialChargesMap_t> PartialCharges_fused;
+    std::vector<IndexSet> indexsets;
+    {
+      JobId_t counter = 0;
+      for (std::set<KeySet>::const_iterator iter = fragments.begin();
+          iter != fragments.end(); ++iter) {
+        const KeySet &currentset = *iter;
+        // place as IndexSet into container
+        {
+          IndexSet tempset;
+          tempset.insert(currentset.begin(), currentset.end());
+          indexsets.push_back(tempset);
+        }
+        const detail::KeysetsToGraph_t::const_iterator graphiter = keyset_graphs.find(currentset);
+        if (graphiter == keyset_graphs.end()) {
+          ELOG(1, "Could not find graph to keyset " << currentset);
+          return Action::failure;
+        }
+        const HomologyGraph &currentgraph = graphiter->second;
+        const detail::GraphFittedChargeMap_t::const_iterator fragmentiter =
+            fittedcharges_per_fragment.find(currentgraph);
+        if (fragmentiter == fittedcharges_per_fragment.end()) {
+          ELOG(1, "Could not find fragment to graph " << currentgraph);
+          return Action::failure;
+        }
+        const PartialNucleiChargeFitter::charges_t &charges = fragmentiter->second;
+        PartialChargesMap_t chargemap;
+        const AtomFragmentsMap::indices_t &full_currentset = atomfragments.getFullKeyset(currentset);
+        IndexedPartialCharges::indices_t indices(full_currentset.begin(), full_currentset.end());
+        IndexedPartialCharges::values_t values(charges.begin(), charges.end());
+        LOG(3, "Inserting " << indices << " with charges " << charges << " into chargemap instance.");
+        boost::fusion::at_key<PartialChargesFused::partial_charges_t>(chargemap) =
+            IndexedPartialCharges(indices, values);
+        PartialCharges_fused.insert( std::make_pair(counter++, chargemap) );
+      }
+    }
+    ASSERT( indexsets.size() == fragments.size(),
+        "PotentialFitPartialChargesAction::performCall() - not all fragments' keysets were created?");
+
+    // prepare index set hierarchy
+    std::map< JobId_t, size_t > IdentityLookup;
+    size_t MaxLevel = std::max_element(fragments.begin(), fragments.end(), keyset_comparator)->size();
+    LOG(3, "The maximum level is " << MaxLevel);
+    size_t max_indices = fragments.size();
+    LOG(3, "There are  " << max_indices << " for the selected atoms.");
+    {
+      for (size_t index = 0; index < max_indices; ++index)
+        IdentityLookup.insert( std::make_pair( (JobId_t)index, index ) );
+    }
+    IndexSetContainer::ptr container(new IndexSetContainer(indexsets));
+    SubsetMap::ptr subsetmap(new SubsetMap(*container));
+
+    // and sum up
+    PartialChargesMap_t ZeroInstances;
+    ZeroInstanceInitializer<PartialChargesMap_t> initZeroInstance(ZeroInstances);
+    boost::mpl::for_each<PartialChargesVector_t>(boost::ref(initZeroInstance));
+    //!> results per level of summed up partial charges
+    std::vector<PartialChargesMap_t> Result_PartialCharges_fused(MaxLevel);
+    //!> results per index set in terms of value and contribution
+    std::map<
+        IndexSet::ptr,
+        std::pair<PartialChargesMap_t, PartialChargesMap_t> > Result_perIndexSet_PartialCharges;
+    AllLevelOrthogonalSummator<PartialChargesMap_t> partialchargeSummer(
+                subsetmap,
+                PartialCharges_fused,
+                container->getContainer(),
+                IdentityLookup,
+                Result_PartialCharges_fused,
+                Result_perIndexSet_PartialCharges,
+                ZeroInstances);
+    boost::mpl::for_each<PartialChargesVector_t>(boost::ref(partialchargeSummer));
+
+    // TODO: place results into fitted_charges
+    const IndexedPartialCharges::indexedvalues_t indexed_partial_charges =
+        boost::fusion::at_key<PartialChargesFused::partial_charges_t>(
+            Result_PartialCharges_fused.back()
+            ).getValues();
+    for (IndexedPartialCharges::indexedvalues_t::const_iterator iter = indexed_partial_charges.begin();
+        iter != indexed_partial_charges.end(); ++iter)
+    fitted_charges.insert( std::make_pair( iter->first, iter->second.charge));
+    LOG(3, "Summation has brought forth the following charges per atom index: " << fitted_charges);
   }