implemented domainType class

diff --git a/src/domains.cpp b/src/domains.cpp
index 44c9545a36fc24aaafa6de5ffc988c1a4edd8db2..eec078bfdcd23b551d4699abf3260129636cee1d 100644 (file)
--- a/src/domains.cpp
+++ b/src/domains.cpp
@@ -178,6 +178,197 @@ void MyFitNew (std::vector< RVec > a, std::vector< RVec > &b, std::vector< std::
     }
 }
 
+class domainsType {
+    public:
+
+        ~domainsType() {}
+
+        domainsType() {}
+
+        domainsType(const unsigned int sliceNum, const std::vector< RVec > reference) {
+            setDefault(sliceNum, reference);
+        }
+
+        void update(const std::vector< std::vector< RVec > > frame, const float epsilon, const int frameNumber) {
+            if (updateCount == window)
+            if ((frameNumber + 1) % ts == 0) {
+                graph.resize(graph.size() + 1);
+                graph.back().resize(graph.front().size());
+                for (unsigned i = 0; i < graph.front().size(); i++) {
+                    setGraph(graph.back()[i], ref);
+                }
+            }
+            for (unsigned int i = 0; i < graph.size(); i++) {
+                for (unsigned int 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 (k1 == k2) {
+                                    graph[i][j][k1][k2].num++;
+                                } else {
+                                    graph[i][j][k1][k2].num++;
+                                    graph[i][j][k2][k1].num++;
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+            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;
+            }
+            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");
+            }
+            std::fprintf(ndxFile,"\n");
+            std::fclose(ndxFile);
+            std::fclose(slFile);
+        }
+
+    private:
+        struct triple {
+            int num;
+            RVec radius;
+            bool check;
+        };
+        std::vector< RVec >                                                     ref;
+        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                                                                     updateCount = 0;
+        int                                                                     dms         = 4;    // selectable
+        int                                                                     dsa         = 0;    // selectable
+        int                                                                     ts          = window / 10;
+
+        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++) {
+                smallGraph[i].resize(reference.size());
+                for (unsigned j = 0; j < reference.size(); j++) {
+                    smallGraph[i][j].num = 0;
+                    smallGraph[i][j].radius = reference[i] - reference[j];
+                    smallGraph[i][j].check = false;
+                }
+            }
+        }
+
+        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++) {
+                    for (unsigned int k = 0; k < graph.front()[i].size(); k++) {
+                        graph.front()[i][domain[j]][k].check = false;
+                        graph.front()[i][k][domain[j]].check = false;
+                    }
+                }
+            }
+        }
+
+        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;
+                    }
+                }
+            }
+            return false;
+        }
+};
+
 struct node {
     short int n;
     RVec r;
@@ -376,11 +567,13 @@ class Domains : public TrajectoryAnalysisModule
 
         std::vector< std::vector< std::vector< std::vector< node > > > >        graph;
 
+        domainsType                                                             testSubject;
+
+
         std::vector< std::vector< unsigned int > >                              domains;
         std::vector< std::vector< int > >                                       domsizes;
 
         std::vector< int >                                                      index;
-        std::vector< int >                                                      numbers;
         std::vector< RVec >                                                     trajectory;
         std::vector< RVec >                                                     reference;
         Selection                                                               selec;
@@ -474,6 +667,9 @@ Domains::initAnalysis(const TrajectoryAnalysisSettings &settings,
     }
 
     trajectory.resize(index.size());
+
+    //testSubject.setDefault(bone, reference);
+    //testSubject.setParams(window, domain_min_size, DomainSearchingAlgorythm, window / tau_jump);
 }
 
 void
@@ -529,6 +725,18 @@ Domains::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, TrajectoryAnal
     }
     #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)
+    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";
 }