From: Anatoly Date: Wed, 29 Jan 2020 15:23:53 +0000 (+0300) Subject: class debugging X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=3aeb99cd02c74482923332caab9fc8887600c3a4;p=alexxy%2Fgromacs-domains.git class debugging --- diff --git a/src/domains.cpp b/src/domains.cpp index eec078b..038460b 100644 --- a/src/domains.cpp +++ b/src/domains.cpp @@ -63,8 +63,6 @@ using namespace gmx; 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) + @@ -185,13 +183,38 @@ class domainsType { 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++) { @@ -199,10 +222,11 @@ class domainsType { } } 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 { @@ -213,109 +237,15 @@ class domainsType { } } } + #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: @@ -324,15 +254,19 @@ class domainsType { 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()); @@ -357,15 +291,99 @@ class domainsType { } } - 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); } }; @@ -509,8 +527,8 @@ void print_domains(std::vector< std::vector< unsigned int > > pd_domains, std::v 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++; @@ -578,15 +596,16 @@ class Domains : public TrajectoryAnalysisModule 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. }; @@ -624,13 +643,21 @@ Domains::initOptions(IOptionsContainer *options, .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); @@ -663,19 +690,18 @@ Domains::initAnalysis(const TrajectoryAnalysisSettings &settings, 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++) { @@ -683,8 +709,9 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnal } } - 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; @@ -706,16 +733,16 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnal 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); @@ -723,21 +750,18 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnal 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