diff CSP2/CSP2_env/env-d9b9114564458d9d-741b3de822f2aaca6c6caa4325c4afce/opt/bbmap-39.01-1/current/sketch/GlocalAligner.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/sketch/GlocalAligner.java	Tue Mar 18 16:23:26 2025 -0400
@@ -0,0 +1,139 @@
+package sketch;
+
+import shared.Tools;
+
+public final class GlocalAligner {
+
+	GlocalAligner(){}
+	
+	public float alignForward(final byte[] query, final byte[] ref){
+		float a=alignForwardInner(query, ref);
+		float b=alignForwardInner(ref, query);
+		return Tools.max(a, b);
+	}
+	
+	/**
+	 * @param query
+	 * @param ref
+	 * @return Identity
+	 */
+	public float alignForwardInner(final byte[] query, final byte[] ref){
+//		if(ref.length<query.length){return alignForward(ref, query);}
+		final int big=ref.length;
+		final int small=query.length;
+//		assert(big>=small);
+		final int arrayLength=big;
+
+		//Stack allocated; faster.
+		final short[] array1=new short[arrayLength+1], array2=new short[arrayLength+1];
+		final short[] consumed1=new short[arrayLength+1], consumed2=new short[arrayLength+1];
+
+		short[] prev=array1, next=array2;
+		short[] prevConsumed=consumed1, nextConsumed=consumed2;
+		
+//		for(int i=0; i<=arrayLength-query.length; i++){prev[i]=0;}
+//		for(short i=(short)(arrayLength-query.length), score=0; i<=arrayLength; i++, score+=pointsDel) {
+//			prev[i]=score;
+//		}
+		final int rmax=ref.length-1;
+		final int qmax=query.length-1;
+		
+		int maxScore=0;
+		int maxConsumed=0;
+		for(short qpos=0; qpos<query.length; qpos++){
+//			prev[0]=(short)(pointsIns*qpos);
+			
+			final byte q=query[qpos];
+			for(int rpos=0, apos=1; rpos<ref.length; rpos++, apos++){
+				final byte r=ref[rpos];
+				final boolean match=(q==r && q!='N');
+				final boolean rEnd=(rpos<1 || rpos>=rmax);
+				final boolean qEnd=(qpos<1 || qpos>=qmax);
+				final short vScore=(short) (prev[apos]+(rEnd ? 0 : pointsIns));
+				final short hScore=(short) (next[apos-1]+(qEnd ? 0 : pointsDel));
+				final short dScore=(short) ((match ? pointsMatch : pointsSub)+prev[apos-1]);
+				
+				short score, consumed;
+				if(dScore>=vScore && dScore>=hScore){
+					score=dScore;
+					consumed=(short)(prevConsumed[apos-1]+1);
+				}else if(vScore>=hScore){
+					score=vScore;
+//					nextConsumed[apos]=(short)(prevConsumed[apos]+(rEnd ? 0 : 1));
+					consumed=(short)(prevConsumed[apos]);
+				}else{
+					score=vScore;
+					consumed=(short)(nextConsumed[apos-1]);
+				}
+
+				if(score<0){
+					score=0;
+					consumed=0;
+				}
+				if(score>maxScore){
+					maxScore=score;
+					maxConsumed=consumed;
+				}
+				next[apos]=score;
+				nextConsumed[apos]=consumed;
+//				//Should be branchless conditional moves
+//				short score=(dScore>=vScore ? dScore : vScore);
+//				score=(hScore>score ? hScore : score);
+//				next[apos]=score;
+			}
+//			iters+=arrayLength;
+			
+			short[] temp=prev;
+			prev=next;
+			next=temp;
+			temp=prevConsumed;
+			prevConsumed=nextConsumed;
+			nextConsumed=temp;
+		}
+
+//		short maxScore=Short.MIN_VALUE;
+//		short maxConsumed=Short.MIN_VALUE;
+//		for(short apos=1; apos<prev.length; apos++){//Grab high score from last iteration
+//			short score=prev[apos];
+//			if(score>=maxScore){
+//				maxScore=score;
+//				maxConsumed=prevConsumed[apos];
+//			}
+//		}
+		
+//		assert(false);
+//		System.err.println("maxScore="+maxScore+"; maxConsumed="+maxConsumed);
+		if(maxConsumed<400){return 0;}
+		
+//		int maxPossibleScore=(small*(pointsMatch-pointsSub));
+//		int rescaledScore=maxScore-small*pointsSub;
+//		final float ratio=rescaledScore/(float)maxPossibleScore;
+//		return ratio;
+		
+		int maxPossibleScore=(maxConsumed*(pointsMatch-pointsSub));
+		int rescaledScore=maxScore-maxConsumed*pointsSub;
+		final float ratio=rescaledScore/(float)maxPossibleScore;
+		return ratio;
+	}
+
+
+	/*--------------------------------------------------------------*/
+	/*----------------           Getters            ----------------*/
+	/*--------------------------------------------------------------*/
+
+	/*--------------------------------------------------------------*/
+	/*----------------            Fields            ----------------*/
+	/*--------------------------------------------------------------*/
+	
+	long iters=0;
+	
+	/*--------------------------------------------------------------*/
+	/*----------------           Constants          ----------------*/
+	/*--------------------------------------------------------------*/
+	
+	public static final short pointsMatch = 1;
+	public static final short pointsSub = -1;
+	public static final short pointsDel = -1;
+	public static final short pointsIns = -1;
+	
+}