#include <omp.h>
#include <thread>
#include <string>
+#include <algorithm>
#include <gromacs/trajectoryanalysis.h>
#include <gromacs/pbcutil/pbc.h>
{
FILE *file;
file = std::fopen(file_name, "w+");
- std::cout << "start printing\n";
for (int i = start; i < correlations.size(); i++) {
if (correlations.size() - start > 1) {
std::fprintf(file, "correlation between 'n' and 'n + %d' frames\n", i);
std::cout << "correlation between 'n' and 'n + " << i << "' frames\n";
}
for (int j = 0; j < correlations[i].size(); j++) {
- for (int k = 0; k < correlations[i][j].size(); k++) {
- //if (std::abs(correlations[i][j][k]) > 0.5) {
- std::fprintf(file, "%3.2f ", correlations[i][j][k]);
- //}
+ for (int f = 0; f < correlations[i][j].size(); f++) {
+ std::fprintf(file, "%3.2f ", correlations[i][j][f]);
}
std::fprintf(file, "\n");
}
}
}
-void make_rout_file(std::vector< int > dist, std::vector< std::vector< int > > rout, const char* file_name)
+void make_rout_file(std::vector< std::vector< int > > rout, const char* file_name)
{
FILE *file;
file = std::fopen(file_name, "w+");
- std::cout << "start printing\n";
for (int i = 0; i < rout.size(); i++) {
- //std::fprintf(file, "sum dist %d\n", dist[i]);
for (int j = 0; j < rout[i].size(); j++) {
std::fprintf(file, "%3d ", rout[i][j]);
}
if (flag) {
path.push_back(std::make_pair(start, depth_go));
}
-
- /*for (int i = 0; i < path.size(); i++) {
- std::cout << path[i].second << " ";
- }
- std::cout << " | " << depth << "\n";*/
-
for (int i = 0; i < graph[start].size(); i++) {
flag = true;
for (int j = 0; j < path.size(); j++) {
return (std::max(steps_old, steps));
}
+bool mysortfunc (std::vector<int> a, std::vector<int> b) { return (a.size() > b.size()); }
+
/*! \brief
* Class used to compute free volume in a simulations box.
*
if (flags[i]) {
int a, b, dist = 90000001;
rvec temp;
- //if (domains.size() != 2) { спорно как ту поступить
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);
}
}
}
- //}
domains_local[a].push_back(i);
flags[i] = false;
}
void
Domains::finishAnalysis(int nframes)
{
- //make_pdb_trajectory(frankenstein_trajectory, "/home/toluk/Develop/samples/reca_rd/for_prez_0.pdb");
+ make_pdb_trajectory(frankenstein_trajectory, "/home/toluk/Develop/samples/reca_rd/10k_frames.pdb");
+
+ //frankenstein_trajectory = trajectory;
+
std::vector< std::vector< std::vector< double > > > crltns;
- int k = 1000, m = 0;
+ int k = 2000, m = 0;
crltns.resize(k + 1);
for (int i = 0; i < crltns.size(); i++) {
crltns[i].resize(index.size());
}
}
- //frankenstein_trajectory = trajectory;
-
#pragma omp parallel
{
#pragma omp for schedule(dynamic)
for (int i = m; i <= k; i += 1) {
std::cout << " \n " << " k = " << i << " \n";
- /*std::vector< float > medx, medy;
- medx.resize(index.size(), 0);
- medy.resize(index.size(), 0);*/
-
rvec *medx, *medy, temp_zero;
snew(medx, index.size());
snew(medy, index.size());
copy_rvec(temp_zero, medy[i]);
}
- for (int j = 0; j < frankenstein_trajectory.size() - k; j++) {
+ for (int j = 0; j < frankenstein_trajectory.size() - i; j++) {
for (int f = 0; f < index.size(); f++) {
rvec temp;
- /*rvec_sub(frankenstein_trajectory[j][f], frankenstein_trajectory[j + 1][f], temp);
- medx[f] += norm(temp);
- rvec_sub(frankenstein_trajectory[j + k - 1][f], frankenstein_trajectory[j + k][f], temp);
- medy[f] += norm(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 + k][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++) {
- /*double temp2 = frankenstein_trajectory.size();
- temp2 -= k;
- medx[j] /= temp2;
- medy[j] /= temp2;*/
-
rvec temp;
real temp2 = frankenstein_trajectory.size();
- temp2 -= k;
+ temp2 -= i;
temp2 = 1 / temp2;
copy_rvec(medx[j], temp);
svmul(temp2, temp, medx[j]);
b[j].resize(index.size(), 0);
c[j].resize(index.size(), 0);
}
- for (int j = 0; j < frankenstein_trajectory.size() - k; j++) {
+ for (int j = 0; j < frankenstein_trajectory.size() - i; j++) {
for (int f1 = 0; f1 < index.size(); f1++) {
for (int f2 = f1; f2 < index.size(); f2++) {
rvec temp1, temp2;
-
- /*rvec_sub(frankenstein_trajectory[j][f1], frankenstein_trajectory[j + 1][f1], temp1);
- rvec_sub(frankenstein_trajectory[j + k - 1][f2], frankenstein_trajectory[j + k][f2], temp2);
- a[f1][f2] += (norm(temp1) - medx[f1]) * (norm(temp2) - medy[f2]);
- b[f1][f2] += (norm(temp1) - medx[f1]) * (norm(temp1) - medx[f1]);
- c[f1][f2] += (norm(temp2) - medy[f2]) * (norm(temp2) - medy[f2]);*/
-
rvec_sub(frankenstein_trajectory[basic_frame][f1], frankenstein_trajectory[j][f1], temp1);
rvec_dec(temp1, medx[f1]);
- rvec_sub(frankenstein_trajectory[basic_frame][f2], frankenstein_trajectory[j + k][f2], temp2);
+ 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]);
sfree(medy);
}
}
-// std::cout << "\noutput file\n";
-// make_correlation_file(crltns, "/home/toluk/Develop/samples/reca_rd/corrs.txt", m);
-// std::cout << "\n\n file done \n\n";
-/*
- std::cout << "\n\n\n\n";
- for (int i2 = 0; i2 < 40; i2++) {
- long int count_here = 0;
- for (int i = 0; i <= k; i++) {
- for (int j = 0; j < index.size(); j++) {
- for (int f = 0; f < index.size(); f++) {
- if (std::abs(crltns[i][j][f]) > (float)i2 / 100) {
- count_here++;
+ //make_correlation_file(crltns, "/home/toluk/Develop/samples/reca_rd/corrs_10k_frames.txt", m);
+
+ for (int i1 = 20; i1 < 70; i1++) {
+ int i5 = 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 (crltns[i2][i3][i4] >= (double)i1 / 100) {
+ i5++;
}
}
}
}
- std::cout << i2 << " " << count_here << "\n";
+ std::cout << i1 << " - " << i5 << "\n";
}
- std::cout << "\n\n\n\n";
-*/
- std::vector< std::vector< std::pair< int, int > > > graph;
+ std::cout << "\n\n\n";
+
+ std::vector< std::vector< int > > graph;
graph.resize(index.size());
- for (int i = 0; i < k; i++) {
+ for (int i = 1; i <= k; i++) {
for (int j = 0; j < index.size(); j++) {
for (int f = 0; f < index.size(); f++) {
- if (std::abs(crltns[i][j][f]) >= 0.05) {
- graph[j].push_back(std::make_pair(f, i));
+ if (std::abs(crltns[i][j][f]) >= 0.30) {
+ bool local_flag = true;
+ for (int f2 = 0; f2 < graph[j].size(); f2++) {
+ if (graph[j][f2] == f) {
+ local_flag = false;
+ }
+ }
+ if (local_flag) {
+ graph[j].push_back(f);
+ }
}
}
}
}
- for (int i = 0; i < index.size(); i++) {
- std::cout << i << "-" << graph[i].size() << " ";
- }
+ /*for (int i = 0; i < index.size(); i++) {
+ if (graph[i].size() > 0) {
+ std::cout << i << " - " << graph[i].size() << "\n";
+ }
+ }*/
+
std::cout << "\n";
std::vector< std::vector < int > > rout_old, rout_new;
- std::vector< int > length_old, length_new;
bool flag = true;
rout_new.resize(0);
- length_new.resize(0);
- length_old.resize(index.size(), 0);
rout_old.resize(index.size());
for (int i = 0; i < index.size(); i++) {
rout_old[i].resize(1, i);
std::cout << "building routs\n" ;
-
- std::vector< std::pair< int, int > > temp_a;
- std::vector< int > howlong;
- int temp_b = 0;
- howlong.resize(index.size(), 0);
-
-/*#pragma omp parallel
-{
- #pragma omp for schedule(dynamic)
- for (int i = 0; i < index.size(); i++) {
- temp_a.resize(1, std::make_pair(i, 0));
- //std::cout << "try out\n";
- temp_b = dive_in(graph, i, 0, frames, 0, 0, temp_a);
- howlong[i] = temp_b;
- std::cout << i << " " << temp_b << "\n";
- }
-}*/
-
- std::vector < int > flag3;
while (flag) {
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]);
- length_new.push_back(length_old[i]);
} else {
for (long int j = 0; j < graph[rout_old[i].back()].size(); j++) {
- if (length_old[i] + graph[rout_old[i].back()][j].second <= frames) {
- flag3 = rout_old[i];
- flag3.push_back(graph[rout_old[i].back()][j].first);
-
- bool flag2 = true;
+ bool flag2 = true, flag22 = true;;
+ flag3 = rout_old[i];
+ for (int f = 0; f < flag3.size(); f++) {
+ if (flag3[f] == graph[rout_old[i].back()][j]) {
+ flag22 = false;
+ }
+ }
+ 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);
- length_new.push_back(length_old[i] + graph[rout_old[i].back()][j].second);
flag = true;
}
}
if (rout_new.size() != 0) {
rout_old = rout_new;
rout_new.resize(0);
- length_old = length_new;
- length_new.resize(0);
- std::cout << length_old.size() << "\n";
+ std::cout << rout_old.size() << "\n";
}
}
+ std::sort(rout_old.begin(), rout_old.end(), mysortfunc);
+
+ std::cout << "\nfinished routs\n" ;
+
+ /*graph.resize(0);
+ graph.resize(index.size());
+ rout_new.resize(0);
+
for (int i = 0; i < rout_old.size(); i++) {
- for (int j = 0; j < rout_old[i].size(); j++) {
- std::cout << rout_old[i][j] << " ";
+ int lflag2 = false;
+ for (int j = 0; j < rout_old[i].size() - 1; j++) {
+ int lflag1 = true;
+ for (int f = 0; f < graph[rout_old[i][j]].size(); f++) {
+ if (graph[rout_old[i][j]][f] == rout_old[i][j + 1]) {
+ lflag1 = false;
+ }
+ }
+ if (lflag1) {
+ graph[rout_old[i][j]].push_back(rout_old[i][j + 1]);
+ lflag2 = true;
+ }
}
- std::cout << "\n";
- }
-
+ if (lflag2) {
+ rout_new.push_back(rout_old[i]);
+ }
+ }*/
+ 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_old[j][f] == rout_old[i][la]) {
+ la++;
+ }
+ }
+ if (la == rout_old[i].size()) {
+ lflag1 = false;
+ }
+ }
+ }
+ if (lflag1) {
+ rout_new.push_back(rout_old[i]);
+ }
+ }
- std::cout << "\nfinished routs\n" ;
- make_rout_file(length_old, rout_old, "/home/toluk/Develop/samples/reca_rd/routs.txt");
+ make_rout_file(rout_new, "routs.txt");
std::cout << "\n\n\n job done \n\n\n";
-
}
void
biod.pnpi.spb.ru / alexxy / gromacs-testing.git / commitdiff