Con muchas actividades 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 solicitando 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 valga la pena leer.
Últimamente he tenido un poco más de tiempo y es hora de afrontar un nuevo desafío. Dejaremos 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 posteriores, 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, permanece muy activo aquí en la Comunidad, respondiendo preguntas de los usuarios, manteniendo su Blog inexplorado y contribuyendo a este blog. La pregunta era sobre la eficiencia de transponer una lista (multidimensional) sin utilizar bucles. Quiero plantear una pregunta más general sobre la que me he preguntado a menudo: "¿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 diferentes enfoques de iteración para ver si también existen diferencias allí. 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 subcadenas específicos. En el caso del ADN, los caracteres de la cadena se componen de sólo 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 humano 1, el más largo, tiene aproximadamente 290 millones de caracteres). ¿Con qué rapidez 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, BAB está contenido dos veces en BABAB , de las posiciones 1 a 3 y nuevamente de 3 a 5. Necesitamos algo que capture ambos sucesos. Empezaremos con substrato usando el siguiente código, donde en cadena 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, tomando 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 se encuentran en Resultados Covid 19 Completos.JMP adjunto a esta publicación. El código utilizado para generar la prueba se encuentra en Test Covid19 Search.JSL.
| Longitud de la subcadena |
substrato |
Í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 substrato es lento. ¿Qué pasaría si “atomizamos” la cadena, la dividimos en caracteres individuales y observamos los rangos de caracteres? La motivación para esto está en el orden de acceder directamente a los valores de cadena mediante la indexación, como es posible en un lenguaje como C. La parte de atomización es necesaria en JSL ya que no se puede indexar una cadena directamente. El código que utilicé 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 este enfoque parece que consumirí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 substrato , 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, si bien 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, cuando 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 sólo es rápido, sino que también es breve.
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 que utiliza 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 las pruebas se realizaron utilizando diferentes subcadenas, alguna indicación de cómo la frecuencia de las subcadenas (# apariciones*1000/29,903) puede afectar el tiempo de búsqueda. El siguiente gráfico 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 de 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 posibles interacciones entre factores, realicé un experimento utilizando un diseño de relleno de espacios, variando la longitud de la cadena buscada (10 000 – 500 000 caracteres), la longitud de la subcadena (2 – 21 caracteres) y la frecuencia ( 0,75 – 3 subcadenas por cada 1.000 caracteres). Se generaron 125 corridas utilizando un diseño Fast Flexible con valores de longitud de cadena fraccionaria redondeados. Para cada ejecución, se generó una cadena aleatoria de caracteres alfabéticos. Dentro de esa cadena, se insertó en ubicaciones aleatorias una subcadena formada por números generados aleatoriamente. Debido a que los valores de frecuencia podrían conducir a un número de inserciones no entero, el número de inserciones se tomó como el límite 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 llenado de espacios 1.JMP adjunto a esta publicación.
Una rápida inspección de los resultados muestra sólo 13 observaciones donde el Contiene El método supera la coincidencia de patrones. De estas, la mayor diferencia es aproximadamente 4.000 microsegundos más rápida (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 fuerte frecuencia mediante la 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 del Profiler para este y los modelos siguientes al público JMP.
Se creó un segundo diseño para llenar el espacio para obtener una mejor representación de los casos en los que Contiene coincidencia de patrones superada. Para este experimento, los rangos para ambas longitudes de cuerda se mantuvieron como antes y la frecuencia varió entre 0 y 1. Los resultados están en Resultados de llenado de espacios 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 favorecieron 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 y a 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 aumente la longitud y la densidad de la cuerda, la combinación de patrones será superior.
Reto del mes
Utilizando 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, analicemos y compartamos. Creo que éste es un problema muy difícil y nos beneficiaremos de un debate abierto.
¿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 proporcionados. Y sería doloroso para datos mucho más largos con subcadenas mucho más largas.
La idea es contar cuántas de cada subcadena ocurren para diferentes longitudes de subcadena, en cada posición 1...longitud del texto. Luego pregunte "¿hay algún recuento > 1?" Utilice una búsqueda binaria en la longitud de la subcadena para encontrar la longitud de subcadena más grande que produzca 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 );
