-/*! \brief
- * Returns an approximation of the geometry factor used for initializing the AWH update size.
- *
- * The geometry factor is defined as the following sum of Gaussians:
- * sum_{k!=0} exp(-0.5*(k*pi*x)^2)/(pi*k)^2,
- * where k is a xArray.size()-dimensional integer vector with k_i in {0,1,..}.
- *
- * \param[in] xArray Array to evaluate.
- * \returns the geometry factor.
- */
-double gaussianGeometryFactor(gmx::ArrayRef<const double> xArray)
-{
- /* For convenience we give the geometry factor function a name: zeta(x) */
- constexpr size_t tableSize = 5;
- std::array<const double, tableSize> xTabulated = { { 1e-5, 1e-4, 1e-3, 1e-2, 1e-1 } };
- std::array<const double, tableSize> zetaTable1d = { { 0.166536811948, 0.16653116886, 0.166250075882,
- 0.162701098306, 0.129272430287 } };
- std::array<const double, tableSize> zetaTable2d = { { 2.31985974274, 1.86307292523, 1.38159772648,
- 0.897554759158, 0.405578211115 } };
-
- gmx::ArrayRef<const double> zetaTable;
-
- if (xArray.size() == 1)
- {
- zetaTable = zetaTable1d;
- }
- else if (xArray.size() == 2)
- { // NOLINT bugprone-branch-clone
- zetaTable = zetaTable2d;
- }
- else
- {
- /* TODO... but this is anyway a rough estimate and > 2 dimensions is not so popular.
- * Remove the above NOLINT when addressing this */
- zetaTable = zetaTable2d;
- }
-
- /* TODO. Really zeta is a function of an ndim-dimensional vector x and we shoudl have a ndim-dimensional lookup-table.
- Here we take the geometric average of the components of x which is ok if the x-components are not very different. */
- double xScalar = 1;
- for (const double& x : xArray)
- {
- xScalar *= x;
- }
-
- GMX_ASSERT(!xArray.empty(), "We should have a non-empty input array");
- xScalar = std::pow(xScalar, 1.0 / xArray.size());
-
- /* Look up zeta(x) */
- size_t xIndex = 0;
- while ((xIndex < xTabulated.size()) && (xScalar > xTabulated[xIndex]))
- {
- xIndex++;
- }
-
- double zEstimate;
- if (xIndex == xTabulated.size())
- {
- /* Take last value */
- zEstimate = zetaTable[xTabulated.size() - 1];
- }
- else if (xIndex == 0)
- {
- zEstimate = zetaTable[xIndex];
- }
- else
- {
- /* Interpolate */
- double x0 = xTabulated[xIndex - 1];
- double x1 = xTabulated[xIndex];
- double w = (xScalar - x0) / (x1 - x0);
- zEstimate = w * zetaTable[xIndex - 1] + (1 - w) * zetaTable[xIndex];
- }
-
- return zEstimate;
-}
-