- minor fixes

diff --git a/src/domains.cpp b/src/domains.cpp
index f8449aba1bf476baf91f11be3b9712123df3f7e6..44c9545a36fc24aaafa6de5ffc988c1a4edd8db2 100644 (file)
--- a/src/domains.cpp
+++ b/src/domains.cpp
@@ -55,6 +55,7 @@
 #include "gromacs/topology/topology.h"
 #include "gromacs/topology/index.h"
 
+//#include "/home/titov_ai/Develop/fittingLib/fittinglib.h"
 
 #define MAX_NTRICVEC 12
 
@@ -115,8 +116,10 @@ void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &mi
     midB[2] = 0;
 
     for (unsigned int i = 0; i < pairs.size(); i++) {
-        rvec_inc(midA, a[pairs[i].first]);
-        rvec_inc(midB, b[pairs[i].second]);
+        midA += a[pairs[i].first];
+        midB += b[pairs[i].second];
+        //rvec_inc(midA, a[pairs[i].first]);
+        //rvec_inc(midB, b[pairs[i].second]);
     }
     midA[0] /= pairs.size();
     midA[1] /= pairs.size();
@@ -131,7 +134,8 @@ void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::
     double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
     RVec ma, mb;
     CalcMid(a, b, ma, mb, pairs);
-    rvec_dec(ma, mb);
+    ma -= mb;
+    //rvec_dec(ma, mb);
     double x = static_cast< double >(ma[0]), y = static_cast< double >(ma[1]), z = static_cast< double >(ma[2]), A = 0, B = 0, C = 0;
     for (unsigned int i = 0; i < pairs.size(); i++) {
         f1 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
@@ -143,7 +147,6 @@ void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::
     if (f1 == 0) {
         return;
     } else {
-        int count = 0;
         while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
             f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0; l = 1;
             for (unsigned int i = 0; i < pairs.size(); i++) {
@@ -159,26 +162,17 @@ void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::
                             static_cast< double >(b[pairs[i].second][0]) + (x - l * fx), static_cast< double >(b[pairs[i].second][1]) + (y - l * fy),
                             static_cast< double >(b[pairs[i].second][2]) + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
                 }
-                count++;
                 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 {
-                    //l *= 0.85;
                     l *= 0.50;
-                    /*if (count >= 14) {
-                        std::cout << count << " " << l << "\n";
-                    }*/
-                    if (l == DBL_MIN) { /*DBL_TRUE_MIN*/
+                    if (l == DBL_MIN) { //DBL_TRUE_MIN
                         break;
                     }
                 }
             }
-            count++;
-            if (count % 100000 == 0) {
-                std::cout << "+100k\n";
-            }
         }
         ApplyFit(b, x, y, z, A, B, C);
     }
@@ -341,6 +335,8 @@ void print_domains(std::vector< std::vector< unsigned int > > pd_domains, std::v
     }
     std::fprintf(fpNdx_,"\n");
     std::fclose(fpNdx_);
+    //std::fprintf(fpNdx2_,"\n");
+    std::fclose(fpNdx2_);
 }
 
 /*! \brief
@@ -376,31 +372,28 @@ class Domains : public TrajectoryAnalysisModule
 
     private:
 
-        std::string                                                 fnNdx_;
-
-        //std::vector< std::vector< std::vector< node > > >           graph;
-        std::vector< std::vector< std::vector< std::vector< node > > > >    graph;
+        std::string                                                             fnNdx_;
 
+        std::vector< std::vector< std::vector< std::vector< node > > > >        graph;
 
-        //std::vector< std::vector< std::vector< unsigned int > > >   domains;
-        std::vector< std::vector< unsigned int > >                  domains;
-        std::vector< std::vector< int > >                           domsizes;
+        std::vector< std::vector< unsigned int > >                              domains;
+        std::vector< std::vector< int > >                                       domsizes;
 
-        std::vector< int >                                          index;
-        std::vector< int >                                          numbers;
-        std::vector< RVec >                                         trajectory;
-        std::vector< RVec >                                         reference;
-        Selection                                                   selec;
-        int                                                         frames              = 0;
-        int                                                         window              = 1000; // selectable
-        int                                                         domain_min_size     = 4; // selectable
+        std::vector< int >                                                      index;
+        std::vector< int >                                                      numbers;
+        std::vector< RVec >                                                     trajectory;
+        std::vector< RVec >                                                     reference;
+        Selection                                                               selec;
+        int                                                                     frames              = 0;
+        int                                                                     window              = 1000; // selectable
+        int                                                                     domain_min_size     = 4; // selectable
 
-        std::vector< std::vector< std::pair< unsigned int, unsigned int > > >           w_rls2;
-        unsigned int                                                bone;
-        double                                                      delta               = 0.95; // selectable
-        double                                                      epsi                = 0.10; // selectable
+        std::vector< std::vector< std::pair< unsigned int, unsigned int > > >   w_rls2;
+        unsigned int                                                            bone;
+        double                                                                  delta               = 0.95; // selectable
+        double                                                                  epsi                = 0.10; // selectable
 
-        int                                                         DomainSearchingAlgorythm = 0; // selectable
+        int                                                                     DomainSearchingAlgorythm = 0; // selectable
         // Copy and assign disallowed by base.
 };
 
@@ -470,7 +463,6 @@ Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
 
     graph.resize(1);
     graph.front().resize(bone);
-    //graph.resize(bone);
 
     w_rls2.resize(bone);
     for (unsigned int i = 0; i < bone; i++) {
@@ -479,20 +471,9 @@ Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
             w_rls2[i].push_back(std::make_pair(j + i, j + i));
         }
         make_graph(index.size(), reference, graph.front()[i]);
-        //make_graph(index.size(), reference, graph[i]);
     }
 
     trajectory.resize(index.size());
-
-    std::vector< int > a;
-
-    /*for (int i = 0; i < 10; i++) {
-        a.push_back(i);
-    }
-    a.erase(a.begin());
-    for (int i = 0; i < 10; i++) {
-        std::cout << a[i] << " ";
-    }*/
 }
 
 void
@@ -529,7 +510,6 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnal
             find_domain_sizes(graph[0], domsizes);
         }
         graph.erase(graph.begin());
-        graph.resize(graph.size() - 1);
         std::cout << (frnr + 1) / (window / tau_jump) - tau_jump << " window's domains have been evaluated\n";
         print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta, window, (frnr + 1) / (window / tau_jump) - tau_jump); // see function for details | numbers from index
     }
@@ -549,52 +529,18 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnal
     }
     #pragma omp barrier
 
-
     std::cout << "frame: " << frnr << " analyzed\n";
 }
 
 void
 Domains::finishAnalysis(int nframes)
 {
-    /*graph.resize(graph.size() - tau_jump);
-    std::cout << "graph size: " << graph.size() << "\n";
-    frames -= 1;
-
-    domains.resize(graph.size());
-
-    std::vector< unsigned int > a;
-    a.resize(0);
-
-    std::cout << "final cheking\n";
 
-    for (unsigned int i = 0; i < graph.size(); i++) {
-        domsizes.resize(0);
-        domains[i].resize(0);
-        check_domains(delta, window, graph[i]);
-        std::cout << "finding domains' sizes\n";
-        find_domain_sizes(graph[i], domsizes);
-
-        while (check_domsizes(domsizes, domain_min_size)) {
-            domains[i].push_back(a);
-            if (DomainSearchingAlgorythm == 0) {
-                get_maxsized_domain(domains[i].back(), graph[i], domsizes);
-            } else if (DomainSearchingAlgorythm == 1) {
-                get_min_domain(domains[i].back(), graph[i], domsizes, domain_min_size);
-            }
-            std::cout << "new domain: " << domains[i].back().size() << " atoms\n";
-            delete_domain_from_graph(graph[i], domains[i].back());
-            domsizes.resize(0);
-            find_domain_sizes(graph[i], domsizes);
-        }
-        std::cout << i * window / tau_jump << " window's domains have been evaluated\n";
-    }*/
 }
 
 void
 Domains::writeOutput()
 {
-    //std::cout << "making output file\n";
-    //print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta, window); // see function for details | numbers from index
     std::cout << "\n END \n";
 }