Mercurial > repos > rliterman > csp2

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