Arrays y árboles de sufijos

Qué cubre este capítulo

Los tres capítulos anteriores buscaban en un texto desde cero cada vez. Ese es el modelo correcto cuando recorres un texto una sola vez, pero es un desperdicio cuando vas a consultar el mismo texto una y otra vez: un genoma en el que buscarás miles de marcadores, un libro al que le harás grep repetidamente, un log que sondearás todo el día. Para eso preprocesas el texto una vez y lo conviertes en un índice buscable, y a partir de ahí cada consulta es rápida. Un suffix array es el índice clásico para esto, y su idea es elegante: ordenar todos los sufijos del texto. Como toda subcadena es el prefijo de algún sufijo, las ocurrencias de cualquier patrón forman un bloque contiguo de los sufijos ordenados, un bloque que encuentras con búsqueda binaria en O(m log n). Este capítulo construye el suffix array con prefix doubling, lo busca, calcula el array LCP que lo acompaña y habilita sus usos más ricos, y muestra la recompensa: un índice que cuesta 12 ms construir se paga solo después de apenas 50 consultas.

Un poco de historia

Los suffix arrays los introdujeron en 1990 Udi Manber y Gene Myers como una alternativa eficiente en memoria al árbol de sufijos, una estructura que Peter Weiner había inventado en 1973 (Donald Knuth llamó a la construcción en tiempo lineal de Weiner "el algoritmo del año"). Los árboles de sufijos son poderosos — responden muchas consultas sobre cadenas en tiempo óptimo — pero son glotones de memoria: en la práctica necesitan de 15 a 20 bytes por carácter, lo que los volvía impracticables para los textos enormes de la genómica. Manber y Myers demostraron que podías obtener casi todo el poder de un árbol de sufijos con un simple array ordenado de posiciones de inicio de sufijos más el array LCP, a una fracción de la memoria (4 bytes por carácter). Esa ganancia de espacio es la razón por la que los suffix arrays, y no los árboles de sufijos, se volvieron la columna vertebral de la búsqueda en ADN. Los algoritmos de construcción fueron mejorando durante dos décadas: Manber-Myers era O(n log n), y en 2003 tres grupos encontraron de forma independiente construcciones verdaderamente lineales en O(n) (el algoritmo DC3/skew y otros), y en 2009 llegó el muy práctico SA-IS. El prefix doubling de este capítulo es la idea de Manber-Myers: la más clara de entender, y suficientemente rápida para que se vea el punto.

La intuición

Toma el texto y lista todos sus sufijos: para banana son banana, anana, nana, ana, na, a. Ahora ordénalos alfabéticamente: a, ana, anana, banana, na, nana. El suffix array solo registra dónde empezaba cada sufijo ordenado en el texto original — aquí [5, 3, 1, 0, 4, 2] — así que son n enteros, no n cadenas. Aquí está la razón por la que ese orden es mágico para buscar: toda ocurrencia de un patrón en el texto es el inicio de algún sufijo que empieza con el patrón. Y en orden ordenado, todos los sufijos que empiezan con un prefijo dado quedan juntos en un solo bloque contiguo. Así que para encontrar todos los lugares donde aparece ana, haces búsqueda binaria sobre los sufijos ordenados buscando el bloque que empieza con ana: encuentras los sufijos ana (posición 3) y anana (posición 1), un bloque de dos, que te da ambas ocurrencias. Dos búsquedas binarias (una para el inicio del bloque y otra para el final) ubican todas las coincidencias en O(m log n).

Lo ingenioso es construir ese orden sin el costo O(n²) de comparar sufijos completos. El prefix doubling lo hace por rondas. Ronda uno: ordena los sufijos por su primer carácter. Ronda dos: ordena por sus primeros dos caracteres, pero no los vuelves a leer, reutilizas el ranking de la ronda uno, porque los primeros dos caracteres del sufijo i son (rank de i por 1 carácter, rank de i+1 por 1 carácter). Ronda tres: primeros cuatro caracteres, a partir de dos ranks traslapados de dos caracteres. Cada ronda duplica cuántos caracteres toma en cuenta el orden, así que después de log n rondas los sufijos quedan completamente ordenados, y cada ronda es un solo sort de pares de enteros. Eso es O(n log² n) — las log n rondas por el sort — y nunca compara cadenas crudas después de la primera ronda.

Complejidad: cómo escala

Construir el suffix array con prefix doubling es $O(n \log^2 n)$ (log n rondas, cada una un sort $O(n \log n)$); el algoritmo especializado SA-IS lo hace en $O(n)$ real. El almacenamiento es $O(n)$ enteros. Cada consulta es $O(m \log n)$: dos búsquedas binarias de $\log n$ pasos, cada una comparando hasta $m$ caracteres. El array LCP agrega $O(n)$ con el algoritmo de Kasai. Toda la propuesta de valor es la amortización — pagas el build una vez y cada consulta posterior sale barata — así que el face-off mide el tiempo acumulado para responder $K$ consultas, el suffix array (construir una vez, luego consultar) contra volver a recorrer el texto en cada consulta:

El suffix array arranca atrás: con una sola consulta ya gastó 12 ms construyendo el índice mientras que el rescaneo apenas recorre el texto una vez. Pero su costo por consulta es minúsculo, así que su línea casi no sube, mientras que la línea del rescaneo trepa linealmente: cada consulta vuelve a leer los 8000 caracteres completos del texto. Las dos se cruzan alrededor de las 50 consultas, y de ahí en adelante el índice se va adelantando cada vez más: para las 2000 consultas el rescaneo ha releído el texto dos mil veces mientras que el índice lo leyó una. Ese cruce es toda la lógica de indexar: el preprocesamiento es un costo fijo que amortizas, y vale la pena exactamente cuando vas a consultar suficientes veces para pagar el build. ¿Buscas en un texto una vez? Recórrelo. ¿Lo buscas cincuenta veces o un millón? Construye el índice.

En qué es bueno y en qué no

Los suffix arrays son la herramienta correcta cuando consultas un texto fijo repetidamente, o cuando necesitas consultas estructurales sobre todo el texto. La genómica es el caso insignia: un genoma de referencia se indexa una vez y se busca para millones de reads (el FM-index detrás de BWA y Bowtie es un suffix array comprimido). Los motores de búsqueda full-text, los detectores de plagio y las herramientas de compresión de datos (la Burrows-Wheeler Transform, el corazón de bzip2, es una rotación ordenada muy emparentada con el suffix array) se apoyan todos en ellos. Y la combinación de suffix array más LCP responde consultas de subcadena repetida más larga, conteo de subcadenas distintas y subcadena común más larga que los matchers de una sola pasada ni pueden tocar. Cuando el texto es fijo y las preguntas son muchas o estructurales, esta es la herramienta.

Donde es la elección equivocada es en una búsqueda de una sola vez: si vas a recorrer un texto nada más una vez, construir un índice es puro overhead; usa Boyer-Moore o str.find y sigue adelante. También le queda mal el texto que cambia con frecuencia: el suffix array se tiene que reconstruir cuando el texto cambia (no hay actualización incremental barata), así que va bien con corpus estáticos o que se actualizan poco, no con documentos que se editan en vivo. Y el build $O(n \log^2 n)$ en Python puro que hacemos aquí es muchísimo más lento que un SA-IS compilado; para textos de verdad a escala genómica usarías pydivsufsort o un índice estilo BWA, no Python hecho a mano. El suffix array es un índice: vale la pena construirlo solo cuando vas a consultar lo suficiente para amortizarlo.

Los datos, o las entradas

El face-off fija un texto aleatorio de 8000 caracteres y responde un número creciente de consultas de seis caracteres de dos maneras: construyendo el suffix array una vez y luego haciendo búsqueda binaria por consulta, contra volver a recorrer todo el texto en cada consulta, midiendo ambos en Python para que la comparación sea justa. La corrección se verifica en serio: sobre cientos de cadenas aleatorias el suffix array por prefix doubling debe ser exactamente igual a la referencia ingenua de ordenar los sufijos, search_all debe coincidir con el str.find de Python para muchos patrones, y el array LCP debe ser igual a un prefijo común más largo calculado por fuerza bruta entre sufijos ordenados adyacentes. La animación hace búsqueda binaria sobre los seis sufijos ordenados de banana buscando el bloque que empieza con ana.

Constrúyelo, una función a la vez

Construyendo el suffix array con prefix doubling — rankeando por 1, luego 2, luego 4, … caracteres:

def build_suffix_array(s):
    """The suffix array: the start indices of all suffixes of `s`, sorted by the suffixes they
    begin. Built by PREFIX DOUBLING — sort suffixes by their first character, then by their first
    2, then 4, 8, ... reusing the previous round's ranks so each doubling is a single sort of
    (rank, next-rank) pairs. O(n log^2 n). After log n rounds every suffix has a unique rank and
    the array is fully sorted."""
    n = len(s)
    sa = list(range(n))
    rank = [ord(c) for c in s]                   # round 0: rank by first character
    tmp = [0] * n
    k = 1
    while True:
        def key(i):
            return (rank[i], rank[i + k] if i + k < n else -1)
        sa.sort(key=key)                          # sort suffixes by (rank, rank k ahead)
        tmp[sa[0]] = 0
        for i in range(1, n):
            tmp[sa[i]] = tmp[sa[i - 1]] + (key(sa[i]) != key(sa[i - 1]))
        rank = tmp[:]
        if rank[sa[-1]] == n - 1:                 # all ranks distinct → fully sorted
            break
        k *= 2
    return sa

Buscándolo — dos búsquedas binarias acotan el bloque contiguo de coincidencias:

def search_all(s, sa, pattern):
    """Every start position of `pattern` in `s`, via two binary searches on the suffix array.
    The occurrences form a contiguous range [lo, hi) of the sorted suffixes — those whose prefix
    is `pattern`. Each comparison looks at up to m characters, and there are log n of them, so
    this is O(m log n). Returns the match positions (sorted)."""
    n, m = len(s), len(pattern)
    if m == 0:
        return sorted(range(n + 1))

    def lower(bound_pattern):
        lo, hi = 0, n
        while lo < hi:
            mid = (lo + hi) // 2
            if s[sa[mid]:sa[mid] + m] < bound_pattern:
                lo = mid + 1
            else:
                hi = mid
        return lo

    start = lower(pattern)                        # first suffix whose prefix ≥ pattern
    # first suffix whose m-prefix is strictly greater than pattern
    lo, hi = start, n
    while lo < hi:
        mid = (lo + hi) // 2
        if s[sa[mid]:sa[mid] + m] <= pattern:
            lo = mid + 1
        else:
            hi = mid
    end = lo
    return sorted(sa[start:end])

El array LCP con el algoritmo de Kasai — el complemento que desbloquea las consultas estructurales:

def kasai_lcp(s, sa):
    """The LCP array: lcp[i] is the length of the longest common prefix of the suffixes ranked i
    and i-1 in the suffix array (adjacent in sorted order). Kasai's algorithm computes it in O(n)
    by scanning suffixes in TEXT order and reusing the fact that dropping the first character can
    lower the LCP by at most one. LCP + suffix array together answer substring-count, longest-
    repeated-substring, and longest-common-substring queries."""
    n = len(s)
    rank = [0] * n
    for i, p in enumerate(sa):
        rank[p] = i
    lcp = [0] * n
    h = 0
    for i in range(n):
        if rank[i] > 0:
            j = sa[rank[i] - 1]
            while i + h < n and j + h < n and s[i + h] == s[j + h]:
                h += 1
            lcp[rank[i]] = h
            if h > 0:
                h -= 1                             # next suffix drops one leading char → LCP falls ≤ 1
        else:
            h = 0
    return lcp

Míralo funcionar

Aquí están los seis sufijos de banana, ordenados: a, ana, anana, banana, na, nana. Para encontrar todas las ocurrencias de ana, hacemos búsqueda binaria sobre esta lista ordenada buscando el bloque que empieza con ana. Naranja es el sufijo que se está comparando en el punto medio; la banda más oscura es el rango que sigue en juego. Míralo estrecharse: los primeros tres caracteres del sufijo del punto medio se comparan contra ana, y la búsqueda se mueve a la izquierda o a la derecha igual que una búsqueda binaria sobre números, porque los sufijos ordenados están ordenados, alfabéticamente. Converge en un bloque contiguo de dos sufijos, ana y anana, cuyas posiciones de inicio (3 y 1) son exactamente donde aparece ana en banana. Cada coincidencia, encontrada con dos búsquedas binarias sobre un índice construido una sola vez:

El código completo

La pestaña from-scratch es el suffix array — build por prefix doubling, consulta por búsqueda binaria y el LCP de Kasai; la pestaña de librería es el suffix array ingenuo (la referencia de corrección) y el contraste de rescanear en cada consulta, con la llamada a pydivsufsort que usarías en producción anotada. Cambia entre ellas.

"""Suffix arrays — index the TEXT once, then find any pattern fast. The last three chapters
searched a text from scratch every time. But if you'll query the same text many times — a book,
a genome, a log file — it pays to preprocess it into a searchable index once and answer every
later query in a flash. A suffix array is that index, and it's beautifully simple: the sorted
order of all the text's suffixes.

Every substring of a text is a prefix of one of its suffixes, so if you sort all n suffixes, the
occurrences of any pattern land in a CONTIGUOUS block of that sorted order — the suffixes that
start with the pattern. Find that block with binary search and you've found every occurrence, in
O(m log n), no matter how many times you search. The suffix array stores just the n starting
positions in sorted-suffix order (n integers, not the suffixes themselves), so it's compact. This
chapter builds it with prefix doubling, searches it by binary search, and computes the companion
LCP array that unlocks its richer uses.
"""


# region: build
def build_suffix_array(s):
    """The suffix array: the start indices of all suffixes of `s`, sorted by the suffixes they
    begin. Built by PREFIX DOUBLING — sort suffixes by their first character, then by their first
    2, then 4, 8, ... reusing the previous round's ranks so each doubling is a single sort of
    (rank, next-rank) pairs. O(n log^2 n). After log n rounds every suffix has a unique rank and
    the array is fully sorted."""
    n = len(s)
    sa = list(range(n))
    rank = [ord(c) for c in s]                   # round 0: rank by first character
    tmp = [0] * n
    k = 1
    while True:
        def key(i):
            return (rank[i], rank[i + k] if i + k < n else -1)
        sa.sort(key=key)                          # sort suffixes by (rank, rank k ahead)
        tmp[sa[0]] = 0
        for i in range(1, n):
            tmp[sa[i]] = tmp[sa[i - 1]] + (key(sa[i]) != key(sa[i - 1]))
        rank = tmp[:]
        if rank[sa[-1]] == n - 1:                 # all ranks distinct → fully sorted
            break
        k *= 2
    return sa
# endregion


# region: search
def search_all(s, sa, pattern):
    """Every start position of `pattern` in `s`, via two binary searches on the suffix array.
    The occurrences form a contiguous range [lo, hi) of the sorted suffixes — those whose prefix
    is `pattern`. Each comparison looks at up to m characters, and there are log n of them, so
    this is O(m log n). Returns the match positions (sorted)."""
    n, m = len(s), len(pattern)
    if m == 0:
        return sorted(range(n + 1))

    def lower(bound_pattern):
        lo, hi = 0, n
        while lo < hi:
            mid = (lo + hi) // 2
            if s[sa[mid]:sa[mid] + m] < bound_pattern:
                lo = mid + 1
            else:
                hi = mid
        return lo

    start = lower(pattern)                        # first suffix whose prefix ≥ pattern
    # first suffix whose m-prefix is strictly greater than pattern
    lo, hi = start, n
    while lo < hi:
        mid = (lo + hi) // 2
        if s[sa[mid]:sa[mid] + m] <= pattern:
            lo = mid + 1
        else:
            hi = mid
    end = lo
    return sorted(sa[start:end])
# endregion


# region: lcp
def kasai_lcp(s, sa):
    """The LCP array: lcp[i] is the length of the longest common prefix of the suffixes ranked i
    and i-1 in the suffix array (adjacent in sorted order). Kasai's algorithm computes it in O(n)
    by scanning suffixes in TEXT order and reusing the fact that dropping the first character can
    lower the LCP by at most one. LCP + suffix array together answer substring-count, longest-
    repeated-substring, and longest-common-substring queries."""
    n = len(s)
    rank = [0] * n
    for i, p in enumerate(sa):
        rank[p] = i
    lcp = [0] * n
    h = 0
    for i in range(n):
        if rank[i] > 0:
            j = sa[rank[i] - 1]
            while i + h < n and j + h < n and s[i + h] == s[j + h]:
                h += 1
            lcp[rank[i]] = h
            if h > 0:
                h -= 1                             # next suffix drops one leading char → LCP falls ≤ 1
        else:
            h = 0
    return lcp
# endregion

"""The library counterpart, the reference, and the contrast. In production you'd build a suffix
array with a compiled library — `pydivsufsort` (SA-IS, true O(n)) or a suffix-automaton library —
never in pure Python for large texts:

    from pydivsufsort import divsufsort
    sa = divsufsort(text.encode())

`naive_suffix_array` below is the O(n^2 log n) reference — sort the actual suffix strings — used to
check the prefix-doubling build. And the real contrast for the *search* is Python's `str.find`:
without an index, every query rescans the whole text (O(n) per query), where the suffix array pays
O(n log n) once and then answers each query in O(m log n). The face-off times "index once, query
many" against "rescan every query."
"""


# region: naive
def naive_suffix_array(s):
    """Sort the suffixes directly by materializing them — simple and obviously correct, but
    O(n^2 log n) time and O(n^2) space. The trusted reference for the prefix-doubling build."""
    return sorted(range(len(s)), key=lambda i: s[i:])
# endregion


# region: rescan
def find_all_builtin(text, pattern):
    """All occurrences via str.find — no index, so it rescans the text each call (O(n) per query).
    Both the correctness reference and the 'no index' contrast in the face-off."""
    if not pattern:
        return sorted(range(len(text) + 1))
    out, i = [], text.find(pattern)
    while i != -1:
        out.append(i)
        i = text.find(pattern, i + 1)
    return out
# endregion

Hecho a mano vs librería

La lección de este capítulo es el indexado como estrategia, y la gráfica del cruce en el face-off es su enunciado más claro: el preprocesamiento cambia un costo fijo por adelantado a cambio de consultas repetidas baratas, y vale la pena exactamente cuando el número de consultas rebasa el cruce. Ese es un patrón de decisión mucho más grande que la búsqueda en cadenas: es la misma lógica detrás de los índices de bases de datos, de ordenar una vez para hacer búsqueda binaria muchas, y de construir una tabla hash antes de un lote de lookups. El suffix array es la versión del mundo de las cadenas: no vuelvas a recorrer un texto que vas a interrogar mil veces; ordena sus sufijos una vez y haz búsqueda binaria para siempre. Los detalles del build también premian el estudio: la reutilización de ranks del prefix doubling es el mismo razonamiento de "reutiliza la ronda anterior" que la programación dinámica, y el LCP lineal de Kasai es el mismo argumento amortizado que KMP. En producción usarías un SA-IS compilado (pydivsufsort) o un FM-index; construir el array tú mismo es lo que convierte "indexa el texto" de una abstracción en la imagen concreta de sufijos ordenados y una búsqueda binaria convergiendo en un bloque.

Dónde te lo vas a encontrar de verdad

Los suffix arrays y sus descendientes comprimidos mueven las búsquedas de texto más grandes del mundo. En genómica, alineadores de reads como BWA y Bowtie indexan genomas de referencia con FM-indexes — suffix arrays comprimidos — para mapear miles de millones de reads de ADN, el caballo de batalla de la secuenciación moderna. La compresión de datos usa la muy emparentada Burrows-Wheeler Transform (rotaciones ordenadas) en bzip2 y en la compresión dentro de esas mismas herramientas genómicas. Los sistemas de búsqueda full-text y los motores de búsqueda de código usan suffix arrays o autómatas de sufijos para consultas de subcadena y regex sobre corpus fijos. Los detectores de plagio y de clones encuentran subcadenas comunes más largas con ellos. La bioinformática usa el array LCP para buscar repeticiones y para ensamblado de secuencias. Donde sea que un texto grande y estático se busque o analice muchas veces, hay un suffix array — normalmente comprimido — como índice por debajo.

Conclusiones

Un suffix array es el orden ordenado de los sufijos de un texto, guardado como n índices de inicio; como las ocurrencias de cualquier patrón forman un bloque contiguo de ese orden, la búsqueda binaria las encuentra todas en O(m log n) después de un build O(n log n). El prefix doubling lo construye rankeando por los primeros 1, 2, 4, … caracteres con reutilización de ranks, y el algoritmo de Kasai agrega el array LCP en O(n) que desbloquea consultas de repetición más larga, subcadenas distintas y subcadena común más larga: casi todo el poder de un árbol de sufijos con una cuarta parte de su memoria. Es un índice: vale la pena construirlo solo cuando vas a consultar un texto fijo lo suficiente para amortizar el costo, que el face-off ubicó alrededor de las 50 consultas.

El bloque de cadenas cierra con un regreso al matching, pero para muchos patrones a la vez. Aho-Corasick construye un solo autómata a partir de todo un diccionario de patrones y encuentra todas sus ocurrencias en una sola pasada lineal sobre el texto: es la función de fallo de KMP generalizada de un patrón a un trie de miles, y es el motor detrás de los antivirus, la detección de intrusiones en redes y cualquier herramienta que tenga que detectar un vocabulario fijo y grande en un stream.