src/domaintype.cpp

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