src/corrtype.cpp

   1 #include "corrtype.h"
   2
   3 // деструктор класса
   4 correlationType::~correlationType() {}
   5
   6 // конструктор класса для инициализации
   7 correlationType::correlationType() {}
   8
   9 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) {
  10     setDefaults(ref, wnd, taau, tau_st, crlUp, effRad, mod, out, sels, rsNames);
  11 }
  12
  13 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) {
  14     if (ref.size() != 0) {reference = ref;}
  15     if (wnd != -1) {window = wnd;}
  16     if (taau != -1) {tau = taau;}
  17     if (tau_st != -1) {tauStep = tau_st;}
  18     if (crlUp != -1) {crlUpBorder = crlUp;}
  19     if (effRad != -1) {effRadius = effRad;}
  20     if (mod != -1) {mode = mod;}
  21     if (out != "") {outputName = out;}
  22     if (sels.size() != 0) {selections = sels;}
  23     if (rsNames.size() > 0) {resNames = rsNames;}
  24     subGraphRouts.resize(0);
  25     trajectoryPartition();
  26 }
  27
  28 void correlationType::update(int frame, std::vector< RVec > curFrame) {
  29     trajectory.push_back(curFrame);
  30     if (trajectory.size() == static_cast< unsigned int>(window + tau + 1)) {
  31         trajectory.erase(trajectory.begin());
  32     }
  33     if ((frame + 1 - static_cast< int >(trajectory.size())) % tauStep == 0) {
  34         correlationEval();
  35         selections.erase(selections.begin());
  36         subGraphRouts.resize(subGraphRouts.size() + 1);
  37         graphCalculations(1, static_cast< unsigned int >(tau));
  38         graphBackBoneEvaluation();
  39     }
  40 }
  41
  42 void correlationType::printData() {
  43     printOutputData();
  44 }
  45
  46 void correlationType::trajectoryPartition() {
  47     std::vector< bool > temp1;
  48     std::pair< unsigned int, unsigned int > temp2;
  49     float temp3;
  50     for (unsigned k = 0; k < selections.size(); k++) {
  51         temp1.resize(0); temp1.resize(index.size(), true);
  52         for (unsigned int i = 0; i < selections[k].size(); i++) {
  53             for (unsigned int j = 0; j < selections[k][i].size(); j++) {
  54                 temp1[selections[k][i][j]] = false;
  55             }
  56         }
  57         temp3 = 9000;
  58         for (unsigned i = 0; i < temp1.size(); i++) {
  59             if (temp1[i]) {
  60                 for (unsigned int i = 0; i < selections[k].size(); i++) {
  61                     for (unsigned int j = 0; j < selections[k][i].size(); j++) {
  62                         if (temp3 > (reference[selections[k][i][j]] - reference[i]).norm()) {
  63                             temp3 = (reference[selections[k][i][j]] - reference[i]).norm();
  64                             temp2 = std::make_pair(i, j);
  65                         }
  66                     }
  67                 }
  68                 selections[k][temp2.first].push_back(i);
  69             }
  70         }
  71     }
  72 }
  73
  74 void correlationType::readWriteCorrelations(int rwMode) {
  75     FILE *file;
  76     if (rwMode == 1) {
  77         file = std::fopen((outputName + "-matrixData").c_str(), "a");
  78         for (unsigned int i = 0; i < matrixes.size(); i++) {
  79             std::fprintf(file, "%d %d\n", count, i);
  80             for (unsigned int j = 0; j < matrixes[i].size(); j++) {
  81                 for (unsigned int f = 0; f < matrixes[i][j].size(); f++) {
  82                     std::fprintf(file, "%.4f ", matrixes[i][j][f]); //~16
  83                 }
  84                 std::fprintf(file, "\n");
  85             }
  86         }
  87         std::fclose(file);
  88       }
  89     if (rwMode == 0) {
  90         file = std::fopen((outputName + "-matrixData").c_str(), "r+");
  91         matrixes.resize(0); matrixes.resize(static_cast< unsigned int >(tau + 1));
  92         for (unsigned int i = 0; i < static_cast< unsigned int >(tau + 1); i++) {
  93             int t0, t1, t2 = std::fscanf(file, "%d %d\n", &t0, &t1);
  94             matrixes[i].resize(0); matrixes[i].resize(index.size());
  95             for (unsigned int j = 0; j < index.size(); j++) {
  96                 matrixes[i][j].resize(index.size());
  97                 for (unsigned int k = 0; k < index.size(); k++) {
  98                     t2 = std::fscanf(file, "%lf ", &matrixes[i][j][k]);
  99                 }
 100             }
 101         }
 102     }
 103 }
 104
 105 void correlationType::correlationEval() {
 106     /*
 107      * fitting block / we add spare atoms to existing domains based on proximity
 108      */
 109     std::vector< std::vector< std::pair< unsigned int, unsigned int > > > pairs;
 110     pairs.resize(0); pairs.resize(selections.front().size());
 111     for (unsigned int i = 0; i < selections.front().size(); i++) {
 112         pairs[i].resize(0);
 113         for (unsigned int j = 0; j < selections.front()[i].size(); j++) {
 114             pairs[i].push_back(std::make_pair(selections.front()[i][j], selections.front()[i][j]));
 115         }
 116     }
 117     fitTrajectory.resize(0); fitTrajectory.resize(selections.front().size(), trajectory);
 118     #pragma omp parallel for schedule(dynamic) firstprivate(reference)
 119     for (unsigned long i = 0; i < selections.front().size(); i++) {
 120         for (unsigned int j = 0; j < fitTrajectory[i].size(); j++) {
 121             MyFitNew(reference, fitTrajectory[i][j], pairs[i], 0);
 122         }
 123     }
 124     #pragma omp barrier
 125     frankensteinTrajectory = trajectory;
 126     for (unsigned int i = 0; i < selections.front().size(); i++) {
 127         for (unsigned int j = 0; j < selections.front()[i].size(); j++) {
 128             for (unsigned int k = 0; k < fitTrajectory[i].size(); k++) {
 129                 frankensteinTrajectory[k][selections.front()[i][j]] = fitTrajectory[i][k][selections.front()[i][j]];
 130             }
 131         }
 132     }
 133     /*
 134      * fitting block
 135      */
 136     matrixes.resize(static_cast< unsigned int >(tau + 1));
 137     for (unsigned int i = 0; i < matrixes.size(); i++) {
 138         matrixes[i].resize(index.size());
 139         for (unsigned int j = 0; j < index.size(); j++) {
 140             matrixes[i][j].resize(index.size());
 141             for (unsigned int k = 0; k < index.size(); k++) {
 142                 matrixes[i][j][k] = 0.;
 143             }
 144         }
 145     }
 146     std::vector< std::vector< double > > a, b, c;
 147     std::vector< double > d;
 148     #pragma omp parallel for ordered schedule(dynamic) shared(matrixes) firstprivate(frankensteinTrajectory, reference)
 149     for (unsigned int i = 0; i <= static_cast< unsigned int >(tau); i++) {
 150         a.resize(0); b.resize(0); c.resize(0); d.resize(0);
 151         for (unsigned int j = 0; j < index.size(); j++) {
 152             a.push_back(d); b.push_back(d); c.push_back(d);
 153             for (unsigned int k = 0; k < index.size(); k++) {
 154                 a[j].push_back(0.); b[j].push_back(0.); c[j].push_back(0.);
 155             }
 156         }
 157         for (unsigned int j = 0; j < static_cast< unsigned int >(window); j++) {
 158             for (unsigned int k1 = 0; k1 < index.size(); k1++) {
 159                 for (unsigned int k2 = 0; k2 < index.size(); k2++) {
 160                     RVec temp1, temp2;
 161                     temp1 = frankensteinTrajectory[j][k1]       - reference[k1];
 162                     temp2 = frankensteinTrajectory[j + i][k2]   - reference[k2];
 163                     a[k1][k2] +=   (static_cast<double>(temp1[0]) * static_cast<double>(temp2[0]) +
 164                                     static_cast<double>(temp1[1]) * static_cast<double>(temp2[1]) +
 165                                     static_cast<double>(temp1[2]) * static_cast<double>(temp2[2]));
 166                     b[k1][k2] +=   (static_cast<double>(temp1[0]) * static_cast<double>(temp1[0]) +
 167                                     static_cast<double>(temp1[1]) * static_cast<double>(temp1[1]) +
 168                                     static_cast<double>(temp1[2]) * static_cast<double>(temp1[2]));
 169                     c[k1][k2] +=   (static_cast<double>(temp2[0]) * static_cast<double>(temp2[0]) +
 170                                     static_cast<double>(temp2[1]) * static_cast<double>(temp2[1]) +
 171                                     static_cast<double>(temp2[2]) * static_cast<double>(temp2[2]));
 172                 }
 173             }
 174         }
 175         for (unsigned int j = 0; j < index.size(); j++) {
 176             for (unsigned int k = 0; k < index.size(); k++) {
 177                 matrixes[i][j][k] = a[j][k] / (std::sqrt(b[j][k] * c[j][k]));
 178             }
 179         }
 180     }
 181     #pragma omp barrier
 182     for (unsigned int i = 0; i < matrixes.size(); i++) {
 183         for (unsigned int j = 0; j < matrixes[i].size(); j++) {
 184             for (unsigned int k = 0; k < matrixes[i][j].size(); k++) {
 185                 matrixes[i][j][k] = std::round(matrixes[i][j][k] * 10000) / 10000;
 186             }
 187         }
 188     }
 189     readWriteCorrelations(1);
 190     count++;
 191 }
 192
 193 void correlationType::graphCalculations(unsigned int tauStart, unsigned int tauEnd) {
 194     graph.resize(0); graph.resize(index.size());
 195     for (unsigned int i = 0; i < index.size(); i++) {
 196         graph[i].resize(index.size(), std::make_pair(0, -1));
 197     }
 198     RVec temp;
 199     for (unsigned int i = tauStart; i <= tauEnd; i++) {
 200         for (unsigned int j = 0; j < index.size(); j++) {
 201             for (unsigned int k = j; k < index.size(); k++) {
 202                 temp = reference[i] - reference[k];
 203                 if (std::max(std::abs(matrixes[i][j][k]), std::abs(matrixes[i][k][j])) >= static_cast< double >(crlUpBorder) &&
 204                     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]))) {
 205                     if (std::abs(matrixes[i][j][k]) > std::abs(matrixes[i][k][j])) {
 206                         graph[j][k].first = matrixes[i][j][k];
 207                     } else {
 208                         graph[j][k].first = matrixes[i][k][j];
 209                     }
 210                     graph[j][k].second = static_cast< int >(i);
 211                 }
 212             }
 213         }
 214     }
 215     std::vector< bool > graph_flags;
 216     graph_flags.resize(0); graph_flags.resize(index.size(), true);
 217     std::vector< unsigned int > a;
 218     std::vector< std::pair< unsigned int, unsigned int > > b;
 219     a.resize(0); b.resize(0);
 220     std::vector< unsigned int > width1, width2, tempSubGraph;
 221     for (unsigned int i = 0; i < index.size(); i++) {
 222         if (graph_flags[i]) {
 223             subGraphPoints.push_back(a);
 224             subGraphRbr.push_back(b);
 225             width1.resize(0); width2.resize(0); tempSubGraph.resize(0);
 226             width1.push_back(i);
 227             tempSubGraph.push_back(i);
 228             graph_flags[i] = false;
 229             while(width1.size() > 0) {
 230                 width2.clear();
 231                 for (unsigned int j = 0; j < width1.size(); j++) {
 232                     for (unsigned int k = 0; k < index.size(); k++) {
 233                         if (graph[width1[j]][k].second > -1 && graph_flags[k]) {
 234                             width2.push_back(k);
 235                             graph_flags[k] = false;
 236                         }
 237                     }
 238                 }
 239                 width1 = width2;
 240                 for (unsigned int j = 0; j < width2.size(); j++) {
 241                     tempSubGraph.push_back(width2[j]);
 242                 }
 243             }
 244             subGraphPoints.back() = tempSubGraph;
 245             for (unsigned int j = 0; j < tempSubGraph.size(); j++) {
 246                 for (unsigned int k = 0; k < index.size(); k++) {
 247                     if (graph[tempSubGraph[j]][k].second > -1) {
 248                         subGraphRbr.back().push_back(std::make_pair(tempSubGraph[j], k));
 249                     }
 250                 }
 251             }
 252         }
 253     }
 254 }
 255
 256 bool correlationType::myComparisonFunction (const std::pair< int, double > i, const std::pair< int, double > j) {
 257     return i.second < j.second;
 258 }
 259
 260 void correlationType::graphBackBoneEvaluation() {
 261     std::vector< double >                                   key;
 262     std::vector< long >                                     path;
 263     std::vector< std::pair< unsigned int, double > >        que;
 264     std::vector< std::pair< unsigned int, unsigned int > >  a;
 265     unsigned int                                            v;
 266     a.resize(0);
 267     subGraphRouts.resize(0);
 268     for (unsigned int i = 0; i < subGraphPoints.size(); i++) {
 269         key.resize(0);
 270         path.resize(0);
 271         que.resize(0);
 272         v = 0;
 273         if (subGraphPoints[i].size() > 2) {
 274             key.resize(index.size(), 2);
 275             path.resize(index.size(), -1);
 276             key[subGraphPoints[i][0]] = 0;
 277             for (unsigned int j = 0; j < subGraphPoints[i].size(); j++) {
 278                 que.push_back(std::make_pair(subGraphPoints[i][j], key[subGraphPoints[i][j]]));
 279             }
 280             std::sort(que.begin(), que.end(), myComparisonFunction);
 281             while (!que.empty()) {
 282                 v = que[0].first;
 283                 que.erase(que.begin());
 284                 for (unsigned int j = 0; j < subGraphRbr[i].size(); j++) {
 285                     long u = -1;
 286                     if (subGraphRbr[i][j].first == v) {
 287                         u = subGraphRbr[i][j].second;
 288                     } else if (subGraphRbr[i][j].second == v) {
 289                         u = subGraphRbr[i][j].first;
 290                     }
 291                     bool flag = false;
 292                     unsigned long pos = 0;
 293                     for (unsigned long k = 0; k < que.size(); k++) {
 294                         if (que[k].first == u) {
 295                             flag = true;
 296                             pos = k;
 297                             k = que.size() + 1;
 298                         }
 299                     }
 300                     if (flag && key[static_cast< unsigned long >(u)] > 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first)) {
 301                         path[static_cast< unsigned long >(u)] = v;
 302                         key[static_cast< unsigned long >(u)] = 1 - std::abs(graph[v][static_cast< unsigned long >(u)].first);
 303                         que[pos].second = key[static_cast< unsigned long >(u)];
 304                         sort(que.begin(), que.end(), myComparisonFunction);
 305                     }
 306                 }
 307             }
 308             subGraphRouts.back().push_back(a);
 309             for (unsigned int j = 0; j < index.size(); j++) {
 310                 if (path[j] != -1) {
 311                     subGraphRouts.back().back().push_back(std::make_pair(j, path[j]));
 312                 }
 313             }
 314         }
 315     }
 316 }
 317
 318 void correlationType::printOutputData() {
 319     FILE *file;
 320     file = std::fopen((outputName + "-arrowsData.txt").c_str(), "w+");
 321     unsigned int same, pre = 0;
 322     std::vector< std::tuple< int, int, std::vector< int > > > table;
 323     std::vector< int > a;
 324     for (unsigned int i = 0; i < subGraphRouts.size(); i++) {
 325         same = i;
 326         for (unsigned int j = i + 1; j < subGraphRouts.size(); j++) {
 327             if (subGraphRouts[j] == subGraphRouts[i]) {
 328                 same = j;
 329             } else {
 330                 break;
 331             }
 332         }
 333         if (i == same) {
 334             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);
 335         } else {
 336             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);
 337         }
 338         for (unsigned int j = 0; j < subGraphRouts[i].size(); j++) {
 339             for (unsigned int k = 0; k < subGraphRouts[i][j].size(); k++) {
 340                 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);
 341             }
 342             std::fprintf(file, "\n");
 343         }
 344         table.resize(0);
 345         for (unsigned int j = 0; j < subGraphRouts[i].size(); j++) {
 346             for (unsigned int k = 0; k < subGraphRouts[i][j].size(); k++) {
 347                 bool flag1 = true, flag2 = true;
 348                 for (unsigned int m = 0; m < table.size(); m++) {
 349                     if (std::get<0>(table[m]) == index[subGraphRouts[i][j][k].first]) {
 350                         std::get<1>(table[m])++;
 351                         std::get<2>(table[m]).push_back(index[subGraphRouts[i][j][k].second]);
 352                         flag1 = false;
 353                     }
 354                     if (std::get<0>(table[m]) == index[subGraphRouts[i][j][k].second]) {
 355                         std::get<1>(table[m])++;
 356                         std::get<2>(table[m]).push_back(index[subGraphRouts[i][j][k].first]);
 357                         flag2 = false;
 358                     }
 359                 }
 360                 if (flag1) {
 361                     a.resize(0);
 362                     a.push_back(index[subGraphRouts[i][j][k].second]);
 363                     table.push_back(std::make_tuple(index[subGraphRouts[i][j][k].first], 1, a));
 364                 }
 365                 if (flag2) {
 366                     a.resize(0);
 367                     a.push_back(index[subGraphRouts[i][j][k].first]);
 368                     table.push_back(std::make_tuple(index[subGraphRouts[i][j][k].second], 1, a));
 369                 }
 370             }
 371         }
 372         for (unsigned int j = 0; j < table.size(); j++) {
 373             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]));
 374             for (unsigned int k = 0; k < std::get<2>(table[j]).size(); k++) {
 375                 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());
 376             }
 377             std::fprintf(file, "\n");
 378         }
 379         if (pre != 0) {
 380             std::vector< std::vector< std::pair< unsigned int, unsigned int > > > temp01, temp02;
 381             std::vector< std::pair< unsigned int, unsigned int > > temp03, temp04, temp05;
 382             temp01 = subGraphRouts[pre];
 383             temp02 = subGraphRouts[same];
 384             temp03.resize(0);
 385             temp04.resize(0);
 386             for (unsigned int j = 0; j < temp01.size(); j++) {
 387                 for (unsigned int k = 0; k < temp01[j].size(); k++) {
 388                     temp03.push_back(temp01[j][k]);
 389                 }
 390             }
 391             for (unsigned int j = 0; j < temp02.size(); j++) {
 392                 for (unsigned int k = 0; k < temp02[j].size(); k++) {
 393                     temp04.push_back(temp02[j][k]);
 394                 }
 395             }
 396             std::sort(temp03.begin(), temp03.end());
 397             std::sort(temp04.begin(), temp04.end());
 398             temp05.resize(0);
 399             std::set_difference(temp03.begin(), temp03.end(), temp04.begin(), temp04.end(), std::inserter(temp05, temp05.begin()));
 400             std::fprintf(file, "minus:\n");
 401             for (unsigned int j = 0; j < temp05.size(); j++) {
 402                 std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.05\n", index[temp05[j].first] + 1, index[temp05[j].second] + 1);
 403             }
 404             temp05.resize(0);
 405             std::set_difference(temp04.begin(), temp04.end(), temp03.begin(), temp03.end(), std::inserter(temp05, temp05.begin()));
 406             std::fprintf(file, "plus:\n");
 407             for (unsigned int j = 0; j < temp05.size(); j++) {
 408                 std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.05\n", index[temp05[j].first] + 1, index[temp05[j].second] + 1);
 409             }
 410         }
 411         pre = same;
 412         i = same;
 413     }
 414     std::fclose(file);
 415 }