From dc294f38ce355c2c8f972b2429648cbd7253e6e9 Mon Sep 17 00:00:00 2001 From: Anatoly Titov Date: Mon, 27 Feb 2017 15:18:11 +0300 Subject: [PATCH] changed the way frames are analyzed --- src/domains.cpp | 181 +++++++++++++++--------------------------------- 1 file changed, 54 insertions(+), 127 deletions(-) diff --git a/src/domains.cpp b/src/domains.cpp index 0ffedb2..ac6e380 100644 --- a/src/domains.cpp +++ b/src/domains.cpp @@ -61,7 +61,7 @@ using namespace gmx; using gmx::RVec; struct node { - int n; + short int n; RVec r; bool yep; }; @@ -242,18 +242,24 @@ class Domains : public TrajectoryAnalysisModule private: std::string fnNdx_; + std::vector< std::vector< std::vector< node > > > graph; + std::vector< std::vector< int > > domains; std::vector< std::vector< int > > domsizes; + std::vector< int > index; std::vector< int > numbers; std::vector< std::vector < RVec > > trajectory; Selection selec; int frames = 0; - int etalon_frame = 0; // should be selectable int domain_min_size = 5; // should be selectable + + real **w_rls; + int bone; double delta = 0.90; //0.95 // should be selectable double epsi = 0.15; //0.3 колебания внутри домена // should be selectable + int domains_ePBC; // Copy and assign disallowed by base. }; @@ -284,10 +290,6 @@ Domains::initOptions(IOptionsContainer *options, .store(&fnNdx_).defaultBasename("domains") .description("Index file from the domains")); // Add option for etalon_frame constant - options->addOption(gmx::IntegerOption("etalon_frame") - .store(&etalon_frame) - .description("basic frame to base evaluation on")); - // Add option for etalon_frame constant options->addOption(gmx::IntegerOption("dms") .store(&domain_min_size) .description("minimum domain size")); @@ -327,40 +329,15 @@ Domains::initAfterFirstFrame(const TrajectoryAnalysisSettings &settings, for (ConstArrayRef::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) { index.push_back(*ai); } - trajectory.resize(frames + 1); - trajectory[frames].resize(selec.atomCount()); + trajectory.resize(2); + trajectory[0].resize(selec.atomCount()); for (int i = 0; i < selec.atomCount(); i++) { - trajectory[frames][i] = fr.x[index[i]]; + trajectory[0][i] = fr.x[index[i]]; } - frames++; - graph.resize(0); -} - -void -Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, - TrajectoryAnalysisModuleData *pdata) -{ - trajectory.resize(frames + 1); - trajectory[frames].resize(index.size()); - for (int i = 0; i < index.size(); i++) { - trajectory[frames][i] = fr.x[index[i]]; - } - frames++; -} -//domains -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA' -//domains -s '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.tpr' -f '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.10th.xtc' -select 'name CA' -//domains -s '/home/toluk/Data/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Data/reca_rd/reca_rd.mono.xtc' - select 'name CA' - -void -Domains::finishAnalysis(int nframes) -{ - frames--; - int bone = index.size() - domain_min_size + 1; + bone = index.size() - domain_min_size + 1; graph.resize(bone); - real **w_rls; - snew(w_rls, bone); for (int i = 0; i < bone; i++) { snew(w_rls[i], index.size()); @@ -371,115 +348,69 @@ Domains::finishAnalysis(int nframes) w_rls[i][j] = 0; } } - } - - /* - rvec **etalon; - snew(etalon, bone); - for (int i = 0; i < bone; i++) { - snew(etalon[i], index.size()); - for (int j = 0; j < index.size(); j++) { - copy_rvec(trajectory[etalon_frame][j].as_vec(), etalon[i][j]); - } - } - */ - - for (int i = 0; i < bone; i++) { rvec *etalon; snew(etalon, index.size()); for (int j = 0; j < index.size(); j++) { - copy_rvec(trajectory[etalon_frame][j].as_vec(), etalon[j]); + copy_rvec(trajectory[0][j].as_vec(), etalon[j]); } reset_x(index.size(), NULL, index.size(), NULL, etalon, w_rls[i]); make_graph(index.size(), etalon, graph[i]); sfree(etalon); } + trajectory[1].resize(index.size()); +} - std::chrono::time_point start1, end1; - start1 = std::chrono::system_clock::now(); - - /* - rvec **traj; - snew(traj, bone); +void +Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, + TrajectoryAnalysisModuleData *pdata) +{ + for (int i = 0; i < index.size(); i++) { + trajectory[1][i] = fr.x[index[i]]; + } + frames++; - for (int i = 1; i < frames; i++) { + #pragma omp parallel + { + #pragma omp for schedule(dynamic) for (int j = 0; j < bone; j++) { - sfree(traj[j]); - snew(traj[j], index.size()); + rvec *etalon, *traj; + snew(etalon, index.size()); for (int k = 0; k < index.size(); k++) { - copy_rvec(trajectory[i][k].as_vec(), traj[j][k]); - } - } - #pragma omp parallel - { - #pragma omp for schedule(dynamic) - for (int j = 0; j < bone; j++) { - //#pragma omp ordered - //{ - reset_x(index.size(), NULL, index.size(), NULL, etalon[j], w_rls[j]); - reset_x(index.size(), NULL, index.size(), NULL, traj[j], w_rls[j]); - do_fit(index.size(), w_rls[j], etalon[j], traj[j]); - update_graph(graph[j], traj[j], epsi); - //} + copy_rvec(trajectory[0][k].as_vec(), etalon[k]); } - } - std::cout << i << " / " << frames << " processed\n"; - if (i % 25 == 0) { - end1 = std::chrono::system_clock::now(); - int seconds = 0, minutes = 0, hours = 0; - seconds = std::chrono::duration_cast(end1-start1).count(); - seconds = (float)seconds * (float)( (float)(frames - i) / (float)i ); - hours = seconds / 3600; - seconds -= hours * 3600; - minutes = seconds / 60; - seconds %= 60; - std::cout << "foretold time: " << hours << " hour(s) " << minutes << " minute(s) " << seconds << "second(s)\n"; - } - } - */ - - for (int i = 1; i < frames; i++) { - #pragma omp parallel - { - #pragma omp for schedule(dynamic) - for (int j = 0; j < bone; j++) { - rvec *etalon, *traj; - snew(etalon, index.size()); - for (int k = 0; k < index.size(); k++) { - copy_rvec(trajectory[etalon_frame][k].as_vec(), etalon[k]); - } - snew(traj, index.size()); - for (int k = 0; k < index.size(); k++) { - copy_rvec(trajectory[i][k].as_vec(), traj[k]); - } - reset_x(index.size(), NULL, index.size(), NULL, etalon, w_rls[j]); - reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[j]); - do_fit(index.size(), w_rls[j], etalon, traj); - update_graph(graph[j], traj, epsi); - sfree(etalon); - sfree(traj); + snew(traj, index.size()); + for (int k = 0; k < index.size(); k++) { + copy_rvec(trajectory[1][k].as_vec(), traj[k]); } - } - std::cout << i << " / " << frames << " processed\n"; - if (i % 25 == 0) { - end1 = std::chrono::system_clock::now(); - int seconds = 0, minutes = 0, hours = 0; - seconds = std::chrono::duration_cast(end1-start1).count(); - seconds = (float)seconds * (float)( (float)(frames - i) / (float)i ); - hours = seconds / 3600; - seconds -= hours * 3600; - minutes = seconds / 60; - seconds %= 60; - std::cout << "foretold time: " << hours << " hour(s) " << minutes << " minute(s) " << seconds << "second(s)\n"; + reset_x(index.size(), NULL, index.size(), NULL, etalon, w_rls[j]); + reset_x(index.size(), NULL, index.size(), NULL, traj, w_rls[j]); + do_fit(index.size(), w_rls[j], etalon, traj); + update_graph(graph[j], traj, epsi); + sfree(etalon); + sfree(traj); } } - std::cout << "finale cheking\n"; + std::cout << "frame: " << frames << " analyzed\n"; +} + +//domains -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -select 'name CA' +//domains -s '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.tpr' -f '/home/toluk/Develop/samples/banana_phone/pgk.md.non-sol.10th.xtc' -select 'name CA' +//domains -s '/home/toluk/Data/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Data/reca_rd/reca_rd.mono.xtc' - select 'name CA' + +void +Domains::finishAnalysis(int nframes) +{ + frames -= 1; + + std::cout << "final cheking\n"; check_domains(delta, frames, graph); + std::cout << "finding domains' sizes\n"; find_domain_sizes(graph, domsizes); + + std::cout << "finding domains\n"; std::vector< int > a; a.resize(0); - std::cout << "finding domains\n"; while (check_domsizes(domsizes, domain_min_size)) { domains.push_back(a); get_maxsized_domain(domains.back(), graph, domsizes); @@ -489,12 +420,8 @@ Domains::finishAnalysis(int nframes) } for (int i = 0; i < bone; i++) { sfree(w_rls[i]); - //sfree(etalon[i]); - //sfree(traj[i]); } sfree(w_rls); - //sfree(etalon); - //sfree(traj); } void -- 2.22.0