Changeset 91e7658 for src/Fragmentation

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:

d56e21

Parents:

06653a

git-author:

Frederik Heber <heber@…> (05/26/16 09:46:06)

git-committer:

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

Message:

SamplingGridUnitTest checks downsample on window.

• FIX: getLengthsOfGrid() asserted that total equals 2^{level which is only true if the window spans the whole grid. We now assert for equal or less.}
• Renamed getLengthsOfGrid -> getLengthsOfWindow.
• Extended downsample() tests to check for level difference of one to three.

Location:

src/Fragmentation/Summation/SetValues

Files:

• SamplingGrid.cpp (modified) (4 diffs)
• unittests/SamplingGridUnitTest.cpp (modified) (4 diffs)
• unittests/SamplingGridUnitTest.hpp (modified) (2 diffs)

src/Fragmentation/Summation/SetValues/SamplingGrid.cpp

@@ -674 +674 @@
 };
 
-static void getLengthsOfGrid(
+static void getLengthsOfWindow(
     int _total[NDIM],
     const SamplingGrid &_grid)
@@ -688 +688 @@
           } else
             _total[index] = 0;
-          ASSERT (_total[index] == ::pow(2, _grid.level),
+          // we can only assert that its atmost the maximum number of grid points
+          ASSERT (_total[index] <= ::pow(2, _grid.level),
               "SamplingGrid::downsample() - total "+toString(_total[index])
-              +" does not match 2^level: "+toString(_grid.level));
+              +" is not equal or less than 2^level: "+toString(_grid.level));
         }
 }
@@ -788 +789 @@
     static_cast<SamplingGridProperties &>(instance) = other;
     instance.setWindowSize(other.begin_window, other.end_window);
-    if (_level >= other.level) {
+    ASSERT( _level <= other.level,
+        "SamplingGrid::downsample() - desired level "+toString(_level)
+        +" larger than level "+toString(other.level)+" of the given values.");
+    if (_level == other.level) {
       instance.sampled_grid = other.sampled_grid;
     } else {
@@ -801 +805 @@
       int length_d[3];
       int length_s[3];
-      getLengthsOfGrid(length_s, other);
+      getLengthsOfWindow(length_s, other);
       for (instance.level = other.level-1; instance.level >= _level; --instance.level) {
-        getLengthsOfGrid(length_d, instance);
+        getLengthsOfWindow(length_d, instance);
         // we always have an eighth of the number of sample points as we stop
+        ASSERT( sourcevalues->size() % 8 == 0,
+            "SamplingGrid::downsample() - at level "+toString( instance.level)
+            +" given grid points "+toString(sourcevalues->size())+" are not even numbered per axis anymore.");
         sampledvalues_t downsampled(sourcevalues->size()/(size_t)8, 0.);
         restrictFullWeight(downsampled, length_d, *sourcevalues, length_s);

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

@@ -47 +47 @@
 
 #include "CodePatterns/Assert.hpp"
+#include "CodePatterns/Log.hpp"
 
 #include <boost/assign.hpp>
@@ -568 +569 @@
       + ::pow(N[2]/(double)_length-.5,2);
 }
+
+#ifdef HAVE_INLINE
+inline
+#endif
+static double getBoundaryCaseFactor(const int N[NDIM], const int &_length)
+{
+  static double third_two=::pow(4., 1./3.);
+  double returnthreshold = 1.;
+  for (size_t i=0;i<NDIM;++i)
+    if ((N[i] == 0) || (N[i] == _length-1))
+      returnthreshold *= third_two;
+  return returnthreshold;
+}
+
 /** UnitTest for downsample()
  */
-void SamplingGridTest::downsampleTest()
+void SamplingGridTest::downsample_gridTest()
 {
   const double begin[NDIM] = { 0., 0., 0. };
@@ -615 +630 @@
   }
 
-  // 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.);
+  // same grid, window equals grids, ever larger largegrids
+  const int base_level = 2;
+  for (int grid_difference = 1; grid_difference <= 3; ++grid_difference) {
+//    LOG(2, "Checking with difference " << grid_difference);
+    const int length_c = ::pow(2, base_level);
+    const int length_f = ::pow(2, base_level+grid_difference);
+    SamplingGrid::sampledvalues_t checkvalues((int)::pow(length_c, (int)NDIM), 0.);
+    SamplingGrid::sampledvalues_t othervalues((int)::pow(length_f, (int)NDIM), 0.);
     int N[NDIM];
     for (N[0]=0; N[0]< length_f; ++N[0]) {
@@ -641 +658 @@
     }
 
-    SamplingGrid largegrid(begin, end, 3, othervalues);
+    SamplingGrid largegrid(begin, end, base_level+grid_difference, 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);
+    SamplingGrid smallgrid(begin, end, base_level);
+    SamplingGrid::downsample(smallgrid, largegrid, base_level);
+    SamplingGrid checkgrid(begin, end, base_level, 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])
+    const double threshold = 3.2e-2;
+    // we skip the first value as checkvalues[0] is always zero, hard to get from downsampling
+    for (N[0]=1; 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;
+              threshold*(1.+grid_difference/4.)*getBoundaryCaseFactor(N, length_c);
           const int index = calculateIndex(N, length_c);
 //          std::cout << "Comparing "
 //              << fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index]))
-//              << " < " << threshold << std::endl;
+//              << " < " << check_threshold << std::endl;
           CPPUNIT_ASSERT(
               fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index])) < check_threshold);
+//          if (fabs(checkgrid.sampled_grid[index]) > 1e-10) {
+//            std::cout << "Comparing "
+//                << fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index])/checkgrid.sampled_grid[index])
+//                << " < " << check_threshold << std::endl;
+//            CPPUNIT_ASSERT(
+//                fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index])/checkgrid.sampled_grid[index]) < check_threshold);
+//          }
         }
   }
 }
+
+/** UnitTest for downsample()
+ */
+void SamplingGridTest::downsample_smallerwindowTest()
+{
+  const double begin[NDIM] = { 0., 0., 0. };
+  const double end[NDIM] = { 2., 2., 2. };
+
+  // same grid, window is half of grid, different level by one
+  const int base_level = 3;
+  for (int grid_difference = 1; grid_difference <= 3; ++grid_difference) {
+//    LOG(2, "Checking with difference " << grid_difference);
+    const int window_length_c = ::pow(2, base_level-1);
+    const int window_length_f = ::pow(2, base_level+grid_difference-1);
+    SamplingGrid::sampledvalues_t checkvalues((int)::pow(window_length_c, (int)NDIM), 0.);
+    SamplingGrid::sampledvalues_t othervalues((int)::pow(window_length_f, (int)NDIM), 0.);
+    int N[NDIM];
+    for (N[0]=0; N[0]< window_length_f; ++N[0]) {
+      for (N[1]=0; N[1]< window_length_f; ++N[1]) {
+        for (N[2]=0; N[2]< window_length_f; ++N[2]) {
+          const int index = calculateIndex(N, window_length_f);
+          const double dist = calculateDistanceSquared(N, window_length_f);
+          othervalues[index] = cos(M_PI*dist/1.5);
+        }
+      }
+    }
+    for (N[0]=0; N[0]< window_length_c; ++N[0]) {
+      for (N[1]=0; N[1]< window_length_c; ++N[1]) {
+        for (N[2]=0; N[2]< window_length_c; ++N[2]) {
+          const int index = calculateIndex(N, window_length_c);
+          const double dist = calculateDistanceSquared(N, window_length_c);
+          checkvalues[index] = cos(M_PI*dist/1.5);
+        }
+      }
+    }
+
+    /** Here, we must initialize the larger grid with zero window first and
+     * then add the values into the smaller window.
+     */
+    SamplingGrid largegrid(begin, end, base_level+grid_difference);
+    const double window_begin[NDIM] = {
+        largegrid.getNearestHigherGridPoint(0.5, 0),
+        largegrid.getNearestHigherGridPoint(0.5, 1),
+        largegrid.getNearestHigherGridPoint(0.5, 2) };
+    const double window_end[NDIM] = {
+        largegrid.getNearestLowerGridPoint(1.5, 0),
+        largegrid.getNearestLowerGridPoint(1.5, 1),
+        largegrid.getNearestLowerGridPoint(1.5, 2) };
+    largegrid.setWindow(window_begin, window_end);
+    largegrid.addOntoWindow(window_begin, window_end, othervalues, +1.);
+//    std::cout << " largegrid value " << largegrid.sampled_grid << std::endl;
+
+    // smallgrid is downsample from full large grid
+    SamplingGrid smallgrid(begin, end, base_level);
+    SamplingGrid::downsample(smallgrid, largegrid, base_level);
+//    std::cout << " smallgrid value " << smallgrid.sampled_grid << std::endl;
+
+    // checkgrid is created in the same way as the large grid (only with smaller level)
+    SamplingGrid checkgrid(begin, end, base_level);
+    checkgrid.setWindow(window_begin, window_end);
+    checkgrid.addOntoWindow(window_begin, window_end, checkvalues, +1.);
+//    std::cout << " checkgrid value " << checkgrid.sampled_grid << std::endl;
+
+    // then we compare over the full length
+    /** Note that the threshold does not get better with increasing grid_difference,
+     * on the contrary! For a grid_difference of 1, the finer grid is exact as it
+     * directly sampled from the function. For a larger grid_difference, the finer
+     * grid (being one level away from the desired coarse grid) is just an
+     * approximation although obtained from a finer sampling. Hence, the larger the
+     * grid_difference, the worse is the smallgrid with respect to checkgrid.
+     */
+    const double threshold = 3.2e-2;
+    // we skip the first value as checkvalues[0] is always zero, hard to get from downsampling
+    for (N[0]=1; N[0]< window_length_c; ++N[0])
+      for (N[1]=0; N[1]< window_length_c; ++N[1])
+        for (N[2]=0; N[2]< window_length_c; ++N[2]) {
+          const double check_threshold =
+              threshold*(1.+grid_difference/4.)*getBoundaryCaseFactor(N, window_length_c);
+          const int index = calculateIndex(N, window_length_c);
+          const double difference = fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index]));
+//          std::cout << "Comparing |"
+//              << scientific
+//              << smallgrid.sampled_grid[index] << " - "
+//              << checkgrid.sampled_grid[index] << "| = "
+//              << fabs((smallgrid.sampled_grid[index] - checkgrid.sampled_grid[index]))
+//              << " < " << check_threshold << " => "
+//              << fixed << setprecision(2)
+//              << 100*difference/threshold
+//              << " %, compare to " << 100*(1.+grid_difference/4.)*getBoundaryCaseFactor(N, window_length_c)
+//              << " %" << std::endl;
+          CPPUNIT_ASSERT( difference < check_threshold );
+        }
+  }
+}

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

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