1 /**** 2 DIAMOND protein aligner 3 Copyright (C) 2013-2017 Benjamin Buchfink <buchfink@gmail.com> 4 5 This program is free software: you can redistribute it and/or modify 6 it under the terms of the GNU General Public License as published by 7 the Free Software Foundation, either version 3 of the License, or 8 (at your option) any later version. 9 10 This program is distributed in the hope that it will be useful, 11 but WITHOUT ANY WARRANTY; without even the implied warranty of 12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13 GNU General Public License for more details. 14 15 You should have received a copy of the GNU General Public License 16 along with this program. If not, see <http://www.gnu.org/licenses/>. 17 ****/ 18 19 #ifndef INTERVAL_PARTITION_H_ 20 #define INTERVAL_PARTITION_H_ 21 22 #include <assert.h> 23 #include <map> 24 #include <limits.h> 25 #include "interval.h" 26 27 struct IntervalNode 28 { IntervalNodeIntervalNode29 IntervalNode(): 30 count(0), 31 min_score(INT_MAX), 32 max_score(0) 33 {} IntervalNodeIntervalNode34 IntervalNode(int count, int min_score, int max_score): 35 count(count), 36 min_score(min_score), 37 max_score(max_score) 38 { 39 } addIntervalNode40 IntervalNode add(int score, int cap) const 41 { 42 return IntervalNode(count + 1, count < cap ? std::min(min_score, score) : min_score, std::max(max_score, score)); 43 } 44 int count, min_score, max_score; 45 }; 46 47 struct IntervalPartition : protected std::map<int, IntervalNode> 48 { 49 50 struct MaxScore {}; 51 struct MinScore {}; 52 53 struct Iterator 54 { IteratorIntervalPartition::Iterator55 Iterator(const_iterator i, const_iterator j, const IntervalPartition &parent) : 56 i_(i), 57 j_(j), 58 parent_(parent) 59 { 60 } goodIntervalPartition::Iterator61 bool good() const 62 { 63 return i_ != parent_.end(); 64 } 65 Iterator& operator++() 66 { 67 i_ = j_; 68 if(j_ != parent_.end()) 69 ++j_; 70 return *this; 71 } 72 std::pair<interval, IntervalNode> operator*() const 73 { 74 return std::make_pair(interval(i_->first, j_ == parent_.end() ? INT_MAX : j_->first), i_->second); 75 } 76 private: 77 const_iterator i_, j_; 78 const IntervalPartition &parent_; 79 }; 80 IntervalPartitionIntervalPartition81 IntervalPartition(int cap): 82 cap(cap) 83 { 84 (*this)[0] = IntervalNode(); 85 } 86 insertIntervalPartition87 void insert(interval k, int score) 88 { 89 iterator i = lower_bound(k.begin_); 90 if (i == end()) 91 i = std::map<int, IntervalNode>::insert(std::make_pair(k.begin_, IntervalNode())).first; 92 else if (i->first != k.begin_) { 93 //assert(i != std::map<int, IntervalNode>::begin()); 94 i--; 95 i = std::map<int, IntervalNode>::insert(std::make_pair(k.begin_, i->second)).first; 96 } 97 IntervalNode last; 98 while (i != end() && i->first < k.end_) { 99 last = i->second; 100 i->second = i->second.add(score, cap); 101 ++i; 102 } 103 if (i == end() || i->first != k.end_) 104 (*this)[k.end_] = last; 105 } 106 coveredIntervalPartition107 int covered(interval k) const 108 { 109 Iterator i = begin(k.begin_); 110 std::pair<interval, IntervalNode> l; 111 int c = 0; 112 while (i.good() && (l = *i).first.begin_ < k.end_) { 113 if (l.second.count >= cap) 114 c += k.overlap(l.first); 115 ++i; 116 } 117 return c; 118 } 119 coveredIntervalPartition120 int covered(interval k, int max_score, const MaxScore&) const 121 { 122 Iterator i = begin(k.begin_); 123 std::pair<interval, IntervalNode> l; 124 int c = 0; 125 while (i.good() && (l = *i).first.begin_ < k.end_) { 126 if (l.second.max_score >= max_score) 127 c += k.overlap(l.first); 128 ++i; 129 } 130 return c; 131 } 132 coveredIntervalPartition133 int covered(interval k, int min_score, const MinScore&) const 134 { 135 Iterator i = begin(k.begin_); 136 std::pair<interval, IntervalNode> l; 137 int c = 0; 138 while (i.good() && (l = *i).first.begin_ < k.end_) { 139 if (l.second.count >= cap && l.second.min_score >= min_score) 140 c += k.overlap(l.first); 141 ++i; 142 } 143 return c; 144 } 145 min_scoreIntervalPartition146 int min_score(interval k) const 147 { 148 Iterator i = begin(k.begin_); 149 std::pair<interval, IntervalNode> l; 150 int s = INT_MAX; 151 while (i.good() && (l = *i).first.begin_ < k.end_) { 152 if (l.second.count < cap) 153 return 0; 154 s = std::min(s, l.second.min_score); 155 ++i; 156 } 157 return s; 158 } 159 max_scoreIntervalPartition160 int max_score(interval k) const 161 { 162 Iterator i = begin(k.begin_); 163 std::pair<interval, IntervalNode> l; 164 int s = INT_MAX; 165 while (i.good() && (l = *i).first.begin_ < k.end_) { 166 s = std::min(s, l.second.max_score); 167 ++i; 168 } 169 assert(s != INT_MAX); 170 return s; 171 } 172 beginIntervalPartition173 Iterator begin(int p) const 174 { 175 const_iterator i = lower_bound(p), j; 176 if (i == end() || i->first != p) { 177 //assert(i != std::map<int, IntervalNode>::begin()); 178 j = i; 179 i--; 180 } 181 else { 182 j = i; 183 ++j; 184 } 185 return Iterator(i, j, *this); 186 } 187 188 const int cap; 189 190 }; 191 192 #endif