Updated exponential fitting to make it robust.
[alexxy/gromacs.git] / src / gromacs / correlationfunctions / tests / expfit.cpp
old mode 100755 (executable)
new mode 100644 (file)
index c9f518b..74fb2eb
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2014, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
 #include <gtest/gtest.h>
 
 #include "gromacs/fileio/xvgr.h"
+#include "gromacs/legacyheaders/oenv.h"
 #include "gromacs/utility/smalloc.h"
 
 #include "testutils/refdata.h"
 #include "testutils/testasserts.h"
 #include "testutils/testfilemanager.h"
 
-//! Number of data files for testing.
-#define expTestNrTypes 3
-
 namespace gmx
 {
 
 namespace
 {
 
+class ExpfitData
+{
+    public:
+        int               nrLines_;
+        std::vector<real> x_, y_;
+        real              startTime_, endTime_, dt_;
+};
+
 class ExpfitTest : public ::testing::Test
 {
 
     protected:
-        static int                 nrLines_;
-        static std::vector<real>   values_[expTestNrTypes];
-        static int                 nrColumns_;
-        static std::vector<real>   standardDev_;
-        static real                startTime_;
-        static real                endTime_;
-        static real                timeDeriv_;
-        test::TestReferenceData    refData_;
-        test::TestReferenceChecker checker_;
+        static std::vector<ExpfitData> data_;
+        test::TestReferenceData        refData_;
+        test::TestReferenceChecker     checker_;
         ExpfitTest( )
             : checker_(refData_.rootChecker())
         {
@@ -85,40 +85,34 @@ class ExpfitTest : public ::testing::Test
         // Static initiation, only run once every test.
         static void SetUpTestCase()
         {
-            double   ** tempValues;
-            std::string fileName[expTestNrTypes];
-            fileName[0] = test::TestFileManager::getInputFilePath("testINVEXP.xvg");
-            fileName[1] = test::TestFileManager::getInputFilePath("testPRES.xvg");
-            fileName[2] = test::TestFileManager::getInputFilePath("testEXP.xvg");
-            for (int i = 0; i < expTestNrTypes; i++)
+            double                ** tempValues = NULL;
+            std::vector<std::string> fileName;
+            fileName.push_back(test::TestFileManager::getInputFilePath("testINVEXP.xvg"));
+            fileName.push_back(test::TestFileManager::getInputFilePath("testPRES.xvg"));
+            fileName.push_back(test::TestFileManager::getInputFilePath("testINVEXP79.xvg"));
+            fileName.push_back(test::TestFileManager::getInputFilePath("testERF.xvg"));
+            fileName.push_back(test::TestFileManager::getInputFilePath("testERREST.xvg"));
+            for (std::vector<std::string>::iterator i = fileName.begin(); i < fileName.end(); ++i)
             {
-                const char * name = fileName[i].c_str();
-                // TODO: this assumes all files have the same length.
-                nrLines_     = read_xvg(name, &tempValues, &nrColumns_);
-
-                // Generating standard deviation
-                if (i == 0)
+                const char * name = i->c_str();
+                int          nrColumns;
+                ExpfitData   ed;
+                ed.nrLines_   = read_xvg(name, &tempValues, &nrColumns);
+                ed.dt_        = tempValues[0][1] - tempValues[0][0];
+                ed.startTime_ = tempValues[0][0];
+                ed.endTime_   = tempValues[0][ed.nrLines_-1];
+                for (int j = 0; j  < ed.nrLines_; j++)
                 {
-                    double fac = 1.0/nrLines_;
-                    for (int j = 0; j < nrLines_; j++)
-                    {
-                        standardDev_.push_back(fac);
-                    }
-                    timeDeriv_ = tempValues[0][1] - tempValues[0][0];
-                    startTime_ = tempValues[0][0];
-                    endTime_   = tempValues[0][nrLines_-1];
-                }
-
-                for (int j = 0; j  < nrLines_; j++)
-                {
-                    values_[i].push_back((real)tempValues[1][j]);
+                    ed.x_.push_back((real)tempValues[0][j]);
+                    ed.y_.push_back((real)tempValues[1][j]);
                 }
+                data_.push_back(ed);
 
                 // Free memory that was allocated in read_xvg
-                for (int i = 0; i < nrColumns_; i++)
+                for (int j = 0; j < nrColumns; j++)
                 {
-                    sfree(tempValues[i]);
-                    tempValues[i] = NULL;
+                    sfree(tempValues[j]);
+                    tempValues[j] = NULL;
                 }
                 sfree(tempValues);
                 tempValues = NULL;
@@ -129,13 +123,26 @@ class ExpfitTest : public ::testing::Test
         {
         }
 
-        void test(int type, double result[], double tolerance, int testType)
+        void test(int type, double result[], double tolerance,
+                  unsigned int testType)
         {
-            int     nfitparm = effnNparams(type);
-
-            do_lmfit(nrLines_, &values_[testType][0], &standardDev_[0], timeDeriv_,
-                     NULL, startTime_, endTime_, NULL, false, type, result, 0);
+            int          nfitparm = effnNparams(type);
+            output_env_t oenv;
 
+            if (testType >= data_.size())
+            {
+                GMX_THROW(InvalidInputError("testType out of range"));
+            }
+            output_env_init_default(&oenv);
+            do_lmfit(data_[testType].nrLines_,
+                     &(data_[testType].y_[0]),
+                     NULL,
+                     data_[testType].dt_,
+                     &(data_[testType].x_[0]),
+                     data_[testType].startTime_,
+                     data_[testType].endTime_,
+                     oenv, false, type, result, 0, NULL);
+            output_env_done(oenv);
             checker_.setDefaultTolerance(test::relativeToleranceAsFloatingPoint(1, tolerance));
             checker_.checkSequenceArray(nfitparm, result, "result");
         }
@@ -143,63 +150,55 @@ class ExpfitTest : public ::testing::Test
 
 
 //static var
-int               ExpfitTest::nrLines_;
-//cppcheck-suppress arrayIndexOutOfBounds fixed in 1.68-dev
-std::vector<real> ExpfitTest::values_[expTestNrTypes];
-int               ExpfitTest::nrColumns_;
-std::vector<real> ExpfitTest::standardDev_;
-real              ExpfitTest::startTime_;
-real              ExpfitTest::endTime_;
-real              ExpfitTest::timeDeriv_;
+std::vector<ExpfitData> ExpfitTest::data_;
 
 TEST_F (ExpfitTest, EffnEXP1) {
     double  param[] = {25};
-    test(effnEXP1, param, 1e-6, 0);
+    test(effnEXP1, param, 1e-5, 0);
 }
 
 TEST_F (ExpfitTest, EffnEXP2) {
     double  param[] = {35, 0.5};
-    test(effnEXP2, param, 1e-6, 0);
+    test(effnEXP2, param, 3e-5, 0);
 }
 
-TEST_F (ExpfitTest, EffnEXP3) {
-    double param[] = {45, 0.5, 5};
-    test(effnEXP3, param, 1e-4, 0);
+TEST_F (ExpfitTest, EffnEXPEXP) {
+    double param[] = {5, 0.5, 45};
+    test(effnEXPEXP, param, 1e-2, 0);
 }
 
 TEST_F (ExpfitTest, EffnEXP5) {
     double  param[] = {0.5, 5, 0.5, 50, 0.002};
-    test(effnEXP5, param, 1e-4, 0);
+    test(effnEXP5, param, 1e-2, 2);
 }
 
 TEST_F (ExpfitTest, EffnEXP7) {
-    double  param[] = {0.5, 5, -0.02, 0.5, 0.5, 50, -0.002};
-    test(effnEXP7, param, 1e-4, 0);
+    double  param[] = {0.1, 2, 0.5, 30, 0.3, 50, -0.002};
+    test(effnEXP7, param, 1e-2, 2);
 }
 
 TEST_F (ExpfitTest, EffnEXP9) {
-    double  param[] = {2, 1200, -1, 300, 0.7, 70, 0.5, 6, -0.5};
-    test(effnEXP9, param, 4e-2, 0);
+    double  param[] = {0.4, 5, 0.2, 30, 0.1, 70, 0.2, 200, -0.05};
+    test(effnEXP9, param, 4e-2, 2);
 }
 
 TEST_F (ExpfitTest, EffnERF) {
-    double  param[] = {0.5, 0.5, 0.5, 1};
-    test(effnERF, param, 1e-2, 0);
+    double  param[] = {80, 120, 180, 5};
+    test(effnERF, param, 1e-1, 3);
 }
 
 TEST_F (ExpfitTest, EffnERREST) {
-    double  param[] = {0.5, 0.7, 0.3};
-    test(effnERREST, param, 1e-4, 2);
+    double  param[] = {1, 0.9, 100};
+    test(effnERREST, param, 5e-3, 4);
 }
 
 TEST_F (ExpfitTest, EffnVAC) {
-    double param[] = {0.5, 0.05};
-    test(effnVAC, param, 1e-4, 0);
+    double param[] = {30, 0.0};
+    test(effnVAC, param, 0.05, 0);
 }
 
-TEST_F (ExpfitTest, DISABLED_EffnPRES) {
-    //TODO: This test is prodocues NaNs and INFs. Fix and then reactivate.
-    double param[] = {0, 10, 4, 1, 0.5, 1};
+TEST_F (ExpfitTest, EffnPRES) {
+    double param[] = {0.6, 10, 7, 1, 0.25, 2};
     test(effnPRES, param, 1e-4, 1);
 }