# https://opensource.org/licenses/MIT

'''
The MIT License (MIT)
Copyright (c) 2016 Philip (Pak Lam) Lee

Permission is hereby granted, free of charge, 
to any person obtaining a copy of this software 
and associated documentation files (the "Software"), 
to deal in the Software without restriction, 
including without limitation the rights to use, 
copy, modify, merge, publish, distribute, sublicense, 
and/or sell copies of the Software, 
and to permit persons to whom the Software is furnished to do so, 
subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, 
INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, 
DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
'''

import array
import itertools as it

def getScore(sequences):
    
    N = len(sequences)
    L = map(len, sequences)
    C = L[:]
    N = len(L)

    for i in range(N): C[i] = L[i] + 1
    for i in range(N-2,0,-1): C[i] = C[i] * C[i+1]
    # the smallest dimension is updated the most...
    
    # previous nearest neighbors, the previous "diagonal" (when all characters match)
    nn = [sum(combo) for i in range(1,N) for combo in it.combinations(C[1:]+[1],i)]
    diag = sum(C[1:]+[1])

    def getIndex(index):
        assert len(index) == N      
        ix = 0
        for i in range(N-1):
            ix += C[i+1] * index[i]
        ix += index[N-1]
        return ix
    
    # indexes the characters from each sequence, it is shifted by 1 [1..L[s]]
    def aligned(index):
        cc = sequences[0][index[0]-1] # current character
        for i in range(1,N):
            if sequences[i][index[i]-1] != cc: 
                return False
        return True

    start = getIndex(N*[1])
    end = C[0]*C[1]
    s = array.array('i', end*[0])
    
    I = N * [1]
    i = start
    while i < end:
        # check the score to see if it aligns
        if aligned(I): 
            s[i] = s[i-diag] + 1
        else:
            s[i] = max([s[i-j] for j in nn])
        
        # get the next set of characters to compare,
        # incrementing a "non-uniform base number"
        if I[0] < L[0]: I[0] += 1
        else:
            d = 0
            while I[d] == L[d]: 
                I[d] = 1
                d += 1
                if d == N: break
            if d == N: break
            I[d] += 1
        i = getIndex(I)

    #return score matrix for backtracking, s[-1] is the score
    return s
    
def backtrackAll(s,sequences):
    # sequences has to be in the same order as when scoring the alignment
    N = len(sequences)
    L = map(len, sequences)
    start = tuple(L) # starts at the end of all strings
    C = L[:]
    N = len(L)

    for i in range(N): C[i] = L[i] + 1
    for i in range(N-2,0,-1): C[i] = C[i] * C[i+1]
    # the smallest dimension is updated the most...
    
    # previous nearest neighbors array indices, 
    B = C[1:]+[1] # the "counting base"
    # nearest neighbors string relative
    def getCoords(combo):
        coords = N*[0]
        for i in combo: coords[i] = 1
        return tuple(coords)
    def getOffset(combo):
        return sum([B[i] for i in combo])
    nr = range(N)
    ns = {getOffset(combo):getCoords(combo) for i in range(1,N) for combo in it.combinations(nr,i)}
    nn = ns.keys()
    
    # get the index for the element in the score matrix (array)
    def getIndex(index):
        assert len(index) == N      
        ix = 0
        for i in range(N-1):
            ix += C[i+1] * index[i]
        ix += index[N-1]
        return ix
    
    # indexes the characters from each sequence, it is shifted by 1 [1..L[s]]
    def aligned(index):
        cc = sequences[0][index[0]-1] # current character
        for i in range(1,N):
            if sequences[i][index[i]-1] != cc: 
                return False
        return True
    
    def getNeighbor(I, neighbor):
        nb = tuple()
        for i in range(N): nb = nb + (I[i] - neighbor[i],)
        return nb
    
    stack = [(start,"")] # iterative (not recursive), depth first search
    aux = set(stack) # auxiliary/helper for dfs
    lcs = set() # store and return longest common sequences
    found = dict()
    
    while not (len(stack) == 0):
        if len(stack) in found: found[len(stack)] += 1
        else: found[len(stack)] = 1
        I, seq = stack.pop() # pop from stack
        if (min(I) <= 0):
            lcs.add(seq) # no more characters will be added to current sequence, it is one of the longest
        elif aligned(I):
            # it agrees, character part of an lcs
            nextnode = (tuple([i-1 for i in I]),sequences[0][I[0]-1]+seq)
            if (not(nextnode in aux)): 
                aux.add(nextnode)
                stack.append(nextnode)
        else:
            # continue with other longest sequence (they have the best current score)
            j = getIndex(I)
            indices = [j-n for n in nn]
            if min(indices) <= 0: continue
            bestscore = max([s[i] for i in indices])
            for n in nn:
                # put into stack the next best scores
                if s[j-n] == bestscore: 
                    nextnode = (getNeighbor(I,ns[n]),seq)
                    if (not(nextnode in aux)): 
                        aux.add(nextnode)
                        stack.append(nextnode)

    return lcs

#test cases, always must sort in length first
strA = 'ACGATACGT'
strB = 'GCCATTAAGT' 
strC = 'GACTATAGAA'
test1 = [strA, strB, strC]
test1.sort(key = len)
score = getScore(test1)
print score[-1]
print backtrackAll(score,test1)

strA = 'ACGATACGT'
strB = 'CCCATTAAGT' 
strC = 'GACTATAGAA'
test2 = [strA, strB, strC]
test2.sort(key = len)
score = getScore(test2)
print score[-1]
print backtrackAll(score,test2)

strA = 'ACGATACGTACCATTCATT'
strB = 'CCCATTAAGTATA' 
strC = 'GACTATAGAATTGACATT'
strD = 'CATTGTTGAATAGA'
test3 = [strA, strB, strC, strD]
test3.sort(key = len)
score = getScore(test3)
print score[-1]
print backtrackAll(score,test3)