src/corrtype.cpp

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