still need to find THE bug

probability is high its in the final part

diff --git a/src/domains.cpp b/src/domains.cpp
index ec4298e5b3d700e1ffe25c0279b4e58ff4663563..11019405e2a69a366d51d30b794681ba25f3de74 100644 (file)
--- a/src/domains.cpp
+++ b/src/domains.cpp
@@ -62,6 +62,8 @@ using namespace gmx;
 
 using gmx::RVec;
 
+const int tau_jump = 10;
+
 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) {
     return  sqrt(   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) +
                     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) +
@@ -316,24 +318,26 @@ bool check_domsizes(std::vector< std::vector< int > > cd_domsizes, int cd_domain
     return false;
 }
 
-void print_domains(std::vector< std::vector< unsigned int > > pd_domains, std::vector< int > index, std::string fnNdx_, int dms, double epsi, double delta) {
+void print_domains(std::vector< std::vector< std::vector< unsigned int > > > pd_domains, std::vector< int > index, std::string fnNdx_, int dms, double epsi, double delta, int window) {
     FILE               *fpNdx_;
     fpNdx_ = std::fopen(fnNdx_.c_str(), "w+");
     int write_count;
-    for (unsigned int i = 0; i < pd_domains.size(); i++) {
-        std::fprintf(fpNdx_, "[domain_№_%3d_%3d_%4.3f_%4.3f]\n", i + 1, dms, epsi, delta);
-        write_count = 0;
-        for (unsigned int j = 0; j < pd_domains[i].size(); j++) {
-            write_count++;
-            if (write_count > 20) {
-                write_count -= 20;
-                std::fprintf(fpNdx_, "\n");
+    for (unsigned int i1 = 0; i1 < pd_domains.size(); i1++) {
+        for (unsigned int i2 = 0; i2 < pd_domains[i1].size(); i2++) {
+            std::fprintf(fpNdx_, "[domain-stPos-%3d-num-%3d-dms-%3d-epsi-%4.3f-delta-%4.3f]\n", static_cast< int >(i1) * window / tau_jump, i2 + 1, dms, epsi, delta);
+            write_count = 0;
+            for (unsigned int j = 0; j < pd_domains[i1][i2].size(); j++) {
+                write_count++;
+                if (write_count > 20) {
+                    write_count -= 20;
+                    std::fprintf(fpNdx_, "\n");
+                }
+                std::fprintf(fpNdx_, "%5d ", index[pd_domains[i1][i2][j]] + 1);
+                std::printf("%5d ", index[pd_domains[i1][i2][j]] + 1);
             }
-            std::fprintf(fpNdx_, "%5d ", index[pd_domains[i][j]] + 1);
-            std::printf("%5d ", index[pd_domains[i][j]] + 1);
+            std::printf("\n\n");
+            std::fprintf(fpNdx_,"\n\n");
         }
-        std::printf("\n\n");
-        std::fprintf(fpNdx_,"\n\n");
     }
     std::fprintf(fpNdx_,"\n");
     std::fclose(fpNdx_);
@@ -378,7 +382,7 @@ class Domains : public TrajectoryAnalysisModule
         std::vector< std::vector< std::vector< std::vector< node > > > >    graph;
 
 
-        std::vector< std::vector< unsigned int > >                  domains;
+        std::vector< std::vector< std::vector< unsigned int > > >   domains;
         std::vector< std::vector< int > >                           domsizes;
 
         std::vector< int >                                          index;
@@ -387,9 +391,10 @@ class Domains : public TrajectoryAnalysisModule
         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;
+        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
@@ -461,7 +466,7 @@ Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
         }
     }
 
-    graph.resize(0);
+    graph.resize(1);
 
     bone = static_cast< unsigned int >(index.size() - static_cast< unsigned int >(domain_min_size) + 1);
     graph[0].resize(bone);
@@ -480,6 +485,14 @@ Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
 void
 Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
 {
+    if ((frnr + 1) % (window / tau_jump) == 0) {
+        graph.resize(graph.size() + 1);
+        graph.back().resize(bone);
+        for (unsigned int i = 0; i < bone; i++) {
+            make_graph(index.size(), reference, graph.back()[i]);
+        }
+    }
+
     for (unsigned int i = 0; i < index.size(); i++) {
         trajectory[i] = fr.x[index[i]];
     }
@@ -489,7 +502,9 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnal
     for (unsigned int j = 0; j < bone; j++) {
         std::vector < RVec > temp = trajectory;
         MyFitNew(reference, temp, w_rls2[j], 0);
-        update_graph(graph[0][j], temp, epsi);
+        for (unsigned int i = static_cast< unsigned int >(std::max(0, (frnr + 1) / (window / tau_jump) - tau_jump + 1)); i < graph.size(); i++) {
+            update_graph(graph[i][j], temp, epsi);
+        }
     }
     #pragma omp barrier
     std::cout << "frame: " << frames << " analyzed\n";
@@ -498,28 +513,44 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnal
 void
 Domains::finishAnalysis(int nframes)
 {
+    graph.resize(graph.size() - tau_jump);
+    std::cout << "graph size: " << graph.size() << "\n";
     frames -= 1;
 
-    std::cout << "final cheking\n";
-    check_domains(delta, frames, graph[0]);
-
-    std::cout << "finding domains' sizes\n";
-    find_domain_sizes(graph[0], domsizes);
+    domains.resize(graph.size());
 
-    std::cout << "finding domains\n";
     std::vector< unsigned int > a;
     a.resize(0);
-    while (check_domsizes(domsizes, domain_min_size)) {
-        domains.push_back(a);
-        if (DomainSearchingAlgorythm == 0) {
-            get_maxsized_domain(domains.back(), graph[0], domsizes);
-        } else if (DomainSearchingAlgorythm == 1) {
-            get_min_domain(domains.back(), graph[0], domsizes, domain_min_size);
-        }
-        std::cout << "new domain: " << domains.back().size() << " atoms\n";
-        delete_domain_from_graph(graph[0], domains.back());
+
+    std::cout << "final cheking\n";
+
+    for (unsigned int i = 0; i < graph.size(); i++) {
         domsizes.resize(0);
-        find_domain_sizes(graph[0], domsizes);
+        domains[i].resize(0);
+        check_domains(delta, frames, graph[i]);
+        std::cout << "finding domains' sizes\n";
+        find_domain_sizes(graph[i], domsizes);
+
+        /*for (int j1 = 0; j1 < domsizes.size(); j1++) {
+            for (int j2 = 0; j2 < domsizes[j1].size(); j2++) {
+                if (domsizes[j1][j2] != 0) {
+                    std::cout << domsizes[j1][j2] << " ";
+                }
+            }
+        }*/
+
+        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);
+        }
     }
 }
 
@@ -527,7 +558,7 @@ void
 Domains::writeOutput()
 {
     std::cout << "making output file\n";
-    print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta); // see function for details | numbers from index
+    print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta, window); // see function for details | numbers from index
     std::cout << "\n END \n";
 }