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 },
 102     gradient_{ referenceDensity.extents() }
 103 {
 104     const auto numVoxels = gradient_.asConstView().mapping().required_span_size();
 105     /* the gradient for the inner product measure of fit is constant and does not
 106      * depend on the compared density, so it is pre-computed here */
 107     std::transform(begin(referenceDensity_), end(referenceDensity_), begin(gradient_),
 108                    [numVoxels](float x) { return x / numVoxels; });
 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 },
 187     gradient_(referenceDensity.extents())
 188 {
 189 }
 190
 191 real DensitySimilarityRelativeEntropy::similarity(density comparedDensity)
 192 {
 193     if (comparedDensity.extents() != referenceDensity_.extents())
 194     {
 195         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 196     }
 197     return std::inner_product(begin(referenceDensity_), end(referenceDensity_),
 198                               begin(comparedDensity), 0., std::plus<>(), relativeEntropyAtVoxel);
 199 }
 200
 201 DensitySimilarityMeasure::density DensitySimilarityRelativeEntropy::gradient(density comparedDensity)
 202 {
 203     if (comparedDensity.extents() != referenceDensity_.extents())
 204     {
 205         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 206     }
 207     std::transform(begin(referenceDensity_), end(referenceDensity_), begin(comparedDensity),
 208                    begin(gradient_), relativeEntropyGradientAtVoxel);
 209     return gradient_.asConstView();
 210 }
 211
 212 std::unique_ptr<DensitySimilarityMeasureImpl> DensitySimilarityRelativeEntropy::clone()
 213 {
 214     return std::make_unique<DensitySimilarityRelativeEntropy>(referenceDensity_);
 215 }
 216
 217 /****************** Cross Correlation *****************************************/
 218
 219 //! Helper values for evaluating the cross correlation
 220 struct CrossCorrelationEvaluationHelperValues
 221 {
 222     //! The mean of the reference density
 223     real meanReference = 0;
 224     //! The mean of the compared density
 225     real meanComparison = 0;
 226     //! The sum of the squared reference density voxel values
 227     real referenceSquaredSum = 0;
 228     //! The sum of the squared compared density voxel values
 229     real comparisonSquaredSum = 0;
 230     //! The covariance of the reference and the compared density
 231     real covariance = 0;
 232 };
 233
 234 /*! \brief Calculate helper values for the cross-correlation.
 235
 236  * Enables numerically stable single-pass cross-correlation evaluation algorithm
 237  * as described in Bennett, J., Grout, R. , Pebay, P., Roe D., Thompson D.
 238  * "Numerically Stable, Single-Pass, Parallel Statistics Algorithms"
 239  * and implemented in boost's correlation coefficient
 240  */
 241 CrossCorrelationEvaluationHelperValues evaluateHelperValues(DensitySimilarityMeasure::density reference,
 242                                                             DensitySimilarityMeasure::density compared)
 243 {
 244     CrossCorrelationEvaluationHelperValues helperValues;
 245
 246     index i = 0;
 247
 248     auto referenceIterator = begin(reference);
 249     for (const real comp : compared)
 250     {
 251         const real refHelper        = *referenceIterator - helperValues.meanReference;
 252         const real comparisonHelper = comp - helperValues.meanComparison;
 253         helperValues.referenceSquaredSum += (i * square(refHelper)) / (i + 1);
 254         helperValues.comparisonSquaredSum += (i * square(comparisonHelper)) / (i + 1);
 255         helperValues.covariance += i * refHelper * comparisonHelper / (i + 1);
 256         helperValues.meanReference += refHelper / (i + 1);
 257         helperValues.meanComparison += comparisonHelper / (i + 1);
 258
 259         ++referenceIterator;
 260         ++i;
 261     }
 262
 263     return helperValues;
 264 }
 265
 266 //! Calculate a single cross correlation gradient entry at a voxel.
 267 class CrossCorrelationGradientAtVoxel
 268 {
 269 public:
 270     //! Set up the gradient calculation with pre-computed values
 271     CrossCorrelationGradientAtVoxel(const CrossCorrelationEvaluationHelperValues& preComputed) :
 272         prefactor_(evaluatePrefactor(preComputed.comparisonSquaredSum, preComputed.referenceSquaredSum)),
 273         comparisonPrefactor_(preComputed.covariance / preComputed.comparisonSquaredSum),
 274         meanReference_(preComputed.meanReference),
 275         meanComparison_(preComputed.meanComparison)
 276     {
 277     }
 278     //! Evaluate the cross correlation gradient at a voxel
 279     real operator()(real reference, real comparison)
 280     {
 281         return prefactor_
 282                * (reference - meanReference_ - comparisonPrefactor_ * (comparison - meanComparison_));
 283     }
 284
 285 private:
 286     static real evaluatePrefactor(real comparisonSquaredSum, real referenceSquaredSum)
 287     {
 288         GMX_ASSERT(comparisonSquaredSum > 0,
 289                    "Squared sum of comparison values needs to be larger than zero.");
 290         GMX_ASSERT(referenceSquaredSum > 0,
 291                    "Squared sum of reference values needs to be larger than zero.");
 292         return 1.0 / (sqrt(comparisonSquaredSum) * sqrt(referenceSquaredSum));
 293     }
 294     const real prefactor_;
 295     const real comparisonPrefactor_;
 296     const real meanReference_;
 297     const real meanComparison_;
 298 };
 299
 300 /*! \internal
 301  * \brief Implementation for DensitySimilarityCrossCorrelation.
 302  *
 303  * The similarity measure itself is documented in DensitySimilarityMeasureMethod::crossCorrelation.
 304  */
 305 class DensitySimilarityCrossCorrelation final : public DensitySimilarityMeasureImpl
 306 {
 307 public:
 308     //! Construct similarity measure by setting the reference density
 309     DensitySimilarityCrossCorrelation(density referenceDensity);
 310     //! The gradient for the cross correlation similarity measure
 311     density gradient(density comparedDensity) override;
 312     //! Clone this
 313     std::unique_ptr<DensitySimilarityMeasureImpl> clone() override;
 314     //! The similarity between reference density and compared density
 315     real similarity(density comparedDensity) override;
 316
 317 private:
 318     //! A view on the reference density
 319     const density referenceDensity_;
 320     //! Stores the gradient of the similarity measure in memory
 321     MultiDimArray<std::vector<float>, dynamicExtents3D> gradient_;
 322 };
 323
 324 DensitySimilarityCrossCorrelation::DensitySimilarityCrossCorrelation(density referenceDensity) :
 325     referenceDensity_{ referenceDensity },
 326     gradient_(referenceDensity.extents())
 327 {
 328 }
 329
 330 real DensitySimilarityCrossCorrelation::similarity(density comparedDensity)
 331 {
 332     if (comparedDensity.extents() != referenceDensity_.extents())
 333     {
 334         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 335     }
 336
 337     CrossCorrelationEvaluationHelperValues helperValues =
 338             evaluateHelperValues(referenceDensity_, comparedDensity);
 339
 340     if ((helperValues.referenceSquaredSum == 0) || (helperValues.comparisonSquaredSum == 0))
 341     {
 342         return 0;
 343     }
 344
 345     // To avoid numerical instability due to large squared density value sums
 346     // division is re-written to avoid multiplying two large numbers
 347     // as product of two separate divisions of smaller numbers
 348     const real covarianceSqrt = sqrt(fabs(helperValues.covariance));
 349     const int  sign           = helperValues.covariance > 0 ? 1 : -1;
 350     return sign * (covarianceSqrt / sqrt(helperValues.referenceSquaredSum))
 351            * (covarianceSqrt / sqrt(helperValues.comparisonSquaredSum));
 352 }
 353
 354 DensitySimilarityMeasure::density DensitySimilarityCrossCorrelation::gradient(density comparedDensity)
 355 {
 356     if (comparedDensity.extents() != referenceDensity_.extents())
 357     {
 358         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 359     }
 360
 361     CrossCorrelationEvaluationHelperValues helperValues =
 362             evaluateHelperValues(referenceDensity_, comparedDensity);
 363
 364     std::transform(begin(referenceDensity_), end(referenceDensity_), begin(comparedDensity),
 365                    begin(gradient_), CrossCorrelationGradientAtVoxel(helperValues));
 366
 367     return gradient_.asConstView();
 368 }
 369
 370 std::unique_ptr<DensitySimilarityMeasureImpl> DensitySimilarityCrossCorrelation::clone()
 371 {
 372     return std::make_unique<DensitySimilarityCrossCorrelation>(referenceDensity_);
 373 }
 374
 375
 376 } // namespace
 377
 378
 379 DensitySimilarityMeasure::DensitySimilarityMeasure(DensitySimilarityMeasureMethod method,
 380                                                    density                        referenceDensity)
 381 {
 382     // chose the implementation depending on the method of density comparison
 383     // throw an error if the method is not known
 384     switch (method)
 385     {
 386         case DensitySimilarityMeasureMethod::innerProduct:
 387             impl_ = std::make_unique<DensitySimilarityInnerProduct>(referenceDensity);
 388             break;
 389         case DensitySimilarityMeasureMethod::relativeEntropy:
 390             impl_ = std::make_unique<DensitySimilarityRelativeEntropy>(referenceDensity);
 391             break;
 392         case DensitySimilarityMeasureMethod::crossCorrelation:
 393             impl_ = std::make_unique<DensitySimilarityCrossCorrelation>(referenceDensity);
 394             break;
 395         default: GMX_THROW(NotImplementedError("Similarity measure not implemented."));
 396     }
 397 }
 398
 399 DensitySimilarityMeasure::density DensitySimilarityMeasure::gradient(density comparedDensity)
 400 {
 401     return impl_->gradient(comparedDensity);
 402 }
 403
 404 real DensitySimilarityMeasure::similarity(density comparedDensity)
 405 {
 406     return impl_->similarity(comparedDensity);
 407 }
 408
 409 DensitySimilarityMeasure::~DensitySimilarityMeasure() = default;
 410
 411 DensitySimilarityMeasure::DensitySimilarityMeasure(const DensitySimilarityMeasure& other) :
 412     impl_(other.impl_->clone())
 413 {
 414 }
 415
 416 DensitySimilarityMeasure& DensitySimilarityMeasure::operator=(const DensitySimilarityMeasure& other)
 417 {
 418     impl_ = other.impl_->clone();
 419     return *this;
 420 }
 421
 422 DensitySimilarityMeasure::DensitySimilarityMeasure(DensitySimilarityMeasure&&) noexcept = default;
 423
 424 DensitySimilarityMeasure& DensitySimilarityMeasure::operator=(DensitySimilarityMeasure&&) noexcept = default;
 425
 426 } // namespace gmx