Avec beaucoup de choses depuis juste après le dernier message, y compris se préparer pour JMP 16 , qui est sorti il y a seulement quelques jours, et sans réponse à mon dernier message demandant des défis de formule de colonne que les lecteurs avaient rencontrés, je ne pensais pas avoir le temps ou le matériel pour publier quelque chose qui vaille la peine d'être lu.
J'ai eu un peu plus de temps ces derniers temps, et il est temps pour un nouveau défi. Nous remettrons les défis de la formule de la colonne à une autre fois et commencerons par un nouveau défi et une nouvelle approche. Pour cet article, et probablement quelques-uns après, je parlerai de une question générale basée sur une réponse de Craige Hales à un message de la communauté . Craig ( @Craige_Hales ), bien qu'il ait pris sa retraite de JMP Development, reste très actif ici dans la Communauté, répondant aux questions des utilisateurs, maintenant son Blog inexploré , et contribuer à ce blog. La question portait sur l'efficacité de la transposition d'une liste (multidimensionnelle) sans utiliser de boucles. Je veux poser une question plus générale que je me suis souvent posée : « Existe-t-il de meilleures alternatives à l'itération et, si oui, quelles sont-elles ? Dans la plupart des cas, mieux signifie plus rapide, mais je considérerai également la compacité du code si cela ne compromet pas trop la vitesse. Nous examinerons également différentes approches de l'itération pour voir s'il existe également des différences. Commençons par rechercher du texte.
Cordes
Une tâche courante en bioinformatique consiste à rechercher une chaîne contenant des informations génétiques pour trouver des modèles de sous-chaîne spécifiques. Pour l'ADN, les caractères de la chaîne sont composés de seulement quatre lettres : A, T, C et G. Le fichier Covid19.TXT (joint à ce post) contient la séquence génétique du génome du Covid19, longue d'environ 30 000 caractères (en revanche, le chromosome 1 humain, le plus long, compte environ 290 millions de caractères). À quelle vitesse peut-on rechercher l'emplacement de toutes les occurrences de certaines sous-chaînes de 2 à 10 caractères ?
Examinons quatre approches, dont trois itératives. Quelle que soit la technique, sachez que les correspondances de sous-chaînes peuvent se chevaucher. Par exemple, BAB est contenu deux fois dans BABAB , des positions 1 à 3 et de nouveau de 3 à 5. Nous avons besoin de quelque chose qui capture ces deux événements. Nous allons commencer par Substr en utilisant le code ci-dessous, où inString est la chaîne à rechercher et inSeq est la sous-chaîne à rechercher.
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);
);
C'était la moins impressionnante des approches, prenant en moyenne 1,2 seconde, quelle que soit la longueur de la sous-chaîne. Les résultats sont présentés ci-dessous (en microsecondes). Le test a généré 10 sous-chaînes aléatoires pour chaque longueur de sous-chaîne et a répété la recherche 20 fois. Les valeurs ont été résumées sur toutes les sous-chaînes d'une longueur donnée. Les détails complets pour les sous-chaînes individuelles sont dans Résultats Covid 19 complets.JMP attaché à ce poste. Le code utilisé pour générer le test se trouve dans Test Covid19 Search.JSL.
| Longueur de la sous-chaîne |
Substr |
Indice |
Contient |
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 |
| sept |
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 |
| dix |
1 208 948 |
89 732 |
106 |
3 742 |
Nous avons vu dans les défis précédents que Substr est lent. Et si nous devions « atomiser » la chaîne, la diviser en caractères individuels et examiner les plages de caractères ? La motivation pour cela est de l'ordre d'accéder directement aux valeurs de chaîne via l'indexation comme vous le pouvez dans un langage comme C. La partie atomisation est nécessaire dans JSL car vous ne pouvez pas indexer directement une chaîne. Le code que j'ai utilisé est ici :
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);
);
Bien que cette approche semble prendre plus de temps car il y a la surcharge initiale de briser la chaîne recherchée en plus de la rechercher, elle est environ 13 à 19 fois plus rapide, comme le montre le tableau des résultats. Contrairement à Substr , cependant, le temps de recherche augmente en fonction de la longueur de la sous-chaîne.
L'approche itérative finale utilise Contient . L'idée est de trouver la sous-chaîne donnée, de marquer son emplacement dans la chaîne recherchée et de répéter le processus en commençant un caractère après l'emplacement jusqu'à ce que la sous-chaîne soit introuvable. L'attrait de cette méthode est que la recherche devient très rapide lorsqu'il y a peu d'occurrences de la sous-chaîne dans la chaîne recherchée. Le code ci-dessous est grâce à 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;
);
Comme vous pouvez le constater, bien que cette approche démarre légèrement plus lentement que la méthode d'indexation pour les sous-chaînes à deux caractères, elle s'accélère rapidement à mesure que les sous-chaînes s'allongent, c'est-à-dire qu'il y a moins de sous-chaînes à rechercher.
L'approche finale utilise ce qui est en train de devenir l'un de mes outils JSL préférés, le pattern matching. Non seulement le code est rapide, mais il est court.
findSequence = Function({inString,inSeq},{i=1,locs={}},
Pat Match(inString,Pat Pos()>>r + inSeq + Pat Test(Insert Into(locs,r+1);0));
locs;
);
Bien qu'il soit légèrement plus lent pour les sous-chaînes à deux caractères, il ne semble pas être affecté par la longueur de la sous-chaîne pour trois caractères et plus.
De toute évidence, l'approche itérative utilisant Contient et l'approche de correspondance de modèles sont des méthodes supérieures pour rechercher de longues chaînes. Mais quelle méthode, et quand ? Étant donné que les tests ont été effectués à l'aide de différentes sous-chaînes, certaines indications sur la manière dont la fréquence des sous-chaînes (nombre d'occurrences * 1000/29 903) peuvent affecter le temps de recherche. Le graphique ci-dessous montre un point de croisement quelque part autour d'une fréquence de 2 sous-chaînes pour 1 000 caractères. Étant donné que tous les tests ont été effectués sur la même chaîne de génome covid19, cependant, on ne sait pas comment le temps de recherche est affecté par la longueur de la chaîne recherchée.
Pour trouver les réponses à ces questions et explorer les interactions possibles entre les facteurs, j'ai mené une expérience en utilisant une conception de remplissage d'espace, en faisant varier la longueur de la chaîne recherchée (10 000 à 500 000 caractères), la longueur de la sous-chaîne (2 à 21 caractères) et la fréquence ( 0,75 - 3 sous-chaînes pour 1 000 caractères). 125 exécutions ont été générées à l'aide d'une conception Fast Flexible avec des valeurs de longueur de chaîne fractionnaires arrondies. Pour chaque exécution, une chaîne aléatoire de caractères alphabétiques a été générée. Dans cette chaîne, une sous-chaîne composée de nombres générés aléatoirement a été insérée à des emplacements aléatoires. Étant donné que les valeurs de fréquence pouvaient conduire à un nombre non entier d'insertions, le nombre d'insertions a été considéré comme le plafond de la fréquence multiplié par la longueur de la chaîne recherchée divisé par 1 000. Chaque exécution a été répétée 10 fois sur la même combinaison chaîne/sous-chaîne recherchée. Le code utilisé pour générer les exécutions est dans Exécuter le script de conception.JSL et les résultats sont là Résultats de remplissage d'espace 1.JMP attaché à ce poste.
Un examen rapide des résultats montre seulement 13 observations où le Contient La méthode surpasse la correspondance de modèle. Parmi ceux-ci, la plus grande différence est d'environ 4 000 microsecondes plus rapide (alors que la correspondance de modèle était jusqu'à 0,72 seconde plus rapide que Contient ). L'ajustement d'un modèle de processus gaussien aux données suggère que la longueur de la chaîne recherchée a eu le plus grand effet avec une fréquence élevée par l'interaction de la longueur de la chaîne recherchée. La longueur de la sous-chaîne avait peu d'effet. Un modèle linéaire a été initialement ajusté aux données, mais l'analyse des résidus a suggéré que l'ajustement pourrait être discutable. j'ai posté une version interactive du Profiler pour cela et les modèles de suivi à JMP Public.
Une deuxième conception de remplissage d'espace a été créée pour obtenir une meilleure représentation des cas où Contient a surpassé la correspondance de motifs. Pour cette expérience, les gammes pour les deux longueurs de cordes sont restées les mêmes et la fréquence a varié entre 0 et 1. Les résultats sont en Résultats de remplissage d'espace 2.JMP . La direction des résultats était similaire à celle d'avant avec des longueurs de chaîne recherchées plus longues et des fréquences plus élevées favorisant la correspondance de motifs.
La combinaison des deux ensembles de résultats et l'ajustement d'un modèle de processus gaussien suggèrent que pour des chaînes de 50 000 caractères, Contient est supérieur si les fréquences sont inférieures à 2. A 100 000 caractères, la fréquence tombe à 1 et à 200 000 caractères, en dessous de 0,5. La conclusion en est que pour les chaînes plus courtes et plus rares, utilisez Contient et à mesure que la longueur et la densité des cordes augmentent, la correspondance des motifs sera supérieure.
Défi du mois
À l'aide du génome de Covid19, trouvez la sous-chaîne la plus longue qui apparaît deux fois ou plus.
Plutôt que de laisser cela sans réponse jusqu'au prochain post, discutons et partageons. Je pense qu'il s'agit d'un problème très difficile, et nous bénéficierons d'une discussion ouverte.
Y compris les chevauchements ? Cela facilite les choses. Je ne sais pas comment je le ferais en excluant les chevauchements. Cette solution est une réponse de 5 à 10 secondes pour les données fournies. Et cela deviendrait pénible pour des données beaucoup plus longues avec des sous-chaînes beaucoup plus longues.
L'idée est de compter combien de chaque sous-chaîne se produisent pour différentes longueurs de sous-chaîne, à chaque position 1...longueur de texte. Demandez ensuite "y a-t-il des nombres > 1 ?" Utilisez une recherche binaire sur la longueur de sous-chaîne pour trouver la plus grande longueur de sous-chaîne qui produit des doublons. Un peu plus de travail devrait être fait à la fin pour voir s'il existe une autre sous-chaîne, plus intéressante, de la même longueur.
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 );
