Con molte cose in corso da subito dopo l'ultimo post, inclusa la preparazione JMP16 , uscito solo pochi giorni fa, e senza risposte al mio ultimo post in cui chiedevo sfide sulla formula delle colonne incontrate dai lettori, non pensavo di avere il tempo o il materiale per pubblicare qualcosa che valesse la pena leggere.
Ho avuto un po' più di tempo ultimamente ed è tempo per una nuova sfida. Presenteremo le sfide della formula delle colonne per un'altra volta e inizieremo con una nuova sfida e un nuovo approccio. Per questo post, e probabilmente per alcuni successivi, parlerò di una domanda generale basata su una risposta di Craige Hales a un post della community . Craige ( @Craige_Hales ), nonostante si sia ritirato dallo sviluppo di JMP, rimane molto attivo qui nella comunità, rispondendo alle domande degli utenti e mantenendo le sue Blog inesplorato e contribuendo a questo blog. La domanda riguardava l'efficienza della trasposizione di un elenco (multidimensionale) senza utilizzare cicli. Voglio porre una domanda più generale su cui mi sono spesso posto: "Esistono alternative migliori all'iterazione e, se sì, quali sono?" Nella maggior parte dei casi migliore significa più veloce, ma prenderò in considerazione anche la compattezza del codice se non compromette troppo la velocità. Daremo anche un'occhiata ai 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 modelli di sottostringhe specifici. 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, lungo circa 30.000 caratteri (al contrario, il cromosoma umano 1, il più lungo, è lungo circa 290 milioni di caratteri). Quanto velocemente è possibile effettuare la ricerca per trovare la posizione di tutte le occorrenze di determinate sottostringhe di lunghezza compresa tra 2 e 10 caratteri?
Esaminiamo quattro approcci, tre dei quali iterativi. Indipendentemente dalla tecnica, tieni presente che le corrispondenze di 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 Sottostr utilizzando il codice seguente, 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 l'approccio meno impressionante, impiegando in media 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. Sono disponibili i dettagli completi per le singole sottostringhe Risultati Covid 19 completi.JMP allegato a questo post. Il codice utilizzato per generare il test è in Test Covid19 Search.JSL.
| Lunghezza della sottostringa |
Sottostr |
Indice |
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 |
Lo abbiamo visto nelle sfide precedenti Sottostr è lento. E se dovessimo “atomizzare” la stringa, suddividerla in singoli caratteri e osservare 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);
);
Anche se questo approccio sembra richiedere più tempo poiché c'è il sovraccarico iniziale di suddividere la stringa cercata oltre a cercarla, è circa 13-19 volte più veloce come mostrato nella tabella dei risultati. A differenza di Sottostr , tuttavia, il tempo di ricerca aumenta in funzione della lunghezza della sottostringa.
L'approccio iterativo finale utilizza Contiene . L'idea è di trovare la sottostringa data, contrassegnare la sua posizione nella stringa cercata e ripetere il processo iniziando un carattere dopo la posizione finché non è possibile trovare la sottostringa. Il vantaggio di questo metodo è che la ricerca diventa molto veloce quando sono presenti poche occorrenze della sottostringa nella stringa cercata. Il codice seguente è 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, sebbene questo approccio inizi leggermente più lentamente rispetto al metodo di indicizzazione per sottostringhe di due caratteri, diventa rapidamente più veloce man mano che le sottostringhe diventano più lunghe, ovvero quando ci sono meno sottostringhe da trovare.
L'approccio finale utilizza quello che sta diventando uno dei miei strumenti JSL preferiti, il pattern match. 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 più.
Chiaramente, l'approccio iterativo utilizzando 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 sottostringhe diverse, alcune indicazioni su come la frequenza delle sottostringhe (# occorrenze*1000/29.903) potrebbero 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 cercata.
Per trovare le risposte a queste domande ed esplorare le possibili interazioni tra fattori, ho eseguito un esperimento utilizzando un design a riempimento di 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). Sono state generate 125 esecuzioni utilizzando un design flessibile e veloce con valori di lunghezza delle stringhe frazionari arrotondati. Per ogni esecuzione è 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 considerato come il limite massimo della frequenza moltiplicato per la lunghezza della stringa cercata divisa per 1.000. Ogni esecuzione è 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 arrivati Risultati riempimento spazio 1.JMP allegato a questo post.
Una rapida ispezione dei risultati mostra solo 13 osservazioni in cui il Contiene il metodo supera la corrispondenza dei modelli. Di questi, la differenza più grande è di circa 4.000 microsecondi più veloce (mentre la corrispondenza dei modelli 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 cercata ha avuto l'effetto maggiore con una forte frequenza mediante l'interazione della lunghezza della stringa cercata. La lunghezza della sottostringa ha avuto scarso effetto. Inizialmente un modello lineare era adatto ai dati, ma l’analisi dei residui ha suggerito che l’adattamento potrebbe essere discutibile. Ho pubblicato una versione interattiva del Profiler per questo e i modelli successivi a JMP pubblico.
È stato creato un secondo progetto che riempie lo spazio per ottenere una migliore rappresentazione dei casi in cui Contiene ha sovraperformato la corrispondenza dei modelli. Per questo esperimento, gli intervalli per entrambe le lunghezze delle corde sono rimasti gli stessi di prima e la frequenza è stata variata tra 0 e 1. I risultati sono arrivati Risultati riempimento spazio 2.JMP . La direzione dei risultati era simile a quella precedente, con lunghezze di corda cercate più lunghe e frequenze più alte che favorivano la corrispondenza dei modelli.
La combinazione di entrambi i set 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 e a 200.000 caratteri, al di sotto di 0,5. La conclusione da ciò è che per stringhe più corte e sparse, utilizzare Contiene e all'aumentare della lunghezza e della densità della corda, la corrispondenza del modello sarà superiore.
Sfida del mese
Utilizzando il genoma del Covid19, trova la sottostringa più lunga che appare due o più volte.
Piuttosto che lasciare questa domanda senza risposta fino al prossimo post, discutiamo e condividiamo. Credo che questo sia un problema molto difficile e trarremo beneficio da una discussione aperta.
Comprese le sovrapposizioni? Questo rende tutto più semplice. 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 è contare quante sottostringhe 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 maggiore che produce duplicati. Alla fine dovrebbe essere fatto un po' più di lavoro per vedere se c'è un'altra sottostringa più interessante 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 );
