X-Git-Url: https://ruin.nu/git/?a=blobdiff_plain;f=src%2Fgenealgorithms.cpp;h=3fb2aaac0c7f180bd73edde5918731c7bcc8c8e0;hb=43b1052d74a207fe667b75e41d0c70c1189c2cc8;hp=d9d86a3dd8fcb651ce497180b062fa3ddc562612;hpb=6724d50f2d36dccd383d95d36f7fe328b2692bc0;p=germs.git diff --git a/src/genealgorithms.cpp b/src/genealgorithms.cpp index d9d86a3..3fb2aaa 100644 --- a/src/genealgorithms.cpp +++ b/src/genealgorithms.cpp @@ -20,9 +20,11 @@ #include "genealgorithms.h" #include "geneorder.h" +#include "componenttree.h" #include #include +#include #include #include using namespace std; @@ -65,54 +67,21 @@ struct FindP{ }; -std::vector findIntervalsAtPoints(const vector& intervals){ - vector points; - for (size_t p = 1; p <= intervals.size(); ++p){ - size_t f = 0; - size_t s = 0; - bool found = false; - size_t n = 0; - for (vector::const_iterator i = intervals.begin(); i != intervals.end(); ++i, ++n){ - if (i->first == p){ - if (!found){ - f = n; - found = true; - }else{ - s = n; - break; - } - } - if (i->second == p){ - if (!found){ - f = n; - found = true; - }else{ - s = n; - break; - } - } - } - points.push_back(Interval(f,s)); - } - return points; - -} - -int countCycles(const GeneOrder& go){ - int cycles = 0; +size_t countCycles(const GeneOrder& go){ + size_t cycles = 0; set marked; vector intervals = findIntervals(go); vector points = findIntervalsAtPoints(intervals); for (size_t p = 1; p < go.size(); ++p){ if (marked.find(p) != marked.end()) continue; - Interval i = intervals[points[p-1].first]; + Interval i = intervals[points[p].first]; while (marked.find(p) == marked.end()){ marked.insert(p); - if (i == intervals[points[p-1].first]) - i = intervals[points[p-1].second]; + if (i == intervals[points[p].first]) + i = intervals[points[p].second]; else - i = intervals[points[p-1].first]; + i = intervals[points[p].first]; if (p == i.first) p = i.second; @@ -124,41 +93,154 @@ int countCycles(const GeneOrder& go){ return cycles; } +size_t inversionDistance(const GeneOrder& go){ + size_t distance = go.size() - 1; + distance -= countCycles(go); + + ComponentTree t(findComponents(go)); + t.makeUnoriented(); + size_t leaves = t.countLeaves(); + distance += leaves; + if (leaves % 2 != 0){ + size_t sb = t.shortBranches(); + if (sb == 0) + distance += 1; + } + + return distance; +} + +int sign(Gene g){ + if (g > 0) + return 1; + if (g < 0) + return -1; + return 0; +} + +struct Abs{ + Gene operator()(Gene x) const{ + return abs(x); + } +}; std::vector findComponents(const GeneOrder& go){ - return vector(); + vector components; + vector os(go.size()-1); + for (size_t i = 0; i < os.size(); ++i) + os[i] = (go[i]*go[i+1] > 0 ? sign(go[i]) : 0); + stack Mdir; + Mdir.push(go.size()-1); + stack Mrev; + Mrev.push(0); + stack Sdir; + Sdir.push(0); + stack Srev; + Srev.push(0); + vector dir; + dir.push_back(go.size()-1); + vector rev; + rev.push_back(0); + size_t s; + vector p(go.list()); + transform(p.begin(),p.end(),p.begin(),Abs()); + for (size_t i = 1; i < go.size(); ++i){ + //Directed + if (p[i-1] > p[i]) + Mdir.push(p[i-1]); + else while (Mdir.top() < p[i]) + Mdir.pop(); + dir.push_back(Mdir.top()); + + s = Sdir.top(); + while(p[Sdir.top()] > p[i] || dir[Sdir.top()] < p[i]){ + Sdir.pop(); + os[Sdir.top()] = (os[Sdir.top()] == os[s] ? os[s] : 0); + s = Sdir.top(); + } + if (go[i] > 0 && dir[i] == dir[s] && static_cast(i - s) == p[i] - p[s]) + components.push_back(Component(p[s],p[i],(s+1 == i ? 0 : os[s]),s,i)); + + //Reverse + if (p[i-1] < p[i]) + Mrev.push(p[i-1]); + else while (Mrev.top() > p[i]) + Mrev.pop(); + rev.push_back(Mrev.top()); + + s = Srev.top(); + while((p[s] < p[i] || rev[s] > p[i]) && s > 0){ + Srev.pop(); + os[Srev.top()] *= (os[Srev.top()] == os[s] ? 1 : 0); + s = Srev.top(); + } + if (go[i] < 0 && rev[i] == rev[s] && static_cast(i - s) == p[s] - p[i]) + components.push_back(Component(-p[s],-p[i],(s+1 == i ? 0 : os[s]),s,i)); + + //Update stacks + if (go[i] > 0) + Sdir.push(i); + else + Srev.push(i); + } + return components; } +int sign2(Gene g){ + if (g < 0) + return -1; + return 1; +} /** - * TODO: Think of a better than O(n^2) implementation - * Possibly move it to GeneOrder too + * */ std::vector findIntervals(const GeneOrder& go){ - vector intervals; - for (size_t i = 0; i < go.size() - 1; ++i){ - size_t f = 0; - size_t s = 0; - bool found = false; - size_t n = 0; - for (GeneOrder::iterator g = go.begin(); g != go.end(); ++g, ++n){ - if (static_cast(abs(*g)) == i){ - f = n; - if (*g >= 0) - ++f; - if (found) - break; - found = true; + const size_t max = go.size(); + vector intervals(go.size()-1,Interval(max,max,false)); + size_t n = 0; + const GeneOrder::iterator end = go.end(); + for (GeneOrder::iterator g = go.begin(); g != end; ++g, ++n){ + size_t i = abs(*g); + if (i < max - 1){ + Interval& curr = intervals[i]; + curr.first = n + (*g >= 0 ? 1 : 0); + + if (curr.second == max) + curr.oriented = *g < 0; + else + curr.oriented ^= *g < 0; } - if(static_cast(abs(*g)) == i+1){ - s = n; - if (*g < 0) - ++s; - if (found) - break; - found = true; + if (i > 0){ + Interval& prev = intervals[i-1]; + prev.second = n + (*g < 0 ? 1 : 0); + + if (prev.first == max) + prev.oriented = *g < 0; + else + prev.oriented ^= *g < 0; } - } - intervals.push_back(Interval(f,s)); + } return intervals; } +/** + * + */ +std::vector findIntervalsAtPoints(const vector& intervals){ + size_t max = intervals.size()+1; + vector points(max,Interval(max,max,false)); + size_t n = 0; + for (vector::const_iterator i = intervals.begin(); i != intervals.end(); ++i, ++n){ + if (points[i->first].first == max){ + points[i->first].first = n; + }else + points[i->first].second = n; + + if (points[i->second].first == max){ + points[i->second].first = n; + }else + points[i->second].second = n; + } + return points; + +}