1 package aligner; 2 3 4 /** FlatAligner with flatter weights */ 5 public final class FlatAligner2 { 6 FlatAligner2()7 public FlatAligner2(){} 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 = -9; 245 public static final int pointsDel = -11; 246 public static final int pointsIns = -11; 247 248 } 249