std::fprintf(file, "correlations >= %0.2f\n\n", crl_border);
for (unsigned int i = 0; i < rout.size(); i++) {
for (unsigned 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, "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_table_routs(double crl_border, std::vector< int > indx, std::vector< std::vector< std::pair< unsigned int, unsigned int > > > rout, const char* file_name)
+{
+ FILE *file;
+ file = std::fopen(file_name, "w+");
+ std::fprintf(file, "correlations >= %0.2f\n\n", crl_border);
+ std::vector< std::tuple< int, int, std::vector< int > > > table;
+
+ for (unsigned int i1 = 0; i1 < rout.size(); i1++) {
+ for (unsigned int i2 = 0; i2 < rout[i1].size(); i2++) {
+ bool flag = true;
+ for (unsigned int i3 = 0; i3 < table.size(); i3++) {
+ if (std::get<0>(table[i3]) == indx[rout[i1][i2].first] + 1) {
+ std::get<1>(table[i3])++;
+ std::get<2>(table[i3]).push_back(indx[rout[i1][i2].second] + 1);
+ flag = false;
+ } else if (std::get<0>(table[i3]) == indx[rout[i1][i2].second] + 1) {
+ std::get<1>(table[i3])++;
+ std::get<2>(table[i3]).push_back(indx[rout[i1][i2].first] + 1);
+ flag = false;
+ }
+ }
+ if (flag) {
+ std::vector< int > a;
+ a.resize(0);
+ a.push_back(indx[rout[i1][i2].second] + 1);
+ table.push_back(std::make_tuple(indx[rout[i1][i2].first] + 1, 1, a));
+ a.resize(0);
+ a.push_back(indx[rout[i1][i2].first] + 1);
+ table.push_back(std::make_tuple(indx[rout[i1][i2].second] + 1, 1, a));
+ }
+ }
+ }
+
+ for (unsigned int i = 0; i < table.size(); i++) {
+ std::fprintf(file, "atomN %d connections %d | ", std::get<0>(table[i]), std::get<1>(table[i]));
+ for (unsigned int j = 0; j < std::get<2>(table[i]).size(); j++) {
+ std::fprintf(file, "%d ", std::get<2>(table[i])[j]);
+ }
+ std::fprintf(file, "\n");
+ }
+ std::fclose(file);
+}
+
void make_best_corrs_graphics(std::vector< std::vector< std::vector< double > > > correlations,
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > rout_pairs,
std::vector< int > indx,
crl.resize(tauE + 1);
for (unsigned int i = 0; i < crl.size(); i++) {
crl[i].resize(traj.front().size());
- for (unsigned int j = 0; j < crl[i].size(); j++) {
- crl[i][j].resize(traj.front().size(), 0);
+ for (unsigned int j = 0; j < traj.front().size(); j++) {
+ //crl[i][j].resize(traj.front().size(), 0);
+ for (unsigned int k = 0; k < traj.front().size(); k++) {
+ crl[i][j].push_back(0.);
+ }
}
}
- #pragma omp parallel for schedule(dynamic) shared(crl, traj, ref)
- for (unsigned int i = 0; i <= tauE; i++) {
- RVec temp1, temp2;
- std::vector< std::vector< double > > a, b, c;
- std::vector< double > d;
- d.resize(traj.front().size(), 0);
- a.resize(traj.front().size(), d);
- b.resize(traj.front().size(), d);
- c.resize(traj.front().size(), d);
- for (unsigned int j = 0; j < traj.size() - i - 1; j++) {
- for (unsigned int f1 = 0; f1 < traj.front().size(); f1++) {
- for (unsigned int f2 = 0; f2 < traj.front().size(); f2++) {
- temp1 = traj[j][f1] - ref[f1];
- temp2 = traj[j + i][f2] - ref[f2];
- a[f1][f2] += (static_cast<double>(temp1[0]) * static_cast<double>(temp2[0]) +
- static_cast<double>(temp1[1]) * static_cast<double>(temp2[1]) +
- static_cast<double>(temp1[2]) * static_cast<double>(temp2[2]));
- b[f1][f2] += (static_cast<double>(temp1[0]) * static_cast<double>(temp1[0]) +
- static_cast<double>(temp1[1]) * static_cast<double>(temp1[1]) +
- static_cast<double>(temp1[2]) * static_cast<double>(temp1[2]));
- c[f1][f2] += (static_cast<double>(temp2[0]) * static_cast<double>(temp2[0]) +
- static_cast<double>(temp2[1]) * static_cast<double>(temp2[1]) +
- static_cast<double>(temp2[2]) * static_cast<double>(temp2[2]));
+ int trjsize = traj.size(), trjsize2 = traj.front().size();
+
+ #pragma omp parallel for ordered schedule(dynamic) shared(crl, traj, ref, trjsize, trjsize2)
+ for (unsigned int i = 0; i <= tauE; i++) {
+
+ #pragma omp ordered
+ {
+
+
+ RVec temp1, temp2;
+ std::vector< std::vector< double > > a, b, c;
+ std::vector< double > d;
+ d.resize(0);
+ for (unsigned int j = 0; j < traj.front().size(); j++) {
+ a.push_back(d);
+ b.push_back(d);
+ c.push_back(d);
+ for (unsigned int k = 0; k < traj.front().size(); k++) {
+ a[j].push_back(0.);
+ b[j].push_back(0.);
+ c[j].push_back(0.);
+ }
+ }
+ /*d.resize(traj.front().size(), 0);
+ a.resize(traj.front().size(), d);
+ b.resize(traj.front().size(), d);
+ c.resize(traj.front().size(), d);*/
+
+ for (unsigned int j = 0; j < trjsize - i - 1; j++) {
+ for (unsigned int f1 = 0; f1 < trjsize2; f1++) {
+ for (unsigned int f2 = 0; f2 < trjsize2; f2++) {
+#pragma omp critical
+{
+ RVec temp3 = traj[j][f1];
+ RVec temp4 = ref[f1];
+ temp1 = temp3 - temp4;
+
+ temp3 = traj[j + i][f2];
+ temp4 = ref[f2];
+
+ temp2 = temp3 - temp4;
+
+ //temp1 = traj[j][f1] - ref[f1];
+ //temp2 = traj[j + i][f2] - ref[f2];
+
+ a[f1][f2] += (static_cast<double>(temp1[0]) * static_cast<double>(temp2[0]) +
+ static_cast<double>(temp1[1]) * static_cast<double>(temp2[1]) +
+ static_cast<double>(temp1[2]) * static_cast<double>(temp2[2]));
+ b[f1][f2] += (static_cast<double>(temp1[0]) * static_cast<double>(temp1[0]) +
+ static_cast<double>(temp1[1]) * static_cast<double>(temp1[1]) +
+ static_cast<double>(temp1[2]) * static_cast<double>(temp1[2]));
+ c[f1][f2] += (static_cast<double>(temp2[0]) * static_cast<double>(temp2[0]) +
+ static_cast<double>(temp2[1]) * static_cast<double>(temp2[1]) +
+ static_cast<double>(temp2[2]) * static_cast<double>(temp2[2]));
+}
+ }
}
}
- }
- #pragma omp critical
- for (unsigned int j = 0; j < traj.front().size(); j++) {
- for (unsigned int f = 0; f < traj.front().size(); f++) {
- crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f]));
+ #pragma omp critical
+ {
+ for (unsigned int j = 0; j < traj.front().size(); j++) {
+ for (unsigned 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 << i << " | ";
+ if (i % 100 == 0) {
+ std::cout << "\n";
+ }
+
+ }
+
}
- std::cout << i << " | ";
- if (i % 100 == 0) {
- std::cout << "\n";
- }
- }
#pragma omp barrier
std::cout << "\n";
return crl;
make_best_corrs_graphics(crltns, rout_new, index, (OutPutName + "_routs_graphics.txt").c_str());
std::cout << "corelation routs' pairs' graphics printed\n";
+ make_table_routs(static_cast< double >(crl_border), index, rout_new, (OutPutName + "_routs_table.txt").c_str());
+
/*std::cout << "extra params\n";
std::vector< double > diffusion;
evaluate_diffusion(trajectory, diffusion);
biod.pnpi.spb.ru / alexxy / gromacs-spacetimecorr.git / blobdiff