Con un sacco di cose in corso da subito dopo l'ultimo post, incluso il prepararsi per JMP 16 , uscito solo pochi giorni fa, e senza risposte al mio ultimo post che chiedeva le sfide della formula della colonna che i lettori avevano incontrato, non pensavo di avere il tempo o il materiale per pubblicare qualcosa che valesse la pena leggere.
Negli ultimi tempi ho avuto un po' più di tempo ed è tempo di una nuova sfida. Esamineremo le sfide della formula della colonna per un'altra volta e inizieremo con una nuova sfida e un nuovo approccio. Per questo post, e probabilmente per alcuni dopo, parlerò una domanda generale basata su una risposta di Craige Hales a un post della community . Craige ( @Craige_Hales ), nonostante si sia ritirato da JMP Development, rimane molto attivo qui nella Community, rispondendo alle domande degli utenti, mantenendo le sue Blog inesplorato , e contribuendo a questo blog. La domanda riguardava l'efficienza della trasposizione di un elenco (multidimensionale) senza l'uso di loop. Voglio porre una domanda più generale che mi sono posto spesso: "Ci sono alternative migliori all'iterazione e, se sì, quali sono?" Nella maggior parte dei casi meglio significa più veloce, ma prenderò in considerazione anche la compattezza del codice se non compromette troppo la velocità. Daremo anche un'occhiata a diversi approcci all'iterazione per vedere se ci sono differenze anche lì. Iniziamo con la ricerca del testo.
stringhe
Un compito comune in bioinformatica è cercare una stringa contenente informazioni genetiche per trovare specifici modelli di sottostringa. Per il DNA, i caratteri della stringa sono composti solo da quattro lettere: A, T, C e G. Il file Covid19.TXT (allegato a questo post) contiene la sequenza genetica del genoma del Covid19, lunga circa 30.000 caratteri (al contrario, il cromosoma 1 umano, il più lungo, è lungo circa 290 milioni di caratteri). Quanto velocemente può essere cercato per trovare la posizione di tutte le occorrenze di determinate sottostringhe di lunghezza compresa tra 2 e 10 caratteri?
Diamo un'occhiata a quattro approcci, tre dei quali iterativi. Indipendentemente dalla tecnica, tieni presente che le corrispondenze delle sottostringhe possono sovrapporsi. Per esempio, BAB è contenuto due volte in BABABÙ , dalle posizioni da 1 a 3 e ancora da 3 a 5. Abbiamo bisogno di qualcosa che catturi entrambi questi eventi. Inizieremo con sost utilizzando il codice qui sotto, dove inString è la stringa da cercare e inSeq è la sottostringa da cercare.
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);
);
Questo è stato il meno impressionante degli approcci, impiegando una media di 1,2 secondi, indipendentemente dalla lunghezza della sottostringa. I risultati sono mostrati di seguito (in microsecondi). Il test ha generato 10 sottostringhe casuali per ciascuna lunghezza di sottostringa e ha ripetuto la ricerca 20 volte. I valori sono stati riepilogati su tutte le sottostringhe di una determinata lunghezza. I dettagli completi per le singole sottostringhe sono disponibili Risultati Covid 19 Full.JMP allegato a questo post. Il codice utilizzato per generare il test è in Test Covid19 Search.JSL.
| Lunghezza della sottostringa |
sost |
Indice |
Contiene |
Pat Match |
| 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 |
Lo abbiamo visto nelle sfide precedenti sost è lento. E se dovessimo "atomizzare" la stringa, suddividerla in singoli caratteri e guardare gli intervalli di caratteri? La motivazione di ciò è nell'ordine di accedere direttamente ai valori di stringa tramite l'indicizzazione come è possibile in un linguaggio come C. La parte di atomizzazione è necessaria in JSL poiché non è possibile indicizzare direttamente una stringa. Il codice che ho usato è qui:
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);
);
Sebbene questo approccio suoni come se richiederebbe più tempo poiché c'è il sovraccarico iniziale di spezzare la stringa cercata oltre a cercarla, è circa 13-19 volte più veloce come si vede nella tabella dei risultati. A differenza di sost , tuttavia, il tempo di ricerca aumenta in funzione della lunghezza della sottostringa.
L'ultimo approccio iterativo utilizza Contiene . L'idea è trovare la sottostringa data, contrassegnare la sua posizione nella stringa cercata e ripetere il processo iniziando un carattere oltre la posizione fino a quando la sottostringa non può essere trovata. Il vantaggio di questo metodo è che la ricerca diventa molto veloce quando ci sono poche occorrenze della sottostringa nella stringa cercata. Il codice qui sotto è grazie 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;
);
Come puoi vedere, mentre questo approccio inizia leggermente più lentamente del metodo di indicizzazione per le sottostringhe di due caratteri, diventa rapidamente più veloce man mano che le sottostringhe si allungano, cioè perché ci sono meno sottostringhe da trovare.
L'approccio finale utilizza quello che sta diventando uno dei miei strumenti JSL preferiti, il pattern matching. Non solo il codice è veloce, ma è breve.
findSequence = Function({inString,inSeq},{i=1,locs={}},
Pat Match(inString,Pat Pos()>>r + inSeq + Pat Test(Insert Into(locs,r+1);0));
locs;
);
Sebbene sia leggermente più lento per le sottostringhe di due caratteri, non sembra essere influenzato dalla lunghezza della sottostringa per tre caratteri e oltre.
Chiaramente, l'approccio iterativo che utilizza Contiene e l'approccio di corrispondenza dei modelli sono metodi superiori per la ricerca di stringhe lunghe. Ma quale metodo e quando? Poiché il test è stato eseguito utilizzando diverse sottostringhe, alcune indicazioni su come la frequenza delle sottostringhe (# occorrenze*1000/29.903) possa influenzare il tempo di ricerca. Il grafico seguente mostra un punto di crossover da qualche parte intorno a una frequenza di 2 sottostringhe per 1.000 caratteri. Poiché tutti i test sono stati eseguiti sulla stessa stringa del genoma covid19, tuttavia, non è chiaro in che modo il tempo di ricerca sia influenzato dalla lunghezza della stringa ricercata.
Per trovare le risposte a queste domande e per esplorare le possibili interazioni tra i fattori, ho eseguito un esperimento utilizzando un design che riempie lo spazio, variando la lunghezza della stringa cercata (10.000 - 500.000 caratteri), la lunghezza della sottostringa (2 - 21 caratteri) e la frequenza ( 0,75 – 3 sottostringhe per 1.000 caratteri). 125 corse sono state generate utilizzando un design Fast Flexible con valori di lunghezza della stringa frazionaria arrotondati. Per ogni corsa è stata generata una stringa casuale di caratteri alfabetici. All'interno di quella stringa, una sottostringa composta da numeri generati casualmente è stata inserita in posizioni casuali. Poiché i valori di frequenza potrebbero portare a un numero non intero di inserimenti, il numero di inserimenti è stato preso come limite massimo della frequenza moltiplicato per la lunghezza della stringa cercata diviso per 1.000. Ogni corsa è stata ripetuta 10 volte sulla stessa combinazione stringa/sottostringa cercata. Il codice utilizzato per generare le esecuzioni è in Esegui Design Script.JSL e i risultati sono in Risultati di riempimento dello spazio 1.JMP allegato a questo post.
Una rapida ispezione dei risultati mostra solo 13 osservazioni in cui il Contiene il metodo supera il pattern matching. Di questi, la differenza maggiore è di circa 4.000 microsecondi più veloce (mentre il pattern matching era fino a 0,72 secondi più veloce di Contiene ). L'adattamento di un modello di processo gaussiano ai dati suggerisce che la lunghezza della stringa ricercata ha avuto l'effetto maggiore con una forte frequenza dall'interazione della lunghezza della stringa ricercata. La lunghezza della sottostringa ha avuto scarso effetto. Inizialmente era stato adattato ai dati un modello lineare, ma l'analisi dei residui ha suggerito che l'adattamento potrebbe essere discutibile. ho postato una versione interattiva del Profiler per questo e per i modelli successivi a JMP Pubblico.
È stato creato un secondo progetto di riempimento dello spazio per ottenere una migliore rappresentazione dei casi in cui Contiene ha sovraperformato il pattern matching. Per questo esperimento, le gamme per entrambe le lunghezze delle corde sono rimaste come prima e la frequenza è stata variata tra 0 e 1. I risultati sono dentro Risultati di riempimento dello spazio 2.JMP . La direzione dei risultati era simile a prima, con lunghezze di corda più lunghe cercate e frequenze più alte che favorivano il pattern matching.
La combinazione di entrambi gli insiemi di risultati e l'adattamento di un modello di processo gaussiano suggerisce che per stringhe di 50.000 caratteri, Contiene è superiore se le frequenze sono inferiori a 2. A 100.000 caratteri, la frequenza scende a 1 ea 200.000 caratteri, inferiore a 0,5. La conclusione da ciò è che per stringhe più corte più sparse, utilizzare Contiene e all'aumentare della lunghezza e della densità delle corde, la corrispondenza del modello sarà superiore.
Sfida del mese
Usando il genoma di Covid19, trova la sottostringa più lunga che appare due o più volte.
Piuttosto che lasciare questo senza risposta fino al prossimo post, discutiamo e condividiamo. Credo che questo sia un problema molto difficile e trarremo vantaggio da una discussione aperta.
Comprese le sovrapposizioni? Questo lo rende facile. Tuttavia, non sono sicuro di come lo farei escludendo le sovrapposizioni. Questa soluzione è una risposta di 5-10 secondi per i dati forniti. E sarebbe doloroso per dati molto più lunghi con sottostringhe molto più lunghe.
L'idea è di contare quante di ogni sottostringa si verificano per diverse lunghezze di sottostringa, in ogni posizione 1...lunghezza del testo. Quindi chiedi "ci sono conteggi > 1?" Utilizzare una ricerca binaria sulla lunghezza della sottostringa per trovare la lunghezza della sottostringa più grande che produce duplicati. Un po' più di lavoro dovrebbe essere fatto alla fine per vedere se c'è un'altra, più interessante, sottostringa della stessa lunghezza.
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 );
