hungarianalgorithm.java

a molecular graph kernel based on iterative graph similarity and optimal assignments

/**
 * ISOAK - Iterative similarity optimal assignment kernel.
 * 
 * Written by Matthias Rupp 2006-2007.
 * Copyright (c) 2006-2007, Matthias Rupp, Ewgenij Proschak, Gisbert Schneider.
 * 
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 *  * The above copyright notice, this permission notice and the following disclaimers
 *    shall be included in all copies or substantial portions of this software.
 *  * The names of the authors may not be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 * 
 * Please cite
 * 
 * Matthias Rupp, Ewgenij Proschak, Gisbert Schneider:
 * A Kernel Approach to Molecular Similarity Based on Iterative Graph Similarity,
 * Journal of Chemical Information and Molecular Modeling, 47(6): 2280-2286, 2007, 
 * DOI http://dx.doi.org/10.1021/ci700274r.
 *
 * if you use this software.
 */

package info.mrupp.isoak1;

import static java.util.Arrays.fill;

/**
 * The Hungarian algorithm (Kuhn-Munkres assignment algorithm).
 * 
 * Computes an assignment of the rows to the columns of a matrix (or vice versa)
 * with minimum or maximum value. Returns the value of the assignment and 
 * optionally the indices of the assignment.
 * 
 * <ul>
 * <li>If the matrix has less rows than columns, rows are assigned to columns.
 *     If the matrix has more rows than columns, columns are assigned to rows.</li>
 * <li>If more than one optimal assignment is to be computed, create one instance
 *     of the HungarianAlgorithm class and reuse it. This will prevent repeated 
 *     reallocation of memory.</li>
 * <li>Not to be used concurrently (use one instance per thread).</li>
 * </ul>
 * 
 * Implementation based on the <a href="http://216.249.163.93/bob.pilgrim/445/munkres.html">
 * exposition by Bob Pilgrim</a>.
 * 
 * @author Matthias Rupp
 * @version 0.1, 2006-08-24
 *
 */

public class HungarianAlgorithm {
	// Internal state of the algorithm object.
	// c is an arbitrary constant outside the range of valid indices; its value changes between computations.
	int[] rowStarred;      // For each row, the column of the starred zero or c if there is none.
	int[] rowPrimed;       // For each row, the column of the primed zero or c if there is none.
	boolean[] rowCovered;  // For each row, true iff row is covered.
	int[] colStarred;      // For each column, the row of the starred zero or c if there is none.
	int[] colPrimed;       // For each column, the row of the primed zero or c if there is none.
	boolean[] colCovered;  // For each column, true iff column is covered.
	int[] pathRows;        // Row coordinates of constructed path.
	int[] pathCols;        // Column coordinates of constructed path.
	int pathSize;          // Number of valid entries in path.
	int[] assignBuf;       // Storage of assignment indices if caller does not provide it.
	XArray2Float X;        // A (possibly modified, e.g. inverted, transposed) copy of the input matrix. 

    /** Initializes the optimizer. 
     * 
     * Only minimal memory is reserved. Additional memory will be reserved depending
     * on passed problem instances. Reuse {@link HungarianAlgorithm HungarianAlgorithm} objects
     * (one per thread).
     */
	public HungarianAlgorithm() {
		rowStarred = new int[0];
		rowPrimed = new int[0];
		rowCovered = new boolean[0];
		colStarred = new int[0];
		colPrimed = new int[0];
		colCovered = new boolean[0];
		pathRows = new int[0];
		pathCols = new int[0];
		pathSize = 0;
		assignBuf = new int[0];
		X = new XArray2Float();
	}

	/** Optimization type (either minimization or maximization of the assignment value). */
	public enum OptimizationType {
        /** The assignment value is minimized. */ MIN, 
        /** The assignment value is maximized. */ MAX };
	
	/** 
	 * The Hungarian algorithm.
	 * 
	 * Computes an optimal assignment for a given matrix.
	 * 
	 * @param matrix     the matrix for which the assignment is computed.
	 * @param type       indicates wether a minimum value or maximum value assignment is computed.
	 * @param assignment a buffer where the (zero-based) indices of the computed assignment are stored. Can be null if the indices are not needed.
	 * @return           the value of the optimal assignment.
	 */
	public float hungarianAlgorithm(Array2Float matrix, OptimizationType type, int[] assignment) {
		// Shortcuts for number of rows and columns.
		final int rows = matrix.rows(); 
		final int cols = matrix.cols();
		
		// Provide memory for assignment indices if not provided by caller.
		if (assignment == null) { 
			if (assignBuf.length < Math.min(rows,cols)) assignBuf = new int[Math.min(rows,cols)]; 
			assignment = assignBuf; 
		}
		
		// Set up a suitable copy of the original input.		
		final boolean inverted = (type == OptimizationType.MAX);	// For maximization, we minimize the negative of the original matrix.
		final boolean transposed = (rows > cols);					// So we can always assign rows to columns.
	
		if (transposed) X.redim(cols,rows); else X.redim(rows,cols);	// Handles memory reallocation if necessary.
		
		if (!inverted && !transposed) for (int r = 0; r < rows; ++r) for (int c = 0; c < cols; ++c) X.set(r,c, matrix.get(r,c));
		else if (!inverted && transposed) for (int r = 0; r < rows; ++r) for (int c = 0; c < cols; ++c) X.set(c,r, matrix.get(r,c));
		else if (inverted && !transposed) for (int r = 0; r < rows; ++r) for (int c = 0; c < cols; ++c) X.set(r,c, -matrix.get(r,c));
		else /* inverted && transposed */ for (int r = 0; r < rows; ++r) for (int c = 0; c < cols; ++c) X.set(c,r, -matrix.get(r,c));
		
		// Compute the assignment via the Hungarian algorithm.
		doHungarianAlgorithm(assignment);
		
		// Compute the assignment's value.
		float result = 0.f;
		for (int i = 0; i < X.rows(); ++i) result += (transposed ? matrix.get(assignment[i],i) : matrix.get(i,assignment[i]));
	
		return result;
	}
		
	/**
	 * Hungarian algorithm (Munkres assignment algorithm).
	 * 
	 * Computes minimum value assignment for matrices with more columns than rows.
	 *
	 * - The internal matrix copy X is used and modified.
	 * 
	 * @param assignment array where the assignment indices are stored.
	 */
	private void doHungarianAlgorithm(int[] assignment) {
		// Initialize algorithm state.
		final int k = X.rows();		// Shortcut for number of rows. 
		final int n = X.cols();		// Shortcut for number of columns.
		final int c = n*k;				// Arbitrary constant indicating special conditions in data arrays.
		if (assignment != null && assignment.length < k) throw new IllegalArgumentException(String.format("Result array too short for Hungarian algorithm (expected %d but got %d).", k, assignment.length));
		rowStarred = initIntArray(rowStarred, k, c);
		rowPrimed = initIntArray(rowPrimed, k, c);
		rowCovered = initBoolArray(rowCovered, k, false);
		colStarred = initIntArray(colStarred, n, c);
		colPrimed = initIntArray(colPrimed, n, c);
		colCovered = initBoolArray(colCovered, n, false);
		pathRows = initIntArray(pathRows, 2*n, 0);
		pathCols = initIntArray(pathCols, 2*n, 0);
		pathSize = 0;
		int rrow, ccol;		// Temporary storage of row and column indices.

		// Step 1: Subtract the row minimum from each row.
		X.subtract_row_mins();

		// Step 2: Find all zeroes; star each zero which does not have a starred zero in it's row or column.
		for (int row = k; --row >= 0; ) {
			for (int col = n; --col >= 0; ) {
				if (X.get(row, col) == 0.f && rowStarred[row] == c && colStarred[col] == c) {
					rowStarred[row] = col;
					colStarred[col] = row;
					break;  // There's a starred zero in this row, so quit searching it.
				}
			}
		}

		// Execute steps 3 to 6.
		for (int step = 3; step != 7; ) {
			switch (step) {
			case 3:
		        // Cover each column containing a starred zero.
				for (int i = n; --i >= 0; ) if (colStarred[i] != c) colCovered[i] = true;

		        // If k columns have been covered, we are done.
				int nc = 0; for (int i = n; --i >= 0; ) if (colCovered[i]) ++nc;  // Count number of covered columns.
				if (nc == k) step = 7; else step = 4;
				
				break;
			
			case 4:
				while (true) {
					// Find a noncovered zero. (It is possible that there is none).
					rrow = c; ccol = c;
					for (int row = k; --row >= 0; ) {
						if (rowCovered[row]) continue;  // If the row is covered, there is no need to search inside it.
						for (int col = n; --col >= 0; ) {
							if (X.get(row, col) == 0.f && !colCovered[col]) {
								rrow = row;
								ccol = col;
								break;
							}
						}
						if (rrow != c) break;
					}
					
					// If no such zero was found, go to step 6.
					if (rrow == c) { step = 6; break; }
					
					// Prime the found zero.
					rowPrimed[rrow] = ccol;
					colPrimed[ccol] = rrow;
					
					// If there is no starred zero in this row, go to step 5.
					if (rowStarred[rrow] == c) {
						pathRows[0] = rrow;
						pathCols[0] = ccol;
						pathSize = 1;
						step = 5;
						break;
					}
					
					// Cover the row and uncover the column with the starred zero.
					rowCovered[rrow] = true;
					colCovered[rowStarred[rrow]] = false;
				}
				
				break;
				
			case 5:
		        // Construct a series of alternating primed and starred zeros.
		        // Primed zero from step 4, starred zero in column, primed zero in row, starred in column, ...
		        ccol = pathCols[0];
		        while (true) {
		            // Add starred zero in primed zeros column.
		        	if (colStarred[ccol] == c) break;  // Exit loop, no starred zero in primed zeros column.
		        	rrow = colStarred[ccol];		// Column is identical.
		        	pathRows[pathSize] = rrow;
		        	pathCols[pathSize] = ccol;
		        	++pathSize;

		            // Add primed zero in starred zeros row.
		            ccol = rowPrimed[rrow];
		            pathRows[pathSize] = rrow;
		            pathCols[pathSize] = ccol;
		            ++pathSize;
		        }

		        // In the series, unstar all starred zeros, star each primed zero. Unprime all zeros, uncover all rows and columns. Go to step 3.

		        // Unstar all zeros in the path.
		        // We do not do this simultaneously with starring the primes since this might lead to two starred zeros in a row or column.
		        for (int i = pathSize; --i >= 0; ) {
		        	// If path element is starred, unstar it.
		        	if (rowStarred[pathRows[i]] == pathCols[i] && colStarred[pathCols[i]] == pathRows[i]) {
		        		rowStarred[pathRows[i]] = c;
		        		colStarred[pathCols[i]] = c;
		        	}
		        }
		                
		        // Star each primed zero in the path.
		        for (int i = pathSize; --i >= 0; ) {
		        	// if Path element is primed, star it.
		        	if (rowPrimed[pathRows[i]] == pathCols[i] && colPrimed[pathCols[i]] == pathRows[i]) {
		        		// Do not bother with unpriming the element since all primes will be removed later on.
		                rowStarred[pathRows[i]] = pathCols[i];
		                colStarred[pathCols[i]] = pathRows[i];		        		
		        	}
		        }

		        // Unprime all zeros in the matrix.
		        fill(rowPrimed, 0, k, c);
		        fill(colPrimed, 0, n, c);

		        // Uncover all matrix entries.
		        fill(rowCovered, 0, k, false);
		        fill(colCovered, 0, n, false);

		        // Go to step 3.
		        step = 3;
		        
		        break;
				
			case 6:
		        // Find smallest uncovered value.
				float min = Float.MAX_VALUE;
				for (int row = k; --row >= 0; ) {
					if (rowCovered[row]) continue;
					for (int col = n; --col >= 0; ) {
						if (!colCovered[col]) min = Math.min(min,X.get(row,col));
					}
				}
		        
		        // Add the value to all covered rows and subtract it from all uncovered columns.
				for (int row = k; --row >= 0; )
					for (int col = n; --col >= 0; )
						if (rowCovered[row] && colCovered[col]) X.set(row, col, X.get(row,col)+min);  // Add where covered by both.
						else if (!rowCovered[row] && !colCovered[col]) X.set(row, col, X.get(row,col)-min); // Subtract where both do not cover.

		        // Go to step 4.
		        step = 4;
		        
		        break;
		        
	        default:
	        	throw new RuntimeException("Invalid step number in Hungarian algorithm.");
			}
		}
		        	
		// Compute assignment indices.
		System.arraycopy(rowStarred, 0, assignment, 0, k); 
	}

	// Helper functions.

	/** Initializes an integer array to given size and constant content. */
	private final static int[] initIntArray(int[] array, int nsize, int c) {
		if (array.length < nsize) array = new int[nsize];	// Allocate new memory only if necessary.
		fill(array, 0, nsize, c);							// Only initialize requested size part of array.
		return array;
	}
	
	/** Initializes a boolean array to given size and constant content. */
	private final static boolean[] initBoolArray(boolean[] array, int nsize, boolean c) {
		if (array.length < nsize) array = new boolean[nsize];	// Allocate new memory only if necessary.
		fill(array, 0, nsize, c);								// Only initialize requested size part of array.
		return array;
	}
	
	/** Helper class which efficiently extends Array2Float by capabilities specific to the Hungarian algorithm. */
	private class XArray2Float extends Array2Float {
		/** Subtracts the minimum of each row from it. */
		public void subtract_row_mins() {
			for (int row = rows; --row >= 0; ) {
				// Determine row minimum.
				float min = data[sub2ind(row,0)]; 
				for (int col = cols; --col > 0; ) { 
					final float value = data[sub2ind(row,col)]; 
					if (value < min) min = value; 
				}
				
				// Subtract row minimum.
				for (int col = cols; --col >= 0; ) { final int index = sub2ind(row,col); data[index] = data[index] - min; }
			}
		}
	}	
}