- new fit
authorAnatoly Titov <toluk@omrb.pnpi.spb.ru>
Mon, 19 Nov 2018 14:17:48 +0000 (17:17 +0300)
committerAnatoly Titov <toluk@omrb.pnpi.spb.ru>
Mon, 19 Nov 2018 14:17:48 +0000 (17:17 +0300)
- outside functions
- new graph constant

src/spacetimecorr.cpp

index cd2eaed44fe7beb4ed8bd971e4c2adae299fa502..f930fb4d8fbaea087dedbd0965e8c2e0df290f37 100644 (file)
@@ -143,42 +143,51 @@ bool isitsubset (std::vector< int > a, std::vector< int > b) {
     }
 }
 
-void correlation_evaluation(std::vector< std::vector< RVec > > traj, int basic_frame, std::vector< std::vector< std::vector< double > > > &crltns, int tauS, int tauE) {
-    #pragma omp parallel shared(crltns)
+void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame, std::vector< std::vector< std::vector< double > > > &crl, int tauS, int tauE) {
+    crl.resize(tauE + 1);
+    for (int i = 0; i < crl.size(); i++) {
+        crl[i].resize(traj.front().size());
+        for (int j = 0; j < crl[i].size(); j++) {
+            crl[i][j].resize(traj.front().size(), 0);
+        }
+    }
+    rvec temp_zero;
+    temp_zero[0] = 0;
+    temp_zero[1] = 0;
+    temp_zero[2] = 0;
+    #pragma omp parallel shared(crl)
     {
         #pragma omp for schedule(dynamic)
         for (int i = tauS; i <= tauE; i += 1) {
-            rvec *medx, *medy, temp_zero;
+            real tmp2;
+            rvec *medx, *medy, temp, temp1, temp2;
             snew(medx, traj.front().size());
             snew(medy, traj.front().size());
-            temp_zero[0] = 0;
-            temp_zero[1] = 0;
-            temp_zero[2] = 0;
             for (int i = 0; i < traj.front().size(); i++) {
                 copy_rvec(temp_zero, medx[i]);
                 copy_rvec(temp_zero, medy[i]);
             }
             for (int j = 0; j < traj.size() - i - 1; j++) {
                 for (int f = 0; f < traj.front().size(); f++) {
-                    rvec temp;
-
-                    rvec_sub(traj[basic_frame][f], traj[j][f], temp);
+                    rvec_sub(traj[b_frame][f], traj[j][f], temp);
                     rvec_inc(medx[f], temp);
 
-                    rvec_sub(traj[basic_frame][f], traj[j + i][f], temp);
+                    rvec_sub(traj[b_frame][f], traj[j + i][f], temp);
                     rvec_inc(medy[f], temp);
                 }
             }
+
             for (int j = 0; j < traj.front().size(); j++) {
-                rvec temp;
-                real temp2 = traj.size() - 1;
+                tmp2 = 1 / (traj.size() - 1 - i);
+
+                /*temp2 = traj.size() - 1;
                 temp2 -= i;
-                temp2 = 1 / temp2;
+                temp2 = 1 / temp2;*/
 
                 copy_rvec(medx[j], temp);
-                svmul(temp2, temp, medx[j]);
+                svmul(tmp2, temp, medx[j]);
                 copy_rvec(medy[j], temp);
-                svmul(temp2, temp, medy[j]);
+                svmul(tmp2, temp, medy[j]);
             }
             std::vector< std::vector< double > > a, b, c;
             a.resize(traj.front().size());
@@ -192,11 +201,11 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int basic_f
             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++) {
-                        rvec temp1, temp2;
-                        rvec_sub(traj[basic_frame][f1], traj[j][f1], temp1);
+
+                        rvec_sub(traj[b_frame][f1], traj[j][f1], temp1);
                         rvec_dec(temp1, medx[f1]);
 
-                        rvec_sub(traj[basic_frame][f2], traj[j + i][f2], temp2);
+                        rvec_sub(traj[b_frame][f2], traj[j + i][f2], temp2);
                         rvec_dec(temp2, medy[f2]);
 
                         a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
@@ -207,20 +216,213 @@ void correlation_evaluation(std::vector< std::vector< RVec > > traj, int basic_f
             }
             for (int j = 0; j < traj.front().size(); j++) {
                 for (int f = 0; f < traj.front().size(); f++) {
-                    crltns[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
+                    crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
                 }
             }
             sfree(medx);
             sfree(medy);
+            std::cout << i << "\n";
+        }
+    }
+    #pragma omp barrier
+}
+
+void domain_chopping(SelectionList seList, std::vector< int > indx, std::vector< std::vector< int > > &filled_domains, std::vector< RVec > frame) {
+    ConstArrayRef< int > atomInd  = seList[0].atomIndices();
+    std::vector< std::vector< int > > init_domains;
+    init_domains.resize(seList.size());
+    for (int i = 0; i < seList.size(); i++) {
+        if (seList.size() != 1 && i == 0) {
+            continue;
+        }
+        atomInd  = seList[i].atomIndices();
+        for (ConstArrayRef<int>::iterator ai = atomInd.begin(); (ai < atomInd.end()); ai++) {
+            init_domains[i].push_back(*ai);
+        }
+    }
+    for (int i = 0; i < init_domains.size(); i++) {
+        for (int j = 0; j < init_domains[i].size(); j++) {
+            int k = 0;
+            while (indx[k] != init_domains[i][j]) {
+                k++;
+            }
+            init_domains[i][j] = k;
+        }
+    }
+    for (int i = 0; i < init_domains.size(); i++) {
+        if (init_domains[i].size() >= 2) {
+            filled_domains.push_back(init_domains[i]);
+            for (int k = 0; k < init_domains[i].size(); k++) {
+                for (int j = i + 1; j < init_domains.size(); j++) {
+                    for (int x = 0; x < init_domains[j].size(); x++) {
+                        if (init_domains[j][x] == init_domains[i][k]) {
+                            init_domains[j].erase(init_domains[j].begin() + x);
+                        }
+                    }
+                }
+            }
+        }
+    }
+    init_domains.resize(0);
+    init_domains = filled_domains;
+    std::vector< bool > flags;
+    flags.resize(indx.size(), true);
+    for (int i = 0; i < init_domains.size(); i++) {
+        for (int j = 0; j < init_domains[i].size(); j++) {
+            flags[init_domains[i][j]] = false;
+        }
+    }
+    int a;
+    rvec temp;
+    for (int i = 0; i < indx.size(); i++) {
+        if (flags[i]) {
+            float dist = 90000001;
+            for (int j = 0; j < init_domains.size(); j++) {
+                for (int k = 0; k < init_domains[j].size(); k++) {
+                    rvec_sub(frame[i], frame[init_domains[j][k]], temp);
+                    if (norm(temp) <= dist) {
+                        dist = norm(temp);
+                        a = j;
+                    }
+                }
+            }
+            filled_domains[a].push_back(i);
+            flags[i] = false;
+        }
+    }
+}
+
+void graph_calculation(std::vector< std::vector< std::pair< double, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr,
+                      std::vector< std::vector< RVec > > traj, int b_frame,
+                      std::vector< std::vector< std::vector< double > > > crl, double crl_b, double e_rad, int tauE) {
+    graph.resize(traj.front().size());
+    for (int i = 0; i < traj.front().size(); i++) {
+        graph[i].resize(traj.front().size(), std::make_pair(0, -1));
+    }
+    rvec temp;
+    for (int i = 1; i <= tauE; i++) {
+        for (int j = 0; j < traj.front().size(); j++) {
+            for (int f = j; f < traj.front().size(); f++) {
+                copy_rvec(traj[b_frame][j], temp);
+                rvec_dec(temp, traj[b_frame][f]);
+                if (std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j])) >= crl_b && norm(temp) <= e_rad && std::abs(graph[j][f].first) < std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j]))) {
+                    if (std::abs(crl[i][j][f]) > std::abs(crl[i][f][j])) {
+                        graph[j][f].first = crl[i][j][f];
+                    } else {
+                        graph[j][f].first = crl[i][f][j];
+                    }
+                    graph[j][f].second = i;
+                }
+            }
+        }
+    }
+    std::vector< bool > graph_flags;
+    graph_flags.resize(traj.front().size(), true);
+    std::vector< int > a;
+    a.resize(0);
+    std::vector< std::pair< int, int > > b;
+    b.resize(0);
+    std::vector< int > que1, que2, que3;
+    for (int i = 0; i < traj.front().size(); i++) {
+        if (graph_flags[i]) {
+            s_graph.push_back(a);
+            s_graph_rbr.push_back(b);
+            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 < traj.front().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]);
+                }
+            }
+            s_graph.back() = que3;
+            for (int j = 0; j < que3.size(); j++) {
+                for (int k = 0; k < traj.front().size(); k++) {
+                    if (graph[que3[j]][k].second > -1) {
+                        s_graph_rbr.back().push_back(std::make_pair(que3[j], k));
+                    }
+                }
+            }
+            //std::cout << s_graph.back().size() << "   ";
         }
     }
-#pragma omp barrier
 }
 
 bool myfunction (const std::pair< int, double > i, const std::pair< int, double > j) {
     return i.second < j.second;
 }
 
+void graph_back_bone_evaluation(std::vector< std::vector < std::pair< int, int > > > &rout_n, int indxSize,
+                                std::vector< std::vector< std::pair< double, int > > > graph, std::vector< std::vector< int > > s_graph, std::vector< std::vector< std::pair< int, int > > > s_graph_rbr) {
+    std::vector< double > key;
+    std::vector< int > path;
+    std::vector< std::pair< int, double > > que;
+    std::vector< std::pair< int, int > > a;
+    int v;
+    for (int i = 0; i < s_graph.size(); i++) {
+        key.resize(0);
+        path.resize(0);
+        que.resize(0);
+        v = 0;
+        if (s_graph[i].size() > 2) {
+            key.resize(indxSize, 2);
+            path.resize(indxSize, -1);
+            key[s_graph[i][0]] = 0;
+            for (int j = 0; j < s_graph[i].size(); j++) {
+                que.push_back(std::make_pair(s_graph[i][j], key[s_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 < s_graph_rbr[i].size(); j++) {
+                    int u = -1;
+                    if (s_graph_rbr[i][j].first == v) {
+                        u = s_graph_rbr[i][j].second;
+                    } else if (s_graph_rbr[i][j].second == v) {
+                        u = s_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);
+                    }
+                }
+            }
+            a.resize(0);
+            rout_n.push_back(a);
+            for (int j = 0; j < indxSize; j++) {
+                if (path[j] != -1) {
+                    rout_n.back().push_back(std::make_pair(j, path[j]));
+                }
+            }
+        }
+    }
+}
+
 /*! \brief
  * \ingroup module_trajectoryanalysis
  */
@@ -261,7 +463,7 @@ class Domains : public TrajectoryAnalysisModule
         std::string                                                 fnNdx_;
         SelectionList                                               sel_;
 
-        std::vector< std::vector< int > >                           domains;
+
         std::vector< std::vector< int > >                           domains_local;
         std::vector< std::vector< RVec > >                          trajectory;
         std::vector< std::vector< RVec > >                          frankenstein_trajectory;
@@ -271,6 +473,7 @@ class Domains : public TrajectoryAnalysisModule
         int                                                         frames              = 0;
         int                                                         basic_frame         = 0;
         int                                                         tau                 = 0;
+        int                                                         graph_tau           = 0;
         double                                                      crl_border          = 0;
         double                                                      eff_rad             = 1.5;
         real                                                        **w_rls;
@@ -301,6 +504,10 @@ Domains::initOptions(IOptionsContainer          *options,
     options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
                             .store(&fnNdx_).defaultBasename("domains")
                             .description("Index file from the domains"));
+    // Add option for graph_tau constant
+    options->addOption(gmx::IntegerOption("Gtau")
+                            .store(&graph_tau)
+                            .description("number of frames for graph to see for time travel"));
     // Add option for tau constant
     options->addOption(gmx::IntegerOption("tau")
                             .store(&tau)
@@ -333,6 +540,7 @@ void
 Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
                              const t_trxframe                       &fr)
 {
+    //std::vector< std::vector< int > > domains;
     ConstArrayRef< int > atomind  = sel_[0].atomIndices();
     index.resize(0);
     for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
@@ -344,85 +552,13 @@ Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings       &settings,
     for (int i = 0; i < index.size(); i++) {
         trajectory.back()[i] = fr.x[index[i]];
     }
-    domains.resize(sel_.size());
-    for (int i = 0; i < sel_.size(); i++) {
-        if (sel_.size() != 1 && i == 0) {
-            continue;
-        }
-        atomind  = sel_[i].atomIndices();
-        for (ConstArrayRef<int>::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) {
-            domains[i].push_back(*ai);
-        }
-    }
-    for (int i = 0; i < domains.size(); i++) {
-        for (int j = 0; j < domains[i].size(); j++) {
-            int k = 0;
-            while (index[k] != domains[i][j]) {
-                k++;
-            }
-            domains[i][j] = k;
-        }
-    }
-    for (int i = 0; i < domains.size(); i++) {
-        if (domains[i].size() >= 2) {
-            domains_local.push_back(domains[i]);
-            for (int k = 0; k < domains[i].size(); k++) {
-                for (int j = i + 1; j < domains.size(); j++) {
-                    for (int x = 0; x < domains[j].size(); x++) {
-                        if (domains[j][x] == domains[i][k]) {
-                            domains[j].erase(domains[j].begin() + x);
-                        }
-                    }
-                }
-            }
-        }
-    }
-    domains.resize(0);
-    domains = domains_local;
-    std::vector< bool >                 flags;
-    flags.resize(index.size(), true);
-    for (int i = 0; i < domains.size(); i++) {
-        for (int j = 0; j < domains[i].size(); j++) {
-            flags[domains[i][j]] = false;
-        }
-    }
-    int a;
-    rvec temp;
-    for (int i = 0; i < index.size(); i++) {
-        if (flags[i]) {
-            float dist = 90000001;
-            for (int j = 0; j < domains.size(); j++) {
-                for (int k = 0; k < domains[j].size(); k++) {
-                    rvec_sub(trajectory.back()[i], trajectory.back()[domains[j][k]], temp);
-                    if (norm(temp) <= dist) {
-                        dist = norm(temp);
-                        a = j;
-                    }
-                }
-            }
-            domains_local[a].push_back(i);
-            flags[i] = false;
-        }
-    }
+
+    domain_chopping(sel_, index, domains_local, trajectory.back());
+
     frankenstein_trajectory.resize(trajectory.size());
     frankenstein_trajectory.back() = trajectory.back();
     trajectory.resize(trajectory.size() + 1);
 
-    /*snew(w_rls, domains_local.size() + 1);
-    for (int i = 0; i < domains_local.size(); i++) {
-        snew(w_rls[i], index.size());
-        for (int j = 0; j < index.size(); j++) {
-            w_rls[i][j] = 0;
-        }
-        for (int j = 0; j < domains_local[i].size(); j++) {
-            w_rls[i][domains_local[i][j]] = 1;
-        }
-    }
-    snew(w_rls[domains_local.size()], index.size());
-    for (int i = 0; i < index.size(); i++) {
-        w_rls[domains_local.size()][i] = 1;
-    }*/
-
     w_rls2.resize(domains_local.size() + 1);
     for (int i = 0; i < domains_local.size(); i++) {
         w_rls2[i].resize(domains_local[i].size());
@@ -458,30 +594,6 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
         for (int j = 0; j < domains_local[i].size(); j++) {
             frankenstein_trajectory.back()[domains_local[i][j]] = traj[domains_local[i][j]];
         }
-        /*rvec *basic, *traj;
-        snew(basic, index.size());
-        for (int k = 0; k < index.size(); k++) {
-            copy_rvec(trajectory[basic_frame][k].as_vec(), basic[k]);
-        }
-        snew(traj, index.size());
-        for (int k = 0; k < index.size(); k++) {
-            copy_rvec(trajectory.back()[k].as_vec(), traj[k]);
-        }
-        reset_x(index.size(), NULL, index.size(), NULL, basic, w_rls[i]);
-        reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[i]);
-        do_fit(index.size(), w_rls[i], basic, traj);
-        for (int j = 0; j < index.size(); j++) {
-            if (w_rls[i][j] == 0) {
-                copy_rvec(basic[j], traj[j]);
-            }
-        }
-        reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[domains_local.size()]);
-
-        for (int j = 0; j < domains_local[i].size(); j++) {
-            frankenstein_trajectory.back()[domains_local[i][j]] = traj[domains_local[i][j]];
-        }
-        sfree(basic);
-        sfree(traj);*/
     }
     std::cout << "frame № " << frames + 1 <<" analysed\n";
     frames++;
@@ -495,87 +607,25 @@ Domains::finishAnalysis(int nframes)
     if (tau > 0) {
         k = tau;
     }
-    crltns.resize(k + 1);
-    for (int i = 0; i < crltns.size(); i++) {
-        crltns[i].resize(index.size());
-        for (int j = 0; j < crltns[i].size(); j++) {
-            crltns[i][j].resize(index.size(), 0);
-        }
-    }
-    std::cout << "Correlation's evaluation - start\n\n";
-    #pragma omp parallel shared(crltns)
-    {
-        #pragma omp for schedule(dynamic)
-        for (int i = m; i <= k; i += 1) {
-            rvec *medx, *medy, temp_zero;
-            snew(medx, index.size());
-            snew(medy, index.size());
-            temp_zero[0] = 0;
-            temp_zero[1] = 0;
-            temp_zero[2] = 0;
-            for (int i = 0; i < index.size(); i++) {
-                copy_rvec(temp_zero, medx[i]);
-                copy_rvec(temp_zero, medy[i]);
-            }
-            for (int j = 0; j < frankenstein_trajectory.size() - i - 1; j++) {
-                for (int f = 0; f < index.size(); f++) {
-                    rvec temp;
 
-                    rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j][f], temp);
-                    rvec_inc(medx[f], temp);
+    /*
+     *
+     *
+     *
+    */
 
-                    rvec_sub(frankenstein_trajectory[basic_frame][f], frankenstein_trajectory[j + i][f], temp);
-                    rvec_inc(medy[f], temp);
-                }
-            }
-            for (int j = 0; j < index.size(); j++) {
-                rvec temp;
-                real temp2 = frankenstein_trajectory.size() - 1;
-                temp2 -= i;
-                temp2 = 1 / temp2;
-
-                copy_rvec(medx[j], temp);
-                svmul(temp2, temp, medx[j]);
-                copy_rvec(medy[j], temp);
-                svmul(temp2, temp, medy[j]);
-            }
-            std::vector< std::vector< double > > a, b, c;
-            a.resize(index.size());
-            b.resize(index.size());
-            c.resize(index.size());
-            for (int j = 0; j < index.size(); j++) {
-                a[j].resize(index.size(), 0);
-                b[j].resize(index.size(), 0);
-                c[j].resize(index.size(), 0);
-            }
-            for (int j = 0; j < frankenstein_trajectory.size() - i - 1; j++) {
-                for (int f1 = 0; f1 < index.size(); f1++) {
-                    for (int f2 = 0; f2 < index.size(); f2++) {
-                        rvec temp1, temp2;
-                        rvec_sub(frankenstein_trajectory[basic_frame][f1], frankenstein_trajectory[j][f1], temp1);
-                        rvec_dec(temp1, medx[f1]);
+    std::cout << "Correlation's evaluation - start\n\n";
+    correlation_evaluation(frankenstein_trajectory, basic_frame, crltns, m, k);
+    std::cout << "Correlation's evaluation - end\n\n";
 
-                        rvec_sub(frankenstein_trajectory[basic_frame][f2], frankenstein_trajectory[j + i][f2], temp2);
-                        rvec_dec(temp2, 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]);
-                    }
-                }
-            }
-            for (int j = 0; j < index.size(); j++) {
-                for (int f = 0; f < index.size(); f++) {
-                    crltns[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
-                }
-            }
-            sfree(medx);
-            sfree(medy);
-        }
-    }
-    #pragma omp barrier
-    std::cout << "Correlation's evaluation - end\n\n";
-    for (int i1 = 0; i1 < 100; i1++) {
+    //number of corelations in the matrixes
+    /*for (int i1 = 0; i1 < 100; i1++) {
         int i5 = 0;
         for (int i2 = 1; i2 < crltns.size(); i2++) {
             for (int i3 = 0; i3 < crltns[i2].size(); i3++) {
@@ -590,128 +640,29 @@ Domains::finishAnalysis(int nframes)
         if (i1 % 10 == 0) {
             std::cout << "\n" ;
         }
-    }
+    }*/
+
+    std::cout << "graph evaluation: start\n";
+
     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() << "   ";
-        }
-    }
+
+
+    graph_calculation(graph, sub_graph, sub_graph_rbr, frankenstein_trajectory, basic_frame, crltns, crl_border, eff_rad, graph_tau/*k*/);
+
+
     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_DomainsNFit_2500_0.70_1.txt");
-    make_best_corrs_graphics(crltns, rout_new, index, "best_graphics_DomainsNFit_2500_0.70_1.txt");
+
+    std::cout << "graph evaluation: end\n";
+    std::cout << "routs evaluation: start\n";
+
+    graph_back_bone_evaluation(rout_new, index.size(), graph, sub_graph, sub_graph_rbr);
+
+    std::cout << "routs evaluation: end\n";
+
+    make_rout_file2(crl_border, index, rout_new, "Routs_DomainsNFit_5000_0.70_1_t.txt");
+    make_best_corrs_graphics(crltns, rout_new, index, "best_graphics_DomainsNFit_5000_0.70_1_t.txt");
     std::cout << "Finish Analysis - end\n\n";
 }