[문자열] 접미사배열(Suffix Array)과 LCP(Longest Common Prefix), Counting Sort 정리

typedef pair<int,int> pii;
vector<pii> sort(vector<pii> &v) {
    int n=v.size();
    vector<pii> ans(n);
    vector<int> cnt(n+1,0);
    vector<int> idx(n,0);
 
    for(int i=0; i<n; i++) cnt[v[i].second]++;
    for(int i=1; i<=n; i++) cnt[i]+=cnt[i-1];
    for(int i=n-1; i>=0; i--) idx[--cnt[v[i].second]]=i;
 
    cnt.clear(); cnt.resize(n+1,0);
    for(int i=0; i<n; i++) cnt[v[i].first]++;
    for(int i=1; i<=n; i++) cnt[i]+=cnt[i-1];
    for(int i=n-1; i>=0; i--) ans[--cnt[v[idx[i]].first]]=v[idx[i]];
 
    return ans;
}

12: cnt배열을 이용해서 v[i].second의 개수를 세준다.
13: 값을 누적시켜준다.
14: idx값은 second 기준으로 정렬했을 때 값이 아닌 인덱스를 저장한다.
만약 second를 기준으로만 정렬하고 말거라면, ans[--cnt[v[i].second]]=v[i].second;일 것이다.
그런데 우리는 값이 아닌 원래 정렬되기 전 상태의 인덱스가 필요하므로 ans[--cnt[v[i].second]]=i;라고 저장한 것이며,이게 사실 답은 아니니까 ans[]배열이 아닌 idx[]배열에 저장한 것 뿐이다.
그럼 14번째 줄까지 하고 for(int i=0; i<n; i++) printf("%d ", idx[i]);를 실행하게 되면,
second기준으로 정렬했을 때 v[]배열의 인덱스를 순서대로 출력하는 것이다. 즉, 이 순서대로 first에 접근하면 second크기 순으로 first에 접근할 수 있으며, 다시 이전과 같은 방법으로 first를 정렬하면 second순서를 이전과 같이 최대한 유지하면서 first를 정렬할 수 있다.

출처: https://plzrun.tistory.com/entry/Counting-Sort-Radix-Sort [plzrun's algorithm]

# Find Longest Common Prefix using suffix array
# first, construct 2 suffix array           -   O(Nlog^N + Mlog^M)
# second,  merge 2 suffix array             -   O(N(a) + M(b))     
a , b  = average of Suffix length compared to each other
# second, find LCP                          -

import time
import math

def suffix_array(txt):
    if not txt:
        return []
    txt += chr(0)
    N = len(txt)
    #print "input string = [", txt, "]"

    equivalence = {t: i for i, t in enumerate(sorted(set(txt)))}
    #print "equi=", equivalence                          #bucket number set
    cls = [equivalence[t] for t in txt]
    #print "cls=", cls                                   # first index sorting

    ns = [(2**i)%N for i in range( int( math.ceil( math.log(N,2))))]
    #print "ns=", ns

    for n in ns[:-1]:
        result = sorted(zip(cls, cls[n:]+cls[:n], range(N)))
        result0, result1, inds = list(zip(*result))

        cls[inds[0]] = 0
        for j in range(1, N):
            cls[inds[j]] = cls[inds[j-1]]
            if (result0[j], result1[j]) != (result0[j-1], result1[j-1]):
                cls[inds[j]] += 1
    n = ns[-1]
    result = sorted(zip(cls, cls[n:]+cls[:n], range(N)))
    return list(list(zip(*result))[2][1:])

def merge(sufarr1, sufarr2 , sufArrIndex1 , sufArrIndex2):      # sort(sumLa)  << time check need
    indMerge = []
    i, j = 0, 0
    len1, len2 = len(sufarr1), len(sufarr2)
    while True:
        if i == len1 and j == len2:             break
        if i == len1:           # if finished sufarr1
            indMerge.append(-sufArrIndex2[j]-1)
            j += 1
        elif j == len2:         # if finished sufarr2
            indMerge.append(sufArrIndex1[i]+1)
            i += 1
        else:
            index1, index2 = sufArrIndex1[i], sufArrIndex2[j]
            str1 , str2 = sufarr1[index1:] , sufarr2[index2:]
            if str1 < str2 :              # standard character-by-character comparison rules
                indMerge.append(index1+1)
                i += 1
            elif str1 > str2:
                indMerge.append(-index2-1)
                j += 1
            else:
                indMerge.append(index1+1)
                indMerge.append(-index2-1)
                i += 1
                j += 1
    return indMerge

def getLCP(suffix1, suffix2, indMerge):
    lcp, where = [0] , []
    length = len(indMerge)
    position = 0 if indMerge[0] > 0 else 1
    where.append(position)
    for x in range(length-1):
        index = indMerge[x]
        if indMerge[x] > 0:
            str1 = suffix1[index - 1:]
        else:
            str1 = suffix2[-index-1:]
        index = indMerge[x+1]
        if indMerge[x+1] > 0:
            str2 = suffix1[index - 1:]
            where.append(0)
        else:
            str2 = suffix2[-index-1:]
            where.append(1)
        size = min(len(str1), len(str2))
        count = 0
        for y in range(size):
            if str1[y] == str2[y]:
                count += 1
            else:
                break
        lcp.append(count)
    return lcp, where

def lcped(suffixMerge,seq1,seq2):
    rank , lcp = [] , [0]
    N = len(suffixMerge)
    for i in range(N):
        rank.insert(suffixMerge[i],i)
    length = 0
    for i in range(N):
        k = rank[i]
        if k:
            j = suffixMerge[k-1]
            while True:
                if j == N :
                    break
                if j+length < N / 2:
                    str1 = seq1[j + length]
                elif j + length < N:
                    str1 = seq2[j + length - N / 2]
                else:
                    break
                if i+length < N / 2:
                    str2 = seq1[i + length]
                elif i + length < N:
                        str2 = seq2[i + length - N / 2]
                else:
                    break
                #str1 = seq1[j + length] if j < N / 2 else seq2[j + length - N / 2]
                #str2 = seq1[i + length] if i < N / 2 else seq2[i + length - N / 2]
                if str1 == str2:
                    length+=1
                else:
                    break
            lcp.insert(k,length)
            if length:
                length-=1
        i+=1
    return lcp

def readSequnce(sequence):
    f = open(sequence)
    for line in f:
        seq = line
    f.close()
    return seq

def main(sequence1,sequence2):
    seq1 = readSequnce(sequence1)
    seq2 = readSequnce(sequence2)
    start_time = time.time()
    la1= suffix_array(seq1)
    la2= suffix_array(seq2)
    print len(la1) , len(la2)
    suffixMerge = merge(seq1,seq2,la1,la2)
    lcp, position = getLCP(seq1,seq2,suffixMerge)
    maxValue, maxStr, compare = lcp[0],'',position[0]
    N = len(suffixMerge)
    for i in range(N):
        if i==0: continue
        if lcp[i] > maxValue and position[i] != position[i-1]:
            maxValue = lcp[i]
            maxStr = seq1[suffixMerge[i]-1:suffixMerge[i]+maxValue-1] if suffixMerge[i] > 0 else seq2[-suffixMerge[i]-1: -suffixMerge[i] + maxValue-1]
    print("--- %s seconds ---" % (time.time() - start_time))
    print maxStr, maxValue

if __name__ == '__main__':
    main('seq819200_1.txt','seq819200_2.txt')