#include <omp.h>
+#include <iostream>
+#include <chrono>
+#include <omp.h>
+#include <thread>
+
+#include <gromacs/utility/gmxomp.h>
#include <gromacs/trajectoryanalysis.h>
#include "new_fit.h"
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;
};
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";
}
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";
}
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";
}
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";
}
biod.pnpi.spb.ru / alexxy / gromacs-rcore.git / blobdiff