From 0f24b1d1feaba85def2c2a060575cb2df8c9e0c6 Mon Sep 17 00:00:00 2001 From: Anatoly Date: Tue, 29 Oct 2019 17:47:39 +0300 Subject: [PATCH] still need to find THE bug probability is high its in the final part --- src/domains.cpp | 99 ++++++++++++++++++++++++++++++++----------------- 1 file changed, 65 insertions(+), 34 deletions(-) diff --git a/src/domains.cpp b/src/domains.cpp index ec4298e..1101940 100644 --- a/src/domains.cpp +++ b/src/domains.cpp @@ -62,6 +62,8 @@ using namespace gmx; using gmx::RVec; +const int tau_jump = 10; + double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) { return sqrt( pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) + pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) + @@ -316,24 +318,26 @@ bool check_domsizes(std::vector< std::vector< int > > cd_domsizes, int cd_domain return false; } -void print_domains(std::vector< std::vector< unsigned int > > pd_domains, std::vector< int > index, std::string fnNdx_, int dms, double epsi, double delta) { +void print_domains(std::vector< std::vector< std::vector< unsigned int > > > pd_domains, std::vector< int > index, std::string fnNdx_, int dms, double epsi, double delta, int window) { FILE *fpNdx_; fpNdx_ = std::fopen(fnNdx_.c_str(), "w+"); int write_count; - for (unsigned int i = 0; i < pd_domains.size(); i++) { - std::fprintf(fpNdx_, "[domain_№_%3d_%3d_%4.3f_%4.3f]\n", i + 1, dms, epsi, delta); - write_count = 0; - for (unsigned int j = 0; j < pd_domains[i].size(); j++) { - write_count++; - if (write_count > 20) { - write_count -= 20; - std::fprintf(fpNdx_, "\n"); + for (unsigned int i1 = 0; i1 < pd_domains.size(); i1++) { + for (unsigned int i2 = 0; i2 < pd_domains[i1].size(); i2++) { + std::fprintf(fpNdx_, "[domain-stPos-%3d-num-%3d-dms-%3d-epsi-%4.3f-delta-%4.3f]\n", static_cast< int >(i1) * window / tau_jump, i2 + 1, dms, epsi, delta); + write_count = 0; + for (unsigned int j = 0; j < pd_domains[i1][i2].size(); j++) { + write_count++; + if (write_count > 20) { + write_count -= 20; + std::fprintf(fpNdx_, "\n"); + } + std::fprintf(fpNdx_, "%5d ", index[pd_domains[i1][i2][j]] + 1); + std::printf("%5d ", index[pd_domains[i1][i2][j]] + 1); } - std::fprintf(fpNdx_, "%5d ", index[pd_domains[i][j]] + 1); - std::printf("%5d ", index[pd_domains[i][j]] + 1); + std::printf("\n\n"); + std::fprintf(fpNdx_,"\n\n"); } - std::printf("\n\n"); - std::fprintf(fpNdx_,"\n\n"); } std::fprintf(fpNdx_,"\n"); std::fclose(fpNdx_); @@ -378,7 +382,7 @@ class Domains : public TrajectoryAnalysisModule std::vector< std::vector< std::vector< std::vector< node > > > > graph; - std::vector< std::vector< unsigned int > > domains; + std::vector< std::vector< std::vector< unsigned int > > > domains; std::vector< std::vector< int > > domsizes; std::vector< int > index; @@ -387,9 +391,10 @@ class Domains : public TrajectoryAnalysisModule std::vector< RVec > reference; Selection selec; int frames = 0; + int window = 1000; // selectable int domain_min_size = 4; // selectable - std::vector< std::vector< std::pair< unsigned int, unsigned int > > > w_rls2; + std::vector< std::vector< std::pair< unsigned int, unsigned int > > > w_rls2; unsigned int bone; double delta = 0.95; // selectable double epsi = 0.10; // selectable @@ -461,7 +466,7 @@ Domains::initAnalysis(const TrajectoryAnalysisSettings &settings, } } - graph.resize(0); + graph.resize(1); bone = static_cast< unsigned int >(index.size() - static_cast< unsigned int >(domain_min_size) + 1); graph[0].resize(bone); @@ -480,6 +485,14 @@ Domains::initAnalysis(const TrajectoryAnalysisSettings &settings, void Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata) { + if ((frnr + 1) % (window / tau_jump) == 0) { + graph.resize(graph.size() + 1); + graph.back().resize(bone); + for (unsigned int i = 0; i < bone; i++) { + make_graph(index.size(), reference, graph.back()[i]); + } + } + for (unsigned int i = 0; i < index.size(); i++) { trajectory[i] = fr.x[index[i]]; } @@ -489,7 +502,9 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnal for (unsigned int j = 0; j < bone; j++) { std::vector < RVec > temp = trajectory; MyFitNew(reference, temp, w_rls2[j], 0); - update_graph(graph[0][j], temp, epsi); + for (unsigned int i = static_cast< unsigned int >(std::max(0, (frnr + 1) / (window / tau_jump) - tau_jump + 1)); i < graph.size(); i++) { + update_graph(graph[i][j], temp, epsi); + } } #pragma omp barrier std::cout << "frame: " << frames << " analyzed\n"; @@ -498,28 +513,44 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnal void Domains::finishAnalysis(int nframes) { + graph.resize(graph.size() - tau_jump); + std::cout << "graph size: " << graph.size() << "\n"; frames -= 1; - std::cout << "final cheking\n"; - check_domains(delta, frames, graph[0]); - - std::cout << "finding domains' sizes\n"; - find_domain_sizes(graph[0], domsizes); + domains.resize(graph.size()); - std::cout << "finding domains\n"; std::vector< unsigned int > a; a.resize(0); - while (check_domsizes(domsizes, domain_min_size)) { - domains.push_back(a); - if (DomainSearchingAlgorythm == 0) { - get_maxsized_domain(domains.back(), graph[0], domsizes); - } else if (DomainSearchingAlgorythm == 1) { - get_min_domain(domains.back(), graph[0], domsizes, domain_min_size); - } - std::cout << "new domain: " << domains.back().size() << " atoms\n"; - delete_domain_from_graph(graph[0], domains.back()); + + std::cout << "final cheking\n"; + + for (unsigned int i = 0; i < graph.size(); i++) { domsizes.resize(0); - find_domain_sizes(graph[0], domsizes); + domains[i].resize(0); + check_domains(delta, frames, graph[i]); + std::cout << "finding domains' sizes\n"; + find_domain_sizes(graph[i], domsizes); + + /*for (int j1 = 0; j1 < domsizes.size(); j1++) { + for (int j2 = 0; j2 < domsizes[j1].size(); j2++) { + if (domsizes[j1][j2] != 0) { + std::cout << domsizes[j1][j2] << " "; + } + } + }*/ + + while (check_domsizes(domsizes, domain_min_size)) { + domains[i].push_back(a); + if (DomainSearchingAlgorythm == 0) { + get_maxsized_domain(domains[i].back(), graph[i], domsizes); + } else if (DomainSearchingAlgorythm == 1) { + get_min_domain(domains[i].back(), graph[i], domsizes, domain_min_size); + } + std::cout << "new domain: " << domains[i].back().size() << " atoms\n"; + delete_domain_from_graph(graph[i], domains[i].back()); + domsizes.resize(0); + find_domain_sizes(graph[i], domsizes); + } } } @@ -527,7 +558,7 @@ void Domains::writeOutput() { std::cout << "making output file\n"; - print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta); // see function for details | numbers from index + print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta, window); // see function for details | numbers from index std::cout << "\n END \n"; } -- 2.22.0