std::fclose(file);
}
-void make_rout_file(float crl_border, std::vector< int > indx, std::vector< std::vector< std::pair< int, int > > > rout, const char* file_name)
-{
- FILE *file;
- file = std::fopen(file_name, "w+");
- std::fprintf(file, "correlations >= %0.2f\n\n", crl_border);
- for (int i = 0; i < rout.size(); i++) {
- for (int j = 0; j < rout[i].size(); j++) {
- std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.15\n", indx[rout[i][j].first]/* + 1*/, indx[rout[i][j].second]/* + 1*/);
- }
- std::fprintf(file, "\n\n");
- }
- std::fclose(file);
-}
-
-void make_best_corrs_graphics(std::vector< std::vector< std::vector< float > > > correlations,
- std::vector< std::vector< std::pair< int, int > > > rout_pairs,
- std::vector< int > indx,
- const char* file_name)
-{
- FILE *file;
- file = std::fopen(file_name, "w+");
- for (int i = 0; i < rout_pairs.size(); i++) {
- for (int j = 0; j < rout_pairs[i].size(); j++) {
- std::fprintf(file, "%3d %3d\n", indx[rout_pairs[i][j].first]/* + 1*/, indx[rout_pairs[i][j].second]/* + 1*/);
- for (int k = 0; k < correlations.size(); k++) {
- std::fprintf(file, "%3.5f ", correlations[k][rout_pairs[i][j].first][rout_pairs[i][j].second]);
- }
- std::fprintf(file, "\n");
- }
- }
- std::fclose(file);
-}
-
-void make_diffusion_file(const char* file_name, std::vector< float > D)
-{
- FILE *file;
- file = std::fopen(file_name, "w+");
- for (int i = 0; i < D.size(); i++) {
- std::fprintf(file, "%f ", D[i]);
- }
- std::fclose(file);
-}
-
-bool mysortfunc (std::vector< int > a, std::vector< int > b) {
- return (a.size() > b.size());
-}
-
-bool isitsubset (std::vector< int > a, std::vector< int > b) {
- if (b.size() == 0) {
- return true;
- }
- std::sort(a.begin(), a.end());
- std::sort(b.begin(), b.end());
- int k = 0;
- for (int i = 0; i < a.size(); i++) {
- if (b[k] == a[i]) {
- k++;
- }
- }
- if (k == b.size()) {
- return true;
- } else {
- return false;
- }
-}
-
-void correlation_evaluation(std::vector< std::vector< RVec > > traj, int b_frame, std::vector< std::vector< std::vector< float > > > &crl, int tauS, int tauE) {
- crl.resize(tauE + 1);
+void correlation_evaluation(std::vector< std::vector< RVec > > traj, int window_size, std::vector< std::vector< std::vector< float > > > &crl) {
+ crl.resize(traj.size() - window_size + 1);
for (int i = 0; i < crl.size(); i++) {
crl[i].resize(traj.front().size());
for (int j = 0; j < crl[i].size(); j++) {
crl[i][j].resize(traj.front().size(), 0);
}
}
+
RVec temp_zero;
temp_zero[0] = 0;
temp_zero[1] = 0;
//changed parrallel
//#pragma omp parallel shared(crl, temp_zero, d)
//{
- #pragma omp parallel for schedule(dynamic)
- for (int i = tauS; i <= tauE; i += 1) {
- std::vector< RVec > medx, medy;
- RVec temp1, temp2;
- float tmp;
- medx.resize(traj.front().size(), temp_zero);
- medy.resize(traj.front().size(), temp_zero);
+ #pragma omp parallel for schedule(dynamic)
+ for (int i = 0; i < crl.size(); i++) {
+ std::vector< RVec > medx, medy;
+ RVec temp1, temp2;
+ float tmp;
+ medx.resize(traj.front().size(), temp_zero);
+ medy.resize(traj.front().size(), temp_zero);
//#pragma omp for schedule(dynamic)
- for (int j = 0; j < traj.size() - i - 1; j++) {
- for (int f = 0; f < traj.front().size(); f++) {
- medx[f] += (traj[b_frame][f] - traj[j][f]);
- medy[f] += (traj[b_frame][f] - traj[j + i][f]);
- }
+ for (int j = i; j < i + window_size - 1; j++) {
+ for (int f = 0; f < traj.front().size(); f++) {
+ medx[f] += (traj[i][f] - traj[j][f]);
+ medy[f] += (traj[i][f] - traj[j][f]);
}
+ }
//#pragma omp barrier
//#pragma omp for schedule(dynamic)
- for (int j = 0; j < traj.front().size(); j++) {
- tmp = traj.size() - 1;
- tmp -= i;
- tmp = 1 / tmp;
+ for (int j = 0; j < traj.front().size(); j++) {
+ tmp = traj.size() - 1;
+ tmp -= i;
+ tmp = 1 / tmp;
//tmp = 1 / (traj.size() - 1 - i);
- medx[j] *= tmp;
- medy[j] *= tmp;
- }
+ medx[j] *= tmp;
+ medy[j] *= tmp;
+ }
//#pragma omp barrier
- std::vector< std::vector< float > > a, b, c;
- a.resize(traj.front().size(), d);
- b.resize(traj.front().size(), d);
- c.resize(traj.front().size(), d);
- for (int j = 0; j < traj.size() - i - 1; j++) {
- for (int f1 = 0; f1 < traj.front().size(); f1++) {
+ std::vector< std::vector< float > > a, b, c;
+ a.resize(traj.front().size(), d);
+ b.resize(traj.front().size(), d);
+ c.resize(traj.front().size(), d);
+ for (int j = i; j < i + window_size - 1; j++) {
+ for (int f1 = 0; f1 < traj.front().size(); f1++) {
//#pragma omp for schedule(dynamic)
- for (int f2 = 0; f2 < traj.front().size(); f2++) {
- temp1 = traj[b_frame][f1] - traj[j][f1] - medx[f1];
- temp2 = traj[b_frame][f2] - traj[j + i][f2] - medy[f2];
- a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
- b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
- c[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
- }
- //#pragma omp barrier
- }
- }
-
- //#pragma omp for schedule(dynamic)
- for (int j = 0; j < traj.front().size(); j++) {
- for (int f = 0; f < traj.front().size(); f++) {
- crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
+ for (int f2 = 0; f2 < traj.front().size(); f2++) {
+ temp1 = traj[i][f1] - traj[j][f1] - medx[f1];
+ temp2 = traj[i][f2] - traj[j][f2] - medy[f2];
+ a[f1][f2] += (temp1[0] * temp2[0] + temp1[1] * temp2[1] + temp1[2] * temp2[2]);
+ b[f1][f2] += (temp1[0] * temp1[0] + temp1[1] * temp1[1] + temp1[2] * temp1[2]);
+ c[f1][f2] += (temp2[0] * temp2[0] + temp2[1] * temp2[1] + temp2[2] * temp2[2]);
}
+ //#pragma omp barrier
}
- //#pragma omp barrier
- medx.resize(0);
- medy.resize(0);
- std::cout << i << " corr done\n";
}
- #pragma omp barrier
- //}
-}
-
-void graph_calculation(std::vector< std::vector< std::pair< float, int > > > &graph, std::vector< std::vector< int > > &s_graph, std::vector< std::vector< std::pair< int, int > > > &s_graph_rbr,
- std::vector< std::vector< RVec > > traj, int b_frame,
- std::vector< std::vector< std::vector< float > > > crl, float crl_b, float e_rad, int tauE) {
- graph.resize(traj.front().size());
- for (int i = 0; i < traj.front().size(); i++) {
- graph[i].resize(traj.front().size(), std::make_pair(0, -1));
- }
- RVec temp;
- for (int i = 1; i <= tauE; i++) {
+ //#pragma omp for schedule(dynamic)
for (int j = 0; j < traj.front().size(); j++) {
- for (int f = j; f < traj.front().size(); f++) {
- temp = traj[b_frame][j] - traj[b_frame][f];
- if (std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j])) >= crl_b && norm(temp) <= e_rad && std::abs(graph[j][f].first) < std::max(std::abs(crl[i][j][f]), std::abs(crl[i][f][j]))) {
- if (std::abs(crl[i][j][f]) > std::abs(crl[i][f][j])) {
- graph[j][f].first = crl[i][j][f];
- } else {
- graph[j][f].first = crl[i][f][j];
- }
- graph[j][f].second = i;
- }
- }
- }
- }
- std::cout << "crl analysed\n";
- std::vector< bool > graph_flags;
- graph_flags.resize(traj.front().size(), true);
- std::vector< int > a;
- a.resize(0);
- std::vector< std::pair< int, int > > b;
- b.resize(0);
- std::vector< int > que1, que2, que3;
- for (int i = 0; i < traj.front().size(); i++) {
- if (graph_flags[i]) {
- s_graph.push_back(a);
- s_graph_rbr.push_back(b);
- que1.resize(0);
- que2.resize(0);
- que3.resize(0);
- que1.push_back(i);
- que3.push_back(i);
- graph_flags[i] = false;
- while(que1.size() > 0) {
- que2.clear();
- for (int k = 0; k < que1.size(); k++) {
- for (int j = 0; j < traj.front().size(); j++) {
- if (graph[que1[k]][j].second > -1 && graph_flags[j]) {
- que2.push_back(j);
- graph_flags[j] = false;
- }
- }
- }
- que1 = que2;
- for (int j = 0; j < que2.size(); j++) {
- que3.push_back(que2[j]);
- }
- }
- s_graph.back() = que3;
- for (int j = 0; j < que3.size(); j++) {
- for (int k = 0; k < traj.front().size(); k++) {
- if (graph[que3[j]][k].second > -1) {
- s_graph_rbr.back().push_back(std::make_pair(que3[j], k));
- }
- }
+ for (int f = 0; f < traj.front().size(); f++) {
+ crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
}
- //std::cout << s_graph.back().size() << " ";
- }
- }
-}
-
-bool myfunction (const std::pair< int, float > i, const std::pair< int, float > j) {
- return i.second < j.second;
-}
-
-void graph_back_bone_evaluation(std::vector< std::vector < std::pair< int, int > > > &rout_n, int indxSize,
- std::vector< std::vector< std::pair< float, int > > > graph, std::vector< std::vector< int > > s_graph, std::vector< std::vector< std::pair< int, int > > > s_graph_rbr) {
- std::vector< float > key;
- std::vector< int > path;
- std::vector< std::pair< int, float > > que;
- std::vector< std::pair< int, int > > a;
- int v;
- for (int i = 0; i < s_graph.size(); i++) {
- key.resize(0);
- path.resize(0);
- que.resize(0);
- v = 0;
- if (s_graph[i].size() > 2) {
- key.resize(indxSize, 2);
- path.resize(indxSize, -1);
- key[s_graph[i][0]] = 0;
- for (int j = 0; j < s_graph[i].size(); j++) {
- que.push_back(std::make_pair(s_graph[i][j], key[s_graph[i][j]]));
- }
- std::sort(que.begin(), que.end(), myfunction);
- while (!que.empty()) {
- v = que[0].first;
- que.erase(que.begin());
- for (int j = 0; j < s_graph_rbr[i].size(); j++) {
- int u = -1;
- if (s_graph_rbr[i][j].first == v) {
- u = s_graph_rbr[i][j].second;
- } else if (s_graph_rbr[i][j].second == v) {
- u = s_graph_rbr[i][j].first;
- }
- bool flag = false;
- int pos = -1;
- for (int k = 0; k < que.size(); k++) {
- if (que[k].first == u) {
- flag = true;
- pos = k;
- k = que.size() + 1;
- }
- }
- if (flag && key[u] > 1 - std::abs(graph[v][u].first)) {
- path[u] = v;
- key[u] = 1 - std::abs(graph[v][u].first);
- que[pos].second = key[u];
- sort(que.begin(), que.end(), myfunction);
- }
- }
- }
- a.resize(0);
- rout_n.push_back(a);
- for (int j = 0; j < indxSize; j++) {
- if (path[j] != -1) {
- rout_n.back().push_back(std::make_pair(j, path[j]));
- }
- }
- }
- }
-}
-
-gmx::RVec evaluate_com(std::vector< RVec > frame) {
- RVec temp;
- temp[0] = 0;
- temp[1] = 0;
- temp[2] = 0;
- for (int i = 0; i < frame.size(); i++) {
- temp += frame[i];
- }
- temp[0] /= frame.size();
- temp[1] /= frame.size();
- temp[2] /= frame.size();
- return temp;
-}
-
-void evaluate_diffusion(std::vector< std::vector< RVec > > trj, std::vector< float > &D/*, int max_frame_depth*/) {
- D.resize(trj.size() * 0.9 - 1/*max_frame_depth*/);
- float temp;
- for (int i = 1; i < trj.size() * 0.9 /*max_frame_depth*/; i++) {
- temp = 0;
- for (int j = 0; j < trj.size() - 1 - i /*trj.size() - 1 - max_frame_depth*/; j++) {
- temp += (evaluate_com(trj[j]) - evaluate_com(trj[j + i])).norm2();
- //D[i][j] = (evaluate_com(trj[j]) - evaluate_com(trj[j + i])).norm2() / (2 * i);
}
- D[i - 1] = temp / (trj.size() - 1 - i) / (2 * i * 0.000001);
+ //#pragma omp barrier
+ medx.resize(0);
+ medy.resize(0);
+ std::cout << i << " corr done\n";
}
+ #pragma omp barrier
+ //}
}
/*! \brief
std::vector< int > index;
int frames = 0;
- int basic_frame = 0;
- int tau = 0; // selectable
- float crl_border = 0; // selectable
- float eff_rad = 1.5; // selectable
+ int W_size = 0;
std::string OutPutName; // selectable
// Copy and assign disallowed by base.
};
//options->addOption(FileNameOption("on").filetype(eftIndex).outputFile()
// .store(&fnNdx_).defaultBasename("domains")
// .description("Index file from the domains"));
- // Add option for tau constant
- options->addOption(gmx::IntegerOption("tau")
- .store(&tau)
+ // Add option for time window size constant
+ options->addOption(gmx::IntegerOption("Window")
+ .store(&W_size)
.description("number of frames for time travel"));
- // Add option for crl_border constant
- options->addOption(FloatOption("crl")
- .store(&crl_border)
- .description("make graph based on corrs > constant"));
- // Add option for eff_rad constant
- options->addOption(FloatOption("ef_rad")
- .store(&eff_rad)
- .description("effective radius for atoms to evaluate corrs"));
// Add option for selection list
options->addOption(SelectionOption("select_domains_and_residue").storeVector(&sel_)
.required().dynamicMask().multiValue()
{
}
-// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test6.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionList5' -tau 5 -crl 0.10
-// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/CorrsTestDomainsNfit.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionListDomainsNFit' -tau 5000 -crl 0.75 -ef_rad 1
-// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/TestCa.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelListCa' -tau 100 -crl 0.75 -ef_rad 1
-// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.tpr' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/testCa.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelListCa' -tau 9000 -crl 0.60 -ef_rad 1
+// -s '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -f '/home/toluk/Develop/samples/reca_rd/reca_rd.mono.xtc' -n '/home/toluk/Develop/samples/reca_rd/test6.ndx' -sf '/home/toluk/Develop/samples/reca_rd/SelectionList5' -Window 1000
-// -s '/home/toluk/Develop/samples/CubeTermal/CUBETermalTest1.pdb' -f '/home/toluk/Develop/samples/CubeTermal/CUBETermalTest.pdb' -n '/home/toluk/Develop/samples/CubeTermal/testCa.ndx' -sf '/home/toluk/Develop/samples/CubeTermal/SelListCa' -tau 900 -crl 0.20 -ef_rad 1
-// cube.000.000.10k.10.3.1stfrm
-// -s '/home/toluk/Develop/samples/JustCube/cube.000.000.10k.10.3.1stfrm.pdb' -f '/home/toluk/Develop/samples/JustCube/cube.000.000.10k.10.3.pdb' -n '/home/toluk/Develop/samples/JustCube/system.ndx' -sf '/home/toluk/Develop/samples/JustCube/SLsystem' -tau 1000 -crl 0.20 -ef_rad 9000 -out_put OLA
void
SpaceTimeCorr::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc,
TrajectoryAnalysisModuleData *pdata)
SpaceTimeCorr::finishAnalysis(int nframes)
{
std::vector< std::vector< std::vector< float > > > crltns;
- std::vector< std::vector< std::pair< float, int > > > graph;
- std::vector< std::vector< int > > sub_graph;
- std::vector< std::vector< std::pair< int, int > > > sub_graph_rbr, rout_new;
- int k = 1000, m = 0;
- if (tau > 0) {
- k = tau;
- }
std::cout << "\nCorrelation's evaluation - start\n";
+ correlation_evaluation(trajectory, W_size, crltns);
+ std::cout << "\nCorrelation's evaluation - end\n";
- correlation_evaluation(trajectory, basic_frame, crltns, m, k);
-
+ std::cout << "corelation matrix - start printing\n";
make_correlation_matrix_file(crltns, (OutPutName + "_matrix.txt").c_str(), 0);
- std::cout << "corelation matrix printed\n";
- make_correlation_pairs_file(crltns, (OutPutName + "_pairs.txt").c_str(), 0);
- std::cout << "corelation pairs printed\n";
-
- std::cout << "Correlation's evaluation - end\n" << "graph evaluation: start\n";
+ std::cout << "corelation matrix - printed\n";
+ std::cout << "corelation pairs - start printing\n";
+ make_correlation_pairs_file(crltns, (OutPutName + "_pairs.txt").c_str(), 0);
+ std::cout << "corelation pairs - printed\n";
-
-
-
- graph_calculation(graph, sub_graph, sub_graph_rbr, trajectory, basic_frame, crltns, crl_border, eff_rad, k);
-
- std::cout << "graph evaluation: end\n" << "routs evaluation: start\n";
-
- graph_back_bone_evaluation(rout_new, index.size(), graph, sub_graph, sub_graph_rbr);
-
- std::cout << "routs evaluation: end\n";
-
- make_rout_file(crl_border, index, rout_new, (OutPutName + "_routs.txt").c_str());
- std::cout << "corelation routs printed\n";
- make_best_corrs_graphics(crltns, rout_new, index, (OutPutName + "_routs_graphics.txt").c_str());
- std::cout << "corelation routs' pairs' graphics printed\n";
-
-
- std::cout << "extra params\n";
- std::vector< float > diffusion;
- evaluate_diffusion(trajectory, diffusion);
- make_diffusion_file((OutPutName + "_diffusion.txt").c_str(), diffusion);
-
- std::cout << "Finish Analysis - end\n\n";
}
void
SpaceTimeCorr::writeOutput()
{
-
+ std::cout << "GAME OVER\n";
}
/*! \brief
biod.pnpi.spb.ru / alexxy / gromacs-testing.git / commitdiff