# 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. "
– user7467
Nov 12, 2020 at 20:56
• Are we only dealing with substitutions here? Or are insertions & deletions allowed? Nov 12, 2020 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, 2020 at 22:02
• @Shaggy insertion and deletions are allowed. The cost is 1 per insert/delete.
– user7467
Nov 12, 2020 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. Nov 13, 2020 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. While adding, discarding, or checking an element in the mapping cost $$\ O(1) \$$ on average.

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

• (str1, str2) not in mapping takes $O(1)$ (a set is a hash table), but mapping | {(str1, str2)}) 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. Nov 13, 2020 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, 2020 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, 2020 at 5:57
• It follows from Hirschhorn’s comment (and proof), with Stirling’s approximation for $\binom{2n}{n}$. Nov 13, 2020 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. Nov 13, 2020 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
delcosts = costs = [(costs + 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?
– user7467
Nov 18, 2020 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)$ Nov 19, 2020 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 Nov 19, 2020 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, 2020 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, 2020 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[i], suffix_sets[j]
substs = frozenset(
t for t in substs if t in suff0 and t 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. Dec 24, 2020 at 6:03
– user7467
Dec 24, 2020 at 9:05
• And how might this be improved?
– user7467
Dec 24, 2020 at 14:00
• @Anush The most obvious way would be to improve the heuristic. Dec 27, 2020 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):
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 in substs1, t), 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[i], suffix_counts[j]
substs = frozenset(
t for t in substs if t in suff0 and t 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?
– user7467
Dec 27, 2020 at 22:42
• @Anush The average case seems to be about O(1.2^n). Dec 30, 2020 at 3:12

# Rust

run it on Rust Playground!

use std::cmp;
use std::collections::HashSet;
use std::time::Instant;

fn f(str1: &str, str2: &str) -> usize {
let mut mapping = HashSet::new();

fn g(p1: usize, p2: usize, str1: &str, str2: &str, mapping: &mut HashSet<(char, char)>) -> usize {
let l1 = str1.len();
let l2 = str2.len();

if l1 == p1 || l2 == p2 {
return l1 + l2 - p1 - p2;
}

let cost1 = g(p1 + 1, p2, str1, str2, mapping) + 1;
let cost2 = g(p1, p2 + 1, str1, str2, mapping) + 1;

let c1 = str1.chars().nth(p1).unwrap();
let c2 = str2.chars().nth(p2).unwrap();

let cost3 = if c1 != c2 && !mapping.contains(&(c1, c2)) {
mapping.insert((c1, c2));
let cost = g(p1 + 1, p2 + 1, str1, str2, mapping) + 1;
mapping.remove(&(c1, c2));
cost
} else {
g(p1 + 1, p2 + 1, str1, str2, mapping)
};

cmp::min(cost1, cmp::min(cost2, cost3))
}

g(0, 0, str1, str2, &mut mapping)
}

fn main() {
let test_cases = [
("apppple", "attttle", 1),
("apxpxple", "atxyxtle", 2),
("apppple", "altttte", 3),
];

for (a, b, expected) in test_cases.iter() {
let start = Instant::now();
let result = f(a, b);
let duration = start.elapsed();
println!("A: {}, B: {}, Output: {}, Expected: {}", a, b, result, expected);
println!("Time elapsed: {:?}", duration);
}
}