src/domaintype.cpp

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