jpayne@68: package icecream;
jpayne@68: 
jpayne@68: import aligner.AlignmentResult;
jpayne@68: 
jpayne@68: public final class IceCreamAlignerJava extends IceCreamAligner {
jpayne@68: 
jpayne@68: 	IceCreamAlignerJava(){}
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * @param query
jpayne@68: 	 * @param ref
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: 	@Override
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: 		final int arrayLength=rstop-rstart+1;
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: 			if(maxScoreThisPass<minScore){//Aggressive early exit
jpayne@68: 				return null;
jpayne@68: 			}
jpayne@68: 
jpayne@68: 			{//Safe early exit
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: 
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: 		maxScore=(maxScore-pointsSub*query.length)/(pointsMatch-pointsSub);//Rescale 0 to length
jpayne@68: 		final float ratio=maxScore/(float)query.length;
jpayne@68: 		
jpayne@68: 		if(ratio<minRatio){return null;}
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: 	@Override
jpayne@68: 	public AlignmentResult alignForwardShort(final byte[] query, final byte[] ref, final int rstart, final int rstop, final int minScore,
jpayne@68: 			final float minRatio) {
jpayne@68: 		
jpayne@68: 		final int arrayLength=query.length;
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);
jpayne@68: 		final int minPassingScore3=minPassingScore-query.length;
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<prev.length; i++) {
jpayne@68: 			prev[i]=pointsIns*i;
jpayne@68: 		}
jpayne@68: 
jpayne@68: 		for(int rpos=rstart; rpos<=rstop; 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: 			
jpayne@68: 			//In Java, the unsplit loop is faster
jpayne@68: 			//-XX:+UseSuperWord seems to have no impact on speed.
jpayne@68: 			//Might be worth packing the query and ref into int arrays and trying again.
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=(dScore>=vScore ? dScore : vScore);
jpayne@68: 				score=(hScore>score ? hScore : score);
jpayne@68: 
jpayne@68: 				next[apos]=score;
jpayne@68: 			}
jpayne@68: 
jpayne@68: //			//Inner DV loop
jpayne@68: //			for(int qpos=0, apos=1; qpos<arrayLength; qpos++, apos++){
jpayne@68: //				final int q=query[qpos];
jpayne@68: //				final boolean match=(q==r);
jpayne@68: //				final int vScore=prev[apos]+pointsIns;
jpayne@68: //				final int dScore=(match ? pointsMatch : pointsSub)+prev[apos-1];
jpayne@68: //
jpayne@68: //				score=(dScore>=vScore ? dScore : vScore);
jpayne@68: //
jpayne@68: //				next[apos]=score;
jpayne@68: //			}
jpayne@68: //			
jpayne@68: //			//Inner I loop
jpayne@68: //			for(int qpos=0, apos=1; qpos<arrayLength; qpos++, apos++){
jpayne@68: //				final int hScore=next[apos-1]+pointsDel;
jpayne@68: //
jpayne@68: //				score=next[apos];
jpayne@68: //				score=(hScore>score ? hScore : score);
jpayne@68: //
jpayne@68: //				next[apos]=score;
jpayne@68: //			}
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: 			int[] temp=prev;
jpayne@68: 			prev=next;
jpayne@68: 			next=temp;
jpayne@68: 		}
jpayne@68: 
jpayne@68: 		maxScore=(maxScore-pointsSub*query.length)/(pointsMatch-pointsSub);//Rescale 0 to length
jpayne@68: 		final float ratio=maxScore/(float)query.length;
jpayne@68: 		
jpayne@68: 		if(ratio<minRatio){return null;}
jpayne@68: 		
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: 	@Override
jpayne@68: 	long iters(){return iters;}
jpayne@68: 
jpayne@68: 	@Override
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 = 1;
jpayne@68: 	public static final int pointsSub = -1;
jpayne@68: 	public static final int pointsDel = -2;
jpayne@68: 	public static final int pointsIns = -2;
jpayne@68: 	
jpayne@68: }