src/gromacs/math/tests/gausstransform.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Tests for Gaussian spreading
  38  *
  39  * \author Christian Blau <blau@kth.se>
  40  * \ingroup module_math
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "gromacs/math/gausstransform.h"
  45
  46 #include <array>
  47 #include <numeric>
  48 #include <vector>
  49
  50 #include <gmock/gmock.h>
  51 #include <gtest/gtest.h>
  52
  53 #include "gromacs/math/vec.h"
  54 #include "gromacs/mdspan/extensions.h"
  55
  56 #include "testutils/testasserts.h"
  57 #include "testutils/testmatchers.h"
  58
  59 namespace gmx
  60 {
  61
  62 namespace test
  63 {
  64
  65 namespace
  66 {
  67
  68 TEST(GaussianOn1DLattice, sumsCloseToOne)
  69 {
  70     const int              spreadWidth = 7;
  71     const real             sigma       = 1.9;
  72     const real             shift       = 0.1;
  73     GaussianOn1DLattice    gauss1d(spreadWidth, sigma);
  74     const real             amplitude = 1;
  75     gauss1d.spread(amplitude, shift);
  76     auto                   sumOverLattice = std::accumulate(std::begin(gauss1d.view()), std::end(gauss1d.view()), 0.);
  77
  78     FloatingPointTolerance tolerance(defaultFloatTolerance());
  79     // The sum over the lattice should be roughly one
  80     EXPECT_FLOAT_EQ_TOL(0.99993300437927246, sumOverLattice, tolerance);
  81 }
  82
  83 TEST(GaussianOn1DLattice, isCorrect)
  84 {
  85     const int                             spreadWidth = 2;
  86     const real                            sigma       = 0.57;
  87     const real                            shift       = -0.2;
  88     GaussianOn1DLattice                   gauss1d(spreadWidth, sigma);
  89     const auto                            viewOnResult = gauss1d.view();
  90     FloatingPointTolerance                tolerance(defaultFloatTolerance());
  91     const real                            amplitude = 1.0;
  92     gauss1d.spread(amplitude, shift);
  93     std::array<float, 2 *spreadWidth + 1> expected =  {
  94         0.0047816522419452667236328125, 0.2613909542560577392578125,
  95         0.65811407566070556640625, 0.07631497085094451904296875, 0.000407583254855126142501831054688
  96     };
  97     EXPECT_THAT(expected, Pointwise(FloatEq(tolerance), viewOnResult));
  98 }
  99
 100 TEST(GaussianOn1DLattice, complementaryAmplitudesSumToZero)
 101 {
 102     const int              spreadWidth = 2;
 103     const real             sigma       = 0.57;
 104     const real             shift       = -0.2;
 105     GaussianOn1DLattice    gauss1d(spreadWidth, sigma);
 106     const auto             viewOnResult = gauss1d.view();
 107     FloatingPointTolerance tolerance(defaultFloatTolerance());
 108     const real             amplitude = 2.0;
 109     gauss1d.spread(amplitude, shift);
 110     std::vector<float>     sumOfComplementaryGaussians;
 111     // keep a copy of the first Gaussian
 112     std::copy(std::begin(viewOnResult), std::end(viewOnResult),
 113               std::back_inserter(sumOfComplementaryGaussians));
 114
 115     gauss1d.spread(-amplitude, shift);
 116     // add the two spread Gaussians
 117     std::transform(std::begin(viewOnResult), std::end(viewOnResult),
 118                    std::begin(sumOfComplementaryGaussians), std::begin(sumOfComplementaryGaussians), std::plus<>());
 119     // Expect all zeros
 120     std::array<float, 2 *spreadWidth + 1> expected = {};
 121     EXPECT_THAT(expected, Pointwise(FloatEq(tolerance), sumOfComplementaryGaussians));
 122 }
 123
 124 TEST(GaussianOn1DLattice, doesNotOverflowForLargeRange)
 125 {
 126     const int           spreadWidth = 200;
 127     const real          sigma       = 1;
 128     const real          shift       = -0.5;
 129     GaussianOn1DLattice gauss1d(spreadWidth, sigma);
 130     const auto          viewOnResult = gauss1d.view();
 131     const real          weightFirst  = 1;
 132     gauss1d.spread(weightFirst, shift);
 133
 134     for (size_t i = 0; i < 187; i++)
 135     {
 136         EXPECT_FLOAT_EQ(0, viewOnResult[i]);
 137     }
 138
 139     std::array<float, 12> expectedResult = {
 140         7.64165518045829568433117268116e-30, 4.57537550091150099862240391475e-25,
 141         1.00779356059942059652096780012e-20, 8.1662360418003390174698785664e-17,
 142         2.43432064870457987026952650922e-13, 2.66955679784075528004905208945e-10,
 143         1.07697594842193211661651730537e-07, 1.59837418323149904608726501465e-05,
 144         0.000872682663612067699432373046875, 0.01752830110490322113037109375,
 145         0.12951759994029998779296875, 0.3520653247833251953125
 146     };
 147
 148     for (size_t i = 188; i < 200; i++)
 149     {
 150         EXPECT_FLOAT_EQ(expectedResult[i-188], viewOnResult[i]);
 151     }
 152
 153     for (size_t i = 200; i < 212; i++)
 154     {
 155         EXPECT_FLOAT_EQ(expectedResult[211-i], viewOnResult[i]);
 156     }
 157
 158     EXPECT_FLOAT_EQ(4.69519506491195688717869977423e-35, viewOnResult[212]);
 159
 160     for (size_t i = 213; i < 2*spreadWidth+1; i++)
 161     {
 162         EXPECT_FLOAT_EQ(0, viewOnResult[i]);
 163     }
 164 }
 165
 166 class GaussTransformTest : public ::testing::Test
 167 {
 168     public:
 169         void isZeroWithinFloatTolerance()
 170         {
 171             for (const auto &x : gaussTransform_.view())
 172             {
 173                 EXPECT_FLOAT_EQ_TOL(0, x, tolerance_);
 174             }
 175         }
 176     protected:
 177         extents<dynamic_extent, dynamic_extent, dynamic_extent> latticeExtent_ = {3, 3, 3 };
 178         DVec                         sigma_          = {1., 1., 1.};
 179         double                       nSigma_         = 5;
 180         GaussTransform3D             gaussTransform_ = {latticeExtent_, {sigma_, nSigma_}};
 181         test::FloatingPointTolerance tolerance_      = test::defaultFloatTolerance();
 182         const RVec                   latticeCenter_  = {1, 1, 1};
 183 };
 184
 185 TEST_F(GaussTransformTest, isZeroUponConstruction)
 186 {
 187     isZeroWithinFloatTolerance();
 188 }
 189
 190 TEST_F(GaussTransformTest, isZeroAddingZeroAmplitudeGauss)
 191 {
 192     gaussTransform_.add({latticeCenter_, 0.});
 193     isZeroWithinFloatTolerance();
 194 }
 195
 196 TEST_F(GaussTransformTest, isZeroAfterSettingZero)
 197 {
 198     gaussTransform_.add({latticeCenter_, 1.});
 199     for (const auto value : gaussTransform_.view())
 200     {
 201         EXPECT_GT(value, 0);
 202     }
 203     gaussTransform_.setZero();
 204     isZeroWithinFloatTolerance();
 205 }
 206
 207 TEST_F(GaussTransformTest, isZeroWhenOutsideRangeinX)
 208 {
 209     DVec coordinateOutsideX(-nSigma_ * sigma_[XX], 0, 0);
 210     gaussTransform_.add({coordinateOutsideX.toRVec(), 1.});
 211     isZeroWithinFloatTolerance();
 212 }
 213
 214 TEST_F(GaussTransformTest, isZeroWhenOutsideRangeinY)
 215 {
 216     DVec coordinateOutsideY(0, -nSigma_ * sigma_[YY], 0);
 217     gaussTransform_.add({coordinateOutsideY.toRVec(), 1.});
 218     isZeroWithinFloatTolerance();
 219 }
 220
 221 TEST_F(GaussTransformTest, isZeroWhenOutsideRangeinZ)
 222 {
 223     RVec coordinateOutsideZ(0, 0, -nSigma_ * sigma_[ZZ]);
 224     gaussTransform_.add({coordinateOutsideZ.toRVec(), 1.});
 225     isZeroWithinFloatTolerance();
 226 }
 227
 228 TEST_F(GaussTransformTest, complementaryGaussAddToZero)
 229 {
 230     gaussTransform_.add({latticeCenter_, -2.0});
 231     gaussTransform_.add({latticeCenter_,  0.8});
 232     gaussTransform_.add({latticeCenter_,  1.2});
 233     isZeroWithinFloatTolerance();
 234 }
 235
 236 TEST_F(GaussTransformTest, centerGaussianInCubeHasExpectedValues)
 237 {
 238     gaussTransform_.add({latticeCenter_, 1.});
 239     const float        center         = 0.0634936392307281494140625;      // center
 240     const float        f              = 0.0385108403861522674560546875;   // face
 241     const float        e              = 0.0233580060303211212158203125;   // edge
 242     const float        c              = 0.014167346991598606109619140625; // corner
 243     std::vector<float> expectedValues = {
 244         c, e, c,
 245         e, f, e,
 246         c, e, c,
 247
 248         e, f, e,
 249         f, center, f,
 250         e, f, e,
 251
 252         c, e, c,
 253         e, f, e,
 254         c, e, c
 255     };
 256     // This assignment to std::vector is needed because the view() on the GaussTransform (aka basic_mdspan) does not provide size() needed by Pointwise
 257     std::vector<float> gaussTransformVector;
 258     gaussTransformVector.assign(gaussTransform_.view().data(), gaussTransform_.view().data() + gaussTransform_.view().mapping().required_span_size());
 259     EXPECT_THAT(expectedValues, testing::Pointwise(FloatEq(tolerance_), gaussTransformVector));
 260 }
 261
 262 } // namespace
 263
 264 } // namespace test
 265
 266 } // namespace gmx