using gmx::RVec;
+// вычисление модуля расстояния между двумя точками, с учётом геометрического преобразования
double F (double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
return sqrt( pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2) +
pow(aiy - biy_plus_y * (cos(A) * cos(C) + sin(A) * sin(B) * sin(C)) + biz_plus_z * (cos(C) * sin(A) - cos(A) * sin(B) * sin(C)) - cos(B) * sin(C) * bix_plus_x, 2) +
pow(aiz + sin(B) * bix_plus_x - cos(A) * cos(B) * biz_plus_z - cos(B) * sin(A) * biy_plus_y, 2) );
}
+// вычисление функции F и её частных производных
void searchF0xyzabc (double &F, double &Fx, double &Fy, double &Fz, double &Fa, double &Fb, double &Fc, double aix, double aiy, double aiz, double bix_plus_x, double biy_plus_y, double biz_plus_z, double A, double B, double C) {
double t01, t02, t03, t04, t05, t06, t07;
t01 = pow(aix + biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) - biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - cos(B) * cos(C) * bix_plus_x, 2);
2 * (biz_plus_z * (sin(A) * sin(C) + cos(A) * cos(C) * sin(B)) - biy_plus_y * (cos(A) * sin(C) - cos(C) * sin(A) * sin(B)) + cos(B) * cos(C) * bix_plus_x) * t07) / (2 * t04);
}
+// применения геометрического преобразования: смещение на вектор (x, y, z) и повороты на градусы A, B, C относительно осей (x, y, z)
void ApplyFit (std::vector< RVec > &b, double x, double y, double z, double A, double B, double C) {
double t0 = 0, t1 = 0, t2 = 0;
for (unsigned int i = 0; i < b.size(); i++) {
}
}
+// подсчёт геометрических центров множеств a и b на основе пар pairs
void CalcMid (std::vector< RVec > a, std::vector< RVec > b, RVec &midA, RVec &midB, std::vector< std::pair< unsigned int, unsigned int > > pairs) {
midA[0] = 0;
midA[1] = 0;
midB[2] /= pairs.size();
}
+// функция фитирования фрейма b на фрейм a на основе пар pairs с "точностью" FtCnst
void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::pair< unsigned int, unsigned int > > pairs, double FtCnst) {
double f1 = 0, f2 = 0, fx = 0, fy = 0, fz = 0, fa = 0, fb = 0, fc = 0, l = 1;
RVec ma, mb;
CalcMid(a, b, ma, mb, pairs);
ma -= mb;
double x = static_cast< double >(ma[0]), y = static_cast< double >(ma[1]), z = static_cast< double >(ma[2]), A = 0, B = 0, C = 0;
+ // поиск стартового отклонения множеств
for (unsigned int i = 0; i < pairs.size(); i++) {
f1 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
static_cast< double >(b[pairs[i].second][0]) + x, static_cast< double >(b[pairs[i].second][1]) + y, static_cast< double >(b[pairs[i].second][2]) + z, A, B, C);
if (f1 == 0) {
return;
} else {
+ // поиск оптимального геометрического преобразования методом градиентного спуска
while (f1 - f2 < 0 || f1 - f2 > FtCnst) {
f1 = 0; fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0; l = 1;
+ // поиск производных и отклонения при заданных параметрах
for (unsigned int i = 0; i < pairs.size(); i++) {
searchF0xyzabc(f1, fx, fy, fz, fa, fb, fc,
static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
}
while (f2 != f1) {
f2 = 0;
+ // поиск отклонения при новых параметрах
for (unsigned int i = 0; i < pairs.size(); i++) {
f2 += F(static_cast< double >(a[pairs[i].first][0]), static_cast< double >(a[pairs[i].first][1]), static_cast< double >(a[pairs[i].first][2]),
static_cast< double >(b[pairs[i].second][0]) + (x - l * fx), static_cast< double >(b[pairs[i].second][1]) + (y - l * fy),
fx = 0; fy = 0; fz = 0; fa = 0; fb = 0; fc = 0;
break;
} else {
- l *= 0.50;
+ // по факту существует минимальное значение переменной типа double
+ // в случае его достижения дважды останавливаем цикл т.к. упёрлись в предел допустимой точности
+ // таким образом гарантируем выход из цикла
if (l == DBL_MIN) { //DBL_TRUE_MIN
break;
}
+ l *= 0.50;
}
}
}
class domainsType {
public:
+ // деструктор класса
~domainsType() {}
+ // конструктор класса для инициализации
domainsType() {}
+ // конструктор класса для полноценной инициализации данных
domainsType(std::vector< int > index, const std::vector< RVec > reference,
const int windowSize, const int domainMinimumSize,
const int domainSearchAlgorythm, const int timeStepBetweenWindowStarts,
}
// set numeric values to "-1" / string value to "" - if you want default settings
+ // функция заполнения необходимыми данными для вычисления структурных доменов
void setDefaults(std::vector< int > index, const std::vector< RVec > reference,
const int windowSize, const int domainMinimumSize,
const int domainSearchAlgorythm, const int timeStepBetweenWindows,
domsizes.resize(sliceNum);
}
+ // фукнция обновления данных для выделения структурных доменов
void update(const std::vector< std::vector< RVec > > frame, const int frameNumber) {
+ // громакс считает с нуля, проверял; инициализируем новое "окно"
if (frameNumber % ts == 0) {
graph.resize(graph.size() + 1);
graph.back().resize(structIndex.size() - static_cast< unsigned long >(dms) + 1);
setGraph(graph.back()[i], ref);
}
}
+ // заполняем структуру данных для структурных доменов
for (unsigned int i = 0; i < graph.size(); i++) {
#pragma omp parallel for
for (unsigned long j = 0; j < graph.front().size(); j++) {
}
#pragma omp barrier
}
+ // считаем число эффективных обновлений
if (graph.size() > 0) {
updateCount++;
}
+ // в случае, если число обновлений равно размеру окна, выделяем на нём домены
if (updateCount == window) {
getDomains();
print(frameNumber);
}
private:
+ // узел матрицы связности: число вхождений в один домен, вектор между элементами в референсной структуре, флаг для проверки
+ // последние два пункта технически можно перевычеслять при каждом обращении, что сократит размер ячейки с 20 байт до 4
+ // но замедлит скорость работы программы, возможно стоит проверить на больших структурах
struct triple {
int num;
RVec radius;
bool check;
};
+ // референстная структура
std::vector< RVec > ref; // must be presented
+ // индекс структуры
std::vector< int > structIndex; // must be presented
+ // рабочие окна -> рассматриваемые основные последовательности -> матрица соотношений в структуре всё на всё
std::vector< std::vector< std::vector< std::vector< triple > > > > graph;
+ // списки доменов
std::vector< std::vector< unsigned int > > domains;
+ // размеры доменов для каждой основной последовательности
std::vector< std::vector< int > > domsizes;
+ // размер рассматриваемого окна в последовательных фреймах траектории
int window = 1000; // selectable
+ // счётчик эффективных последовательных обновлений структуры graph
int updateCount = 0;
+ // минимальный рассматриваемый размер домена - осмысленно брать 4+
int dms = 4; // selectable
+ // какой из алгоритмов выделения использовать 0 - максимальный домен / 1 - минимальный домен
int dsa = 0; // selectable
+ // константа допустимых тепловых колебаний относительно центра домена, находящегося в конкретном атоме структуры
double eps = 0.2; // selectable
+ // константа допустимой входимости атома в выделяемый домен, призвана сгладить резкую границу константы eps
double dlt = 0.95; // selectable
+ // шаг между началами рассматриваемых окон на траектории в фреймах
int ts = window / 10; // selectable
+ // название выходного файла для записи доменов, так же происходит создание selectionList для выделенных доменов
std::string outPut = "default.ndx"; // selectable
+ // инициализация матриц соотношений в структуре
void setGraph(std::vector< std::vector< triple > > &smallGraph, const std::vector< RVec > reference) {
smallGraph.resize(reference.size());
for (unsigned i = 0; i < reference.size(); i++) {
}
}
+ // "зануление" элементов матриц вошедших в домен
void deleteDomainFromGraph(std::vector< unsigned int > domain) {
for (unsigned int i = 0; i < graph.front().size(); i++) {
for (unsigned int j = 0; j < domain.size(); j++) {
}
}
+ // подсчёт размеров всех потенциально возможных доменов и проверка на наличие домена для выделения
bool searchDomainSizes() {
bool flag = false;
for (unsigned int i = 0; i < graph.front().size(); i++) {
return flag;
}
+ // функция выделения доменов
void getDomains() {
+ // проверка элементов матриц на создание на их основе доменов
for (unsigned int i = 0; i < graph.front().size(); i++) {
for (unsigned int j = 0; j < graph.front()[i].size(); j++) {
for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
}
}
domains.resize(0);
+ // итеративное выделение доменов одним из двух(пока что) способов
+ // в случае если цикл запустился, то гарантированно имеется домен для выделения
while (searchDomainSizes()) {
unsigned int t1 = 0, t2 = 0;
for (unsigned int i = 0; i < domsizes.size(); i++) {
for (unsigned int j = 0; j < domsizes[i].size(); j++) {
if ((dsa == 0) && (domsizes[i][j] > domsizes[t1][t2])) {
- // подумать как понять какой домен лучше при одинаковом размере
+ // из-за физических ограничений по памяти в общем случае нельзя хранить всю подноготную данных на
+ // основе которых было полученно "вхождение" атомов в домен, потому мы только знаем, что размер
+ // совпадает, а сказать который из двух стабильнее нельзя, либо нужно придумать что-то новое
t1 = i;
t2 = j;
}
}
}
}
- if (domsizes[t1][t2] >= dms) {
- domains.resize(domains.size() + 1);
- for (unsigned int i = 0; i < graph.front()[t1][t2].size(); i++) {
- if (graph.front()[t1][t2][i].check) {
- domains.back().push_back(i);
- }
+ // выделяем найдённый домен
+ domains.resize(domains.size() + 1);
+ for (unsigned int i = 0; i < graph.front()[t1][t2].size(); i++) {
+ if (graph.front()[t1][t2][i].check) {
+ domains.back().push_back(i);
}
- deleteDomainFromGraph(domains.back());
- } else {
- break;
}
+ // удаляем его из матриц
+ deleteDomainFromGraph(domains.back());
}
+ // удаляем рассматриваемое окно из вектора
graph.erase(graph.begin());
+ // смещаем счётчик обновлений для следующего окна, std::max для случая если окна не пересекаются
updateCount = std::max(0, window - ts);
}
+ // функция печати ndx и selectionList файлов
void print(int currentFrame) {
if (domains.size() == 0) {
return;
}
FILE *ndxFile, *slFile;
ndxFile = std::fopen(outPut.c_str(), "a+");
- slFile = std::fopen(("SelectionList-" + outPut.substr(0, outPut.size() - 4)).c_str(), "a+");
+ slFile = std::fopen(("selectionList-" + outPut.substr(0, outPut.size() - 4)).c_str(), "a+");
int writeCount;
for (unsigned int i = 0; i < domains.size(); i++) {
- std::fprintf(ndxFile, "[domain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", currentFrame - window, i + 1, dms, eps, dlt);
- std::fprintf(slFile, "group %cdomain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f%c;\n", '"', currentFrame - window, i + 1, dms, eps, dlt, '"');
- writeCount = 0;
- for (unsigned int j = 0; j < domains[i].size(); j++) {
- writeCount++;
- if (writeCount > 20) {
- writeCount -= 20;
- std::fprintf(ndxFile, "\n");
- }
- std::fprintf(ndxFile, "%5d ", structIndex[domains[i][j]] + 1);
- std::printf("%5d ", structIndex[domains[i][j]] + 1);
+ // domain - стартовая позиция в фреймах - номер домена - минимальный размер домена -
+ // константа тепловых колебаний (отклонения) - константа входимости (отклонения)
+ std::fprintf(ndxFile, "[domain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", currentFrame - window, i + 1, dms, eps, dlt);
+ std::fprintf(slFile, "group %cdomain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f%c;\n", '"', currentFrame - window, i + 1, dms, eps, dlt, '"');
+ writeCount = 0;
+ for (unsigned int j = 0; j < domains[i].size(); j++) {
+ writeCount++;
+ if (writeCount > 20) {
+ writeCount -= 20;
+ std::fprintf(ndxFile, "\n");
}
- std::fprintf(ndxFile,"\n\n");
- std::printf("\n\n");
+ std::fprintf(ndxFile, "%5d ", structIndex[domains[i][j]] + 1);
+ std::printf("%5d ", structIndex[domains[i][j]] + 1);
+ }
+ std::fprintf(ndxFile,"\n\n");
+ std::printf("\n\n");
}
std::fprintf(ndxFile,"\n");
std::fclose(ndxFile);
Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top)
{
+ // считывание индекса
index.resize(0);
for (ArrayRef< const int >::iterator ai = selec.atomIndices().begin(); (ai < selec.atomIndices().end()); ai++) {
index.push_back(*ai);
}
+ // считывание референсной структуры
reference.resize(0);
if (top.hasFullTopology()) {
for (unsigned int i = 0; i < index.size(); i++) {
reference.push_back(top.x().at(index[i]));
}
}
+
+ // вычисление числа основных (от основаб основание) последовательностей
+ // в идеале нужно брать все сочетания, но ограничиваемся последовательными наборами
bone = static_cast< unsigned int >(index.size() - static_cast< unsigned int >(domain_min_size) + 1);
+ // создание пар для фитирования структур
fitPairs.resize(bone);
for (unsigned int i = 0; i < bone; i++) {
fitPairs[i].resize(0);
trajectory.resize(index.size());
+ // задание необходимых параметров для вычисления доменов
testSubject.setDefaults(index, reference, window, domain_min_size, DomainSearchingAlgorythm, twStep, bone, epsi, delta, fnNdx_);
}
void
Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnalysisModuleData *pdata)
{
+ // считывания текущего фрейма траектории
for (unsigned int i = 0; i < index.size(); i++) {
trajectory[i] = fr.x[index[i]];
}
+ // создание копий фрейма для каждой основной последовательности
std::vector< std::vector< RVec > > tsTemp;
tsTemp.resize(0);
tsTemp.resize(bone, trajectory);
- #pragma omp parallel for ordered schedule(dynamic) shared(fitPairs) firstprivate(reference)
+ // фитирование каждой копии
+ #pragma omp parallel for ordered schedule(dynamic) firstprivate(reference)
for (unsigned int i = 0; i < bone; i++) {
MyFitNew(reference, tsTemp[i], fitPairs[i], 0);
}
#pragma omp barrier
+ // (обновление/потенциальное выделение) доменов
testSubject.update(tsTemp, frnr);
}
biod.pnpi.spb.ru / alexxy / gromacs-domains.git / commitdiff