Edit distance where a substitution only costs the first time

Question

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.

Answer 1

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:
      mapping.add((str1[p1], str2[p2]))
      cost3 = g(p1 + 1, p2 + 1) + 1
      mapping.discard((str1[p1], str2[p2]))
    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...

Answer 2

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!

Answer 3

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:
            closed_ij.add((g, 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)]))

Answer 4

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)

Answer 5

Rust

Port of @tsh's answer

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);
    }
}