- changed the way routs are found

- added potential multithreading to the last part

diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp
index d7645fa5c0bd78b4aa95830be9c4199762752c68..0b0a8922fc43b01c6c241e787e186408db4fd625 100644 (file)
--- a/src/spacetimecorr.cpp
+++ b/src/spacetimecorr.cpp
@@ -120,6 +120,25 @@ void make_best_corrs_graphics(std::vector< std::vector< std::vector< double > >
 
 bool mysortfunc (std::vector< int > a, std::vector< int > b) { return (a.size() > b.size()); }
 
+bool isitsubset (std::vector< int > a, std::vector< int > b) {
+    if (b.size() == 0) {
+        return true;
+    }
+    std::sort(a.begin(), a.end());
+    std::sort(b.begin(), b.end());
+    int k = 0;
+    for (int i = 0; i < a.size(); i++) {
+        if (b[k] == a[i]) {
+            k++;
+        }
+    }
+    if (k == b.size()) {
+        return true;
+    } else {
+        return false;
+    }
+}
+
 /*! \brief
  * \ingroup module_trajectoryanalysis
  */
@@ -487,7 +506,7 @@ Domains::finishAnalysis(int nframes)
         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)i1 / 100) {
+                    if (std::abs(crltns[i2][i3][i4]) >= (double)i1 / 100 && i3 != i4) {
                         i5++;
                     }
                 }
@@ -510,7 +529,7 @@ Domains::finishAnalysis(int nframes)
     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) {
+                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));
                 }
@@ -549,11 +568,16 @@ 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;
+    std::vector< std::vector < int > > rout_old, rout_new, rout_out, rout_stack;
     bool flag = true;
     rout_new.resize(0);
     rout_old.resize(index.size());
+    rout_stack.resize(0);
     for (int i = 0; i < index.size(); i++) {
         rout_old[i].resize(1, i);
     }
@@ -562,12 +586,11 @@ Domains::finishAnalysis(int nframes)
         flag = false;
         std::vector < int > flag3;
         for (long int i = 0; i < rout_old.size(); i++) {
-            if (graph[rout_old[i].back()].size() == 0 && rout_old[i].size() > 2) {
-                rout_new.push_back(rout_old[i]);
+            if (isitsubset(rout_old[i], graph[rout_old[i].back()]) && rout_old[i].size() > 2) {
+                rout_stack.push_back(rout_old[i]);
             } else {
-                bool flag222 = true;
-                for (long int j = 0; j < graph[rout_old[i].back()].size(); j++) {
-                    bool flag2 = true, flag22 = true;
+                for (int j = 0; j < graph[rout_old[i].back()].size(); j++) {
+                    bool flag22 = true;
                     flag3 = rout_old[i];
                     for (int f = 0; f < flag3.size(); f++) {
                         if (flag3[f] == graph[rout_old[i].back()][j]) {
@@ -576,22 +599,10 @@ Domains::finishAnalysis(int nframes)
                     }
                     if (flag22) {
                         flag3.push_back(graph[rout_old[i].back()][j]);
-                        for (int f = 0; f < rout_new.size(); f++) {
-                            if (rout_new[f] == flag3) {
-                                flag2 = false;
-                            }
-                        }
-                        if (flag2) {
-                            rout_new.push_back(flag3);
-                            flag = true;
-                            flag222 = false;
-
-                        }
+                        rout_new.push_back(flag3);
+                        flag = true;
                     }
                 }
-                if (flag222 && rout_old[i].size() > 2) {
-                    rout_new.push_back(rout_old[i]);
-                }
             }
         }
         if (rout_new.size() != 0) {
@@ -601,33 +612,54 @@ Domains::finishAnalysis(int nframes)
         }
     }
 
+    for (int i = 0; i < rout_stack.size(); i++) {
+        rout_old.push_back(rout_stack[i]);
+    }
+    rout_stack.resize(0);
+
     std::sort(rout_old.begin(), rout_old.end(), mysortfunc);
 
+    std::cout << "\n\nClearing subrouts " << rout_old.size() << "\n\n";
+
     rout_out = rout_old;
     for (int i = 0; i < rout_out.size(); i++) {
         std::sort(rout_out[i].begin(), rout_out[i].end());
     }
 
+    std::vector< bool > rout_flags;
+    rout_flags.resize(rout_out.size(), false);
     rout_new.resize(0);
-    for (int i = rout_old.size() - 1; i >= 0; i--) {
-        bool lflag1 = true;
-        for (int j = 0; j < i; j++) {
-            if (rout_old[i].size() <= rout_old[j].size()) {
-                int la = 0;
-                for (int f = 0; f < rout_old[j].size(); f++) {
-                    if (rout_out[j][f] == rout_out[i][la]) {
-                        la++;
+    //#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++) {
+                if (rout_old[i].size() <= rout_old[j].size()) {
+                    int la = 0;
+                    for (int f = 0; f < rout_old[j].size(); f++) {
+                        if (rout_out[j][f] == rout_out[i][la]) {
+                            la++;
+                        }
+                    }
+                    if (la == rout_out[i].size()) {
+                        lflag1 = false;
                     }
                 }
-                if (la == rout_out[i].size()) {
-                    lflag1 = false;
-                }
+            }
+            if (lflag1) {
+                //rout_flags[i] = true;
+                rout_new.push_back(rout_old[i]);
             }
         }
-        if (lflag1) {
+    //}
+
+    /*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";
 
     //