content: package com.actelion.research.chem.descriptor.pharmacophoretree;
  
  
  	/*
  	 * Created on Apr 25, 2005
  	 * 
  	 * Munkres-Kuhn (Hungarian) Algorithm Clean Version: 0.11
  	 * 
  	 * Konstantinos A. Nedas                     
  	 * Department of Spatial Information Science & Engineering
  	 * University of Maine, Orono, ME 04469-5711, USA
  	 * kostas@spatial.maine.edu
  	 * http://www.spatial.maine.edu/~kostas       
  	 *
  	 * This Java class implements the Hungarian algorithm [a.k.a Munkres' algorithm,
  	 * a.k.a. Kuhn algorithm, a.k.a. Assignment problem, a.k.a. Marriage problem,
  	 * a.k.a. Maximum Weighted Maximum Cardinality Bipartite Matching].
  	 *
  	 * [It can be used as a method call from within any main (or other function).]
  	 * It takes 2 arguments:
  	 * a. A 2-D array (could be rectangular or square).
  	 * b. A string ("min" or "max") specifying whether you want the min or max assignment.
  	 * [It returns an assignment matrix[array.length][2] that contains the row and col of
  	 * the elements (in the original inputted array) that make up the optimum assignment.]
  	 *  
  	 * [This version contains only scarce comments. If you want to understand the 
  	 * inner workings of the algorithm, get the tutorial version of the algorithm
  	 * from the same website you got this one (http://www.spatial.maine.edu/~kostas/dev/soft/munkres.htm)]
  	 * 
  	 * Any comments, corrections, or additions would be much appreciated. 
  	 * Credit due to professor Bob Pilgrim for providing an online copy of the
  	 * pseudocode for this algorithm (http://216.249.163.93/bob.pilgrim/445/munkres.html)
  	 * 
  	 * Feel free to redistribute this source code, as long as this header--with
  	 * the exception of sections in brackets--remains as part of the file.
  	 * 
  	 * Requirements: JDK 1.5.0_01 or better.
  	 * [Created in Eclipse 3.1M6 (www.eclipse.org).]
  	 * 
  	 */
  
  	import static java.lang.Math.floor;
  	import static java.lang.Math.round;
  
  	import java.util.Random;
  	import java.util.Scanner;
  
  	public class HungarianAlgorithm {
  
  		//********************************//
  		//METHODS FOR CONSOLE INPUT-OUTPUT//
  		//********************************//
  		
  		public static int readInput(String prompt)	//Reads input,returns double.
  		{
  			Scanner in = new Scanner(System.in);
  			System.out.print(prompt);
  			int input = in.nextInt();
  			return input;
  		}
  		public static void printTime(double time)	//Formats time output.
  		{
  			String timeElapsed = "";
  			int days = (int)floor(time)/(24 * 3600);
  			int hours = (int)floor(time%(24*3600))/(3600);
  			int minutes = (int)floor((time%3600)/60);
  			int seconds = (int)round(time%60);
  			
  			if (days > 0)
  				timeElapsed = Integer.toString(days) + "d:";
  			if (hours > 0)
  				timeElapsed = timeElapsed + Integer.toString(hours) + "h:";
  			if (minutes > 0)
  				timeElapsed = timeElapsed + Integer.toString(minutes) + "m:";
  			
  			timeElapsed = timeElapsed + Integer.toString(seconds) + "s";
  			System.out.print("\nTotal time required: " + timeElapsed + "\n\n");
  		}
  		
  		//*******************************************//
  		//METHODS THAT PERFORM ARRAY-PROCESSING TASKS//
  		//*******************************************//
  		
  		public static void generateRandomArray	//Generates random 2-D array.
  		(double[][] array, String randomMethod)	
  		{
  			Random generator = new Random();
  			for (int i=0; i<array.length; i++)
  			{
  				for (int j=0; j<array[i].length; j++)
  				{
  					if (randomMethod.equals("random"))
  						{array[i][j] = generator.nextDouble();}
  					if (randomMethod.equals("gaussian"))
  					{
  							array[i][j] = generator.nextGaussian()/4;		//range length to 1.
  							if (array[i][j] > 0.5) {array[i][j] = 0.5;}		//eliminate outliers.
  							if (array[i][j] < -0.5) {array[i][j] = -0.5;}	//eliminate outliers.
  							array[i][j] = array[i][j] + 0.5;				//make elements positive.
  					}
  				}
  			}			
  		}
  		public static double findLargest		//Finds the largest element in a positive array.
  		(double[][] array)
  		//works for arrays where all values are >= 0.
  		{
  			double largest = 0;
  			for (int i=0; i<array.length; i++)
  			{
  				for (int j=0; j<array[i].length; j++)
  				{
  					if (array[i][j] > largest)
  					{
  						largest = array[i][j];
  					}
  				}
  			}
  				
  			return largest;
  		}
  		public static double[][] transpose		//Transposes a double[][] array.
  		(double[][] array)	
  		{
  			double[][] transposedArray = new double[array[0].length][array.length];
  			for (int i=0; i<transposedArray.length; i++)
  			{
  				for (int j=0; j<transposedArray[i].length; j++)
  				{transposedArray[i][j] = array[j][i];}
  			}
  			return transposedArray;
  		}
  		public static double[][] copyOf			//Copies all elements of an array to a new array.
  		(double[][] original)	
  		{
  			double[][] copy = new double[original.length][original[0].length];
  			for (int i=0; i<original.length; i++)
  			{
  				//Need to do it this way, otherwise it copies only memory location
  				System.arraycopy(original[i], 0, copy[i], 0, original[i].length);
  			}
  			
  			return copy;
  		}
  		
  		//**********************************//
  		//METHODS OF THE HUNGARIAN ALGORITHM//
  		//**********************************//
  		
  		public static int[][] hgAlgorithm (double[][] array, String sumType)
  		{
  			double[][] cost = copyOf(array);	//Create the cost matrix
  			
  			if (sumType.equalsIgnoreCase("max"))	//Then array is weight array. Must change to cost.
  			{
  				double maxWeight = findLargest(cost);
  				for (int i=0; i<cost.length; i++)		//Generate cost by subtracting.
  				{
  					for (int j=0; j<cost[i].length; j++)
  					{
  						cost [i][j] = (maxWeight - cost [i][j]);
  					}
  				}
  			}
  			double maxCost = findLargest(cost);		//Find largest cost matrix element (needed for step 6).
  			
  			int[][] mask = new int[cost.length][cost[0].length];	//The mask array.
  			int[] rowCover = new int[cost.length];					//The row covering vector.
  			int[] colCover = new int[cost[0].length];				//The column covering vector.
  			int[] zero_RC = new int[2];								//Position of last zero from Step 4.
  			int step = 1;											
  			boolean done = false;
  			while (done == false)	//main execution loop
  			{ 
  				switch (step)
  			    {
  					case 1:
  						step = hg_step1(step, cost);     
  			    	    break;
  			    	case 2:
  			    	    step = hg_step2(step, cost, mask, rowCover, colCover);
  						break;
  			    	case 3:
  			    	    step = hg_step3(step, mask, colCover);
  						break;
  			    	case 4:
  			    	    step = hg_step4(step, cost, mask, rowCover, colCover, zero_RC);
  						break;
  			    	case 5:
  						step = hg_step5(step, mask, rowCover, colCover, zero_RC);
  						break;
  			    	case 6:
  			    	   	step = hg_step6(step, cost, rowCover, colCover, maxCost);
  						break;
  			  	    case 7:
  			    	    done=true;
  			    	    break;
  			    }
  			}//end while
  			
  			int[][] assignment = new int[array.length][2];	//Create the returned array.
  			for (int i=0; i<mask.length; i++)
  			{
  				for (int j=0; j<mask[i].length; j++)
  				{
  					if (mask[i][j] == 1)
  					{
  						assignment[i][0] = i;
  						assignment[i][1] = j;
  					}
  				}
  			}
  			
  			//If you want to return the min or max sum, in your own main method
  			//instead of the assignment array, then use the following code:
  			/*
  			double sum = 0; 
  			for (int i=0; i<assignment.length; i++)
  			{
  				sum = sum + array[assignment[i][0]][assignment[i][1]];
  			}
  			return sum;
  			*/
  			//Of course you must also change the header of the method to:
  			//public static double hgAlgorithm (double[][] array, String sumType)
  			
  			return assignment;
  		}
  		public static int hg_step1(int step, double[][] cost)
  		{
  			//What STEP 1 does:
  			//For each row of the cost matrix, find the smallest element
  			//and subtract it from from every other element in its row. 
  		    
  		   	double minval;
  		   	
  			for (int i=0; i<cost.length; i++)	
  		   	{									
  		   	    minval=cost[i][0];
  		   	    for (int j=0; j<cost[i].length; j++)	//1st inner loop finds min val in row.
  		   	    {
  		   	        if (minval>cost[i][j])
  		   	        {
  		   	            minval=cost[i][j];
  		   	        }
  				}
  				for (int j=0; j<cost[i].length; j++)	//2nd inner loop subtracts it.
  		   	    {
  		   	        cost[i][j]=cost[i][j]-minval;
  		   	    }
  			}
  		   			    
  			step=2;
  			return step;
  		}
  		public static int hg_step2(int step, double[][] cost, int[][] mask, int[]rowCover, int[] colCover)
  		{
  			//What STEP 2 does:
  			//Marks uncovered zeros as starred and covers their row and column.
  			
  			for (int i=0; i<cost.length; i++)
  		    {
  		        for (int j=0; j<cost[i].length; j++)
  		        {
  		            if ((cost[i][j]==0) && (colCover[j]==0) && (rowCover[i]==0))
  		            {
  		                mask[i][j]=1;
  						colCover[j]=1;
  		                rowCover[i]=1;
  					}
  		        }
  		    }
  								
  			clearCovers(rowCover, colCover);	//Reset cover vectors.
  				    
  			step=3;
  			return step;
  		}
  		public static int hg_step3(int step, int[][] mask, int[] colCover)
  		{
  			//What STEP 3 does:
  			//Cover columns of starred zeros. Check if all columns are covered.
  			
  			for (int i=0; i<mask.length; i++)	//Cover columns of starred zeros.
  		    {
  		        for (int j=0; j<mask[i].length; j++)
  		        {
  		            if (mask[i][j] == 1)
  		            {
  		                colCover[j]=1;
  					}
  		        }
  		    }
  		    
  			int count=0;						
  			for (int j=0; j<colCover.length; j++)	//Check if all columns are covered.
  		    {
  		        count=count+colCover[j];
  		    }
  			
  			if (count>=mask.length)	//Should be cost.length but ok, because mask has same dimensions.	
  		    {
  				step=7;
  			}
  		    else
  			{
  				step=4;
  			}
  		    	
  			return step;
  		}
  		public static int hg_step4(int step, double[][] cost, int[][] mask, int[] rowCover, int[] colCover, int[] zero_RC)
  		{
  			//What STEP 4 does:
  			//Find an uncovered zero in cost and prime it (if none go to step 6). Check for star in same row:
  			//if yes, cover the row and uncover the star's column. Repeat until no uncovered zeros are left
  			//and go to step 6. If not, save location of primed zero and go to step 5.
  			
  			int[] row_col = new int[2];	//Holds row and col of uncovered zero.
  			boolean done = false;
  			while (done == false)
  			{
  				row_col = findUncoveredZero(row_col, cost, rowCover, colCover);
  				if (row_col[0] == -1)
  				{
  					done = true;
  					step = 6;
  				}
  				else
  				{
  					mask[row_col[0]][row_col[1]] = 2;	//Prime the found uncovered zero.
  					
  					boolean starInRow = false;
  					for (int j=0; j<mask[row_col[0]].length; j++)
  					{
  						if (mask[row_col[0]][j]==1)		//If there is a star in the same row...
  						{
  							starInRow = true;
  							row_col[1] = j;		//remember its column.
  						}
  					}
  								
  					if (starInRow==true)	
  					{
  						rowCover[row_col[0]] = 1;	//Cover the star's row.
  						colCover[row_col[1]] = 0;	//Uncover its column.
  					}
  					else
  					{
  						zero_RC[0] = row_col[0];	//Save row of primed zero.
  						zero_RC[1] = row_col[1];	//Save column of primed zero.
  						done = true;
  						step = 5;
  					}
  				}
  			}
  			
  			return step;
  		}
  		public static int[] findUncoveredZero	//Aux 1 for hg_step4.
  		(int[] row_col, double[][] cost, int[] rowCover, int[] colCover)
  		{
  			row_col[0] = -1;	//Just a check value. Not a real index.
  			row_col[1] = 0;
  			
  			int i = 0; boolean done = false;
  			while (done == false)
  			{
  				int j = 0;
  				while (j < cost[i].length)
  				{
  					if (cost[i][j]==0 && rowCover[i]==0 && colCover[j]==0)
  					{
  						row_col[0] = i;
  						row_col[1] = j;
  						done = true;
  					}
  					j = j+1;
  				}//end inner while
  				i=i+1;
  				if (i >= cost.length)
  				{
  					done = true;
  				}
  			}//end outer while
  			
  			return row_col;
  		}
  		public static int hg_step5(int step, int[][] mask, int[] rowCover, int[] colCover, int[] zero_RC)
  		{
  			//What STEP 5 does:	
  			//Construct series of alternating primes and stars. Start with prime from step 4.
  			//Take star in the same column. Next take prime in the same row as the star. Finish
  			//at a prime with no star in its column. Unstar all stars and star the primes of the
  			//series. Erasy any other primes. Reset covers. Go to step 3.
  			
  			int count = 0;												//Counts rows of the path matrix.
  			int[][] path = new int[(mask[0].length*mask.length)][2];	//Path matrix (stores row and col).
  			path[count][0] = zero_RC[0];								//Row of last prime.
  			path[count][1] = zero_RC[1];								//Column of last prime.
  			
  			boolean done = false;
  			while (done == false)
  			{ 
  				int r = findStarInCol(mask, path[count][1]);
  				if (r>=0)
  				{
  					count = count+1;
  					path[count][0] = r;					//Row of starred zero.
  					path[count][1] = path[count-1][1];	//Column of starred zero.
  				}
  				else
  				{
  					done = true;
  				}
  				
  				if (done == false)
  				{
  					int c = findPrimeInRow(mask, path[count][0]);
  					count = count+1;
  					path[count][0] = path [count-1][0];	//Row of primed zero.
  					path[count][1] = c;					//Col of primed zero.
  				}
  			}//end while
  			
  			convertPath(mask, path, count);
  			clearCovers(rowCover, colCover);
  			erasePrimes(mask);
  			
  			step = 3;
  			return step;
  			
  		}
  		public static int findStarInCol			//Aux 1 for hg_step5.
  		(int[][] mask, int col)
  		{
  			int r=-1;	//Again this is a check value.
  			for (int i=0; i<mask.length; i++)
  			{
  				if (mask[i][col]==1)
  				{
  					r = i;
  				}
  			}
  					
  			return r;
  		}
  		public static int findPrimeInRow		//Aux 2 for hg_step5.
  		(int[][] mask, int row)
  		{
  			int c = -1;
  			for (int j=0; j<mask[row].length; j++)
  			{
  				if (mask[row][j]==2)
  				{
  					c = j;
  				}
  			}
  			
  			return c;
  		}
  		public static void convertPath			//Aux 3 for hg_step5.
  		(int[][] mask, int[][] path, int count)
  		{
  			for (int i=0; i<=count; i++)
  			{
  				if (mask[(path[i][0])][(path[i][1])]==1)
  				{
  					mask[(path[i][0])][(path[i][1])] = 0;
  				}
  				else
  				{
  					mask[(path[i][0])][(path[i][1])] = 1;
  				}
  			}
  		}
  		public static void erasePrimes			//Aux 4 for hg_step5.
  		(int[][] mask)
  		{
  			for (int i=0; i<mask.length; i++)
  			{
  				for (int j=0; j<mask[i].length; j++)
  				{
  					if (mask[i][j]==2)
  					{
  						mask[i][j] = 0;
  					}
  				}
  			}
  		}
  		public static void clearCovers			//Aux 5 for hg_step5 (and not only).
  		(int[] rowCover, int[] colCover)
  		{
  			for (int i=0; i<rowCover.length; i++)
  			{
  				rowCover[i] = 0;
  			}
  			for (int j=0; j<colCover.length; j++)
  			{
  				colCover[j] = 0;
  			}
  		}
  		public static int hg_step6(int step, double[][] cost, int[] rowCover, int[] colCover, double maxCost)
  		{
  			//What STEP 6 does:
  			//Find smallest uncovered value in cost: a. Add it to every element of covered rows
  			//b. Subtract it from every element of uncovered columns. Go to step 4.
  			
  			double minval = findSmallest(cost, rowCover, colCover, maxCost);
  			
  			for (int i=0; i<rowCover.length; i++)
  			{
  				for (int j=0; j<colCover.length; j++)
  				{
  					if (rowCover[i]==1)
  					{
  						cost[i][j] = cost[i][j] + minval;
  					}
  					if (colCover[j]==0)
  					{
  						cost[i][j] = cost[i][j] - minval;
  					}
  				}
  			}
  				
  			step = 4;
  			return step;
  		}
  		public static double findSmallest		//Aux 1 for hg_step6.
  		(double[][] cost, int[] rowCover, int[] colCover, double maxCost)
  		{
  			double minval = maxCost;				//There cannot be a larger cost than this.
  			for (int i=0; i<cost.length; i++)		//Now find the smallest uncovered value.
  			{
  				for (int j=0; j<cost[i].length; j++)
  				{
  					if (rowCover[i]==0 && colCover[j]==0 && (minval > cost[i][j]))
  					{
  						minval = cost[i][j];
  					}
  				}
  			}
  			
  			return minval;
  		}
  		
  		//***********//
  		//MAIN METHOD//
  		//***********//
  		/*
  		public static void main(String[] args) 
  		{
  			//Below enter "max" or "min" to find maximum sum or minimum sum assignment.
  			String sumType = "max";		
  					
  			/*
  			double[][] array =
  			{
  					{1.0},{1.5}
  			};
  			*/
  			
  			//<UNCOMMENT> BELOW AND COMMENT BLOCK ABOVE TO USE A RANDOMLY GENERATED MATRIX
  			//int numOfRows = readInput("How many rows for the matrix? ");
  			//int numOfCols = readInput("How many columns for the matrix? ");
  			//double[][] array = new double[numOfRows][numOfCols];
  			//generateRandomArray(array, "random");	//All elements within [0,1].
  			//</UNCOMMENT>
  			/*
  			if (array.length > array[0].length)
  			{
  				System.out.println("Array transposed (because rows>columns).\n");	//Cols must be >= Rows.
  				array = transpose(array);
  			}
  					
  			//<COMMENT> TO AVOID PRINTING THE MATRIX FOR WHICH THE ASSIGNMENT IS CALCULATED
  			System.out.println("\n(Printing out only 2 decimals)\n");
  			System.out.println("The matrix is:");
  			for (int i=0; i<array.length; i++)
  			{
  				for (int j=0; j<array[i].length; j++)
  					{System.out.printf("%.2f\t", array[i][j]);}
  				System.out.println();
  			}
  			System.out.println();
  			//</COMMENT>*/
  			/*
  			double startTime = System.nanoTime();	
  			int[][] assignment = new int[array.length][2];
  			assignment = hgAlgorithm(array, sumType);	//Call Hungarian algorithm.
  			double endTime = System.nanoTime();
  							
  			System.out.println("The winning assignment (" + sumType + " sum) is:\n");	
  			double sum = 0;
  			for (int i=0; i<assignment.length; i++)
  			{
  				//<COMMENT> to avoid printing the elements that make up the assignment
  				System.out.printf("array(%d,%d) = %.2f\n", (assignment[i][0]+1), (assignment[i][1]+1),
  						array[assignment[i][0]][assignment[i][1]]);
  				sum = sum + array[assignment[i][0]][assignment[i][1]];
  				//</COMMENT>
  			}
  			
  			System.out.printf("\nThe %s is: %.2f\n", sumType, sum);
  			printTime((endTime - startTime)/1000000000.0);
  			*/
  		//}
  	
  	}
  
  
  