Changeset 06653a for src/Fragmentation/Summation

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:

91e7658

Parents:

cb30d9

git-author:

Frederik Heber <heber@…> (06/08/16 13:22:11)

git-committer:

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

Message:

SamplingGrid has downsample() function.

• we need this when small fragment grids are calculated at a higher resolution than full_sample grid is supposed to be calculated at.
• extended SamplingGridUnitTest to downsample() function, needed special attention on boundary of grid.

Location:

src/Fragmentation/Summation/SetValues

Files:

• SamplingGrid.cpp (modified) (3 diffs)
• SamplingGrid.hpp (modified) (1 diff)
• unittests/SamplingGridUnitTest.cpp (modified) (2 diffs)
• unittests/SamplingGridUnitTest.hpp (modified) (2 diffs)

src/Fragmentation/Summation/SetValues/SamplingGrid.cpp

@@ -43 +43 @@
 #include "Fragmentation/Summation/SetValues/SamplingGrid.hpp"
 
+#include <boost/assign.hpp>
 #include <boost/bind.hpp>
 #include <algorithm>
@@ -610 +611 @@
             "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
         ASSERT( sourceiter != source_sampled_grid.end(),
-            "SamplingGrid::addWindowOntoWindow() - destiter is already at end of window.");
+            "SamplingGrid::addWindowOntoWindow() - sourceiter is already at end of window.");
         op(*destiter, *sourceiter);
         ++destiter;
@@ -658 +659 @@
 }
 
+/** Struct contains a single point with displacements from the
+ * central point and the weight in the restriction.
+ */
+struct PointWeight_t {
+  PointWeight_t(const int &d1, const int &d2, const int &d3, const double &_weight) :
+    displacement(NDIM),
+    weight(_weight)
+  {
+    displacement[0] = d1; displacement[1] = d2; displacement[2] = d3;
+  }
+  typedef std::vector<int>  displacement_t;
+  displacement_t displacement;
+  double weight;
+};
+
+static void getLengthsOfGrid(
+    int _total[NDIM],
+    const SamplingGrid &_grid)
+{
+        const size_t gridpoints_axis = _grid.getGridPointsPerAxis();
+        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
+        for (size_t index=0; index<NDIM; ++index) {
+          if (fabs(_grid.end[index] - _grid.begin[index]) > std::numeric_limits<double>::epsilon()*1e4) {
+            const double delta = (double)gridpoints_axis/(_grid.end[index] - _grid.begin[index]);
+            _total[index] = delta*(_grid.end_window[index] - _grid.begin_window[index])+round_offset;
+          } else
+            _total[index] = 0;
+          ASSERT (_total[index] == ::pow(2, _grid.level),
+              "SamplingGrid::downsample() - total "+toString(_total[index])
+              +" does not match 2^level: "+toString(_grid.level));
+        }
+}
+
+//!> stencil for full weight restriction, see vmg's stencils.hpp
+static const std::vector< PointWeight_t > FullWeightNearestNeighbor =
+    boost::assign::list_of
+    ( PointWeight_t( 0,  0,  0, 0.125) )
+    ( PointWeight_t( 1,  0,  0, 0.0625) )
+    ( PointWeight_t(-1,  0,  0, 0.0625) )
+    ( PointWeight_t( 0,  1,  0, 0.0625) )
+    ( PointWeight_t( 0, -1,  0, 0.0625) )
+    ( PointWeight_t( 0,  0,  1, 0.0625) )
+    ( PointWeight_t( 0,  0, -1, 0.0625) )
+    ( PointWeight_t( 1,  1,  0, 0.03125) )
+    ( PointWeight_t( 1, -1,  0, 0.03125) )
+    ( PointWeight_t(-1,  1,  0, 0.03125) )
+    ( PointWeight_t(-1, -1,  0, 0.03125) )
+    ( PointWeight_t( 0,  1,  1, 0.03125) )
+    ( PointWeight_t( 0,  1, -1, 0.03125) )
+    ( PointWeight_t( 0, -1,  1, 0.03125) )
+    ( PointWeight_t( 0, -1, -1, 0.03125) )
+    ( PointWeight_t( 1,  0,  1, 0.03125) )
+    ( PointWeight_t( 1,  0, -1, 0.03125) )
+    ( PointWeight_t(-1,  0,  1, 0.03125) )
+    ( PointWeight_t(-1,  0, -1, 0.03125) )
+    ( PointWeight_t( 1,  1,  1, 0.015625) )
+    ( PointWeight_t( 1,  1, -1, 0.015625) )
+    ( PointWeight_t( 1, -1,  1, 0.015625) )
+    ( PointWeight_t(-1,  1,  1, 0.015625) )
+    ( PointWeight_t( 1, -1, -1, 0.015625) )
+    ( PointWeight_t(-1,  1, -1, 0.015625) )
+    ( PointWeight_t(-1, -1,  1, 0.015625) )
+    ( PointWeight_t(-1, -1, -1, 0.015625) )
+;
+
+int getValidIndex(
+    const PointWeight_t::displacement_t &_disp,
+    const int N[NDIM],
+    const int length[NDIM])
+{
+  int index = 0;
+  // we simply truncate in case of out of bounds access
+  if ((N[2]+_disp[2] >= 0) && (N[2]+_disp[2] < length[2]))
+    index += _disp[2];
+  if ((N[1]+_disp[1] >= 0) && (N[1]+_disp[1] < length[1]))
+    index += _disp[1]*length[2];
+  if ((N[0]+_disp[0] >= 0) && (N[0]+_disp[0] < length[0]))
+    index += _disp[0]*length[1]*length[2];
+  return index;
+}
+
+void restrictFullWeight(
+    SamplingGrid::sampledvalues_t &_coarse_level,
+    const int length_c[NDIM],
+    const SamplingGrid::sampledvalues_t &_fine_level,
+    const int length_f[NDIM])
+{
+  int N_c[NDIM];
+  int N_f[NDIM];
+  SamplingGrid::sampledvalues_t::iterator coarseiter = _coarse_level.begin();
+  for(N_c[0]=0, N_f[0]=0; (N_c[0] < length_c[0]) && (N_f[0] < length_f[0]); ++N_c[0], N_f[0] +=2) {
+    for(N_c[1]=0, N_f[1]=0; (N_c[1] < length_c[1]) && (N_f[1] < length_f[1]); ++N_c[1], N_f[1] +=2) {
+      for(N_c[2]=0, N_f[2]=0; (N_c[2] < length_c[2]) && (N_f[2] < length_f[2]); ++N_c[2], N_f[2] +=2) {
+        const int index_base = N_f[2] + (N_f[1] + N_f[0]*length_f[1])*length_f[2];
+        // go through stencil and add each point relative to displacement with weight
+        for (std::vector< PointWeight_t >::const_iterator weightiter = FullWeightNearestNeighbor.begin();
+            weightiter != FullWeightNearestNeighbor.end(); ++weightiter) {
+          const PointWeight_t::displacement_t disp = weightiter->displacement;
+          const int index_disp = getValidIndex(disp, N_f, length_f);
+          *coarseiter += _fine_level[index_base+index_disp]*weightiter->weight;
+        }
+        ++coarseiter;
+      }
+      ASSERT ( (N_c[2] == length_c[2]) && (N_f[2] == length_f[2]),
+          "restrictFullWeight() - N_c "+toString(N_c[2])+" != length_c "+toString(length_c[2])
+          +" or N_f "+toString(N_f[2])+" != length_f "+toString(length_f[2]));
+    }
+    ASSERT ( (N_c[1] == length_c[1]) && (N_f[1] == length_f[1]),
+        "restrictFullWeight() - N_c "+toString(N_c[1])+" != length_c "+toString(length_c[1])
+        +" or N_f "+toString(N_f[1])+" != length_f "+toString(length_f[1]));
+  }
+  ASSERT ( (N_c[0] == length_c[0]) && (N_f[0] == length_f[0]),
+      "restrictFullWeight() - N_c "+toString(N_c[0])+" != length_c "+toString(length_c[0])
+      +" or N_f "+toString(N_f[0])+" != length_f "+toString(length_f[0]));
+  ASSERT( coarseiter ==  _coarse_level.end(),
+      "restrictFullWeight() - coarseiter is not at end of values.");
+}
+
+void SamplingGrid::downsample(
+    SamplingGrid& instance,
+    const SamplingGrid& other,
+    const int _level)
+{
+  if (&instance != &other) {
+    // take over properties
+    static_cast<SamplingGridProperties &>(instance) = other;
+    instance.setWindowSize(other.begin_window, other.end_window);
+    if (_level >= other.level) {
+      instance.sampled_grid = other.sampled_grid;
+    } else {
+      // if desired level is smaller we need to downsample
+      // we do this similarly to vmg::RestrictionFullWeight (i.e. a full nearest
+      // neighbor interpolation) and always one grid level at a time till we
+      // have reached the desired one
+
+      // the reference such that we never have to copy the full grid but only
+      // downsampled ones
+      const sampledvalues_t * sourcevalues = &other.sampled_grid;
+      int length_d[3];
+      int length_s[3];
+      getLengthsOfGrid(length_s, other);
+      for (instance.level = other.level-1; instance.level >= _level; --instance.level) {
+        getLengthsOfGrid(length_d, instance);
+        // we always have an eighth of the number of sample points as we stop
+        sampledvalues_t downsampled(sourcevalues->size()/(size_t)8, 0.);
+        restrictFullWeight(downsampled, length_d, *sourcevalues, length_s);
+        // then copy the downsampled values
+        instance.sampled_grid = downsampled;
+        sourcevalues = &instance.sampled_grid;
+        // and exchange lengths
+        for (size_t i=0;i<3;++i) {
+          length_s[i] = length_d[i];
+        }
+      }
+      instance.level = _level;
+
+      // and finally, renormalize downsampled grid to old value
+//      instance *= other.integral()/instance.integral();
+    }
+  }
+}
+
 std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other)
 {

src/Fragmentation/Summation/SetValues/SamplingGrid.hpp

@@ -166 +166 @@
   }
 
+  /** Sample a given grid \a other down to grid level \a _level and store
+   * in \a instance.
+   *
+   * \param instance given instance to store downsampled grid in
+   * \param other instance to get grid from
+   * \param _level level to sample down to
+   */
+  static void downsample(SamplingGrid& instance, const SamplingGrid& other, const int _level);
+
   /** Returns the numeric integral over the grid.
    *

src/Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.cpp

@@ -49 +49 @@
 
 #include <boost/assign.hpp>
+
 #include <cmath>
 #include <numeric>
@@ -549 +550 @@
 }
 
+#ifdef HAVE_INLINE
+inline
+#endif
+static int calculateIndex(const int N[NDIM], const int &_length)
+{
+  return N[2] + N[1]*_length + N[0]*_length*_length;
+}
+
+#ifdef HAVE_INLINE
+inline
+#endif
+static double calculateDistanceSquared(const int N[NDIM], const int &_length)
+{
+  return
+      ::pow(N[0]/(double)_length-.5,2)
+      + ::pow(N[1]/(double)_length-.5,2)
+      + ::pow(N[2]/(double)_length-.5,2);
+}
+/** UnitTest for downsample()
+ */
+void SamplingGridTest::downsampleTest()
+{
+  const double begin[NDIM] = { 0., 0., 0. };
+  const double end[NDIM] = { 1., 1., 1. };
+
+  // simple test, one level difference, same value everywhere
+  {
+    SamplingGrid::sampledvalues_t checkvalues;
+    SamplingGrid::sampledvalues_t othervalues;
+    const double othergrid_value = 1.5;
+    for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
+      checkvalues += othergrid_value;
+    for (size_t i=0; i< NUMBEROFSAMPLES(3); ++i)
+      othervalues += othergrid_value;
+
+    SamplingGrid largegrid(begin, end, 3, othervalues);
+//    std::cout << " largegrid value " << largegrid.sampled_grid << std::endl;
+    SamplingGrid smallgrid(begin, end, 2);
+    SamplingGrid::downsample(smallgrid, largegrid, 2);
+    SamplingGrid checkgrid(begin, end, 2, checkvalues);
+//    std::cout << " smallgrid value " << smallgrid.sampled_grid << std::endl;
+//    std::cout << " checkgrid value " << checkgrid.sampled_grid << std::endl;
+    CPPUNIT_ASSERT_EQUAL( smallgrid, checkgrid );
+  }
+
+  // simple test, two level difference, same value everywhere
+  {
+    SamplingGrid::sampledvalues_t checkvalues;
+    SamplingGrid::sampledvalues_t othervalues;
+    const double othergrid_value = 1.5;
+    for (size_t i=0; i< NUMBEROFSAMPLES(2); ++i)
+      checkvalues += othergrid_value;
+    for (size_t i=0; i< NUMBEROFSAMPLES(4); ++i)
+      othervalues += othergrid_value;
+
+    SamplingGrid largegrid(begin, end, 4, othervalues);
+//    std::cout << " largegrid value " << largegrid.sampled_grid << std::endl;
+    SamplingGrid smallgrid(begin, end, 2);
+    SamplingGrid::downsample(smallgrid, largegrid, 2);
+    SamplingGrid checkgrid(begin, end, 2, checkvalues);
+//    std::cout << " smallgrid value " << smallgrid.sampled_grid << std::endl;
+//    std::cout << " checkgrid value " << checkgrid.sampled_grid << std::endl;
+    CPPUNIT_ASSERT_EQUAL( smallgrid, checkgrid );
+  }
+
+  // same grid, window equals grids, different values
+  {
+    const int length_c = ::pow(2,2);
+    const int length_f = ::pow(2,3);
+    SamplingGrid::sampledvalues_t checkvalues((int)::pow(length_c, NDIM), 0.);
+    SamplingGrid::sampledvalues_t othervalues((int)::pow(length_f, NDIM), 0.);
+    int N[NDIM];
+    for (N[0]=0; N[0]< length_f; ++N[0]) {
+      for (N[1]=0; N[1]< length_f; ++N[1]) {
+        for (N[2]=0; N[2]< length_f; ++N[2]) {
+          const int index = calculateIndex(N, length_f);
+          const double dist = calculateDistanceSquared(N, length_f);
+          othervalues[index] = cos(M_PI*dist/1.5);
+        }
+      }
+    }
+    for (N[0]=0; N[0]< length_c; ++N[0]) {
+      for (N[1]=0; N[1]< length_c; ++N[1]) {
+        for (N[2]=0; N[2]< length_c; ++N[2]) {
+          const int index = calculateIndex(N, length_c);
+          const double dist = calculateDistanceSquared(N, length_c);
+          checkvalues[index] = cos(M_PI*dist/1.5);
+        }
+      }
+    }
+
+    SamplingGrid largegrid(begin, end, 3, othervalues);
+//    std::cout << " largegrid value " << largegrid.sampled_grid << std::endl;
+    SamplingGrid smallgrid(begin, end, 2);
+    SamplingGrid::downsample(smallgrid, largegrid, 2);
+    SamplingGrid checkgrid(begin, end, 2, checkvalues);
+//    std::cout << " smallgrid value " << smallgrid.sampled_grid << std::endl;
+//    std::cout << " checkgrid value " << checkgrid.sampled_grid << std::endl;
+    static const double threshold = 1.2e-1;
+    for (N[0]=0; N[0]< length_c; ++N[0])
+      for (N[1]=0; N[1]< length_c; ++N[1])
+        for (N[2]=0; N[2]< length_c; ++N[2]) {
+          const double check_threshold =
+                    (((N[0] == 0) || (N[0] == length_c-1))
+                  && ((N[1] == 0) || (N[1] == length_c-1))
+                  && ((N[2] == 0) || (N[2] == length_c-1)))
+              ? 2.*threshold : threshold;
+          const int index = calculateIndex(N, length_c);
+//          std::cout << "Comparing "
+//              << fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index]))
+//              << " < " << threshold << std::endl;
+          CPPUNIT_ASSERT(
+              fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index])) < check_threshold);
+        }
+  }
+}

src/Fragmentation/Summation/SetValues/unittests/SamplingGridUnitTest.hpp

@@ -38 +38 @@
     CPPUNIT_TEST ( equality_Test );
     CPPUNIT_TEST ( serializeTest );
+    CPPUNIT_TEST ( downsampleTest );
     CPPUNIT_TEST_SUITE_END();
 
@@ -57 +58 @@
       void equality_Test();
       void serializeTest();
+      void downsampleTest();
 
 private: