src/gromacs/math/densityfit.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020,2021, 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  * Implements density similarity measures and their derivatives.
  38  *
  39  * \author Christian Blau <blau@kth.se>
  40  * \ingroup module_math
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "densityfit.h"
  45
  46 #include <algorithm>
  47 #include <numeric>
  48
  49 #include "gromacs/math/multidimarray.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/utility/exceptions.h"
  52
  53 namespace gmx
  54 {
  55
  56 class DensitySimilarityMeasureImpl
  57 {
  58 public:
  59     virtual ~DensitySimilarityMeasureImpl();
  60     //! convenience typedef
  61     using density = DensitySimilarityMeasure::density;
  62     //! \copydoc DensitySimilarityMeasure::gradient(DensitySimilarityMeasure::density comparedDensity)
  63     virtual density gradient(density comparedDensity) = 0;
  64     //! \copydoc DensitySimilarityMeasure::similarity(density comparedDensity)
  65     virtual real similarity(density comparedDensity) = 0;
  66     //! clone to allow copy operations
  67     virtual std::unique_ptr<DensitySimilarityMeasureImpl> clone() = 0;
  68 };
  69 DensitySimilarityMeasureImpl::~DensitySimilarityMeasureImpl() = default;
  70
  71 namespace
  72 {
  73
  74 /****************** Inner Product *********************************************/
  75
  76 /*! \internal
  77  * \brief Implementation for DensitySimilarityInnerProduct.
  78  *
  79  * The similarity measure itself is documented in DensitySimilarityMeasureMethod::innerProduct.
  80  */
  81 class DensitySimilarityInnerProduct final : public DensitySimilarityMeasureImpl
  82 {
  83 public:
  84     //! Construct similarity measure by setting the reference density
  85     DensitySimilarityInnerProduct(density referenceDensity);
  86     //! The gradient for the inner product similarity measure is the reference density divided by the number of voxels
  87     density gradient(density comparedDensity) override;
  88     //! Clone this
  89     std::unique_ptr<DensitySimilarityMeasureImpl> clone() override;
  90     //! The similarity between reference density and compared density
  91     real similarity(density comparedDensity) override;
  92
  93 private:
  94     //! A view on the reference density
  95     const density referenceDensity_;
  96     //! Stores the gradient of the similarity measure in memory
  97     MultiDimArray<std::vector<float>, dynamicExtents3D> gradient_;
  98 };
  99
 100 DensitySimilarityInnerProduct::DensitySimilarityInnerProduct(density referenceDensity) :
 101     referenceDensity_{ referenceDensity }, gradient_{ referenceDensity.extents() }
 102 {
 103     const auto numVoxels = gradient_.asConstView().mapping().required_span_size();
 104     /* the gradient for the inner product measure of fit is constant and does not
 105      * depend on the compared density, so it is pre-computed here */
 106     std::transform(begin(referenceDensity_), end(referenceDensity_), begin(gradient_), [numVoxels](float x) {
 107         return x / numVoxels;
 108     });
 109 }
 110
 111 real DensitySimilarityInnerProduct::similarity(density comparedDensity)
 112 {
 113     if (comparedDensity.extents() != referenceDensity_.extents())
 114     {
 115         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 116     }
 117     /* the similarity measure uses the gradient instead of the reference,
 118      * here, because it is the reference density divided by the number of voxels */
 119     return std::inner_product(begin(gradient_), end(gradient_), begin(comparedDensity), 0.);
 120 }
 121
 122 DensitySimilarityMeasure::density DensitySimilarityInnerProduct::gradient(density comparedDensity)
 123 {
 124     /* even though the gradient density does not depend on the compad density,
 125      * still checking the extents to make sure we're consistent */
 126     if (comparedDensity.extents() != referenceDensity_.extents())
 127     {
 128         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 129     }
 130
 131     return gradient_.asConstView();
 132 }
 133
 134 std::unique_ptr<DensitySimilarityMeasureImpl> DensitySimilarityInnerProduct::clone()
 135 {
 136     return std::make_unique<DensitySimilarityInnerProduct>(referenceDensity_);
 137 }
 138
 139 /****************** Relative Entropy *****************************************/
 140
 141 //! Calculate a single summand in the relative entropy sum.
 142 real relativeEntropyAtVoxel(real reference, real comparison)
 143 {
 144     if ((reference > 0) && (comparison > 0))
 145     {
 146         return reference * (std::log(comparison / reference));
 147     }
 148     return 0.;
 149 }
 150
 151 //! Calculate a single relative entropy gradient entry at a voxel.
 152 real relativeEntropyGradientAtVoxel(real reference, real comparison)
 153 {
 154     if ((reference > 0) && (comparison > 0))
 155     {
 156         return reference / comparison;
 157     }
 158     return 0.;
 159 }
 160
 161 /*! \internal
 162  * \brief Implementation for DensitySimilarityRelativeEntropy.
 163  *
 164  * The similarity measure itself is documented in DensitySimilarityMeasureMethod::RelativeEntropy.
 165  */
 166 class DensitySimilarityRelativeEntropy final : public DensitySimilarityMeasureImpl
 167 {
 168 public:
 169     //! Construct similarity measure by setting the reference density
 170     DensitySimilarityRelativeEntropy(density referenceDensity);
 171     //! The gradient for the relative entropy similarity measure
 172     density gradient(density comparedDensity) override;
 173     //! Clone this
 174     std::unique_ptr<DensitySimilarityMeasureImpl> clone() override;
 175     //! The similarity between reference density and compared density
 176     real similarity(density comparedDensity) override;
 177
 178 private:
 179     //! A view on the reference density
 180     const density referenceDensity_;
 181     //! Stores the gradient of the similarity measure in memory
 182     MultiDimArray<std::vector<float>, dynamicExtents3D> gradient_;
 183 };
 184
 185 DensitySimilarityRelativeEntropy::DensitySimilarityRelativeEntropy(density referenceDensity) :
 186     referenceDensity_{ referenceDensity }, gradient_(referenceDensity.extents())
 187 {
 188 }
 189
 190 real DensitySimilarityRelativeEntropy::similarity(density comparedDensity)
 191 {
 192     if (comparedDensity.extents() != referenceDensity_.extents())
 193     {
 194         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 195     }
 196     return std::inner_product(begin(referenceDensity_),
 197                               end(referenceDensity_),
 198                               begin(comparedDensity),
 199                               0.,
 200                               std::plus<>(),
 201                               relativeEntropyAtVoxel);
 202 }
 203
 204 DensitySimilarityMeasure::density DensitySimilarityRelativeEntropy::gradient(density comparedDensity)
 205 {
 206     if (comparedDensity.extents() != referenceDensity_.extents())
 207     {
 208         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 209     }
 210     std::transform(begin(referenceDensity_),
 211                    end(referenceDensity_),
 212                    begin(comparedDensity),
 213                    begin(gradient_),
 214                    relativeEntropyGradientAtVoxel);
 215     return gradient_.asConstView();
 216 }
 217
 218 std::unique_ptr<DensitySimilarityMeasureImpl> DensitySimilarityRelativeEntropy::clone()
 219 {
 220     return std::make_unique<DensitySimilarityRelativeEntropy>(referenceDensity_);
 221 }
 222
 223 /****************** Cross Correlation *****************************************/
 224
 225 //! Helper values for evaluating the cross correlation
 226 struct CrossCorrelationEvaluationHelperValues
 227 {
 228     //! The mean of the reference density
 229     real meanReference = 0;
 230     //! The mean of the compared density
 231     real meanComparison = 0;
 232     //! The sum of the squared reference density voxel values
 233     real referenceSquaredSum = 0;
 234     //! The sum of the squared compared density voxel values
 235     real comparisonSquaredSum = 0;
 236     //! The covariance of the reference and the compared density
 237     real covariance = 0;
 238 };
 239
 240 /*! \brief Calculate helper values for the cross-correlation.
 241
 242  * Enables numerically stable single-pass cross-correlation evaluation algorithm
 243  * as described in Bennett, J., Grout, R. , Pebay, P., Roe D., Thompson D.
 244  * "Numerically Stable, Single-Pass, Parallel Statistics Algorithms"
 245  * and implemented in boost's correlation coefficient
 246  */
 247 CrossCorrelationEvaluationHelperValues evaluateHelperValues(DensitySimilarityMeasure::density reference,
 248                                                             DensitySimilarityMeasure::density compared)
 249 {
 250     CrossCorrelationEvaluationHelperValues helperValues;
 251
 252     index i = 0;
 253
 254     const auto* referenceIterator = begin(reference);
 255     for (const real comp : compared)
 256     {
 257         const real refHelper        = *referenceIterator - helperValues.meanReference;
 258         const real comparisonHelper = comp - helperValues.meanComparison;
 259         helperValues.referenceSquaredSum += (i * square(refHelper)) / (i + 1);
 260         helperValues.comparisonSquaredSum += (i * square(comparisonHelper)) / (i + 1);
 261         helperValues.covariance += i * refHelper * comparisonHelper / (i + 1);
 262         helperValues.meanReference += refHelper / (i + 1);
 263         helperValues.meanComparison += comparisonHelper / (i + 1);
 264
 265         ++referenceIterator;
 266         ++i;
 267     }
 268
 269     return helperValues;
 270 }
 271
 272 //! Calculate a single cross correlation gradient entry at a voxel.
 273 class CrossCorrelationGradientAtVoxel
 274 {
 275 public:
 276     //! Set up the gradient calculation with pre-computed values
 277     CrossCorrelationGradientAtVoxel(const CrossCorrelationEvaluationHelperValues& preComputed) :
 278         prefactor_(evaluatePrefactor(preComputed.comparisonSquaredSum, preComputed.referenceSquaredSum)),
 279         comparisonPrefactor_(preComputed.covariance / preComputed.comparisonSquaredSum),
 280         meanReference_(preComputed.meanReference),
 281         meanComparison_(preComputed.meanComparison)
 282     {
 283     }
 284     //! Evaluate the cross correlation gradient at a voxel
 285     real operator()(real reference, real comparison) const
 286     {
 287         return prefactor_
 288                * (reference - meanReference_ - comparisonPrefactor_ * (comparison - meanComparison_));
 289     }
 290
 291 private:
 292     static real evaluatePrefactor(real comparisonSquaredSum, real referenceSquaredSum)
 293     {
 294         GMX_ASSERT(comparisonSquaredSum > 0,
 295                    "Squared sum of comparison values needs to be larger than zero.");
 296         GMX_ASSERT(referenceSquaredSum > 0,
 297                    "Squared sum of reference values needs to be larger than zero.");
 298         return 1.0 / (sqrt(comparisonSquaredSum) * sqrt(referenceSquaredSum));
 299     }
 300     const real prefactor_;
 301     const real comparisonPrefactor_;
 302     const real meanReference_;
 303     const real meanComparison_;
 304 };
 305
 306 /*! \internal
 307  * \brief Implementation for DensitySimilarityCrossCorrelation.
 308  *
 309  * The similarity measure itself is documented in DensitySimilarityMeasureMethod::crossCorrelation.
 310  */
 311 class DensitySimilarityCrossCorrelation final : public DensitySimilarityMeasureImpl
 312 {
 313 public:
 314     //! Construct similarity measure by setting the reference density
 315     DensitySimilarityCrossCorrelation(density referenceDensity);
 316     //! The gradient for the cross correlation similarity measure
 317     density gradient(density comparedDensity) override;
 318     //! Clone this
 319     std::unique_ptr<DensitySimilarityMeasureImpl> clone() override;
 320     //! The similarity between reference density and compared density
 321     real similarity(density comparedDensity) override;
 322
 323 private:
 324     //! A view on the reference density
 325     const density referenceDensity_;
 326     //! Stores the gradient of the similarity measure in memory
 327     MultiDimArray<std::vector<float>, dynamicExtents3D> gradient_;
 328 };
 329
 330 DensitySimilarityCrossCorrelation::DensitySimilarityCrossCorrelation(density referenceDensity) :
 331     referenceDensity_{ referenceDensity }, gradient_(referenceDensity.extents())
 332 {
 333 }
 334
 335 real DensitySimilarityCrossCorrelation::similarity(density comparedDensity)
 336 {
 337     if (comparedDensity.extents() != referenceDensity_.extents())
 338     {
 339         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 340     }
 341
 342     CrossCorrelationEvaluationHelperValues helperValues =
 343             evaluateHelperValues(referenceDensity_, comparedDensity);
 344
 345     if ((helperValues.referenceSquaredSum == 0) || (helperValues.comparisonSquaredSum == 0))
 346     {
 347         return 0;
 348     }
 349
 350     // To avoid numerical instability due to large squared density value sums
 351     // division is re-written to avoid multiplying two large numbers
 352     // as product of two separate divisions of smaller numbers
 353     const real covarianceSqrt = sqrt(fabs(helperValues.covariance));
 354     const int  sign           = helperValues.covariance > 0 ? 1 : -1;
 355     return sign * (covarianceSqrt / sqrt(helperValues.referenceSquaredSum))
 356            * (covarianceSqrt / sqrt(helperValues.comparisonSquaredSum));
 357 }
 358
 359 DensitySimilarityMeasure::density DensitySimilarityCrossCorrelation::gradient(density comparedDensity)
 360 {
 361     if (comparedDensity.extents() != referenceDensity_.extents())
 362     {
 363         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 364     }
 365
 366     CrossCorrelationEvaluationHelperValues helperValues =
 367             evaluateHelperValues(referenceDensity_, comparedDensity);
 368
 369     std::transform(begin(referenceDensity_),
 370                    end(referenceDensity_),
 371                    begin(comparedDensity),
 372                    begin(gradient_),
 373                    CrossCorrelationGradientAtVoxel(helperValues));
 374
 375     return gradient_.asConstView();
 376 }
 377
 378 std::unique_ptr<DensitySimilarityMeasureImpl> DensitySimilarityCrossCorrelation::clone()
 379 {
 380     return std::make_unique<DensitySimilarityCrossCorrelation>(referenceDensity_);
 381 }
 382
 383
 384 } // namespace
 385
 386
 387 DensitySimilarityMeasure::DensitySimilarityMeasure(DensitySimilarityMeasureMethod method,
 388                                                    density                        referenceDensity)
 389 {
 390     // chose the implementation depending on the method of density comparison
 391     // throw an error if the method is not known
 392     switch (method)
 393     {
 394         case DensitySimilarityMeasureMethod::innerProduct:
 395             impl_ = std::make_unique<DensitySimilarityInnerProduct>(referenceDensity);
 396             break;
 397         case DensitySimilarityMeasureMethod::relativeEntropy:
 398             impl_ = std::make_unique<DensitySimilarityRelativeEntropy>(referenceDensity);
 399             break;
 400         case DensitySimilarityMeasureMethod::crossCorrelation:
 401             impl_ = std::make_unique<DensitySimilarityCrossCorrelation>(referenceDensity);
 402             break;
 403         default: GMX_THROW(NotImplementedError("Similarity measure not implemented."));
 404     }
 405 }
 406
 407 DensitySimilarityMeasure::density DensitySimilarityMeasure::gradient(density comparedDensity)
 408 {
 409     return impl_->gradient(comparedDensity);
 410 }
 411
 412 real DensitySimilarityMeasure::similarity(density comparedDensity)
 413 {
 414     return impl_->similarity(comparedDensity);
 415 }
 416
 417 DensitySimilarityMeasure::~DensitySimilarityMeasure() = default;
 418
 419 DensitySimilarityMeasure::DensitySimilarityMeasure(const DensitySimilarityMeasure& other) :
 420     impl_(other.impl_->clone())
 421 {
 422 }
 423
 424 DensitySimilarityMeasure& DensitySimilarityMeasure::operator=(const DensitySimilarityMeasure& other)
 425 {
 426     impl_ = other.impl_->clone();
 427     return *this;
 428 }
 429
 430 DensitySimilarityMeasure::DensitySimilarityMeasure(DensitySimilarityMeasure&&) noexcept = default;
 431
 432 DensitySimilarityMeasure& DensitySimilarityMeasure::operator=(DensitySimilarityMeasure&&) noexcept = default;
 433
 434 } // namespace gmx