/*
* Copyright 2010 Google Inc.
*
* Licensed under the Apache License, Version 2.0 (the "License"); you may not
* use this file except in compliance with the License. You may obtain a copy of
* the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
* WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the
* License for the specific language governing permissions and limitations under
* the License.
*/
package com.google.gwt.dev.util.editdistance;
/**
* Computes Levenshtein string-edit distance using the
* algorithm of Eugene Myers (see "A fast bit-vector algorithm for
* approximate string matching based on dynamic progamming",
* Journal of the ACM, 46(3): 395-415, 1999), using the reformulations
* of Heikki Hyyro (see "Explaining and Extending the Bit-parallel
* Approximate String Matching Algorithm of Myers") along with
* Gonzalo Navarro ("Faster Bit-Parallel Approximate String Matching",
* Proc. 13th Combinatorial Pattern Matching (CPM'2002)).
*/
public abstract class MyersBitParallelEditDistance
implements GeneralEditDistance, Cloneable {
/*
* The Myers algorithm is based on an adaptation of the traditional
* dynamic programming algorithm, in which a matrix of distances
* is computed column-wise (or row-wise). For strings of length
* M and N, an MxN matrix (D) is filled with distances for substrings:
* D[i,j] holds the distance for string1[0..i]=>string2[0..j].
* The matrix is initially populated with the trivial values
* D[0,j]=j and D[i,0]=i; and then expanded with the rule:
* D[i,j] = min( D[i-1,j]+1, // insertion
* D[i,j-1]+1, // deletion
* (D[i-1,j-1]
* + (string1[i]==string2[j])
* ? 0 // match
* : 1 // substitution ) )
*
* The Myers algorithm takes advantage of the observation that
* the difference from one matrix position to its higher diagonal
* is either 0 or 1, and consequently row-to-row and column-to-column
* differences are 1, 0, or -1. It organizes whole columns of
* differences into bit arrays, computing new columns with
* arithmetic on whole words. The edit distance is computed
* by tracking the column changes in the last row.
*
* The description below follows the naming from the Hyrro papers above.
* The pattern string (length M) is arranged vertically; the
* string to be matched (length N) is horizontal. Columns are
* represented by bit-arrays with (1<<j) set for row i. The
* inter-cell differences are grouped into column-variables:
* VN has bits set when the transition from row i-1 to row i is
* negative; VP has bits set when it is positive; HN and HP
* have bits set for column j-1 to j transitions that are
* negative or positive, respectively. DZ has bits set when
* the diagonal from (i-1,j-1) to (i,j) is zero.
*
* The computation proceeds as follows:
* We start at column 0, where all VP bits are 1. The distance
* is accumulated in row M, with initial value M in the matrix.
*
* For each new column, we compute the diagonal-zero (DZ) first.
* This happens under three circumstances:
*
* when the pattern matches in this position ("match" above):
* We make use of a bitmask, one for each character
* in the alphabet, that has bits set wherever that character
* appears in the pattern. These bitmasks are pre-computed
* when the pattern is established (in the class constructor).
* In the column loop, we simply look up the bitmask for the
* character in this column. We call this PM.
*
* when the prior column has a negative vertical transition
* in this row ("deletion"):
* This is simply the VN bit array from the prior column.
*
* when the row above has a negative transition in this column
* ("insertion").
* This happens when the diagonal coming into the cell above
* (in *this* column) is zero but the vertical transition
* in the preceding column for that row is positive. That is:
*
* D[i-2,j-1]
* | \
* +1 0
* | \
* D[i-1,j-1] -1 -> D[i-1,j]
* | \ |
* ? 0 +1
* | \ |
* D[i,j-1] D[i,j]
*
* That is DZ[i,j] = DZ[i-1,j] & VP[i-1,j-1].
* Note that this is recursive in the DZ[*,j] bit array, but we
* want to compute the whole array at once. We need an expression
* that will propagate a 1 bit down the DZ column when the associated
* VP bit is set. Noting that the DZ[*,0] is the least significant
* bit, this means propagating to more significant bits. The carry bit
* of an addition has this behavior. If we set
* DZ = (PM | VN) // first two cases
* then we can propagate bits up iff the corresponding VP bit
* is set with:
* DZ |= (((DZ & VP) + VP) ^ VP);
* The (DZ&VP) piece covers our required condition on VP. Adding
* the VP bits causes any 1 bits from (DZ&VP) to result in a carry,
* but also introduces a 1 bit where VP=1 but DZ=0, which are then
* cleared with the ^VP piece. The propagated bits are |ed into the
* base case to get a final result.
*
* Once we have DZ, we can compute HP and HN. Note that VP and VN have
* not yet been updated, and so refer to the preceding column; DZ refers
* to the current column.
* HN is true when DZ is true in this column but VP was true
* in the prior column (HN = DZ & VP)
* HP is true when: the prior VN is true (otherwise, the diagonal would
* be negative, and that cannot happen); or, DZ is false (meaning
* the diagonal increased by 1) and the prior-column VP is also false
* (HP = VN | ~(DZ|VP))
*
* Once we have HP and HN, we can tally the effect on the last row (adding
* or subtracting one, respectively).
*
* Finally, we can compute VP and VN for this column, in preparation
* for the next round. These computations are analagous to the ones for
* HN and HP above, except that there, the incoming V values already
* reflected the prior column and DZ this one, but now the H and DZ values
* reflect the same column. So, we first shift the H values to associate
* with the proper bits from DZ.
*/
/**
* A trivial implementation for the zero-length string.
*/
static class Empty extends MyersBitParallelEditDistance {
Empty(CharSequence s) {
super(s);
}
@Override
public GeneralEditDistance duplicate() {
return this; /* thread-safe */
}
@Override
public int getDistance(CharSequence s, int k) {
return s.length();
}
}
/**
* An implementation using multiple integers for each column.
* This follows the pattern for the single-word implementation, with
* these changes:
* sums/shifts must propagate from one word to the next;
* to make that easier (and cheaper), we use only (SIZE-1) bits per word;
*/
static class Multi extends MyersBitParallelEditDistance {
/* How many integers we use per column */
int count;
/**
* Where the last-row bit lives -- only in the last array slot, though.
*/
final int lastBitPosition;
/**
* Bitmaps for pattern string: [pattern_index][column_index].
*/
final int[][] positions;
int[] verticalNegativesReusable;
/* Reusable arrays for vertical changes */
int[] verticalPositivesReusable;
/**
* A mask with those bits set.
*/
final int wordMask = (-1 >>> 1);
/**
* How many bits we use per word.
*/
final int wordSize = Integer.SIZE - 1;
/**
* Constructs a multi-word engine.
*/
Multi(CharSequence s) {
super(s);
/* Compute number of words to use (rounding up) */
count = (m + wordSize - 1) / wordSize;
/* Precompute bitmaps for pattern string */
positions = PatternBitmap.map(s, idx, new int [idx.size()][], wordSize);
/* Where in last word does last row appear */
lastBitPosition = (1 << ((m - 1) % wordSize));
/* Initialize scratchpad items */
perThreadInit();
}
@Override
public int getDistance(CharSequence s, int k) {
indices = idx.map(s, indices);
/* Initialize verticalPositive to all bits on, verticalNegative all off */
int[] verticalPositives = verticalPositivesReusable;
java.util.Arrays.fill(verticalPositives, wordMask);
int[] verticalNegatives = verticalNegativesReusable;
java.util.Arrays.fill(verticalNegatives, 0);
int distance = m;
int len = s.length();
/* We can only miss the distance-- below this many times: */
int maxMisses = k + len - m;
if (maxMisses < 0) {
maxMisses = Integer.MAX_VALUE;
}
outer:
for (int j = 0; j < len; j++) {
int[] position = positions[indices[j]];
/* Carry bits from one word to the next */
int sum = 0;
int horizontalPositiveShift = 1;
int horizontalNegativeShift = 0;
/* Iterate through words for this column */
for (int i = 0; i < count; i++) {
int verticalNegative = verticalNegatives[i];
int patternMatch = (position[i] | verticalNegative);
int verticalPositive = verticalPositives[i];
sum = (verticalPositive & patternMatch)
+ (verticalPositive) + (sum >>> wordSize);
int diagonalZero = ((sum & wordMask) ^ verticalPositive)
| patternMatch;
int horizontalPositive = (verticalNegative
| ~(diagonalZero | verticalPositive));
int horizontalNegative = diagonalZero & verticalPositive;
if (i == (count - 1)) { /* only last bit in last word */
if ((horizontalNegative & lastBitPosition) != 0) {
distance--;
} else if ((horizontalPositive & lastBitPosition) != 0) {
distance++;
if ((maxMisses -= 2) < 0) {
break outer;
}
} else if (--maxMisses < 0) {
break outer;
}
}
horizontalPositive = ((horizontalPositive << 1)
| horizontalPositiveShift);
horizontalPositiveShift = (horizontalPositive >>> wordSize);
horizontalNegative = ((horizontalNegative << 1)
| horizontalNegativeShift);
horizontalNegativeShift = (horizontalNegative >>> wordSize);
verticalPositives[i] = (horizontalNegative
| ~(diagonalZero | horizontalPositive))
& wordMask;
verticalNegatives[i] = (diagonalZero & horizontalPositive) & wordMask;
}
}
return distance;
}
@Override
protected void perThreadInit() {
super.perThreadInit();
/* Allocate verticalPositive/verticalNegative arrays */
verticalPositivesReusable = new int[count];
verticalNegativesReusable = new int[count];
}
}
/*
* The following code is duplicated with both "int" and "long"
* as the word primitive type. To improve maintainability, all instances
* of the word type are marked with "WORD" in comments. This
* allows the "long" version to be regenerated from the "int" version
* through a mechanical replacement. To make this work, the
* class name also needs to include the primitive type (without
* altered capitalization) and the WORD marker -- this is a quite
* intentional violation of the usual class naming rules.
*/
/**
* An implementation using "int" as the word size.
*/
// Replace "int/*WORD*/" with "long/*WORD*/" to generate the long version.
static class TYPEint/*WORD*/ extends MyersBitParallelEditDistance {
final int/*WORD*/ lastBitPosition;
final int/*WORD*/[] map;
@SuppressWarnings("cast")
TYPEint/*WORD*/(CharSequence s) {
super(s);
/* Precompute bitmaps for this pattern */
map = PatternBitmap.map(s, idx, new int/*WORD*/[idx.size()]);
/* Compute the bit that represents a change in the last row */
lastBitPosition = (((int/*WORD*/) 1) << (m - 1));
}
@Override
public int getDistance(CharSequence s, int k) {
int len = s.length();
/* Quick check based on length */
if (((len - m) > k) || ((m - len) > k)) {
return k + 1;
}
/* Map characters to their integer positions in the bitmap array */
indices = idx.map(s, indices);
/* Initially, vertical change is all positive (none negative) */
int/*WORD*/ verticalPositive = -1;
int/*WORD*/ verticalNegative = 0;
int distance = m;
/* We can only miss the "distance--" below this many times: */
int maxMisses = k + len - m;
if (maxMisses < 0) {
maxMisses = Integer.MAX_VALUE;
}
for (int j = 0; j < len; j++) {
/* Where is diagonal zero: matches, or prior VN; plus recursion */
int/*WORD*/ diagonalZero = map[indices[j]] | verticalNegative;
diagonalZero |= (((diagonalZero & verticalPositive) + verticalPositive)
^ verticalPositive);
/* Compute horizontal changes */
int/*WORD*/ horizontalPositive = verticalNegative
| ~(diagonalZero | verticalPositive);
int/*WORD*/ horizontalNegative = diagonalZero & verticalPositive;
/* Update final distance based on horizontal changes */
if ((horizontalNegative & lastBitPosition) != 0) {
distance--;
} else if ((horizontalPositive & lastBitPosition) != 0) {
distance++;
if ((maxMisses -= 2) < 0) {
break;
}
} else if (--maxMisses < 0) {
break;
}
/* Shift Hs to next row, compute new Vs analagously to Hs above */
horizontalPositive = (horizontalPositive << 1) | 1;
verticalPositive = (horizontalNegative << 1)
| ~(diagonalZero | horizontalPositive);
verticalNegative = diagonalZero & horizontalPositive;
}
return distance;
}
}
/**
* An implementation using "long" as the word size.
* GENERATED MECHANICALLY FROM THE "int" VERSION ABOVE.
* DO NOT EDIT THIS ONE -- EDIT ABOVE AND REAPPLY QUERY/REPLACE.
*/
static class TYPElong/*WORD*/ extends MyersBitParallelEditDistance {
final long/*WORD*/ lastBitPosition;
final long/*WORD*/[] map;
TYPElong/*WORD*/(CharSequence s) {
super(s);
/* Precompute bitmaps for this pattern */
map = PatternBitmap.map(s, idx, new long/*WORD*/[idx.size()]);
/* Compute the bit that represents a change in the last row */
lastBitPosition = (((long/*WORD*/) 1) << (m - 1));
}
@Override
public int getDistance(CharSequence s, int k) {
int len = s.length();
/* Quick check based on length */
if (((len - m) > k) || ((m - len) > k)) {
return k + 1;
}
/* Map characters to their integer positions in the bitmap array */
indices = idx.map(s, indices);
/* Initially, vertical change is all positive (none negative) */
long/*WORD*/ verticalPositive = -1;
long/*WORD*/ verticalNegative = 0;
int distance = m;
/* We can only miss the "distance--" below this many times: */
int maxMisses = k + len - m;
if (maxMisses < 0) {
maxMisses = Integer.MAX_VALUE;
}
for (int j = 0; j < len; j++) {
/* Where is diagonal zero: matches, or prior VN; plus recursion */
long/*WORD*/ diagonalZero = map[indices[j]] | verticalNegative;
diagonalZero |= (((diagonalZero & verticalPositive) + verticalPositive)
^ verticalPositive);
/* Compute horizontal changes */
long/*WORD*/ horizontalPositive = verticalNegative
| ~(diagonalZero | verticalPositive);
long/*WORD*/ horizontalNegative = diagonalZero & verticalPositive;
/* Update final distance based on horizontal changes */
if ((horizontalNegative & lastBitPosition) != 0) {
distance--;
} else if ((horizontalPositive & lastBitPosition) != 0) {
distance++;
if ((maxMisses -= 2) < 0) {
break;
}
} else if (--maxMisses < 0) {
break;
}
/* Shift Hs to next row, compute new Vs analagously to Hs above */
horizontalPositive = (horizontalPositive << 1) | 1;
verticalPositive = (horizontalNegative << 1)
| ~(diagonalZero | horizontalPositive);
verticalNegative = diagonalZero & horizontalPositive;
}
return distance;
}
}
/**
* Chooses an appropriate implementation for a given pattern string.
* @param s pattern string
* @return distance calculator appropriate for the pattern
*/
public static MyersBitParallelEditDistance getInstance(CharSequence s) {
int m = s.length();
return (m <= Integer.SIZE) ?
((m == 0) ? new Empty(s) : new TYPEint(s)) :
(s.length() <= Long.SIZE) ?
new TYPElong(s) :
new Multi(s);
}
/**
* Tests a computation manually.
*/
public static void main(String[] args) {
MyersBitParallelEditDistance b = getInstance(args[0]);
int k = args.length > 2 ? Integer.parseInt(args[2]) : 0;
System.out.println("Result: " + b.getDistance(args[1], k));
}
/**
* Index mapping for pattern string.
*/
final CharIndex idx;
/**
* Reusable array of indices for target strings.
*/
int[] indices = new int[0];
/**
* Length of pattern.
*/
final int m;
/**
* Constructs a distance calculator for a given string.
*/
protected MyersBitParallelEditDistance(CharSequence s) {
m = s.length();
idx = CharIndex.getInstance(s);
}
public GeneralEditDistance duplicate() {
try {
return (MyersBitParallelEditDistance) clone();
} catch (CloneNotSupportedException x) { /*IMPOSSIBLE */
throw new IllegalStateException("Cloneable object would not clone");
}
}
/**
* Computes distance from the pattern to a given string, bounded by
* a limiting distance @see(GeneralEditDistance.getDistance(CharSequence,int)).
*/
public abstract int getDistance(CharSequence s, int k);
@Override
protected Object clone() throws CloneNotSupportedException {
Object obj = super.clone();
/* Re-initialize any non-thread-safe parts */
((MyersBitParallelEditDistance) obj).perThreadInit();
return obj;
}
/**
* Initializes items that cannot be shared among threads.
*/
protected void perThreadInit() {
indices = new int[0];
}
}