- implemented new algorythm

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index 0b0a8922fc43b01c6c241e787e186408db4fd625..220de09249e463ca609fb18276ac5904d3ffa19e 100644 (file)
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -101,6 +101,20 @@ void make_rout_file(double crl_border, std::vector< int > indx, std::vector< std
     std::fclose(file);
 }
 
+void make_rout_file2(double crl_border, std::vector< int > indx, std::vector< std::vector< std::pair< int, int > > > rout, const char* file_name)
+{
+    FILE *file;
+    file = std::fopen(file_name, "w+");
+    std::fprintf(file, "correlations >= %0.2f\n\n", crl_border);
+    for (int i = 0; i < rout.size(); i++) {
+        for (int j = 0; j < rout[i].size() - 1; j++) {
+            std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.15\n", indx[rout[i][j].first] + 1, indx[rout[i][j].second] + 1);
+        }
+        std::fprintf(file, "\n\n");
+    }
+    std::fclose(file);
+}
+
 void make_best_corrs_graphics(std::vector< std::vector< std::vector< double > > > correlations,
                               std::vector< std::pair< int, int > > corr_pairs,
                               std::vector< int > indx,
@@ -139,6 +153,10 @@ bool isitsubset (std::vector< int > a, std::vector< int > b) {
     }
 }
 
+bool myfunction (const std::pair< int, double > i, const std::pair< int, double > j) {
+    return i.second < j.second;
+}
+
 /*! \brief
  * \ingroup module_trajectoryanalysis
  */
@@ -420,7 +438,7 @@ Domains::finishAnalysis(int nframes)
     }
 
     std::cout << "Correlation's evaluation - start\n\n";
-    #pragma omp parallel
+    #pragma omp parallel shared(crltns)
     {
         #pragma omp for schedule(dynamic)
         for (int i = m; i <= k; i += 1) {
@@ -500,7 +518,6 @@ Domains::finishAnalysis(int nframes)
     //
 
     std::cout << "Correlation's overview - start\n\n";
-    int imax = 0;
     for (int i1 = 0; i1 < 100; i1++) {
         int i5 = 0;
         for (int i2 = 1; i2 < crltns.size(); i2++) {
@@ -512,38 +529,15 @@ Domains::finishAnalysis(int nframes)
                 }
             }
         }
-        if (i5 > 0) {
-            imax = i1;
-        }
         std::cout << i1 << " - " << i5 << " | ";
         if (i1 % 10 == 0) {
             std::cout << "\n" ;
         }
     }
-    std::cout << "\n";
+    std::cout << "\nCorrelation's overview - end \n\n";
 
-    std::set< std::pair< int, int > > crr_prs;
-    crr_prs.clear();
-    std::vector< std::pair< int, int > > corr_prs;
-    corr_prs.resize(0);
-    for (int i2 = 1; i2 < crltns.size(); i2++) {
-        for (int i3 = 0; i3 < crltns[i2].size(); i3++) {
-            for (int i4 = 0; i4 < crltns[i2][i3].size(); i4++) {
-                if (std::abs(crltns[i2][i3][i4]) >= (double)(imax - 5) / 100 && i3 != i4) {
-                    std::cout << " found it \n";
-                    crr_prs.insert(std::make_pair(i3, i4));
-                }
-            }
-        }
-    }
-    for (std::set< std::pair< int, int > >::iterator it = crr_prs.begin(); it != crr_prs.end(); it++) {
-        corr_prs.push_back(*it);
-    }
-    std::cout << "Correlation's overview - end \n\n";
-    std::cout << imax << " " << crr_prs.size() << " " << corr_prs.size() << "\n\n";
-    make_best_corrs_graphics(crltns, corr_prs, index, "best_graphics.txt");
 
-    std::vector< std::vector< int > > graph;
+    /*std::vector< std::vector< int > > graph;
     graph.resize(index.size());
     rvec temp;
 
@@ -568,10 +562,6 @@ Domains::finishAnalysis(int nframes)
         }
     }
 
-    //kross correlation
-    //corelations for pairs in the routs
-    //check on self to self connections
-
     std::cout << "\n";
     std::vector< std::vector < int > > rout_old, rout_new, rout_out, rout_stack;
     bool flag = true;
@@ -629,9 +619,9 @@ Domains::finishAnalysis(int nframes)
     std::vector< bool > rout_flags;
     rout_flags.resize(rout_out.size(), false);
     rout_new.resize(0);
-    //#pragma omp parallel
-    //{
-        //#pragma omp for schedule(dynamic)
+    #pragma omp parallel
+    {
+        #pragma omp for schedule(dynamic)
         for (int i = rout_old.size() - 1; i >= 0; i--) {
             bool lflag1 = true;
             for (int j = 0; j < i; j++) {
@@ -648,22 +638,168 @@ Domains::finishAnalysis(int nframes)
                 }
             }
             if (lflag1) {
-                //rout_flags[i] = true;
-                rout_new.push_back(rout_old[i]);
+                rout_flags[i] = true;
+                //rout_new.push_back(rout_old[i]);
             }
         }
-    //}
+    }
 
-    /*for (int i = 0; i < rout_flags.size(); i++) {
+    for (int i = 0; i < rout_flags.size(); i++) {
         if (rout_flags[i]) {
             rout_new.push_back(rout_old[i]);
             std::cout << " " << i << "\n";
         }
-    }*/
+    }
     std::cout << "Path's finding - end\n\n";
 
+    std::set< std::pair< int, int > > crr_prs;
+    crr_prs.clear();
+    std::vector< std::pair< int, int > > corr_prs;
+    corr_prs.resize(0);
+    for (int i = 0; i < rout_new.size(); i++) {
+        for (int j = 0; j < rout_new[i].size() - 1; j++) {
+            crr_prs.insert(std::make_pair(rout_new[i][j], rout_new[i][j + 1]));
+        }
+    }
+    for (std::set< std::pair< int, int > >::iterator it = crr_prs.begin(); it != crr_prs.end(); it++) {
+        corr_prs.push_back(*it);
+    }
+    make_best_corrs_graphics(crltns, corr_prs, index, "best_graphics.txt");
+    */
+
+    std::vector< std::vector< std::pair< double, int > > > graph;
+    graph.resize(index.size());
+    for (int i = 0; i < index.size(); i++) {
+        graph[i].resize(index.size(), std::make_pair(0, -1));
+    }
+    rvec temp;
+    for (int i = 1; i <= k; i++) {
+        for (int j = 0; j < index.size(); j++) {
+            for (int f = j; f < index.size(); f++) {
+                copy_rvec(frankenstein_trajectory[basic_frame][j], temp);
+                rvec_dec(temp, frankenstein_trajectory[basic_frame][f]);
+                if (std::max(std::abs(crltns[i][j][f]), std::abs(crltns[i][f][j])) >= crl_border && norm(temp) <= eff_rad && std::abs(graph[j][f].first) < std::max(std::abs(crltns[i][j][f]), std::abs(crltns[i][f][j]))) {
+                    if (std::abs(crltns[i][j][f]) > std::abs(crltns[i][f][j])) {
+                        graph[j][f].first = crltns[i][j][f];
+                    } else {
+                        graph[j][f].first = crltns[i][f][j];
+                    }
+                    graph[j][f].second = i;
+                }
+            }
+        }
+    }
+    std::vector< bool > graph_flags;
+    graph_flags.resize(index.size(), true);
+    std::vector< std::vector< int > > sub_graph;
+    std::vector< std::vector< std::pair< int, int > > > sub_graph_rbr;
+    std::vector< int > a;
+    std::vector< std::pair< int, int > > b;
+    a.resize(0);
+    b.resize(0);
+    for (int i = 0; i < index.size(); i++) {
+        if (graph_flags[i]) {
+            sub_graph.push_back(a);
+            sub_graph_rbr.push_back(b);
+            std::vector< int > que1, que2, que3;
+            que1.resize(0);
+            que2.resize(0);
+            que3.resize(0);
+            que1.push_back(i);
+            que3.push_back(i);
+            graph_flags[i] = false;
+            while(que1.size() > 0) {
+                que2.clear();
+                for (int k = 0; k < que1.size(); k++) {
+                    for (int j = 0; j < index.size(); j++) {
+                        if (graph[que1[k]][j].second > -1 && graph_flags[j]) {
+                            que2.push_back(j);
+                            graph_flags[j] = false;
+                        }
+                    }
+                }
+                que1 = que2;
+                for (int j = 0; j < que2.size(); j++) {
+                    que3.push_back(que2[j]);
+                }
+            }
+            sub_graph.back() = que3;
+            for (int j = 0; j < que3.size(); j++) {
+                for (int k = 0; k < index.size(); k++) {
+                    if (graph[que3[j]][k].second > -1) {
+                        sub_graph_rbr.back().push_back(std::make_pair(que3[j], k));
+                    }
+                }
+            }
+            std::cout << sub_graph.back().size() << "   ";
+        }
+    }
+    std::cout << "sub graphs - end\n\n";
+
+    std::vector< std::vector < std::pair< int, int > > > rout_new;
+    for (int i = 0; i < sub_graph.size(); i++) {
+        if (sub_graph[i].size() > 2) {
+            std::vector< double > key;
+            std::vector< int > path;
+            std::vector< std::pair< int, double > > que;
+            int v;
+            key.resize(index.size(), 2);
+            path.resize(index.size(), -1);
+            key[sub_graph[i][0]] = 0;
+
+            for (int j = 0; j < sub_graph[i].size(); j++) {
+                que.push_back(std::make_pair(sub_graph[i][j], key[sub_graph[i][j]]));
+            }
+            std::sort(que.begin(), que.end(), myfunction);
+
+            while (!que.empty()) {
+                v = que[0].first;
+                que.erase(que.begin());
+                for (int j = 0; j < sub_graph_rbr[i].size(); j++) {
+                    int u = -1;
+                    if (sub_graph_rbr[i][j].first == v) {
+                        u = sub_graph_rbr[i][j].second;
+                    } else if (sub_graph_rbr[i][j].second == v) {
+                        u = sub_graph_rbr[i][j].first;
+                    }
+                    bool flag = false;
+                    int pos = -1;
+                    for (int k = 0; k < que.size(); k++) {
+                        if (que[k].first == u) {
+                            flag = true;
+                            pos = k;
+                            k = que.size() + 1;
+                        }
+                    }
+                    if (flag && key[u] > 1 - std::abs(graph[v][u].first)) {
+                        path[u] = v;
+                        key[u] = 1 - std::abs(graph[v][u].first);
+                        que[pos].second = key[u];
+                        sort(que.begin(), que.end(), myfunction);
+                    }
+                }
+            }
+            std::vector < std::pair< int, int > > a;
+            a.resize(0);
+            rout_new.push_back(a);
+            for (int j = 0; j < index.size(); j++) {
+                if (path[j] != -1) {
+                    rout_new.back().push_back(std::make_pair(j, path[j]));
+                }
+            }
+        }
+    }
+    make_rout_file2(crl_border, index, rout_new, "routs_new.txt");
+
+
+
+
+
+
+
+
     //
-    make_rout_file(crl_border, index, rout_new, "routs.txt");
+    //make_rout_file(crl_border, index, rout_new, "routs.txt");
     //
 
     std::cout << "Finish Analysis - end\n\n";