//% Copyright(c) 2006, Shinichi Morishita. All Rights Reserved.

public class Chapter6_ApproximateStringMatching {

	final static double matchScore		= 1;
	final static double mismatchPenalty = -0.5;
	final static double gapPenalty 	 	= -0.5;
	final static double startupGapPenalty = -2;	
	
	public static void NeedlemanWunsch(int[] stringX, int[] stringY){
		// Initialize constants.
		int stringXlen = stringX.length;
		int stringYlen = stringY.length; 
		double[][] score = new double[stringXlen+1][stringYlen+1];
		// Initialize the score array.
		for(int x=0; x < stringXlen+1; x++) score[x][0] = gapPenalty*x;
		for(int y=0; y < stringYlen+1; y++) score[0][y] = gapPenalty*y;
		// Inductive steps.
		double oneScore; 
		for(int x=1; x < stringXlen+1; x++){
			for(int y=1; y < stringYlen+1; y++){
				if(stringX[x-1] == stringY[y-1]) 
					 oneScore = matchScore; 
				else oneScore = mismatchPenalty;
				score[x][y] = max(	score[x-1][y-1] + oneScore, 
									score[x-1][y]   + gapPenalty,
									score[x][y-1]   + gapPenalty);
			}
		}	
		// Print out the score.
		System.out.println("Alignment score = "+score[stringX.length][stringY.length]);
		// Restore the alignment from the bottom-right to the top-left.
		int[][] alignment = new int[stringXlen+stringYlen+1][2];
		int i = 0; int x,y;
		for(x = stringX.length, y = stringY.length; !(x == 0 && y == 0); i++){
			if(0 < x && 0 < y && stringX[x-1] == stringY[y-1]) 
				 oneScore = matchScore; 
			else oneScore = mismatchPenalty;
			if(0 < x && 0 < y && score[x][y] == score[x-1][y-1] + oneScore){
				alignment[i][0] = stringX[x-1]; alignment[i][1] = stringY[y-1];
				x = x-1; y = y-1;
			}else if(0 < x && 0 <= y && score[x][y] == score[x-1][y] + gapPenalty){
				alignment[i][0] = stringX[x-1]; alignment[i][1] = -1; x = x-1; 
			}else if(0 <= x && 0 < y && score[x][y] == score[x][y-1] + gapPenalty){
				alignment[i][0] = -1; alignment[i][1] = stringY[y-1]; y = y-1;
			}
		}
		printAlignment(alignment, i-1);
	}

	public static void NeedlemanWunschHirschberg(int[] stringX, int[] stringY){
		double[][] tmpScore = new double[3][stringY.length+1];
		int[] footprints = new int[stringX.length+1];  
		double[] footprintScores = new double[stringX.length+1]; 
		footprints[0] = 0;
		footprintScores[0] = 0;
		footprints[stringX.length] = stringY.length;
		NWHrecursion(footprints, footprintScores, 
				tmpScore, stringX, stringY, 0, 0, stringX.length, stringY.length);

		// Print the score of the optimal alignment.
		System.out.println("Alignment score = "+footprintScores[stringX.length]);
		// Print out the alignment by parsing footprints and footprintScores.
		int[][] alignment = new int[stringX.length+stringY.length+1][2];
		int x, j;   // x and j scan footprints and alignment.
		for( j=0, x=1; x < footprints.length; x++){
			if( footprints[x-1] == footprints[x]){ // Insert a gap in stringY.
				alignment[j][0] = stringX[x-1]; alignment[j][1] = -1; j++;
			}else{ // Insert gaps in stringX.
				double diffScores = footprintScores[x]-footprintScores[x-1];
				int diffFootprints = footprints[x]-footprints[x-1];	
				if(diffScores == matchScore + gapPenalty*(diffFootprints-1)){
					boolean match = false;
					for(int y=footprints[x-1]+1; y <= footprints[x]; y++, j++){
						if(stringX[x-1] == stringY[y-1] && !match){ 	
								alignment[j][0] = stringX[x-1]; match=true;
						}else	alignment[j][0] = -1; 
						alignment[j][1] = stringY[y-1];
					}
				}else if(diffScores == gapPenalty*(diffFootprints+1)){
					for(int y=footprints[x-1]+1; y <= footprints[x]; y++, j++){
						alignment[j][0] = -1; 
						alignment[j][1] = stringY[y-1];
					}
					alignment[j][0] = stringX[x-1]; alignment[j][1] = -1; j++;		 					
				}else if(diffScores == mismatchPenalty + gapPenalty*(diffFootprints-1)){
					for(int y=footprints[x-1]+1; y <= footprints[x]; y++, j++){
						if(y == footprints[x]) 	
								alignment[j][0] = stringX[x-1];
						else	alignment[j][0] = -1; 
						alignment[j][1] = stringY[y-1];
					}
				}			
			}	
		}
		for(int k=0; k<j; k++)	System.out.print(int2char(alignment[k][0])+" ");
		System.out.println();
		for(int k=0; k<j; k++)
			if(alignment[k][0] == -1 || alignment[k][1] == -1) System.out.print("  ");
			else if(alignment[k][0] == alignment[k][1]) System.out.print("| ");
			else System.out.print("x ");
		System.out.println();		
		for(int k=0; k<j; k++)	System.out.print(int2char(alignment[k][1])+" ");
		System.out.println();		
	}
	
	public static void NWHrecursion( 
				int[] footprints, double[] footprintScores, 
				double[][] tmpScore, int[] stringX, int[] stringY, 
				int topLeftX, int topLeftY, int bottomRightX, int bottomRightY){
		int middleX = (topLeftX + bottomRightX)/2;	// Probe the middle vertical line.
		if(middleX == topLeftX || middleX == bottomRightX) return;

		// Scan from the top-left. Initialize the score array.
		int prev = 0; int current = 1;
		for(int y = topLeftY; y <= bottomRightY; y++) 
			tmpScore[current][y] = gapPenalty*(y-topLeftY);
		// Inductive steps.
		double oneScore; 
		for( int x = topLeftX+1; x <= middleX; x++ ){	
			prev = (prev+1)%2; current = (current+1)%2; // Exchange prev and current.
			tmpScore[current][topLeftY] = gapPenalty*(x-topLeftX); 			
			for( int y = topLeftY+1; y <= bottomRightY; y++){
				if(stringX[x-1] == stringY[y-1]) oneScore = matchScore; 
				else 							 oneScore = mismatchPenalty;
				tmpScore[current][y] = 
					max(tmpScore[prev][y-1]   + oneScore,	// Diagonal.
						tmpScore[prev][y]     + gapPenalty,	// Horizontal.
						tmpScore[current][y-1]+ gapPenalty);	// Vertical.
			}
		}			

		// Set aside the scores of the middle vertical line.
		int firstHalf = 2;
		for(int y = topLeftY; y <= bottomRightY; y++) 
			tmpScore[firstHalf][y] = tmpScore[current][y];		

		// Scan from the bottom-right. Initialize the score array.
		prev = 0; current = 1;
		for(int y = bottomRightY; topLeftY <= y; y--) 
			tmpScore[current][y] = gapPenalty*(bottomRightY-y);
		// Inductive steps. 
		for( int x = bottomRightX-1; middleX <= x; x-- ){
			prev = (prev+1)%2; current = (current+1)%2; // Exchange prev and current.	
			tmpScore[current][bottomRightY] = gapPenalty*(bottomRightX-x); 			
			for( int y = bottomRightY-1; topLeftY <= y; y--){
				// Indexes of stringX and stringY are changed slightly.
				if(stringX[x] == stringY[y]) oneScore = matchScore; 
				else 						 oneScore = mismatchPenalty;
				tmpScore[current][y] = 
					max(tmpScore[prev][y+1]   + oneScore,	// Diagonal.
						tmpScore[prev][y]     + gapPenalty,	// Horizontal.
						tmpScore[current][y+1]+ gapPenalty);	// Vertical.
			}
		}

		// Set the footprint to the maximum sum in the middle column.
		int y = topLeftY; int middleY = topLeftY;
		double maxSum = tmpScore[firstHalf][y]+tmpScore[current][y]; 
		for( y++; y <= bottomRightY; y++){
			double currentSum = tmpScore[firstHalf][y]+tmpScore[current][y];
			if(maxSum < currentSum){ maxSum = currentSum; middleY = y; }
		}
		footprints[middleX] = middleY;
		footprintScores[middleX] = footprintScores[topLeftX]+tmpScore[firstHalf][middleY];
		footprintScores[bottomRightX] = footprintScores[topLeftX]+maxSum;
	
		// Recursive application of this procedure to the upper-left half and lower-right half blocks.
		NWHrecursion(footprints, footprintScores, tmpScore, 
				stringX, stringY, topLeftX, topLeftY, middleX, middleY);
		NWHrecursion(footprints, footprintScores, tmpScore, 
				stringX, stringY, middleX, middleY, bottomRightX, bottomRightY);
	}
	
	
	public static void SmithWaterman( int[] stringX, int[] stringY){
		// Initialize constants.
		int stringXlen = stringX.length;
		int stringYlen = stringY.length; 
		double[][] score = new double[stringXlen+1][stringYlen+1];	
		// Initialize the score array by setting the boundary elements to 0.
		for(int x=0; x < stringXlen+1; x++) score[x][0] = 0;
		for(int y=0; y < stringYlen+1; y++) score[0][y] = 0;	
		// Inductive steps.
		int xMax, yMax; xMax = yMax = 0;
		double maxScore = 0; double oneScore; 
		for(int x=1; x < stringXlen+1; x++){
			for(int y=1; y < stringYlen+1; y++){
				if(0 < x && 0 < y && stringX[x-1] == stringY[y-1]) 
					 oneScore = matchScore; 
				else oneScore = mismatchPenalty;
				score[x][y] = max(	0,  // This allows generation of a new path.
									score[x-1][y-1] + oneScore, 
									score[x-1][y]   + gapPenalty,
									score[x][y-1]   + gapPenalty);
				if(maxScore < score[x][y]){
					maxScore = score[x][y]; xMax=x; yMax=y;
				}
			}
		}
//		printScoreTable(score);
		// Restore the alignment from the maximum value until it hits the zero.		
		int[][] alignment = new int[stringXlen+stringYlen][2];
		int i = 0; int x,y;
		for(x = xMax, y = yMax; !(score[x][y] == 0); i++){ // Changed to hit the zero.
			if(0 < x && 0 < y && stringX[x-1] == stringY[y-1]) 
				 oneScore = matchScore; 
			else oneScore = mismatchPenalty;
			if(0 < x && 0 < y && score[x][y] == score[x-1][y-1] + oneScore){
				alignment[i][0] = stringX[x-1]; alignment[i][1] = stringY[y-1];
				x = x-1; y = y-1;
			}else if(0 < x && 0 <= y && score[x][y] == score[x-1][y] + gapPenalty){
				alignment[i][0] = stringX[x-1]; alignment[i][1] = -1; x = x-1; 
			}else if(0 <= x && 0 < y && score[x][y] == score[x][y-1] + gapPenalty){
				alignment[i][0] = -1;  alignment[i][1] = stringY[y-1]; y = y-1;
			}
		}
		printAlignment(alignment, i-1);
	}

	public static void overlapAlignment(int[] stringX, int[] stringY){
		// Initialize constants.
		int stringXlen = stringX.length;
		int stringYlen = stringY.length; 
		double[][] score = new double[stringXlen+1][stringYlen+1];
		// Initilize the score array.
		for(int x=0; x < stringXlen+1; x++) score[x][0] = 0;
		for(int y=0; y < stringYlen+1; y++) score[0][y] = 0;
		// Inductive steps.
		double oneScore; 
		for(int x=1; x < stringXlen+1; x++){
			for(int y=1; y < stringYlen+1; y++){
				if(stringX[x-1] == stringY[y-1]) 
					 oneScore = matchScore; 
				else oneScore = mismatchPenalty;
				score[x][y] = max(	score[x-1][y-1] + oneScore, 
									score[x-1][y]   + gapPenalty,
									score[x][y-1]   + gapPenalty);
			}
		}	
//		printScoreTable(score);	
	     // The traceback procesure for restoring the best overlap alignment. 
		int x,y; int xMax, yMax; xMax = yMax = 0; double maxScore = 0; 		
		for(x=0; x<stringXlen+1; x++)  // Search the bottom row.
			if(maxScore < score[x][stringYlen]){ 
				maxScore = score[x][stringYlen]; xMax=x; yMax=stringYlen; }
		for(y=0; y<stringYlen+1; y++)  // Search the rightmost column.
			if(maxScore < score[stringXlen][y]){ 
				maxScore = score[stringXlen][y]; xMax=stringXlen; yMax=y; }	
		// Store the alignment of the bottom row or the rightmost column.
		int[][] alignment = new int[stringXlen+stringYlen][2];
		int i = 0;
		for(x=stringXlen, y=stringYlen; !(x == xMax && y == yMax); i++){  
			if(y == yMax){ alignment[i][0] = stringX[x-1]; alignment[i][1] = -1; x = x-1; }
			else{	       alignment[i][0] = -1; alignment[i][1] = stringY[y-1]; y = y-1; }			
		}
		//	Until hitting the top row or leftmost column.
		for(x = xMax, y = yMax; !(x == 0 || y == 0); i++){ 
			if(0 < x && 0 < y && stringX[x-1] == stringY[y-1]) 
				 oneScore = matchScore; 
			else oneScore = mismatchPenalty;
			if(0 <= x && 0 < y && score[x][y] == score[x][y-1] + gapPenalty){
				alignment[i][0] = -1;  alignment[i][1] = stringY[y-1]; y = y-1;
			}else if(0 < x && 0 < y && score[x][y] == score[x-1][y-1] + oneScore){
				alignment[i][0] = stringX[x-1]; alignment[i][1] = stringY[y-1];
				x = x-1; y = y-1;
			}else if(0 < x && 0 <= y && score[x][y] == score[x-1][y] + gapPenalty){
				alignment[i][0] = stringX[x-1]; alignment[i][1] = -1; x = x-1; 
			}
		}
		// Store the alignment of the top row or the leftmost column.
		for( ; 0 < x || 0 < y; i++){  // Store remaining characters in stringX.
			if(0 < x){ alignment[i][0] = stringX[x-1]; alignment[i][1] = -1; x = x-1; }
			else{	   alignment[i][0] = -1; alignment[i][1] = stringY[y-1]; y = y-1; }			
		}		
		printAlignment(alignment, i-1);
	}	

	public static void Gotoh(int[] stringX, int[] stringY){
		// Initialize constants.
		int stringXlen = stringX.length;
		int stringYlen = stringY.length; 
		double lowestScoreSum = Math.min(startupGapPenalty, 
				Math.min(mismatchPenalty, gapPenalty))*(stringXlen + stringYlen);
		double[][] diagonal   = new double[stringXlen+1][stringYlen+1];
		double[][] vertical   = new double[stringXlen+1][stringYlen+1];
		double[][] horizontal = new double[stringXlen+1][stringYlen+1];

		// Initilize the score array.
		for(int x=0; x < stringXlen+1; x++){ 
			diagonal[x][0]   = lowestScoreSum;
			horizontal[x][0] = x*gapPenalty + startupGapPenalty;
			vertical[x][0]   = lowestScoreSum;
		}
		for(int y=0; y < stringYlen+1; y++){
			diagonal[0][y]   = lowestScoreSum;			
			horizontal[0][y] = lowestScoreSum;
			vertical[0][y]   = y*gapPenalty + startupGapPenalty;
		}
		diagonal[0][0] = 0;
		
		// Inductive steps.
		double oneScore; 
		for(int x=1; x < stringXlen+1; x++){
			for(int y=1; y < stringYlen+1; y++){
				if(stringX[x-1] == stringY[y-1]) 
					 oneScore = matchScore; 
				else oneScore = mismatchPenalty;
				diagonal[x][y] = oneScore 
					+ max(	diagonal[x-1][y-1],
							horizontal[x-1][y-1], 
							vertical[x-1][y-1] );	
				horizontal[x][y] = gapPenalty 
					+ max( 	diagonal[x-1][y] + startupGapPenalty,
							horizontal[x-1][y], 
							vertical[x-1][y] + startupGapPenalty);
				vertical[x][y] = gapPenalty  
					+ max( 	diagonal[x][y-1]   + startupGapPenalty,
							horizontal[x][y-1] + startupGapPenalty, 
							vertical[x][y-1]);					
			}
		}
		
		// Restore the alignment from the bottom-right to the top-left.
		int x = stringXlen; int y = stringYlen;
		char traceType;
		double trace = max(diagonal[x][y], horizontal[x][y], vertical[x][y]);
		// Move in the vertical or horizontal direction if there are multiple choices.
		if(trace == vertical[x][y]) 	    traceType = 'v';
		else if(trace == horizontal[x][y])	traceType = 'h';
		else								traceType = 'd';		
		int[][] alignment = new int[stringXlen+stringYlen][2];
		int i = 0; 		
		for( ; !(x == 0 && y == 0); i++){
			// No need to perform an array boundary check when x-1 and y-1 are accessed.
			switch(traceType){ 
			case 'd':
				if(stringX[x-1] == stringY[y-1]) 
					 oneScore = matchScore; 
				else oneScore = mismatchPenalty;
				if(diagonal[x][y] == oneScore + diagonal[x-1][y-1])    		
					 traceType = 'd';
				else if(diagonal[x][y] == oneScore + horizontal[x-1][y-1])	
					 traceType = 'h';
				else traceType = 'v';
				alignment[i][0] = stringX[x-1]; alignment[i][1] = stringY[y-1]; 
				x = x-1; y = y-1; 
				break;
			case 'h':
				if(horizontal[x][y] == gapPenalty + diagonal[x-1][y] + startupGapPenalty)	
					 traceType = 'd';
				else if(horizontal[x][y] == gapPenalty + horizontal[x-1][y]) 				
					 traceType = 'h';
				else traceType = 'v';
				alignment[i][0] = stringX[x-1]; alignment[i][1] = -1; x = x-1; 
				break;
			case 'v':
				if(vertical[x][y] == gapPenalty + diagonal[x][y-1] + startupGapPenalty)    		
					 traceType = 'd';
				else if(vertical[x][y] == gapPenalty + horizontal[x][y-1] + startupGapPenalty) 	
					 traceType = 'h';
				else traceType = 'v';
				alignment[i][0] = -1; alignment[i][1] = stringY[y-1]; y = y-1; 
				break;
			}
		}
		printAlignment(alignment, i-1);
	}
	
	
	// Subroutines
	public static double max(double x, double y, double z){
		return Math.max(Math.max(x,y),z);
	}
	public static double max(double x, double y, double z, double u){
		return Math.max(Math.max(Math.max(x,y),z), u);
	}	
	public static void printScoreTable(double[][] score){
		int Xlen = score.length;
		int Ylen = score[0].length;
		for(int y=0; y < Ylen; y++){
			for(int x=0; x < Xlen; x++)
				System.out.print(score[x][y]+"\t");
			System.out.println();
		}
	}
	public static int[] characters2IntArray(String targetString){
		int targetLen = targetString.length();
		int[] targetIntArray = new int[targetLen];
		for(int i = 0; i < targetLen; i++)
			targetIntArray[i] = (int)targetString.charAt(i);
		return targetIntArray;
	}
	public static char int2char( int i){
		if(i == -1) return '-'; else return (char)i;
	}
	public static int[] generateRandomNucleotides(int targetLen){
		// Generate an array in which elements are selected from A,C,G,T at random.
		String nucleotides = "ACGT";
		int[] target = new int[targetLen];
		for(int i=0; i<targetLen; i++) target[i] = nucleotides.charAt((int)(Math.random()*4));
		return target;
	}
	public static void printAlignment(int[][] a, int lastPos){
		for(int j=lastPos; 0<=j; j--) System.out.print(int2char(a[j][0])+" ");
		System.out.println();
		for(int j=lastPos; 0<=j; j--)
			if(a[j][0] == -1 || a[j][1] == -1) System.out.print("  ");
			else if(a[j][0] == a[j][1]) System.out.print("| ");
			else System.out.print("x ");
		System.out.println();
		for(int j=lastPos; 0<=j; j--) System.out.print(int2char(a[j][1])+" ");	
		System.out.println();
	}
	
	public static void main( String args[] ){
		int[] stringX, stringY;
	
		// Example in the textbook.
		stringX = characters2IntArray("ATCGAT");		
		stringY = characters2IntArray("ATGCGT");		
		System.out.println("\nNeedleman Wunsch");
		NeedlemanWunsch(stringX, stringY);
		System.out.println("\nNeedleman Wunsch Hirschberg");	
		NeedlemanWunschHirschberg(stringX, stringY);	

		// Handling random sequences using the Needleman-Wunsch algorithm.
		for(int i=0; i< 3; i++){
			stringX = generateRandomNucleotides(30);
			stringY = generateRandomNucleotides(20);	
			System.out.println("\nNeedleman Wunsch");
			NeedlemanWunsch(stringX, stringY);
			System.out.println("\nNeedleman Wunsch Hirschberg");	
			NeedlemanWunschHirschberg(stringX, stringY);
		}

		// Example in the textbook.
		// Comparison between Needleman-Wunsch and Smith-Waterman.
		stringX = characters2IntArray("GATGCGATTA");	
		stringY = characters2IntArray("AGGATG");	
		System.out.println("\nNeedleman Wunsch");
		NeedlemanWunsch(stringX, stringY);
		System.out.println("\nSmith Waterman");
		SmithWaterman(stringX, stringY);
		
		// Another example of comparing Needleman-Wunsch and Smith-Waterman.	
		stringX = characters2IntArray("ATGCGATTAG");	
		stringY = characters2IntArray("TCAACGTTCT");	
		System.out.println("\nNeedleman Wunsch");
		NeedlemanWunsch(stringX, stringY);
		System.out.println("\nSmith Waterman");
		SmithWaterman(stringX, stringY);
	
		// Example in the textbook.
		// Comparison between overlap alignment, Needleman-Wunsch, and Smith-Waterman.
		stringX = characters2IntArray("ATGCGATACGATGCTATGCTAG");	
		stringY = characters2IntArray("GATGCTAGCTAGCATGCCCGG");				
		System.out.println("\nOverlap alignment");
		overlapAlignment(stringX, stringY);
		System.out.println("\nOverlap alignment (The two arguments are exchanged.)");
		overlapAlignment(stringY, stringX);	// Arguments are exchanged.	
		System.out.println("\nNeedleman Wunsch");
		NeedlemanWunsch(stringX, stringY);
		System.out.println("\nSmith Waterman");
		SmithWaterman(stringX, stringY);

		// Example in the textbook.
		// Comparison between Needleman-Wunsch and Gotoh.		
		stringX = characters2IntArray("AAGCATTGGTCATCCATGCT");	
		stringY = characters2IntArray("ATTGCCAT");			
		System.out.println("\nNeedleman Wunsch");
		NeedlemanWunsch(stringX, stringY);
		System.out.println("\nGotoh");
		Gotoh(stringX, stringY);
		System.out.println("\nGotoh@(The two arguments are exchanged.)");
		Gotoh(stringY, stringX); 	// Arguments are exchanged.	

	}	
}