src/domaintype.cpp

   1 #include "domaintype.h"
   2
   3 // деструктор класса
   4 domainType::~domainType() {}
   5
   6 // конструктор класса для инициализации
   7 domainType::domainType() {}
   8
   9 // конструктор класса для полноценной инициализации данных
  10 domainType::domainType(const std::vector< size_t > &index, const std::vector< RVec > &reference,
  11                     const int windowSize, const int domainMinimumSize,
  12                     const int domainSearchAlgorythm, const int timeStepBetweenWindowStarts,
  13                     const unsigned int sliceNum, const double epsilon, const double delta,
  14                     const std::string outPutFileName) {
  15     setDefaults(index, reference, windowSize, domainMinimumSize, domainSearchAlgorythm, timeStepBetweenWindowStarts, sliceNum, epsilon, delta, outPutFileName);
  16 }
  17
  18 // set numeric values to "-1" / string value to "" - if you want default settings
  19 // функция заполнения необходимыми данными для вычисления структурных доменов
  20 void domainType::setDefaults(const std::vector< size_t > &index, const std::vector< RVec > &reference,
  21                         const int windowSize, const int domainMinimumSize,
  22                         const int domainSearchAlgorythm, const int timeStepBetweenWindows,
  23                         const unsigned int sliceNum, const double epsilon, const double delta,
  24                         const std::string outPutFileName) {
  25     graph.resize(0);
  26     structIndex = index;
  27     refTable.resize(reference.size());
  28     for (size_t i {0}; i < reference.size(); i++) {
  29         refTable[i].resize(0);
  30         for (size_t j {0}; j < reference.size(); j++) {
  31             refTable[i].push_back(reference[i] - reference[j]);
  32         }
  33     }
  34     // слава шиве, слава индусам... зачем столько ифов?
  35     if (epsilon != -1) { eps = epsilon; }
  36     if (delta != -1) { dlt = delta; }
  37     if (windowSize != -1) { window = windowSize; }
  38     if (domainMinimumSize != -1) { dms = domainMinimumSize; }
  39     if (domainSearchAlgorythm != -1) { dsa = domainSearchAlgorythm; }
  40     if (timeStepBetweenWindows != -1) { ts = timeStepBetweenWindows; }
  41     if (outPutFileName != "") { outPut = outPutFileName; }
  42     domsizes.resize(0);
  43     domsizes.resize(sliceNum);
  44     updateCount = 0;
  45 }
  46
  47 // фукнция обновления данных для выделения структурных доменов
  48 void domainType::update(const std::vector< std::vector< RVec > > &frame, const int frameNumber) {
  49     // громакс считает с нуля, проверял; инициализируем новое "окно"
  50     if (frameNumber % ts == 0) {
  51         graph.resize(graph.size() + 1);
  52         graph.back().resize(structIndex.size() - static_cast< size_t >(dms) + 1);
  53         for (auto &i : graph.back()) {
  54             setGraph(i);
  55         }
  56     }
  57     // заполняем структуру данных для структурных доменов
  58     for (auto &i : graph) {
  59         #pragma omp parallel for
  60         for (size_t j = 0; j < i.size(); j++) {
  61             for (size_t k1 {0}; k1 < i[j].size(); k1++) {
  62                 for (size_t k2 {k1}; k2 < i[j].size(); k2++) {
  63                     if ((frame[j][k1] - frame[j][k2] - refTable[k1][k2]).norm() <= static_cast< float >(eps)) {
  64                         if (k1 != k2) {
  65                             i[j][k1][k2].num++;
  66                             i[j][k2][k1].num++;
  67                         } else {
  68                             i[j][k1][k2].num++;
  69                         }
  70                     }
  71                 }
  72             }
  73         }
  74         #pragma omp barrier
  75     }
  76     // считаем число эффективных обновлений
  77     if (graph.size() > 0) {
  78         updateCount++;
  79     }
  80     // в случае, если число обновлений равно размеру окна, выделяем на нём домены
  81     if (updateCount == window) {
  82         getDomains();
  83         print(frameNumber);
  84     }
  85 }
  86
  87 // инициализация матриц соотношений в структуре
  88 void domainType::setGraph(std::vector< std::vector< node > > &smallGraph) {
  89     smallGraph.resize(0);
  90     smallGraph.resize(refTable.size());
  91     for (size_t i {0}; i < refTable.size(); i++) {
  92         smallGraph[i].resize(0);
  93         smallGraph[i].resize(refTable.size());
  94         for (size_t j {0}; j < refTable.size(); j++) {
  95             smallGraph[i][j].num = 0;
  96             smallGraph[i][j].check = false;
  97         }
  98     }
  99 }
 100
 101 // "зануление" элементов матриц вошедших в домен
 102 void domainType::deleteDomainFromGraph(const std::vector< size_t > &domain) {
 103     for (auto &i : graph.front()) {
 104         for (const auto &j : domain) {
 105             for (size_t k {0}; k < i.size(); k++) {
 106                 i[j][k].check = false;
 107                 i[k][j].check = false;
 108             }
 109         }
 110     }
 111 }
 112
 113 // подсчёт размеров всех потенциально возможных доменов и проверка на наличие домена для выделения
 114  bool domainType::searchDomainSizes() {
 115     bool flag {false};
 116     for (size_t i {0}; i < graph.front().size(); i++) {
 117         domsizes[i].resize(0); // необходимо переопределять память
 118         domsizes[i].resize(graph.front()[i].size(), 0);
 119         for (size_t j {0}; j < graph.front()[i].size(); j++) {
 120             for (const auto &k : graph.front()[i][j]) {
 121                 if (k.check) {
 122                     domsizes[i][j]++;
 123                 }
 124             }
 125             if ((!flag) && (domsizes[i][j] >= dms)) {
 126                 flag = true;
 127             }
 128         }
 129     }
 130     return flag;
 131 }
 132
 133 // функция выделения доменов
 134 void domainType::getDomains() {
 135     // проверка элементов матриц на создание на их основе доменов
 136     for (auto &i : graph.front()) {
 137         for (auto &j : i) {         //  -> матрица соотношений в структуре всё на всё
 138             for (auto &k : j) {     //  -> матрица соотношений в структуре всё на всё
 139                 if (k.num >= window * dlt) {
 140                     k.check = true;
 141                 }
 142             }
 143         }
 144     }
 145     domains.resize(0);
 146     // итеративное выделение доменов одним из двух(пока что) способов
 147     // в случае если цикл запустился, то гарантированно имеется домен для выделения
 148     while (searchDomainSizes()) {
 149         size_t t1 {0}, t2 {0};
 150         for (size_t i {0}; i < domsizes.size(); i++) {
 151             for (size_t j {0}; j < domsizes[i].size(); j++) {
 152                 if ((dsa == 0) && (domsizes[i][j] > domsizes[t1][t2])) {
 153                     // из-за физических ограничений по памяти в общем случае нельзя хранить всю подноготную данных на
 154                     // основе которых было полученно "вхождение" атомов в домен, потому мы только знаем, что размер
 155                     // совпадает, а сказать который из двух стабильнее нельзя, либо нужно придумать что-то новое
 156                     t1 = i;
 157                     t2 = j;
 158                 }
 159                 if ((dsa == 2) && ((domsizes[i][j] < domsizes[t1][t2]) || ((domsizes[i][j] >= dms) && (domsizes[t1][t2] < dms)))) {
 160                     t1 = i;
 161                     t2 = j;
 162                 }
 163             }
 164         }
 165         // выделяем найдённый домен
 166         domains.resize(domains.size() + 1);
 167         for (size_t i {0}; i < graph.front()[t1][t2].size(); i++) {
 168             if (graph.front()[t1][t2][i].check) {
 169                 domains.back().push_back(i);
 170             }
 171         }
 172         //выход из поиска доменов при данном режиме обработки
 173         if (dsa == 1) {
 174             break;
 175         }
 176         // удаляем его из матриц
 177         deleteDomainFromGraph(domains.back());
 178     }
 179     // удаляем рассматриваемое окно из вектора
 180     graph.erase(graph.begin());
 181     // смещаем счётчик обновлений для следующего окна, std::max для случая если окна не пересекаются
 182     updateCount = std::max(0, window - ts);
 183 }
 184
 185 // функция печати ndx и selectionList файлов
 186 void domainType::print(int currentFrame) {
 187     if (domains.size() == 0) {
 188         return;
 189     }
 190     FILE *ndxFile {std::fopen(outPut.c_str(), "a+")}, *slFile {std::fopen(("selectionList-" + outPut.substr(0, outPut.size() - 4)).c_str(), "a+")};
 191     short int writeCount {0};
 192     for (size_t i {0}; i < domains.size(); i++) {
 193         // domain - стартовая позиция в фреймах - номер домена - минимальный размер домена -
 194         // константа тепловых колебаний (отклонения) - константа входимости (отклонения)
 195         std::fprintf(ndxFile, "[domain-stPos-%06d-num-%03lu-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", currentFrame - window + 1, i + 1, dms, eps, dlt);
 196         std::fprintf(slFile, "group %cdomain-stPos-%06d-num-%03lu-dms-%03d-epsi-%04.3f-delta-%04.3f%c;\n", '"', currentFrame - window + 1, i + 1, dms, eps, dlt, '"');
 197         writeCount = 0;
 198         std::printf("\n");
 199         for (size_t j {0}; j < domains[i].size(); j++) {
 200             writeCount++;
 201             if (writeCount > 20) {
 202                 writeCount -= 20;
 203                 std::fprintf(ndxFile, "\n");
 204             }
 205             std::fprintf(ndxFile, "%5lu ", structIndex[domains[i][j]] + 1);
 206             std::printf("%5lu ", structIndex[domains[i][j]] + 1);
 207         }
 208         std::fprintf(ndxFile,"\n\n");
 209         std::printf("\n\n");
 210     }
 211     std::fprintf(ndxFile,"\n");
 212     std::fclose(ndxFile);
 213     std::fclose(slFile);
 214 }