1 package aligner;
2 
3 
4 
5 public final class FlatAligner {
6 
FlatAligner()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 	 */
alignForward(final byte[] query, final byte[] ref, final int rstart, final int rstop, final int minScore, final float minRatio)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 	 */
alignForwardShort(final byte[] query, final byte[] ref, final int rstart, final int rstop, final float minRatio)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 
iters()229 	long iters(){return iters;}
itersShort()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 }
254