src/gromacs/math/tests/gausstransform.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020, 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 = { 0.0047816522419452667236328125,
  94                                                         0.2613909542560577392578125,
  95                                                         0.65811407566070556640625,
  96                                                         0.07631497085094451904296875,
  97                                                         0.000407583254855126142501831054688 };
  98     EXPECT_THAT(expected, Pointwise(FloatEq(tolerance), viewOnResult));
  99 }
 100
 101 TEST(GaussianOn1DLattice, complementaryAmplitudesSumToZero)
 102 {
 103     const int              spreadWidth = 2;
 104     const real             sigma       = 0.57;
 105     const real             shift       = -0.2;
 106     GaussianOn1DLattice    gauss1d(spreadWidth, sigma);
 107     const auto             viewOnResult = gauss1d.view();
 108     FloatingPointTolerance tolerance(defaultFloatTolerance());
 109     const real             amplitude = 2.0;
 110     gauss1d.spread(amplitude, shift);
 111     std::vector<float> sumOfComplementaryGaussians;
 112     // keep a copy of the first Gaussian
 113     std::copy(std::begin(viewOnResult), std::end(viewOnResult), std::back_inserter(sumOfComplementaryGaussians));
 114
 115     gauss1d.spread(-amplitude, shift);
 116     // add the two spread Gaussians
 117     std::transform(std::begin(viewOnResult),
 118                    std::end(viewOnResult),
 119                    std::begin(sumOfComplementaryGaussians),
 120                    std::begin(sumOfComplementaryGaussians),
 121                    std::plus<>());
 122     // Expect all zeros
 123     std::array<float, 2 * spreadWidth + 1> expected = {};
 124     EXPECT_THAT(expected, Pointwise(FloatEq(tolerance), sumOfComplementaryGaussians));
 125 }
 126
 127 TEST(GaussianOn1DLattice, doesNotOverflowForLargeRange)
 128 {
 129     const int           spreadWidth = 200;
 130     const real          sigma       = 1;
 131     const real          shift       = -0.5;
 132     GaussianOn1DLattice gauss1d(spreadWidth, sigma);
 133     const auto          viewOnResult = gauss1d.view();
 134     const real          weightFirst  = 1;
 135     gauss1d.spread(weightFirst, shift);
 136
 137     for (size_t i = 0; i < 187; i++)
 138     {
 139         EXPECT_FLOAT_EQ(0, viewOnResult[i]);
 140     }
 141
 142     std::array<float, 12> expectedResult = {
 143         7.64165518045829568433117268116e-30, 4.57537550091150099862240391475e-25,
 144         1.00779356059942059652096780012e-20, 8.1662360418003390174698785664e-17,
 145         2.43432064870457987026952650922e-13, 2.66955679784075528004905208945e-10,
 146         1.07697594842193211661651730537e-07, 1.59837418323149904608726501465e-05,
 147         0.000872682663612067699432373046875, 0.01752830110490322113037109375,
 148         0.12951759994029998779296875,        0.3520653247833251953125
 149     };
 150
 151     for (size_t i = 188; i < 200; i++)
 152     {
 153         EXPECT_FLOAT_EQ(expectedResult[i - 188], viewOnResult[i]);
 154     }
 155
 156     for (size_t i = 200; i < 212; i++)
 157     {
 158         EXPECT_FLOAT_EQ(expectedResult[211 - i], viewOnResult[i]);
 159     }
 160
 161     EXPECT_FLOAT_EQ(4.69519506491195688717869977423e-35, viewOnResult[212]);
 162
 163     for (size_t i = 213; i < 2 * spreadWidth + 1; i++)
 164     {
 165         EXPECT_FLOAT_EQ(0, viewOnResult[i]);
 166     }
 167 }
 168
 169 class GaussTransformTest : public ::testing::Test
 170 {
 171 public:
 172     void isZeroWithinFloatTolerance()
 173     {
 174         for (const auto& x : gaussTransform_.constView())
 175         {
 176             EXPECT_FLOAT_EQ_TOL(0, x, tolerance_);
 177         }
 178     }
 179
 180 protected:
 181     extents<dynamic_extent, dynamic_extent, dynamic_extent> latticeExtent_ = { 3, 3, 3 };
 182     DVec                                                    sigma_         = { 1., 1., 1. };
 183     double                                                  nSigma_        = 5;
 184     GaussTransform3D             gaussTransform_ = { latticeExtent_, { sigma_, nSigma_ } };
 185     test::FloatingPointTolerance tolerance_      = test::defaultFloatTolerance();
 186     const RVec                   latticeCenter_  = { 1, 1, 1 };
 187 };
 188
 189 TEST_F(GaussTransformTest, isZeroUponConstruction)
 190 {
 191     isZeroWithinFloatTolerance();
 192 }
 193
 194 TEST_F(GaussTransformTest, isZeroAddingZeroAmplitudeGauss)
 195 {
 196     gaussTransform_.add({ latticeCenter_, 0. });
 197     isZeroWithinFloatTolerance();
 198 }
 199
 200 TEST_F(GaussTransformTest, isZeroAfterSettingZero)
 201 {
 202     gaussTransform_.add({ latticeCenter_, 1. });
 203     for (const auto value : gaussTransform_.constView())
 204     {
 205         EXPECT_GT(value, 0);
 206     }
 207     gaussTransform_.setZero();
 208     isZeroWithinFloatTolerance();
 209 }
 210
 211 TEST_F(GaussTransformTest, isZeroWhenOutsideRangeinX)
 212 {
 213     DVec coordinateOutsideX(-nSigma_ * sigma_[XX], 0, 0);
 214     gaussTransform_.add({ coordinateOutsideX.toRVec(), 1. });
 215     isZeroWithinFloatTolerance();
 216 }
 217
 218 TEST_F(GaussTransformTest, isZeroWhenOutsideRangeinY)
 219 {
 220     DVec coordinateOutsideY(0, -nSigma_ * sigma_[YY], 0);
 221     gaussTransform_.add({ coordinateOutsideY.toRVec(), 1. });
 222     isZeroWithinFloatTolerance();
 223 }
 224
 225 TEST_F(GaussTransformTest, isZeroWhenOutsideRangeinZ)
 226 {
 227     RVec coordinateOutsideZ(0, 0, -nSigma_ * sigma_[ZZ]);
 228     gaussTransform_.add({ coordinateOutsideZ.toRVec(), 1. });
 229     isZeroWithinFloatTolerance();
 230 }
 231
 232 TEST_F(GaussTransformTest, complementaryGaussAddToZero)
 233 {
 234     gaussTransform_.add({ latticeCenter_, -2.0 });
 235     gaussTransform_.add({ latticeCenter_, 0.8 });
 236     gaussTransform_.add({ latticeCenter_, 1.2 });
 237     isZeroWithinFloatTolerance();
 238 }
 239
 240 TEST_F(GaussTransformTest, centerGaussianInCubeHasExpectedValues)
 241 {
 242     gaussTransform_.add({ latticeCenter_, 1. });
 243     const float        center         = 0.0634936392307281494140625;      // center
 244     const float        f              = 0.0385108403861522674560546875;   // face
 245     const float        e              = 0.0233580060303211212158203125;   // edge
 246     const float        c              = 0.014167346991598606109619140625; // corner
 247     std::vector<float> expectedValues = { c, e, c, e, f,      e, c, e, c,
 248
 249                                           e, f, e, f, center, f, e, f, e,
 250
 251                                           c, e, c, e, f,      e, c, e, c };
 252     // This assignment to std::vector is needed because the view() on the GaussTransform (aka basic_mdspan) does not provide size() needed by Pointwise
 253     std::vector<float> gaussTransformVector;
 254     gaussTransformVector.assign(gaussTransform_.constView().data(),
 255                                 gaussTransform_.constView().data()
 256                                         + gaussTransform_.constView().mapping().required_span_size());
 257     EXPECT_THAT(expectedValues, testing::Pointwise(FloatEq(tolerance_), gaussTransformVector));
 258 }
 259
 260 TEST_F(GaussTransformTest, view)
 261 {
 262     gaussTransform_.add({ latticeCenter_, 1. });
 263     const float        f              = 0.0385108403861522674560546875;   // face
 264     const float        e              = 0.0233580060303211212158203125;   // edge
 265     const float        c              = 0.014167346991598606109619140625; // corner
 266     std::vector<float> expectedValues = { c, e, c, e, f,  e, c, e, c,
 267
 268                                           e, f, e, f, 0., f, e, f, e,
 269
 270                                           c, e, c, e, f,  e, c, e, c };
 271
 272     gaussTransform_.view()[1][1][1] = 0;
 273     // This assignment to std::vector is needed because the view() on the
 274     // GaussTransform (aka basic_mdspan) does not provide size() that is
 275     // needed by Pointwise
 276     // \todo Update when mdspan functionality is extended
 277     std::vector<float> gaussTransformVector;
 278     gaussTransformVector.assign(
 279             gaussTransform_.view().data(),
 280             gaussTransform_.view().data() + gaussTransform_.view().mapping().required_span_size());
 281     EXPECT_THAT(expectedValues, testing::Pointwise(FloatEq(tolerance_), gaussTransformVector));
 282 }
 283
 284 } // namespace
 285
 286 } // namespace test
 287
 288 } // namespace gmx