Multiple longest common subsequence (another algorithm)

Question

Here is the algorithm:

Input: a set of \$M\$ strings \$S=\{s_i\}\$ of length \$N\$ with alphabet \$\Sigma\$

Make an index of each string by chars. Maintain a position \$p=(p_i)\$ where \$p_i\$ is an index into string \$s_i\$. Set the initial position such that \$\forall i: p_i = |s_i|+1\$ (1-based indexing).

While chars left in strings do:

• Select the next point among matches (same char in every string) by looping through the alphabet shared by the strings and picking the point with the minimal m-dimensional euclidean distance to the current point.

Is this \$O(|\Sigma|*M*\log(N))\$?

Here is (quite bad) sample python code:

from collections import defaultdict
import sys, bisect

def mlcs(strings):
    strings = list(set(strings))
    indexes = dict()

    sigma = set()
    for x in strings:
        sigma = sigma.union(set(x))

    for x in strings:
        tx = defaultdict(list)
        for i in range(len(x)):
            tx[x[i]] += [i]
        indexes[x] = tx

    def match(pos):
        s = True
        for i in range(len(strings)-1):
            x = strings[i]
            y = strings[i+1]
            s &= x[pos[x]] == y[pos[y]]
        return s

    def distance(ind1, ind2):
        d = 0
        for x in strings:
            d += (ind1[x] - ind2[x])**2
        return d

    def find(xs, y):
        i = bisect.bisect_right(xs, y)
        return i-1

    pos = {}
    for x in strings:
        pos[x] = len(x)
    result = ""
    loop = True
    while loop:
        for x in strings:
            if pos[x] < 0:
                loop = False
        if not loop:
            break
        indr = None
        dr = 12777216
        for c in sigma:
            ind = defaultdict(int)
            for x in strings:
                indxc = indexes[x][c]
                if indxc:
                    i = find(indxc, pos[x] - 1)
                    ind[x] = indxc[i]
                    if ind[x] > pos[x] - 1:
                        ind[x] = -128                    
                else:
                    ind[x] = -128
            d = distance(ind, pos)
            if d < dr:
                indr = ind
                dr = d
        x = strings[0]
        def notnegative(p):
            r = True
            for x in p.values():
                if x < 0:
                    r = False
            return r
        if notnegative(indr):
            result = x[indr[x]] + result
        pos = indr
    return result

def check(string, seq):
    i = 0
    j = 0
    while i < len(string) and j < len(seq):
        if string[i] == seq[j]:
            j += 1
        i += 1
    return len(seq) - j

def checkall(strings, seq):
    for x in strings:
        a = check(x, seq)
        if not a == 0:
            print(x, seq, a)
            return False
    return True

strings = []
sigma = ["a", "b", "c", "d"]
#sigma = ["a", "b"]
sigmax = ["e", "f", "g"]
import random
Nx = 67
N = 128
M = 13 # Number of strings
x = ""
for i in range(Nx):
    x += random.choice(sigmax)
for i in range(M):    
    string = ""
    for j in range(N):
        string += random.choice(sigma)
    indexes = list(range(N))
    random.shuffle(indexes)
    indexes = sorted(indexes[:len(x)])
    for j in range(len(x)):
        string = string[:indexes[j]]+x[j]+string[indexes[j]+1:]
    strings += [string]

#strings = ["abbab", "ababa", "abbba"]
#strings = ["abab", "baba", "abaa"]
#strings = ["bacda", "abcde", "decac"]
#strings = ["babbabbb", "bbbaabaa", "abbbabab", "abbababa"]
print("Strings:", strings)
l = mlcs(strings)
print("Lcs:", l, len(l), checkall(strings, l))

Answer 1

1. Analysis

This is a greedy algorithm: at each step it prefixes one letter to the result, choosing the letter than minimizes the change in the indexes into the strings.

The main loop of the algorithm adds one letter to the result, so can execute at most \$ N \$ times. Each iteration of the loop considers each letter in the alphabet \$ Σ \$, and each of the \$ M \$ strings, and searches (by bisection) the list of occurrences of that letter in that string, taking \$ O(\log N) \$. Thus the overall runtime is \$ O(\left|Σ\right|NM\log N) \$. (You missed a factor of \$ N \$.)

This is polynomial time, but the longest common subsequence problem is NP-hard, so the algorithm does not solve the problem in the general case.

Here's a simple case where mlcs returns the wrong result:

>>> mlcs(('ab', 'aba'))
'a'

2. Review

1. There's no docstring for mlcs. What does this function do? How do I call it? What does it return?

2. The case of no strings is not handled:

   >>> mlcs(())
   Traceback (most recent call last):
       ...
   IndexError: list index out of range
   

   This should probably be a ValueError, like max(()).

3. The case where all strings are empty is not handled:

   >>> mlcs(('', ''))
   Traceback (most recent call last):
       ....
   AttributeError: 'NoneType' object has no attribute 'values'
   

   This should return the empty string.

4. The code takes the union of the letters in the strings. But the only letters that can appear in the longest common subsequence are those letters that appear in all the strings, so the intersection is needed here here, not the union.

5. The code accumulates a list by repeated addition:

   for i in range(len(x)):
       tx[x[i]] += [i]
   

   This wastes space because Python might have to allocate a new list each time and copy the old list across. It's more efficient to accumulate a list using the append method.

6. Iteration over the indexes of a sequence can often be simplified using enumerate, for example the loop above can be written:

   for i, letter in enumerate(x):
       tx[letter].append(i)
   

7. pos is a dictionary mapping input string to the current position in that string in the search. Since strings is a list of strings, it would make more sense for pos to be a list of positions, so that pos[i] was the current position in strings[i]. Then match simplifies to:

   def match(pos):
       c = strings[0][pos[0]]
       return all(c == s[p] for s, p in zip(strings, pos))
   

8. The function match is not called, so could be omitted.

9. The initialization of pos:

   pos = {}
   for x in strings:
       pos[x] = len(x)
   

   Could be written as a dictionary comprehension:

   pos = {s: len(s) for s in strings}
   

   or, if you follow my advice above about making pos into a list, a list comprehension:

   pos = [len(s) for s in strings]
   

10. This loop:

    for x in strings:
        if pos[x] < 0:
            loop = False
    if not loop:
        break
    

    always tests all elements of pos, so that even if the first one is negative, it will go on and uselessly test all the others. It would be better to fail as soon as possible:

    if any(pos[s] < 0 for s in strings):
        break
    

    or, if you follow my advice above about making pos into a list:

    if any(p < 0 for p in pos):
        break
    

11. Note that the change discussed above gets rid of the only assignment to loop, so that we have:

    loop = True
    while loop:
        if any(p < 0 for p in pos):
            break
    

    which simplifies to:

    while all(p >= 0 for p in pos):
    

12. Instead of setting ind[x] and then setting it back to the default value based on a condition:

    ind = defaultdict(int)
    for x in strings:
        indxc = indexes[x][c]
        if indxc:
            i = find(indxc, pos[x] - 1)
            ind[x] = indxc[i]
            if ind[x] > pos[x] - 1:
                ind[x] = -128                    
        else:
            ind[x] = -128
    

    initialize it to the default value and set it only if the condition does not apply:

    ind = defaultdict(lambda:-128)
    for x in strings:
        indxc = indexes[x][c]
        if indxc:
            i = find(indxc, pos[x] - 1)
            if indxc[i] <= pos[x] - 1:
                ind[x] = indxc[i]
    

13. Just as I recommended making pos into a list, I also recommend making ind into a list. If you did, then distance would become:

    def distance(v, w):
        return sum((i - j)**2 for i, j in zip(v, w)):
    

    and since this is only called from one place, you could easily inline it.

14. Similarly, indexes would be more conveniently organized so that indexes[letter][i] is a list of the indexes of the occurrences of letter in strings[i].

15. A magic number like this:

    dr = 12777216
    

    needs to be explained. It looks to me as though this needs to be some number larger than the biggest possible distance between pos and ind. But for very long strings (tens of thousands of letters), this won't be the case. It would be more reliable to start at infinity:

    min_distance = float('inf')
    

    But it would be better to reorganize the code to call the built-in min function, and catch the ValueError that is raised when there are no elements. See the revised code in §3 below for how this is done.

16. Similarly, the magic number -128 needs explanation. The idea seems to be that this is so far away from the index of any character in the string that the distance will be large and the candidate will be rejected. But is this really true? If the strings are longer than 128 characters then a candidate might be included erroneously. It would be better to reject candidates where bisect_right returns 0.

17. The handling of the case where bisect_right returns 0 seems incorrect to me. In this case find returns -1, and then this leads to the assignment ind[x] = indxc[-1] which succeeds (getting the last element in indxc) but this is wrong.

3. Revised code

from bisect import bisect

def mlcs(strings):
    """Return a long common subsequence of the strings.
    Uses a greedy algorithm, so the result is not necessarily the
    longest common subsequence.

    """
    if not strings:
        raise ValueError("mlcs() argument is an empty sequence")
    strings = list(set(strings)) # deduplicate
    alphabet = set.intersection(*(set(s) for s in strings))

    # indexes[letter][i] is list of indexes of letter in strings[i].
    indexes = {letter:[[] for _ in strings] for letter in alphabet}
    for i, s in enumerate(strings):
        for j, letter in enumerate(s):
            if letter in alphabet:
                indexes[letter][i].append(j)

    # pos[i] is current position of search in strings[i].
    pos = [len(s) for s in strings]

    # Generate candidate positions for next step in search.
    def candidates():
        for letter, letter_indexes in indexes.items():
            distance, candidate = 0, []
            for ind, p in zip(letter_indexes, pos):
                i = bisect.bisect_right(ind, p - 1) - 1
                q = ind[i]
                if i < 0 or q > p - 1:
                    break
                candidate.append(q)
                distance += (p - q)**2
            else:
                yield distance, letter, candidate

    result = []
    while True:
        try:
            # Choose the closest candidate position, if any.
            _, letter, pos = min(candidates())
        except ValueError:
            return ''.join(reversed(result))
        result.append(letter)

Answer 2

The question has been answered above, but here are some further thoughts:

The problem with solving LCS exactly comes from the fact that there are an exponential number of matches in the worst case, $O(N^M)$. This is the case when no deduplication is used and all the strings are equal. The LCS grid is a DAG, where each grid cell is a vertex, so a longest path and longest sequence could be found in time linear in the number of matches. But for this to be efficient (polynomial time) only a polynomial number of matches must be considered, which may lead to an inexact solution.

The number of matches is exactly $\Sigma_{k}\Pi_{i}|s_i|_{a_k}$, where $a_k$ is the kth letter of the (common) alphabet and $s_i$ is the ith string, and $|s_i|_{a_k}$ is the number of occurrences of a letter in a string. This is exponential in the number of strings.