后缀数组的自底向上（bottom-up）遍历算法
后缀数组自底向上遍历等价于后缀树的自底向上遍历。由于后缀数组不是树型结构，在遍历时除了SA本身之外还需要额外的信息，这时Suffix Array就是一个增强的后缀数组(Enhanced Suffix Array)了。该算法使用后缀数组的一个增强信息---LCP表，并通过堆栈模拟自底向上的遍历。遍历的结果就是一颗虚拟的lcp-interval 树，其中每一个结点对应后缀树的一个内部结点。有些应用中，遍历时需要知道每个结点的孩子信息，因此在下面的实现中提供了两个版本 bottomUpTraverseWithoutChildren和bottomUpTraverseWithChildren。
需要说明 的是，树中每一个lcp-interval结点表示为：lcp-[i..j]，其中lcp相当于后缀树的pathlen，i和j分别是以该lcp- interval结点为根结点表示的子树中的最小和最大后缀数组下标，例如mississippi$的一个lcp-interval结点是 1-[1..4]，表示后缀数组中第1个后缀到第4个后缀这总共四个后缀作为该lcp-interval结点为根的子树的四个叶子结点。
实现：
import java.util.ArrayList;
import java.util.List;
import java.util.Stack;
/**
*
* Bottom-Up Traversal of a suffix array (with LCP table)
* (The suffix array is constructed with prefix doubling algorithm)
*
*
* Copyright (c) 2011 ljs (http://blog.csdn.net/ljsspace/)
* Licensed under GPL (http://www.opensource.org/licenses/gpl-license.php)
*
* @author ljs
* 2011-07-23
*
*/
public class ESA_BottomUpTraversal {
public static final char MAX_CHAR = '\u00FF';
private String text;
private int[] suffixarray;
private int[] ranktable;
private int[] lcptable;
public ESA_BottomUpTraversal(String text){
this.text = text;
}
class Suffix{
int[] sa;
//Note: the p-th suffix in sa: SA[rank[p]-1]];
//p is the index of the array "rank", start with 0;
//a text S's p-th suffix is S[p..n], n=S.length-1.
int[] rank;
boolean done;
public Suffix(int[] sa,int[] rank){
this.sa = sa;
this.rank = rank;
}
}
//a prefix of suffix[isuffix] represented with digits
class Tuple{
int isuffix; //the p-th suffix
int[] digits;
public Tuple(int suffix,int[] digits){
this.isuffix = suffix;
this.digits = digits;
}
public String toString(){
StringBuffer sb = new StringBuffer();
sb.append(isuffix);
sb.append("(");
for(int i=0;i<digits.length;i++){
sb.append(digits[i]);
if(i<digits.length-1)
sb.append("-");
}
sb.append(")");
return sb.toString();
}
}
//d: the digit to do countingsort
//max: A value's range is 0...max
private void countingSort(int d,Tuple[] tA,Tuple[] tB,int max){
//init the counter array
int[] C = new int[max+1];
for(int i=0;i<=max;i++){
C[i] = 0;
}
//stat the count
for(int j=0;j<tA.length;j++){
C[tA[j].digits[d]]++;
}
//process the counter array C
for(int i=1;i<=max;i++){
C[i]+=C[i-1];
}
//distribute the values
for(int j=tA.length-1;j>=0;j--){
//C[A[j]] <= A.length
tB[--C[tA[j].digits[d]]]=tA[j];
}
}
//tA: input
//tB: output for rank caculation
private void radixSort(Tuple[] tA,Tuple[] tB,int max,int digitsLen){
int len = tA.length;
int digitsTotalLen = tA[0].digits.length;
for(int d=digitsTotalLen-1,j=0;j<digitsLen;d--,j++){
this.countingSort(d, tA, tB, max);
//assign tB to tA
if(j<digitsLen-1){
for(int i=0;i<len;i++){
tA[i] = tB[i];
}
}
}
}
//max is the maximum value in any digit of TA.digits[], used for counting sort
//tA: input
//tB: the place holder, reused between iterations
private Suffix rank(Tuple[] tA,Tuple[] tB,int max,int digitsLen){
int len = tA.length;
radixSort(tA,tB,max,digitsLen);
int digitsTotalLen = tA[0].digits.length;
//caculate rank and sa
int[] sa = new int[len];
sa[0] = tB[0].isuffix;
int[] rank = new int[len+2]; //add 2 for sentinel
rank[len]=1;rank[len+1] = 1;
int r = 1; //rank starts with 1
rank[tB[0].isuffix] = r;
for(int i=1;i<len;i++){
sa[i] = tB[i].isuffix;
boolean equalLast = true;
for(int j=digitsTotalLen-digitsLen;j<digitsTotalLen;j++){
if(tB[i].digits[j]!=tB[i-1].digits[j]){
equalLast = false;
break;
}
}
if(!equalLast){
r++;
}
rank[tB[i].isuffix] = r;
}
Suffix suffix = new Suffix(sa,rank);
//judge if we are done
if(r==len){
suffix.done = true;
}else{
suffix.done = false;
}
return suffix;
}
private int[] orderSuffixes(Tuple[] tA,Tuple[] tB,int max,int digitsLen){
int len = tA.length;
radixSort(tA,tB,max,digitsLen);
//caculate rank and sa
int[] sa = new int[len];
for(int i=0;i<len;i++){
sa[i] = tB[i].isuffix;
}
return sa;
}
//rank needs sentinel: len+2
private Suffix reduce(int[] rank,int max){
int len = rank.length - 2;
int n1 = (len+1)/3;
int n2 = len/3;
Tuple[] tA = new Tuple[n1+n2];
Tuple[] tB = new Tuple[n1+n2];
for(int i=0,j=1;i<n1;i++,j+=3){
int r1 = rank[j];
int r2 = rank[j+1];
int r3 = rank[j+2];
tA[i] = new Tuple(i,new int[]{r1,r2,r3});
}
for(int i=n1,j=2;i<n1+n2;i++,j+=3){
int r1 = rank[j];
int r2 = rank[j+1];
int r3 = rank[j+2];
tA[i] = new Tuple(i,new int[]{r1,r2,r3});
}
return rank(tA,tB,max,3);
}
private int[] skew(int[] rank,int max){
int len = rank.length - 2;
//step 1: caculate sa12
Suffix suffixT12 = reduce(rank,max);
int[] sa12 = null;
if(!suffixT12.done){
int[] rankT12 = suffixT12.rank;
int maxT12 = rankT12[suffixT12.sa[suffixT12.sa.length-1]];
sa12 = skew(rankT12,maxT12);
// debug for string: GACCCACCACC#
//s12 = new Suffix();
//s12.rank = new int[]{3,6,5,4,7,2,1,1,1};
//s12.sa = new int[]{7,6,5,0,3,2,1,4};
//s12.done =true;
}else{
sa12 = suffixT12.sa;
}
//index conversion for sa12
int n1 = (len+1)/3;
for(int j=0;j<sa12.length;j++){
if(sa12[j]<n1){
sa12[j] = 1 + 3*sa12[j];
}else{
sa12[j] = 2 + 3*(sa12[j]-n1);
}
}
//recaculate rank for sa12
int[] rank12 = new int[len+2];
rank12[len] = 1;rank12[len+1] = 1;
for(int k=0;k<sa12.length;k++){
rank12[sa12[k]] = k+1;
}
//step 2: caculate sa0
int n0=(len+2)/3;
Tuple[] tA = new Tuple[n0];
Tuple[] tB = new Tuple[n0];
for(int i=0,j=0;i<n0;i++,j+=3){
int r1 = rank[j];
int r2 = rank12[j+1];
tA[i] = new Tuple(i,new int[]{r1,r2});
}
int max12 = rank12[sa12[sa12.length-1]];
int[] sa0 = orderSuffixes(tA,tB,max<max12?max12:max,2);
//index conversion for sa0
for(int j=0;j<n0;j++){
sa0[j] = 3*sa0[j];
}
//step 3: merge sa12 and sa0
int[] sa = new int[len];
int i=0,j=0;
int k=0;
while(i<sa12.length && j<sa0.length){
int p = sa12[i];
int q = sa0[j];
if(p%3==1){
//case 1
if(rank[p]<rank[q]){
sa[k++] = p;i++;
}else if(rank[p]>rank[q]){
sa[k++] = q;j++;
}else{
if(rank12[p+1]<rank12[q+1]){
sa[k++] = p;i++;
}else{
sa[k++] = q;j++;
}
}
}else{
//case 2
if(rank[p]<rank[q]){
sa[k++] = p;i++;
}else if(rank[p]>rank[q]){
sa[k++] = q;j++;
}else{
if(rank[p+1]<rank[q+1]){
sa[k++] = p;i++;
}else if(rank[p+1]>rank[q+1]){
sa[k++] = q;j++;
}else{
if(rank12[p+2]<rank12[q+2]){
sa[k++] = p;i++;
}else{
sa[k++] = q;j++;
}
}
}
}
}
for(int m=i;m<sa12.length;m++){
sa[k++] = sa12[m];
}
for(int m=j;m<sa0.length;m++){
sa[k++] = sa0[m];
}
return sa;
}
//Precondition: the last char in text must be less than other chars.
public void computeSuffixArray(){
if(text == null)return;
int len = text.length();
if(len == 0) return ;
char base = text.charAt(len-1); //the smallest char
Tuple[] tA = new Tuple[len];
Tuple[] tB = new Tuple[len]; //placeholder
for(int i=0;i<len;i++){
tA[i] = new Tuple(i,new int[]{0,text.charAt(i)-base});
}
Suffix suffix = rank(tA,tB,MAX_CHAR-base,1);
int max = suffix.rank[suffix.sa[len-1]];
int[] sa = skew(suffix.rank,max);
//caculate rank for result suffix array
int[] r = new int[len];
for(int k=0;k<sa.length;k++){
r[sa[k]] = k+1;
}
this.suffixarray = sa;
this.ranktable = r;
}
//rank[p]'s index starts with 1 (not 0)
public void computeLCPtable(){
if(text == null)return;
int len = text.length();
if(len == 0) return ;
int[] sa = this.suffixarray;
int[] rank = this.ranktable;
int[] lcpz = new int[len];
//base case: p=0
//caculate LCP of suffix[0]
int lcp = 0;
int r = rank[0]-1;
if(r>0){
int q=sa[r-1];
//caculate LCP by definition
for(int i=0,j=q;i<len && j<len;i++,j++){
if(text.charAt(i) != text.charAt(j)){
lcp=i;
break;
}
}
}
lcpz[0] = lcp;
//other cases: p>=1
//ignore p == sa[0] because LCP=0 for suffix[p] where rank[p]=0
for(int p=1;p<len && p != sa[0];p++){
int h = lcpz[p-1];
int q=sa[rank[p]-2];
lcp = 0;
if(h>1){ //for h<=1, caculate LCP by definition (i.e. start with lcp=0)
//jump h-1 chars for suffix[p] and suffix[q]
lcp = h-1;
}
for(int i=p+lcp,j=q+lcp,k=0;i<len && j<len;i++,j++,k++){
if(text.charAt(i) != text.charAt(j)){
lcp+=k;
break;
}
}
lcpz[p] = lcp;
}
//caculate LCP
int[] LCP = new int[len];
for(int i=0;i<len;i++){
LCP[i] = lcpz[sa[i]];
}
this.lcptable = LCP;
}
private void reportLCP(){
System.out.format("Text: %s%n",text);
int len = this.text.length();
System.out.println("suffix array:");
for(int i=0;i<len;i++){
System.out.format(" %s", this.suffixarray[i]);
}
System.out.println();
System.out.println("rank table:");
for(int i=0;i<len;i++){
System.out.format(" %s", this.ranktable[i]);
}
System.out.println();
System.out.println("lcp table:");
for(int i=0;i<len;i++){
System.out.format(" %s", this.lcptable[i]);
}
}
class LCPInterval{
int lcp; //the lcp value of the lcp-interval
int lb; //the left boundary suffix index
int rb; //the right boundary suffix index
List<LCPInterval> children = new ArrayList<LCPInterval>();
public LCPInterval(int lcp,int lb,int rb){
this.lcp = lcp;
this.lb = lb;
this.rb = rb;
}
public String toString(){
return String.format("%d-[%d..%d]",
this.lcp,this.lb,this.rb);
}
}
private void reportLCPInterval(LCPInterval interval){
if(interval.children.size()>0){
StringBuilder sb = new StringBuilder();
for(LCPInterval child:interval.children){
sb.append(child.toString());
sb.append(",");
}
sb.deleteCharAt(sb.length()-1);
System.out.format("%s, children={%s}%n",
interval,sb.toString());
}else{
System.out.format("%s%n", interval);
}
}
//traverse the corresponding suffix tree with a bottom-up approach
//each internal node is equivalent to an lcp-interval.
public void bottomUpTraverseWithoutChildren(){
int len = text.length();
Stack<LCPInterval> stack = new Stack<LCPInterval>();
int lb = -1;
//push root's first child
stack.push(new LCPInterval(0,0,-1));
for(int i=1;i<len;i++){
//save lb
lb = i - 1;
//pop all the deeper suffixes
while(lcptable[i]<stack.peek().lcp){
LCPInterval interval = stack.pop();
//update the popped interval's rb
interval.rb = i-1;
reportLCPInterval(interval);
//update lb
lb = interval.lb;
}
if(lcptable[i]>stack.peek().lcp){
stack.push(new LCPInterval(lcptable[i],lb,-1));
} //if lcptable[i]==interval.lcp, no push because rb is updated when popped
}
while(!stack.isEmpty()){
LCPInterval interval = stack.pop();
//update the popped interval's rb
interval.rb = len-1;
reportLCPInterval(interval);
}
}
public void bottomUpTraverseWithChildren(){
int len = text.length();
Stack<LCPInterval> stack = new Stack<LCPInterval>();
int lb = -1;
LCPInterval lastInterval = null;
//push root's first child
stack.push(new LCPInterval(0,0,-1));
for(int i=1;i<len;i++){
//save lb
lb = i - 1;
//pop all the deeper suffixes
while(lcptable[i]<stack.peek().lcp){
lastInterval = stack.pop();
//update the popped interval's rb
lastInterval.rb = i-1;
//Note: when lastInterval is popped out, if it has children then its
//children are already known!
reportLCPInterval(lastInterval);
//update lb
lb = lastInterval.lb;
if(lcptable[i]<=stack.peek().lcp){
//case 1: top>next>i, then top is next's child
LCPInterval next = stack.peek();
next.children.add(lastInterval);
lastInterval = null;
}
}
if(lcptable[i]>stack.peek().lcp){
LCPInterval curr=new LCPInterval(lcptable[i],lb,-1);
if(lastInterval != null){
//case 2: top>i>next, then top is i's child
curr.children.add(lastInterval);
lastInterval = null;
}
stack.push(curr);
} //if lcptable[i]==interval.lcp, no push because rb is updated when popped
}
while(!stack.isEmpty()){
lastInterval = stack.pop();
//update the popped interval's rb
lastInterval.rb = len-1;
reportLCPInterval(lastInterval);
if(!stack.isEmpty()){//case 1: top > next, i.e. top is next's child
LCPInterval next = stack.peek();
next.children.add(lastInterval);
}
}
}
public void solve(){
this.computeSuffixArray();
this.computeLCPtable();
this.reportLCP();
System.out.format("%nbottom-up traversal with no children list: %n");
this.bottomUpTraverseWithoutChildren();
System.out.format("%nbottom-up traversal with children list: %n");
this.bottomUpTraverseWithChildren();
}
public static void main(String[] args) {
String text = "mississippi#";
ESA_BottomUpTraversal esa = new ESA_BottomUpTraversal(text);
esa.solve();
System.out.format("%n********************************%n");
text = "GACCCACCACC#";
esa = new ESA_BottomUpTraversal(text);
esa.solve();
System.out.format("%n********************************%n");
}
}
测试：
Text: mississippi#
suffix array:
11 10 7 4 1 0 9 8 6 3 5 2
rank table:
6 5 12 10 4 11 9 3 8 7 2 1
lcp table:
0 0 1 1 4 0 0 1 0 2 1 3
bottom-up traversal with no children list:
4-[3..4]
1-[1..4]
1-[6..7]
2-[8..9]
3-[10..11]
1-[8..11]
0-[0..11]
bottom-up traversal with children list:
4-[3..4]
1-[1..4], children={4-[3..4]}
1-[6..7]
2-[8..9]
3-[10..11]
1-[8..11], children={2-[8..9],3-[10..11]}
0-[0..11], children={1-[1..4],1-[6..7],1-[8..11]}
********************************
Text: GACCCACCACC#
suffix array:
11 8 5 1 10 7 4 9 6 3 2 0
rank table:
12 4 11 10 7 3 9 6 2 8 5 1
lcp table:
0 0 3 3 0 1 4 1 2 5 2 0
bottom-up traversal with no children list:
3-[1..3]
4-[5..6]
5-[8..9]
2-[7..10]
1-[4..10]
0-[0..11]
bottom-up traversal with children list:
3-[1..3]
4-[5..6]
5-[8..9]
2-[7..10], children={5-[8..9]}
1-[4..10], children={4-[5..6],2-[7..10]}
0-[0..11], children={3-[1..3],1-[4..10]}
********************************
参考：
Mohamed Ibrahim Abouelhoda, Stefan Kurtz, Enno Ohlebusch: Replacing suffix tree with enhanced suffix arrays (2004)
