Mercurial > repos > rliterman > csp2

comparison CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/aligner/FlatAligner.java @ 68:5028fdace37b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author jpayne
date Tue, 18 Mar 2025 16:23:26 -0400
parents
children
comparison
equal deleted inserted replaced
67:0e9998148a16 68:5028fdace37b
1 package aligner;
2
3
4
5 public final class FlatAligner {
6
7 public FlatAligner(){}
8
9 /**
10 * @param query
11 * @param ref
12 * @param qstart
13 * @param rstart
14 * @param rstop
15 * @param minScore Quit early if score drops below this
16 * @param minRatio Don't return results if max score is less than this fraction of max possible score
17 * @return
18 */
19 public AlignmentResult alignForward(final byte[] query, final byte[] ref, final int rstart, final int rstop, final int minScore,
20 final float minRatio) {
21 // if(ref.length<16300){return ica16.alignForward(query, ref, qstart, rstart, rstop, minScore, minRatio);}
22 final int arrayLength=rstop-rstart+1;
23 // if(array1==null || arrayLength+1>array1.length){
24 // array1=new int[(int)((arrayLength+2)*1.2)];
25 // array2=new int[(int)((arrayLength+2)*1.2)];
26 // }
27
28 //Stack allocated; faster.
29 final int[] array1=new int[arrayLength+1], array2=new int[arrayLength+1];
30
31 final int minPassingScore=(int)(query.length*minRatio*pointsMatch);
32 final int minPassingScore3=minPassingScore-query.length*pointsMatch;
33
34 int[] prev=array1, next=array2;
35 int maxScore=-999999;
36 int maxQpos=-1;
37 int maxRpos=-1;
38
39 for(int i=0; i<=arrayLength-query.length; i++){prev[i]=0;}
40 for(int i=arrayLength-query.length, score=0; i<=arrayLength; i++, score+=pointsDel) {
41 prev[i]=score;
42 }
43
44 int currentRstart=rstart;
45 for(int qpos=0; qpos<query.length; qpos++){
46 prev[0]=pointsIns*qpos;
47
48 final byte q=query[qpos];
49 int maxScoreThisPass=-9999;
50 final int remainingBases=(query.length-qpos);
51 final int remainingPoints=remainingBases*pointsMatch;
52 final int minViableScore=minPassingScore3-remainingPoints;
53 // int minViableRstart=rstop+1;
54 for(int rpos=currentRstart, apos=1+currentRstart-rstart; rpos<=rstop; rpos++, apos++){
55 final byte r=ref[rpos];
56 final boolean match=(q==r);
57 final int vScore=prev[apos]+pointsIns;
58 final int hScore=next[apos-1]+pointsDel;
59 final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1];
60
61 //Slow branchy code
62 // final int score=Tools.max(vScore, hScore, dScore);
63 // next[apos]=score;
64 // if(score>=maxScoreThisPass){
65 // maxScoreThisPass=score;
66 // if(score>=maxScore){
67 // maxScore=score;
68 // maxQpos=qpos;
69 // maxRpos=rpos;
70 // }
71 // }
72
73 //Should be branchless conditional moves
74 int score=(dScore>=vScore ? dScore : vScore);
75 score=(hScore>score ? hScore : score);
76 next[apos]=score;
77 maxScoreThisPass=(score>maxScoreThisPass ? score : maxScoreThisPass);
78
79 // minViableRstart=((score<minViableScore || rpos>minViableRstart) ? minViableRstart : rpos);
80 }
81 iters+=arrayLength;
82 // currentRstart=minViableRstart;
83 // System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore);
84 // System.err.println(Arrays.toString(prev));
85 // System.err.println(Arrays.toString(next));
86 if(maxScoreThisPass<minScore){//Aggressive early exit
87 // System.err.print('.');
88 // System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore);
89 return null;
90 }
91
92 {//Safe early exit
93 // if((maxScoreThisPass+remaining+query.length)/2<minPassingScore){return null;}
94 // if((maxScoreThisPass+remaining+query.length)<minPassingScore2){return null;}
95 if(maxScoreThisPass<minViableScore){return null;}
96 }
97
98 int[] temp=prev;
99 prev=next;
100 next=temp;
101 }
102 // System.err.print('*');
103
104 maxScore=-999999;
105 for(int rpos=rstart, apos=1; rpos<=rstop; rpos++, apos++){//Grab high score from last iteration
106 int score=prev[apos];
107 if(score>=maxScore){
108 maxScore=score;
109 maxQpos=query.length;
110 maxRpos=rpos;
111 }
112 }
113
114
115 // maxScore=(maxScore+query.length)/2;//Rescale 0 to length
116 maxScore=(maxScore-pointsSub*query.length)/(pointsMatch-pointsSub);//Rescale 0 to length
117 final float ratio=maxScore/(float)query.length;
118 // System.err.println("maxqPos="+maxQpos+", maxScore="+maxScore+", minScore="+(minRatio*query.length));
119 // if(minRatio*query.length>maxScore){return null;}
120 if(ratio<minRatio){return null;}
121 // System.err.print('^');
122 return new AlignmentResult(maxScore, maxQpos, maxRpos, query.length, ref.length, rstart, rstop, ratio);
123 }
124
125 /**
126 * @param query
127 * @param ref
128 * @param qstart
129 * @param rstart
130 * @param rstop
131 * @param minScore Quit early if score drops below this
132 * @param minRatio Don't return results if max score is less than this fraction of max possible score
133 * @return
134 */
135 public AlignmentResult alignForwardShort(final byte[] query, final byte[] ref, final int rstart, final int rstop,
136 final float minRatio) {
137 // if(ref.length<16300){return ica16.alignForwardShort(query, ref, qstart, rstart, rstop, minScore, minRatio);}
138
139 final int arrayLength=query.length;
140 // if(array1==null || arrayLength+1>array1.length){
141 // array1=new int[(int)((arrayLength+2)*1.2)];
142 // array2=new int[(int)((arrayLength+2)*1.2)];
143 // }
144
145 //Stack allocated; faster.
146 final int[] array1=new int[arrayLength+1], array2=new int[arrayLength+1];
147
148 final int minPassingScore=(int)(pointsMatch*query.length*minRatio);
149
150 int[] prev=array1, next=array2;
151 int maxScore=-1;
152 int maxQpos=-1;
153 int maxRpos=-1;
154
155 for(int i=0; i<prev.length; i++) {
156 prev[i]=pointsIns*i;
157 }
158
159 for(int rpos=rstart; rpos<=rstop && maxScore<minPassingScore; rpos++){
160 if(rstop-rpos<arrayLength){prev[0]=next[0]+pointsDel;}
161
162 final byte r=ref[rpos];
163 int score=0;
164 for(int qpos=0, apos=1; qpos<arrayLength; qpos++, apos++){
165 final byte q=query[qpos];
166 final boolean match=(q==r);
167 final int vScore=prev[apos]+pointsIns;
168 final int hScore=next[apos-1]+pointsDel;
169 final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1];
170
171 // score=Tools.max(vScore, hScore, dScore);
172 score=(dScore>=vScore ? dScore : vScore);
173 score=(hScore>score ? hScore : score);
174
175 next[apos]=score;
176 }
177 itersShort+=arrayLength;
178
179 if(score>=maxScore){
180 maxScore=score;
181 maxQpos=arrayLength-1;
182 maxRpos=rpos;
183 }
184
185 // System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore);
186 // System.err.println(Arrays.toString(prev));
187 // System.err.println(Arrays.toString(next));
188 // if(maxScoreThisPass<minScore){//Aggressive early exit
189 //// System.err.print('.');
190 //// System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore);
191 // return null;
192 // }
193
194 // {//Safe early exit
195 // int remaining=query.length-qpos;
196 //// if((maxScoreThisPass+remaining+query.length)/2<minPassingScore){return null;}
197 //// if((maxScoreThisPass+remaining+query.length)<minPassingScore2){return null;}
198 // if((maxScoreThisPass+remaining)<minPassingScore3){return null;}
199 // }
200
201 int[] temp=prev;
202 prev=next;
203 next=temp;
204 }
205 // System.err.print('*');
206
207 // maxScore=-999999;
208 // for(int rpos=rstart, apos=1; rpos<=rstop; rpos++, apos++){//Grab high score from last iteration
209 // int score=prev[apos];
210 // if(score>=maxScore){
211 // maxScore=score;
212 // maxQpos=query.length;
213 // maxRpos=rpos;
214 // }
215 // }
216 final float ratio=maxScore/(float)(pointsMatch*query.length);
217 // System.err.println(ratio+"\t"+maxScore+"\t"+query.length);
218 // System.err.println("maxqPos="+maxQpos+", maxScore="+maxScore+", minScore="+(minRatio*query.length));
219 // if(minRatio*query.length>maxScore){return null;}
220 if(ratio<minRatio){return null;}
221 // System.err.print('^');
222 return new AlignmentResult(maxScore, maxQpos, maxRpos, query.length, ref.length, rstart, rstop, ratio);
223 }
224
225 /*--------------------------------------------------------------*/
226 /*---------------- Getters ----------------*/
227 /*--------------------------------------------------------------*/
228
229 long iters(){return iters;}
230 long itersShort(){return itersShort;}
231
232 /*--------------------------------------------------------------*/
233 /*---------------- Fields ----------------*/
234 /*--------------------------------------------------------------*/
235
236 long iters = 0;
237 long itersShort = 0;
238
239 /*--------------------------------------------------------------*/
240 /*---------------- Constants ----------------*/
241 /*--------------------------------------------------------------*/
242
243 // public static final int pointsMatch = 10;
244 // public static final int pointsSub = -11;
245 // public static final int pointsDel = -9;
246 // public static final int pointsIns = -9;
247
248 public static final int pointsMatch = 120;
249 public static final int pointsSub = -123;
250 public static final int pointsDel = -187;
251 public static final int pointsIns = -243;
252
253 }