From 903e450211b68c3e83ff7bec18852634c7c9b223 Mon Sep 17 00:00:00 2001 From: Anatoly Titov Date: Thu, 6 Oct 2016 03:14:29 +0300 Subject: [PATCH] added new algorythm --- src/rcore.cpp | 135 +++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 122 insertions(+), 13 deletions(-) diff --git a/src/rcore.cpp b/src/rcore.cpp index d67df70..a80e6e7 100644 --- a/src/rcore.cpp +++ b/src/rcore.cpp @@ -38,6 +38,12 @@ #include +#include +#include +#include +#include + +#include #include #include "new_fit.h" @@ -79,11 +85,13 @@ class AnalysisTemplate : public TrajectoryAnalysisModule Selection selec; std::vector< std::pair< int, int > > fitting_pairs; std::vector< std::vector < RVec > > trajectory; - std::vector< int > noise; + std::vector< std::vector < RVec > > fitting_temp; + std::vector< long double > noise; std::vector< int > index; - const long double epsi = 0.15; - const long double fitting_prec = 0.0000001; + long double epsi = 0.15; + long double fitting_prec = 0.0001; int frames = 0; + int iter_count = 0; }; @@ -136,19 +144,20 @@ AnalysisTemplate::initOptions(IOptionsContainer *options, settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);*/ } -//domains -s '/home/toluk/Рабочий стол/reca_rd_2008/reca_rd.md.non-sol.tpr' -f '/home/toluk/Рабочий стол/reca_rd_2008/reca_rd.md.non-sol.xtc' -select 'name CA' -dt 100 +// -s '/home/toluk/Рабочий стол/reca_rd_2008/reca_rd.md.non-sol.tpr' -f '/home/toluk/Рабочий стол/reca_rd_2008/reca_rd.md.non-sol.xtc' -select 'name CA' -dt 500 +// -s '/home/toluk/Data/dusc_trna/EcDusC_tRNA_FMN.md_npt.non-sol.tpr' -f '/home/toluk/Data/dusc_trna/EcDusC_tRNA_FMN.md_npt.non-sol.xtc' -select 'name CA' -dt 500 void AnalysisTemplate::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation & /*top*/) { - std::cout << "select start\n"; + //std::cout << "select start\n"; index.resize(0); ConstArrayRef< int > atomind = selec.atomIndices(); for (ConstArrayRef< int >::iterator ai = atomind.begin(); (ai < atomind.end()); ai++) { index.push_back(*ai); } - std::cout << "select finish\n"; + //std::cout << "select finish\n"; } @@ -156,14 +165,14 @@ void AnalysisTemplate::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata) { - std::cout << "trajectory start\n"; + //std::cout << "trajectory start\n"; 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++; - std::cout << "trajectory finish\n"; + //std::cout << "trajectory finish\n"; } @@ -171,30 +180,125 @@ void AnalysisTemplate::finishAnalysis(int /*nframes*/) { std::cout << "analys start\n"; + //std::cout << "\n"; + fitting_pairs.resize(0); for (int i = 0; i < index.size(); i++) { fitting_pairs.push_back(std::make_pair (i, i)); } - long double noise_mid = 9000; - while (noise_mid > epsi) { + //std::cout << "001\n"; + long double noise_mid = 1; + bool flag = true; + std::vector< std::pair< int, int > > get_out; + //freopen("/home/toluk/Data/dusc_trna/output.txt", "w+", stdout); + + /* + + while (flag) { + std::cout << "start iteration\n"; + noise.resize(index.size(), 0); + fitting_temp = trajectory; + for (int i = 1; i < frames; i++) { - new_fit(trajectory[0], trajectory[i], fitting_pairs, index.size(), fitting_prec); + new_fit(trajectory[0], fitting_temp[i], fitting_pairs, index.size(), fitting_prec); + if (i % 50 == 0) { + std::cout << i << " / " << frames << "\n"; + } for (int j = 0; j < fitting_pairs.size(); j++) { - noise[fitting_pairs[j].first] += norm(trajectory[0][fitting_pairs[j].first] - trajectory[i][fitting_pairs[j].first]); + noise[fitting_pairs[j].first] += norm(trajectory[0][fitting_pairs[j].first] - fitting_temp[i][fitting_pairs[j].first]); } } + + flag = noise_mid; noise_mid = 0; - for (int i = 0; i < fitting_pairs.size(); i++) { + for (int i = 0; i < index.size(); i++) { noise[fitting_pairs[i].first] /= (frames - 1); noise_mid += noise[fitting_pairs[i].first]; } noise_mid /= fitting_pairs.size(); + + get_out = fitting_pairs; + fitting_pairs.resize(0); + for (int i = 0; i < index.size(); i++) { + fitting_pairs.push_back(std::make_pair (i, i)); + } + for (int i = 0; i < fitting_pairs.size(); i++) { if (noise[fitting_pairs[i].first] > noise_mid) { fitting_pairs.erase(fitting_pairs.begin() + i); i--; } } + + if (get_out == fitting_pairs) { + flag = false; + } + + std::cout << noise_mid << " | " << fitting_pairs.size() << "\n"; + for (int i = 0; i < fitting_pairs.size(); i++) { + std::cout << fitting_pairs[i].first << " "; + } + std::cout << "\n\n"; + for (int i = 0; i < get_out.size(); i++) { + std::cout << get_out[i].first << " "; + } + for (int i = 0; i < fitting_pairs.size(); i++) { + std::cout << noise[get_out[i].first] << " "; + } + std::cout << "\n\n"; + std::cout << "finish iteration\n"; + iter_count++; } + + */ + + noise_mid = 9999; + + while (noise_mid > epsi) { + std::cout << "start iteration\n"; + noise.resize(index.size(), 0); + fitting_temp = trajectory; + + for (int i = 1; i < frames; i++) { + new_fit(trajectory[0], fitting_temp[i], fitting_pairs, index.size(), fitting_prec); + if (i % 50 == 0) { + std::cout << i << " / " << frames << "\n"; + } + for (int j = 0; j < index.size(); j++) { + noise[j] += norm(trajectory[0][j] - fitting_temp[i][j]); + } + } + + noise_mid = 0; + for (int i = 0; i < index.size(); i++) { + noise[i] /= (frames - 1); + } + for (int i = 0; i < fitting_pairs.size(); i++) { + noise_mid += noise[fitting_pairs[i].first]; + } + + noise_mid /= fitting_pairs.size(); + + for (int i = 0; i < fitting_pairs.size(); i++) { + if (noise[fitting_pairs[i].first] > noise_mid) { + fitting_pairs.erase(fitting_pairs.begin() + i); + i--; + } + } + + std::cout << noise_mid << " | " << fitting_pairs.size() << "\n"; + std::cout << "finish iteration\n"; + iter_count++; + } + + fitting_pairs.resize(0); + for (int i = 0; i < noise.size(); i++) { + if (noise[i] <= epsi) { + fitting_pairs.push_back(std::make_pair (i, i)); + } + } + + + std::cout << "analys finish\n"; } @@ -203,9 +307,14 @@ void AnalysisTemplate::writeOutput() { std::cout << "output start\n"; + std::cout << "number of iterations: " << iter_count << "\n"; for (int i = 0; i < fitting_pairs.size(); i++) { std::cout << fitting_pairs[i].first << " "; } + std::cout << "\n"; + for (int i = 0; i < fitting_pairs.size(); i++) { + std::cout << noise[fitting_pairs[i].first] << " "; + } std::cout << "output finish\n"; } -- 2.22.0