Changeset 57fda0

Timestamp:

Sep 14, 2016, 6:42:53 PM (8 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_StructOpt_integration_tests, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, 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_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, GeometryObjects, Gui_displays_atomic_force_velocity, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, RotateToPrincipalAxisSystem_UndoRedo, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, Ubuntu_1604_changes, stable

Children:

b839fb

Parents:

6a5eb2

git-author:

Frederik Heber <heber@…> (06/06/16 14:46:53)

git-committer:

Frederik Heber <heber@…> (09/14/16 18:42:53)

Message:

SamplingGrid::operator*=(), ::superposeOtherGrids() and ::integral() now allow non-equivalent grids.

File:

• src/Fragmentation/Summation/SetValues/SamplingGrid.cpp (modified) (4 diffs)

src/Fragmentation/Summation/SetValues/SamplingGrid.cpp

@@ -45 +45 @@
 #include <boost/assign.hpp>
 #include <boost/bind.hpp>
+#include <boost/shared_ptr.hpp>
+
 #include <algorithm>
 #include <limits>
@@ -158 +160 @@
 }
 
+static void getDownsampledGrid(
+    const SamplingGrid &_reference_grid,
+    const SamplingGrid &_grid,
+    boost::shared_ptr<SamplingGrid> &_weight_downsampled)
+{
+  static const double round_offset(
+      (std::numeric_limits<size_t>::round_style == std::round_toward_zero) ?
+          0.5 : 0.); // need offset to get to round_toward_nearest behavior
+  const int surplus_level = _reference_grid.getSurplusLevel(_grid)+round_offset;
+  // need to downsample first
+  _weight_downsampled.reset(new SamplingGrid); // may use default cstor as downsample does all settings
+  SamplingGrid::downsample(*_weight_downsampled, _grid, _reference_grid.level-surplus_level);
+}
+
 SamplingGrid& SamplingGrid::operator*=(const SamplingGrid& other)
 {
   // check that grids are compatible
-  if (isEquivalent(other)) {
-    /// get minimum of window
-    double min_begin_window[NDIM];
-    double min_end_window[NDIM];
-    bool doShrink = false;
-    for (size_t index=0; index<NDIM;++index) {
-      if (begin_window[index] <= other.begin_window[index]) {
-        min_begin_window[index] = other.begin_window[index];
-        doShrink = true;
-      } else {
-        min_begin_window[index] = begin_window[index];
-      }
-      if (end_window[index] <= other.end_window[index]) {
-        min_end_window[index] = end_window[index];
-      } else {
-        min_end_window[index] = other.end_window[index];
-        doShrink = true;
-      }
-    }
-    LOG(4, "DEBUG: min begin is " << min_begin_window[0] << "," << min_begin_window[1] << "," << min_begin_window[2] << ".");
-    LOG(4, "DEBUG: min end is " << min_end_window[0] << "," << min_end_window[1] << "," << min_end_window[2] << ".");
-    if (doShrink)
-      shrinkWindow(min_begin_window, min_end_window);
-    addWindowOntoWindow(
-        other.begin_window,
-        other.end_window,
-        begin_window,
-        end_window,
-        sampled_grid,
-        other.sampled_grid,
-        boost::bind(multiplyElements, _1, _2, 1.),
-        sourcewindow);
-  } else {
-    ASSERT(0, "SamplingGrid::operator*=() - multiplying incongruent grids is so far not in the cards.");
-  }
+  ASSERT(isCompatible(other),
+      "SamplingGrid::operator*=() - multiplying incomatible grids is so far not in the cards.");
+  const SamplingGrid *other_grid = &other;
+  boost::shared_ptr<SamplingGrid> other_downsampled;
+  if (!isEquivalent(other)) {
+    getDownsampledGrid(*this, other, other_downsampled);
+    other_grid = other_downsampled.get();
+  }
+  /// get minimum of window
+  double min_begin_window[NDIM];
+  double min_end_window[NDIM];
+  bool doShrink = false;
+  for (size_t index=0; index<NDIM;++index) {
+    if (begin_window[index] <= other.begin_window[index]) {
+      min_begin_window[index] = other.begin_window[index];
+      doShrink = true;
+    } else {
+      min_begin_window[index] = begin_window[index];
+    }
+    if (end_window[index] <= other.end_window[index]) {
+      min_end_window[index] = end_window[index];
+    } else {
+      min_end_window[index] = other.end_window[index];
+      doShrink = true;
+    }
+  }
+  LOG(4, "DEBUG: min begin is " << min_begin_window[0] << "," << min_begin_window[1] << "," << min_begin_window[2] << ".");
+  LOG(4, "DEBUG: min end is " << min_end_window[0] << "," << min_end_window[1] << "," << min_end_window[2] << ".");
+  if (doShrink)
+    shrinkWindow(min_begin_window, min_end_window);
+  addWindowOntoWindow(
+      other_grid->begin_window,
+      other_grid->end_window,
+      begin_window,
+      end_window,
+      sampled_grid,
+      other_grid->sampled_grid,
+      boost::bind(multiplyElements, _1, _2, 1.),
+      sourcewindow);
+
   return *this;
 }
@@ -201 +222 @@
 void SamplingGrid::superposeOtherGrids(const SamplingGrid &other, const double prefactor)
 {
-  /// check that grids are compatible
-  if (isEquivalent(other)) {
-    /// get maximum of window
-    double max_begin_window[NDIM];
-    double max_end_window[NDIM];
-    bool doExtend = false;
-    for (size_t index=0; index<NDIM;++index) {
-      if (begin_window[index] >= other.begin_window[index]) {
-        max_begin_window[index] = other.begin_window[index];
-        doExtend = true;
-      } else {
-        max_begin_window[index] = begin_window[index];
-      }
-      if (end_window[index] >= other.end_window[index]) {
-        max_end_window[index] = end_window[index];
-      } else {
-        max_end_window[index] = other.end_window[index];
-        doExtend = true;
-      }
-    }
-    LOG(4, "DEBUG: max begin is " << max_begin_window[0] << "," << max_begin_window[1] << "," << max_begin_window[2] << ".");
-    LOG(4, "DEBUG: max end is " << max_end_window[0] << "," << max_end_window[1] << "," << max_end_window[2] << ".");
-    if (doExtend)
-      extendWindow(max_begin_window, max_end_window);
-    /// and copy other into larger window, too
-    addOntoWindow(other.begin_window, other.end_window, other.sampled_grid, prefactor);
-  } else {
-    ASSERT(0, "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards.");
-  }
+  // check that grids are compatible
+  ASSERT(isCompatible(other),
+      "SamplingGrid::superposeOtherGrids() - superposing incompatible grids is so far not in the cards.");
+  const SamplingGrid *other_grid = &other;
+  boost::shared_ptr<SamplingGrid> other_downsampled;
+  if (!isEquivalent(other)) {
+    getDownsampledGrid(*this, other, other_downsampled);
+    other_grid = other_downsampled.get();
+  }
+  /// get maximum of window
+  double max_begin_window[NDIM];
+  double max_end_window[NDIM];
+  bool doExtend = false;
+  for (size_t index=0; index<NDIM;++index) {
+    if (begin_window[index] >= other.begin_window[index]) {
+      max_begin_window[index] = other.begin_window[index];
+      doExtend = true;
+    } else {
+      max_begin_window[index] = begin_window[index];
+    }
+    if (end_window[index] >= other.end_window[index]) {
+      max_end_window[index] = end_window[index];
+    } else {
+      max_end_window[index] = other.end_window[index];
+      doExtend = true;
+    }
+  }
+  LOG(4, "DEBUG: max begin is " << max_begin_window[0] << "," << max_begin_window[1] << "," << max_begin_window[2] << ".");
+  LOG(4, "DEBUG: max end is " << max_end_window[0] << "," << max_end_window[1] << "," << max_end_window[2] << ".");
+  if (doExtend)
+    extendWindow(max_begin_window, max_end_window);
+  /// and copy other into larger window, too
+  addOntoWindow(other_grid->begin_window, other_grid->end_window, other_grid->sampled_grid, prefactor);
 }
 
@@ -262 +287 @@
 double SamplingGrid::integral(const SamplingGrid &weight) const
 {
-  if (isEquivalent(weight)) {
-    const double volume_element = getVolume()/(double)getTotalGridPoints();
-    double int_value = 0.;
-    sampledvalues_t::const_iterator iter = sampled_grid.begin();
-    sampledvalues_t::const_iterator weightiter = weight.sampled_grid.begin();
-    for (;iter != sampled_grid.end();++iter,++weightiter)
-      int_value += (*weightiter) * (*iter);
-    int_value *= volume_element;
-    //LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
-    return int_value;
-  } else
-    return 0.;
-}
+  // check that grids are compatible
+  ASSERT(isCompatible(weight),
+      "SamplingGrid::integral() - integrating with weights from incompatible grids is so far not in the cards.");
+  const SamplingGrid *weight_grid = &weight;
+  boost::shared_ptr<SamplingGrid> weight_downsampled;
+  if (!isEquivalent(weight)) {
+    getDownsampledGrid(*this, weight, weight_downsampled);
+    weight_grid = weight_downsampled.get();
+  }
+  const double volume_element = getVolume()/(double)getTotalGridPoints();
+  double int_value = 0.;
+  sampledvalues_t::const_iterator iter = sampled_grid.begin();
+  sampledvalues_t::const_iterator weightiter = weight_grid->sampled_grid.begin();
+  for (;iter != sampled_grid.end();++iter,++weightiter)
+    int_value += (*weightiter) * (*iter);
+  int_value *= volume_element;
+  //LOG(2, "DEBUG: SamplingGrid::integral() is " << scientific << setprecision(13) << int_value << ".");
+  return int_value;
+}
+
 void SamplingGrid::setWindowSize(
     const double _begin_window[NDIM],