src/gromacs/math/densityfit.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   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/vec.h"
  50 #include "gromacs/utility/exceptions.h"
  51
  52 #include "multidimarray.h"
  53
  54 namespace gmx
  55 {
  56
  57 class DensitySimilarityMeasureImpl
  58 {
  59     public:
  60         virtual ~DensitySimilarityMeasureImpl();
  61         //! convenience typedef
  62         using density = DensitySimilarityMeasure::density;
  63         //! \copydoc DensitySimilarityMeasure::gradient(DensitySimilarityMeasure::density comparedDensity)
  64         virtual density gradient(density comparedDensity) = 0;
  65         //! \copydoc DensitySimilarityMeasure::gradient(density comparedDensity)
  66         virtual float similarity(density comparedDensity) = 0;
  67         //! clone to allow copy operations
  68         virtual std::unique_ptr<DensitySimilarityMeasureImpl> clone() = 0;
  69 };
  70 DensitySimilarityMeasureImpl::~DensitySimilarityMeasureImpl() = default;
  71
  72 namespace
  73 {
  74
  75 /*! \internal \brief
  76  * Implementation for DensitySimilarityInnerProduct.
  77  * The similarity measure itself is documented in DensitySimilarityMeasureMethod::innerProduct.
  78  */
  79 class DensitySimilarityInnerProduct final : public DensitySimilarityMeasureImpl
  80 {
  81     public:
  82         //! Construct similarity measure by setting the reference density
  83         DensitySimilarityInnerProduct(density referenceDensity);
  84         //! The gradient for the inner product similarity measure is the reference density divided by the number of voxels
  85         density gradient(density comparedDensity) override;
  86         //! Clone this
  87         std::unique_ptr<DensitySimilarityMeasureImpl> clone() override;
  88         //! The similarity between reference density and compared density
  89         float similarity(density comparedDensity) override;
  90     private:
  91         //! A view on the reference density
  92         const density referenceDensity_;
  93         //! Stores the gradient of the similarity measure in memory
  94         MultiDimArray<std::vector<float>, dynamicExtents3D> gradient_;
  95 };
  96
  97 DensitySimilarityInnerProduct::DensitySimilarityInnerProduct(density referenceDensity) :
  98     referenceDensity_ {referenceDensity },
  99 gradient_ {
 100     referenceDensity.extent(XX), referenceDensity.extent(YY), referenceDensity.extent(ZZ)
 101 }
 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_),
 107                    [numVoxels](float x){return x/numVoxels; });
 108 }
 109
 110 float DensitySimilarityInnerProduct::similarity(density comparedDensity)
 111 {
 112     if (comparedDensity.extents() != referenceDensity_.extents())
 113     {
 114         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 115     }
 116     /* the similarity measure uses the gradient instead of the reference,
 117      * here, because it is the reference density divided by the number of voxels */
 118     return std::inner_product(begin(gradient_), end(gradient_), begin(comparedDensity), 0.);
 119 }
 120
 121 DensitySimilarityMeasure::density DensitySimilarityInnerProduct::gradient(density comparedDensity)
 122 {
 123     /* even though the gradient density does not depend on the compad density,
 124      * still checking the extents to make sure we're consistent */
 125     if (comparedDensity.extents() != referenceDensity_.extents())
 126     {
 127         GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
 128     }
 129
 130     return gradient_.asConstView();
 131 }
 132
 133 std::unique_ptr<DensitySimilarityMeasureImpl> DensitySimilarityInnerProduct::clone()
 134 {
 135     return std::make_unique<DensitySimilarityInnerProduct>(referenceDensity_);
 136 }
 137
 138 }   // namespace
 139
 140
 141 DensitySimilarityMeasure::DensitySimilarityMeasure(DensitySimilarityMeasureMethod method, density referenceDensity)
 142 {
 143     // chose the implementation depending on the method of density comparison
 144     switch (method)
 145     {
 146         case DensitySimilarityMeasureMethod::innerProduct:
 147             impl_ = std::make_unique<DensitySimilarityInnerProduct>(referenceDensity);
 148             break;
 149
 150         default:
 151             break;
 152     }
 153 }
 154
 155 DensitySimilarityMeasure::density DensitySimilarityMeasure::gradient(density comparedDensity)
 156 {
 157     return impl_->gradient(comparedDensity);
 158 }
 159
 160 float DensitySimilarityMeasure::similarity(density comparedDensity)
 161 {
 162     return impl_->similarity(comparedDensity);
 163 }
 164
 165 DensitySimilarityMeasure::~DensitySimilarityMeasure()  = default;
 166
 167 DensitySimilarityMeasure::DensitySimilarityMeasure(const DensitySimilarityMeasure &other)
 168     : impl_(other.impl_->clone())
 169 {
 170 }
 171
 172 DensitySimilarityMeasure &DensitySimilarityMeasure::operator=(const DensitySimilarityMeasure &other)
 173 {
 174     impl_ = other.impl_->clone();
 175     return *this;
 176 }
 177
 178 DensitySimilarityMeasure::DensitySimilarityMeasure(DensitySimilarityMeasure &&) noexcept = default;
 179
 180 DensitySimilarityMeasure &DensitySimilarityMeasure::operator=(DensitySimilarityMeasure &&) noexcept = default;
 181
 182 } // namespace gmx