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