Wu-Manber algorithm in Python
Bit-parallel approximate matching with k errors, using one bit vector per error level.
- File
wu_manber_approximate_matching.py, 51 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/wu-manber-approximate-matching/. The authors are listed in prov/QT4WGM.jsonl.
"""The Wu-Manber bit-parallel algorithm for approximate string matching.
Edit model: unit-cost insertion, deletion and substitution (Levenshtein),
no transpositions. Indices are 0-based. The result is the sorted list of end
positions j in 0..len(text) such that some substring text[h:j] has edit
distance at most k to the pattern.
Bit i - 1 of the vector R[e] (for 1 <= i <= m) is set after column j when the
pattern prefix of length i matches a substring that ends at j with at most e
errors. This is the shift-and form, with 1 for an active state.
"""
def wu_manber(pattern: str, text: str, k: int) -> list[int]:
"""Return the end positions of the approximate matches of pattern in text.
Precondition: k >= 0. Cost: O((k + 1) n) operations on m-bit vectors,
O(sigma + k) vectors of extra space.
"""
return [j for j, _ in wu_manber_states(pattern, text, k) if j >= 0]
def wu_manber_states(pattern: str, text: str, k: int) -> list[tuple[int, list[int]]]:
"""Return, for each column j = 0..n, the pair (j or -1, [R[0], ..., R[k]]),
where j marks a reported end and -1 a column without a match."""
if k < 0:
raise ValueError("k must be at least 0")
m = len(pattern)
full = (1 << m) - 1
high = 1 << (m - 1) if m > 0 else 0
masks: dict[str, int] = {} # @line L1
for i, ch in enumerate(pattern):
masks[ch] = masks.get(ch, 0) | (1 << i)
r = [(1 << min(e, m)) - 1 for e in range(k + 1)] # @line L2
hit = m == 0 or (r[k] & high) != 0 # @line L3
out = [(0 if hit else -1, r[:])]
for j in range(1, len(text) + 1): # @line L4
b = masks.get(text[j - 1], 0)
old = r[0]
r[0] = ((old << 1) | 1) & b # @line L5
for e in range(1, k + 1): # @line L6
prev_old = old
old = r[e]
r[e] = ((((old << 1) | 1) & b) # @line L7
| (prev_old << 1) | 1 | prev_old | (r[e - 1] << 1)) & full
hit = m == 0 or (r[k] & high) != 0 # @line L8
out.append((j if hit else -1, r[:]))
return out
test_wu_manber_approximate_matching.py 70 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/wu-manber-approximate-matching/. The authors are listed in prov/QT4WGM.jsonl.
import random
import pytest
from wu_manber_approximate_matching import wu_manber, wu_manber_states
def edit_distance(a: str, b: str) -> int:
d = [[i + j if i * j == 0 else 0 for j in range(len(b) + 1)] for i in range(len(a) + 1)]
for i in range(1, len(a) + 1):
for j in range(1, len(b) + 1):
d[i][j] = min(d[i - 1][j] + 1, d[i][j - 1] + 1, d[i - 1][j - 1] + (a[i - 1] != b[j - 1]))
return d[len(a)][len(b)]
def sellers_table(p: str, t: str) -> list[list[int]]:
c = [[i if j == 0 else 0 for j in range(len(t) + 1)] for i in range(len(p) + 1)]
for i in range(1, len(p) + 1):
for j in range(1, len(t) + 1):
c[i][j] = min(c[i - 1][j] + 1, c[i][j - 1] + 1, c[i - 1][j - 1] + (p[i - 1] != t[j - 1]))
return c
def brute(pattern: str, text: str, k: int) -> list[int]:
return [
j for j in range(len(text) + 1)
if min(edit_distance(pattern, text[h:j]) for h in range(j + 1)) <= k
]
def test_example() -> None:
assert wu_manber("GATAA", "CAGATAAGAGAA", 1) == [6, 7, 8, 12]
assert wu_manber("GATAA", "CAGATAAGAGAA", 0) == [7]
assert wu_manber("survey", "surgery", 2) == [5, 6, 7]
def test_edges() -> None:
assert wu_manber("", "", 0) == [0]
assert wu_manber("", "ab", 0) == [0, 1, 2]
assert wu_manber("ab", "", 1) == []
assert wu_manber("ab", "", 2) == [0]
assert wu_manber("ab", "xyz", 5) == [0, 1, 2, 3]
with pytest.raises(ValueError):
wu_manber("a", "a", -1)
def test_vectors_are_the_active_sets() -> None:
rng = random.Random(1992)
for alphabet in ("ab", "acgt"):
for _ in range(200):
p = "".join(rng.choice(alphabet) for _ in range(rng.randint(0, 7)))
t = "".join(rng.choice(alphabet) for _ in range(rng.randint(0, 10)))
k = rng.randint(0, 4)
c = sellers_table(p, t)
for j, (_, r) in enumerate(wu_manber_states(p, t, k)):
for e in range(k + 1):
want = sum(1 << (i - 1) for i in range(1, len(p) + 1) if c[i][j] <= e)
assert r[e] == want
def test_random_against_brute_force() -> None:
rng = random.Random(35)
for alphabet in ("ab", "acgt"):
for _ in range(400):
p = "".join(rng.choice(alphabet) for _ in range(rng.randint(0, 6)))
t = "".join(rng.choice(alphabet) for _ in range(rng.randint(0, 12)))
k = rng.randint(0, 4)
assert wu_manber(p, t, k) == brute(p, t, k)
What the next level needs
- the tests did not run in this build