Myers Bit-Vector
## Introducción y fundamentos previos
El enfoque clásico para el problema de la búsqueda aproximada (buscar el patron _P_ en un texto _T_
con no mas de _k_ errores) se basa en la programación dinámica, que construye una matriz D de tamaño
*(M+1)*×*(N+1)* (donde _M_ es la longitud de _P_ y _N_ es la longitud de _T_), resultando en una
complejidad temporal de _O(M×N)_. Si bien es robusto, este costo es prohibitivo para textos
extensos.
El algoritmo de Eugene Myers, publicado en 1999, introdujo un paradigma diferente. Utiliza
paralelismo a nivel de bits para codificar el estado de la computación en los registros de la CPU.
En lugar de calcular explícitamente cada celda de la matriz _D_, el algoritmo de Myers representa
columnas enteras de la matriz (o sus diferencias) como vectores de bits y las actualiza en tiempo
constante (para _M≤W_) o _O(M/W)_ (donde _W_ es el tamaño de palabra del procesador), logrando una
aceleración masiva.
### Solución tradicional con programación dinámica
El enfoque tradicional para la búsqueda aproximada se basa en el algoritmo de Wagner-Fischer, que
calcula la distancia de Levenshtein (o distancia de edición). Resumidamente, utiliza la matriz de
programación dinámica _D_ de tamaño *(M+1)*×*(N+1)* mencionada previamente. El valor en cada celda
_D[i,j]_ representa el costo (distancia de edición) mínimo para encontrar una coincidencia entre los
primeros i caracteres del patrón (_P[1..i]_) y una subcadena del texto que termina en _T[j]_.
### Ejemplo de uso (Patrón "MAR" vs texto "PAR")
**Caso base:** La primera fila y primera columna se rellenan de la siguiente forma:
_D[0,j]=0_, _D[i,0]=i_
**Siguientes elementos:** Para el resto de elementos d[i,j] en la matriz, se utiliza el siguiente
criterio, fila por fila:
$$
d[i,j] = min\begin{cases}
D[i-1,j]+1 \\
D[i,j-1]+1 \\
D[i-1,j-1] + costo\_sust
\end{cases}
$$
donde:
$$
costo\_sust=\begin{cases}
0, \text{ si } P[i]=T[j] \\
1, \text{ si } P[i]\neq T[j]
\end{cases}
$$
### Ejemplo ilustrado
Si queremos por ejemplo, las coincidencias de la palabra "PAR" con no mas de un solo error en la
palabra "MAR", bastaría con revisar los elementos de la ultima fila que tengan un número menor o
igual a uno. Lo que significa que si la palabra PAR termina en ese elemento, entonces se cumplirá
con el criterio.
Como se ha mencionado previamente, el principal problema es la complejidad. Para encontrar las
coincidencias, debe calcular el valor de cada una de las M×N celdas de la matriz. El cálculo de cada
celda toma un tiempo constante _O(1)_, pero el número total de celdas es el producto de ambas
longitudes, resultando en una complejidad temporal de _O(M×N)_.
Esto es problemático debido a que el costo es multiplicativo. Si el texto es muy largo (ej. = 1.000
millones, como un genoma) y el patrón es moderadamente largo (ej. = 100), el número de operaciones
es prohibitivo (100 × 1.000 millones = 100 mil millones de operaciones). El rendimiento se degrada
linealmente con la longitud del patrón.
La solución que propone el algoritmo de Myers implica una reducción completa del calculo de la
matriz. El algoritmo no calcula los valores _D[i,j]_ directamente. En su lugar, explota la propiedad
de que los valores en celdas adyacentes de la matriz _D_ difieren en una cantidad pequeña.
### Preliminares
Antes de describir el algoritmo, se definen las **operaciones binarias** fundamentales utilizadas en
las ecuaciones:
| Símbolo | Nombre | Descripción | Ejemplo |
| ---------- | ----------------------------- | -------------------------------------------------------------------- | ------------------------- |
| `∨` | OR bit a bit | Devuelve 1 si al menos un bit de las dos posiciones comparadas es 1. | `1010 ∨ 1100 = 1110` |
| `∧` | AND bit a bit | Devuelve 1 solo si ambos bits son 1. | `1010 ∧ 1100 = 1000` |
| `¬` | NOT bit a bit | Invierte cada bit (0→1, 1→0). | `¬1010 = 0101` |
| `⊕` | XOR bit a bit | Devuelve 1 si los bits son diferentes. | `1010 ⊕ 1100 = 0110` |
| `+` | Suma binaria | Opera con acarreo entre bits (como una suma entera). | `0111 + 0001 = 1000` |
| `<<` | Desplazamiento a la izquierda | Mueve todos los bits una posición hacia los más significativos. | `0011 << 1 = 0110` |
| `bit_m(X)` | Bit más significativo | Bit en la posición del último carácter del patrón. | Si `X = 101`, `bit₃(X)=1` |
## Definición del algoritmo
El algoritmo propone su solución partiendo de la siguiente premisa; la diferencia entre entradas
adyacentes en la matriz de programación dinámica $D[i,j]$ ($\Delta$) es siempre un valor
$\in \{-1, 0, 1\}$ Así, se definen principalmente tres tipos de propiedades:
- Adyacencia horizontal $\Delta h_{i,j}=D_{i,j}-D_{i,j-1}\in \{-1,0,+1\}$
- Adyacencia vertical $\Delta v_{i,j}=D_{i,j}-D_{i-1,j}\in \{-1,0,+1\}$
- Diagonal $\Delta d_{i,j}=D_{i,j}-D_{i-1,j-1}\in \{0,+1\}$
A través de estas propiedades, se definen los siguientes vectores de bits o máscaras:
- Delta vertical positivo: $VP_{ij}\equiv (\Delta v_{i,j}=+1$)
- Delta vertical negativo: $VN_{ij}\equiv (\Delta v_{i,j}=-1$)
- Delta horizontal positivo: $HP_{ij}\equiv (\Delta h_{i,j}=+1$)
- Delta horizontal negativo: $HN_{ij}\equiv (\Delta h_{i,j}=-1$)
- Delta diagonal cero: $D0_{ij}\equiv (\Delta d_{i,j}=0$)
Estas cinco máscaras codifican completamente la evolución de la matriz $D[i,j]$ a través de
relaciones booleanas entre sí. A continuación se detallaran cuales son estas relaciones y su
significado. Las siguientes equivalencias definen cómo se relacionan las deltas verticales,
horizontales y diagonales entre posiciones adyacentes de la matriz, Estas fórmulas son el núcleo
lógico del algoritmo y permiten calcular toda una columna mediante operaciones de bits.
### Relación entre `HN`, `VP` y `D0`
$$HN_{i,j} \Leftrightarrow VP_{i,j-1} \land D0_{i,j}$$
Si el valor de una celda es menor que el de su izquierda (Δh = −1), entonces la celda anterior en la
misma fila debió haber tenido Δv = +1 y su diagonal inmediata no cambió (D0 = 1).
### Relación entre `VN`, `HP` y `D0`
$$VN_{i,j} \Leftrightarrow HP_{i-1,j} \land D0_{i,j}$$
Una disminución vertical (Δv = −1) ocurre cuando arriba hay una celda con incremento horizontal (Δh
= +1) y la diagonal actual es igual (D0 = 1).
### Relación entre `HP`, `VN`, `VP` y `D0`
$$HP_{i,j} \Leftrightarrow VN_{i,j-1} \lor \neg(VP_{i,j-1} \lor D0_{i,j})$$
Hay un aumento horizontal (Δh = +1) si:
- la celda anterior tenía Δv = −1, o
- la celda anterior **no tenía ni** Δv = +1 **ni** D0 = 1 (es decir, si el valor a la izquierda era
menor).
### Relación entre `VP`, `HN`, `HP` y `D0`
$$VP_{i,j} \Leftrightarrow HN_{i-1,j} \lor \neg(HP_{i-1,j} \lor D0_{i,j})$$
Se produce un aumento vertical (Δv = +1) si:
- la celda de arriba tenía Δh = −1, o
- la celda de arriba no tenía ni Δh = +1 ni D0 = 1.
### Relación entre `D0`, coincidencias y deltas
$$D0_{i,j} \Leftrightarrow (p_i = t_j) \lor VN_{i,j-1} \lor HN_{i-1,j}$$
Dos celdas consecutivas tienen el mismo valor (Δd = 0) si los caracteres coinciden ($p_i = t_j$), o
el valor se propagó desde la izquierda con Δh = −1, o desde arriba con Δv = −1.
### Codificación de caracteres
Adicionalmente, el algoritmo propone una codificación para cada carácter $c\in P$ del patrón $P$ de
largo $m$: Cada codificación tendrá $m$ valores binarios que se determinarán, de derecha a izquierda
de la siguiente forma para cada $i\in \{1,...m\}$:
- 1, si $c_i=p_i$
- 0, si $c_i \neq p_i$
Ejemplo con el patrón MAR:
- $M\rightarrow 001$
- $A\rightarrow 010$
- $R\rightarrow 100$
A esta codificación se le conoce como máscara de coincidencia ($Eq[t_j]$). Notar que si existe algún
carácter en el texto pero no en el patrón, entonces su máscara de coincidencia es un vector de ceros
de tamaño $m$.
### ¿Como utilizamos esta información para calcular equivalentemente la distancia de edición?
#### Resolución de dependencias
Como puede observarse en las relaciones anteriores, las variables `VP`, `VN`, `HP`, `HN` y `D0`
están interdependientes. Por ejemplo, `D0` depende de `HN`, que depende de `VP`, y esta última
depende nuevamente de `D0`. Para romper este ciclo y permitir el cálculo directo de toda una
columna, Myers propuso una formulación cerrada para `D0`.
Sea:
- `X` el vector que marca coincidencias o propagaciones negativas desde la izquierda:
$$
X = Eq[t_j] \lor VN_{i,j-1}
$$
- `Y` el vector que marca incrementos verticales desde la posición anterior:
$$
Y = VP_{i,j-1}
$$
El valor de `D0` puede calcularse **bit a bit** como:
$$
D0 = (((X \land Y) + Y) \oplus Y) \lor X
$$
Esta ecuación constituye el **núcleo algebraico del algoritmo de Myers**, y se basa en el hecho de
que la operación de suma binaria `(+)` propaga los bits de acarreo hacia arriba, simulando así la
propagación de coincidencias diagonales en la matriz de programación dinámica.
En términos intuitivos:
- `X` representa las coincidencias o “correcciones” que reducen la distancia.
- `Y` representa los posibles incrementos acumulados en la columna anterior.
- La suma `(X ∧ Y) + Y` permite extender coincidencias en secuencia a través de bits adyacentes, de
la misma manera que una cadena de ceros se propaga en una resta binaria.
### Cálculo de los vectores auxiliares
Una vez obtenido `D0`, el resto de los vectores puede actualizarse en una sola pasada mediante
operaciones lógicas a nivel de bit:
$$
\begin{aligned}
HP &= VN \lor \neg(VP \lor D0) \\
HN &= VP \land D0 \\
VP' &= HN \lor \neg(HP \lor X) \\
VN' &= HP \land X
\end{aligned}
$$
Cada línea de este conjunto de ecuaciones describe cómo los cambios horizontales y verticales se
propagan hacia la siguiente columna (`VP'`, `VN'`).
Así, cada paso del algoritmo avanza una columna completa en el texto `T`, manteniendo el estado del
patrón `P` comprimido en una palabra de bits.
### Actualización del puntaje de edición
Una vez que se actualizan `VP` y `VN`, el algoritmo puede determinar la diferencia en el valor de la
última fila (`D[M,j]`), correspondiente al costo total acumulado hasta el carácter `t_j`.
Para ello, se utiliza el bit más significativo (correspondiente a `p_M`) de los vectores `HP` y
`HN`:
$$
\text{si } HP[M] = 1 \Rightarrow \text{score} \leftarrow \text{score} + 1
$$
$$
\text{si } HN[M] = 1 \Rightarrow \text{score} \leftarrow \text{score} - 1
$$
Este puntaje (`score`) representa la **distancia de edición** entre el patrón `P` y el prefijo del
texto que termina en `t_j`. Así, los pasos de ejecución del algoritmo son los siguientes:
#### Paso 1: Inicialización
Entradas del algoritmo:
- Patron P[i,..m]
- Texto T[1,...n]
- tolerancia $k$
#### Paso 2: Búsqueda
### Ejemplo ilustrado (Búsqueda de patrón "MAR" en el texto "PAR")
## Implementación del algoritmo
## Aplicaciones del algoritmo
Una de las aplicaciones más relevantes del algoritmo de Myers es la identificación de variaciones
genéticas y mutaciones puntuales en secuencias de ADN. El genoma humano contiene aproximadamente 3
mil millones de pares de bases, y la búsqueda eficiente de patrones con errores es fundamental en
genómica.
**Ejemplo:** Consideremos la búsqueda de la secuencia del gen BRCA1 (implicado en cáncer de mama)
con hasta 2 mutaciones permitidas:
En la primera posición encontramos una mutación (C→C en lugar de T en la posición 4), y en la
segunda encontramos dos mutaciones (C→G en la posición 4 y A→A en la posición 5). Con el algoritmo
Myers, estas variaciones se identifican en tiempo lineal O(n) en lugar de O(m×n).
## Referencias
1. **Eugene W. Myers.**
_A Fast Bit-Vector Algorithm for Approximate String Matching._
_Journal of the ACM (JACM)_, Vol. 46, No. 3, 1999, pp. 395–415.
DOI: [10.1145/316542.316550](https://doi.org/10.1145/316542.316550)
2. **Eugene W. Myers.**
_Bit-Parallel string matching— Verification Presentation._
Documento técnico complementario
Presentación [https://www.mi.fu-berlin.de/wiki/pub/ABI/RnaSeqP4/myers-bitvector-verification.pdf](https://www.mi.fu-berlin.de/wiki/pub/ABI/RnaSeqP4/myers-bitvector-verification.pdf)
Para cada carácter c ∈ Σ:
Eq[c] ← vector binario de longitud m
Para i = 1 hasta m:
Eq[c][i] = 1 si P[i] = c
Eq[c][i] = 0 en caso contrario
VP ← (1...1)_m // todos los bits en 1
VN ← (0...0)_m // todos los bits en 0
score ← m // distancia inicial = longitud del patrón
Para cada carácter del texto $T[j]$ (con j desde 1 hasta n):
Eqj ← Eq[T[j]]
X ← Eqj ∨ VN
D0 ← (((X ∧ VP) + VP) ⊕ VP) ∨ X
HN ← VP ∧ D0
HP ← VN ∨ ¬(VP ∨ D0)
si bit_m(HP) = 1 → score ← score + 1
si bit_m(HN) = 1 → score ← score - 1
VP' ← (HN << 1) ∨ ¬((HP << 1) ∨ X)
VN' ← (HP << 1) ∧ X
VP ← VP'
VN ← VN'
si score ≤ k → agregar j a lista de coincidencias
#include <bits/stdc++.h>
using namespace std;
using u64 = uint64_t;
static constexpr int W = 64;
string bits_to_string(u64 x, int m) {
string s;
for (int i = m - 1; i >= 0; --i)
s.push_back((x >> i) & 1ULL ? '1' : '0');
return s;
}
inline pair<u64, int> shl_with_carry(u64 x, int carry_in) {
int carry_out = (int)((x >> (W - 1)) & 1ULL);
u64 res = (x << 1) | (u64)carry_in;
return {res, carry_out};
}
vector<u64> build_masks(int m) {
int nb = (m + W - 1) / W;
vector<u64> mask(nb, ~0ULL);
if (m % W) {
int lastbits = m % W;
mask[nb - 1] = (1ULL << lastbits) - 1ULL;
}
return mask;
}
vector<vector<u64>> build_Eq(const string &pattern) {
int m = (int)pattern.size();
int nb = (m + W - 1) / W;
vector<vector<u64>> Eq(256, vector<u64>(nb, 0ULL));
for (int i = 0; i < m; ++i) {
unsigned char c = (unsigned char)pattern[i];
int b = i / W;
int pos = i % W;
Eq[c][b] |= (1ULL << pos);
}
return Eq;
}
vector<int> myers_search_debug(const string &text, const string &pattern, int k) {
int n = (int)text.size();
int m = (int)pattern.size();
int nb = (m + W - 1) / W;
auto mask = build_masks(m);
auto Eq = build_Eq(pattern);
// === Paso 1: Inicialización ===
vector<u64> VP(nb), VN(nb);
for (int b = 0; b < nb; ++b) {
VP[b] = mask[b]; // VP ← (1...1)
VN[b] = 0ULL; // VN ← (0...0)
}
int score = m;
vector<int> occ;
vector<u64> HP(nb), HN(nb), D0(nb), Eqj(nb), X(nb);
vector<u64> VP_next(nb), VN_next(nb);
cout << "=== Algoritmo Myers Bit-Vector ===\n";
cout << "Patrón: \"" << pattern << "\" (" << m << " chars)\n";
cout << "Texto : \"" << text << "\" (" << n << " chars)\n\n";
// === Paso 2: Iteración principal ===
for (int j = 0; j < n; ++j) {
unsigned char tc = (unsigned char)text[j];
cout << "Iteración " << j + 1 << "\n";
cout << string(70, '-') << "\n";
for (int b = 0; b < nb; ++b) {
Eqj[b] = Eq[tc][b];
// === 1. Calcular X y D0 ===
X[b] = Eqj[b] | VN[b];
u64 tmp = X[b] & VP[b];
u64 add = tmp + VP[b];
D0[b] = (add ^ VP[b]) | X[b];
D0[b] &= mask[b];
// === 2. Calcular HN y HP ===
HN[b] = VP[b] & D0[b];
HP[b] = VN[b] | (~(VP[b] | D0[b]) & mask[b]);
}
// === 3. Actualizar el puntaje ===
int last_block = nb - 1;
int last_pos = (m - 1) % W;
int bit_hp = (int)((HP[last_block] >> last_pos) & 1ULL);
int bit_hn = (int)((HN[last_block] >> last_pos) & 1ULL);
if (bit_hp) score++;
else if (bit_hn) score--;
// === 4. Mostrar información ===
int bits_to_show = (m < 64) ? m : 64;
cout << "Codificación del carácter '" << text[j] << "' (Eq): "
<< bits_to_string(Eqj[0], bits_to_show) << "\n";
cout << "D0 = " << bits_to_string(D0[0], bits_to_show) << "\n";
cout << "HN = " << bits_to_string(HN[0], bits_to_show) << "\n";
cout << "HP = " << bits_to_string(HP[0], bits_to_show) << "\n";
cout << "VN = " << bits_to_string(VN[0], bits_to_show) << "\n";
cout << "VP = " << bits_to_string(VP[0], bits_to_show) << "\n";
cout << "Score actual = " << score << "\n";
// === 5. Calcular VN' y VP' ===
int carry_in_hp = 0, carry_in_hn = 0;
for (int b = 0; b < nb; ++b) {
auto [HP_shift, c1] = shl_with_carry(HP[b], carry_in_hp);
auto [HN_shift, c2] = shl_with_carry(HN[b], carry_in_hn);
carry_in_hp = c1;
carry_in_hn = c2;
VN_next[b] = (HP_shift & X[b]) & mask[b];
VP_next[b] = (HN_shift | (~(HP_shift | X[b]) & mask[b])) & mask[b];
}
// === 6. Actualizar los vectores ===
VN = VN_next;
VP = VP_next;
// === 7. Registrar coincidencias si aplica ===
if (score <= k) occ.push_back(j);
cout << "------------------------------------------\n\n";
}
// === Resultado final ===
cout << "\n=== Resultado ===\n";
cout << "Score final: " << score << "\n";
cout << "Coincidencias (índices finales): ";
for (int p : occ)
cout << p << " ";
cout << "\n";
return occ;
}