Seitdem kurz nach dem letzten Beitrag viel los ist, einschließlich der Vorbereitungen JMP 16 , das erst vor ein paar Tagen herauskam, und da auf meinen letzten Beitrag mit der Frage nach den Herausforderungen bei der Spaltenformel, auf die die Leser gestoßen waren, keine Antworten gab, glaubte ich nicht, dass ich die Zeit oder das Material hatte, etwas zu posten, das es wert wäre, gelesen zu werden.
Ich hatte in letzter Zeit etwas mehr Zeit und es ist Zeit für eine neue Herausforderung. Wir werden die Herausforderungen der Spaltenformel für ein anderes Mal verschieben und mit einer neuen Herausforderung und einem neuen Ansatz beginnen. In diesem Beitrag und wahrscheinlich auch in einigen darauffolgenden werde ich darüber sprechen eine allgemeine Frage, die auf einer Antwort von Craige Hales auf einen Community-Beitrag basiert . Craige ( @Craige_Hales ), obwohl er sich aus der JMP-Entwicklung zurückgezogen hat, bleibt hier in der Community sehr aktiv, beantwortet Benutzerfragen und pflegt seine eigenen Uncharted-Blog , und Beitrag zu diesem Blog. Die Frage betraf die Effizienz der Transponierung einer (mehrdimensionalen) Liste ohne Verwendung von Schleifen. Ich möchte eine allgemeinere Frage stellen, die ich mir oft gestellt habe: „Gibt es bessere Alternativen zur Iteration und wenn ja, welche?“ In den meisten Fällen bedeutet besser schneller, aber ich berücksichtige auch die Kompaktheit des Codes, wenn dadurch die Geschwindigkeit nicht zu stark beeinträchtigt wird. Wir werfen auch einen Blick auf verschiedene Iterationsansätze, um festzustellen, ob es auch dort Unterschiede gibt. Beginnen wir mit der Textsuche.
Saiten
Eine häufige Aufgabe in der Bioinformatik besteht darin, eine Zeichenfolge mit genetischen Informationen zu durchsuchen, um bestimmte Teilzeichenfolgenmuster zu finden. Bei DNA bestehen die Zeichenfolgen nur aus vier Buchstaben: A, T, C und G. Die Datei Covid19.TXT (diesem Beitrag beigefügt) enthält die genetische Sequenz für das Covid19-Genom, etwa 30.000 Zeichen lang (im Gegensatz dazu ist das menschliche Chromosom 1, das längste, etwa 290 Millionen Zeichen lang). Wie schnell kann gesucht werden, um den Ort aller Vorkommen bestimmter Teilzeichenfolgen mit einer Länge von 2 bis 10 Zeichen zu finden?
Schauen wir uns vier Ansätze an, drei davon iterativ. Beachten Sie unabhängig von der Technik, dass Teilzeichenfolgenübereinstimmungen überlappen können. Zum Beispiel, BAB ist zweimal darin enthalten BABAB , von Position 1 auf 3 und erneut von 3 auf 5. Wir brauchen etwas, das diese beiden Ereignisse erfasst. Wir beginnen mit Substr Verwenden Sie den folgenden Code, wo inString ist die zu durchsuchende Zeichenfolge und inSeq ist die Teilzeichenfolge, nach der gesucht werden soll.
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);
);
Dies war der am wenigsten beeindruckende Ansatz und dauerte unabhängig von der Länge des Teilstrings durchschnittlich 1,2 Sekunden. Die Ergebnisse werden unten angezeigt (in Mikrosekunden). Der Test generierte 10 zufällige Teilzeichenfolgen für jede Teilzeichenfolgenlänge und wiederholte die Suche 20 Mal. Die Werte wurden über alle Teilzeichenfolgen einer bestimmten Länge zusammengefasst. Vollständige Details zu einzelnen Teilzeichenfolgen finden Sie in Covid-19-Ergebnisse vollständig.JMP diesem Beitrag beigefügt. Der zum Generieren des Tests verwendete Code befindet sich in Test Covid19 Search.JSL.
| Teilstringlänge |
Substr |
Index |
Enthält |
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 |
Das haben wir bei früheren Herausforderungen gesehen Substr ist langsam. Was wäre, wenn wir die Zeichenfolge „atomisieren“, in einzelne Zeichen aufteilen und uns die Zeichenbereiche ansehen würden? Die Motivation hierfür liegt darin, über die Indizierung direkt auf Zeichenfolgenwerte zuzugreifen, wie dies in einer Sprache wie C möglich ist. Der Atomisierungsteil wird in JSL benötigt, da Sie einen String nicht direkt indizieren können. Der Code, den ich verwendet habe, ist hier:
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);
);
Dieser Ansatz scheint zwar zeitaufwändiger zu sein, da neben der Suche zunächst auch der Aufwand für das Aufteilen der gesuchten Zeichenfolge anfällt, er ist jedoch etwa 13-19-mal schneller, wie aus der Ergebnistabelle hervorgeht. nicht wie Substr Allerdings nimmt die Suchzeit mit der Länge der Teilzeichenfolge zu.
Der letzte iterative Ansatz verwendet Enthält . Die Idee besteht darin, die angegebene Teilzeichenfolge zu finden, ihre Position in der gesuchten Zeichenfolge zu markieren und den Vorgang beginnend mit einem Zeichen hinter der Position zu wiederholen, bis die Teilzeichenfolge nicht mehr gefunden werden kann. Der Vorteil dieser Methode besteht darin, dass die Suche sehr schnell erfolgt, wenn die Teilzeichenfolge in der gesuchten Zeichenfolge nur wenige Vorkommen aufweist. Der folgende Code ist Brady Brady zu verdanken ( @brady_brady 
findSequence = Function({inString,inSeq},{pos=0,ans={}},
while(1,
pos = contains(inString, inSeq, pos);
if(pos, insertinto(ans, pos++), break())
);
ans;
);
Wie Sie sehen, beginnt dieser Ansatz zwar etwas langsamer als die Indizierungsmethode für Teilzeichenfolgen mit zwei Zeichen, wird aber schnell schneller, je länger die Teilzeichenfolgen werden, d. h. je weniger Teilzeichenfolgen zu finden sind.
Der letzte Ansatz verwendet eines meiner Lieblings-JSL-Tools, den Mustervergleich. Der Code ist nicht nur schnell, sondern auch kurz.
findSequence = Function({inString,inSeq},{i=1,locs={}},
Pat Match(inString,Pat Pos()>>r + inSeq + Pat Test(Insert Into(locs,r+1);0));
locs;
);
Während es bei Teilzeichenfolgen mit zwei Zeichen etwas langsamer ist, scheint die Länge der Teilzeichenfolge ab drei Zeichen keinen Einfluss darauf zu haben.
Ganz klar, der iterative Ansatz verwendet Enthält und der Mustervergleichsansatz sind überlegene Methoden zum Durchsuchen langer Zeichenfolgen. Aber welche Methode und wann? Da die Tests mit verschiedenen Teilzeichenfolgen durchgeführt wurden, gibt es einige Hinweise darauf, wie sich die Häufigkeit der Teilzeichenfolgen (Anzahl Vorkommen*1000/29.903) auf die Suchzeit auswirken kann. Das folgende Diagramm zeigt einen Kreuzungspunkt bei etwa einer Häufigkeit von 2 Teilzeichenfolgen pro 1.000 Zeichen. Da jedoch alle Tests mit demselben Covid19-Genomstring durchgeführt wurden, ist unklar, wie sich die Suchzeit auf die Länge des gesuchten Strings auswirkt.
Um Antworten auf diese Fragen zu finden und mögliche Wechselwirkungen zwischen Faktoren zu untersuchen, habe ich ein Experiment mit einem raumfüllenden Design durchgeführt und dabei die Länge der gesuchten Zeichenfolge (10.000 – 500.000 Zeichen), die Länge der Teilzeichenfolge (2 – 21 Zeichen) und die Häufigkeit variiert ( 0,75 – 3 Teilzeichenfolgen pro 1.000 Zeichen). 125 Läufe wurden mit einem Fast-Flexible-Design mit gerundeten gebrochenen Stringlängenwerten generiert. Für jeden Lauf wurde eine zufällige Folge alphabetischer Zeichen generiert. Innerhalb dieser Zeichenfolge wurde an zufälligen Stellen eine Teilzeichenfolge aus zufällig generierten Zahlen eingefügt. Da Häufigkeitswerte zu einer nicht ganzzahligen Anzahl von Einfügungen führen können, wurde die Anzahl der Einfügungen als Obergrenze der Häufigkeit multipliziert mit der Länge der gesuchten Zeichenfolge dividiert durch 1.000 angenommen. Jeder Lauf wurde zehnmal mit derselben gesuchten Kombination aus Zeichenfolge und Teilzeichenfolge wiederholt. Der zum Generieren der Läufe verwendete Code befindet sich in Design Script.JSL ausführen und die Ergebnisse liegen vor Ergebnisse der Raumfüllung 1.JMP diesem Beitrag beigefügt.
Eine schnelle Überprüfung der Ergebnisse zeigt nur 13 Beobachtungen, bei denen die Enthält Die Methode übertrifft den Mustervergleich. Von diesen ist der größte Unterschied etwa 4.000 Mikrosekunden schneller (während der Mustervergleich sogar 0,72 Sekunden schneller war). Enthält ). Die Anpassung eines Gaußschen Prozessmodells an die Daten legt nahe, dass die gesuchte Zeichenfolgenlänge den größten Effekt mit einer starken Häufigkeit durch die Interaktion der gesuchten Zeichenfolgenlänge hatte. Die Länge der Teilzeichenfolge hatte kaum Auswirkungen. Zunächst wurde ein lineares Modell an die Daten angepasst, die Residuenanalyse deutete jedoch darauf hin, dass die Anpassung möglicherweise fraglich ist. Ich habe gepostet eine interaktive Version des Profilers für dieses und die Folgemodelle an JMP Public.
Ein zweites raumfüllendes Design wurde erstellt, um Fälle besser darzustellen, in denen Enthält übertraf den Mustervergleich. Für dieses Experiment blieben die Bereiche für beide Saitenlängen wie zuvor und die Frequenz wurde zwischen 0 und 1 variiert. Die Ergebnisse liegen vor Ergebnisse der Raumfüllung 2.JMP . Die Richtung der Ergebnisse war ähnlich wie zuvor, wobei längere gesuchte Zeichenfolgenlängen und höhere Frequenzen den Mustervergleich begünstigten.
Die Kombination beider Ergebnissätze und die Anpassung eines Gaußschen Prozessmodells legen nahe, dass für Zeichenfolgen mit 50.000 Zeichen Enthält ist besser, wenn die Frequenzen unter 2 liegen. Bei 100.000 Zeichen sinkt die Häufigkeit auf 1 und bei 200.000 Zeichen auf unter 0,5. Die Schlussfolgerung daraus ist, dass für spärlichere, kürzere Zeichenfolgen verwendet werden sollte Enthält Und mit zunehmender Saitenlänge und -dichte wird die Musteranpassung besser.
Herausforderung für den Monat
Finden Sie mithilfe des Covid19-Genoms den längsten Teilstring, der zwei- oder mehrmals vorkommt.
Anstatt dies bis zum nächsten Beitrag unbeantwortet zu lassen, lasst uns diskutieren und teilen. Ich glaube, das ist ein sehr schwieriges Problem, und wir werden von einer offenen Diskussion profitieren.
Einschließlich Überschneidungen? Das macht es einfach. Ich bin mir jedoch nicht sicher, wie ich es machen würde, ohne die Überlappungen. Bei dieser Lösung handelt es sich um eine 5-10-sekündige Antwort auf die bereitgestellten Daten. Und bei viel längeren Daten mit viel längeren Teilzeichenfolgen würde es schmerzhaft werden.
Die Idee besteht darin, zu zählen, wie viele von jedem Teilstring für unterschiedliche Teilstringlängen an jeder Position 1...Textlänge vorkommen. Fragen Sie dann: „Gibt es Zählungen > 1?“ Verwenden Sie eine binäre Suche nach der Teilzeichenfolgenlänge, um die größte Teilzeichenfolgenlänge zu finden, die Duplikate erzeugt. Am Ende sollte noch etwas mehr Arbeit geleistet werden, um zu sehen, ob es einen weiteren, interessanteren Teilstring derselben Länge gibt.
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 );
