There are multiple valid similarity measures between the reference density and
the simulated density, each motivated by the experimental source of the
-reference density data. For the density-guided simulations in gromacs, two
-measures are provided, the inner product of the simulated density,
+reference density data. For the density-guided simulations in |Gromacs|, the following
+measures are provided:
+
+The inner product of the simulated density,
.. math:: S_{\mathrm{inner-product}}[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}] =
- \frac{1}{N_\mathrm{voxel}}\sum_{v=1}^{N_\mathrm{voxel}} \rho^{\mathrm{ref}}_v \rho^{\mathrm{sim}}_v\,\mathrm{,}
+ \frac{1}{N_\mathrm{voxel}}\sum_{v=1}^{N_\mathrm{voxel}} \rho^{\mathrm{ref}}_v \rho^{\mathrm{sim}}_v\,\mathrm{.}
-as well as the negative relative entropy between two densities
+The negative relative entropy between two densities,
.. math:: S_{\mathrm{relative-entropy}}[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}] =
\sum_{v=1, \rho^{\mathrm{ref}}>0, \rho^{\mathrm{sim}}>0}^{N_\mathrm{voxel}} \rho^\mathrm{ref} [\log(\rho^\mathrm{sim}_v)-\log(\rho^\mathrm{ref}_v)]\,\mathrm{.}
+
+The cross correlation between two densities,
+
+.. math:: S_{\mathrm{cross-correlation}}[\rho^{\mathrm{ref}},\rho^{\mathrm{sim}}] =
+ \frac{\sum_{v}\left((\rho_v^{\mathrm{ref}} - \bar{\rho}^{\mathrm{ref}})(\rho_v^{\mathrm{sim}} - \bar{\rho}^{\mathrm{sim}})\right)}
+ {\sqrt{\sum_v(\rho_v^{\mathrm{ref}} - \bar{\rho}^{\mathrm{ref}})^2 \sum_v(\rho_v^{\mathrm{sim}} - \bar{\rho}^{\mathrm{sim}})^2}}\mathrm{.}
+
Declaring regions to fit
{
addDensityFittingMdpOutputValue(builder, groupString_, c_groupTag_);
- addDensityFittingMdpOutputValueComment(builder, "; Similarity measure between densities: inner-product, or relative-entropy", c_similarityMeasureTag_);
+ addDensityFittingMdpOutputValueComment(builder, "; Similarity measure between densities: inner-product, relative-entropy, or cross-correlation", c_similarityMeasureTag_);
addDensityFittingMdpOutputValue<std::string>(builder,
c_densitySimilarityMeasureMethodNames[parameters_.similarityMeasureMethod_],
c_similarityMeasureTag_);
writeKeyValueTreeAsMdp(&writer, builder.build());
}
stream.close();
- std::string expected
+ std::string expectedString
= {
"\n"
"; Density guided simulation\n"
"density-guided-simulation-active = true\n"
"density-guided-simulation-group = protein\n"
- "; Similarity measure between densities: inner-product, or relative-entropy\n"
+ "; Similarity measure between densities: inner-product, relative-entropy, or cross-correlation\n"
"density-guided-simulation-similarity-measure = inner-product\n"
"; Atom amplitude for spreading onto grid: unity, mass, or charge\n"
"density-guided-simulation-atom-spreading-weight = unity\n"
"density-guided-simulation-normalize-densities = true\n"
};
- EXPECT_EQ(expected, stream.toString());
+ EXPECT_EQ(expectedString, stream.toString());
}
return std::make_unique<DensitySimilarityRelativeEntropy>(referenceDensity_);
}
+/****************** Cross Correlation *****************************************/
+
+//! Helper values for evaluating the cross correlation
+struct CrossCorrelationEvaluationHelperValues
+{
+ //! The mean of the reference density
+ real meanReference = 0;
+ //! The mean of the compared density
+ real meanComparison = 0;
+ //! The sum of the squared reference density voxel values
+ real referenceSquaredSum = 0;
+ //! The sum of the squared compared density voxel values
+ real comparisonSquaredSum = 0;
+ //! The covariance of the refernce and the compared density
+ real covariance = 0;
+};
+
+/*! \brief Calculate helper values for the cross-correlation.
+
+ * Enables numerically stable single-pass cross-correlation evaluation algorithm
+ * as described in Bennett, J., Grout, R. , Pebay, P., Roe D., Thompson D.
+ * "Numerically Stable, Single-Pass, Parallel Statistics Algorithms"
+ * and implemented in boost's correlation coefficient
+ */
+CrossCorrelationEvaluationHelperValues
+evaluateHelperValues(DensitySimilarityMeasure::density reference,
+ DensitySimilarityMeasure::density compared)
+{
+ CrossCorrelationEvaluationHelperValues helperValues;
+
+ index i = 0;
+
+ auto referenceIterator = begin(reference);
+ for (const real comp : compared)
+ {
+ const real refHelper = *referenceIterator - helperValues.meanReference;
+ const real comparisonHelper = comp - helperValues.meanComparison;
+ helperValues.referenceSquaredSum += (i * square(refHelper)) / (i + 1);
+ helperValues.comparisonSquaredSum += (i * square(comparisonHelper)) / (i + 1);
+ helperValues.covariance += i * refHelper * comparisonHelper / (i + 1);
+ helperValues.meanReference += refHelper / (i + 1);
+ helperValues.meanComparison += comparisonHelper / (i + 1);
+
+ ++referenceIterator;
+ ++i;
+ }
+
+ return helperValues;
+}
+
+//! Calculate a single cross correlation gradient entry at a voxel.
+class CrossCorrelationGradientAtVoxel
+{
+ public:
+
+ //! Set up the gradident calculation with pre-computed values
+ CrossCorrelationGradientAtVoxel(const CrossCorrelationEvaluationHelperValues &preComputed)
+ : prefactor_(evaluatePrefactor(preComputed.comparisonSquaredSum, preComputed.referenceSquaredSum)),
+ comparisonPrefactor_(preComputed.covariance / preComputed.comparisonSquaredSum),
+ meanReference_(preComputed.meanReference), meanComparison_(preComputed.meanComparison)
+ {
+ }
+ //! Evaluate the cross correlation gradient at a voxel
+ real operator()(real reference, real comparison)
+ {
+ return prefactor_ * (reference - meanReference_ - comparisonPrefactor_ * (comparison-meanComparison_));
+ }
+ private:
+ real evaluatePrefactor(real comparisonSquaredSum, real referenceSquaredSum)
+ {
+ GMX_ASSERT(comparisonSquaredSum > 0,
+ "Squared sum of comparison values needs to be larger than zero.");
+ GMX_ASSERT(referenceSquaredSum > 0,
+ "Squared sum of reference values needs to be larger than zero.");
+ return 1.0 / (sqrt(comparisonSquaredSum) * sqrt(referenceSquaredSum));
+ }
+ const real prefactor_;
+ const real comparisonPrefactor_;
+ const real meanReference_;
+ const real meanComparison_;
+};
+
+/*! \internal
+ * \brief Implementation for DensitySimilarityCrossCorrelation.
+ *
+ * The similarity measure itself is documented in DensitySimilarityMeasureMethod::crossCorrelation.
+ */
+class DensitySimilarityCrossCorrelation final : public DensitySimilarityMeasureImpl
+{
+ public:
+ //! Construct similarity measure by setting the reference density
+ DensitySimilarityCrossCorrelation(density referenceDensity);
+ //! The gradient for the cross correlation similarity measure
+ density gradient(density comparedDensity) override;
+ //! Clone this
+ std::unique_ptr<DensitySimilarityMeasureImpl> clone() override;
+ //! The similarity between reference density and compared density
+ real similarity(density comparedDensity) override;
+
+ private:
+ //! A view on the reference density
+ const density referenceDensity_;
+ //! Stores the gradient of the similarity measure in memory
+ MultiDimArray<std::vector<float>, dynamicExtents3D> gradient_;
+};
+
+DensitySimilarityCrossCorrelation::DensitySimilarityCrossCorrelation(density referenceDensity)
+ : referenceDensity_ { referenceDensity }, gradient_(referenceDensity.extents())
+{
+}
+
+real DensitySimilarityCrossCorrelation::similarity(density comparedDensity)
+{
+ if (comparedDensity.extents() != referenceDensity_.extents())
+ {
+ GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
+ }
+
+ CrossCorrelationEvaluationHelperValues helperValues = evaluateHelperValues(referenceDensity_, comparedDensity);
+
+ if ((helperValues.referenceSquaredSum == 0) || (helperValues.comparisonSquaredSum == 0))
+ {
+ return 0;
+ }
+ real covarianceSqrt = sqrt(helperValues.covariance);
+ return (covarianceSqrt / sqrt(helperValues.referenceSquaredSum)) * (covarianceSqrt / sqrt(helperValues.comparisonSquaredSum));
+}
+
+DensitySimilarityMeasure::density DensitySimilarityCrossCorrelation::gradient(density comparedDensity)
+{
+ if (comparedDensity.extents() != referenceDensity_.extents())
+ {
+ GMX_THROW(RangeError("Reference density and compared density need to have same extents."));
+ }
+
+ CrossCorrelationEvaluationHelperValues helperValues = evaluateHelperValues(referenceDensity_, comparedDensity);
+
+ std::transform(begin(referenceDensity_), end(referenceDensity_),
+ begin(comparedDensity),
+ begin(gradient_),
+ CrossCorrelationGradientAtVoxel(helperValues));
+
+ return gradient_.asConstView();
+}
+
+std::unique_ptr<DensitySimilarityMeasureImpl> DensitySimilarityCrossCorrelation::clone()
+{
+ return std::make_unique<DensitySimilarityCrossCorrelation>(referenceDensity_);
+}
+
+
} // namespace
DensitySimilarityMeasure::DensitySimilarityMeasure(DensitySimilarityMeasureMethod method, density referenceDensity)
{
// chose the implementation depending on the method of density comparison
+ // throw an error if the method is not known
switch (method)
{
case DensitySimilarityMeasureMethod::innerProduct:
case DensitySimilarityMeasureMethod::relativeEntropy:
impl_ = std::make_unique<DensitySimilarityRelativeEntropy>(referenceDensity);
break;
- default:
+ case DensitySimilarityMeasureMethod::crossCorrelation:
+ impl_ = std::make_unique<DensitySimilarityCrossCorrelation>(referenceDensity);
break;
+ default:
+ GMX_THROW(NotImplementedError("Similarity measure not implemented."));
}
}
*
* \f[
* \mathrm{Similarity}(\rho_{\mathrm{r}},\rho_{\mathrm{c}}) =
- * \frac{1}{N_\mathrm{voxel}}/\sum_{v=1}^{N_\mathrm{voxel}} \rho^v_{\mathrm{r}} \rho^v_{\mathrm{c}}
+ * \frac{1}{N_\mathrm{voxel}}/\sum_{v=1}^{N_\mathrm{voxel}} \rho^v_{\mathrm{r}}
+ * \rho^v_{\mathrm{c}}
* \f]
*/
innerProduct,
- /*! \brief Measure similarity between densities, using the relative entropy between them.
+
+ /*! \brief Measure similarity between densities by negative relative entropy.
* \note Voxels with negative values are ignored.
*
* \f[
* \f]
*/
relativeEntropy,
+
+ /*! \brief Measure similarity between densities by cross correlation.
+ *
+ * \f[
+ * \mathrm{Similarity}(\rho_{\mathrm{r}},\rho_{\mathrm{c}}) =
+ * \frac{\sum_{v}\left((\rho_{\mathrm{r}} - \bar{\rho}_{\mathrm{r}})(\rho_{\mathrm{c}} - \bar{\rho}_{\mathrm{c}})\right)}
+ * {\sqrt{\sum_v(\rho_{\mathrm{r}} - \bar{\rho}_{\mathrm{r}})^2 \sum_v (\rho_{\mathrm{c}} - \bar{\rho}_{\mathrm{c}})^2}}
+ * \f]
+ */
+ crossCorrelation,
Count,
};
//! Name the methods that may be used to evaluate similarity between densities
const EnumerationArray<DensitySimilarityMeasureMethod, const char *const>
-c_densitySimilarityMeasureMethodNames = {{ "inner-product", "relative-entropy"}};
+c_densitySimilarityMeasureMethodNames = {{ "inner-product", "relative-entropy", "cross-correlation"}};
/* Forward declaration of implementation class outside class to allow
* choose implementation class during construction of the DensitySimilarityMeasure*/
/*! \brief Chose comparison method and set reference density.
* \param[in] method defines how densities are compared to one another
* \param[in] referenceDensity
+ * \throws NotImplementedError if method is not known
*/
DensitySimilarityMeasure(DensitySimilarityMeasureMethod method, density referenceDensity);
~DensitySimilarityMeasure();
MultiDimArray<std::vector<float>, dynamicExtents3D> referenceDensity(3, 3, 3);
std::iota(begin(referenceDensity), end(referenceDensity), -1);
- DensitySimilarityMeasure measure(DensitySimilarityMeasureMethod::innerProduct, referenceDensity.asConstView());
+ DensitySimilarityMeasure measure(DensitySimilarityMeasureMethod::relativeEntropy, referenceDensity.asConstView());
MultiDimArray<std::vector<float>, dynamicExtents3D> comparedDensity(3, 3, 3);
std::iota(begin(comparedDensity), end(comparedDensity), -2);
checker.checkSequence(gradientView.begin(), gradientView.end(), "relative-entropy-gradient");
}
+TEST(DensitySimilarityTest, CrossCorrelationIsOne)
+{
+ MultiDimArray<std::vector<float>, dynamicExtents3D> referenceDensity(100, 100, 100);
+ std::iota(begin(referenceDensity), end(referenceDensity), 10000);
+
+ DensitySimilarityMeasure measure(DensitySimilarityMeasureMethod::crossCorrelation,
+ referenceDensity.asConstView());
+
+ MultiDimArray<std::vector<float>, dynamicExtents3D> comparedDensity(100, 100, 100);
+ std::iota(begin(comparedDensity), end(comparedDensity), -10000);
+
+ const real expectedSimilarity = 1;
+ EXPECT_REAL_EQ(expectedSimilarity, measure.similarity(comparedDensity.asConstView()));
+}
+
+TEST(DensitySimilarityTest, CrossCorrelationGradientIsZeroWhenCorrelated)
+{
+ MultiDimArray<std::vector<float>, dynamicExtents3D> referenceDensity(30, 30, 30);
+ std::iota(begin(referenceDensity), end(referenceDensity), -1);
+
+ DensitySimilarityMeasure measure(DensitySimilarityMeasureMethod::crossCorrelation,
+ referenceDensity.asConstView());
+
+ MultiDimArray<std::vector<float>, dynamicExtents3D> comparedDensity(30, 30, 30);
+ std::iota(begin(comparedDensity), end(comparedDensity), -2);
+
+ // Need this conversion to ArrayRef, because Pointwise requires size()
+ // member function not provided by basic_mdspan
+ const basic_mdspan<const float, dynamicExtents3D> gradient
+ = measure.gradient(comparedDensity.asConstView());
+ ArrayRef<const float> gradientView(gradient.data(),
+ gradient.data() + gradient.mapping().required_span_size());
+
+ std::array<float, 27000> expectedSimilarityGradient = {};
+
+ EXPECT_THAT(expectedSimilarityGradient, Pointwise(FloatEq(defaultFloatTolerance()), gradientView));
+}
+
+TEST(DensitySimilarityTest, CrossCorrelationGradientIsCorrect)
+{
+ MultiDimArray<std::vector<float>, dynamicExtents3D> referenceDensity(3, 3, 3);
+ std::iota(begin(referenceDensity), end(referenceDensity), -1);
+
+ DensitySimilarityMeasure measure(DensitySimilarityMeasureMethod::crossCorrelation,
+ referenceDensity.asConstView());
+
+ MultiDimArray<std::vector<float>, dynamicExtents3D> comparedDensity(3, 3, 3);
+ std::iota(begin(comparedDensity), end(comparedDensity), -2);
+
+ // some non-linear transformation, so that we break the correlation
+ for (float &valueToCompare : comparedDensity)
+ {
+ valueToCompare *= valueToCompare;
+ }
+
+ // Need this conversion to ArrayRef, because Pointwise requires size()
+ // member function not provided by basic_mdspan
+ const basic_mdspan<const float, dynamicExtents3D> gradient
+ = measure.gradient(comparedDensity.asConstView());
+ ArrayRef<const float> gradientView(gradient.data(),
+ gradient.data() + gradient.mapping().required_span_size());
+
+ TestReferenceData refData;
+ TestReferenceChecker checker(refData.rootChecker());
+ checker.setDefaultTolerance(defaultFloatTolerance());
+ checker.checkSequence(gradientView.begin(), gradientView.end(), "cross-correlation-gradient");
+}
+
} // namespace test
} // namespace gmx
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Sequence Name="cross-correlation-gradient">
+ <Int Name="Length">27</Int>
+ <Real>-0.00014969852</Real>
+ <Real>-0.00011995764</Real>
+ <Real>-9.2403832e-05</Real>
+ <Real>-6.7037101e-05</Real>
+ <Real>-4.3857435e-05</Real>
+ <Real>-2.2864861e-05</Real>
+ <Real>-4.0593445e-06</Real>
+ <Real>1.25591e-05</Real>
+ <Real>2.699047e-05</Real>
+ <Real>3.9234772e-05</Real>
+ <Real>4.9291986e-05</Real>
+ <Real>5.7162139e-05</Real>
+ <Real>6.2845218e-05</Real>
+ <Real>6.6341221e-05</Real>
+ <Real>6.7650144e-05</Real>
+ <Real>6.6772009e-05</Real>
+ <Real>6.3706793e-05</Real>
+ <Real>5.8454505e-05</Real>
+ <Real>5.101514e-05</Real>
+ <Real>4.138871e-05</Real>
+ <Real>2.9575202e-05</Real>
+ <Real>1.557462e-05</Real>
+ <Real>-6.130212e-07</Real>
+ <Real>-1.8987734e-05</Real>
+ <Real>-3.9549544e-05</Real>
+ <Real>-6.2298401e-05</Real>
+ <Real>-8.7234337e-05</Real>
+ </Sequence>
+</ReferenceData>
<ReferenceData>
<Sequence Name="relative-entropy-gradient">
<Int Name="Length">27</Int>
- <Real>-0.037037037</Real>
<Real>0</Real>
- <Real>0.037037037</Real>
- <Real>0.074074075</Real>
- <Real>0.11111111</Real>
- <Real>0.14814815</Real>
- <Real>0.18518518</Real>
- <Real>0.22222222</Real>
- <Real>0.25925925</Real>
- <Real>0.2962963</Real>
- <Real>0.33333334</Real>
- <Real>0.37037036</Real>
- <Real>0.4074074</Real>
- <Real>0.44444445</Real>
- <Real>0.48148149</Real>
- <Real>0.51851851</Real>
- <Real>0.55555558</Real>
- <Real>0.5925926</Real>
- <Real>0.62962961</Real>
- <Real>0.66666669</Real>
- <Real>0.7037037</Real>
- <Real>0.74074072</Real>
- <Real>0.77777779</Real>
- <Real>0.81481481</Real>
- <Real>0.85185188</Real>
- <Real>0.8888889</Real>
- <Real>0.92592591</Real>
+ <Real>0</Real>
+ <Real>0</Real>
+ <Real>2</Real>
+ <Real>1.5</Real>
+ <Real>1.3333334</Real>
+ <Real>1.25</Real>
+ <Real>1.2</Real>
+ <Real>1.1666666</Real>
+ <Real>1.1428572</Real>
+ <Real>1.125</Real>
+ <Real>1.1111112</Real>
+ <Real>1.1</Real>
+ <Real>1.0909091</Real>
+ <Real>1.0833334</Real>
+ <Real>1.0769231</Real>
+ <Real>1.0714285</Real>
+ <Real>1.0666667</Real>
+ <Real>1.0625</Real>
+ <Real>1.0588236</Real>
+ <Real>1.0555556</Real>
+ <Real>1.0526316</Real>
+ <Real>1.05</Real>
+ <Real>1.0476191</Real>
+ <Real>1.0454545</Real>
+ <Real>1.0434783</Real>
+ <Real>1.0416666</Real>
</Sequence>
</ReferenceData>
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff