- changed to prefix "++"
[alexxy/gromacs-spacetimecorr.git] / src / corrtype.cpp
index b5739fa179de254f5dfa7891b42b37cbc24d09ce..1722fe909f1670802b63f1de5b318ccf9920293e 100644 (file)
@@ -3,36 +3,49 @@
 // деструктор класса
 correlationType::~correlationType() {}
 
-// конструктор класса для инициализации
-correlationType::correlationType() {}
-
-correlationType::correlationType(const std::vector< RVec > &ref, int wnd, int taau, int tau_st, float crlUp, float effRad, int mod,
-                                 const std::string &out, const std::vector< int > &indx,
-                                 const std::vector< std::vector < std::vector < size_t > > > &sels,
-                                 const std::vector< std::string > &rsNames) {
-    setDefaults(ref, wnd, taau, tau_st, crlUp, effRad, mod, out, indx, sels, rsNames);
+// конструктор класса с необходимыми параметрами
+correlationType::correlationType(const std::vector< RVec > &inputReference, const int inputWindow, const int inputTau,
+                                 const int inputTauStep, const float inputCrlUpBorder, const float inputEffRadius, const int inputMode,
+                                 const std::string &inputOutputName, const std::vector< size_t > &inputIndex,
+                                 const std::vector< std::vector < std::vector < size_t > > > &inputSelections,
+                                 const std::vector< std::string > &inputResNames) {
+    setDefaults(inputReference, inputWindow, inputTau, inputTauStep, inputCrlUpBorder, inputEffRadius, inputMode,
+                inputOutputName, inputIndex, inputSelections, inputResNames);
 }
 
-void correlationType::setDefaults(const std::vector< RVec > &ref, int wnd, int taau, int tau_st, float crlUp, float effRad, int mod,
-                                  const std::string &out, const std::vector< int > &indx,
-                                  const std::vector< std::vector < std::vector < size_t > > > &sels,
-                                  const 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;}
-    if (indx.size() > 0) {index = indx;}
+// функция заполнения необходимых параметров
+void correlationType::setDefaults(const std::vector< RVec > &inputReference, const int inputWindow, const int inputTau,
+                                  const int inputTauStep, const float inputCrlUpBorder, const float inputEffRadius, const int inputMode,
+                                  const std::string &inputOutputName, const std::vector< size_t > &inputIndex,
+                                  const std::vector< std::vector < std::vector < size_t > > > &inputSelections,
+                                  const std::vector< std::string > &inputResNames) {
+    if (inputReference.size() == 0 || inputIndex.size() == 0 || inputSelections.size() == 0 || inputResNames.size() == 0) {
+        throw "\nempty vectors\n";
+    }
+    reference = inputReference;
+    inputWindow > 0 ? window = inputWindow : window = 1000;
+    inputTau > 0 ? tau = inputTau : tau = window / 2;
+    inputTauStep > 0 ? tauStep = inputTauStep : tauStep = window / 10;
+    inputCrlUpBorder > 0 ? crlUpBorder = inputCrlUpBorder : crlUpBorder = 0.5;
+    inputEffRadius > 0 ? effRadius = inputEffRadius : effRadius = 1.5;
+    inputMode == 0 || inputMode == 1 ? mode = inputMode : mode = 1;
+    inputOutputName.size() > 0 ? outputName = inputOutputName : outputName = std::chrono::system_clock::to_time_t(std::chrono::system_clock::now());
+    index = inputIndex;
+    selections = inputSelections;
+    if (selections.front().size() == 0) {
+        selections.front().push_back(index);
+    }
+    for (size_t i {1}; i < selections.size(); ++i) {
+        if (selections[i].size() == 0) {
+            selections[i] = selections[i - 1];
+        }
+    }
+    resNames = inputResNames;
     subGraphRouts.resize(0);
     trajectoryPartition();
 }
 
+// функция обновления данных для подсчёта корреляций
 void correlationType::update(const int frameNum, const std::vector< RVec > &curFrame) {
     trajectory.push_back(curFrame);
     int temp = window + tau;
@@ -42,14 +55,17 @@ void correlationType::update(const int frameNum, const std::vector< RVec > &curF
     if (mode == 1) {
         if ( ((frameNum - temp + 1) >= 0) && ((frameNum - temp + 1) % tauStep == 0)) {
             correlationEval();
-            selections.erase(selections.begin());
+            if (selections.size() > 1) {
+                selections.erase(selections.begin());
+            }
         }
     }
 }
 
+// функция считывания подсчитанных данных из класса
 void correlationType::readEval(size_t frameNum) {
     if (mode == 0) {
-        for (size_t i {0}; i <= (frameNum + 1 - tau - window) / tauStep; i++) {
+        for (size_t i {0}; i <= (frameNum + 1 - tau - window) / tauStep; ++i) {
             readWriteCorrelations(0);
             subGraphRouts.resize(subGraphRouts.size() + 1);
             graphCalculations(1, static_cast< size_t >(tau));
@@ -58,22 +74,24 @@ void correlationType::readEval(size_t frameNum) {
     }
 }
 
+// функция вывода графов в виде стрелок в тектовом формате в файл
 void correlationType::printData() {
     printOutputData();
 }
 
+// функция дополнения SELECTIONS до полного покрытия объекта
 void correlationType::trajectoryPartition() {
     std::vector< bool > temp1;
     std::pair< size_t, size_t > temp2;
     float temp3;
     std::vector< std::vector < std::vector < size_t > > > selectionsTemp;
     selectionsTemp.resize(selections.size());
-    for (size_t i1 {0}; i1 < selections.size(); i1++) {
+    for (size_t i1 {0}; i1 < selections.size(); ++i1) {
         selectionsTemp[i1].resize(selections[i1].size());
-        for (size_t i2 {0}; i2 < selections[i1].size(); i2++) {
+        for (size_t i2 {0}; i2 < selections[i1].size(); ++i2) {
             selectionsTemp[i1][i2].resize(0);
-            for (size_t i3 {0}; i3 < selections[i1][i2].size(); i3++) {
-                for (size_t i4 {0}; i4 < index.size(); i4++) {
+            for (size_t i3 {0}; i3 < selections[i1][i2].size(); ++i3) {
+                for (size_t i4 {0}; i4 < index.size(); ++i4) {
                     if (selections[i1][i2][i3] == index[i4]) {
                         selectionsTemp[i1][i2].push_back(i4);
                         break;
@@ -91,11 +109,11 @@ void correlationType::trajectoryPartition() {
                 temp1[j] = false;
             }
         }
-        for (size_t i {0}; i < temp1.size(); i++) {
+        for (size_t i {0}; i < temp1.size(); ++i) {
             if (temp1[i]) {
                 temp3 = std::numeric_limits<float>::max();
-                for (size_t f {0}; f < k.size(); f++) {
-                    for (size_t j {0}; j < k[f].size(); j++) {
+                for (size_t f {0}; f < k.size(); ++f) {
+                    for (size_t j {0}; j < k[f].size(); ++j) {
                         if (float temp4 {(reference[k[f][j]] - reference[i]).norm()}; temp3 > temp4) {
                             temp3 = temp4;
                             temp2 = std::make_pair(f, j);
@@ -108,14 +126,15 @@ void correlationType::trajectoryPartition() {
     }
 }
 
+// функция чтения/записи корреляций
 void correlationType::readWriteCorrelations(int rwMode) {
     FILE *file;
     if (rwMode == 1) {
         file = std::fopen((outputName + "-matrixData").c_str(), "a");
-        for (size_t i {0}; i < matrixes.size(); i++) {
+        for (size_t i {0}; i < matrixes.size(); ++i) {
             std::fprintf(file, "%d %lu\n", count, i);
-            for (size_t j {0}; j < matrixes[i].size(); j++) {
-                for (size_t f {0}; f < matrixes[i][j].size(); f++) {
+            for (size_t j {0}; j < matrixes[i].size(); ++j) {
+                for (size_t f {0}; f < matrixes[i][j].size(); ++f) {
                     std::fprintf(file, "%.4f ", matrixes[i][j][f]); //~16
                 }
                 std::fprintf(file, "\n");
@@ -127,13 +146,13 @@ void correlationType::readWriteCorrelations(int rwMode) {
         file = std::fopen((outputName + "-matrixData").c_str(), "r+");
         matrixes.resize(0);
         matrixes.resize(static_cast< unsigned int >(tau + 1));
-        for (size_t i {0}; i < static_cast< size_t >(tau + 1); i++) {
+        for (size_t i {0}; i < static_cast< size_t >(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 (size_t j {0}; j < index.size(); j++) {
+            for (size_t j {0}; j < index.size(); ++j) {
                 matrixes[i][j].resize(index.size(), 0.);
-                for (size_t k {0}; k < index.size(); k++) {
+                for (size_t k {0}; k < index.size(); ++k) {
                     t2 = std::fscanf(file, "%lf ", &matrixes[i][j][k]);
                 }
             }
@@ -141,34 +160,36 @@ void correlationType::readWriteCorrelations(int rwMode) {
     }
 }
 
+// фитирование траектории на основе SELECTIONS и объеденение её воедино
 inline void correlationType::trajectoryFitting() {
     std::vector< std::vector< std::pair< size_t, size_t > > > pairs;
     pairs.resize(0);
     pairs.resize(selections.front().size());
-    for (size_t i {0}; i < selections.front().size(); i++) {
+    for (size_t i {0}; i < selections.front().size(); ++i) {
         pairs[i].resize(0);
-        for (size_t j {0}; j < selections.front()[i].size(); j++) {
+        for (size_t 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 (size_t i = 0; i < selections.front().size(); i++) {
-        for (size_t j {0}; j < fitTrajectory[i].size(); j++) {
+    for (size_t i = 0; i < selections.front().size(); ++i) {
+        for (size_t j {0}; j < fitTrajectory[i].size(); ++j) {
             MyFitNew(reference, fitTrajectory[i][j], pairs[i], 0);
         }
     }
     #pragma omp barrier
-    for (size_t i {0}; i < selections.front().size(); i++) {
-        for (size_t j {0}; j < selections.front()[i].size(); j++) {
-            for (size_t k {0}; k < fitTrajectory[i].size(); k++) {
+    for (size_t i {0}; i < selections.front().size(); ++i) {
+        for (size_t j {0}; j < selections.front()[i].size(); ++j) {
+            for (size_t k {0}; k < fitTrajectory[i].size(); ++k) {
                 trajectory[k][selections.front()[i][j]] = fitTrajectory[i][k][selections.front()[i][j]];
             }
         }
     }
 }
 
+// инициализация матрицы корреляций
 inline void correlationType::matrixNullFitting() {
     matrixes.resize(0);
     matrixes.resize(static_cast< size_t >(tau + 1));
@@ -185,11 +206,12 @@ inline void correlationType::matrixNullFitting() {
     }
 }
 
+// функция подсчёта корреляций
 void correlationType::correlationEval() {
     trajectoryFitting();
     matrixNullFitting();
     #pragma omp parallel for ordered schedule(dynamic) shared(matrixes) firstprivate(trajectory, reference)
-    for (size_t i = 0; i <= static_cast< size_t >(tau); i++) {
+    for (size_t i = 0; i <= static_cast< size_t >(tau); ++i) {
         std::vector< std::vector< double > > a, b, c;
         std::vector< double > d;
         d.resize(index.size(), 0.);
@@ -199,9 +221,9 @@ void correlationType::correlationEval() {
         a.resize(index.size(), d);
         b.resize(index.size(), d);
         c.resize(index.size(), d);
-        for (size_t j {0}; j < static_cast< size_t >(window); j++) {
-            for (size_t k1 {0}; k1 < index.size(); k1++) {
-                for (size_t k2 {0}; k2 < index.size(); k2++) {
+        for (size_t j {0}; j < static_cast< size_t >(window); ++j) {
+            for (size_t k1 {0}; k1 < index.size(); ++k1) {
+                for (size_t k2 {0}; k2 < index.size(); ++k2) {
                     RVec temp1, temp2;
                     temp1 = trajectory[j][k1]     - reference[k1];
                     temp2 = trajectory[j + i][k2] - reference[k2];
@@ -217,36 +239,37 @@ void correlationType::correlationEval() {
                 }
             }
         }
-        for (size_t j {0}; j < index.size(); j++) {
-            for (size_t k {0}; k < index.size(); k++) {
+        for (size_t j {0}; j < index.size(); ++j) {
+            for (size_t 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 (size_t i {0}; i < matrixes.size(); i++) {
-        for (size_t j {0}; j < matrixes[i].size(); j++) {
-            for (size_t k {0}; k < matrixes[i][j].size(); k++) {
+    for (size_t i {0}; i < matrixes.size(); ++i) {
+        for (size_t j {0}; j < matrixes[i].size(); ++j) {
+            for (size_t k {0}; k < matrixes[i][j].size(); ++k) {
                 matrixes[i][j][k] = std::round(matrixes[i][j][k] * 10000) / 10000;
             }
         }
     }
     readWriteCorrelations(1);
-    count++;
+    ++count;
 }
 
+// функция нахождения графов взаимодействия
 void correlationType::graphCalculations(size_t tauStart, size_t tauEnd) {
     graph.resize(0);
     graph.resize(index.size());
     subGraphPoints.resize(0);
     subGraphRbr.resize(0);
-    for (size_t i {0}; i < index.size(); i++) {
+    for (size_t i {0}; i < index.size(); ++i) {
         graph[i].resize(index.size(), std::make_pair(0, -1));
     }
     RVec temp;
-    for (size_t i {tauStart}; i <= tauEnd; i++) {
-        for (size_t j {0}; j < index.size(); j++) {
-            for (size_t k {j}; k < index.size(); k++) {
+    for (size_t i {tauStart}; i <= tauEnd; ++i) {
+        for (size_t j {0}; j < index.size(); ++j) {
+            for (size_t k {j}; k < index.size(); ++k) {
                 temp = reference[j] - reference[k];
                 if (double tempIf {std::max(std::abs(matrixes[i][j][k]), std::abs(matrixes[i][k][j]))}; (tempIf >= static_cast< double >(crlUpBorder)) &&
                     (static_cast< float >(norm(temp)) <= effRadius) && (std::abs(graph[j][k].first) < tempIf)) {
@@ -264,7 +287,7 @@ void correlationType::graphCalculations(size_t tauStart, size_t tauEnd) {
     a.resize(0);
     b.resize(0);
     std::vector< size_t > width1, width2, tempSubGraph;
-    for (size_t i {0}; i < index.size(); i++) {
+    for (size_t i {0}; i < index.size(); ++i) {
         if (graph_flags[i]) {
             subGraphPoints.push_back(a);
             subGraphRbr.push_back(b);
@@ -276,8 +299,8 @@ void correlationType::graphCalculations(size_t tauStart, size_t tauEnd) {
             graph_flags[i] = false;
             while(width1.size() > 0) {
                 width2.resize(0);
-                for (size_t j {0}; j < width1.size(); j++) {
-                    for (size_t k {0}; k < index.size(); k++) {
+                for (size_t j {0}; j < width1.size(); ++j) {
+                    for (size_t k {0}; k < index.size(); ++k) {
                         if ((graph[width1[j]][k].second > -1) && graph_flags[k]) {
                             width2.push_back(k);
                             graph_flags[k] = false;
@@ -285,13 +308,13 @@ void correlationType::graphCalculations(size_t tauStart, size_t tauEnd) {
                     }
                 }
                 width1 = width2;
-                for (size_t j {0}; j < width2.size(); j++) {
+                for (size_t j {0}; j < width2.size(); ++j) {
                     tempSubGraph.push_back(width2[j]);
                 }
             }
             subGraphPoints.back() = tempSubGraph;
-            for (size_t j {0}; j < tempSubGraph.size(); j++) {
-                for (size_t k {0}; k < index.size(); k++) {
+            for (size_t j {0}; j < tempSubGraph.size(); ++j) {
+                for (size_t k {0}; k < index.size(); ++k) {
                     if (graph[tempSubGraph[j]][k].second > -1) {
                         subGraphRbr.back().push_back(std::make_pair(tempSubGraph[j], k));
                     }
@@ -301,10 +324,12 @@ void correlationType::graphCalculations(size_t tauStart, size_t tauEnd) {
     }
 }
 
+// вспомогательная функция сравнения
 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;
@@ -313,7 +338,7 @@ void correlationType::graphBackBoneEvaluation() {
     size_t                                      v;
     a.resize(0);
     subGraphRouts.back().resize(0);
-    for (size_t i {0}; i < subGraphPoints.size(); i++) {
+    for (size_t i {0}; i < subGraphPoints.size(); ++i) {
         key.resize(0);
         path.resize(0);
         que.resize(0);
@@ -322,14 +347,14 @@ void correlationType::graphBackBoneEvaluation() {
             key.resize(index.size(), 2);
             path.resize(index.size(), -1);
             key[subGraphPoints[i][0]] = 0;
-            for (size_t j {0}; j < subGraphPoints[i].size(); j++) {
+            for (size_t 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.front().first;
                 que.erase(que.begin());
-                for (size_t j {0}; j < subGraphRbr[i].size(); j++) {
+                for (size_t j {0}; j < subGraphRbr[i].size(); ++j) {
                     long u {-1};
                     if (subGraphRbr[i][j].first == v) {
                         u = subGraphRbr[i][j].second;
@@ -338,7 +363,7 @@ void correlationType::graphBackBoneEvaluation() {
                     }
                     bool flag {false};
                     size_t pos {0};
-                    for (size_t k {0}; k < que.size(); k++) {
+                    for (size_t k {0}; k < que.size(); ++k) {
                         if (que[k].first == u) {
                             flag = true;
                             pos = k;
@@ -355,7 +380,7 @@ void correlationType::graphBackBoneEvaluation() {
                 }
             }
             subGraphRouts.back().push_back(a);
-            for (size_t j {0}; j < index.size(); j++) {
+            for (size_t j {0}; j < index.size(); ++j) {
                 if (path[j] != -1) {
                     subGraphRouts.back().back().push_back(std::make_pair(path[j], j));
                 }
@@ -364,15 +389,16 @@ void correlationType::graphBackBoneEvaluation() {
     }
 }
 
+// функция вывода графов в текстовом виде в файл
 void correlationType::printOutputData() {
     FILE *file {std::fopen((outputName + "-arrowsData.txt").c_str(), "w+")};
     size_t same, pre {0};
     std::vector< std::tuple< int, int, std::vector< int > > > table;
     table.resize(0);
     std::vector< int > a;
-    for (size_t i {0}; i < subGraphRouts.size(); i++) {
+    for (size_t i {0}; i < subGraphRouts.size(); ++i) {
         same = i;
-        for (size_t j {i + 1}; j < subGraphRouts.size(); j++) {
+        for (size_t j {i + 1}; j < subGraphRouts.size(); ++j) {
             if (subGraphRouts[j] == subGraphRouts[i]) {
                 same = j;
             } else {
@@ -384,9 +410,9 @@ void correlationType::printOutputData() {
         } 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 (size_t j {0}; j < subGraphRouts[i].size(); j++) {
-            for (size_t 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);
+        for (size_t j {0}; j < subGraphRouts[i].size(); ++j) {
+            for (size_t k {0}; k < subGraphRouts[i][j].size(); ++k) {
+                std::fprintf(file, "cgo_arrow (id %3lu), (id %3lu), radius=0.05\n", index[subGraphRouts[i][j][k].first] + 1, index[subGraphRouts[i][j][k].second] + 1);
             }
             std::fprintf(file, "\n");
         }
@@ -394,14 +420,14 @@ void correlationType::printOutputData() {
         for (auto &j : subGraphRouts[i]) {
             for (auto &k : j) {
                 bool flag1 {true}, flag2 {true};
-                for (size_t m {0}; m < table.size(); m++) {
+                for (size_t m {0}; m < table.size(); ++m) {
                     if (std::get<0>(table[m]) == index[k.first]) {
-                        std::get<1>(table[m])++;
+                        ++std::get<1>(table[m]);
                         std::get<2>(table[m]).push_back(index[k.second]);
                         flag1 = false;
                     }
                     if (std::get<0>(table[m]) == index[k.second]) {
-                        std::get<1>(table[m])++;
+                        ++std::get<1>(table[m]);
                         std::get<2>(table[m]).push_back(index[k.first]);
                         flag2 = false;
                     }
@@ -418,9 +444,9 @@ void correlationType::printOutputData() {
                 }
             }
         }
-        for (size_t j {0}; j < table.size(); j++) {
+        for (size_t j {0}; j < table.size(); ++j) {
             std::fprintf(file, "residue %s connections %d | ", (resNames[static_cast< size_t >(std::find(index.begin(), index.end(), std::get<0>(table[j])) - index.begin())]).c_str(), std::get<1>(table[j]));
-            for (size_t k {0}; k < std::get<2>(table[j]).size(); k++) {
+            for (size_t k {0}; k < std::get<2>(table[j]).size(); ++k) {
                 std::fprintf(file, "%s ", (resNames[static_cast< size_t >(std::find(index.begin(), index.end(), std::get<2>(table[j])[k]) - index.begin())]).c_str());
             }
             std::fprintf(file, "\n");
@@ -448,13 +474,13 @@ void correlationType::printOutputData() {
             std::set_difference(temp03.begin(), temp03.end(), temp04.begin(), temp04.end(), std::inserter(temp05, temp05.begin()));
             std::fprintf(file, "minus:\n");
             for (auto &j : temp05) {
-                std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.05\n", index[j.first] + 1, index[j.second] + 1);
+                std::fprintf(file, "cgo_arrow (id %3lu), (id %3lu), radius=0.05\n", index[j.first] + 1, index[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 (auto &j : temp05) {
-                std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.05\n", index[j.first] + 1, index[j.second] + 1);
+                std::fprintf(file, "cgo_arrow (id %3lu), (id %3lu), radius=0.05\n", index[j.first] + 1, index[j.second] + 1);
             }
         }
         pre = same;