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