Edit distance where a substitution only costs the first time

This challenge is about the following variant of edit distance. Say we have a cost of 1 for inserts, deletes and substitutions as usual with one exception. A substitution for a given letter x for a letter y only costs 1 the first time. Any further substitutions of x for y cost 0.

As simple examples:

A = apppple
B = attttle


cost 1 (not 4) to transform A into B. This is because we change p for a t four times but we only charge for the first one.

A = apxpxple
B = atxyxtle


cost 2 (not 3) as the p to t substitution only costs 1 even though we do it twice.

A = apppple
B = altttte


cost 3. p -> l costs 1. p -> t costs 1 as does l -> t.

If we assume the total length of the input is n, you must give the big O() of your algorithm. The best time complexity wins.

We should assume the alphabet size can be as large as n.

• @SunnyMoon It isn't because it's a different distance function. The important twist is "A substitution for a given letter x for a letter y only costs 1 the first time. " – Anush Nov 12 '20 at 20:56
• Are we only dealing with substitutions here? Or are insertions & deletions allowed? – Shaggy Nov 12 '20 at 21:48
• If we can assume that there are only 26 possible letters (or some constant), it seems one can "cheat" the time complexity by iterating over all possible sets of letter-to-letter replacements that are used. There are an astronomical number of these, but it's a constant independent of the length of the input so doesn't affect the algorithm's big-O. – xnor Nov 12 '20 at 22:02
• @Shaggy insertion and deletions are allowed. The cost is 1 per insert/delete. – Anush Nov 12 '20 at 22:38
• @SunnyMoon Note that being a dupe on Code Golf requires that a serious solution to one challenge is expected to be a serious contender (with minor changes) to the other. Even if the task is 100% identical, a code-golf and a fastest-algorithm are not likely to be dupes, unless the shortest algorithm is also the fastest. – Bubbler Nov 13 '20 at 5:40

Naive recursive algorithm $$\ O(\frac{(1+\sqrt2)^n}{\sqrt{n}}) \$$

The basic idea is try adding, removing and replacing one character and recursive.

def f(str1, str2):
l1, l2 = len(str1), len(str2)
mapping = set()
def g(p1, p2):
if l1 == p1 or l2 == p2:
return l1 + l2 - p1 - p2
cost1 = g(p1 + 1, p2) + 1
cost2 = g(p1, p2 + 1) + 1
if str1[p1] != str2[p2] and (str1[p1], str2[p2]) not in mapping:
cost3 = g(p1 + 1, p2 + 1) + 1
else:
cost3 = g(p1 + 1, p2 + 1)
return min(cost1, cost2, cost3)
return g(0, 0)


The function g will be invoked $$\ O(\frac{(1+\sqrt2)^n}{\sqrt{n}}) \$$ times[1]. While adding, discarding, or checking an element in the mapping cost $$\ O(1) \$$ on average.

[1]: Based on https://oeis.org/A027618 ; I don't know why actually...

• (str1[0], str2[0]) not in mapping takes $O(1)$ (a set is a hash table), but mapping | {(str1[0], str2[0])}) takes $O(n)$ (it constructs a new set), as do str1[1:] and str2[1:] (they construct a new str). But A001850(n/2) is $O{\left(\frac{(1 + \sqrt 2)^n}{\sqrt n}\right)}$, so you can take off a $\sqrt n$ factor from the number of invocations. – Anders Kaseorg Nov 13 '20 at 4:38
• @AndersKaseorg str1[1:] can be replaced by passing another variable pos1; mapping | {...} can be replaced by mapping |= {...} The O(n) is not necessary. – tsh Nov 13 '20 at 5:31
• @AndersKaseorg Sorry but I didn't find the formula $\frac{(1+\sqrt2)^n}{\sqrt{n}}$ in the page you linked. And also cannot get it by myself. Could you tell me more about where I can find this formula? – tsh Nov 13 '20 at 5:57
• It follows from Hirschhorn’s comment (and proof), with Stirling’s approximation for $\binom{2n}{n}$. – Anders Kaseorg Nov 13 '20 at 6:17
• Since every non-leaf node in the call tree has degree three, just over 2/3 of the nodes are leaves, so it suffices just to count those. (Or see A027618.) Your argument for $\Omega((1 + \sqrt 2)^n)$ pessimistically assumes that the recursion terminates at p1 + p2 == l1 + l2, while it actually terminates at p1 == l1 or p2 == l2. – Anders Kaseorg Nov 13 '20 at 6:58

Based on the Wagner–Fischer algorithm, each cell (here I just store a single row rather than the whole matrix) contains a list of tuples of edit distances and their associated transpositions. The matrix itself requires $$\ O(n^2) \$$ to traverse, but I don't know the cost of tracking the sets of transpositions, which varies across the grid roughly in line with Delannoy numbers. The code also merges distances with the same transpositions; I don't know whether that increases or reduces the complexity.

def dist(fr, to):
costs = [[(i, frozenset())] for i in range(len(fr) + 1)]
for j in range(len(to)):
editcosts = costs[0]
delcosts = costs[0] = [(costs[0][0][0] + 1, frozenset())]
for i in range(len(fr)):
inscosts = costs[i + 1]
allcosts = [(c + 1, s) for c, s in delcosts + inscosts] + [(c, s) if fr[i] == to[j] or fr[i] + to[j] in s else (c + 1, s | {fr[i] + to[j]}) for c, s in editcosts]
allcosts = [min((c, s) for c, s in allcosts if s == e) for e in {s for c, s in allcosts}]
costs[i + 1] = allcosts
editcosts = inscosts
delcosts = allcosts
return min(costs.pop())

print(dist("apxpxple", "atxyxtle"))


Try it online!

• I'll try to run an experiment with larger inputs to see what the running time might be. Unless someone can analyse this answer? – Anush Nov 18 '20 at 19:30
• @Anush if no distances have the same transpositions (which for an alphabet of size n is not far from average case of 1), then for each cell (i,j) from the Wagner-Fischer algorithm, this forms a list with one tuple (c,s) for every path from (0,0) to (i,j). The number of paths are the Delannoy numbers D(i,j), so the total would be $O(\sum_{i=0}^n\sum_{j=0}^n D(m,n))$. ias.sabanciuniv.edu/documents/kiselmanpaper.pdf proves that $3^{min(i,j)}\leq D(i,j)\leq(\sqrt{2}+1)^{i+j}$. Summing those across (i,j), we can show that the algorithm is $O((\sqrt{2}+1)^{2n})$ and $\Omega(3^n)$ – fireflame241 Nov 19 '20 at 0:28
• The main two differences between this answer and tsh's "Naive recursive algorithm" is that (1) Wagner-Fischer proceeds forward from (0,0) rather than backwards from (n,n), and (2) this answer generates a list in O(D(i,j)) for each cell rather than having O(D(n,n)) total runtime – fireflame241 Nov 19 '20 at 0:36
• Maybe rewrite the two lines about allcosts may reduce complexity. But I'm not sure since I actually didn't understand what happened there. – tsh Nov 19 '20 at 2:53
• The time complexity may be $O(e^{1.2n})$ ~ $O(e^{1.4n})$. Guessed based on first few values – tsh Nov 19 '20 at 3:13

A* Search, unknown complexity. Seems to at least achieve $$\o(2^{n})\$$ and is quite fast in practice. No reason to think it's optimal, but making it better seemed like too much work.

import collections
import heapq

def dist(s0, s1):
l0, l1 = len(s0), len(s1)
suffix_sets = tuple([frozenset(s[i:]) for i in range(len(s) + 1)] for s in (s0, s1))
q = []
def push_node(g, i, j, substs):
suff0, suff1 = suffix_sets[0][i], suffix_sets[1][j]
substs = frozenset(
t for t in substs if t[0] in suff0 and t[1] in suff1
)
substs0, substs1 = {c0 for c0, _ in substs}, {c1 for _, c1 in substs}
h0 = abs((l0 - i) - (l1 - j))
h1 = max(len((suff0 - suff1) - substs0), len((suff1 - suff0) - substs1))
h = max(h0, h1)
heapq.heappush(q, (g + h, -g, -len(substs), i, j, substs))
closed = collections.defaultdict(lambda: set())
push_node(0, 0, 0, frozenset())
while True:
(f, neg_g, _, i, j, substs) = heapq.heappop(q)
g = -neg_g
if (i == l0) or (j == l1):
return f
closed_ij = closed[(i, j)]
if (g, substs) in closed_ij:
continue
subsumed = False
for (g1, substs1) in closed_ij:
if g1 <= g and substs.issubset(substs1):
subsumed = True
break
if subsumed:
continue
if substs:
c0, c1 = s0[i], s1[j]
if (c0 == c1) or ((c0, c1) in substs):
push_node(g, i + 1, j + 1, substs)
continue
push_node(g + 1, i + 1, j, substs)
push_node(g + 1, i, j + 1, substs)
push_node(g + 1, i + 1, j + 1, substs | frozenset([(c0, c1)]))

• I also have a version which uses a Trie for the subsumption testing, and is a bit faster. – user1502040 Dec 24 '20 at 6:03
• Please add another answer with the trie version too! – Anush Dec 24 '20 at 9:05
• And how might this be improved? – Anush Dec 24 '20 at 14:00
• @Anush The most obvious way would be to improve the heuristic. – user1502040 Dec 27 '20 at 0:35

I worked a bit more on this, against my better judgement :). A* search, more complicated than my other answer.

import collections
import heapq
import itertools

class SubsumptionTrie:
def __init__(self):
self.max_len = 0
self.min_v = float('inf')
self.children = {}
def insert(self, k, v, i=0):
self.max_len = max(self.max_len, len(k) - i)
self.min_v = min(self.min_v, v)
if i >= len(k):
return
c = k[i]
if c not in self.children:
self.children[c] = SubsumptionTrie()
self.children[c].insert(k, v, i + 1)
def subsumption_test(self, k, v, i=0):
if i >= len(k):
return self.min_v <= v
if (self.min_v > v) or (self.max_len < (len(k) - i)):
return False
c = k[i]
for b, child in self.children.items():
if b < c and self.children[b].subsumption_test(k, v, i):
return True
if c in self.children:
return self.children[c].subsumption_test(k, v, i + 1)
return False

def heuristic(s0, s1, n0, n1, substs):
substs0 = {a for a, _ in substs} | set(s1.keys())
substs1 = {b for _, b in substs} | set(s0.keys())
if n0 > n1:
s0, s1 = s1, s0
n0, n1 = n1, n0
substs0, substs1 = substs1, substs0
h0 = 0
m1 = sum(n for c, n in s1.items() if c in substs1)
for c, n in sorted(s0.items(), key=lambda t: t[1]):
if (c not in substs0) or (m1 < 0):
h0 += 1
else:
m1 -= n
h1 = 0
for c, n in sorted(s1.items(), key=lambda t: (t[0] in substs1, t[1]), reverse=True):
n0 -= n
if n0 < 0:
excess = min(-n0, n)
h1 += excess
n -= excess
if n and (c not in substs1):
h1 += 1
return max(h0, h1)

def dist(s0, s1):
labels = itertools.count()
relabel = collections.defaultdict(lambda: next(labels))
s0, s1 = [relabel[c] for c in s0], [relabel[c] for c in s1]
label_count = next(labels)
l0, l1 = len(s0), len(s1)
suffix_counts = tuple([collections.Counter(s[i:]) for i in range(len(s) + 1)] for s in (s0, s1))
q = []
def push_node(g, i, j, substs):
suff0, suff1 = suffix_counts[0][i], suffix_counts[1][j]
substs = frozenset(
t for t in substs if t[0] in suff0 and t[1] in suff1
)
h = heuristic(suff0, suff1, l0 - i, l1 - j, substs)
heapq.heappush(q, (g + h, -g, -len(substs), i, j, substs))
push_node(0, 0, 0, frozenset())
closed = collections.defaultdict(lambda: SubsumptionTrie())
while True:
(f, neg_g, _, i, j, substs) = heapq.heappop(q)
g = -neg_g
if (i == l0) or (j == l1):
return f
closed_ij = closed[(i, j)]
k = sorted(a * label_count + b for a, b in substs)
if closed_ij.subsumption_test(k, g):
continue
closed_ij.insert(k, g)
c0, c1 = s0[i], s1[j]
if (c0 == c1) or ((c0, c1) in substs):
push_node(g, i + 1, j + 1, substs)
continue
push_node(g + 1, i + 1, j + 1, substs | frozenset(((c0, c1),)))
push_node(g + 1, i + 1, j, substs)
push_node(g + 1, i, j + 1, substs)

• Can you say anything about its performance? – Anush Dec 27 '20 at 22:42
• @Anush The average case seems to be about O(1.2^n). – user1502040 Dec 30 '20 at 3:12