Beaucoup de choses se sont passées depuis juste après le dernier post, y compris la préparation pour JMP16 , qui est sorti il y a seulement quelques jours, et sans réponse à mon dernier message demandant les 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 valait la peine d'être lu.
J'ai eu un peu plus de temps ces derniers temps, et il est temps de relever un nouveau défi. Nous reporterons les défis de la formule des colonnes à une autre fois et commencerons avec 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é . Craige ( @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 en contribuant à ce blog. La question portait sur l'efficacité de la transposition d'une liste (multidimensionnelle) sans utiliser de boucles. Je souhaite 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, lesquelles ? Dans la plupart des cas, mieux signifie plus rapide, mais je tiendrai également compte de la compacité du code si cela ne compromet pas trop la vitesse. Nous examinerons également différentes approches d'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înes spécifiques. Pour l'ADN, les chaînes de caractères sont composées de seulement quatre lettres : A, T, C et G. Le fichier Covid19.TXT (ci-joint à ce message) 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 effectuer une recherche pour trouver 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 encore de 3 à 5. Nous avons besoin de quelque chose qui capture ces deux événements. Nous commencerons par Sous-titre 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 l'approche la moins impressionnante, 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 des sous-chaînes individuelles se trouvent dans Résultats complets de Covid 19.JMP joint à ce post. Le code utilisé pour générer le test est dans Test Covid19 Search.JSL.
| Longueur de la sous-chaîne |
Sous-titre |
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 |
| 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 |
| dix |
1 208 948 |
89 732 |
106 |
3 742 |
Nous avons vu dans les défis précédents que Sous-titre 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 en est de l'ordre d'accéder directement aux valeurs de chaîne via l'indexation, comme vous pouvez le faire dans un langage comme C. La partie atomisation est nécessaire en JSL puisque 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 liée à la rupture de la chaîne recherchée en plus de sa recherche, elle est environ 13 à 19 fois plus rapide, comme le montre le tableau des résultats. Contrairement à Sous-titre , cependant, le temps de recherche augmente en fonction de la longueur de la sous-chaîne.
La dernière approche itérative 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 ne puisse plus être trouvée. L'intérêt 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 dû à 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 devient rapidement plus rapide à mesure que les sous-chaînes s'allongent, c'est-à-dire qu'il y a moins de sous-chaînes à trouver.
L'approche finale utilise ce qui est en train de devenir l'un de mes outils JSL préférés, la correspondance de modèles. 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 les 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 en utilisant différentes sous-chaînes, certaines indications sur la manière dont la fréquence des sous-chaînes (nombre d'occurrences*1 000/29 903) peut 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. Cependant, étant donné que tous les tests ont été effectués sur la même chaîne du génome covid19, 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 utilisant un modèle de remplissage d'espace, 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 rapide et 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 peuvent 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 en Exécuter Design Script.JSL et les résultats sont là Résultats de remplissage d'espace 1.JMP joint à ce post.
Une inspection rapide des résultats montre seulement 13 observations où le Contient La méthode surpasse la correspondance de modèles. Parmi ceux-ci, la plus grande différence est environ 4 000 microsecondes plus rapide (alors que la correspondance de modèles é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 chaîne recherchée a eu l'effet le plus important avec une forte fréquence par interaction entre la longueur de chaîne recherchée. La longueur des sous-chaînes 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 ce modèle et les modèles de suivi à JMP Public.
Un deuxième modèle remplissant l'espace a été créé pour obtenir une meilleure représentation des cas où Contient a surpassé la correspondance de modèles. Pour cette expérience, les plages pour les deux longueurs de cordes sont restées comme avant 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 précédente, avec des longueurs de chaîne recherchées plus longues et des fréquences plus élevées favorisant la correspondance de modèles.
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 chute à 1 et à 200 000 caractères, en dessous de 0,5. La conclusion est que pour les cordes plus courtes et plus clairsemées, utilisez Contient et à mesure que la longueur et la densité des chaînes augmentent, la correspondance des motifs sera supérieure.
Défi du mois
À l’aide du génome 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 que nous bénéficierons d’un débat ouvert.
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 sous-chaînes apparaissent pour différentes longueurs de sous-chaîne, à chaque position 1...longueur du texte. Demandez ensuite « y a-t-il des comptes > 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 );
