I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

Ukkonen’s suffix tree construction

From Rosetta Code
Ukkonen’s suffix tree construction is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Suffix Trees are very useful in numerous string processing and computational biology problems.

The task is to create a function which implements Ukkonen’s algorithm to create a useful Suffix Tree as described:

Part 1
Part 2
Part 3
Part 4
Part 5
Part 6

Using Arithmetic-geometric mean/Calculate Pi generate the first 1000, 10000, and 100000 decimal places of pi. Using your implementation with an alphabet of 0 through 9 (plus $ say to make the tree explicit) find the longest repeated string in each list. Time your results and demonstrate that your implementation is linear (i.e. that 10000 takes approx. 10 times as long as 1000). You may vary the size of the lists of decimal places of pi to give reasonable answers.

Go[edit]

This is a translation of the C code here which is an extended form of the code in Part 6 of the task description for finding the longest repeated substring of a given string. In the interests of brevity, the extensive comments in the C version have been largely omitted. The C code doesn't compile as it stands but I have added a fix in the Talk Page.

For convenience I have included the code from the Arithmetic-geometric_mean/Calculate_Pi#Go task in the same package.

It takes around 25 seconds on my machine (Celeron @1.6GHz) to calculate the first 100,000 (or so) decimal places of Pi. Having done that, the timings for extracting the longest repeated sequence of digits are quick and fairly linear as expected.

As the task doesn't say whether overlapping sequences are to be counted, I've assumed that they are as this is what the algorithm naturally produces.

package main
 
import (
"fmt"
"math/big"
"time"
)
 
var maxChar = 128
 
type Node struct {
children []*Node
suffixLink *Node
start int
end *int
suffixIndex int
}
 
var (
text string
root *Node
lastNewNode *Node
activeNode *Node
activeEdge = -1
activeLength = 0
remainingSuffixCount = 0
leafEnd = -1
rootEnd *int
splitEnd *int
size = -1
)
 
func newNode(start int, end *int) *Node {
node := new(Node)
node.children = make([]*Node, maxChar)
node.suffixLink = root
node.start = start
node.end = end
node.suffixIndex = -1
return node
}
 
func edgeLength(n *Node) int {
if n == root {
return 0
}
return *(n.end) - n.start + 1
}
 
func walkDown(currNode *Node) bool {
if activeLength >= edgeLength(currNode) {
activeEdge += edgeLength(currNode)
activeLength -= edgeLength(currNode)
activeNode = currNode
return true
}
return false
}
 
func extendSuffixTree(pos int) {
leafEnd = pos
remainingSuffixCount++
lastNewNode = nil
for remainingSuffixCount > 0 {
if activeLength == 0 {
activeEdge = pos
}
if activeNode.children[text[activeEdge]] == nil {
activeNode.children[text[activeEdge]] = newNode(pos, &leafEnd)
if lastNewNode != nil {
lastNewNode.suffixLink = activeNode
lastNewNode = nil
}
} else {
next := activeNode.children[text[activeEdge]]
if walkDown(next) {
continue
}
if text[next.start+activeLength] == text[pos] {
if lastNewNode != nil && activeNode != root {
lastNewNode.suffixLink = activeNode
lastNewNode = nil
}
activeLength++
break
}
temp := next.start + activeLength - 1
splitEnd = &temp
split := newNode(next.start, splitEnd)
activeNode.children[text[activeEdge]] = split
split.children[text[pos]] = newNode(pos, &leafEnd)
next.start += activeLength
split.children[text[next.start]] = next
if lastNewNode != nil {
lastNewNode.suffixLink = split
}
lastNewNode = split
}
remainingSuffixCount--
if activeNode == root && activeLength > 0 {
activeLength--
activeEdge = pos - remainingSuffixCount + 1
} else if activeNode != root {
activeNode = activeNode.suffixLink
}
}
}
 
func setSuffixIndexByDFS(n *Node, labelHeight int) {
if n == nil {
return
}
if n.start != -1 {
// Uncomment line below to print suffix tree
// fmt.Print(text[n.start: *(n.end) +1])
}
leaf := 1
for i := 0; i < maxChar; i++ {
if n.children[i] != nil {
// Uncomment the 3 lines below to print suffix index
//if leaf == 1 && n.start != -1 {
// fmt.Printf(" [%d]\n", n.suffixIndex)
//}
leaf = 0
setSuffixIndexByDFS(n.children[i], labelHeight+edgeLength(n.children[i]))
}
}
if leaf == 1 {
n.suffixIndex = size - labelHeight
// Uncomment line below to print suffix index
//fmt.Printf(" [%d]\n", n.suffixIndex)
}
}
 
func buildSuffixTree() {
size = len(text)
temp := -1
rootEnd = &temp
root = newNode(-1, rootEnd)
activeNode = root
for i := 0; i < size; i++ {
extendSuffixTree(i)
}
labelHeight := 0
setSuffixIndexByDFS(root, labelHeight)
}
 
func doTraversal(n *Node, labelHeight int, maxHeight, substringStartIndex *int) {
if n == nil {
return
}
if n.suffixIndex == -1 {
for i := 0; i < maxChar; i++ {
if n.children[i] != nil {
doTraversal(n.children[i], labelHeight+edgeLength(n.children[i]),
maxHeight, substringStartIndex)
}
}
} else if n.suffixIndex > -1 && (*maxHeight < labelHeight-edgeLength(n)) {
*maxHeight = labelHeight - edgeLength(n)
*substringStartIndex = n.suffixIndex
}
}
 
func getLongestRepeatedSubstring(s string) {
maxHeight := 0
substringStartIndex := 0
doTraversal(root, 0, &maxHeight, &substringStartIndex)
// Uncomment line below to print maxHeight and substringStartIndex
// fmt.Printf("maxHeight %d, substringStartIndex %d\n", maxHeight, substringStartIndex)
if s == "" {
fmt.Printf("  %s is: ", text)
} else {
fmt.Printf("  %s is: ", s)
}
k := 0
for ; k < maxHeight; k++ {
fmt.Printf("%c", text[k+substringStartIndex])
}
if k == 0 {
fmt.Print("No repeated substring")
}
fmt.Println()
}
 
func calculatePi() *big.Float {
one := big.NewFloat(1)
two := big.NewFloat(2)
four := big.NewFloat(4)
prec := uint(325 * 1024) // enough to calculate Pi to 100,182 decimal digits
 
a := big.NewFloat(1).SetPrec(prec)
g := new(big.Float).SetPrec(prec)
 
// temporary variables
t := new(big.Float).SetPrec(prec)
u := new(big.Float).SetPrec(prec)
 
g.Quo(a, t.Sqrt(two))
sum := new(big.Float)
pow := big.NewFloat(2)
 
for a.Cmp(g) != 0 {
t.Add(a, g)
t.Quo(t, two)
g.Sqrt(u.Mul(a, g))
a.Set(t)
pow.Mul(pow, two)
t.Sub(t.Mul(a, a), u.Mul(g, g))
sum.Add(sum, t.Mul(t, pow))
}
 
t.Mul(a, a)
t.Mul(t, four)
pi := t.Quo(t, u.Sub(one, sum))
return pi
}
 
func main() {
tests := []string{
"GEEKSFORGEEKS$",
"AAAAAAAAAA$",
"ABCDEFG$",
"ABABABA$",
"ATCGATCGA$",
"banana$",
"abcpqrabpqpq$",
"pqrpqpqabab$",
}
fmt.Println("Longest Repeated Substring in:\n")
for _, test := range tests {
text = test
buildSuffixTree()
getLongestRepeatedSubstring("")
}
fmt.Println()
 
pi := calculatePi()
piStr := fmt.Sprintf("%v", pi)
piStr = piStr[2:] // remove initial 3.
numbers := []int{1e3, 1e4, 1e5}
maxChar = 58
for _, number := range numbers {
start := time.Now()
text = piStr[0:number] + "$"
buildSuffixTree()
getLongestRepeatedSubstring(fmt.Sprintf("first %d d.p. of Pi", number))
elapsed := time.Now().Sub(start)
fmt.Printf(" (this took %s)\n\n", elapsed)
}
}
Output:

Sample run:

Longest Repeated Substring in:

  GEEKSFORGEEKS$ is: GEEKS
  AAAAAAAAAA$ is: AAAAAAAAA
  ABCDEFG$ is: No repeated substring
  ABABABA$ is: ABABA
  ATCGATCGA$ is: ATCGA
  banana$ is: ana
  abcpqrabpqpq$ is: ab
  pqrpqpqabab$ is: ab

  first 1000 d.p. of Pi is: 23846
  (this took 7.728858ms)

  first 10000 d.p. of Pi is: 7111369
  (this took 57.524478ms)

  first 100000 d.p. of Pi is: 041021944
  (this took 599.770281ms)

Julia[edit]

Translation of: Go

Uses array indices instead of the Go version's node pointers.

const oo = typemax(Int)
 
"""The suffix-tree's node."""
mutable struct Node
children::Dict{Char, Int}
start::Int
ending::Int
suffixlink::Int
suffixindex::Int
end
 
Node() = Node(Dict(), 0, oo, 0, -1)
Node(start, ending) = Node(Dict(), start, ending, 0, -1)
 
""" Ukkonen Suffix-Tree """
mutable struct SuffixTree
nodes::Vector{Node}
text::Vector{Char}
root::Int
position::Int
currentnode::Int
needsuffixlink::Int
remainder::Int
activenode::Int
activelength::Int
activeedge::Int
end
 
edgelength(st, n::Node) = min(n.ending, st.position + 1) - n.start
 
function newnode(st, start, ending)
st.currentnode += 1
st.nodes[st.currentnode] = Node(start, ending)
return st.currentnode
end
 
function SuffixTree(str::String)
nodes = [Node() for _ in 1:length(str) * 2]
st = SuffixTree(nodes, [c for c in str], 1, 0, 0, 0, 0, 1, 1, 1)
st.root = newnode(st, 0, 0)
st.activenode = st.root
for i in 1:length(st.text)
extendsuffixtree(st, i)
end
setsuffixindexbyDFS(st, st.nodes[st.root], 0)
return st
end
 
function addsuffixlink(st, nodenum::Int)
if st.needsuffixlink > 0
st.nodes[st.needsuffixlink].suffixlink = nodenum
end
st.needsuffixlink = nodenum
end
 
activeedge(st) = st.text[st.activeedge]
 
function walkdown!(st, currnode::Int)
len = edgelength(st, st.nodes[currnode])
st.activelength < len && return false
st.activeedge += len
st.activelength -= len
st.activenode = currnode
return true
end
 
function extendsuffixtree(st, pos)
st.position = pos
st.needsuffixlink = 0
st.remainder += 1
while st.remainder > 0
st.activelength == 0 && (st.activeedge = st.position)
if !haskey(st.nodes[st.activenode].children, activeedge(st))
nodenum = newnode(st, st.position, oo)
st.nodes[st.activenode].children[activeedge(st)] = nodenum
addsuffixlink(st, st.activenode)
else
next = st.nodes[st.activenode].children[activeedge(st)]
walkdown!(st, next) && continue
if st.text[st.nodes[next].start + st.activelength] == st.text[pos]
addsuffixlink(st, st.activenode)
st.activelength += 1
break
end
splt = newnode(st, st.nodes[next].start, st.nodes[next].start + st.activelength)
st.nodes[st.activenode].children[activeedge(st)] = splt
nodenum = newnode(st, st.position, oo)
st.nodes[splt].children[st.text[pos]] = nodenum
st.nodes[next].start += st.activelength
st.nodes[splt].children[st.text[st.nodes[next].start]] = next
addsuffixlink(st, splt)
end
st.remainder -= 1
if st.activenode == st.root && st.activelength > 0
st.activelength -= 1
st.activeedge = st.position - st.remainder + 1
elseif st.activenode != st.root
st.activenode = st.nodes[st.activenode].suffixlink
end
end
end
 
function setsuffixindexbyDFS(st, node, labelheight, verbose=false)
verbose && node.start > 0 && print(st.text[node.start:min(node.ending, length(st.text))])
isleaf = true
for child in map(v -> st.nodes[v], collect(values(node.children)))
verbose && isleaf && node.start > 0 && println(" [", node.suffixindex, "]")
isleaf = false
setsuffixindexbyDFS(st, child, labelheight + edgelength(st, child))
end
if isleaf
idx = length(st.text) - labelheight
node.suffixindex = idx
verbose && println(" [$idx]")
end
end
 
function dotraversal(st)
maxheight, substringstartindices = 0, [0]
function traversal(node::Node, labelheight)
if node.suffixindex == -1
for child in map(v -> st.nodes[v], collect(values(node.children)))
traversal(child, labelheight + edgelength(st, child))
end
elseif maxheight < labelheight - edgelength(st, node)
maxheight = labelheight - edgelength(st, node)
substringstartindices = [node.suffixindex + 1]
elseif maxheight == labelheight - edgelength(st, node)
push!(substringstartindices, node.suffixindex + 1)
end
end
traversal(st.nodes[st.root], 0)
return maxheight, substringstartindices
end
 
function getlongestrepeatedsubstring(st::SuffixTree, label="", printresult=true)
len, starts = dotraversal(st)
substring = len == 0 ? "" :
join(unique(map(x -> String(st.text[x:x+len-1]), starts)), " (or) ")
if printresult
print(" ", label == "" ? String(st.text) : label, ": ")
println(len == 0 ? "No repeated substring." : substring)
end
return substring
end
 
function testsuffixtree()
tests = [
"CAAAABAAAABD\$",
"GEEKSFORGEEKS\$",
"AAAAAAAAAA\$",
"ABCDEFG\$",
"ABABABA\$",
"ATCGATCGA\$",
"banana\$",
"abcpqrabpqpq\$",
"pqrpqpqabab\$",
]
println("Longest Repeated Substring in:\n")
for test in tests
st = SuffixTree(test)
getlongestrepeatedsubstring(st)
end
println()
sπ = ""
setprecision(4000000) do
sπ = string(BigFloat(π))[3:end]
end
for number in [1000, 10000, 100000, 1000000]
text = sπ[1:number] * "\$"
@time begin
st = SuffixTree(text)
getlongestrepeatedsubstring(st, "first $number d.p. of π")
end
end
end
 
testsuffixtree()
 
Output:
Longest Repeated Substring in:

  CAAAABAAAABD: AAAAB
  GEEKSFORGEEKS: GEEKS
  AAAAAAAAAA: AAAAAAAAA
  ABCDEFG: No repeated substring.
  ABABABA: ABABA
  ATCGATCGA: ATCGA
  banana: ana
  abcpqrabpqpq: ab (or) pq
  pqrpqpqabab: ab (or) pq

  first 1000 d.p. of π: 60943 (or) 42019 (or) 82534 (or) 99999 (or) 93751 (or) 23846 (or) 33446
  0.003336 seconds (34.86 k allocations: 4.252 MiB)
  first 10000 d.p. of π: 7111369 (or) 8530614
  0.038749 seconds (351.60 k allocations: 42.460 MiB, 16.54% gc time)
  first 100000 d.p. of π: 134926158 (or) 041021944 (or) 201890888
  0.533892 seconds (3.52 M allocations: 425.035 MiB, 22.42% gc time)
  first 1000000 d.p. of π: 756130190263
  6.008879 seconds (35.25 M allocations: 4.152 GiB, 23.20% gc time)

Phix[edit]

Translation of: Go
-- demo/rosetta/Ukkonens_Suffix_Tree.exw
integer maxChar = 'z'
 
sequence children = {},
suffixLinks = {},
starts = {},
ends = {},
suffixIndices = {}
string text
atom pSplitEnd,
pRootEnd,
pLeafEnd = allocate_word(0,true)
integer root = NULL,
lastNewNode,
activeNode,
activeEdge = -1,
activeLength = 0,
remainingSuffixCount = 0,
size = -1
 
function newNode(integer start, atom pFinish, bool bKids=false)
children = append(children,iff(bKids?repeat(NULL,maxChar):0))
suffixLinks = append(suffixLinks,root)
starts = append(starts,start)
ends = append(ends,pFinish)
suffixIndices = append(suffixIndices,-1)
return length(children)
end function
 
function edgeLength(integer n)
if n == root then
return 0
end if
return peekns(ends[n]) - starts[n] + 1
end function
 
function walkDown(integer currNode)
integer l = edgeLength(currNode)
if activeLength >= l then
activeEdge += l
activeLength -= l
activeNode = currNode
return true
end if
return false
end function
 
procedure extendSuffixTree(integer pos)
poken(pLeafEnd,pos)
remainingSuffixCount += 1
lastNewNode = NULL
while remainingSuffixCount > 0 do
if activeLength == 0 then
activeEdge = pos
end if
integer ta = text[activeEdge]
object ca = children[activeNode]
integer next = iff(ca=0?NULL:ca[ta])
if next == null then
if ca=0 then
children[activeNode] = repeat(NULL,maxChar)
end if
children[activeNode][ta] = newNode(pos, pLeafEnd)
if lastNewNode != NULL then
suffixLinks[lastNewNode] = activeNode
lastNewNode = NULL
end if
else
if walkDown(next) then
continue
end if
integer tp = text[pos]
if text[starts[next]+activeLength] == tp then
if lastNewNode != NULL and activeNode != root then
suffixLinks[lastNewNode] = activeNode
lastNewNode = NULL
end if
activeLength += 1
exit
end if
integer temp = starts[next] + activeLength - 1
pSplitEnd = allocate_word(temp,true)
integer splitnode = newNode(starts[next], pSplitEnd,true)
ta = text[activeEdge]
children[activeNode][ta] = splitnode
children[splitnode][tp] = newNode(pos, pLeafEnd)
starts[next] += activeLength
children[splitnode][text[starts[next]]] = next
if lastNewNode != NULL then
suffixLinks[lastNewNode] = splitnode
end if
lastNewNode = splitnode
end if
remainingSuffixCount -= 1
if activeNode == root and activeLength > 0 then
activeLength -= 1
activeEdge = pos - remainingSuffixCount + 1
elsif activeNode != root then
activeNode = suffixLinks[activeNode]
end if
end while
end procedure
 
procedure setSuffixIndexByDFS(integer n, integer labelHeight)
if n!=NULL then
-- Uncomment the 3 lines below to print suffix tree
--if starts[n]!=-1 then
-- printf(1,text[starts[n]..peekns(ends[n])])
--end if
if children[n]=0 then
suffixIndices[n] = size - labelHeight
-- Uncomment line below to print suffix index
--printf(1," [%d]\n", suffixIndices[n])
else
bool leaf = true
for i=1 to maxChar do
integer nci = children[n][i]
if nci != NULL then
-- Uncomment the 3 lines below to print suffix index
--if leaf and starts[n] != -1 then
-- printf(1," [%d]\n", suffixIndices[n])
--end if
leaf = false
setSuffixIndexByDFS(nci, labelHeight+edgeLength(nci))
end if
end for
if leaf then ?9/0 end if -- (sanity check)
end if
end if
end procedure
 
procedure buildSuffixTree()
size = length(text)
pRootEnd = allocate_word(-1,true)
root = newNode(-1, pRootEnd)
activeNode = root
for i=1 to size do
extendSuffixTree(i)
end for
integer labelHeight = 0
setSuffixIndexByDFS(root, labelHeight)
end procedure
 
procedure doTraversal(integer n, integer labelHeight, atom pMaxHeight, pSubstringStartIndex)
if n!=NULL then
integer nsi = suffixIndices[n], newHeight
if nsi == -1 then
for i=1 to maxChar do
integer nci = children[n][i]
if nci != NULL then
newHeight = labelHeight+edgeLength(nci)
doTraversal(nci, newHeight, pMaxHeight, pSubstringStartIndex)
end if
end for
elsif nsi > -1 then
newHeight = labelHeight-edgeLength(n)
if peekns(pMaxHeight)<newHeight then
poken(pMaxHeight,newHeight)
poken(pSubstringStartIndex,nsi)
end if
end if
end if
end procedure
 
procedure getLongestRepeatedSubstring(string s)
atom pMaxHeight = allocate_word(0,true),
pSubstringStartIndex = allocate_word(0,true)
doTraversal(root, 0, pMaxHeight, pSubstringStartIndex)
integer maxHeight = peekns(pMaxHeight),
start = peekns(pSubstringStartIndex)
-- Uncomment line below to print maxHeight and substringStartIndex
--printf(1,"maxHeight %d, substringStartIndex %d\n", {maxHeight, start})
string t = iff(maxHeight=0?"No repeated substring"
 :text[start+1..start+maxHeight])
printf(1,"  %s is: %s\n", {iff(s=""?text:s),t})
end procedure
 
constant tests = {"GEEKSFORGEEKS$",
"AAAAAAAAAA$",
"ABCDEFG$",
"ABABABA$",
"ATCGATCGA$",
"banana$",
"abcpqrabpqpq$",
"pqrpqpqabab$"}
printf(1,"Longest Repeated Substring in:\n")
for i=1 to length(tests) do
text = tests[i]
buildSuffixTree()
getLongestRepeatedSubstring("")
end for
printf(1,"\n")
 
-- Sadly gmp crashes when I try 100,000 dp, but I suspect it would be fine with more than my paltry 4GB ram
include mpfr.e
mpfr pi = mpfr_init(0,-10001) -- (set precision to 10,000 dp, plus the "3.")
mpfr_const_pi(pi)
string piStr = mpfr_sprintf("%.10000Rf", pi), s="Pi"
piStr = piStr[3..$] -- discard leading "3."
maxChar = '9'
for i=3 to 6 do
atom t0 = time()
integer n = power(10,i)
if i>=5 then
text = repeat('0',n)&"$"
for c=1 to n do
text[c] = rand(10)+'0'-1
end for
s = "a made up number"
else
text = piStr[1..n] & "$"
end if
buildSuffixTree()
getLongestRepeatedSubstring(sprintf("first %,d d.p. of %s", {n,s}))
printf(1," (this took %gs)\n\n", {time()-t0})
end for
Output:
Longest Repeated Substring in:
  GEEKSFORGEEKS$ is: GEEKS
  AAAAAAAAAA$ is: AAAAAAAAA
  ABCDEFG$ is: No repeated substring
  ABABABA$ is: ABABA
  ATCGATCGA$ is: ATCGA
  banana$ is: ana
  abcpqrabpqpq$ is: ab
  pqrpqpqabab$ is: ab

  first 1,000 d.p. of Pi is: 23846
  (this took 0s)

  first 10,000 d.p. of Pi is: 7111369
  (this took 0.047s)

  first 100,000 d.p. of a made up number is: 022105155
  (this took 0.422s)

  first 1,000,000 d.p. of a made up number is: 09014346795
  (this took 4.734s)