from collections import defaultdict
from math import isfinite
from numbers import Real
class StringEdit1:
"""Distance between strings.
Parameters
----------
insertion_cost : float
Cost of inserting a character (default 1.0).
deletion_cost : float
Cost of deleting a character (default 1.0).
substitution_cost : float | None
Cost of substituting a character. Defaults to the sum
of insertion and deletion costs.
"""
def __init__(
self,
insertion_cost: float = 1.,
deletion_cost: float = 1.,
substitution_cost: float | None = None,
) -> None:
self._insertion_cost = self._validate_cost(
insertion_cost,
'insertion_cost',
)
self._deletion_cost = self._validate_cost(
deletion_cost,
'deletion_cost',
)
if substitution_cost is None:
self._substitution_cost = (
self._insertion_cost + self._deletion_cost
)
else:
self._substitution_cost = self._validate_cost(
substitution_cost,
'substitution_cost',
)
self._call_counter: defaultdict[tuple[str, str], int] = defaultdict(int)
@staticmethod
def _validate_cost(cost: float, name: str) -> float:
if isinstance(cost, bool) or not isinstance(cost, Real):
raise TypeError(f'{name} must be a real number')
if not isfinite(cost) or cost < 0:
raise ValueError(f'{name} must be finite and nonnegative')
return float(cost)
def __call__(self, source: str, target: str) -> float:
"""Compute the edit distance between two strings.
Parameters
----------
source : str
The source string.
target : str
The target string.
Returns
-------
float
The minimum edit distance.
"""
self._call_counter = defaultdict(int)
return self._naive_levenshtein(source, target)
def _naive_levenshtein(self, source: str, target: str) -> float:
self._call_counter[(source, target)] += 1
# base case
if len(source) == 0:
return self._insertion_cost*len(target)
if len(target) == 0:
return self._deletion_cost*len(source)
# test if last characters of the strings match
if source[len(source) - 1] == target[len(target) - 1]:
sub_cost = 0.
else:
sub_cost = self._substitution_cost
# Minimum over deletion, insertion, and substitution/match.
return min(
self._naive_levenshtein(source[:-1], target)
+ self._deletion_cost,
self._naive_levenshtein(source, target[:-1])
+ self._insertion_cost,
self._naive_levenshtein(source[:-1], target[:-1]) + sub_cost,
)
@property
def call_counter(self) -> dict[tuple[str, str], int]:
"""The number of times each subproblem was computed."""
return self._call_counterLevenshtein distance
Okay, so how do we compute edit (or Levenshtein) distance (Levenshtein 1966)? The basic idea is to define it recursively, deciding at each point in the string whether we want to insert/delete an element (each at some cost \(c\)) or whether we want to try matching the string.
\[d_\text{lev}(\mathbf{a}, \mathbf{b}) = \begin{cases} c_\text{ins} \times |\mathbf{b}| & \text{if } |\mathbf{a}| = 0 \\ c_\text{del} \times |\mathbf{a}| & \text{if } |\mathbf{b}| = 0 \\ \min \begin{cases}d_\text{lev}(a_1\ldots a_{|\mathbf{a}|-1}, \mathbf{b}) + c_\text{del} \\ d_\text{lev}(\mathbf{a}, b_1\ldots b_{|\mathbf{b}|-1}) + c_\text{ins} \\ d_\text{lev}(a_1\ldots a_{|\mathbf{a}|-1}, b_1\ldots b_{|\mathbf{b}|-1}) + c_\text{sub} \times \mathbb{1}[a_{|\mathbf{a}|} \neq b_{|\mathbf{b}|}]\end{cases} & \text{otherwise}\end{cases}\]
where \(c_\text{sub}\) defaults to \(c_\text{del} + c_\text{ins}\).
The recurrence computes a minimum weighted edit cost when all three costs are nonnegative. Calling that cost a metric requires more: insertion and deletion must have the same positive cost, and substitution must have the same positive cost in both directions. Under those assumptions, every edit has an inverse of the same cost, and the shortest-path interpretation gives symmetry and the triangle inequality. If insertion and deletion costs differ, the code still computes a directed edit cost, but that cost is not a metric.
Why does this recurrence return the minimum cost rather than merely the cost of some edit sequence? We can prove the claim using the last-edit decomposition (LED).
Fix arbitrary prefixes \(a_1\ldots a_i\) and \(b_1\ldots b_j\). Represent an edit sequence as a path from \((0,0)\) to \((i,j)\) in a grid. A vertical step from \((i-1,j)\) deletes \(a_i\), a horizontal step from \((i,j-1)\) inserts \(b_j\), and a diagonal step from \((i-1,j-1)\) matches the final symbols or substitutes \(b_j\) for \(a_i\). Every edit path into \((i,j)\) must end with exactly one of these three step types. Thus, every candidate belongs to one of three exhaustive cases.
Let \(D[i,j]\) be the value assigned by the displayed recurrence to the prefixes \(a_1\ldots a_i\) and \(b_1\ldots b_j\). We show that \(D[i,j]\) equals the minimum cost among all edit paths between those prefixes. Start with the chart boundaries. If \(i=0\), the only possible path inserts all \(j\) target symbols, so \(D[0,j]=j\,c_\text{ins}\). If \(j=0\), the only possible path deletes all \(i\) source symbols, so \(D[i,0]=i\,c_\text{del}\).
Now take \(i,j>0\) and suppose that the claim holds for cells whose coordinate sum is smaller than \(i+j\). Consider an arbitrary path into \((i,j)\).
- If its final step is a deletion, its earlier part ends at \((i-1,j)\). By the inductive hypothesis, that earlier part costs at least \(D[i-1,j]\), so the whole path costs at least \(D[i-1,j]+c_\text{del}\).
- If its final step is an insertion, the same reasoning gives the lower bound \(D[i,j-1]+c_\text{ins}\).
- If its final step is diagonal, the lower bound is \(D[i-1,j-1]\) when \(a_i=b_j\) and \(D[i-1,j-1]+c_\text{sub}\) otherwise.
Since every path falls into one of these cases, no path can cost less than the minimum in the recurrence. For the other direction, choose the predecessor that attains that minimum. The inductive hypothesis supplies a path with the predecessor’s chart cost, and appending the corresponding final edit constructs a path whose cost is exactly the new chart value. The value is thus both a lower bound on all edit paths and the cost of an available edit path. It is the optimum.
Compare kat and kæt. At the cell for the full strings, the diagonal predecessor already aligns ka with kæ; adding a zero-cost match for the final t completes that alignment. In contrast, an insertion or deletion would add its positive cost. The recurrence chooses the diagonal path. Thus, LED explains both the recursive implementation below and the dynamic-programming chart later in the section.
editdist = StringEdit1(1, 1)
editdist('æbstɹækt', 'æbstɹækt'), editdist('æbstɹækt', 'æbstɹækʃən'), editdist('æbstɹækʃən', 'æbstɹækt'), editdist('æbstɹækt', '')Okay. So here’s the thing. This looks nice, but it’s actually not that efficient because we’re actually redoing a whole ton of work.
editdist('æbstɹækʃən', 'æbstɹækt')
editdist.call_counterWe could try to get around this by memoizing using the lru_cache decorator.
from functools import lru_cache
class StringEdit2(StringEdit1):
"""Distance between strings with LRU cache memoization.
Parameters
----------
insertion_cost : float
Cost of inserting a character (default 1.0).
deletion_cost : float
Cost of deleting a character (default 1.0).
substitution_cost : float | None
Cost of substituting a character.
"""
def __call__(self, source: str, target: str) -> float:
"""Compute one memoized distance with a fresh per-call cache."""
self._call_counter = defaultdict(int)
self._naive_levenshtein.cache_clear()
return self._naive_levenshtein(source, target)
@lru_cache(maxsize=None)
def _naive_levenshtein(self, source: str, target: str) -> float:
self._call_counter[(source, target)] += 1
# base case
if len(source) == 0:
return self._insertion_cost*len(target)
if len(target) == 0:
return self._deletion_cost*len(source)
# test if last characters of the strings match
if source[len(source) - 1] == target[len(target) - 1]:
sub_cost = 0
else:
sub_cost = self._substitution_cost
# Minimum over deletion, insertion, and substitution/match.
return min(
self._naive_levenshtein(source[:-1], target)
+ self._deletion_cost,
self._naive_levenshtein(source, target[:-1])
+ self._insertion_cost,
self._naive_levenshtein(source[:-1], target[:-1]) + sub_cost,
)%%timeit
editdist = StringEdit1(1, 1)
editdist('æbstɹækt', 'æbstɹækt'), editdist('æbstɹækt', 'æbstɹækʃən'), editdist('æbstɹækʃən', 'æbstɹækt'), editdist('æbstɹækt', '')%%timeit
editdist = StringEdit2(1, 1)
editdist('æbstɹækt', 'æbstɹækt'), editdist('æbstɹækt', 'æbstɹækʃən'), editdist('æbstɹækʃən', 'æbstɹækt'), editdist('æbstɹækt', '')editdist = StringEdit2(1, 1)
editdist('æbstɹækʃən', 'æbstɹækt')
editdist.call_counterThat helps a lot. Why? Because it computes the distance for each pair of substrings only once. This is effectively what the Wagner–Fischer algorithm (Wagner and Fischer 1974) that you read about is doing. This is our first instance of a dynamic programming algorithm. The basic idea for Wagner–Fischer (and other algorithms we’ll use later in the class) is to cache the values for a function within a chart whose rows correspond to positions in the source string and whose columns correspond to positions in the target string.
import numpy as np
class StringEdit3(StringEdit2):
"""Distance between strings using the Wagner-Fischer algorithm.
Parameters
----------
insertion_cost : float
Cost of inserting a character (default 1.0).
deletion_cost : float
Cost of deleting a character (default 1.0).
substitution_cost : float | None
Cost of substituting a character.
"""
def __call__(self, source: str, target: str) -> float:
"""Compute the edit distance using a dynamic programming chart.
Parameters
----------
source : str
The source string.
target : str
The target string.
Returns
-------
float
The minimum edit distance.
"""
return self._wagner_fisher(source, target)
def _wagner_fisher(self, source: str, target: str) -> float:
n, m = len(source), len(target)
source, target = '#'+source, '#'+target
distance = np.zeros([n+1, m+1], dtype=float)
for i in range(1,n+1):
distance[i,0] = distance[i-1,0]+self._deletion_cost
for j in range(1,m+1):
distance[0,j] = distance[0,j-1]+self._insertion_cost
for i in range(1,n+1):
for j in range(1,m+1):
if source[i] == target[j]:
substitution_cost = 0.
else:
substitution_cost = self._substitution_cost
costs = np.array([distance[i-1,j]+self._deletion_cost,
distance[i-1,j-1]+substitution_cost,
distance[i,j-1]+self._insertion_cost])
distance[i,j] = costs.min()
return float(distance[n, m])editdist = StringEdit3(1, 1)
editdist('æbstɹækt', 'æbstɹækʃən')So why use Wagner-Fisher when we can just use memoization on the naive algorithm? The reason is that the chart used in Wagner-Fisher allows us to very easily store information about the implicit alignment between string elements. This notion of alignment is the same as the one we saw above when talking about how best to match up a square to the face of a cube when discussing boolean vectors.
So what do we need to do add to our previous implementation of Wagner-Fisher to store backtraces? Importantly, note that you will need to return a list of backtraces because there could be multiple equally good ones. (This point will come up for all of the dynamic programming algorithms we look at and, as we’ll see, is actually abstractly related to syntactic ambiguity.)
class StringEdit4(StringEdit3):
"""Distance, alignment, and edit paths between strings.
Parameters
----------
insertion_cost : float
Cost of inserting a character.
deletion_cost : float
Cost of deleting a character.
substitution_cost : float | None
Cost of substituting a character.
"""
def __call__(self, source: str,
target: str) -> tuple[
float,
list[list[tuple[int | None, int | None]]],
]:
"""Compute edit distance and all optimal alignments.
Parameters
----------
source : str
The source string.
target : str
The target string.
Returns
-------
tuple[float, list[list[tuple[int | None, int | None]]]]
The minimum distance and all optimal alignments. `None`
marks a gap on the corresponding side.
"""
return self._wagner_fisher(source, target)
def _wagner_fisher(
self,
source: str,
target: str,
) -> tuple[
float,
list[list[tuple[int | None, int | None]]],
]:
"""Compute minimum edit distance and alignment."""
n, m = len(source), len(target)
source, target = self._add_sentinel(source, target)
distance = np.zeros([n+1, m+1], dtype=float)
pointers = np.zeros([n+1, m+1], dtype=list)
pointers[0,0] = []
for i in range(1,n+1):
distance[i,0] = distance[i-1,0]+self._deletion_cost
pointers[i,0] = [(i-1,0)]
for j in range(1,m+1):
distance[0,j] = distance[0,j-1]+self._insertion_cost
pointers[0,j] = [(0,j-1)]
for i in range(1,n+1):
for j in range(1,m+1):
if source[i] == target[j]:
substitution_cost = 0.
else:
substitution_cost = self._substitution_cost
costs = np.array([distance[i-1,j]+self._deletion_cost,
distance[i-1,j-1]+substitution_cost,
distance[i,j-1]+self._insertion_cost])
distance[i,j] = costs.min()
best_edits = np.where(costs==distance[i,j])[0]
indices = [(i-1,j), (i-1,j-1), (i,j-1)]
pointers[i,j] = [indices[edit] for edit in best_edits]
pointer_backtrace = self._construct_backtrace(pointers,
idx=(n,m))
return distance[n,m], pointer_backtrace
def _construct_backtrace(
self,
pointers: np.ndarray,
idx: tuple[int, int],
) -> list[list[tuple[int | None, int | None]]]:
if idx == (0,0):
return [[]]
else:
pointer_backtrace = [backtrace+[self._alignment_step(prev_idx, idx)]
for prev_idx in pointers[idx]
for backtrace in self._construct_backtrace(pointers,
prev_idx)]
return pointer_backtrace
@staticmethod
def _alignment_step(
previous: tuple[int, int],
current: tuple[int, int],
) -> tuple[int | None, int | None]:
"""Translate one chart edge into source/target coordinates."""
old_i, old_j = previous
i, j = current
if (i-old_i, j-old_j) == (1, 0):
return i-1, None
if (i-old_i, j-old_j) == (0, 1):
return None, j-1
if (i-old_i, j-old_j) == (1, 1):
return i-1, j-1
raise ValueError('backtrace edge must be horizontal, vertical, or diagonal')
def _add_sentinel(
self,
source: str,
target: str,
) -> tuple[str, str]:
return '#' + source, '#' + targeteditdist = StringEdit4(1, 1)
editdist('æbstɹækʃən', 'æbstɹækt')The boundary cases expose gap-coordinate errors that nonempty examples can hide.
assert StringEdit4(1, 1)('', 'ab') == (
2, [[(None, 0), (None, 1)]],
)
assert StringEdit4(1, 1)('ab', '') == (
2, [[(0, None), (1, None)]],
)
assert StringEdit4(1, 1, 1)('a', 'b') == (
1, [[(0, 0)]],
)This isn’t particularly interpretable, so we can postprocess the output slightly to better see what’s going on.
class StringEdit5(StringEdit4):
"""Distance, alignment, and human-readable edit paths between strings.
Parameters
----------
insertion_cost : float
Cost of inserting a character.
deletion_cost : float
Cost of deleting a character.
substitution_cost : float | None
Cost of substituting a character.
"""
def __call__(self, source: str,
target: str) -> tuple[float, list[list[tuple[str, str]]]]:
distance, alignment = self._wagner_fisher(source, target)
return distance, [[(
'ε' if source_idx is None else source[source_idx],
'ε' if target_idx is None else target[target_idx],
) for source_idx, target_idx in path]
for path in alignment]editdist = StringEdit5(1, 1)
editdist('æbstɹækʃən', 'æbstɹækt')We now have a way of computing the distance between any two strings and recovering the alignment that achieves that distance. In the next section, we’ll put this to work: given a novel string and an entire lexicon, we’ll compute the distance from the novel string to every word in the lexicon and use the resulting distribution of distances to predict how wordlike the novel string sounds—implementing the “language as a region around known strings” idea from the module overview.
References
Levenshtein, Vladimir I. 1966. “Binary Codes Capable of Correcting Deletions, Insertions, and Reversals.” Soviet Physics Doklady 10 (8): 707–10.
Wagner, Robert A., and Michael J. Fischer. 1974. “The String-to-String Correction Problem.” Journal of the ACM 21 (1): 168–73. https://doi.org/10.1145/321796.321811.