jpayne@68: package aligner;
jpayne@68: 
jpayne@68: 
jpayne@68: /** FlatAligner with flatter weights */
jpayne@68: public final class FlatAligner2 {
jpayne@68: 
jpayne@68: 	public FlatAligner2(){}
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * @param query
jpayne@68: 	 * @param ref
jpayne@68: 	 * @param qstart
jpayne@68: 	 * @param rstart
jpayne@68: 	 * @param rstop
jpayne@68: 	 * @param minScore Quit early if score drops below this
jpayne@68: 	 * @param minRatio Don't return results if max score is less than this fraction of max possible score
jpayne@68: 	 * @return
jpayne@68: 	 */
jpayne@68: 	public AlignmentResult alignForward(final byte[] query, final byte[] ref, final int rstart, final int rstop, final int minScore,
jpayne@68: 			final float minRatio) {
jpayne@68: //		if(ref.length<16300){return ica16.alignForward(query, ref, qstart, rstart, rstop, minScore, minRatio);}
jpayne@68: 		final int arrayLength=rstop-rstart+1;
jpayne@68: 		//		if(array1==null || arrayLength+1>array1.length){
jpayne@68: 		//			array1=new int[(int)((arrayLength+2)*1.2)];
jpayne@68: 		//			array2=new int[(int)((arrayLength+2)*1.2)];
jpayne@68: 		//		}
jpayne@68: 
jpayne@68: 		//Stack allocated; faster.
jpayne@68: 		final int[] array1=new int[arrayLength+1], array2=new int[arrayLength+1];
jpayne@68: 
jpayne@68: 		final int minPassingScore=(int)(query.length*minRatio*pointsMatch);
jpayne@68: 		final int minPassingScore3=minPassingScore-query.length*pointsMatch;
jpayne@68: 
jpayne@68: 		int[] prev=array1, next=array2;
jpayne@68: 		int maxScore=-999999;
jpayne@68: 		int maxQpos=-1;
jpayne@68: 		int maxRpos=-1;
jpayne@68: 
jpayne@68: 		for(int i=0; i<=arrayLength-query.length; i++){prev[i]=0;}
jpayne@68: 		for(int i=arrayLength-query.length, score=0; i<=arrayLength; i++, score+=pointsDel) {
jpayne@68: 			prev[i]=score;
jpayne@68: 		}
jpayne@68: 
jpayne@68: 		int currentRstart=rstart;
jpayne@68: 		for(int qpos=0; qpos<query.length; qpos++){
jpayne@68: 			prev[0]=pointsIns*qpos;
jpayne@68: 
jpayne@68: 			final byte q=query[qpos];
jpayne@68: 			int maxScoreThisPass=-9999;
jpayne@68: 			final int remainingBases=(query.length-qpos);
jpayne@68: 			final int remainingPoints=remainingBases*pointsMatch;
jpayne@68: 			final int minViableScore=minPassingScore3-remainingPoints;
jpayne@68: 			//			int minViableRstart=rstop+1;
jpayne@68: 			for(int rpos=currentRstart, apos=1+currentRstart-rstart; rpos<=rstop; rpos++, apos++){
jpayne@68: 				final byte r=ref[rpos];
jpayne@68: 				final boolean match=(q==r);
jpayne@68: 				final int vScore=prev[apos]+pointsIns;
jpayne@68: 				final int hScore=next[apos-1]+pointsDel;
jpayne@68: 				final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1];
jpayne@68: 
jpayne@68: 				//Slow branchy code
jpayne@68: 				//				final int score=Tools.max(vScore, hScore, dScore);
jpayne@68: 				//				next[apos]=score;
jpayne@68: 				//				if(score>=maxScoreThisPass){
jpayne@68: 				//					maxScoreThisPass=score;
jpayne@68: 				//					if(score>=maxScore){
jpayne@68: 				//						maxScore=score;
jpayne@68: 				//						maxQpos=qpos;
jpayne@68: 				//						maxRpos=rpos;
jpayne@68: 				//					}
jpayne@68: 				//				}
jpayne@68: 
jpayne@68: 				//Should be branchless conditional moves
jpayne@68: 				int score=(dScore>=vScore ? dScore : vScore);
jpayne@68: 				score=(hScore>score ? hScore : score);
jpayne@68: 				next[apos]=score;
jpayne@68: 				maxScoreThisPass=(score>maxScoreThisPass ? score : maxScoreThisPass);
jpayne@68: 
jpayne@68: 				//				minViableRstart=((score<minViableScore || rpos>minViableRstart) ? minViableRstart : rpos);
jpayne@68: 			}
jpayne@68: 			iters+=arrayLength;
jpayne@68: 			//			currentRstart=minViableRstart;
jpayne@68: 			//			System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore);
jpayne@68: 			//			System.err.println(Arrays.toString(prev));
jpayne@68: 			//			System.err.println(Arrays.toString(next));
jpayne@68: 			if(maxScoreThisPass<minScore){//Aggressive early exit
jpayne@68: 				//				System.err.print('.');
jpayne@68: 				//				System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore);
jpayne@68: 				return null;
jpayne@68: 			}
jpayne@68: 
jpayne@68: 			{//Safe early exit
jpayne@68: 				//				if((maxScoreThisPass+remaining+query.length)/2<minPassingScore){return null;}
jpayne@68: 				//				if((maxScoreThisPass+remaining+query.length)<minPassingScore2){return null;}
jpayne@68: 				if(maxScoreThisPass<minViableScore){return null;}
jpayne@68: 			}
jpayne@68: 
jpayne@68: 			int[] temp=prev;
jpayne@68: 			prev=next;
jpayne@68: 			next=temp;
jpayne@68: 		}
jpayne@68: 		//		System.err.print('*');
jpayne@68: 
jpayne@68: 		maxScore=-999999;
jpayne@68: 		for(int rpos=rstart, apos=1; rpos<=rstop; rpos++, apos++){//Grab high score from last iteration
jpayne@68: 			int score=prev[apos];
jpayne@68: 			if(score>=maxScore){
jpayne@68: 				maxScore=score;
jpayne@68: 				maxQpos=query.length;
jpayne@68: 				maxRpos=rpos;
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 
jpayne@68: 
jpayne@68: 		//		maxScore=(maxScore+query.length)/2;//Rescale 0 to length
jpayne@68: 		maxScore=(maxScore-pointsSub*query.length)/(pointsMatch-pointsSub);//Rescale 0 to length
jpayne@68: 		final float ratio=maxScore/(float)query.length;
jpayne@68: 		//		System.err.println("maxqPos="+maxQpos+", maxScore="+maxScore+", minScore="+(minRatio*query.length));
jpayne@68: 		//					if(minRatio*query.length>maxScore){return null;}
jpayne@68: 		if(ratio<minRatio){return null;}
jpayne@68: 		//		System.err.print('^');
jpayne@68: 		return new AlignmentResult(maxScore, maxQpos, maxRpos, query.length, ref.length, rstart, rstop, ratio);
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	/**
jpayne@68: 	 * @param query
jpayne@68: 	 * @param ref
jpayne@68: 	 * @param qstart
jpayne@68: 	 * @param rstart
jpayne@68: 	 * @param rstop
jpayne@68: 	 * @param minScore Quit early if score drops below this
jpayne@68: 	 * @param minRatio Don't return results if max score is less than this fraction of max possible score
jpayne@68: 	 * @return
jpayne@68: 	 */
jpayne@68: 	public AlignmentResult alignForwardShort(final byte[] query, final byte[] ref, final int rstart, final int rstop,
jpayne@68: 			final float minRatio) {
jpayne@68: //		if(ref.length<16300){return ica16.alignForwardShort(query, ref, qstart, rstart, rstop, minScore, minRatio);}
jpayne@68: 
jpayne@68: 		final int arrayLength=query.length;
jpayne@68: 		//		if(array1==null || arrayLength+1>array1.length){
jpayne@68: 		//			array1=new int[(int)((arrayLength+2)*1.2)];
jpayne@68: 		//			array2=new int[(int)((arrayLength+2)*1.2)];
jpayne@68: 		//		}
jpayne@68: 
jpayne@68: 		//Stack allocated; faster.
jpayne@68: 		final int[] array1=new int[arrayLength+1], array2=new int[arrayLength+1];
jpayne@68: 
jpayne@68: 		final int minPassingScore=(int)(pointsMatch*query.length*minRatio);
jpayne@68: 
jpayne@68: 		int[] prev=array1, next=array2;
jpayne@68: 		int maxScore=-1;
jpayne@68: 		int maxQpos=-1;
jpayne@68: 		int maxRpos=-1;
jpayne@68: 
jpayne@68: 		for(int i=0; i<prev.length; i++) {
jpayne@68: 			prev[i]=pointsIns*i;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		for(int rpos=rstart; rpos<=rstop && maxScore<minPassingScore; rpos++){
jpayne@68: 			if(rstop-rpos<arrayLength){prev[0]=next[0]+pointsDel;}
jpayne@68: 
jpayne@68: 			final byte r=ref[rpos];
jpayne@68: 			int score=0;
jpayne@68: 			for(int qpos=0, apos=1; qpos<arrayLength; qpos++, apos++){
jpayne@68: 				final byte q=query[qpos];
jpayne@68: 				final boolean match=(q==r);
jpayne@68: 				final int vScore=prev[apos]+pointsIns;
jpayne@68: 				final int hScore=next[apos-1]+pointsDel;
jpayne@68: 				final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1];
jpayne@68: 
jpayne@68: 				//				score=Tools.max(vScore, hScore, dScore);
jpayne@68: 				score=(dScore>=vScore ? dScore : vScore);
jpayne@68: 				score=(hScore>score ? hScore : score);
jpayne@68: 
jpayne@68: 				next[apos]=score;
jpayne@68: 			}
jpayne@68: 			itersShort+=arrayLength;
jpayne@68: 
jpayne@68: 			if(score>=maxScore){
jpayne@68: 				maxScore=score;
jpayne@68: 				maxQpos=arrayLength-1;
jpayne@68: 				maxRpos=rpos;
jpayne@68: 			}
jpayne@68: 
jpayne@68: 			//			System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore);
jpayne@68: 			//			System.err.println(Arrays.toString(prev));
jpayne@68: 			//			System.err.println(Arrays.toString(next));
jpayne@68: 			//			if(maxScoreThisPass<minScore){//Aggressive early exit
jpayne@68: 			////				System.err.print('.');
jpayne@68: 			////				System.err.println("qPos="+qpos+", maxScoreThisPass="+maxScoreThisPass+", maxScore="+maxScore+", minScore="+minScore);
jpayne@68: 			//				return null;
jpayne@68: 			//			}
jpayne@68: 
jpayne@68: 			//			{//Safe early exit
jpayne@68: 			//				int remaining=query.length-qpos;
jpayne@68: 			////				if((maxScoreThisPass+remaining+query.length)/2<minPassingScore){return null;}
jpayne@68: 			////				if((maxScoreThisPass+remaining+query.length)<minPassingScore2){return null;}
jpayne@68: 			//				if((maxScoreThisPass+remaining)<minPassingScore3){return null;}
jpayne@68: 			//			}
jpayne@68: 
jpayne@68: 			int[] temp=prev;
jpayne@68: 			prev=next;
jpayne@68: 			next=temp;
jpayne@68: 		}
jpayne@68: 		//		System.err.print('*');
jpayne@68: 
jpayne@68: 		//		maxScore=-999999;
jpayne@68: 		//		for(int rpos=rstart, apos=1; rpos<=rstop; rpos++, apos++){//Grab high score from last iteration
jpayne@68: 		//			int score=prev[apos];
jpayne@68: 		//			if(score>=maxScore){
jpayne@68: 		//				maxScore=score;
jpayne@68: 		//				maxQpos=query.length;
jpayne@68: 		//				maxRpos=rpos;
jpayne@68: 		//			}
jpayne@68: 		//		}
jpayne@68: 		final float ratio=maxScore/(float)(pointsMatch*query.length);
jpayne@68: //		System.err.println(ratio+"\t"+maxScore+"\t"+query.length);
jpayne@68: 		//		System.err.println("maxqPos="+maxQpos+", maxScore="+maxScore+", minScore="+(minRatio*query.length));
jpayne@68: 		//					if(minRatio*query.length>maxScore){return null;}
jpayne@68: 		if(ratio<minRatio){return null;}
jpayne@68: 		//		System.err.print('^');
jpayne@68: 		return new AlignmentResult(maxScore, maxQpos, maxRpos, query.length, ref.length, rstart, rstop, ratio);
jpayne@68: 	}
jpayne@68: 
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------           Getters            ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	long iters(){return iters;}
jpayne@68: 	long itersShort(){return itersShort;}
jpayne@68: 
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------            Fields            ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	long iters = 0;
jpayne@68: 	long itersShort = 0;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------           Constants          ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	public static final int pointsMatch = 10;
jpayne@68: 	public static final int pointsSub = -9;
jpayne@68: 	public static final int pointsDel = -11;
jpayne@68: 	public static final int pointsIns = -11;
jpayne@68: 	
jpayne@68: }