using gmx::RVec;
-const int tau_jump = 10;
-
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) +
domainsType() {}
- domainsType(const unsigned int sliceNum, const std::vector< RVec > reference) {
- setDefault(sliceNum, reference);
+ domainsType(std::vector< int > index, const std::vector< RVec > reference,
+ const int windowSize, const int domainMinimumSize,
+ const int domainSearchAlgorythm, const int timeStepBetweenWindowStarts,
+ const unsigned int sliceNum, const double epsilon, const double delta,
+ const std::string outPutFileName) {
+ setDefaults(index, reference, windowSize, domainMinimumSize, domainSearchAlgorythm, timeStepBetweenWindowStarts, sliceNum, epsilon, delta, outPutFileName);
}
- void update(const std::vector< std::vector< RVec > > frame, const float epsilon, const int frameNumber) {
- if (updateCount == window)
- if ((frameNumber + 1) % ts == 0) {
+ // 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,
+ const unsigned int sliceNum, const double epsilon, const double delta,
+ const std::string outPutFileName) {
+ graph.resize(1);
+ graph.front().resize(sliceNum);
+ for (unsigned i = 0; i < sliceNum; i++) {
+ setGraph(graph.front()[i], reference);
+ }
+ structIndex = index;
+ ref = reference;
+ if (epsilon != -1) { eps = epsilon; }
+ if (delta != -1) { dlt = delta; }
+ if (windowSize != -1) { window = windowSize; }
+ if (domainMinimumSize != -1) { dms = domainMinimumSize; }
+ if (domainSearchAlgorythm != -1) { dsa = domainSearchAlgorythm; }
+ if (timeStepBetweenWindows != -1) { ts = timeStepBetweenWindows; }
+ if (outPutFileName != "") { outPut = outPutFileName; }
+ }
+
+ void update(const std::vector< std::vector< RVec > > frame, const int frameNumber) {
+ if (frameNumber % ts == 0) {
graph.resize(graph.size() + 1);
graph.back().resize(graph.front().size());
for (unsigned i = 0; i < graph.front().size(); i++) {
}
}
for (unsigned int i = 0; i < graph.size(); i++) {
- for (unsigned int j = 0; j < graph.front().size(); j++) {
+ #pragma omp parallel for
+ for (unsigned long j = 0; j < graph.front().size(); j++) {
for (unsigned int k1 = 0; k1 < graph[i][j].size(); k1++) {
for (unsigned int k2 = k1; k2 < graph[i][j].size(); k2++) {
- if ((frame[j][k1] - frame[j][k2] - graph[i][j][k1][k2].radius).norm() <= epsilon) {
+ if ((frame[j][k1] - frame[j][k2] - graph[i][j][k1][k2].radius).norm() <= static_cast< float >(eps)) {
if (k1 == k2) {
graph[i][j][k1][k2].num++;
} else {
}
}
}
+ #pragma omp barrier
}
- updateCount++;
- }
-
- void getDomains(const float delta) {
- if (updateCount != window) {
- return;
- } else {
- //присмотреться к этому моменту
- updateCount -= ts;
- 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++) {
- if (graph.front()[i][j][k].num >= window * delta) {
- graph.front()[i][j][k].check = true;
- }
- }
- }
- }
- domains.resize(0);
- while (checkDomainSizes()) {
- domsizes.resize(graph.front().size());
- for (unsigned int i = 0; i < graph.front().size(); i++) {
- domsizes[i].resize(graph.front()[i].size(), 0);
- for (unsigned int j = 0; j < graph.front()[i].size(); j++) {
- for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
- if (graph.front()[i][j][k].check) {
- domsizes[i][j]++;
- }
- }
- }
- }
- 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 ((dsa == 1) && ((domsizes[i][j] < domsizes[t1][t2]) || ((domsizes[i][j] >= dms) && (domsizes[t1][t2] < dms)))) {
- t1 = i;
- t2 = j;
- }
- }
- }
- 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());
- }
- if (graph.size() > static_cast< unsigned int >(window / ts)) {
- graph.erase(graph.begin());
- }
- }
- }
-
- void setParams(const int windowSize, const int domainMinimumSize, const int domainSearchAlgorythm, const int timeStepBetweenWindows) {
- window = windowSize;
- dms = domainMinimumSize;
- dsa = domainSearchAlgorythm;
- ts = timeStepBetweenWindows;
- }
-
- void setDefault(const unsigned int sliceNum, const std::vector< RVec > reference) {
- graph.resize(1);
- graph.front().resize(sliceNum);
- for (unsigned i = 0; i < sliceNum; i++) {
- setGraph(graph.front()[i], reference);
- }
- ref = reference;
- }
-
- void print(std::vector< int > index, std::string fileName, double epsilon, double delta, int startingPos) {
- if (domains.size() == 0) {
- return;
+ if (graph.size() > 0) {
+ updateCount++;
}
- FILE *ndxFile, *slFile;
- ndxFile = std::fopen(fileName.c_str(), "a+");
- slFile = std::fopen(("SelectionList-" + fileName.substr(0, fileName.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", static_cast< int >(startingPos) * ts, i + 1, dms, epsilon, delta);
- std::fprintf(slFile, "group %cdomain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f%c;\n", '"', static_cast< int >(startingPos) * ts, i + 1, dms, epsilon, delta, '"');
- 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 ", index[domains[i][j]] + 1);
- std::printf("%5d ", index[domains[i][j]] + 1);
- }
- std::fprintf(ndxFile,"\n\n");
- std::printf("\n\n");
+ if (updateCount == window) {
+ getDomains();
+ print(frameNumber);
}
- std::fprintf(ndxFile,"\n");
- std::fclose(ndxFile);
- std::fclose(slFile);
}
private:
RVec radius;
bool check;
};
- std::vector< RVec > ref;
+ 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
+ int window = 1000; // selectable
int updateCount = 0;
- int dms = 4; // selectable
- int dsa = 0; // selectable
- int ts = window / 10;
+ int dms = 4; // selectable
+ int dsa = 0; // selectable
+ double eps = 0.2; // selectable
+ double dlt = 0.95; // selectable
+ int ts = window / 10; // selectable
+ std::string outPut = "default.ndx"; // selectable
void setGraph(std::vector< std::vector< triple > > &smallGraph, const std::vector< RVec > reference) {
smallGraph.resize(reference.size());
}
}
- bool checkDomainSizes() {
- for (unsigned int i = 0; i < domsizes.size(); i++) {
- for (unsigned int j = 0; j < domsizes[i].size(); j++) {
- if (domsizes[i][j] >= dms) {
- return true;
+ bool searchDomainSizes() {
+ // it goes infinite without these 2 lines, curious even if it was resized before
+ // resizing inside vectors with 0 is not enough
+ domsizes.resize(0);
+ domsizes.resize(graph.front().size());
+ bool flag = false;
+ for (unsigned int i = 0; i < graph.front().size(); i++) {
+ domsizes[i].resize(graph.front()[i].size(), 0);
+ for (unsigned int j = 0; j < graph.front()[i].size(); j++) {
+ for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
+ if (graph.front()[i][j][k].check) {
+ domsizes[i][j]++;
+ }
+ if (!flag) {
+ if (domsizes[i][j] >= dms) {
+ flag = true;
+ }
+ }
+ }
+ }
+ }
+ 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++) {
+ if (graph.front()[i][j][k].num >= window * dlt) {
+ graph.front()[i][j][k].check = true;
+ }
+ }
+ }
+ }
+ 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 ((dsa == 1) && ((domsizes[i][j] < domsizes[t1][t2]) || ((domsizes[i][j] >= dms) && (domsizes[t1][t2] < dms)))) {
+ 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);
+ }
}
+ deleteDomainFromGraph(domains.back());
+ } else {
+ break;
}
}
- return false;
+ graph.erase(graph.begin());
+ updateCount = std::max(0, window - ts);
+ }
+
+ 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+");
+ 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);
+ }
+ std::fprintf(ndxFile,"\n\n");
+ std::printf("\n\n");
+ }
+ std::fprintf(ndxFile,"\n");
+ std::fclose(ndxFile);
+ std::fclose(slFile);
}
};
fpNdx2_ = std::fopen(("SelectionList-" + fnNdx_.substr(0, fnNdx_.size() - 4)).c_str(), "a+");
int write_count;
for (unsigned int i1 = 0; i1 < pd_domains.size(); i1++) {
- std::fprintf(fpNdx_, "[domain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", static_cast< int >(st_pos) * window / tau_jump, i1 + 1, dms, epsi, delta);
- std::fprintf(fpNdx2_, "group %cdomain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f%c;\n", '"', static_cast< int >(st_pos) * window / tau_jump, i1 + 1, dms, epsi, delta, '"');
+ std::fprintf(fpNdx_, "[domain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f]\n", static_cast< int >(st_pos) * window / 10, i1 + 1, dms, epsi, delta);
+ std::fprintf(fpNdx2_, "group %cdomain-stPos-%06d-num-%03d-dms-%03d-epsi-%04.3f-delta-%04.3f%c;\n", '"', static_cast< int >(st_pos) * window / 10, i1 + 1, dms, epsi, delta, '"');
write_count = 0;
for (unsigned int j = 0; j < pd_domains[i1].size(); j++) {
write_count++;
std::vector< RVec > reference;
Selection selec;
int frames = 0;
- int window = 1000; // selectable
- int domain_min_size = 4; // selectable
+ int window = -1; // selectable
+ int domain_min_size = -1; // selectable
std::vector< std::vector< std::pair< unsigned int, unsigned int > > > w_rls2;
unsigned int bone;
- double delta = 0.95; // selectable
- double epsi = 0.10; // selectable
+ double delta = -1; // selectable
+ double epsi = -1; // selectable
+ int twStep = -1;
- int DomainSearchingAlgorythm = 0; // selectable
+ int DomainSearchingAlgorythm = -1; // selectable
// Copy and assign disallowed by base.
};
.store(&DomainSearchingAlgorythm)
.description("Domain's Searching Algorythm: 0 == default (from bigger to smaller) | 1 == (from smaller to bigger)"));
// Add option for epsi constant
- options->addOption(DoubleOption("epsilon")
+ options->addOption(DoubleOption("eps")
.store(&epsi)
.description("thermal vibrations' constant"));
// Add option for delta constant
- options->addOption(DoubleOption("delta")
+ options->addOption(DoubleOption("dlt")
.store(&delta)
.description("domain membership probability"));
+ // Add option for domain min size constant
+ options->addOption(gmx::IntegerOption("wSize")
+ .store(&window)
+ .description("flowing window to evaluate domains from"));
+ // Add option for domain min size constant
+ options->addOption(gmx::IntegerOption("twStep")
+ .store(&twStep)
+ .description("time step between windows' starting positions"));
// Control input settings
settings->setFlags(TrajectoryAnalysisSettings::efNoUserPBC);
settings->setFlag(TrajectoryAnalysisSettings::efUseTopX);
for (unsigned int j = 0; j < static_cast< unsigned int >(domain_min_size); j++) {
w_rls2[i].push_back(std::make_pair(j + i, j + i));
}
- make_graph(index.size(), reference, graph.front()[i]);
+ //make_graph(index.size(), reference, graph.front()[i]);
}
trajectory.resize(index.size());
- //testSubject.setDefault(bone, reference);
- //testSubject.setParams(window, domain_min_size, DomainSearchingAlgorythm, window / tau_jump);
+ 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)
{
- if ((frnr + 1) % (window / tau_jump) == 0) {
+ /*if (frnr % twStep == 0) {
graph.resize(graph.size() + 1);
graph.back().resize(bone);
for (unsigned int i = 0; i < bone; i++) {
}
}
- int temp = graph.size() - tau_jump;
+ int temp = graph.size() - window / twStep;
+ // most probably == 0
if (temp > 0) {
domains.resize(0);
std::vector< unsigned int > a;
find_domain_sizes(graph[0], domsizes);
}
graph.erase(graph.begin());
- std::cout << (frnr + 1) / (window / tau_jump) - tau_jump << " window's domains have been evaluated\n";
- print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta, window, (frnr + 1) / (window / tau_jump) - tau_jump); // see function for details | numbers from index
- }
+ std::cout << (frnr + 1) / twStep - (window / twStep)<< " window's domains have been evaluated\n";
+ print_domains(domains, index, fnNdx_, domain_min_size, epsi, delta, window, (frnr + 1) / twStep - (window / twStep)); // see function for details | numbers from index
+ }*/
for (unsigned int i = 0; i < index.size(); i++) {
trajectory[i] = fr.x[index[i]];
}
frames++;
- #pragma omp parallel for ordered schedule(dynamic) shared(w_rls2, graph) firstprivate(trajectory, reference, epsi, frnr, window)
+ /*#pragma omp parallel for ordered schedule(dynamic) shared(w_rls2, graph) firstprivate(trajectory, reference, epsi, frnr, window)
for (unsigned int j = 0; j < bone; j++) {
std::vector < RVec > temp = trajectory;
MyFitNew(reference, temp, w_rls2[j], 0);
update_graph(graph[i][j], temp, epsi);
}
}
- #pragma omp barrier
+ #pragma omp barrier*/
std::vector< std::vector< RVec > > tsTemp;
tsTemp.resize(bone, trajectory);
- #pragma omp parallel for ordered schedule(dynamic) shared(w_rls2, graph) firstprivate(trajectory, reference, epsi, frnr, window)
+ #pragma omp parallel for ordered schedule(dynamic) shared(w_rls2) firstprivate(trajectory, reference)
for (unsigned int i = 0; i < bone; i++) {
MyFitNew(reference, tsTemp[i], w_rls2[i], 0);
}
#pragma omp barrier
- testSubject.update(tsTemp, epsi, frnr);
- testSubject.getDomains(delta);
- testSubject.print(index, fnNdx_, epsi, delta, (frnr + 1) / (window / tau_jump) - tau_jump);
-
- std::cout << "frame: " << frnr << " analyzed\n";
+ testSubject.update(tsTemp, frnr);
+ //std::cout << "frame: " << frnr << " analyzed\n";
}
void
biod.pnpi.spb.ru / alexxy / gromacs-domains.git / commitdiff