Mercurial > repos > rliterman > csp2

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