#include "corrtype.h"
-corrtype::corrtype()
-{
+// деструктор класса
+correlationType::~correlationType() {}
+// конструктор класса для инициализации
+correlationType::correlationType() {}
+
+correlationType::correlationType(std::vector< RVec > ref, int wnd, int taau, int tau_st, float crlUp, float effRad, int mod, std::string out, std::vector< std::vector < std::vector < unsigned int > > > sels, std::vector< std::string > rsNames) {
+ setDefaults(ref, wnd, taau, tau_st, crlUp, effRad, mod, out, sels, rsNames);
+}
+
+void correlationType::setDefaults(std::vector< RVec > ref, int wnd, int taau, int tau_st, float crlUp, float effRad, int mod, std::string out, std::vector< std::vector < std::vector < unsigned int > > > sels, std::vector< std::string > rsNames) {
+ if (ref.size() != 0) {reference = ref;}
+ if (wnd != -1) {window = wnd;}
+ if (taau != -1) {tau = taau;}
+ if (tau_st != -1) {tauStep = tau_st;}
+ if (crlUp != -1) {crlUpBorder = crlUp;}
+ if (effRad != -1) {effRadius = effRad;}
+ if (mod != -1) {mode = mod;}
+ if (out != "") {outputName = out;}
+ if (sels.size() != 0) {selections = sels;}
+ if (rsNames.size() > 0) {resNames = rsNames;}
+ subGraphRouts.resize(0);
+ trajectoryPartition();
+}
+
+void correlationType::update(int frame, std::vector< RVec > curFrame) {
+ trajectory.push_back(curFrame);
+ if (trajectory.size() == static_cast< unsigned int>(window + tau + 1)) {
+ trajectory.erase(trajectory.begin());
+ }
+ if ((frame + 1 - static_cast< int >(trajectory.size())) % tauStep == 0) {
+ correlationEval();
+ selections.erase(selections.begin());
+ subGraphRouts.resize(subGraphRouts.size() + 1);
+ graphCalculations(1, static_cast< unsigned int >(tau));
+ graphBackBoneEvaluation();
+ }
+}
+
+void correlationType::printData() {
+ printOutputData();
+}
+
+void correlationType::trajectoryPartition() {
+ std::vector< bool > temp1;
+ std::pair< unsigned int, unsigned int > temp2;
+ float temp3;
+ for (unsigned k = 0; k < selections.size(); k++) {
+ temp1.resize(0); temp1.resize(index.size(), true);
+ for (unsigned int i = 0; i < selections[k].size(); i++) {
+ for (unsigned int j = 0; j < selections[k][i].size(); j++) {
+ temp1[selections[k][i][j]] = false;
+ }
+ }
+ temp3 = 9000;
+ for (unsigned i = 0; i < temp1.size(); i++) {
+ if (temp1[i]) {
+ for (unsigned int i = 0; i < selections[k].size(); i++) {
+ for (unsigned int j = 0; j < selections[k][i].size(); j++) {
+ if (temp3 > (reference[selections[k][i][j]] - reference[i]).norm()) {
+ temp3 = (reference[selections[k][i][j]] - reference[i]).norm();
+ temp2 = std::make_pair(i, j);
+ }
+ }
+ }
+ selections[k][temp2.first].push_back(i);
+ }
+ }
+ }
+}
+
+void correlationType::readWriteCorrelations(int rwMode) {
+ FILE *file;
+ if (rwMode == 1) {
+ file = std::fopen((outputName + "-matrixData").c_str(), "a");
+ for (unsigned int i = 0; i < matrixes.size(); i++) {
+ std::fprintf(file, "%d %d\n", count, i);
+ for (unsigned int j = 0; j < matrixes[i].size(); j++) {
+ for (unsigned int f = 0; f < matrixes[i][j].size(); f++) {
+ std::fprintf(file, "%.4f ", matrixes[i][j][f]); //~16
+ }
+ std::fprintf(file, "\n");
+ }
+ }
+ std::fclose(file);
+ }
+ if (rwMode == 0) {
+ file = std::fopen((outputName + "-matrixData").c_str(), "r+");
+ matrixes.resize(0); matrixes.resize(static_cast< unsigned int >(tau + 1));
+ for (unsigned int i = 0; i < static_cast< unsigned int >(tau + 1); i++) {
+ int t0, t1, t2 = std::fscanf(file, "%d %d\n", &t0, &t1);
+ matrixes[i].resize(0); matrixes[i].resize(index.size());
+ for (unsigned int j = 0; j < index.size(); j++) {
+ matrixes[i][j].resize(index.size());
+ for (unsigned int k = 0; k < index.size(); k++) {
+ t2 = std::fscanf(file, "%lf ", &matrixes[i][j][k]);
+ }
+ }
+ }
+ }
+}
+
+void correlationType::correlationEval() {
+ /*
+ * fitting block / we add spare atoms to existing domains based on proximity
+ */
+ std::vector< std::vector< std::pair< unsigned int, unsigned int > > > pairs;
+ pairs.resize(0); pairs.resize(selections.front().size());
+ for (unsigned int i = 0; i < selections.front().size(); i++) {
+ pairs[i].resize(0);
+ for (unsigned int j = 0; j < selections.front()[i].size(); j++) {
+ pairs[i].push_back(std::make_pair(selections.front()[i][j], selections.front()[i][j]));
+ }
+ }
+ fitTrajectory.resize(0); fitTrajectory.resize(selections.front().size(), trajectory);
+ #pragma omp parallel for schedule(dynamic) firstprivate(reference)
+ for (unsigned long i = 0; i < selections.front().size(); i++) {
+ for (unsigned int j = 0; j < fitTrajectory[i].size(); j++) {
+ MyFitNew(reference, fitTrajectory[i][j], pairs[i], 0);
+ }
+ }
+ #pragma omp barrier
+ frankensteinTrajectory = trajectory;
+ for (unsigned int i = 0; i < selections.front().size(); i++) {
+ for (unsigned int j = 0; j < selections.front()[i].size(); j++) {
+ for (unsigned int k = 0; k < fitTrajectory[i].size(); k++) {
+ frankensteinTrajectory[k][selections.front()[i][j]] = fitTrajectory[i][k][selections.front()[i][j]];
+ }
+ }
+ }
+ /*
+ * fitting block
+ */
+ matrixes.resize(static_cast< unsigned int >(tau + 1));
+ for (unsigned int i = 0; i < matrixes.size(); i++) {
+ matrixes[i].resize(index.size());
+ for (unsigned int j = 0; j < index.size(); j++) {
+ matrixes[i][j].resize(index.size());
+ for (unsigned int k = 0; k < index.size(); k++) {
+ matrixes[i][j][k] = 0.;
+ }
+ }
+ }
+ std::vector< std::vector< double > > a, b, c;
+ std::vector< double > d;
+ #pragma omp parallel for ordered schedule(dynamic) shared(matrixes) firstprivate(frankensteinTrajectory, reference)
+ for (unsigned int i = 0; i <= static_cast< unsigned int >(tau); i++) {
+ a.resize(0); b.resize(0); c.resize(0); d.resize(0);
+ for (unsigned int j = 0; j < index.size(); j++) {
+ a.push_back(d); b.push_back(d); c.push_back(d);
+ for (unsigned int k = 0; k < index.size(); k++) {
+ a[j].push_back(0.); b[j].push_back(0.); c[j].push_back(0.);
+ }
+ }
+ for (unsigned int j = 0; j < static_cast< unsigned int >(window); j++) {
+ for (unsigned int k1 = 0; k1 < index.size(); k1++) {
+ for (unsigned int k2 = 0; k2 < index.size(); k2++) {
+ RVec temp1, temp2;
+ temp1 = frankensteinTrajectory[j][k1] - reference[k1];
+ temp2 = frankensteinTrajectory[j + i][k2] - reference[k2];
+ a[k1][k2] += (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[k1][k2] += (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[k1][k2] += (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]));
+ }
+ }
+ }
+ for (unsigned int j = 0; j < index.size(); j++) {
+ for (unsigned int k = 0; k < index.size(); k++) {
+ matrixes[i][j][k] = a[j][k] / (std::sqrt(b[j][k] * c[j][k]));
+ }
+ }
+ }
+ #pragma omp barrier
+ for (unsigned int i = 0; i < matrixes.size(); i++) {
+ for (unsigned int j = 0; j < matrixes[i].size(); j++) {
+ for (unsigned int k = 0; k < matrixes[i][j].size(); k++) {
+ matrixes[i][j][k] = std::round(matrixes[i][j][k] * 10000) / 10000;
+ }
+ }
+ }
+ readWriteCorrelations(1);
+ count++;
+}
+
+void correlationType::graphCalculations(unsigned int tauStart, unsigned int tauEnd) {
+ graph.resize(0); graph.resize(index.size());
+ for (unsigned int i = 0; i < index.size(); i++) {
+ graph[i].resize(index.size(), std::make_pair(0, -1));
+ }
+ RVec temp;
+ for (unsigned int i = tauStart; i <= tauEnd; i++) {
+ for (unsigned int j = 0; j < index.size(); j++) {
+ for (unsigned int k = j; k < index.size(); k++) {
+ temp = reference[i] - reference[k];
+ if (std::max(std::abs(matrixes[i][j][k]), std::abs(matrixes[i][k][j])) >= static_cast< double >(crlUpBorder) &&
+ static_cast< float >(norm(temp)) <= effRadius && std::abs(graph[j][k].first) < std::max(std::abs(matrixes[i][j][k]), std::abs(matrixes[i][k][j]))) {
+ if (std::abs(matrixes[i][j][k]) > std::abs(matrixes[i][k][j])) {
+ graph[j][k].first = matrixes[i][j][k];
+ } else {
+ graph[j][k].first = matrixes[i][k][j];
+ }
+ graph[j][k].second = static_cast< int >(i);
+ }
+ }
+ }
+ }
+ std::vector< bool > graph_flags;
+ graph_flags.resize(0); graph_flags.resize(index.size(), true);
+ std::vector< unsigned int > a;
+ std::vector< std::pair< unsigned int, unsigned int > > b;
+ a.resize(0); b.resize(0);
+ std::vector< unsigned int > width1, width2, tempSubGraph;
+ for (unsigned int i = 0; i < index.size(); i++) {
+ if (graph_flags[i]) {
+ subGraphPoints.push_back(a);
+ subGraphRbr.push_back(b);
+ width1.resize(0); width2.resize(0); tempSubGraph.resize(0);
+ width1.push_back(i);
+ tempSubGraph.push_back(i);
+ graph_flags[i] = false;
+ while(width1.size() > 0) {
+ width2.clear();
+ for (unsigned int j = 0; j < width1.size(); j++) {
+ for (unsigned int k = 0; k < index.size(); k++) {
+ if (graph[width1[j]][k].second > -1 && graph_flags[k]) {
+ width2.push_back(k);
+ graph_flags[k] = false;
+ }
+ }
+ }
+ width1 = width2;
+ for (unsigned int j = 0; j < width2.size(); j++) {
+ tempSubGraph.push_back(width2[j]);
+ }
+ }
+ subGraphPoints.back() = tempSubGraph;
+ for (unsigned int j = 0; j < tempSubGraph.size(); j++) {
+ for (unsigned int k = 0; k < index.size(); k++) {
+ if (graph[tempSubGraph[j]][k].second > -1) {
+ subGraphRbr.back().push_back(std::make_pair(tempSubGraph[j], k));
+ }
+ }
+ }
+ }
+ }
+}
+
+bool correlationType::myComparisonFunction (const std::pair< int, double > i, const std::pair< int, double > j) {
+ return i.second < j.second;
+}
+
+void correlationType::graphBackBoneEvaluation() {
+ std::vector< double > key;
+ std::vector< long > path;
+ std::vector< std::pair< unsigned int, double > > que;
+ std::vector< std::pair< unsigned int, unsigned int > > a;
+ unsigned int v;
+ a.resize(0);
+ subGraphRouts.resize(0);
+ for (unsigned int i = 0; i < subGraphPoints.size(); i++) {
+ key.resize(0);
+ path.resize(0);
+ que.resize(0);
+ v = 0;
+ if (subGraphPoints[i].size() > 2) {
+ key.resize(index.size(), 2);
+ path.resize(index.size(), -1);
+ key[subGraphPoints[i][0]] = 0;
+ for (unsigned int j = 0; j < subGraphPoints[i].size(); j++) {
+ que.push_back(std::make_pair(subGraphPoints[i][j], key[subGraphPoints[i][j]]));
+ }
+ std::sort(que.begin(), que.end(), myComparisonFunction);
+ while (!que.empty()) {
+ v = que[0].first;
+ que.erase(que.begin());
+ for (unsigned int j = 0; j < subGraphRbr[i].size(); j++) {
+ long u = -1;
+ if (subGraphRbr[i][j].first == v) {
+ u = subGraphRbr[i][j].second;
+ } else if (subGraphRbr[i][j].second == v) {
+ u = subGraphRbr[i][j].first;
+ }
+ bool flag = false;
+ unsigned long pos = 0;
+ for (unsigned long k = 0; k < que.size(); k++) {
+ if (que[k].first == u) {
+ flag = true;
+ pos = k;
+ k = que.size() + 1;
+ }
+ }
+ if (flag && key[static_cast< unsigned long >(u)] > 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first)) {
+ path[static_cast< unsigned long >(u)] = v;
+ key[static_cast< unsigned long >(u)] = 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first);
+ que[pos].second = key[static_cast< unsigned long >(u)];
+ sort(que.begin(), que.end(), myComparisonFunction);
+ }
+ }
+ }
+ subGraphRouts.back().push_back(a);
+ for (unsigned int j = 0; j < index.size(); j++) {
+ if (path[j] != -1) {
+ subGraphRouts.back().back().push_back(std::make_pair(j, path[j]));
+ }
+ }
+ }
+ }
+}
+
+void correlationType::printOutputData() {
+ FILE *file;
+ file = std::fopen((outputName + "-arrowsData.txt").c_str(), "w+");
+ unsigned int same, pre = 0;
+ std::vector< std::tuple< int, int, std::vector< int > > > table;
+ std::vector< int > a;
+ for (unsigned int i = 0; i < subGraphRouts.size(); i++) {
+ same = i;
+ for (unsigned int j = i + 1; j < subGraphRouts.size(); j++) {
+ if (subGraphRouts[j] == subGraphRouts[i]) {
+ same = j;
+ } else {
+ break;
+ }
+ }
+ if (i == same) {
+ std::fprintf(file, "\n Starting time point = %d | correlations >= %0.2f | tau = %d | window = %d\n\n", static_cast< int >(i) * tauStep, static_cast< double >(crlUpBorder), tau, window);
+ } else {
+ std::fprintf(file, "\n Starting time point = [%d ; %d] | correlations >= %0.2f | tau = %d | window = %d\n\n", static_cast< int >(i) * tauStep, static_cast< int >(same) * tauStep, static_cast< double >(crlUpBorder), tau, window);
+ }
+ for (unsigned int j = 0; j < subGraphRouts[i].size(); j++) {
+ for (unsigned int k = 0; k < subGraphRouts[i][j].size(); k++) {
+ std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.05\n", index[subGraphRouts[i][j][k].first] + 1, index[subGraphRouts[i][j][k].second] + 1);
+ }
+ std::fprintf(file, "\n");
+ }
+ table.resize(0);
+ for (unsigned int j = 0; j < subGraphRouts[i].size(); j++) {
+ for (unsigned int k = 0; k < subGraphRouts[i][j].size(); k++) {
+ bool flag1 = true, flag2 = true;
+ for (unsigned int m = 0; m < table.size(); m++) {
+ if (std::get<0>(table[m]) == index[subGraphRouts[i][j][k].first]) {
+ std::get<1>(table[m])++;
+ std::get<2>(table[m]).push_back(index[subGraphRouts[i][j][k].second]);
+ flag1 = false;
+ }
+ if (std::get<0>(table[m]) == index[subGraphRouts[i][j][k].second]) {
+ std::get<1>(table[m])++;
+ std::get<2>(table[m]).push_back(index[subGraphRouts[i][j][k].first]);
+ flag2 = false;
+ }
+ }
+ if (flag1) {
+ a.resize(0);
+ a.push_back(index[subGraphRouts[i][j][k].second]);
+ table.push_back(std::make_tuple(index[subGraphRouts[i][j][k].first], 1, a));
+ }
+ if (flag2) {
+ a.resize(0);
+ a.push_back(index[subGraphRouts[i][j][k].first]);
+ table.push_back(std::make_tuple(index[subGraphRouts[i][j][k].second], 1, a));
+ }
+ }
+ }
+ for (unsigned int j = 0; j < table.size(); j++) {
+ std::fprintf(file, "residue %s connections %d | ", (resNames[static_cast< unsigned int >(std::find(index.begin(), index.end(), std::get<0>(table[j])) - index.begin())]).c_str(), std::get<1>(table[j]));
+ for (unsigned int k = 0; k < std::get<2>(table[j]).size(); k++) {
+ std::fprintf(file, "%s ", (resNames[static_cast< unsigned int >(std::find(index.begin(), index.end(), std::get<2>(table[j])[k]) - index.begin())]).c_str());
+ }
+ std::fprintf(file, "\n");
+ }
+ if (pre != 0) {
+ std::vector< std::vector< std::pair< unsigned int, unsigned int > > > temp01, temp02;
+ std::vector< std::pair< unsigned int, unsigned int > > temp03, temp04, temp05;
+ temp01 = subGraphRouts[pre];
+ temp02 = subGraphRouts[same];
+ temp03.resize(0);
+ temp04.resize(0);
+ for (unsigned int j = 0; j < temp01.size(); j++) {
+ for (unsigned int k = 0; k < temp01[j].size(); k++) {
+ temp03.push_back(temp01[j][k]);
+ }
+ }
+ for (unsigned int j = 0; j < temp02.size(); j++) {
+ for (unsigned int k = 0; k < temp02[j].size(); k++) {
+ temp04.push_back(temp02[j][k]);
+ }
+ }
+ std::sort(temp03.begin(), temp03.end());
+ std::sort(temp04.begin(), temp04.end());
+ temp05.resize(0);
+ std::set_difference(temp03.begin(), temp03.end(), temp04.begin(), temp04.end(), std::inserter(temp05, temp05.begin()));
+ std::fprintf(file, "minus:\n");
+ for (unsigned int j = 0; j < temp05.size(); j++) {
+ std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.05\n", index[temp05[j].first] + 1, index[temp05[j].second] + 1);
+ }
+ temp05.resize(0);
+ std::set_difference(temp04.begin(), temp04.end(), temp03.begin(), temp03.end(), std::inserter(temp05, temp05.begin()));
+ std::fprintf(file, "plus:\n");
+ for (unsigned int j = 0; j < temp05.size(); j++) {
+ std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.05\n", index[temp05[j].first] + 1, index[temp05[j].second] + 1);
+ }
+ }
+ pre = same;
+ i = same;
+ }
+ std::fclose(file);
}
* \ingroup module_trajectoryanalysis
*/
-
#include <iostream>
#include <vector>
#include <cfloat>
-
-
#include <gromacs/trajectoryanalysis.h>
// #include <gromacs/trajectoryanalysis/topologyinformation.h>
#include "/home/titov_ai/Develop/gromacs-original/src/gromacs/trajectoryanalysis/topologyinformation.h"
+#include "newfit.h"
+#include "corrtype.h"
+
using namespace gmx;
using gmx::RVec;
-// вычисление модуля расстояния между двумя точками, с учётом геометрического преобразования
-double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
- return sqrt( pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) +
- pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) +
- pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2) );
-}
-
-// вычисление функции F и её частных производных
-void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
- double t01, t02, t03, t04, t05, t06, t07;
- t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
- t02 = pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2);
- t03 = pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2);
- t04 = sqrt(t01 + t02 + t03);
- t05 = (aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x);
- t06 = (aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y);
- t07 = (aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x);
- F += t04;
- Fx += -(2 * cos(B) * cos(C) * t05 - 2 * sin(B) * t06 + 2 * cos(B) * sin(C) * t07) / (2 * t04);
- Fy += -(2 * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) * t07 - 2 * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) * t05 + 2 * cos(B) * sin(A) * t06) / (2 * t04);
- Fz += -(2 * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) * t05 - 2 * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) * t07 + 2 * cos(A) * cos(B) * t06) / (2 * t04);
- Fa += -(2 * (cos(A) * cos(B) * biy_plus_y - cos(B) * sin(A) * biz_plus_z) * t06 -
- 2 * (biy_plus_y * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + biz_plus_z * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C))) * t07 +
- 2 * (biy_plus_y * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) + biz_plus_z * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B))) * t05) / (2 * t04);
- Fb += -(2 * (cos(A) * cos(B) * sin(C) * biz_plus_z - sin(B) * sin(C) * bix_plus_x + cos(B) * sin(A) * sin(C) * biy_plus_y) * t07 +
- 2 * (cos(A) * cos(B) * cos(C) * biz_plus_z - cos(C) * sin(B) * bix_plus_x + cos(B) * cos(C) * sin(A) * biy_plus_y) * t05 -
- 2 * (cos(B) * bix_plus_x + sin(A) * sin(B) * biy_plus_y + cos(A) * sin(B) * biz_plus_z) * t06) / (2 * t04);
- Fc += (2 * (biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * bix_plus_x) * t05 -
- 2 * (biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * bix_plus_x) * t07) / (2 * t04);
-}
-
-// применения геометрического преобразования: смещение на вектор (x, y, z) и повороты на градусы A, B, C относительно осей (x, y, z)
-void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
- double t0 = 0, t1 = 0, t2 = 0;
- for (unsigned int i = 0; i < b.size(); i++) {
- t0 = static_cast< double >(b[i][0]);
- t1 = static_cast< double >(b[i][1]);
- t2 = static_cast< double >(b[i][2]);
- b[i][0] = static_cast< float >((t2 + z) * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - (t1 + y) * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * (t0 + x));
- b[i][1] = static_cast< float >((t1 + y) * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) - (t2 + z) * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) + cos(B) * sin(C) * (t0 + x));
- b[i][2] = static_cast< float >(cos(A) * cos(B) * (t2 + z) - sin(B) * (t0 + x) + cos(B) * sin(A) * (t1 + y));
- }
-}
-
-// подсчёт геометрических центров множеств a и b на основе пар pairs
-void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< unsigned int, unsigned int > > pairs) {
- midA[0] = 0;
- midA[1] = 0;
- midA[2] = 0;
-
- midB[0] = 0;
- midB[1] = 0;
- midB[2] = 0;
-
- for (unsigned int i = 0; i < pairs.size(); i++) {
- midA += a[pairs[i].first];
- midB += b[pairs[i].second];
- }
- midA[0] /= pairs.size();
- midA[1] /= pairs.size();
- midA[2] /= pairs.size();
-
- midB[0] /= pairs.size();
- midB[1] /= pairs.size();
- midB[2] /= pairs.size();
-}
-
-// функция фитирования фрейма b на фрейм a на основе пар pairs с "точностью" FtCnst
-void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, double FtCnst) {
- double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
- RVec ma, mb;
- CalcMid(a, b, ma, mb, pairs);
- ma -= mb;
- double x = static_cast< double >(ma[0]), y = static_cast< double >(ma[1]), z = static_cast< double >(ma[2]), A = 0, B = 0, C = 0;
- // поиск стартового отклонения множеств
- for (unsigned int i = 0; i < pairs.size(); i++) {
- f1 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
- static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y, static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
- }
- if (FtCnst == 0) {
- FtCnst = 0.00001;
- }
- if (f1 == 0) {
- return;
- } else {
- // поиск оптимального геометрического преобразования методом градиентного спуска
- while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
- f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0; l = 1;
- // поиск производных и отклонения при заданных параметрах
- for (unsigned int i = 0; i < pairs.size(); i++) {
- searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc,
- static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
- static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y,
- static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
- }
- while (f2 != f1) {
- f2 = 0;
- // поиск отклонения при новых параметрах
- for (unsigned int i = 0; i < pairs.size(); i++) {
- f2 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
- static_cast< double >(b[pairs[i].second][0]) + (x - l * fx), static_cast< double >(b[pairs[i].second][1]) + (y - l * fy),
- static_cast< double >(b[pairs[i].second][2]) + (z - l * fz), A - l * fa, B - l * fb, C - l * fc);
- }
- if (f2 < f1) {
- x -= l * fx; y -= l * fy; z -= l * fz; A -= l * fa; B -= l * fb; C -= l * fc;
- fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
- break;
- } else {
- // по факту существует минимальное значение переменной типа double
- // в случае его достижения дважды останавливаем цикл т.к. упёрлись в предел допустимой точности
- // таким образом гарантируем выход из цикла
- if (l == DBL_MIN) { //DBL_TRUE_MIN
- break;
- }
- l *= 0.50;
- }
- }
- }
- ApplyFit(b, x, y, z, A, B, C);
- }
-}
-
-class correlationType {
-
- public:
-
- // деструктор класса
- ~correlationType() {}
-
- // конструктор класса для инициализации
- correlationType() {}
-
- correlationType(std::vector< RVec > ref, int wnd, int taau, int tau_st, float crlUp, float effRad, int mod, std::string out, std::vector< std::vector < std::vector < unsigned int > > > sels, std::vector< std::string > rsNames) {
- setDefaults(ref, wnd, taau, tau_st, crlUp, effRad, mod, out, sels, rsNames);
- }
-
- void setDefaults(std::vector< RVec > ref, int wnd, int taau, int tau_st, float crlUp, float effRad, int mod, std::string out, std::vector< std::vector < std::vector < unsigned int > > > sels, std::vector< std::string > rsNames) {
- if (ref.size() != 0) {reference = ref;}
- if (wnd != -1) {window = wnd;}
- if (taau != -1) {tau = taau;}
- if (tau_st != -1) {tauStep = tau_st;}
- if (crlUp != -1) {crlUpBorder = crlUp;}
- if (effRad != -1) {effRadius = effRad;}
- if (mod != -1) {mode = mod;}
- if (out != "") {outputName = out;}
- if (sels.size() != 0) {selections = sels;}
- if (rsNames.size() > 0) {resNames = rsNames;}
- subGraphRouts.resize(0);
- trajectoryPartition();
- }
-
- void update(int frame, std::vector< RVec > curFrame) {
- trajectory.push_back(curFrame);
- if (trajectory.size() == static_cast< unsigned int>(window + tau + 1)) {
- trajectory.erase(trajectory.begin());
- }
- if ((frame + 1 - static_cast< int >(trajectory.size())) % tauStep == 0) {
- correlationEval();
- selections.erase(selections.begin());
- subGraphRouts.resize(subGraphRouts.size() + 1);
- graphCalculations(1, static_cast< unsigned int >(tau));
- graphBackBoneEvaluation();
- }
- }
-
- void printData() {
- printOutputData();
- }
-
- private:
- std::vector< std::vector< std::vector< double > > > matrixes;
- std::vector< std::vector< RVec > > trajectory;
- std::vector< std::vector< RVec > > frankensteinTrajectory;
- std::vector< std::vector< std::vector< RVec > > > fitTrajectory;
- std::vector< RVec > reference;
- int window = 1000; // selectable
- int tau = window / 2; // selectable
- int tauStep = window / 10; // selectable
- float crlUpBorder = 0.5; // selectable
- float effRadius = 1.5; // selectable
- int mode = 0; // selectable
- std::string outputName = ""; // selectable
- std::vector< int > index;
- std::vector< std::string > resNames;
- std::vector< std::vector < std::vector < unsigned int > > > selections;
- int count = 0;
-
- std::vector< std::vector< std::pair< double, int > > > graph;
- std::vector< std::vector< unsigned int > > subGraphPoints;
- std::vector< std::vector< std::pair< unsigned int, unsigned int > > > subGraphRbr;
- std::vector< std::vector< std::vector< std::pair< unsigned int, unsigned int > > > > subGraphRouts;
-
- void trajectoryPartition() {
- std::vector< bool > temp1;
- std::pair< unsigned int, unsigned int > temp2;
- float temp3;
- for (unsigned k = 0; k < selections.size(); k++) {
- temp1.resize(0); temp1.resize(index.size(), true);
- for (unsigned int i = 0; i < selections[k].size(); i++) {
- for (unsigned int j = 0; j < selections[k][i].size(); j++) {
- temp1[selections[k][i][j]] = false;
- }
- }
- temp3 = 9000;
- for (unsigned i = 0; i < temp1.size(); i++) {
- if (temp1[i]) {
- for (unsigned int i = 0; i < selections[k].size(); i++) {
- for (unsigned int j = 0; j < selections[k][i].size(); j++) {
- if (temp3 > (reference[selections[k][i][j]] - reference[i]).norm()) {
- temp3 = (reference[selections[k][i][j]] - reference[i]).norm();
- temp2 = std::make_pair(i, j);
- }
- }
- }
- selections[k][temp2.first].push_back(i);
- }
- }
- }
- }
-
- void readWriteCorrelations(int rwMode) {
- FILE *file;
- if (rwMode == 1) {
- file = std::fopen((outputName + "-matrixData").c_str(), "a");
- for (unsigned int i = 0; i < matrixes.size(); i++) {
- std::fprintf(file, "%d %d\n", count, i);
- for (unsigned int j = 0; j < matrixes[i].size(); j++) {
- for (unsigned int f = 0; f < matrixes[i][j].size(); f++) {
- std::fprintf(file, "%.4f ", matrixes[i][j][f]); //~16
- }
- std::fprintf(file, "\n");
- }
- }
- std::fclose(file);
- }
- if (rwMode == 0) {
- file = std::fopen((outputName + "-matrixData").c_str(), "r+");
- matrixes.resize(0); matrixes.resize(static_cast< unsigned int >(tau + 1));
- for (unsigned int i = 0; i < static_cast< unsigned int >(tau + 1); i++) {
- int t0, t1, t2 = std::fscanf(file, "%d %d\n", &t0, &t1);
- matrixes[i].resize(0); matrixes[i].resize(index.size());
- for (unsigned int j = 0; j < index.size(); j++) {
- matrixes[i][j].resize(index.size());
- for (unsigned int k = 0; k < index.size(); k++) {
- t2 = std::fscanf(file, "%lf ", &matrixes[i][j][k]);
- }
- }
- }
- }
-
- }
-
- void correlationEval() {
- /*
- * fitting block / we add spare atoms to existing domains based on proximity
- */
- std::vector< std::vector< std::pair< unsigned int, unsigned int > > > pairs;
- pairs.resize(0); pairs.resize(selections.front().size());
- for (unsigned int i = 0; i < selections.front().size(); i++) {
- pairs[i].resize(0);
- for (unsigned int j = 0; j < selections.front()[i].size(); j++) {
- pairs[i].push_back(std::make_pair(selections.front()[i][j], selections.front()[i][j]));
- }
- }
- fitTrajectory.resize(0); fitTrajectory.resize(selections.front().size(), trajectory);
- #pragma omp parallel for schedule(dynamic) firstprivate(reference)
- for (unsigned long i = 0; i < selections.front().size(); i++) {
- for (unsigned int j = 0; j < fitTrajectory[i].size(); j++) {
- MyFitNew(reference, fitTrajectory[i][j], pairs[i], 0);
- }
- }
- #pragma omp barrier
- frankensteinTrajectory = trajectory;
- for (unsigned int i = 0; i < selections.front().size(); i++) {
- for (unsigned int j = 0; j < selections.front()[i].size(); j++) {
- for (unsigned int k = 0; k < fitTrajectory[i].size(); k++) {
- frankensteinTrajectory[k][selections.front()[i][j]] = fitTrajectory[i][k][selections.front()[i][j]];
- }
- }
- }
- /*
- * fitting block
- */
- matrixes.resize(static_cast< unsigned int >(tau + 1));
- for (unsigned int i = 0; i < matrixes.size(); i++) {
- matrixes[i].resize(index.size());
- for (unsigned int j = 0; j < index.size(); j++) {
- matrixes[i][j].resize(index.size());
- for (unsigned int k = 0; k < index.size(); k++) {
- matrixes[i][j][k] = 0.;
- }
- }
- }
- std::vector< std::vector< double > > a, b, c;
- std::vector< double > d;
- #pragma omp parallel for ordered schedule(dynamic) shared(matrixes) firstprivate(frankensteinTrajectory, reference)
- for (unsigned int i = 0; i <= static_cast< unsigned int >(tau); i++) {
- a.resize(0); b.resize(0); c.resize(0); d.resize(0);
- for (unsigned int j = 0; j < index.size(); j++) {
- a.push_back(d); b.push_back(d); c.push_back(d);
- for (unsigned int k = 0; k < index.size(); k++) {
- a[j].push_back(0.); b[j].push_back(0.); c[j].push_back(0.);
- }
- }
- for (unsigned int j = 0; j < static_cast< unsigned int >(window); j++) {
- for (unsigned int k1 = 0; k1 < index.size(); k1++) {
- for (unsigned int k2 = 0; k2 < index.size(); k2++) {
- RVec temp1, temp2;
- temp1 = frankensteinTrajectory[j][k1] - reference[k1];
- temp2 = frankensteinTrajectory[j + i][k2] - reference[k2];
- a[k1][k2] += (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[k1][k2] += (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[k1][k2] += (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]));
- }
- }
- }
- for (unsigned int j = 0; j < index.size(); j++) {
- for (unsigned int k = 0; k < index.size(); k++) {
- matrixes[i][j][k] = a[j][k] / (std::sqrt(b[j][k] * c[j][k]));
- }
- }
- }
- #pragma omp barrier
- for (unsigned int i = 0; i < matrixes.size(); i++) {
- for (unsigned int j = 0; j < matrixes[i].size(); j++) {
- for (unsigned int k = 0; k < matrixes[i][j].size(); k++) {
- matrixes[i][j][k] = std::round(matrixes[i][j][k] * 10000) / 10000;
- }
- }
- }
- readWriteCorrelations(1);
- count++;
- }
-
- void graphCalculations(unsigned int tauStart, unsigned int tauEnd) {
- graph.resize(0); graph.resize(index.size());
- for (unsigned int i = 0; i < index.size(); i++) {
- graph[i].resize(index.size(), std::make_pair(0, -1));
- }
- RVec temp;
- for (unsigned int i = tauStart; i <= tauEnd; i++) {
- for (unsigned int j = 0; j < index.size(); j++) {
- for (unsigned int k = j; k < index.size(); k++) {
- temp = reference[i] - reference[k];
- if (std::max(std::abs(matrixes[i][j][k]), std::abs(matrixes[i][k][j])) >= static_cast< double >(crlUpBorder) &&
- static_cast< float >(norm(temp)) <= effRadius && std::abs(graph[j][k].first) < std::max(std::abs(matrixes[i][j][k]), std::abs(matrixes[i][k][j]))) {
- if (std::abs(matrixes[i][j][k]) > std::abs(matrixes[i][k][j])) {
- graph[j][k].first = matrixes[i][j][k];
- } else {
- graph[j][k].first = matrixes[i][k][j];
- }
- graph[j][k].second = static_cast< int >(i);
- }
- }
- }
- }
- std::vector< bool > graph_flags;
- graph_flags.resize(0); graph_flags.resize(index.size(), true);
- std::vector< unsigned int > a;
- std::vector< std::pair< unsigned int, unsigned int > > b;
- a.resize(0); b.resize(0);
- std::vector< unsigned int > width1, width2, tempSubGraph;
- for (unsigned int i = 0; i < index.size(); i++) {
- if (graph_flags[i]) {
- subGraphPoints.push_back(a);
- subGraphRbr.push_back(b);
- width1.resize(0); width2.resize(0); tempSubGraph.resize(0);
- width1.push_back(i);
- tempSubGraph.push_back(i);
- graph_flags[i] = false;
- while(width1.size() > 0) {
- width2.clear();
- for (unsigned int j = 0; j < width1.size(); j++) {
- for (unsigned int k = 0; k < index.size(); k++) {
- if (graph[width1[j]][k].second > -1 && graph_flags[k]) {
- width2.push_back(k);
- graph_flags[k] = false;
- }
- }
- }
- width1 = width2;
- for (unsigned int j = 0; j < width2.size(); j++) {
- tempSubGraph.push_back(width2[j]);
- }
- }
- subGraphPoints.back() = tempSubGraph;
- for (unsigned int j = 0; j < tempSubGraph.size(); j++) {
- for (unsigned int k = 0; k < index.size(); k++) {
- if (graph[tempSubGraph[j]][k].second > -1) {
- subGraphRbr.back().push_back(std::make_pair(tempSubGraph[j], k));
- }
- }
- }
- }
- }
- }
-
- static bool myComparisonFunction (const std::pair< int, double > i, const std::pair< int, double > j) {
- return i.second < j.second;
- }
-
- void graphBackBoneEvaluation() {
- std::vector< double > key;
- std::vector< long > path;
- std::vector< std::pair< unsigned int, double > > que;
- std::vector< std::pair< unsigned int, unsigned int > > a;
- unsigned int v;
- a.resize(0);
- subGraphRouts.resize(0);
- for (unsigned int i = 0; i < subGraphPoints.size(); i++) {
- key.resize(0);
- path.resize(0);
- que.resize(0);
- v = 0;
- if (subGraphPoints[i].size() > 2) {
- key.resize(index.size(), 2);
- path.resize(index.size(), -1);
- key[subGraphPoints[i][0]] = 0;
- for (unsigned int j = 0; j < subGraphPoints[i].size(); j++) {
- que.push_back(std::make_pair(subGraphPoints[i][j], key[subGraphPoints[i][j]]));
- }
- std::sort(que.begin(), que.end(), myComparisonFunction);
- while (!que.empty()) {
- v = que[0].first;
- que.erase(que.begin());
- for (unsigned int j = 0; j < subGraphRbr[i].size(); j++) {
- long u = -1;
- if (subGraphRbr[i][j].first == v) {
- u = subGraphRbr[i][j].second;
- } else if (subGraphRbr[i][j].second == v) {
- u = subGraphRbr[i][j].first;
- }
- bool flag = false;
- unsigned long pos = 0;
- for (unsigned long k = 0; k < que.size(); k++) {
- if (que[k].first == u) {
- flag = true;
- pos = k;
- k = que.size() + 1;
- }
- }
- if (flag && key[static_cast< unsigned long >(u)] > 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first)) {
- path[static_cast< unsigned long >(u)] = v;
- key[static_cast< unsigned long >(u)] = 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first);
- que[pos].second = key[static_cast< unsigned long >(u)];
- sort(que.begin(), que.end(), myComparisonFunction);
- }
- }
- }
- subGraphRouts.back().push_back(a);
- for (unsigned int j = 0; j < index.size(); j++) {
- if (path[j] != -1) {
- subGraphRouts.back().back().push_back(std::make_pair(j, path[j]));
- }
- }
- }
- }
- }
-
- void printOutputData() {
- FILE *file;
- file = std::fopen((outputName + "-arrowsData.txt").c_str(), "w+");
-
- unsigned int same, pre = 0;
- std::vector< std::tuple< int, int, std::vector< int > > > table;
- std::vector< int > a;
-
- for (unsigned int i = 0; i < subGraphRouts.size(); i++) {
- same = i;
- for (unsigned int j = i + 1; j < subGraphRouts.size(); j++) {
- if (subGraphRouts[j] == subGraphRouts[i]) {
- same = j;
- } else {
- break;
- }
- }
- if (i == same) {
- std::fprintf(file, "\n Starting time point = %d | correlations >= %0.2f | tau = %d | window = %d\n\n", static_cast< int >(i) * tauStep, static_cast< double >(crlUpBorder), tau, window);
- } else {
- std::fprintf(file, "\n Starting time point = [%d ; %d] | correlations >= %0.2f | tau = %d | window = %d\n\n", static_cast< int >(i) * tauStep, static_cast< int >(same) * tauStep, static_cast< double >(crlUpBorder), tau, window);
- }
- for (unsigned int j = 0; j < subGraphRouts[i].size(); j++) {
- for (unsigned int k = 0; k < subGraphRouts[i][j].size(); k++) {
- std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.05\n", index[subGraphRouts[i][j][k].first] + 1, index[subGraphRouts[i][j][k].second] + 1);
- }
- std::fprintf(file, "\n");
- }
- table.resize(0);
-
- for (unsigned int j = 0; j < subGraphRouts[i].size(); j++) {
- for (unsigned int k = 0; k < subGraphRouts[i][j].size(); k++) {
- bool flag1 = true, flag2 = true;
- for (unsigned int m = 0; m < table.size(); m++) {
- if (std::get<0>(table[m]) == index[subGraphRouts[i][j][k].first]) {
- std::get<1>(table[m])++;
- std::get<2>(table[m]).push_back(index[subGraphRouts[i][j][k].second]);
- flag1 = false;
- }
- if (std::get<0>(table[m]) == index[subGraphRouts[i][j][k].second]) {
- std::get<1>(table[m])++;
- std::get<2>(table[m]).push_back(index[subGraphRouts[i][j][k].first]);
- flag2 = false;
- }
- }
- if (flag1) {
- a.resize(0);
- a.push_back(index[subGraphRouts[i][j][k].second]);
- table.push_back(std::make_tuple(index[subGraphRouts[i][j][k].first], 1, a));
- }
- if (flag2) {
- a.resize(0);
- a.push_back(index[subGraphRouts[i][j][k].first]);
- table.push_back(std::make_tuple(index[subGraphRouts[i][j][k].second], 1, a));
- }
- }
- }
-
- for (unsigned int j = 0; j < table.size(); j++) {
- std::fprintf(file, "residue %s connections %d | ", (resNames[static_cast< unsigned int >(std::find(index.begin(), index.end(), std::get<0>(table[j])) - index.begin())]).c_str(), std::get<1>(table[j]));
- for (unsigned int k = 0; k < std::get<2>(table[j]).size(); k++) {
- std::fprintf(file, "%s ", (resNames[static_cast< unsigned int >(std::find(index.begin(), index.end(), std::get<2>(table[j])[k]) - index.begin())]).c_str());
- }
- std::fprintf(file, "\n");
- }
- if (pre != 0) {
- std::vector< std::vector< std::pair< unsigned int, unsigned int > > > temp01, temp02;
- std::vector< std::pair< unsigned int, unsigned int > > temp03, temp04, temp05;
- temp01 = subGraphRouts[pre];
- temp02 = subGraphRouts[same];
- temp03.resize(0);
- temp04.resize(0);
- for (unsigned int j = 0; j < temp01.size(); j++) {
- for (unsigned int k = 0; k < temp01[j].size(); k++) {
- temp03.push_back(temp01[j][k]);
- }
- }
- for (unsigned int j = 0; j < temp02.size(); j++) {
- for (unsigned int k = 0; k < temp02[j].size(); k++) {
- temp04.push_back(temp02[j][k]);
- }
- }
- std::sort(temp03.begin(), temp03.end());
- std::sort(temp04.begin(), temp04.end());
- temp05.resize(0);
- std::set_difference(temp03.begin(), temp03.end(), temp04.begin(), temp04.end(), std::inserter(temp05, temp05.begin()));
- std::fprintf(file, "minus:\n");
- for (unsigned int j = 0; j < temp05.size(); j++) {
- std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.05\n", index[temp05[j].first] + 1, index[temp05[j].second] + 1);
- }
- temp05.resize(0);
- std::set_difference(temp04.begin(), temp04.end(), temp03.begin(), temp03.end(), std::inserter(temp05, temp05.begin()));
- std::fprintf(file, "plus:\n");
- for (unsigned int j = 0; j < temp05.size(); j++) {
- std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.05\n", index[temp05[j].first] + 1, index[temp05[j].second] + 1);
- }
- }
- pre = same;
- i = same;
- }
- std::fclose(file);
- }
-};
-
/*! \brief
* \ingroup module_trajectoryanalysis
*/
biod.pnpi.spb.ru / alexxy / gromacs-spacetimecorr.git / commitdiff