changed the way frames are analyzed

diff --git a/src/domains.cpp b/src/domains.cpp
index 0ffedb2459d7264505fe7081b6c6189fd553f207..ac6e3808d2876929d896bb5b3bbae72902977c87 100644 (file)
--- a/src/domains.cpp
+++ b/src/domains.cpp
@@ -61,7 +61,7 @@ using namespace gmx;
 using gmx::RVec;
 
 struct node {
-    int  n;
+    short int n;
     RVec r;
     bool yep;
 };
@@ -242,18 +242,24 @@ class Domains : public TrajectoryAnalysisModule
     private:
 
         std::string                                                 fnNdx_;
+
         std::vector< std::vector< std::vector< node > > >           graph;
+
         std::vector< std::vector< int > >                           domains;
         std::vector< std::vector< int > >                           domsizes;
+
         std::vector< int >                                          index;
         std::vector< int >                                          numbers;
         std::vector< std::vector < RVec > >                         trajectory;
         Selection                                                   selec;
         int                                                         frames              = 0;
-        int                                                         etalon_frame        = 0; // should be selectable
         int                                                         domain_min_size     = 5; // should be selectable
+
+        real                                                        **w_rls;
+        int                                                         bone;
         double                                                      delta               = 0.90; //0.95 // should be selectable
         double                                                      epsi                = 0.15; //0.3 колебания внутри домена // should be selectable
+
         int                                                         domains_ePBC;
         // Copy and assign disallowed by base.
 };
@@ -284,10 +290,6 @@ Domains::initOptions(IOptionsContainer          *options,
                             .store(&fnNdx_).defaultBasename("domains")
                             .description("Index file from the domains"));
     // Add option for etalon_frame constant
-    options->addOption(gmx::IntegerOption("etalon_frame")
-                            .store(&etalon_frame)
-                            .description("basic frame to base evaluation on"));
-    // Add option for etalon_frame constant
     options->addOption(gmx::IntegerOption("dms")
                             .store(&domain_min_size)
                             .description("minimum domain size"));
@@ -327,40 +329,15 @@ Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
     for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
         index.push_back(*ai);
     }
-    trajectory.resize(frames + 1);
-    trajectory[frames].resize(selec.atomCount());
+    trajectory.resize(2);
+    trajectory[0].resize(selec.atomCount());
 
     for (int i = 0; i < selec.atomCount(); i++) {
-        trajectory[frames][i] = fr.x[index[i]];
+        trajectory[0][i] = fr.x[index[i]];
     }
-    frames++;
-    graph.resize(0);
-}
-
-void
-Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
-                      TrajectoryAnalysisModuleData *pdata)
-{
-    trajectory.resize(frames + 1);
-    trajectory[frames].resize(index.size());
-    for (int i = 0; i < index.size(); i++) {
-        trajectory[frames][i] = fr.x[index[i]];
-    }
-    frames++;
-}
 
-//domains -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA'
-//domains -s '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.tpr' -f '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.10th.xtc' -select 'name CA'
-//domains -s '/home/toluk/Data/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Data/reca_rd/reca_rd.mono.xtc' - select 'name CA'
-
-void
-Domains::finishAnalysis(int nframes)
-{
-    frames--;
-    int bone = index.size() - domain_min_size + 1;
+    bone = index.size() - domain_min_size + 1;
     graph.resize(bone);
-    real **w_rls;
-
     snew(w_rls, bone);
     for (int i = 0; i < bone; i++) {
         snew(w_rls[i], index.size());
@@ -371,115 +348,69 @@ Domains::finishAnalysis(int nframes)
                 w_rls[i][j] = 0;
             }
         }
-    }
-
-    /*
-    rvec **etalon;
-    snew(etalon, bone);
-    for (int i = 0; i < bone; i++) {
-        snew(etalon[i], index.size());
-        for (int j = 0; j < index.size(); j++) {
-            copy_rvec(trajectory[etalon_frame][j].as_vec(), etalon[i][j]);
-        }
-    }
-    */
-
-    for (int i = 0; i < bone; i++) {
         rvec *etalon;
         snew(etalon, index.size());
         for (int j = 0; j < index.size(); j++) {
-            copy_rvec(trajectory[etalon_frame][j].as_vec(), etalon[j]);
+            copy_rvec(trajectory[0][j].as_vec(), etalon[j]);
         }
         reset_x(index.size(), NULL, index.size(), NULL, etalon, w_rls[i]);
         make_graph(index.size(), etalon, graph[i]);
         sfree(etalon);
     }
+    trajectory[1].resize(index.size());
+}
 
-    std::chrono::time_point<std::chrono::system_clock> start1, end1;
-    start1 = std::chrono::system_clock::now();
-
-    /*
-    rvec **traj;
-    snew(traj, bone);
+void
+Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
+                      TrajectoryAnalysisModuleData *pdata)
+{
+    for (int i = 0; i < index.size(); i++) {
+        trajectory[1][i] = fr.x[index[i]];
+    }
+    frames++;
 
-    for (int i = 1; i < frames; i++) {
+    #pragma omp parallel
+    {
+        #pragma omp for schedule(dynamic)
         for (int j = 0; j < bone; j++) {
-            sfree(traj[j]);
-            snew(traj[j], index.size());
+            rvec *etalon, *traj;
+            snew(etalon, index.size());
             for (int k = 0; k < index.size(); k++) {
-                copy_rvec(trajectory[i][k].as_vec(), traj[j][k]);
-            }
-        }
-        #pragma omp parallel
-        {
-            #pragma omp for schedule(dynamic)
-            for (int j = 0; j < bone; j++) {
-                //#pragma omp ordered
-                //{
-                reset_x(index.size(), NULL, index.size(), NULL, etalon[j], w_rls[j]);
-                reset_x(index.size(), NULL, index.size(), NULL, traj[j], w_rls[j]);
-                do_fit(index.size(), w_rls[j], etalon[j], traj[j]);
-                update_graph(graph[j], traj[j], epsi);
-                //}
+                copy_rvec(trajectory[0][k].as_vec(), etalon[k]);
             }
-        }
-        std::cout << i << " / " << frames << " processed\n";
-        if (i % 25 == 0) {
-            end1 = std::chrono::system_clock::now();
-            int seconds = 0, minutes = 0, hours = 0;
-            seconds = std::chrono::duration_cast<std::chrono::seconds>(end1-start1).count();
-            seconds = (float)seconds * (float)( (float)(frames - i) / (float)i );
-            hours = seconds / 3600;
-            seconds -= hours * 3600;
-            minutes = seconds / 60;
-            seconds %= 60;
-            std::cout << "foretold time: " << hours << " hour(s) " << minutes << " minute(s) " << seconds << "second(s)\n";
-        }
-    }
-    */
-
-    for (int i = 1; i < frames; i++) {
-        #pragma omp parallel
-        {
-            #pragma omp for schedule(dynamic)
-            for (int j = 0; j < bone; j++) {
-                rvec *etalon, *traj;
-                snew(etalon, index.size());
-                for (int k = 0; k < index.size(); k++) {
-                    copy_rvec(trajectory[etalon_frame][k].as_vec(), etalon[k]);
-                }
-                snew(traj, index.size());
-                for (int k = 0; k < index.size(); k++) {
-                    copy_rvec(trajectory[i][k].as_vec(), traj[k]);
-                }
-                reset_x(index.size(), NULL, index.size(), NULL, etalon, w_rls[j]);
-                reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[j]);
-                do_fit(index.size(), w_rls[j], etalon, traj);
-                update_graph(graph[j], traj, epsi);
-                sfree(etalon);
-                sfree(traj);
+            snew(traj, index.size());
+            for (int k = 0; k < index.size(); k++) {
+                copy_rvec(trajectory[1][k].as_vec(), traj[k]);
             }
-        }
-        std::cout << i << " / " << frames << " processed\n";
-        if (i % 25 == 0) {
-            end1 = std::chrono::system_clock::now();
-            int seconds = 0, minutes = 0, hours = 0;
-            seconds = std::chrono::duration_cast<std::chrono::seconds>(end1-start1).count();
-            seconds = (float)seconds * (float)( (float)(frames - i) / (float)i );
-            hours = seconds / 3600;
-            seconds -= hours * 3600;
-            minutes = seconds / 60;
-            seconds %= 60;
-            std::cout << "foretold time: " << hours << " hour(s) " << minutes << " minute(s) " << seconds << "second(s)\n";
+            reset_x(index.size(), NULL, index.size(), NULL, etalon, w_rls[j]);
+            reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[j]);
+            do_fit(index.size(), w_rls[j], etalon, traj);
+            update_graph(graph[j], traj, epsi);
+            sfree(etalon);
+            sfree(traj);
         }
     }
-    std::cout << "finale cheking\n";
+    std::cout << "frame: " << frames << " analyzed\n";
+}
+
+//domains -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA'
+//domains -s '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.tpr' -f '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.10th.xtc' -select 'name CA'
+//domains -s '/home/toluk/Data/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Data/reca_rd/reca_rd.mono.xtc' - select 'name CA'
+
+void
+Domains::finishAnalysis(int nframes)
+{
+    frames -= 1;
+
+    std::cout << "final cheking\n";
     check_domains(delta, frames, graph);
+
     std::cout << "finding domains' sizes\n";
     find_domain_sizes(graph, domsizes);
+
+    std::cout << "finding domains\n";
     std::vector< int > a;
     a.resize(0);
-    std::cout << "finding domains\n";
     while (check_domsizes(domsizes, domain_min_size)) {
         domains.push_back(a);
         get_maxsized_domain(domains.back(), graph, domsizes);
@@ -489,12 +420,8 @@ Domains::finishAnalysis(int nframes)
     }
     for (int i = 0; i < bone; i++) {
         sfree(w_rls[i]);
-        //sfree(etalon[i]);
-        //sfree(traj[i]);
     }
     sfree(w_rls);
-    //sfree(etalon);
-    //sfree(traj);
 }
 
 void