Mercurial > repos > rliterman > csp2
diff CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/aligner/FlatAligner2.java @ 68:5028fdace37b
planemo upload commit 2e9511a184a1ca667c7be0c6321a36dc4e3d116d
author | jpayne |
---|---|
date | Tue, 18 Mar 2025 16:23:26 -0400 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/aligner/FlatAligner2.java Tue Mar 18 16:23:26 2025 -0400 @@ -0,0 +1,248 @@ +package aligner; + + +/** FlatAligner with flatter weights */ +public final class FlatAligner2 { + + public FlatAligner2(){} + + /** + * @param query + * @param ref + * @param qstart + * @param rstart + * @param rstop + * @param minScore Quit early if score drops below this + * @param minRatio Don't return results if max score is less than this fraction of max possible score + * @return + */ + public AlignmentResult alignForward(final byte[] query, final byte[] ref, final int rstart, final int rstop, final int minScore, + final float minRatio) { +// if(ref.length<16300){return ica16.alignForward(query, ref, qstart, rstart, rstop, minScore, minRatio);} + final int arrayLength=rstop-rstart+1; + // if(array1==null || arrayLength+1>array1.length){ + // array1=new int[(int)((arrayLength+2)*1.2)]; + // array2=new int[(int)((arrayLength+2)*1.2)]; + // } + + //Stack allocated; faster. + final int[] array1=new int[arrayLength+1], array2=new int[arrayLength+1]; + + final int minPassingScore=(int)(query.length*minRatio*pointsMatch); + final int minPassingScore3=minPassingScore-query.length*pointsMatch; + + int[] prev=array1, next=array2; + int maxScore=-999999; + int maxQpos=-1; + int maxRpos=-1; + + for(int i=0; i<=arrayLength-query.length; i++){prev[i]=0;} + for(int i=arrayLength-query.length, score=0; i<=arrayLength; i++, score+=pointsDel) { + prev[i]=score; + } + + int currentRstart=rstart; + for(int qpos=0; qpos<query.length; qpos++){ + prev[0]=pointsIns*qpos; + + final byte q=query[qpos]; + int maxScoreThisPass=-9999; + final int remainingBases=(query.length-qpos); + final int remainingPoints=remainingBases*pointsMatch; + final int minViableScore=minPassingScore3-remainingPoints; + // int minViableRstart=rstop+1; + for(int rpos=currentRstart, apos=1+currentRstart-rstart; rpos<=rstop; rpos++, apos++){ + final byte r=ref[rpos]; + final boolean match=(q==r); + final int vScore=prev[apos]+pointsIns; + final int hScore=next[apos-1]+pointsDel; + final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1]; + + //Slow branchy code + // final int score=Tools.max(vScore, hScore, dScore); + // next[apos]=score; + // if(score>=maxScoreThisPass){ + // maxScoreThisPass=score; + // if(score>=maxScore){ + // maxScore=score; + // maxQpos=qpos; + // maxRpos=rpos; + // } + // } + + //Should be branchless conditional moves + int score=(dScore>=vScore ? dScore : vScore); + score=(hScore>score ? hScore : score); + next[apos]=score; + maxScoreThisPass=(score>maxScoreThisPass ? score : maxScoreThisPass); + + // minViableRstart=((score<minViableScore || rpos>minViableRstart) ? minViableRstart : rpos); + } + iters+=arrayLength; + // currentRstart=minViableRstart; + // System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore); + // System.err.println(Arrays.toString(prev)); + // System.err.println(Arrays.toString(next)); + if(maxScoreThisPass<minScore){//Aggressive early exit + // System.err.print('.'); + // System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore); + return null; + } + + {//Safe early exit + // if((maxScoreThisPass+remaining+query.length)/2<minPassingScore){return null;} + // if((maxScoreThisPass+remaining+query.length)<minPassingScore2){return null;} + if(maxScoreThisPass<minViableScore){return null;} + } + + int[] temp=prev; + prev=next; + next=temp; + } + // System.err.print('*'); + + maxScore=-999999; + for(int rpos=rstart, apos=1; rpos<=rstop; rpos++, apos++){//Grab high score from last iteration + int score=prev[apos]; + if(score>=maxScore){ + maxScore=score; + maxQpos=query.length; + maxRpos=rpos; + } + } + + + // maxScore=(maxScore+query.length)/2;//Rescale 0 to length + maxScore=(maxScore-pointsSub*query.length)/(pointsMatch-pointsSub);//Rescale 0 to length + final float ratio=maxScore/(float)query.length; + // System.err.println("maxqPos="+maxQpos+", maxScore="+maxScore+", minScore="+(minRatio*query.length)); + // if(minRatio*query.length>maxScore){return null;} + if(ratio<minRatio){return null;} + // System.err.print('^'); + return new AlignmentResult(maxScore, maxQpos, maxRpos, query.length, ref.length, rstart, rstop, ratio); + } + + /** + * @param query + * @param ref + * @param qstart + * @param rstart + * @param rstop + * @param minScore Quit early if score drops below this + * @param minRatio Don't return results if max score is less than this fraction of max possible score + * @return + */ + public AlignmentResult alignForwardShort(final byte[] query, final byte[] ref, final int rstart, final int rstop, + final float minRatio) { +// if(ref.length<16300){return ica16.alignForwardShort(query, ref, qstart, rstart, rstop, minScore, minRatio);} + + final int arrayLength=query.length; + // if(array1==null || arrayLength+1>array1.length){ + // array1=new int[(int)((arrayLength+2)*1.2)]; + // array2=new int[(int)((arrayLength+2)*1.2)]; + // } + + //Stack allocated; faster. + final int[] array1=new int[arrayLength+1], array2=new int[arrayLength+1]; + + final int minPassingScore=(int)(pointsMatch*query.length*minRatio); + + int[] prev=array1, next=array2; + int maxScore=-1; + int maxQpos=-1; + int maxRpos=-1; + + for(int i=0; i<prev.length; i++) { + prev[i]=pointsIns*i; + } + + for(int rpos=rstart; rpos<=rstop && maxScore<minPassingScore; rpos++){ + if(rstop-rpos<arrayLength){prev[0]=next[0]+pointsDel;} + + final byte r=ref[rpos]; + int score=0; + for(int qpos=0, apos=1; qpos<arrayLength; qpos++, apos++){ + final byte q=query[qpos]; + final boolean match=(q==r); + final int vScore=prev[apos]+pointsIns; + final int hScore=next[apos-1]+pointsDel; + final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1]; + + // score=Tools.max(vScore, hScore, dScore); + score=(dScore>=vScore ? dScore : vScore); + score=(hScore>score ? hScore : score); + + next[apos]=score; + } + itersShort+=arrayLength; + + if(score>=maxScore){ + maxScore=score; + maxQpos=arrayLength-1; + maxRpos=rpos; + } + + // System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore); + // System.err.println(Arrays.toString(prev)); + // System.err.println(Arrays.toString(next)); + // if(maxScoreThisPass<minScore){//Aggressive early exit + //// System.err.print('.'); + //// System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore); + // return null; + // } + + // {//Safe early exit + // int remaining=query.length-qpos; + //// if((maxScoreThisPass+remaining+query.length)/2<minPassingScore){return null;} + //// if((maxScoreThisPass+remaining+query.length)<minPassingScore2){return null;} + // if((maxScoreThisPass+remaining)<minPassingScore3){return null;} + // } + + int[] temp=prev; + prev=next; + next=temp; + } + // System.err.print('*'); + + // maxScore=-999999; + // for(int rpos=rstart, apos=1; rpos<=rstop; rpos++, apos++){//Grab high score from last iteration + // int score=prev[apos]; + // if(score>=maxScore){ + // maxScore=score; + // maxQpos=query.length; + // maxRpos=rpos; + // } + // } + final float ratio=maxScore/(float)(pointsMatch*query.length); +// System.err.println(ratio+"\t"+maxScore+"\t"+query.length); + // System.err.println("maxqPos="+maxQpos+", maxScore="+maxScore+", minScore="+(minRatio*query.length)); + // if(minRatio*query.length>maxScore){return null;} + if(ratio<minRatio){return null;} + // System.err.print('^'); + return new AlignmentResult(maxScore, maxQpos, maxRpos, query.length, ref.length, rstart, rstop, ratio); + } + + /*--------------------------------------------------------------*/ + /*---------------- Getters ----------------*/ + /*--------------------------------------------------------------*/ + + long iters(){return iters;} + long itersShort(){return itersShort;} + + /*--------------------------------------------------------------*/ + /*---------------- Fields ----------------*/ + /*--------------------------------------------------------------*/ + + long iters = 0; + long itersShort = 0; + + /*--------------------------------------------------------------*/ + /*---------------- Constants ----------------*/ + /*--------------------------------------------------------------*/ + + public static final int pointsMatch = 10; + public static final int pointsSub = -9; + public static final int pointsDel = -11; + public static final int pointsIns = -11; + +}