/usr/include/blasr/algorithms/alignment/sdp/NonoverlappingSparseDynamicProgramming.h is in libblasr-dev 0~20151014+gitbe5d1bf-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 | #ifndef NONOVERLAPPING_SPARSE_DYNAMIC_PROGRAMMING_H_
#define NONOVERLAPPING_SPARSE_DYNAMIC_PROGRAMMING_H_
/*
*
* Compute the subset of fragmentSet of largest weight that is not overlapping.
*
*/
#include "FragmentSort.h"
#include "SDPSet.h"
#include "SDPFragment.h"
#include "SDPColumn.h"
#include "../AlignmentUtils.h"
#include "../../../datastructures/alignment/Alignment.h"
template<typename T_Fragment>
int SDPHeaviestSubsequence(int queryLength,
vector<T_Fragment> &fragmentSet,
int indel, int match,
vector<int> &maxFragmentChain,
AlignmentType alignType=Global) {
maxFragmentChain.clear();
if (fragmentSet.size() < 1)
return 0;
std::sort(fragmentSet.begin(), fragmentSet.end(), LexicographicFragmentSort<T_Fragment>());
SDPSet<SDPColumn> activeSet;
int sweepX;
int trailRow;
int fSweep;
int fi;
for (fi = 0; fi < fragmentSet.size(); fi++) {
fragmentSet[fi].index = fi;
}
sweepX = fragmentSet[0].x;
fSweep = 0;
int maxChainLength = 0;
int maxChainFragment = -1;
int minFragmentCost, minFragmentIndex;
minFragmentCost = INF_INT;
minFragmentIndex = -1;
for (sweepX = 0; sweepX < queryLength; sweepX++) {
//
// For each fragment ending in sweepX, attempt to compute its cost
// by looking for the best previous match.
//
int fSweepStart = fSweep;
while (fSweep < fragmentSet.size() and
fragmentSet[fSweep].x == sweepX) {
//
// Compute the cost of every fragment in the sweep.
//
int cp = INF_INT, cl = INF_INT, ca = INF_INT;
SDPColumn curCol, predCol;
curCol.col = fragmentSet[fSweep].y;
//
// Search preceeding fragments.
//
//
// Compute the cost of fragment_f
int foundPrev = 0;
T_Fragment fragStartPoint;
fragStartPoint.x = fragmentSet[fSweep].x - fragmentSet[fSweep].GetLength();
fragStartPoint.y = fragmentSet[fSweep].y - fragmentSet[fSweep].GetLength();
T_Fragment predFrag;
int predFragScore;
if (activeSet.Predecessor(fragStartPoint, predFrag)) {
//
// predCol points to the fragment with greatest value less than curCol.
//
predFragScore = fragmentSet[predCol.optFragment].cost +
abs((fragmentSet[fSweep].x - fragmentSet[fSweep].y)
- (fragmentSet[predCol.optFragment].x - fragmentSet[predCol.optFragment].y)) * indel;
foundPrev = 1;
}
//
// Now assign the score for the fragment. When doing a local alignment
// start a new chain if the score is not good.
if (foundPrev and
(alignType == Global or
(alignType == Local and predFragScore < 0))) {
fragmentSet[fSweep].chainPrev = predCol.optFragment;
fragmentSet[fSweep].cost = predFragScore - fragmentSet[fSweep].weight;
fragmentSet[fSweep].chainLength = fragmentSet[fragmentSet[fSweep].chainPrev].chainLength + 1;
}
else {
fragmentSet[fSweep].chainPrev = -1;
fragmentSet[fSweep].cost = fragmentSet[fSweep].GetLength() * match - fragmentSet[fSweep].weight;
fragmentSet[fSweep].chainLength = 1;
}
if (minFragmentCost > fragmentSet[fSweep].cost) {
minFragmentCost = fragmentSet[fSweep].cost;
minFragmentIndex = fSweep;
// maxChainLength = fragmentSet[fSweep].chainLength;
}
if (fragmentSet[fSweep].chainLength > maxChainLength) {
maxChainLength = fragmentSet[fSweep].chainLength;
maxChainFragment = fSweep;
}
fSweep++;
}
//
// Each column must contain the highest scoring hit in that column.
//
fSweep = fSweepStart;
while (fSweep < fragmentSet.size() and
fragmentSet[fSweep].x == sweepX) {
//
// These elements are removed from the sweep set since they are done being processed.
// If they are the lowest cost in the value, update activeSet
//
SDPColumn col;
int storeCol = 0;
col.col = fragmentSet[fSweep].y;
if (activeSet.Member(col)) {
if (fragmentSet[col.optFragment].cost < fragmentSet[fSweep].cost) {
storeCol = 1;
}
}
else {
storeCol = 1;
}
if (storeCol) {
col.col = fragmentSet[fSweep].y;
col.optFragment = fSweep;
//
// Insert new column or replace col with a more optimal one.
//
activeSet.Insert(col);
//
// The invariant structure of the activeSet is that
// after inserting a fragment of score S at column col,
// the score of all columns greater than 'col' in col set
// must be less than col.
//
// To preserve this invariant, when an element is inserted
// at 'col', look to columns greater. As long as any columns
// have scores that are greater than col, remove them.
// Once a column col_next has been found that has a score less than S
// by the structure of the loop invariant, all columns greater than col_next
// are guaranteed to have lower score than S, so we can continue searching
// through this loop.
//
// Since fragments are processed at most once, this remains O(M).
SDPColumn successorCol = col;
while (activeSet.Successor(col, successorCol) and
fragmentSet[successorCol.optFragment].cost > fragmentSet[fSweep].cost) {
activeSet.Delete(successorCol);
}
}
//
// Now remove this fragment, it is at the end of the sweep line.
//
int deleted;
deleted = activeSet.Delete(fragmentSet[fSweep]);
assert(deleted);
++fSweep;
}
}
if (alignType == Local) {
maxChainFragment = minFragmentIndex;
}
while (maxChainFragment != -1) {
maxFragmentChain.push_back(maxChainFragment);
maxChainFragment = fragmentSet[maxChainFragment].chainPrev;
}
std::reverse(maxFragmentChain.begin(), maxFragmentChain.end());
return maxFragmentChain.size();
}
#endif
|