Con muchas cosas desde justo después de la última publicación, incluida la preparación para JMP 16 , que salió hace solo unos días, y sin respuestas a mi última publicación en la que pedía los desafíos de fórmula de columna que los lectores habían encontrado, no pensé que tenía el tiempo o el material para publicar algo que valiera la pena leer.
He tenido un poco más de tiempo últimamente, y es hora de un nuevo desafío. Presentaremos los desafíos de la fórmula de la columna para otro momento y comenzaremos con un nuevo desafío y un nuevo enfoque. Para esta publicación, y probablemente algunas después de esa, hablaré sobre una pregunta general basada en una respuesta de Craige Hales a una publicación de la comunidad . craig ( @Craige_Hales ), a pesar de haberse retirado de JMP Development, sigue muy activo aquí en la Comunidad, respondiendo preguntas de los usuarios, manteniendo su Blog desconocido , y contribuyendo a este blog. La pregunta era sobre la eficiencia de transponer una lista (multidimensional) sin usar bucles. Quiero plantear una pregunta más general que a menudo me he preguntado: "¿Existen mejores alternativas a la iteración y, de ser así, cuáles son?" En la mayoría de los casos, mejor significa más rápido, pero también consideraré la compacidad del código si no compromete demasiado la velocidad. También veremos los diferentes enfoques de la iteración para ver si también hay diferencias. Comencemos con la búsqueda de texto.
Instrumentos de cuerda
Una tarea común en bioinformática es buscar una cadena que contenga información genética para encontrar patrones de subcadena específicos. Para el ADN, los caracteres de la cadena se componen de solo cuatro letras: A, T, C y G. El archivo covid19.txt (adjunto a esta publicación) contiene la secuencia genética del genoma de Covid19, de aproximadamente 30,000 caracteres (en contraste, el cromosoma 1 humano, el más largo, tiene aproximadamente 290 millones de caracteres). ¿Qué tan rápido se puede buscar para encontrar la ubicación de todas las apariciones de ciertas subcadenas de 2 a 10 caracteres de longitud?
Veamos cuatro enfoques, tres de ellos iterativos. Independientemente de la técnica, tenga en cuenta que las coincidencias de subcadenas pueden superponerse. Por ejemplo, bebé está contenido dos veces en BABAB , de la posición 1 a la 3 y nuevamente de la 3 a la 5. Necesitamos algo que capture ambas ocurrencias. Empezaremos con Sust. utilizando el siguiente código, donde enCadena es la cadena a buscar y inSeq es la subcadena a buscar.
findSequence = Function({inString,inSeq},{i,seqLen,checkTo,outVec},
seqLen = Length(inSeq);
checkTo = Length(covid19)-seqLen+1; //Only need to go until the substring includes the last character
outVec = J(checkTo,1,0);
For(i=1,i<=checkTo,i++,
outVec[i] = Substr(inString,i,seqLen) == inSeq;
);
Loc(outVec);
);
Este fue el menos impresionante de los enfoques, con un promedio de 1,2 segundos, independientemente de la longitud de la subcadena. Los resultados se muestran a continuación (en microsegundos). La prueba generó 10 subcadenas aleatorias para cada longitud de subcadena y repitió la búsqueda 20 veces. Los valores se resumieron en todas las subcadenas de una longitud determinada. Los detalles completos de las subcadenas individuales están en Resultados Covid 19 Completo.JMP adjunto a esta publicación. El código utilizado para generar la prueba está en Test Covid19 Search.JSL.
| Longitud de la subcadena |
Sust. |
Índice |
Contiene |
PatMatch |
| 2 |
1,227,342 |
60,312 |
79,640 |
4,046 |
| 3 |
1,245,909 |
62,650 |
22,143 |
3,103 |
| 4 |
1,215,641 |
66,045 |
4,938 |
2,881 |
| 5 |
1,218,983 |
72,337 |
1,293 |
2,880 |
| 6 |
1,220,846 |
76,705 |
357 |
2,956 |
| 7 |
1,249,843 |
77,127 |
234 |
3,010 |
| 8 |
1,212,982 |
84,413 |
132 |
3,078 |
| 9 |
1,227,892 |
87,479 |
113 |
3,331 |
| 10 |
1,208,948 |
89,732 |
106 |
3,742 |
Hemos visto en desafíos anteriores que Sust. es lento. ¿Qué pasaría si tuviéramos que "atomizar" la cadena, dividirla en caracteres individuales y observar los rangos de caracteres? La motivación para esto está en el orden de acceder directamente a los valores de cadena a través de la indexación que pueda en un lenguaje como C. La parte de atomización es necesaria en JSL ya que no puede indexar una cadena directamente. El código que usé está aquí:
findSequence = Function({inString,inSeq},{i,seqLen,checkTo,outVec,inList,seqList},
seqLen = Length(inSeq);
checkTo = Length(inString)-Length(inSeq)+1; //Only need to go until the substring includes the last character
outVec = J(checkTo,1,0);
inList = Items(inString,"");
seqList = Items(inSeq,"");
For(i=1,i<=checkTo,i++,
outVec[i] = inList[i::(i+seqLen-1)] == seqList;
);
Loc(outVec);
);
Si bien parece que este enfoque llevaría más tiempo, ya que existe la sobrecarga inicial de dividir la cadena buscada además de buscarla, es entre 13 y 19 veces más rápido, como se ve en la tabla de resultados. A diferencia de Sust. , sin embargo, el tiempo de búsqueda aumenta en función de la longitud de la subcadena.
El enfoque iterativo final utiliza Contiene . La idea es encontrar la subcadena dada, marcar su ubicación en la cadena buscada y repetir el proceso comenzando un carácter más allá de la ubicación hasta que no se pueda encontrar la subcadena. El atractivo de este método es que la búsqueda se vuelve muy rápida cuando hay pocas apariciones de la subcadena en la cadena buscada. El siguiente código es gracias a Brady Brady ( @brady_brady 
findSequence = Function({inString,inSeq},{pos=0,ans={}},
while(1,
pos = contains(inString, inSeq, pos);
if(pos, insertinto(ans, pos++), break())
);
ans;
);
Como puede ver, aunque este enfoque comienza un poco más lento que el método de indexación para subcadenas de dos caracteres, rápidamente se vuelve más rápido a medida que las subcadenas se hacen más largas, es decir, a medida que hay menos subcadenas para encontrar.
El enfoque final utiliza lo que se está convirtiendo en una de mis herramientas JSL favoritas, la coincidencia de patrones. El código no solo es rápido, sino que es corto.
findSequence = Function({inString,inSeq},{i=1,locs={}},
Pat Match(inString,Pat Pos()>>r + inSeq + Pat Test(Insert Into(locs,r+1);0));
locs;
);
Si bien es un poco más lento para las subcadenas de dos caracteres, no parece verse afectado por la longitud de la subcadena para tres caracteres o más.
Claramente, el enfoque iterativo usando Contiene y el enfoque de coincidencia de patrones son métodos superiores para buscar cadenas largas. Pero, ¿qué método y cuándo? Dado que la prueba se realizó con diferentes subcadenas, alguna indicación de cómo la frecuencia de la subcadena (# ocurrencias*1000/29,903) puede afectar el tiempo de búsqueda. La siguiente gráfica muestra un punto de cruce en algún lugar alrededor de una frecuencia de 2 subcadenas por cada 1000 caracteres. Sin embargo, debido a que todas las pruebas se realizaron en la misma cadena del genoma covid19, no está claro cómo el tiempo de búsqueda se ve afectado por la longitud de la cadena buscada.
Para encontrar las respuestas a estas preguntas y explorar las posibles interacciones entre los factores, realicé un experimento usando un diseño de relleno de espacio, variando la longitud de la cadena buscada (10 000 a 500 000 caracteres), la longitud de la subcadena (2 a 21 caracteres) y la frecuencia ( 0,75 – 3 subcadenas por 1000 caracteres). Se generaron 125 ejecuciones utilizando un diseño Fast Flexible con valores de longitud de cadena fraccionarios redondeados. Para cada ejecución, se generó una cadena aleatoria de caracteres alfabéticos. Dentro de esa cadena, se insertó una subcadena formada por números generados aleatoriamente en ubicaciones aleatorias. Debido a que los valores de frecuencia podrían conducir a un número no entero de inserciones, el número de inserciones se tomó como el techo de la frecuencia multiplicada por la longitud de la cadena buscada dividida por 1000. Cada ejecución se repitió 10 veces en la misma combinación de cadena/subcadena buscada. El código utilizado para generar las ejecuciones está en Ejecutar script de diseño.JSL y los resultados están en Resultados de relleno de espacio 1.JMP adjunto a esta publicación.
Una inspección rápida de los resultados muestra solo 13 observaciones donde el Contiene método supera a la coincidencia de patrones. De estos, la mayor diferencia es de unos 4.000 microsegundos más rápido (mientras que la coincidencia de patrones fue hasta 0,72 segundos más rápida que Contiene ). Ajustar un modelo de proceso gaussiano a los datos sugiere que la longitud de la cadena buscada tuvo el mayor efecto con una frecuencia fuerte por interacción de la longitud de la cadena buscada. La longitud de la subcadena tuvo poco efecto. Inicialmente se ajustó un modelo lineal a los datos, pero el análisis de residuos sugirió que el ajuste podría ser cuestionable. he publicado una versión interactiva de Profiler para este y los modelos de seguimiento a JMP Público.
Se creó un segundo diseño que llena el espacio para obtener una mejor representación de los casos en los que Contiene superó la coincidencia de patrones. Para este experimento, los rangos para ambas longitudes de cuerda se mantuvieron como antes y la frecuencia se varió entre 0 y 1. Los resultados están en Resultados de relleno de espacio 2.JMP . La dirección de los resultados fue similar a la anterior, con longitudes de cadena buscadas más largas y frecuencias más altas que favorecían la coincidencia de patrones.
La combinación de ambos conjuntos de resultados y el ajuste de un modelo de proceso gaussiano sugiere que para cadenas de 50 000 caracteres, Contiene es superior si las frecuencias están por debajo de 2. A 100.000 caracteres, la frecuencia baja a 1 ya 200.000 caracteres, por debajo de 0,5. La conclusión de esto es que para cadenas más cortas y dispersas, use Contiene y a medida que aumenta la longitud y la densidad de la cuerda, la coincidencia de patrones será superior.
Reto del mes
Usando el genoma de Covid19, encuentre la subcadena más larga que aparece dos o más veces.
En lugar de dejar esto sin respuesta hasta la próxima publicación, hablemos y compartamos. Creo que este es un problema muy difícil, y nos beneficiaremos de una discusión abierta.
¿Incluyendo superposiciones? Eso lo hace fácil. Sin embargo, no estoy seguro de cómo lo haría excluyendo las superposiciones. Esta solución es una respuesta de 5 a 10 segundos para los datos suministrados. Y sería doloroso para datos mucho más largos con subcadenas mucho más largas.
La idea es contar cuántos de cada subcadena ocurren para diferentes longitudes de subcadena, en cada posición 1...longitud de texto. Luego pregunte "¿hay cuentas > 1?" Utilice una búsqueda binaria en la longitud de la subcadena para encontrar la mayor longitud de la subcadena que produce duplicados. Se debe trabajar un poco más al final para ver si hay otra subcadena más interesante de la misma longitud.
txt = Load Text File( "z:/covid-19.txt" );
// finding the longest duplicated substring, including overlaps.
// if there is one, then there will be shorter ones of every possible
// shorter length. So, as the possible length of the substring grows,
// there is a boundary at some length where all shorter lengths have
// duplicates and no longer lengths have duplicates. find the boundary.
// do a binary search for the longest string that occurs 2 or more times
// testlength is the length during a pass
// see https://en.wikipedia.org/wiki/Binary_search_algorithm#Procedure
// for a *real* binary search. this code is potentially wrong and very
// poorly tested!
// Changing the maxlength guess below from 1...1000 will produce different times.
// 20 is a good starting point, bigger than 200 or so is very slow because
// length(txt) * maxlength strings are being manipulated, a lot.
start = Tick Seconds();
minlength = 1;
maxlength = 27; // initial guess, will (slowly) auto-adjust bigger as needed
firstpass = 1; // needed to detect guess too small
minkey = .;
minval = .;
while( 1,
testlength = if( firstpass == 1,
Min( Length( txt ), maxlength )// this binary search starts at the big end, then splits on subsequent passes. This lets code
, // else
Floor( (minlength + maxlength) / 2 )// below test that the big end is big enough on the first pass.
);
Write( "testing strings of length ", testlength, "\!n" );
Wait( 0 );
aa = [=> 0];
patmatch( txt,
Pat Len( testlength ) >> sequence + //
pattest( // build the associative array of key=sequence, val=count
aa[sequence]++;
// usually I use 1; on the next line. that always succeeds.
// here, I want pattest to make the matcher back up and retry...
0; // fail! try again at next position
),
NULL, // unused positional replacement
FULLSCAN // needed to bypass a heuristic and get the last little bit
);
keys = aa << getkeys;
vals = aa << getvalues;
idx = Loc Max( Matrix( vals ) );
bigN = vals[idx];
//
// the binary search is for the biggest testlength that has bigN > 1.
// bigN is the maximum repeats of a string of testlength characters.
//
if( bigN == 1,
maxlength = testlength;// maxlength promises to be 1-too-big at the end
, //else
minlength = testlength;
Write( " found some, minlength=", minlength, " maxlength=", maxlength, "\!n" );
minkey = keys[idx];
minval = vals[idx];
if( firstpass == 1, // looking for a big-enough high end, still
maxlength = Min( 2 * maxlength, Length( txt ) );
Write( "Increasing max, minlength=", minlength, " maxlength=", maxlength, "\!n" );
Continue();
);
);
firstpass = 0;
if( minlength >= maxlength - 1,
if( Is Missing( minkey ),
Continue();// needs one more pass
, // else
Break();// got it
)
);
);
Write( "\!n", minkey, " occurs ", minval, " times at " );
patmatch(
txt,
Pat Pos() >> location + //
minkey + //
pattest(
Write( location + 1, " " );
1;
) + //
Pat Fail() // this will also retry the match, 1 character further along
);
Write( "\!n" );
stop = Tick Seconds();
Show( stop - start );
