Mercurial > repos > rliterman > csp2

annotate CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/sketch/GlocalAligner.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 sketch;
jpayne@68 2
jpayne@68 3 import shared.Tools;
jpayne@68 4
jpayne@68 5 public final class GlocalAligner {
jpayne@68 6
jpayne@68 7 GlocalAligner(){}
jpayne@68 8
jpayne@68 9 public float alignForward(final byte[] query, final byte[] ref){
jpayne@68 10 float a=alignForwardInner(query, ref);
jpayne@68 11 float b=alignForwardInner(ref, query);
jpayne@68 12 return Tools.max(a, b);
jpayne@68 13 }
jpayne@68 14
jpayne@68 15 /**
jpayne@68 16 * @param query
jpayne@68 17 * @param ref
jpayne@68 18 * @return Identity
jpayne@68 19 */
jpayne@68 20 public float alignForwardInner(final byte[] query, final byte[] ref){
jpayne@68 21 // if(ref.length<query.length){return alignForward(ref, query);}
jpayne@68 22 final int big=ref.length;
jpayne@68 23 final int small=query.length;
jpayne@68 24 // assert(big>=small);
jpayne@68 25 final int arrayLength=big;
jpayne@68 26
jpayne@68 27 //Stack allocated; faster.
jpayne@68 28 final short[] array1=new short[arrayLength+1], array2=new short[arrayLength+1];
jpayne@68 29 final short[] consumed1=new short[arrayLength+1], consumed2=new short[arrayLength+1];
jpayne@68 30
jpayne@68 31 short[] prev=array1, next=array2;
jpayne@68 32 short[] prevConsumed=consumed1, nextConsumed=consumed2;
jpayne@68 33
jpayne@68 34 // for(int i=0; i<=arrayLength-query.length; i++){prev[i]=0;}
jpayne@68 35 // for(short i=(short)(arrayLength-query.length), score=0; i<=arrayLength; i++, score+=pointsDel) {
jpayne@68 36 // prev[i]=score;
jpayne@68 37 // }
jpayne@68 38 final int rmax=ref.length-1;
jpayne@68 39 final int qmax=query.length-1;
jpayne@68 40
jpayne@68 41 int maxScore=0;
jpayne@68 42 int maxConsumed=0;
jpayne@68 43 for(short qpos=0; qpos<query.length; qpos++){
jpayne@68 44 // prev[0]=(short)(pointsIns*qpos);
jpayne@68 45
jpayne@68 46 final byte q=query[qpos];
jpayne@68 47 for(int rpos=0, apos=1; rpos<ref.length; rpos++, apos++){
jpayne@68 48 final byte r=ref[rpos];
jpayne@68 49 final boolean match=(q==r && q!='N');
jpayne@68 50 final boolean rEnd=(rpos<1 || rpos>=rmax);
jpayne@68 51 final boolean qEnd=(qpos<1 || qpos>=qmax);
jpayne@68 52 final short vScore=(short) (prev[apos]+(rEnd ? 0 : pointsIns));
jpayne@68 53 final short hScore=(short) (next[apos-1]+(qEnd ? 0 : pointsDel));
jpayne@68 54 final short dScore=(short) ((match ? pointsMatch : pointsSub)+prev[apos-1]);
jpayne@68 55
jpayne@68 56 short score, consumed;
jpayne@68 57 if(dScore>=vScore && dScore>=hScore){
jpayne@68 58 score=dScore;
jpayne@68 59 consumed=(short)(prevConsumed[apos-1]+1);
jpayne@68 60 }else if(vScore>=hScore){
jpayne@68 61 score=vScore;
jpayne@68 62 // nextConsumed[apos]=(short)(prevConsumed[apos]+(rEnd ? 0 : 1));
jpayne@68 63 consumed=(short)(prevConsumed[apos]);
jpayne@68 64 }else{
jpayne@68 65 score=vScore;
jpayne@68 66 consumed=(short)(nextConsumed[apos-1]);
jpayne@68 67 }
jpayne@68 68
jpayne@68 69 if(score<0){
jpayne@68 70 score=0;
jpayne@68 71 consumed=0;
jpayne@68 72 }
jpayne@68 73 if(score>maxScore){
jpayne@68 74 maxScore=score;
jpayne@68 75 maxConsumed=consumed;
jpayne@68 76 }
jpayne@68 77 next[apos]=score;
jpayne@68 78 nextConsumed[apos]=consumed;
jpayne@68 79 // //Should be branchless conditional moves
jpayne@68 80 // short score=(dScore>=vScore ? dScore : vScore);
jpayne@68 81 // score=(hScore>score ? hScore : score);
jpayne@68 82 // next[apos]=score;
jpayne@68 83 }
jpayne@68 84 // iters+=arrayLength;
jpayne@68 85
jpayne@68 86 short[] temp=prev;
jpayne@68 87 prev=next;
jpayne@68 88 next=temp;
jpayne@68 89 temp=prevConsumed;
jpayne@68 90 prevConsumed=nextConsumed;
jpayne@68 91 nextConsumed=temp;
jpayne@68 92 }
jpayne@68 93
jpayne@68 94 // short maxScore=Short.MIN_VALUE;
jpayne@68 95 // short maxConsumed=Short.MIN_VALUE;
jpayne@68 96 // for(short apos=1; apos<prev.length; apos++){//Grab high score from last iteration
jpayne@68 97 // short score=prev[apos];
jpayne@68 98 // if(score>=maxScore){
jpayne@68 99 // maxScore=score;
jpayne@68 100 // maxConsumed=prevConsumed[apos];
jpayne@68 101 // }
jpayne@68 102 // }
jpayne@68 103
jpayne@68 104 // assert(false);
jpayne@68 105 // System.err.println("maxScore="+maxScore+"; maxConsumed="+maxConsumed);
jpayne@68 106 if(maxConsumed<400){return 0;}
jpayne@68 107
jpayne@68 108 // int maxPossibleScore=(small*(pointsMatch-pointsSub));
jpayne@68 109 // int rescaledScore=maxScore-small*pointsSub;
jpayne@68 110 // final float ratio=rescaledScore/(float)maxPossibleScore;
jpayne@68 111 // return ratio;
jpayne@68 112
jpayne@68 113 int maxPossibleScore=(maxConsumed*(pointsMatch-pointsSub));
jpayne@68 114 int rescaledScore=maxScore-maxConsumed*pointsSub;
jpayne@68 115 final float ratio=rescaledScore/(float)maxPossibleScore;
jpayne@68 116 return ratio;
jpayne@68 117 }
jpayne@68 118
jpayne@68 119
jpayne@68 120 /*--------------------------------------------------------------*/
jpayne@68 121 /*---------------- Getters ----------------*/
jpayne@68 122 /*--------------------------------------------------------------*/
jpayne@68 123
jpayne@68 124 /*--------------------------------------------------------------*/
jpayne@68 125 /*---------------- Fields ----------------*/
jpayne@68 126 /*--------------------------------------------------------------*/
jpayne@68 127
jpayne@68 128 long iters=0;
jpayne@68 129
jpayne@68 130 /*--------------------------------------------------------------*/
jpayne@68 131 /*---------------- Constants ----------------*/
jpayne@68 132 /*--------------------------------------------------------------*/
jpayne@68 133
jpayne@68 134 public static final short pointsMatch = 1;
jpayne@68 135 public static final short pointsSub = -1;
jpayne@68 136 public static final short pointsDel = -1;
jpayne@68 137 public static final short pointsIns = -1;
jpayne@68 138
jpayne@68 139 }