Myers bit-parallel algorithm in Python
Encodes DP columns as bit vectors, giving O(nm/w) time for edit distance and approximate matching.
- File
myers_bit_parallel_algorithm.py, 81 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/myers-bit-parallel-algorithm/. The authors are listed in prov/MXPHKN.jsonl.
"""Myers' bit-parallel algorithm for approximate string matching, in the form
given by Hyyrö.
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 pv (of mv) is set when C[i][j] - C[i-1][j] is +1 (is -1), where C
is the Sellers table of the current column j.
"""
def myers(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(n) operations on m-bit vectors.
"""
return [j for j, score, _, _ in myers_columns(pattern, text) if score <= k]
def myers_columns(pattern: str, text: str) -> list[tuple[int, int, int, int]]:
"""Return (j, C[m][j], pv, mv) for each column j = 0..n."""
m = len(pattern)
if m == 0:
return [(j, 0, 0, 0) for j in range(len(text) + 1)]
full = (1 << m) - 1
high = 1 << (m - 1)
peq: dict[str, int] = {} # @line L1
for i, ch in enumerate(pattern):
peq[ch] = peq.get(ch, 0) | (1 << i)
pv, mv, score = full, 0, m # @line L2
out = [(0, score, pv, mv)] # @line L3
for j in range(1, len(text) + 1): # @line L4
eq = peq.get(text[j - 1], 0) # @line L5
xv = eq | mv
xh = ((((eq & pv) + pv) & full) ^ pv) | eq # @line L6
ph = mv | (~(xh | pv) & full) # @line L7
mh = pv & xh
if ph & high: # @line L8
score += 1
elif mh & high:
score -= 1
ph = (ph << 1) & full # @line L9
mh = (mh << 1) & full
pv = mh | (~(xv | ph) & full) # @line L10
mv = ph & xv
out.append((j, score, pv, mv)) # @line L11
return out
def edit_distance(a: str, b: str) -> int:
"""Return the Levenshtein distance of a and b with the same vectors, where
the first row of the table is C[0][j] = j instead of 0."""
m = len(a)
if m == 0:
return len(b)
full = (1 << m) - 1
high = 1 << (m - 1)
peq: dict[str, int] = {}
for i, ch in enumerate(a):
peq[ch] = peq.get(ch, 0) | (1 << i)
pv, mv, score = full, 0, m
for ch in b:
eq = peq.get(ch, 0)
xv = eq | mv
xh = ((((eq & pv) + pv) & full) ^ pv) | eq
ph = mv | (~(xh | pv) & full)
mh = pv & xh
if ph & high:
score += 1
elif mh & high:
score -= 1
ph = ((ph << 1) | 1) & full # row 0 grows by 1 in each column
mh = (mh << 1) & full
pv = mh | (~(xv | ph) & full)
mv = ph & xv
return score
test_myers_bit_parallel_algorithm.py 74 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/myers-bit-parallel-algorithm/. The authors are listed in prov/MXPHKN.jsonl.
import random
import pytest
from myers_bit_parallel_algorithm import edit_distance, myers, myers_columns
def lev(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(lev(pattern, text[h:j]) for h in range(j + 1)) <= k
]
def test_example() -> None:
assert myers("GATAA", "CAGATAAGAGAA", 1) == [6, 7, 8, 12]
assert myers("GATAA", "CAGATAAGAGAA", 0) == [7]
assert myers("survey", "surgery", 2) == [5, 6, 7]
assert [s for _, s, _, _ in myers_columns("survey", "surgery")] == [6, 5, 4, 3, 3, 2, 2, 2]
assert edit_distance("kitten", "sitting") == 3
def test_edges() -> None:
assert myers("", "", 0) == [0]
assert myers("", "ab", 0) == [0, 1, 2]
assert myers("ab", "", 1) == []
assert myers("ab", "", 2) == [0]
assert myers("ab", "xyz", 5) == [0, 1, 2, 3]
assert edit_distance("", "abc") == 3
assert edit_distance("abc", "") == 3
def test_vectors_are_the_column_differences() -> None:
rng = random.Random(1999)
for alphabet in ("ab", "acgt"):
for _ in range(200):
p = "".join(rng.choice(alphabet) for _ in range(rng.randint(1, 8)))
t = "".join(rng.choice(alphabet) for _ in range(rng.randint(0, 10)))
c = sellers_table(p, t)
for j, score, pv, mv in myers_columns(p, t):
assert score == c[len(p)][j]
for i in range(1, len(p) + 1):
v = c[i][j] - c[i - 1][j]
assert (pv >> (i - 1)) & 1 == (v == 1)
assert (mv >> (i - 1)) & 1 == (v == -1)
def test_random_against_brute_force() -> None:
rng = random.Random(2001)
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 myers(p, t, k) == brute(p, t, k)
assert edit_distance(p, t) == lev(p, t)
What the next level needs
- the tests did not run in this build