- changed "double" to "long double" for fit const

// enabled better precision
- added tests to check on fitting algo

diff --git a/src/domaintests.cpp b/src/domaintests.cpp
index 192361a2543b2b0b9f872010bba1c9990e3a41a6..c22150783f84968e9bf42c8a8c0eccd223fae816 100644 (file)
--- a/src/domaintests.cpp
+++ b/src/domaintests.cpp
@@ -44,17 +44,179 @@ TEST( fitTests, fitTest_searchF0xyzabc)
 
 TEST( fitTests, fitTest_ApplyFit)
 {
-    // шаблон - применять готовое наложение
+    RVec mid1, mid2, x1, y1, z1, x2, y2, z2;
+
+    x1[0] = 1; x1[1] = 0; x1[2] = 0;
+    y1[0] = 0; y1[1] = 1; y1[2] = 0;
+    z1[0] = 0; z1[1] = 0; z1[2] = 1;
+
+    x2[0] = -1; x2[1] = 0; x2[2] = 0;
+    y2[0] = 0; y2[1] = -1; y2[2] = 0;
+    z2[0] = 0; z2[1] = 0; z2[2] = -1;
+
+    std::vector< RVec > testFrame1, testFrame2, testFrame3;
+    testFrame1.resize(0);
+    testFrame2.resize(0);
+    std::vector< std::pair< unsigned int, unsigned int > > testPairs;
+    testPairs.resize(0);
+
+    for (unsigned int i = 0; i < 27; i++) {
+        testFrame1.push_back(x1 * static_cast< float >(i % 3) + y1 * static_cast< float >((i % 9) / 3) + z1 * static_cast< float >(i / 9));
+        testFrame2.push_back(x2 * static_cast< float >(i % 3) + y2 * static_cast< float >((i % 9) / 3) + z2 * static_cast< float >(i / 9));
+        testPairs.push_back(std::make_pair(i, i));
+    }
+
+    testFrame3.resize(0);
+    testFrame3 = testFrame2;
+    ApplyFit(testFrame3, 2, 2, 2, 0, 0, 0);
+    CalcMid(testFrame1, testFrame3, mid1, mid2, testPairs);
+    ASSERT_NEAR((mid1 - mid2).norm(), 0, 0.000001);
+
+    testFrame3.resize(0);
+    testFrame3 = testFrame1;
+    ApplyFit(testFrame3, 0, 0, 0, M_PI, 0, -M_PI / 2);
+    CalcMid(testFrame2, testFrame3, mid1, mid2, testPairs);
+    ASSERT_NEAR((mid1 - mid2).norm(), 0, 0.000001);
+
+    testFrame3.resize(0);
+    testFrame3 = testFrame1;
+    ApplyFit(testFrame3, 0, 0, 0, 0, M_PI, M_PI / 2);
+    CalcMid(testFrame2, testFrame3, mid1, mid2, testPairs);
+    ASSERT_NEAR((mid1 - mid2).norm(), 0, 0.000001);
+
+    testFrame3.resize(0);
+    testFrame3 = testFrame1;
+    ApplyFit(testFrame3, 0, 0, 0, M_PI, M_PI / 2, 0);
+    CalcMid(testFrame2, testFrame3, mid1, mid2, testPairs);
+    ASSERT_NEAR((mid1 - mid2).norm(), 0, 0.000001);
+
+    testFrame3.resize(0);
+    testFrame3 = testFrame1;
+    ApplyFit(testFrame3, 0, 0, 0, 3 * M_PI / 2, M_PI / 2, M_PI / 2);
+    CalcMid(testFrame2, testFrame3, mid1, mid2, testPairs);
+    ASSERT_NEAR((mid1 - mid2).norm(), 0, 0.000001);
+
+    testFrame3.resize(0);
+    testFrame3 = testFrame1;
+    ApplyFit(testFrame1, -1, -1, -1, M_PI, 0, -M_PI / 2);
+    ApplyFit(testFrame3, 0, 0, 0, M_PI, 0, -M_PI / 2);
+    ApplyFit(testFrame3, 1, 1, 1, 0, 0, 0);
+    ApplyFit(testFrame1, 2, 4, 8, M_PI / 6, M_PI / 6, M_PI / 6);
+    ApplyFit(testFrame3, 2, 4, 8, M_PI / 6, M_PI / 6, M_PI / 6);
+    CalcMid(testFrame1, testFrame3, mid1, mid2, testPairs);
+    ASSERT_NEAR((mid1 - mid2).norm(), 0, 0.000001);
+    double testF = 0;
+    for (unsigned int i = 0; i < testFrame1.size(); i++) {
+        testF += F(testFrame1[i][0], testFrame1[i][1], testFrame1[i][2], testFrame3[i][0], testFrame3[i][1], testFrame3[i][2], 0, 0, 0);
+    }
+    ASSERT_NEAR(testF, 0, 0.000001);
 }
 
 TEST( fitTests, fitTest_CalcMid)
 {
-    // шаблон - проверка функции
+    RVec mid1, mid2, x1, y1, z1, x2, y2, z2;
+
+    x1[0] = 1; x1[1] = 0; x1[2] = 0;
+    y1[0] = 0; y1[1] = 1; y1[2] = 0;
+    z1[0] = 0; z1[1] = 0; z1[2] = 1;
+
+    x2[0] = -1; x2[1] = 0; x2[2] = 0;
+    y2[0] = 0; y2[1] = -1; y2[2] = 0;
+    z2[0] = 0; z2[1] = 0; z2[2] = -1;
+
+    std::vector< RVec > testFrame1, testFrame2;
+    testFrame1.resize(0);
+    testFrame2.resize(0);
+    std::vector< std::pair< unsigned int, unsigned int > > testPairs;
+    testPairs.resize(0);
+
+    for (unsigned int i = 0; i < 27; i++) {
+        testFrame1.push_back(x1 * static_cast< float >(i % 3) + y1 * static_cast< float >((i % 9) / 3) + z1 * static_cast< float >(i / 9));
+        testFrame2.push_back(x2 * static_cast< float >(i % 3) + y2 * static_cast< float >((i % 9) / 3) + z2 * static_cast< float >(i / 9));
+        testPairs.push_back(std::make_pair(i, i));
+    }
+
+    CalcMid(testFrame1, testFrame1, mid1, mid2, testPairs);
+    ASSERT_NEAR((mid1 - mid2).norm(), 0, 0.000001);
+
+    CalcMid(testFrame2, testFrame2, mid1, mid2, testPairs);
+    ASSERT_NEAR((mid1 - mid2).norm(), 0, 0.000001);
+
+    CalcMid(testFrame1, testFrame2, mid1, mid2, testPairs);
+    ASSERT_NEAR((mid1 - mid2).norm(), 2 * std::sqrt(3), 0.000001);
 }
 
 TEST( fitTests, fitTest_MyFitNew)
 {
-    // шаблон - глобальное применение
+    RVec a, x1, y1, z1;
+
+    x1[0] = 1; x1[1] = 0; x1[2] = 0;
+    y1[0] = 0; y1[1] = 1; y1[2] = 0;
+    z1[0] = 0; z1[1] = 0; z1[2] = 1;
+
+    a[0] = -2; a[1] = -2; a[2] = -2;
+
+    std::vector< RVec > testFrame1, testFrame2, testFrame3;
+    testFrame1.resize(0);
+    testFrame2.resize(0);
+    testFrame3.resize(0);
+    std::vector< std::pair< unsigned int, unsigned int > > testPairs;
+    testPairs.resize(0);
+
+    for (unsigned int i = 0; i < 27; i++) {
+        testFrame1.push_back(x1 * static_cast< float >(i % 3) + y1 * static_cast< float >((i % 9) / 3) + z1 * static_cast< float >(i / 9));
+        testFrame2.push_back(testFrame1.back() + a);
+        testPairs.push_back(std::make_pair(i, i));
+    }
+
+    testFrame3 = testFrame2;
+    MyFitNew(testFrame1, testFrame3, testPairs, 0);
+    double testF = 0;
+    for (unsigned int i = 0; i < testFrame1.size(); i++) {
+        testF += F(testFrame1[i][0], testFrame1[i][1], testFrame1[i][2], testFrame3[i][0], testFrame3[i][1], testFrame3[i][2], 0, 0, 0);
+    }
+    ASSERT_NEAR(testF, 0, 0.000001);
+
+    testFrame3.resize(0);
+    testFrame3 = testFrame2;
+    ApplyFit(testFrame3, 7.1534, 0.5591, -3.1415, M_PI / 23, M_PI / 17, M_PI / 2.2);
+    MyFitNew(testFrame1, testFrame3, testPairs, 0.000001);
+    testF = 0;
+    for (unsigned int i = 0; i < testFrame1.size(); i++) {
+        testF += F(testFrame1[i][0], testFrame1[i][1], testFrame1[i][2], testFrame3[i][0], testFrame3[i][1], testFrame3[i][2], 0, 0, 0);
+    }
+    ASSERT_NEAR(testF, 0, 0.01);
+
+    testFrame3.resize(0);
+    testFrame3 = testFrame2;
+    ApplyFit(testFrame3, 7.1534, 0.5591, -3.1415, M_PI / 23, M_PI / 17, M_PI / 2.2);
+    MyFitNew(testFrame1, testFrame3, testPairs, 0.0000001);
+    testF = 0;
+    for (unsigned int i = 0; i < testFrame1.size(); i++) {
+        testF += F(testFrame1[i][0], testFrame1[i][1], testFrame1[i][2], testFrame3[i][0], testFrame3[i][1], testFrame3[i][2], 0, 0, 0);
+    }
+    ASSERT_NEAR(testF, 0, 0.001);
+
+    testFrame3.resize(0);
+    testFrame3 = testFrame2;
+    ApplyFit(testFrame3, 7.1534, 0.5591, -3.1415, M_PI / 23, M_PI / 17, M_PI / 2.2);
+    MyFitNew(testFrame1, testFrame3, testPairs, 0.0000001);
+    testF = 0;
+    for (unsigned int i = 0; i < testFrame1.size(); i++) {
+        testF += F(testFrame1[i][0], testFrame1[i][1], testFrame1[i][2], testFrame3[i][0], testFrame3[i][1], testFrame3[i][2], 0, 0, 0);
+    }
+    ASSERT_NEAR(testF, 0, 0.0001);
+
+    testFrame3.resize(0);
+    testFrame3 = testFrame2;
+    ApplyFit(testFrame3, 7.1534, 0.5591, -3.1415, M_PI / 23, M_PI / 17, M_PI / 2.2);
+    MyFitNew(testFrame1, testFrame3, testPairs, 0.00000001);
+    testF = 0;
+    for (unsigned int i = 0; i < testFrame1.size(); i++) {
+        testF += F(testFrame1[i][0], testFrame1[i][1], testFrame1[i][2], testFrame3[i][0], testFrame3[i][1], testFrame3[i][2], 0, 0, 0);
+    }
+    ASSERT_NEAR(testF, 0, 0.00001);
+
 }
 
 // т.к. домены это класс с кучей приватных функций ещё под вопросом как их тестить
@@ -77,6 +239,44 @@ TEST( domainTests, domainTest_update)
 }
 
 int main(int argc, char *argv[]) {
+
+    /*RVec a, x1, y1, z1;
+
+    x1[0] = 1; x1[1] = 0; x1[2] = 0;
+    y1[0] = 0; y1[1] = 1; y1[2] = 0;
+    z1[0] = 0; z1[1] = 0; z1[2] = 1;
+
+    a[0] = -2; a[1] = -2; a[2] = -2;
+
+    std::vector< RVec > testFrame1, testFrame2;
+    testFrame1.resize(0);
+    testFrame2.resize(0);
+    std::vector< std::pair< unsigned int, unsigned int > > testPairs;
+    testPairs.resize(0);
+
+    for (unsigned int i = 0; i < 27; i++) {
+        testFrame1.push_back(x1 * static_cast< float >(i % 3) + y1 * static_cast< float >((i % 9) / 3) + z1 * static_cast< float >(i / 9));
+        testFrame2.push_back(testFrame1.back() + a);
+        testPairs.push_back(std::make_pair(i, i));
+    }
+
+    ApplyFit(testFrame2, 7.1534, 0.5591, -3.1415, M_PI / 23, M_PI / 17, M_PI / 2.2);
+    MyFitNew(testFrame1, testFrame2, testPairs, 0.000000001);
+
+    std::cout << "\n";
+    for (unsigned int i = 0; i < 27; i++) {
+        std::cout << testFrame1[i][0] << " " << testFrame1[i][1] << " " << testFrame1[i][2] << " ";
+        std::cout << testFrame2[i][0] << " " << testFrame2[i][1] << " " << testFrame2[i][2] << "\n";
+    }
+    std::cout << "\n";
+
+    double testF = 0;
+    for (unsigned int i = 0; i < testFrame1.size(); i++) {
+        testF += F(testFrame1[i][0], testFrame1[i][1], testFrame1[i][2], testFrame2[i][0], testFrame2[i][1], testFrame2[i][2], 0, 0, 0);
+    }
+
+    std::cout << testF << "\n\n";*/
+
     ::testing::InitGoogleTest( &argc, argv );
     return RUN_ALL_TESTS();
 }

diff --git a/src/newfit.cpp b/src/newfit.cpp
index eee47daca2de1e6cec16a4e8dd22bb996c12e0ad..0343ac12088a1bb9b657ee238c25cf0dd3c0adb2 100644 (file)
--- a/src/newfit.cpp
+++ b/src/newfit.cpp
@@ -8,7 +8,7 @@ double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus
 }
 
 // вычисление функции F и её частных производных
-void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
+void searchF0xyzabc (long double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
     double t01, t02, t03, t04, t05, t06, t07;
     t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
     t02 = pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2);
@@ -68,8 +68,9 @@ void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &mi
 }
 
 // функция фитирования фрейма b на фрейм a на основе пар pairs с "точностью" FtCnst
-void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, double FtCnst) {
-    double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
+void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, long double FtCnst) {
+    double fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0;
+    long double f1 = 0, f2 = 0, l = 1;
     RVec ma, mb;
     CalcMid(a, b, ma, mb, pairs);
     ma -= mb;
@@ -80,9 +81,10 @@ void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::
                 static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y, static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
     }
     if (FtCnst == 0) {
-        FtCnst = 0.00001;
+        FtCnst = 0.0000001;
     }
     if (f1 == 0) {
+        ApplyFit(b, x, y, z, A, B, C);
         return;
     } else {
         // поиск оптимального геометрического преобразования методом градиентного спуска
@@ -105,7 +107,6 @@ void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::
                 }
                 if (f2 < f1) {
                     x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
-                    fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
                     break;
                 } else {
                     // по факту существует минимальное значение переменной типа double

diff --git a/src/newfit.h b/src/newfit.h
index b065bf54e42ecb0a403585f79f4db8a1cb103de4..474dbe812e7109b1de48459f9162cbeddf3e2b39 100644 (file)
--- a/src/newfit.h
+++ b/src/newfit.h
@@ -3,6 +3,8 @@
 
 #include <gromacs/trajectoryanalysis.h>
 #include "/home/titov_ai/Develop/gromacs-original/src/gromacs/trajectoryanalysis/topologyinformation.h"
+#include <gromacs/math/vectypes.h>
+#include <gromacs/math/vec.h>
 
 #include <iostream>
 #include <cfloat>
@@ -17,14 +19,15 @@ using namespace gmx;
 double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C);
 
 // вычисление функции F и её частных производных
-void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C);
+void searchF0xyzabc (long double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C);
 
 // применения геометрического преобразования: смещение на вектор (x, y, z) и повороты на градусы A, B, C относительно осей (x, y, z)
 void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C);
 
 // подсчёт геометрических центров множеств a и b на основе пар pairs
 void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< unsigned int, unsigned int > > pairs);
+
 // функция фитирования фрейма b на фрейм a на основе пар pairs с "точностью" FtCnst
-void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, double FtCnst);
+void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, long double FtCnst);
 
 #endif // NEWFIT_H