#include <iostream>
#include <omp.h>
#include <stdio.h>
-#include <chrono>
-
-
#include <iostream>
-#include <chrono>
-#include <omp.h>
-#include <thread>
#include <cstdlib>
#include <gromacs/utility/gmxomp.h>
#include <gromacs/trajectoryanalysis.h>
#include "gromacs/math/do_fit.h"
#include "gromacs/utility/smalloc.h"
-#include "newfit.h"
using namespace gmx;
Selection selec;
- std::vector< std::pair< int, int > > fitting_pairs;
std::vector< std::vector < RVec > > trajectory;
- std::vector< std::vector < RVec > > fitting_temp;
std::vector< double > noise;
std::vector< int > index;
- double epsi = 0.1;
- double fitting_prec = 0.003;
+ double epsi = 0.1; // should be selectable
+ int etalon_frame = 0; // should be selectable
int frames = 0;
- int iter_count = 0;
-
-
};
options->addOption(SelectionOption("select")
.store(&selec).required()
.description("Atoms that are considered as part of the excluded volume"));
- /*options->addOption(SelectionOption("fitting_prec")
- .store(&fitting_prec).required()
- .description("Recommended prec is 10^(-5)"));*/
-
-
/* options->addOption(DoubleOption("cutoff").store(&cutoff_)
.description("Cutoff for distance calculation (0 = no cutoff)"));
settings->setFlag(TrajectoryAnalysisSettings::efRequireTop);*/
}
-
-// -s '/home/toluk/Develop/samples/dusc_trna/EcDusC_tRNA_FMN.md_npt.non-sol.tpr' -f '/home/toluk/Develop/samples/dusc_trna/EcDusC_tRNA_FMN.md_npt.non-sol.xtc' -select 'name CA' -dt 500
-// -s '/home/toluk/Develop/samples/reca_rd_2008/reca_rd.md.non-sol.tpr' -f '/home/toluk/Develop/samples/reca_rd_2008/reca_rd.md.non-sol.xtc' -select 'name CA' -dt 500
-// -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'
-// -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'
-// /home/toluk/Develop/samples/banana_phone/
void
AnalysisTemplate::initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation & /*top*/)
void
AnalysisTemplate::finishAnalysis(int /*nframes*/)
{
-
real *w_rls;
rvec *x1, *x2, temp;
-
- snew(x1, trajectory[0].size());
- snew(x2, trajectory[0].size());
- snew(w_rls, trajectory[0].size());
-
std::vector< bool > flag;
- flag.resize(index.size(), true);
-
- for (int j = 0; j < index.size(); j++) {
- w_rls[j] = 1;
- }
-
- for (int i = 0; i < trajectory[0].size(); i++) {
- copy_rvec(trajectory[0][i].as_vec(), x1[i]);
- }
-
double noise_mid = 9999;
+ int count = 0;
- int count;
-
- FILE * pFile;
- pFile = fopen ("/home/toluk/Develop/testing/reca_old_core.out","w+");
-
- std::chrono::time_point<std::chrono::system_clock> start1, end1;
- start1 = std::chrono::system_clock::now();
+ snew(x1, index.size());
+ snew(x2, index.size());
+ snew(w_rls, index.size());
+ flag.resize(index.size(), true);
+ noise.resize(index.size(), 0);
while (noise_mid > epsi) {
- noise.resize(index.size(), 0);
noise_mid = 0;
+ count = 0;
+ noise.assign(index.size(), 0);
- for (int j = 0; j < index.size(); j++) {
- if (flag[j]) {
- w_rls[j] = 1;
+ for (int i = 0; i < index.size(); i++) {
+ if (flag[i]) {
+ w_rls[i] = 1;
} else {
- w_rls[j] = 0;
+ w_rls[i] = 0;
}
}
- for (int i = 1; i < frames; i++) {
+ for (int i = 0; i < frames; i++) {
for (int j = 0; j < trajectory[i].size(); j++) {
copy_rvec(trajectory[i][j].as_vec(), x2[j]);
- copy_rvec(trajectory[0][j].as_vec(), x1[j]);
+ copy_rvec(trajectory[etalon_frame][j].as_vec(), x1[j]); // i think its nessesary to do because of lots of reset_x
}
reset_x(index.size(), NULL, index.size(), NULL, x1, w_rls);
reset_x(index.size(), NULL, index.size(), NULL, x2, w_rls);
noise[j] += norm(temp);
}
}
+
for (int i = 0; i < index.size(); i++) {
- noise[i] /= (frames - 1);
+ noise[i] /= frames;
}
- count = 0;
-
for (int i = 0; i < index.size(); i++) {
if (flag[i]) {
count++;
}
noise_mid /= count;
- count = 0;
for (int i = 0; i < index.size(); i++) {
if (noise[i] > noise_mid) {
flag[i] = false;
}
- if (flag[i]) {
- count++;
- }
}
- std::cout << "noise_mid: " << noise_mid << " | number of atoms: " << count << "\n";
- std::cout << "\n\n finish iteration\n\n";
- iter_count++;
}
- end1 = std::chrono::system_clock::now();
- std::cout << "taken time: " << std::chrono::duration_cast<std::chrono::seconds>(end1-start1).count() << " seconds\n";
- fitting_pairs.resize(0);
for (int i = 0; i < noise.size(); i++) {
if (noise[i] <= epsi) {
- fprintf(pFile, "%d ", i + 1);
+ //fprintf(pFile, "%d ", i + 1);
+ // numbers from index, +1 for starting from 1 not 0
}
}
- fclose(pFile);
sfree(x1);
sfree(x2);
sfree(w_rls);
-
-
-
-
-/*
- // USING NEW FIT
-
- std::cout << "analys start\n";
- fitting_pairs.resize(index.size());
- for (int i = 0; i < index.size(); i++) {
- fitting_pairs[i] = std::make_pair (i, i);
- }
- double noise_mid = 9999;
-
- FILE * pFile;
- pFile = fopen ("/home/toluk/Develop/testing/reca_old_core.out","w+");
-
- std::chrono::time_point<std::chrono::system_clock> start1, end1;
- start1 = std::chrono::system_clock::now();
- while (noise_mid > epsi) {
- noise.resize(index.size(), 0);
- fitting_temp = trajectory;
- noise_mid = 0;
-
- for (int i = 1; i < frames; i++) {
- NewFit::new_fit(trajectory.front(), fitting_temp[i], fitting_pairs, fitting_prec);
- std::cout << i << " / " << frames << "\n";
- for (int j = 0; j < index.size(); j++) {
- noise[j] += norm(trajectory.front()[j] - fitting_temp[i][j]);
- }
- }
- 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: " << noise_mid << " | number of atoms: " << fitting_pairs.size() << "\n";
- std::cout << "\n\n finish iteration\n\n";
- iter_count++;
- }
- end1 = std::chrono::system_clock::now();
- std::cout << "taken time: " << std::chrono::duration_cast<std::chrono::seconds>(end1-start1).count() << " seconds\n";
-
- fitting_pairs.resize(0);
- for (int i = 0; i < noise.size(); i++) {
- if (noise[i] <= epsi) {
- fprintf(pFile, "%d ", i + 1);
- }
- }
- fclose(pFile);
-
- //60 61 62 73 74 75 76 188 189 215 217
- // 24 35 38 39 44 45 46 53 54 55 56 59 60 61 62 63 73 75 76 77 78 79 80 81 82 90 91 92 93 98 114 115 116 120 121 123 124 127 138 139 140 141 142 143 147 148 149 171 172 173 174 176 185 186 187 188 189 190 191 213 214 215 216 217 218 219 220 221 222 223 224 225 226 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
- // fitting constant chart calculation
-
-*/
-
-
- // NEW FIT CONSTANT TESTING
-/*
- fitting_pairs.resize(index.size());
- for (int i = 0; i < index.size(); i++) {
- fitting_pairs[i] = std::make_pair (i, i);
- }
- std::chrono::time_point<std::chrono::system_clock> start1, end1;
- FILE * pFile;
- pFile = fopen ("/home/toluk/Develop/testing/output.out","w+");
- double fitting_prec = 1;
- double temp1, temp2;
- std::vector < RVec > fit_temp;
- fprintf (pFile, "starting dist %f\n", NewFit::dist(trajectory[0], trajectory[900], fitting_pairs)) / fitting_pairs.size();
- fprintf (pFile, "fitting prec constant | new dist | taken time\n");
- for (int i = 1; i <= 10; i++) {
- temp1 = fitting_prec / 10;
- for (int j = 0; j <= 9; j++) {
- temp2 = fitting_prec - temp1 * j;
- fit_temp.resize(0);
- fit_temp = trajectory[900];
- start1 = std::chrono::system_clock::now();
- NewFit::new_fit(trajectory[0], fit_temp, fitting_pairs, temp2);
- end1 = std::chrono::system_clock::now();
- std::cout << NewFit::dist(trajectory[0], fit_temp, fitting_pairs) / fitting_pairs.size() << " | " << std::chrono::duration_cast<std::chrono::milliseconds>(end1-start1).count() << " | " << temp2 << "\n";
- fprintf (pFile, "%25.10f %25.10f %10ld millisecconds\n", temp2, NewFit::dist(trajectory[0], fit_temp, fitting_pairs) / fitting_pairs.size(), std::chrono::duration_cast<std::chrono::milliseconds>(end1-start1).count());
- }
- fitting_prec /= 10;
- }
- fclose (pFile);
-*/
-
-
-
-
-
-
-
-
-
- // NEW / OLD FIT COMPARISON
-/*
-
-
- std::chrono::time_point<std::chrono::system_clock> start1, end1, start2, end2;
- FILE * pFile1, * pFile2, * pFile3;
- pFile1 = fopen ("/home/toluk/Develop/testing/reca_old_old_fit.out","w+");
- pFile2 = fopen ("/home/toluk/Develop/testing/reca_old_new_fit.out","w+");
- pFile3 = fopen ("/home/toluk/Develop/testing/reca_old_comparison.out","w+");
- double fitting_prec = 0.00001;
- real *w_rls;
- rvec *x1, *x2;
- bool flag;
- snew(w_rls, index.size());
- std::vector < std::pair< int, int > > fit_pairs;
- std::vector < int > a;
- snew(x1, trajectory[0].size());
- snew(x2, trajectory[0].size());
- fit_pairs.resize(0);
- a.resize(0);
-
- for (int i1 = 0; i1 < trajectory.size(); i1++) {
- for (int i2 = 0; i2 < trajectory[i1].size(); i2++) {
- trajectory[i1][i2] = trajectory[i1][i2].scale(10);
- }
- NewFit::center_frame(trajectory[i1], a);
- }
-
- for (int i = 0; i < trajectory[0].size(); i++) {
- fit_pairs.push_back(std::make_pair(i, i));
- copy_rvec(trajectory[0][i].as_vec(), x1[i]);
- }
-
- for (int i = 0; i < index.size(); i++) {
- w_rls[i] = 1;
- }
-
-
-
- double start_dist = 0, d1, d2, d3;
- fprintf (pFile3, "| frame | dist stt | dist old | dist new | dist dif |\n");
- for (int i = 0; i < frames; i++) {
- start_dist = NewFit::dist(trajectory[0], trajectory[i], fit_pairs) / index.size();
-
- for (int i1 = 0; i1 < trajectory[0].size(); i1++) {
- copy_rvec(trajectory[0][i1].as_vec(), x1[i1]);
- copy_rvec(trajectory[i][i1].as_vec(), x2[i1]);
- }
-
- start1 = std::chrono::system_clock::now();
- NewFit::new_fit(trajectory[0], trajectory[i], fit_pairs, fitting_prec);
- end1 = std::chrono::system_clock::now();
-
- std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(end1-start1).count() << "\n";
-
- flag = true;
- start2 = std::chrono::system_clock::now();
-
- try {
- reset_x(index.size(), NULL, index.size(), NULL, x1, w_rls);
- reset_x(index.size(), NULL, index.size(), NULL, x2, w_rls);
- do_fit(index.size(), w_rls, x1, x2);
- }
- catch (...) {
- flag = false;
- }
- end2 = std::chrono::system_clock::now();
-
- std::cout << std::chrono::duration_cast<std::chrono::milliseconds>(end2-start2).count() << "\n\n";
-
- d1 = NewFit::dist_old_format(x1, x2, fit_pairs) / index.size();
- d2 = NewFit::dist(trajectory[0], trajectory[i], fit_pairs) / index.size();
- d3 = d2 - d1;
-
-
-
- fprintf(pFile3, "%14d %16f %16f %16f %16f\n", i, start_dist, d1, d2, d3);
- printf("%5d %13f %13f %13f %13f\n\n", i, start_dist, d1, d2, d3);
- if (i != frames - 1) {
- fprintf(pFile1, "MODEL %d\n", i + 1);
- fprintf(pFile2, "MODEL %d\n", i + 1);
- for (int i1 = 0; i1 < index.size(); i1++) {
- fprintf(pFile1, "ATOM %5d CA CA %5d %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s \n", i1, i1, x2[i1][0], x2[i1][1], x2[i1][2], 1.0, 20.0, " ", " ");
- fprintf(pFile2, "ATOM %5d CA CA %5d %8.3f%8.3f%8.3f%6.2f%6.2f %4s%2s \n", i1, i1, trajectory[i][i1][0], trajectory[i][i1][1], trajectory[i][i1][2], 1.0, 20.0, " ", " ");
- }
- fprintf(pFile1, "END\n");
- fprintf(pFile2, "END\n");
- }
- }
- fclose (pFile1);
- fclose (pFile2);
- fclose (pFile3);
-*/
-
- // RAND TESTING
-/*
- printf("\n\n\n\n\n\n\n\n\n\n");
-
-
- std::vector< RVec > a, b;
- std::vector < std::pair< int, int > > pairs;
- rvec *x1, *x2;
- real *w_rls;
- int k = 100;
- a.resize(k);
- b.resize(k);
- snew(x1, k);
- snew(x2, k);
- snew(w_rls, k);
- for (int i = 0; i < k; i++) {
- pairs.push_back(std::make_pair(i, i));
- w_rls[i] = 1;
- }
- RVec temp;
- while (true) {
- for (int i = 0; i < k; i++) {
- temp[0] = (float)(rand() % 10000) / 100;
- temp[1] = (float)(rand() % 10000) / 100;
- temp[2] = (float)(rand() % 10000) / 100;
- a[i] = temp;
- copy_rvec(a[i].as_vec(), x1[i]);
-
- temp[0] = (float)(rand() % 10000) / 100;
- temp[1] = (float)(rand() % 10000) / 100;
- temp[2] = (float)(rand() % 10000) / 100;
- b[i] = temp;
- copy_rvec(b[i].as_vec(), x2[i]);
- }
-
- NewFit::new_fit(a, b, pairs, 0.00001);
- std::cout << NewFit::dist(a, b, pairs) / k << " ";
-
- reset_x(k, NULL, k, NULL, x1, w_rls);
- reset_x(k, NULL, k, NULL, x2, w_rls);
- do_fit(k, w_rls, x1, x2);
- std::cout << NewFit::dist_old_format(x1, x2, pairs) / k << "\n";
- }
-
- printf("\n\n\n\n\n\n\n\n\n\n");
-
-*/
-
-
-
- std::cout << "analys finish\n";
}
void
AnalysisTemplate::writeOutput()
{
- /*std::cout << "number of iterations: " << iter_count << "\n" << "noise constant: " << epsi << "\n" << "atoms' numbers:\n";
- for (int i = 0; i < fitting_pairs.size(); i++) {
- std::cout << fitting_pairs[i].first << " ";
- }
- std::cout << "\n\n";
- std::cout << "atoms' noise values:\n";
- for (int i = 0; i < fitting_pairs.size(); i++) {
- std::cout << noise[fitting_pairs[i].first] << " ";
- }*/
+
}
/*! \brief
biod.pnpi.spb.ru / alexxy / gromacs-rcore.git / commitdiff