2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * Implements Gaussian function evaluations on lattices and related functionality
39 * \author Christian Blau <blau@kth.se>
41 * \ingroup module_math
45 #include "gausstransform.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/multidimarray.h"
54 #include "gromacs/math/utilities.h"
59 /********************************************************************
60 * GaussianOn1DLattice::Impl
63 class GaussianOn1DLattice::Impl
66 Impl(int numGridPointsForSpreadingHalfWidth
, real sigma
);
68 Impl(const Impl
&other
) = default;
69 Impl
&operator=(const Impl
&other
) = default;
71 /*! \brief evaluate Gaussian function at all lattice points
72 * \param[in] amplitude the amplitude of the Gaussian
73 * \param[in] dx distance from the center
75 void spread(double amplitude
, real dx
);
76 //! Largest distance in number of gridpoints from 0
77 int numGridPointsForSpreadingHalfWidth_
;
78 /*! \brief Avoid overflow for E2^offset and underflow for E3(i).
80 * Occurs when sigma is much smaller than numGridPointsForSpreadingHalfWidth_.
82 * E2^offset smaller than maximum float requires
83 * \f$exp(dx / (2*square(sigma))^numGridPointsForSpreadingHalfWidth_ \leq max_float \f$
84 * The maximum expected distance of the Gaussian center to the next lattice point is dx = 0.5,
85 * thus the maximum spread distance here is \f$4 * sigma^2 * \log(\mathrm{maxfloat})\f$ .
87 * E3(i) larger than minmium float requires
88 * exp(i^2 / 2*(sigma)^2) > min_float
89 * Thus the maximum spread distance here is \f$\sigma \sqrt(-2\log(\mathrm{minfloat}))\f$
91 int maxEvaluatedSpreadDistance_
;
92 //! Width of the Gaussian function
94 //! The result of the spreading calculation
95 std::vector
<float> spreadingResult_
;
96 //! Pre-calculated exp(-gridIndex^2/2 * (sigma^2)) named as in Greengard2004
97 std::vector
<float> e3_
;
98 /*! \brief Equal to std::floor(std::log(std::numeric_limits<float>::max())).
99 * Above expression is not constexpr and a const variable would implicitly delete default copy assignment.
100 * Therefore resorting to setting number manually.
102 static constexpr double c_logMaxFloat
= 88.72284;
103 static constexpr double c_logMinFloat
= -87.33654;
106 GaussianOn1DLattice::Impl::Impl(int numGridPointsForSpreadingHalfWidth
, real sigma
) :
107 numGridPointsForSpreadingHalfWidth_(numGridPointsForSpreadingHalfWidth
),
109 spreadingResult_(2 * numGridPointsForSpreadingHalfWidth
+ 1)
111 maxEvaluatedSpreadDistance_
= std::min(numGridPointsForSpreadingHalfWidth_
, static_cast<int>(std::floor(4 * square(sigma
) * c_logMaxFloat
)) - 1);
112 maxEvaluatedSpreadDistance_
= std::min(maxEvaluatedSpreadDistance_
, static_cast<int>(std::floor(sigma
* sqrt(-2.0 * c_logMinFloat
))) - 1);
114 std::generate_n(std::back_inserter(e3_
), maxEvaluatedSpreadDistance_
+ 1,
115 [sigma
, latticeIndex
= 0]() mutable {
116 return std::exp(-0.5 * square(latticeIndex
++ / sigma
));
119 std::fill(std::begin(spreadingResult_
), std::end(spreadingResult_
), 0.);
122 void GaussianOn1DLattice::Impl::spread(double amplitude
, real dx
)
124 /* The spreading routine implements the fast gaussian gridding as in
126 * Leslie Greengard and June-Yub Lee,
127 * "Accelerating the Nonuniform Fast Fourier Transform"
128 * SIAM REV 2004 Vol. 46, No. 3, pp. 443-454 DOI. 10.1137/S003614450343200X
130 * Following the naming conventions for e1, e2 and e3, nu = 1, m = numGridPointsForSpreadingHalfWidth_.
132 * Speed up is achieved by factorization of the exponential that is evaluted
133 * at regular lattice points i, where the distance from the
134 * Gaussian center is \f$x-i\f$:
137 * a * \exp(-(x^2-2*i*x+ i^2)/(2*\sigma^2)) =
138 * a * \exp(-x^2/2*\sigma^2) * \exp(x/\sigma^2)^i * \exp(i/2*sigma^2) =
139 * e_1(x) * e_2(x)^i * e_3(i)
142 * Requiring only two exp evaluations per spreading operation.
145 const double e1
= amplitude
* exp(-0.5 * dx
* dx
/ square(sigma_
)) / (sqrt(2 * M_PI
) * sigma_
);
146 spreadingResult_
[numGridPointsForSpreadingHalfWidth_
] = e1
;
148 const double e2
= exp(dx
/ square(sigma_
));
150 double e2pow
= e2
; //< powers of e2, e2^offset
152 // Move outwards from mid-point, using e2pow value for both points simultaneously
153 // o o o<----O---->o o o
154 for (int offset
= 1; offset
< maxEvaluatedSpreadDistance_
; offset
++)
156 const double e1_3
= e1
* e3_
[offset
];
157 spreadingResult_
[numGridPointsForSpreadingHalfWidth_
+ offset
] = e1_3
* e2pow
;
158 spreadingResult_
[numGridPointsForSpreadingHalfWidth_
- offset
] = e1_3
/ e2pow
;
161 // separate statement for gridpoints at the end of the range avoids
162 // overflow for large sigma and saves one e2 multiplication operation
163 spreadingResult_
[numGridPointsForSpreadingHalfWidth_
- maxEvaluatedSpreadDistance_
] = (e1
/ e2pow
) * e3_
[maxEvaluatedSpreadDistance_
];
164 spreadingResult_
[numGridPointsForSpreadingHalfWidth_
+ maxEvaluatedSpreadDistance_
] = (e1
* e2pow
) * e3_
[maxEvaluatedSpreadDistance_
];
167 /********************************************************************
168 * GaussianOn1DLattice
171 GaussianOn1DLattice::GaussianOn1DLattice(int numGridPointsForSpreadingHalfWidth_
, real sigma
) : impl_(new Impl(numGridPointsForSpreadingHalfWidth_
, sigma
))
175 GaussianOn1DLattice::~GaussianOn1DLattice () {}
177 void GaussianOn1DLattice::spread(double amplitude
, real dx
)
179 impl_
->spread(amplitude
, dx
);
182 ArrayRef
<const float> GaussianOn1DLattice::view()
184 return impl_
->spreadingResult_
;
187 GaussianOn1DLattice::GaussianOn1DLattice(const GaussianOn1DLattice
&other
)
188 : impl_(new Impl(*other
.impl_
))
192 GaussianOn1DLattice
&GaussianOn1DLattice::operator=(const GaussianOn1DLattice
&other
)
194 *impl_
= *other
.impl_
;
198 GaussianOn1DLattice::GaussianOn1DLattice(GaussianOn1DLattice
&&) noexcept
= default;
200 GaussianOn1DLattice
&GaussianOn1DLattice::operator=(GaussianOn1DLattice
&&) noexcept
= default;
205 //! rounds real-valued coordinate to the closest integer values
206 IVec
closestIntegerPoint(const RVec
&coordinate
)
209 roundToInt(coordinate
[XX
]),
210 roundToInt(coordinate
[YY
]),
211 roundToInt(coordinate
[ZZ
])
215 /*! \brief Substracts a range from a three-dimensional integer coordinate and ensures
216 * the resulting coordinate is within a lattice.
217 * \param[in] index point in lattice
218 * \param[in] range to be shifted
219 * \returns Shifted index or zero if shifted index is smaller than zero.
221 IVec
rangeBeginWithinLattice(const IVec
&index
, const IVec
&range
)
223 return elementWiseMax({0, 0, 0}, index
- range
);
226 /*! \brief Adds a range from a three-dimensional integer coordinate and ensures
227 * the resulting coordinate is within a lattice.
228 * \param[in] index point in lattice
229 * \param[in] extents extent of the lattice
230 * \param[in] range to be shifted
231 * \returns Shifted index or the lattice extent if shifted index is larger than the extent
233 IVec
rangeEndWithinLattice(const IVec
&index
, const dynamicExtents3D
&extents
, const IVec
&range
)
235 IVec
extentAsIvec(static_cast<int>(extents
.extent(ZZ
)), static_cast<int>(extents
.extent(YY
)), static_cast<int>(extents
.extent(XX
)));
236 return elementWiseMin(extentAsIvec
, index
+ range
);
242 /********************************************************************
243 * OuterProductEvaluator
246 mdspan
<const float, dynamic_extent
, dynamic_extent
>
247 OuterProductEvaluator::operator()(ArrayRef
<const float> x
, ArrayRef
<const float> y
)
249 data_
.resize(ssize(x
), ssize(y
));
250 for (gmx::index xIndex
= 0; xIndex
< ssize(x
); ++xIndex
)
252 const auto xValue
= x
[xIndex
];
253 std::transform(std::begin(y
), std::end(y
), begin(data_
.asView()[xIndex
]),
254 [xValue
](float yValue
) { return xValue
* yValue
; });
256 return data_
.asConstView();
259 /********************************************************************
263 IntegerBox::IntegerBox(const IVec
&begin
, const IVec
&end
) : begin_
{begin
}, end_
{
268 const IVec
&IntegerBox::begin() const{return begin_
; }
269 const IVec
&IntegerBox::end() const { return end_
; }
271 bool IntegerBox::empty() const { return !((begin_
[XX
] < end_
[XX
] ) && (begin_
[YY
] < end_
[YY
]) && (begin_
[ZZ
] < end_
[ZZ
])); }
273 IntegerBox
spreadRangeWithinLattice(const IVec
¢er
, dynamicExtents3D extent
, IVec range
)
275 const IVec begin
= rangeBeginWithinLattice(center
, range
);
276 const IVec end
= rangeEndWithinLattice(center
, extent
, range
);
281 /********************************************************************
282 * GaussianSpreadKernel
285 IVec
GaussianSpreadKernelParameters::Shape::latticeSpreadRange() const
287 DVec
range(std::ceil(sigma_
[XX
] * spreadWidthMultiplesOfSigma_
), std::ceil(sigma_
[YY
] * spreadWidthMultiplesOfSigma_
), std::ceil(sigma_
[ZZ
] * spreadWidthMultiplesOfSigma_
));
288 return range
.toIVec();
291 /********************************************************************
292 * GaussTransform3D::Impl
296 * Private implementation class for GaussTransform3D.
298 class GaussTransform3D::Impl
301 //! Construct from extent and spreading width and range
302 Impl(const dynamicExtents3D
&extent
,
303 const GaussianSpreadKernelParameters::Shape
&kernelShapeParameters
);
306 Impl(const Impl
&other
) = default;
308 Impl
&operator=(const Impl
&other
) = default;
309 //! Add another gaussian
310 void add(const GaussianSpreadKernelParameters::PositionAndAmplitude
&localParamters
);
311 //! The width of the Gaussian in lattice spacing units
312 BasicVector
<double> sigma_
;
313 //! The spread range in lattice points
315 //! The result of the Gauss transform
316 MultiDimArray
<std::vector
<float>, dynamicExtents3D
> data_
;
317 //! The outer product of a Gaussian along the z and y dimension
318 OuterProductEvaluator outerProductZY_
;
319 //! The three one-dimensional Gaussians, whose outer product is added to the Gauss transform
320 std::array
<GaussianOn1DLattice
, DIM
> gauss1d_
;
323 GaussTransform3D::Impl::Impl(const dynamicExtents3D
&extent
,
324 const GaussianSpreadKernelParameters::Shape
&kernelShapeParameters
)
325 : sigma_
{kernelShapeParameters
.sigma_
},
327 kernelShapeParameters
.latticeSpreadRange()
332 gauss1d_( {GaussianOn1DLattice(spreadRange_
[XX
], sigma_
[XX
]),
333 GaussianOn1DLattice(spreadRange_
[YY
], sigma_
[YY
]),
334 GaussianOn1DLattice(spreadRange_
[ZZ
], sigma_
[ZZ
]) } )
338 void GaussTransform3D::Impl::add(const GaussianSpreadKernelParameters::PositionAndAmplitude
&localParameters
)
340 const IVec closestLatticePoint
= closestIntegerPoint(localParameters
.coordinate_
);
341 const auto spreadRange
= spreadRangeWithinLattice(closestLatticePoint
, data_
.asView().extents(), spreadRange_
);
343 // do nothing if the added Gaussian will never reach the lattice
344 if (spreadRange
.empty())
349 for (int dimension
= XX
; dimension
<= ZZ
; ++dimension
)
351 // multiply with amplitude so that Gauss3D = (amplitude * Gauss_x) * Gauss_y * Gauss_z
352 const float gauss1DAmplitude
= dimension
> XX
? 1.0 : localParameters
.amplitude_
;
353 gauss1d_
[dimension
].spread(gauss1DAmplitude
, localParameters
.coordinate_
[dimension
] - closestLatticePoint
[dimension
]);
356 const auto spreadZY
= outerProductZY_(gauss1d_
[ZZ
].view(), gauss1d_
[YY
].view());
357 const auto spreadX
= gauss1d_
[XX
].view();
358 const IVec spreadGridOffset
= spreadRange_
- closestLatticePoint
;
360 // \todo optimize these loops if performance critical
361 // The looping strategy uses that the last, x-dimension is contiguous in the memory layout
362 for (int zLatticeIndex
= spreadRange
.begin()[ZZ
]; zLatticeIndex
< spreadRange
.end()[ZZ
]; ++zLatticeIndex
)
364 const auto zSlice
= data_
.asView()[zLatticeIndex
];
366 for (int yLatticeIndex
= spreadRange
.begin()[YY
]; yLatticeIndex
< spreadRange
.end()[YY
]; ++yLatticeIndex
)
368 const auto ySlice
= zSlice
[yLatticeIndex
];
369 const float zyPrefactor
= spreadZY(zLatticeIndex
+ spreadGridOffset
[ZZ
], yLatticeIndex
+ spreadGridOffset
[YY
]);
371 for (int xLatticeIndex
= spreadRange
.begin()[XX
]; xLatticeIndex
< spreadRange
.end()[XX
]; ++xLatticeIndex
)
373 const float xPrefactor
= spreadX
[xLatticeIndex
+ spreadGridOffset
[XX
]];
374 ySlice
[xLatticeIndex
] += zyPrefactor
* xPrefactor
;
380 /********************************************************************
384 GaussTransform3D::GaussTransform3D(const dynamicExtents3D
&extent
,
385 const GaussianSpreadKernelParameters::Shape
&kernelShapeParameters
) : impl_(new Impl(extent
, kernelShapeParameters
))
389 void GaussTransform3D::add(const GaussianSpreadKernelParameters::PositionAndAmplitude
&localParameters
)
391 impl_
->add(localParameters
);
394 void GaussTransform3D::setZero()
396 std::fill(begin(impl_
->data_
), end(impl_
->data_
), 0.);
399 basic_mdspan
<const float, dynamicExtents3D
> GaussTransform3D::view()
401 return impl_
->data_
.asConstView();
404 GaussTransform3D::~GaussTransform3D()
407 GaussTransform3D::GaussTransform3D(const GaussTransform3D
&other
)
408 : impl_(new Impl(*other
.impl_
))
412 GaussTransform3D
&GaussTransform3D::operator=(const GaussTransform3D
&other
)
414 *impl_
= *other
.impl_
;
418 GaussTransform3D::GaussTransform3D(GaussTransform3D
&&) noexcept
= default;
420 GaussTransform3D
&GaussTransform3D::operator=(GaussTransform3D
&&) noexcept
= default;