- changed correlation evaluation

- added reference frame as reference

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index bc6bba86fdabb75e1c56dc8ff70204ade8453beb..348f7d72dfc419800c5d908cf45a74d88d42d1c7 100644 (file)
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -51,6 +51,7 @@
 #include <gromacs/trajectoryanalysis.h>
 #include <gromacs/utility/smalloc.h>
 #include <gromacs/math/do_fit.h>
+#include <gromacs/trajectoryanalysis/topologyinformation.h>
 
 using namespace gmx;
 
@@ -163,7 +164,7 @@ bool isitsubset (std::vector< int > a, std::vector< int > b) {
     }
 }
 
-void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame, std::vector< std::vector< std::vector< float > > > &crl, int tauS, int tauE) {
+void correlation_evaluation(std::vector< RVec > ref, std::vector< std::vector< RVec > > traj, int b_frame, std::vector< std::vector< std::vector< float > > > &crl, int tauS, int tauE) {
     crl.resize(tauE + 1);
     for (int i = 0; i < crl.size(); i++) {
         crl[i].resize(traj.front().size());
@@ -178,69 +179,33 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame
 
     std::vector< float > d;
     d.resize(traj.front().size(), 0);
-    //int nth = omp_get_thread_num();
-
-
-    //changed parrallel
-    //#pragma omp parallel shared(crl, temp_zero, d)
-    //{
-        #pragma omp parallel for schedule(dynamic)
-        for (int i = tauS; i <= tauE; i += 1) {
-            std::vector< RVec > medx, medy;
-            RVec temp1, temp2;
-            float tmp;
-            medx.resize(traj.front().size(), temp_zero);
-            medy.resize(traj.front().size(), temp_zero);
-            //#pragma omp for schedule(dynamic)
-            for (int j = 0; j < traj.size() - i - 1; j++) {
-                for (int f = 0; f < traj.front().size(); f++) {
-                    medx[f] += (traj[b_frame][f] - traj[j][f]);
-                    medy[f] += (traj[b_frame][f] - traj[j + i][f]);
-                }
-            }
-            //#pragma omp barrier
-            //#pragma omp for schedule(dynamic)
-            for (int j = 0; j < traj.front().size(); j++) {
-                tmp = traj.size() - 1;
-                tmp -= i;
-                tmp = 1 / tmp;
-                //tmp = 1 / (traj.size() - 1 - i);
-                medx[j] *= tmp;
-                medy[j] *= tmp;
-            }
-            //#pragma omp barrier
-            std::vector< std::vector< float > > a, b, c;
-            a.resize(traj.front().size(), d);
-            b.resize(traj.front().size(), d);
-            c.resize(traj.front().size(), d);
-            for (int j = 0; j < traj.size() - i - 1; j++) {
-                for (int f1 = 0; f1 < traj.front().size(); f1++) {
-                    //#pragma omp for schedule(dynamic)
-                    for (int f2 = 0; f2 < traj.front().size(); f2++) {
-                        temp1 = traj[b_frame][f1] - traj[j][f1] - medx[f1];
-                        temp2 = traj[b_frame][f2] - traj[j + i][f2] - medy[f2];
-                        a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
-                        b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
-                        c[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
-                    }
-                    //#pragma omp barrier
-                }
-            }
 
-            //#pragma omp for schedule(dynamic)
-            for (int j = 0; j < traj.front().size(); j++) {
-                for (int f = 0; f < traj.front().size(); f++) {
-                    crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
+    #pragma omp parallel for schedule(dynamic)
+    for (int i = tauS; i <= tauE; i += 1) {
+        RVec temp1, temp2;
+        std::vector< std::vector< float > > a, b, c;
+        a.resize(traj.front().size(), d);
+        b.resize(traj.front().size(), d);
+        c.resize(traj.front().size(), d);
+        for (int j = 0; j < traj.size() - i - 1; j++) {
+            for (int f1 = 0; f1 < traj.front().size(); f1++) {
+                for (int f2 = 0; f2 < traj.front().size(); f2++) {
+                    temp1 = traj[j][f1] - ref[f1];
+                    temp2 = traj[j + i][f2] - ref[f2];
+                    a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
+                    b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
+                    c[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
                 }
             }
-            //#pragma omp barrier
-            medx.resize(0);
-            medy.resize(0);
-            std::cout << i << " corr done\n";
         }
-        #pragma omp barrier
-    //}
-
+        for (int j = 0; j < traj.front().size(); j++) {
+            for (int f = 0; f < traj.front().size(); f++) {
+                crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
+            }
+        }
+        std::cout << i << " corr done\n";
+    }
+    #pragma omp barrier
 }
 
 void graph_calculation(std::vector< std::vector< std::pair< float, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr,
@@ -441,6 +406,7 @@ class SpaceTimeCorr : public TrajectoryAnalysisModule
         SelectionList                                               sel_;
 
         std::vector< std::vector< RVec > >                          trajectory;
+        std::vector< RVec >                                         reference;
 
         std::vector< int >                                          index;
         int                                                         frames              = 0;
@@ -494,6 +460,7 @@ SpaceTimeCorr::initOptions(IOptionsContainer          *options,
                             .description("<your name here> + <local file tag>.txt"));
     // Control input settings
     settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
+    settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
     settings->setPBC(true);
 }
 
@@ -507,6 +474,13 @@ SpaceTimeCorr::initAnalysis(const TrajectoryAnalysisSettings &settings,
         index.push_back(*ai);
     }
     trajectory.resize(0);
+
+    reference.resize(0);
+    if (top.hasFullTopology()) {
+        for (int i = 0; i < index.size(); i++) {
+            reference.push_back(top.x().at(index[i]));
+        }
+    }
 }
 
 void
@@ -523,6 +497,10 @@ SpaceTimeCorr::initAfterFirstFrame(const TrajectoryAnalysisSettings       &setti
 // -s '/home/toluk/Develop/samples/CubeTermal/CUBETermalTest1.pdb' -f '/home/toluk/Develop/samples/CubeTermal/CUBETermalTest.pdb' -n '/home/toluk/Develop/samples/CubeTermal/testCa.ndx' -sf '/home/toluk/Develop/samples/CubeTermal/SelListCa' -tau 900 -crl 0.20 -ef_rad 1
 // cube.000.000.10k.10.3.1stfrm
 // -s '/home/toluk/Develop/samples/JustCube/cube.000.000.10k.10.3.1stfrm.pdb' -f '/home/toluk/Develop/samples/JustCube/cube.000.000.10k.10.3.pdb' -n '/home/toluk/Develop/samples/JustCube/system.ndx' -sf '/home/toluk/Develop/samples/JustCube/SLsystem' -tau 1000 -crl 0.20 -ef_rad 9000 -out_put OLA
+
+// -s '*.pdb' -f '*.pdb' -n '*.ndx' -sf 'name' -tau 1000 -crl 0.20 -ef_rad 9000 -out_put OLA
+// -s '*.tpr' -f '*.xtc' -n '*.ndx' -sf 'name' -tau 1000 -crl 0.30 -ef_rad 20 -out_put 'test_run'
+
 void
 SpaceTimeCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
                       TrajectoryAnalysisModuleData *pdata)
@@ -549,21 +527,17 @@ SpaceTimeCorr::finishAnalysis(int nframes)
 
     std::cout << "\nCorrelation's evaluation - start\n";
 
-    correlation_evaluation(trajectory, basic_frame, crltns, m, k);
+    correlation_evaluation(reference, trajectory, basic_frame, crltns, m, k);
 
     make_correlation_matrix_file(crltns, (OutPutName + "_matrix.txt").c_str(), 0);
     std::cout << "corelation matrix printed\n";
+
     make_correlation_pairs_file(crltns, (OutPutName + "_pairs.txt").c_str(), 0);
     std::cout << "corelation pairs printed\n";
 
     std::cout << "Correlation's evaluation - end\n" << "graph evaluation: start\n";
 
-
-
-
-
     graph_calculation(graph, sub_graph, sub_graph_rbr, trajectory, basic_frame, crltns, crl_border, eff_rad, k);
-
     std::cout << "graph evaluation: end\n" << "routs evaluation: start\n";
 
     graph_back_bone_evaluation(rout_new, index.size(), graph, sub_graph, sub_graph_rbr);
@@ -572,14 +546,15 @@ SpaceTimeCorr::finishAnalysis(int nframes)
 
     make_rout_file(crl_border, index, rout_new, (OutPutName + "_routs.txt").c_str());
     std::cout << "corelation routs printed\n";
+
     make_best_corrs_graphics(crltns, rout_new, index, (OutPutName + "_routs_graphics.txt").c_str());
     std::cout << "corelation routs' pairs' graphics printed\n";
 
 
-    std::cout << "extra params\n";
+    /*std::cout << "extra params\n";
     std::vector< float > diffusion;
     evaluate_diffusion(trajectory, diffusion);
-    make_diffusion_file((OutPutName + "_diffusion.txt").c_str(), diffusion);
+    make_diffusion_file((OutPutName + "_diffusion.txt").c_str(), diffusion);*/
 
     std::cout << "Finish Analysis - end\n\n";
 }