optimized fitting algo

diff --git a/include/newfit.h b/include/newfit.h
index 0cd4ee08f2edbd334b8e7d17fbb9fcf93ed33d18..ec1e76d68b075fcf2eeb45673fe996dbf60ada1d 100644 (file)
--- a/include/newfit.h
+++ b/include/newfit.h
@@ -8,30 +8,36 @@ using gmx::RVec;
 
 namespace NewFit {
 
-    void implement_fit(std::vector< RVec > &if_frame, RVec if_step, double if_A, double if_B, double if_C) {
-        for (int i = 0; i < if_frame.size(); i++) {
-            if_frame[i] += if_step;
+    void center_frame(std::vector< RVec > &cf_frame, std::vector< int > cf_atoms) {
+        if (cf_atoms.empty()) {
+            cf_atoms.resize(0);
+            for (int i = 0; i < cf_frame.size(); i++) {
+                cf_atoms.push_back(i);
+            }
         }
-        for (int i = 0; i < if_frame.size(); i++) {
-            if_step[0] = if_frame[i][0];
-            if_step[1] = if_frame[i][1];
-            if_step[2] = if_frame[i][2];
-            if_frame[i][1] = if_step[1] * cos(if_A) - if_step[2] * sin(if_A);
-            if_frame[i][2] = if_step[1] * sin(if_A) + if_step[2] * cos(if_A);
+        RVec cf_temp;
+        cf_temp[0] = 0;
+        cf_temp[1] = 0;
+        cf_temp[2] = 0;
+        for (int i = 0; i < cf_atoms.size(); i++) {
+            cf_temp += cf_frame[cf_atoms[i]];
         }
-        for (int i = 0; i < if_frame.size(); i++) {
-            if_step[0] = if_frame[i][0];
-            if_step[1] = if_frame[i][1];
-            if_step[2] = if_frame[i][2];
-            if_frame[i][0] = if_step[0] * cos(if_B) + if_step[2] * sin(if_B);
-            if_frame[i][2] = - if_step[0] * sin(if_B) + if_step[2] * cos(if_B);
+        cf_temp[0] /= cf_atoms.size();
+        cf_temp[1] /= cf_atoms.size();
+        cf_temp[2] /= cf_atoms.size();
+        for (int i = 0; i < cf_atoms.size(); i++) {
+            cf_frame[cf_atoms[i]] -= cf_temp;
         }
+    }
+
+    void implement_fit(std::vector< RVec > &if_frame, RVec if_step, double if_A, double if_B, double if_C) {
+        RVec temp;
         for (int i = 0; i < if_frame.size(); i++) {
-            if_step[0] = if_frame[i][0];
-            if_step[1] = if_frame[i][1];
-            if_step[2] = if_frame[i][2];
-            if_frame[i][0] = if_step[0] * cos(if_C) - if_step[1] * sin(if_C);
-            if_frame[i][1] = if_step[0] * sin(if_C) + if_step[1] * cos(if_C);
+            if_frame[i] += if_step;
+            temp = if_frame[i];
+            if_frame[i][0] = cos(if_C) * cos(if_B) * temp[0] + (-sin(if_C)* cos(if_A) + sin(if_B) * cos(if_C) * sin(if_A)) * temp[1] + (sin(if_C) * sin(if_A) + sin(if_B) * cos(if_C) * cos(if_A)) * temp[2];
+            if_frame[i][1] = sin(if_C) * cos(if_B) * temp[0] + (cos(if_C) * cos(if_A) + sin(if_C) * sin(if_B) * sin(if_A)) * temp[1] + (-cos(if_C) * sin(if_A) + sin(if_C) * sin(if_B) * cos(if_A)) * temp[2];
+            if_frame[i][2] = -sin(if_B) * temp[0] + cos(if_B) * sin(if_A) * temp[1] + cos(if_B) * cos(if_A) * temp[2];
         }
     }
 
@@ -84,17 +90,24 @@ namespace NewFit {
         return distance;
     }
 
+    double dist_old_format(rvec *alfa, rvec *beta, std::vector< std::pair< int, int > > cros) {
+        double distance = 0;
+        for (int i = 0; i < cros.size(); i++) {
+            distance += F(alfa[cros[i].first][0], alfa[cros[i].first][1], alfa[cros[i].first][2], beta[cros[i].second][0], beta[cros[i].second][1], beta[cros[i].second][2], 1, 0, 1, 0, 1, 0);
+        }
+        return distance;
+    }
+
     void new_fit(std::vector< RVec > alfa, std::vector< RVec > &beta, std::vector< std::pair< int, int > > cros, double prec) {
         double bx, by, bz;
         RVec step;
         double x1 = 0, y1 = 0, z1 = 0, A1 = 0, B1 = 0, C1 = 0;
         double x2 = 0, y2 = 0, z2 = 0, A2 = 0, B2 = 0, C2 = 0;
         double cA, sA, cB, sB, cC, sC, temp1, temp2, temp3;
-        double H1 = 0.00001, H2 = 0.00001, F1 = 0, F2 = 0;
+        long double H1 = 0.0001, H2 = 0.0001, F1 = 0, F2 = 0, H1_temp, H2_temp, H1_const = 1, H2_const = M_PI/180;
+        F1 = dist(alfa, beta, cros) / cros.size();
+        bool flag = true;
         do {
-            if (F1 > F2 || F1 == 0) {
-                F1 = F2;
-            }
             F2 = 0;
             for (int i = 0; i < cros.size(); i++) {
                 bx = beta[cros[i].second][0] + x1;
@@ -116,26 +129,21 @@ namespace NewFit {
                 B2 += FB(temp1, temp2, temp3, bx, by, bz, cA, sA, cB, sB, cC, sC);
                 C2 += FC(temp1, temp2,        bx, by, bz, cA, sA, cB, sB, cC, sC);
             }
-            if (std::max(std::fabs(H1*x2), std::max(std::fabs(H1*y2), std::fabs(H1*z2))) > 1) {
-                H1 = 1 / std::max(std::fabs(H1*x2), std::max(std::fabs(H1*y2), std::fabs(H1*z2)));
-            }
-            x2 = x1 - H1 * x2;
-            y2 = y1 - H1 * y2;
-            z2 = z1 - H1 * z2;
-            if (std::max(std::fabs(H2*A2), std::max(std::fabs(H2*B2), std::fabs(H2*C2))) > M_PI/60) {
-                H2 = M_PI/60 / std::max(std::fabs(H2*A2), std::max(std::fabs(H2*B2), std::fabs(H2*C2)));
-            }
-            A2 = A1 - H2 * A2;
-            if (std::fabs(A2) > M_PI) {
-                A2 = 0;
-            }
-            B2 = B1 - H2 * B2;
-            if (std::fabs(B2) > M_PI) {
-                B2 = 0;
+            H1_temp = std::max(std::fabs(H1*x2), std::max(std::fabs(H1*y2), std::fabs(H1*z2)));
+            if (H1_temp > H1_const) {
+                H1 /= 10;
+            } else {
+                x2 = x1 - H1 * x2;
+                y2 = y1 - H1 * y2;
+                z2 = z1 - H1 * z2;
             }
-            C2 = C1 - H2 * C2;
-            if (std::fabs(C2) > M_PI) {
-                C2 = 0;
+            H2_temp = std::max(std::fabs(H2*A2), std::max(std::fabs(H2*B2), std::fabs(H2*C2)));
+            if (H2_temp > H2_const) {
+                H2 /= 10;
+            } else {
+                A2 = A1 - H2 * A2;
+                B2 = B1 - H2 * B2;
+                C2 = C1 - H2 * C2;
             }
             for (int i = 0; i < cros.size(); i++) {
                 cA = cos(A2);
@@ -146,10 +154,8 @@ namespace NewFit {
                 sC = sin(C2);
                 F2 += F(alfa[cros[i].first][0], alfa[cros[i].first][1], alfa[cros[i].first][2], beta[cros[i].second][0] + x2, beta[cros[i].second][1] + y2, beta[cros[i].second][2] + z2, cA, sA, cB, sB, cC, sC);
             }
-            if (F2 > F1 && F1 != 0) {
-                H1 /= 10;
-                H2 /= 10;
-            } else {
+            F2 /= cros.size();
+            if (F2 <= F1) {
                 x1 = x2;
                 y1 = y2;
                 z1 = z2;
@@ -162,8 +168,12 @@ namespace NewFit {
                 A2 = 0;
                 B2 = 0;
                 C2 = 0;
+                if (std::fabs(F1 - F2) <= prec) {
+                    flag = false;
+                }
+                F1 = F2;
             }
-        } while (std::fabs(F1 - F2) > prec/*0.00001*/);
+        } while (flag);
         step[0] = x1;
         step[1] = y1;
         step[2] = z1;