hoy me gustaría volver un poco sobre el tema de los suffix array que ya ha sido tratado alguna vez en el blog.
En 1990, Manber y Myers publicaban Suffix arrays: a new method for on-line string searches, como motivo del primer ACM-SIAM symposium on Discrete algorithms. En éste fantástico e ilustrativo artículo, introducen la comparación de esta estructura con los suffix trees y proponen definiciones, notación y algoritmos para su construcción, búsqueda y optimización.
La mayor desventaja de suffix trees reside en el espacio requerido, y con los suffix arrays se consigue reducir considerablemente, siendo aún así competitivos frente a los suffix trees en términos de búsqueda y construcción.
Aquí veremos la primera parte, la creación de un suffix array ordenado (sección 3 del artículo original, sorprendentemente), con un algoritmo que requiere O(NlogN) como peor caso, si bien con optimizaciones se podría conseguir un tiempo esperado O(N); lo cual ya se consiguió de manera formal en 2003.
Quizás hoy día los suffix array van perdiendo protagonismo de primera plana, pero siguen siendo la base para muchas aplicaciones, por ejemplo implicados en métodos de obtención de la secuencia Barrows-Wheeler. Posiblemente en entradas posteriores hablemos de mejoras en el uso de suffix arrays y de otras estructuras relacionadas como FM-index y BWT.
Por ahora, aquí os dejo el código en Python (-V 2.6.6). Podéis cambiar la secuencia a indexar en la variable inicial A.
Y nos podemos preguntar ¿por qué usar éste código en lugar de la simple línea de ordenamiento MergeSort que aparecía en el código Perl de la entrada anterior sobre suffix arrays? Adjunto una tabla que lo explica en términos de tiempo requerido para la construcción del índice.
| Long texto | Myers | MergeSort |
| 20 | 0''019 | 0''004 |
| 100 | 0''023 | 0''005 |
| 1,000 | 0''055 | 0''013 |
| 11,000 | 0''717 | 0''401 |
| 110,000 | 22''200 | 41''075 |
| 220,000 | 1'36''004 | 2'42''189 |
| 660,000 | 9'31''297 | 25'54''598 |
Por supuesto estamos espectantes ante cualquier comentario, crítica o sugerencia.
¡Un saludo!
#!/usr/bin/python
# Simple algorithm for Suffix Array creation,
# based on Manber & Myers, 1990
# by Cantalapiedra, 2011
############### FUNCTIONS ##############
## Function to print a Suffix Array
def printSA(text, sa = [], sa_name = "Suffix Array"):
if len(sa) == 0:
print "No records in {0}".format(sa_name)
else:
print sa_name
for i in sa:
print "{0}\t{1}".format(i, text[i:])
print "\n"
return
## Function to print one or several lists
def printTemp(things = [], title = "##", post_title = "#"):
if len(things) == 0:
print "No things to print"
else:
print title
for a in things:
print a
print post_title
print "\n"
return
############### VARIABLES ##############
A = 'TTTTAGATCGATCGACTAGA$'
pos = [] # Contains the SA based on suffix numbers (SN)
prm = [] # Inverse of pos, each position is a SN containing its position in SA
bH = [] # Boolean array marking the buckets
b2H = [] # Temporary array for marking the 2H-buckets
count = [] # Counts the number of SN in a 2H-bucket
## All above variables have length = |A|
################ START ###################
# Some verbose
print "Your text is {0}\n".format(A[:len(A) - 1])
# Initial SA (POS)
pos = [(a,i) for i,a in enumerate(A)]
pos = sorted(pos)
pos = [a[1] for a in pos]
printSA(A, pos, "INITIAL SA")
# Initialize PRM
prm = [0]*len(A)
# Initial positions and buckets
prev_letter = ""
for i, pos_i in enumerate(pos):
prm[pos_i] = i
letter = A[pos_i]
if letter != prev_letter:
bH.append(1)
else:
bH.append(0)
prev_letter = letter
######### H-STAGES ##########
count_changed = []
h_stage = 1
while h_stage < len(A): # max log2(|A| + 1)
#print "stage {1} init {0}".format(time.time(), h_stage) IMPORT time module to use this
# FIRST: re-assign prm as of buckets starting points
# and create count and b2H arrays
prev_one = 0
for i, bh_i in enumerate(bH):
if bh_i == 1:
prev_one = i # Is a bucket starting position
if bh_i == 0:
prm[pos[i]] = prev_one # Points to its bucket starting pos
count.append(0)
b2H.append(0)
# Reset count[] to 0s
for a in range(len(count)):
count[a] = 0
# THEN: for every bucket...
# Here the entire SA is scanned, but in fact takes account
# of the current bucket being covered
prev_count = [a for a in count]
for i, pos_i in enumerate(pos):
if bH[i] == 1:
b2H[i] = 1
# We are moving the H-prefix (Ti) of suffix in i position
ti = pos_i - h_stage
if ti < 0: continue
### This is the heart of the algorithm
# If its the first change of this bucket during this H-bucket run...
if prev_count[prm[ti]] == count[prm[ti]]:
change = True
else:
change = False
count[prm[ti]] += 1
count_changed.append(prm[ti])
prm[ti] = prm[ti] + count[prm[ti]] - 1
# ... that means this is a new 2H-bucket
if change:
b2H[prm[ti]] = 1
#printTemp([prm, count, b2H], "After moving suff: {0}".format(ti), "")
# If next position is another H-bucket, or is the last bucket
if (i < len(bH) - 1 and bH[i + 1] == 1) or i == len(bH) - 1:
#printTemp([prm, count, b2H], "", "NEXT H-BUCKET {0}".format(i + 1))
for j in count_changed:
prev_count[j] = count[j]
count_changed = []
# If there are no buckets left
if 0 not in b2H:
break
#print "end stage {0}".format(time.time())
# If there are no buckets left
if 0 not in b2H:
break
#printTemp([prm, count, b2H], "", "NEXT ROUND {0}".format(h_stage*2))
# Update POS from PRM
for j, prm_j in enumerate(prm):
pos[prm_j] = j
# Update BH from B2H
bH = b2H
# Reset vars for next H-STAGE
count = []
b2H = []
h_stage *= 2
# Update POS from PRM
for j, prm_j in enumerate(prm):
pos[prm_j] = j
printSA(A, pos, "FINAL SA")
#
La salida del ejemplo es:
Your text is TTTTAGATCGATCGACTAGA
INITIAL SA
$
AGATCGATCGACTAGA$
ATCGATCGACTAGA$
ATCGACTAGA$
ACTAGA$
AGA$
A$
CGATCGACTAGA$
CGACTAGA$
CTAGA$
GATCGATCGACTAGA$
GATCGACTAGA$
GACTAGA$
GA$
TTTTAGATCGATCGACTAGA$
TTTAGATCGATCGACTAGA$
TTAGATCGATCGACTAGA$
TAGATCGATCGACTAGA$
TCGATCGACTAGA$
TCGACTAGA$
TAGA$
FINAL SA
$
A$
ACTAGA$
AGA$
AGATCGATCGACTAGA$
ATCGACTAGA$
ATCGATCGACTAGA$
CGACTAGA$
CGATCGACTAGA$
CTAGA$
GA$
GACTAGA$
GATCGACTAGA$
GATCGATCGACTAGA$
TAGA$
TAGATCGATCGACTAGA$
TCGACTAGA$
TCGATCGACTAGA$
TTAGATCGATCGACTAGA$
TTTAGATCGATCGACTAGA$
TTTTAGATCGATCGACTAGA$
