Compare trajectories to reference data
authorPascal Merz <pascal.merz@me.com>
Sun, 27 Sep 2020 11:23:20 +0000 (11:23 +0000)
committerArtem Zhmurov <zhmurov@gmail.com>
Sun, 27 Sep 2020 11:23:20 +0000 (11:23 +0000)
This introduces functionality to compare trajectories (x/v/f/box)
to reference data. It also allows to define the maximum number of
frames which should be checked. The latter allows to have a single
reference data, but decide based on the simulation setup (e.g.
parallelization scheme) how many frames are expected to match
before numerical divergence gets too large. To allow this, the
reference data functionality needs to give the option to ignore
unused reference frames.

src/programs/mdrun/tests/trajectorycomparison.cpp
src/programs/mdrun/tests/trajectorycomparison.h
src/testutils/refdata.cpp
src/testutils/refdata.h

index 16a8a0980c40ec534efc29944221e65db4e76081..72aea1dcef002cf2ed413ea73111962854be1ad4 100644 (file)
@@ -38,6 +38,7 @@
  * produced by mdrun.
  *
  * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \author Pascal Merz <pascal.merz@me.com>
  * \ingroup module_mdrun_integration_tests
  */
 #include "gmxpre.h"
 #include "trajectorycomparison.h"
 
 #include <gmock/gmock.h>
+#include <gtest/gtest-spi.h>
 
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/trajectory/trajectoryframe.h"
+#include "gromacs/utility/strconvert.h"
 
+#include "testutils/refdata.h"
 #include "testutils/testasserts.h"
 #include "testutils/testmatchers.h"
 
+#include "trajectoryreader.h"
+
 namespace gmx
 {
 namespace test
@@ -268,5 +274,161 @@ void TrajectoryComparison::operator()(const TrajectoryFrame& reference, const Tr
     numComparedFrames_++;
 }
 
+/*! \brief Return whether the \c comparisonConditions and emptiness of
+ * the test frame means that a comparison should be attempted.
+ *
+ * This allows the framework to determine whether it is an error if a
+ * comparison cannot be made. */
+static bool shouldDoComparison(ArrayRef<const RVec> values, ComparisonConditions comparisonConditions)
+{
+    if (comparisonConditions == ComparisonConditions::NoComparison)
+    {
+        return false;
+    }
+    if (values.empty())
+    {
+        // Empty array ref is only acceptable if comparison is not required
+        if (comparisonConditions != ComparisonConditions::CompareIfTestFound
+            && comparisonConditions != ComparisonConditions::CompareIfBothFound)
+        {
+            ADD_FAILURE() << "Test frame lacked quantity for required comparison";
+        }
+        return false;
+    }
+    return true;
+}
+
+/*! \brief Compares the positions from \c frame to reference data
+ * according to the \c matchSettings and \c tolerance.
+ */
+static void checkPositionsAgainstReference(const TrajectoryFrame&              frame,
+                                           const TrajectoryFrameMatchSettings& matchSettings,
+                                           const TrajectoryTolerances&         tolerance,
+                                           TestReferenceChecker*               checker)
+{
+    SCOPED_TRACE("Comparing positions");
+    if (shouldDoComparison(frame.x(), matchSettings.coordinatesComparison))
+    {
+        auto positions = frame.x();
+        if (frame.hasBox() && (matchSettings.handlePbcIfPossible || matchSettings.requirePbcHandling))
+        {
+            positions = putAtomsInBox(frame);
+        }
+        if (!frame.hasBox() && matchSettings.requirePbcHandling)
+        {
+            ADD_FAILURE() << "Comparing positions required PBC handling, "
+                             "but the test frame did not have a box";
+            return;
+        }
+        checker->setDefaultTolerance(tolerance.coordinates);
+        std::string id = frame.frameName() + " X";
+        checker->checkSequence(std::begin(positions), std::end(positions), id.c_str());
+    }
+}
+
+/*! \brief Compares the velocities from \c frame to reference data
+ * according to the \c matchSettings and \c tolerance.
+ */
+static void checkVelocitiesAgainstReference(const TrajectoryFrame&              frame,
+                                            const TrajectoryFrameMatchSettings& matchSettings,
+                                            const TrajectoryTolerances&         tolerance,
+                                            TestReferenceChecker*               checker)
+{
+    SCOPED_TRACE("Comparing velocities");
+    if (shouldDoComparison(frame.v(), matchSettings.velocitiesComparison))
+    {
+        checker->setDefaultTolerance(tolerance.velocities);
+        std::string id = frame.frameName() + " V";
+        checker->checkSequence(std::begin(frame.v()), std::end(frame.v()), id.c_str());
+    }
+}
+
+/*! \brief Compares the forces from \c frame to reference data
+ * according to the \c matchSettings and \c tolerance.
+ */
+static void checkForcesAgainstReference(const TrajectoryFrame&              frame,
+                                        const TrajectoryFrameMatchSettings& matchSettings,
+                                        const TrajectoryTolerances&         tolerance,
+                                        TestReferenceChecker*               checker)
+{
+    SCOPED_TRACE("Comparing forces");
+    if (shouldDoComparison(frame.f(), matchSettings.forcesComparison))
+    {
+        checker->setDefaultTolerance(tolerance.forces);
+        std::string id = frame.frameName() + " F";
+        checker->checkSequence(std::begin(frame.f()), std::end(frame.f()), id.c_str());
+    }
+}
+
+/*! \brief Compares the box from \c frame to reference data
+ * according to the \c matchSettings and \c tolerance.
+ */
+static void checkBoxAgainstReference(const TrajectoryFrame&              frame,
+                                     const TrajectoryFrameMatchSettings& matchSettings,
+                                     const TrajectoryTolerances&         tolerance,
+                                     TestReferenceChecker*               checker)
+{
+    SCOPED_TRACE("Comparing box");
+    if (!matchSettings.mustCompareBox)
+    {
+        return;
+    }
+    if (!frame.hasBox())
+    {
+        ADD_FAILURE() << "Comparing the box was required, "
+                         "but the test frame did not have one";
+        return;
+    }
+    checker->setDefaultTolerance(tolerance.box);
+    // Do the comparison
+    for (int d = 0; d < DIM; ++d)
+    {
+        for (int dd = 0; dd < DIM; ++dd)
+        {
+            checker->checkReal(
+                    frame.box()[d][dd],
+                    (frame.frameName() + " Box[" + toString(d) + "][" + toString(dd) + "]").c_str());
+        }
+    }
+}
+
+void TrajectoryComparison::operator()(const TrajectoryFrame& frame, TestReferenceChecker* checker) const
+{
+    SCOPED_TRACE("Comparing trajectory frame " + frame.frameName() + " against reference data");
+    checkBoxAgainstReference(frame, matchSettings_, tolerances_, checker);
+    checkPositionsAgainstReference(frame, matchSettings_, tolerances_, checker);
+    checkVelocitiesAgainstReference(frame, matchSettings_, tolerances_, checker);
+    checkForcesAgainstReference(frame, matchSettings_, tolerances_, checker);
+}
+
+void checkTrajectoryAgainstReferenceData(const std::string&          trajectoryFilename,
+                                         const TrajectoryComparison& trajectoryComparison,
+                                         TestReferenceChecker*       checker)
+{
+    checkTrajectoryAgainstReferenceData(trajectoryFilename, trajectoryComparison, checker,
+                                        MaxNumFrames::compareAllFrames());
+}
+
+void checkTrajectoryAgainstReferenceData(const std::string&          trajectoryFilename,
+                                         const TrajectoryComparison& trajectoryComparison,
+                                         TestReferenceChecker*       checker,
+                                         MaxNumFrames                maxNumFrames)
+{
+    ;
+    TrajectoryFrameReader reader(trajectoryFilename);
+    unsigned int          counter = 0;
+    for (counter = 0; counter < maxNumFrames && reader.readNextFrame(); counter++)
+    {
+        auto frame = reader.frame();
+        trajectoryComparison(frame, checker);
+    }
+
+    if (counter == maxNumFrames && reader.readNextFrame())
+    {
+        // There would have been at least one more frame!
+        checker->disableUnusedEntriesCheck();
+    }
+}
+
 } // namespace test
 } // namespace gmx
index 9e2847f7388589046f595aef4d311df1620b3d7e..d2cc51495f0547c0b92f144478ee3a20f4116b09 100644 (file)
@@ -38,6 +38,7 @@
  * produced by mdrun.
  *
  * \author Mark Abraham <mark.j.abraham@gmail.com>
+ * \author Pascal Merz <pascal.merz@me.com>
  * \ingroup module_mdrun_integration_tests
  */
 #ifndef GMX_PROGRAMS_MDRUN_TESTS_TRAJECTORYCOMPARISON_H
@@ -55,6 +56,8 @@ class TrajectoryFrame;
 namespace test
 {
 
+class TestReferenceChecker;
+
 /*! \internal
  * \brief Helper struct for testing different trajectory components with different tolerances. */
 struct TrajectoryTolerances
@@ -104,10 +107,10 @@ struct TrajectoryFrameMatchSettings
 };
 
 /*! \internal
- * \brief Function object to compare the fields of the two frames for
- * equality given the \c matchSettings_ and \c tolerances_.
+ * \brief Function object to compare the fields of two frames vs each others or one
+ * frame vs reference data for equality given the \c matchSettings_ and \c tolerances_.
  *
- * The two frames are required to have valid and matching values for
+ * If using two frames, they are required to have valid and matching values for
  * time and step. According to \c matchSettings_, box, position coordinates,
  * velocities and/or forces will be compared between frames, using the
  * \c tolerances_. Comparisons will only occur when both frames have
@@ -116,7 +119,16 @@ struct TrajectoryFrameMatchSettings
  * GoogleTest expectation failure will be given. If a comparison is
  * required by \c matchSettings_ but cannot be done because either (or
  * both) frames lack the requisite data, descriptive expectation
- * failures will be given. */
+ * failures will be given.
+ *
+ * When comparing a frame to reference data, according to \c matchSettings_,
+ * box, position coordinates, velocities and/or forces will be compared to
+ * reference data, using the \c tolerances_. If a comparison fails, a
+ * GoogleTest expectation failure will be given. If a comparison is
+ * required by \c matchSettings_ but cannot be done because the
+ * frame lacks the requisite data, descriptive expectation
+ * failures will be given.
+ */
 class TrajectoryComparison
 {
 public:
@@ -128,6 +140,9 @@ public:
     /*! \brief Compare reference with test given the \c
      * matchSettings_ within \c tolerances_ */
     void operator()(const TrajectoryFrame& reference, const TrajectoryFrame& test) const;
+    /*! \brief Compare frame to reference given the \c
+     * matchSettings_ within \c tolerances_ */
+    void operator()(const TrajectoryFrame& frame, TestReferenceChecker* checker) const;
 
 private:
     //! Specifies expected behavior in comparisons
@@ -143,6 +158,27 @@ private:
     mutable unsigned int numComparedFrames_ = 0;
 };
 
+/*! \internal
+ * \brief Helper function comparing a trajectory to reference.
+ *
+ * This opens a trajectory file, loops over its frames, and uses the
+ * \c trajectoryComparison and \c checker arguments to check each frame
+ * to reference data.
+ *
+ * Optionally, the max number of trajectory frames to be compared can
+ * be specified to allow to only test the first few frames of a
+ * simulation, e.g. to avoid failures due to numerical divergence.
+ */
+//! \{
+void checkTrajectoryAgainstReferenceData(const std::string&          trajectoryFilename,
+                                         const TrajectoryComparison& trajectoryComparison,
+                                         TestReferenceChecker*       checker);
+void checkTrajectoryAgainstReferenceData(const std::string&          trajectoryFilename,
+                                         const TrajectoryComparison& trajectoryComparison,
+                                         TestReferenceChecker*       checker,
+                                         MaxNumFrames                maxNumFrames);
+//! \}
+
 } // namespace test
 } // namespace gmx
 
index e2dc6ed6fd027e117de6343e9620d58ef181c1cd..8c92d6754f5e4256d2bdd0ab53dee71690222523 100644 (file)
@@ -53,6 +53,7 @@
 
 #include <gtest/gtest.h>
 
+#include "gromacs/math/vectypes.h"
 #include "gromacs/options/basicoptions.h"
 #include "gromacs/options/ioptionscontainer.h"
 #include "gromacs/utility/any.h"
@@ -953,7 +954,7 @@ void TestReferenceChecker::checkRealFromString(const std::string& value, const c
 }
 
 
-void TestReferenceChecker::checkVector(const int value[3], const char* id)
+void TestReferenceChecker::checkVector(const BasicVector<int>& value, const char* id)
 {
     TestReferenceChecker compound(checkCompound(Impl::cVectorType, id));
     compound.checkInteger(value[0], "X");
@@ -962,7 +963,7 @@ void TestReferenceChecker::checkVector(const int value[3], const char* id)
 }
 
 
-void TestReferenceChecker::checkVector(const float value[3], const char* id)
+void TestReferenceChecker::checkVector(const BasicVector<float>& value, const char* id)
 {
     TestReferenceChecker compound(checkCompound(Impl::cVectorType, id));
     compound.checkReal(value[0], "X");
@@ -971,7 +972,7 @@ void TestReferenceChecker::checkVector(const float value[3], const char* id)
 }
 
 
-void TestReferenceChecker::checkVector(const double value[3], const char* id)
+void TestReferenceChecker::checkVector(const BasicVector<double>& value, const char* id)
 {
     TestReferenceChecker compound(checkCompound(Impl::cVectorType, id));
     compound.checkReal(value[0], "X");
@@ -980,6 +981,24 @@ void TestReferenceChecker::checkVector(const double value[3], const char* id)
 }
 
 
+void TestReferenceChecker::checkVector(const int value[3], const char* id)
+{
+    checkVector(BasicVector<int>(value), id);
+}
+
+
+void TestReferenceChecker::checkVector(const float value[3], const char* id)
+{
+    checkVector(BasicVector<float>(value), id);
+}
+
+
+void TestReferenceChecker::checkVector(const double value[3], const char* id)
+{
+    checkVector(BasicVector<double>(value), id);
+}
+
+
 void TestReferenceChecker::checkAny(const Any& any, const char* id)
 {
     if (any.isType<bool>())
index 21b8a2c96019280e4ee6403fb205ae7385a3f815..bc5323acfcc03816152c19a4d38a871a2e7b8ddb 100644 (file)
@@ -61,6 +61,9 @@ class KeyValueTreeObject;
 class KeyValueTreeValue;
 class Any;
 
+template<typename ValueType>
+class BasicVector;
+
 namespace test
 {
 
@@ -168,6 +171,7 @@ class TestReferenceDataImpl;
        compound.checkInteger(function2ToTest(3), "ValueWith3");
        compound.checkInteger(function2ToTest(5), "ValueWith5");
        checker.checkInteger(functionToTest(4), "ValueWith4");
+       checker.checkVector(functionProducingRVec(), "Describe The RVec");
    }
    } // namespace test
    } // namespace gmx
@@ -391,6 +395,12 @@ public:
     void checkVector(const float value[3], const char* id);
     //! Check a vector of three double-precision floating point values.
     void checkVector(const double value[3], const char* id);
+    //! Check a BasicVector of ints, ie. IVec
+    void checkVector(const BasicVector<int>& value, const char* id);
+    //! Check a BasicVector of floats, ie. RVec
+    void checkVector(const BasicVector<float>& value, const char* id);
+    //! Check a BasicVector of doubles, ie. DVec
+    void checkVector(const BasicVector<double>& value, const char* id);
     //! Check a single floating-point value from a string.
     void checkRealFromString(const std::string& value, const char* id);
     //! Checks a any value that contains a supported simple type.
@@ -458,6 +468,12 @@ public:
     void checkValue(const float value[3], const char* id) { checkVector(value, id); }
     //! Check a vector of three double-precision floating point values.
     void checkValue(const double value[3], const char* id) { checkVector(value, id); }
+    //! Check a BasicVector of integer values, ie. IVec.
+    void checkValue(const BasicVector<int>& value, const char* id) { checkVector(value, id); }
+    //! Check a BasicVector of float values, ie. RVec.
+    void checkValue(const BasicVector<float>& value, const char* id) { checkVector(value, id); }
+    //! Check a BasicVector of double values, ie. DVec.
+    void checkValue(const BasicVector<double>& value, const char* id) { checkVector(value, id); }
     //! Check a generic key-value tree value.
     void checkValue(const KeyValueTreeValue& value, const char* id)
     {