自從上一篇文章之後，還有很多事情要做，包括準備好JMP 16 ，僅在幾天前發布，並且沒有回复我上一篇詢問讀者遇到的列公式挑戰的帖子，我認為我沒有時間或材料來發布值得您花時間閱讀的內容。
我最近有更多的時間，是時候迎接新的挑戰了。 我們將再次列出列公式挑戰，並從新挑戰和新方法開始。 對於這篇文章，以及之後的幾篇文章，我將討論基於 Craige Hales 對社區帖子的回复的一般性問題. 克雷格 ( @Craige_Hales )，儘管已從 JMP Development 退休，但在社區中仍然非常活躍，回答用戶問題，維護他的未知博客，並為本博客做出貢獻。 問題是關於在不使用循環的情況下轉置（多維）列表的效率。 我想提出一個我經常想知道的更普遍的問題：“有沒有更好的迭代替代方案，如果有，它們是什麼？”在大多數情況下，更好意味著更快，但如果它不會對速度造成太大影響，我也會考慮代碼緊湊性。 我們還將研究不同的迭代方法，看看是否也存在差異。 讓我們從搜索文本開始。
字符串
生物信息學中的一項常見任務是搜索包含遺傳信息的字符串以找到特定的子字符串模式。 對於 DNA，字符串字符僅由四個字母組成：A、T、C 和 G。 文件Covid19.TXT （附於這篇文章）包含 Covid19 基因組的基因序列，長約 30,000 個字符（相比之下，最長的人類染色體 1 長約 2.9 億個字符）。 搜索長度為 2 到 10 個字符的特定子字符串的所有出現位置的速度有多快？
讓我們看看四種方法，其中三種是迭代的。 無論採用何種技術，請注意子字符串匹配可以重疊。 例如， BAB包含兩次巴巴，從位置 1 到 3，再從位置 3 到 5。 我們需要能夠捕捉這兩種情況的東西。 我們將從子串使用下面的代碼，其中字符串是要搜索的字符串序列是要搜索的子字符串。
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);
);
這是這些方法中最不令人印象深刻的方法，平均需要 1.2 秒，無論子字符串長度如何。 結果如下所示（以微秒為單位）。 測試為每個子串長度生成 10 個隨機子串，並重複搜索 20 次。 對給定長度的所有子字符串進行匯總。 各個子字符串的完整詳細信息在Covid 19 結果 Full.JMP附在這個帖子上。 用於生成測試的代碼位於 Test Covid19 Search.JSL 中。
| 子串長度 |
子串 |
指數 |
包含 |
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 |
我們在之前的挑戰中看到子串是緩慢的。 如果我們要將字符串“原子化”，將其分解為單個字符並查看字符範圍會怎樣？ 這樣做的動機是通過索引直接訪問字符串值，就像在 C 語言中一樣。 JSL 中需要原子化部分，因為您不能直接索引字符串。 我使用的代碼在這裡：
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);
);
雖然這種方法聽起來會更耗時，因為除了搜索之外還有分解搜索字符串的初始開銷，但從結果表中可以看出它的速度大約快 13-19 倍。 不像子串，但是，搜索時間隨著子字符串長度的增加而增加。
最後的迭代方法使用包含. 這個想法是找到給定的子字符串，在搜索的字符串中標記它的位置，然後重複這個過程，從該位置開始一個字符，直到找不到子字符串。 這種方法的吸引力在於，當搜索到的字符串中出現的子字符串很少時，搜索變得非常快。 下面的代碼感謝 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;
);
如您所見，雖然這種方法開始時比雙字符子字符串的索引方法稍慢，但隨著子字符串變長（即要查找的子字符串越來越少），它很快變得更快。
最後一種方法使用了我最喜歡的 JSL 工具之一，即模式匹配。 代碼不僅速度快，而且很短。
findSequence = Function({inString,inSeq},{i=1,locs={}},
Pat Match(inString,Pat Pos()>>r + inSeq + Pat Test(Insert Into(locs,r+1);0));
locs;
);
雖然兩個字符的子字符串稍微慢一些，但它似乎不受三個字符或更大字符的子字符串長度的影響。
顯然，迭代方法使用包含和模式匹配方法是搜索長字符串的優越方法。 但是哪種方法，何時？ 由於測試是使用不同的子字符串完成的，因此子字符串頻率（# 出現次數*1000/29,903）如何影響搜索時間的一些跡象。 下圖顯示了一個交叉點，其頻率約為每 1,000 個字符 2 個子串的頻率。 但是，由於所有測試都是在同一個 covid19 基因組字符串上完成的，因此尚不清楚搜索時間如何受搜索字符串長度的影響。
為了找到這些問題的答案，並探索因素之間可能的相互作用，我使用空間填充設計進行了一個實驗，改變了搜索字符串長度（10,000 - 500,000 個字符）、子字符串長度（2 - 21 個字符）和頻率（每 1,000 個字符 0.75 – 3 個子字符串）。 使用快速靈活設計生成了 125 次運行，其中小數字符串長度值四捨五入。 對於每次運行，都會生成一個隨機的字母字符串。 在該字符串中，由隨機生成的數字組成的子字符串被插入到隨機位置。 因為頻率值可能導致非整數插入次數，所以插入次數取為頻率上限乘以搜索字符串長度除以 1,000。 每次運行在相同的搜索字符串/子字符串組合上重複 10 次。 用於生成運行的代碼位於執行設計腳本.JSL結果在空間填充結果 1.JMP附在這個帖子上。
對結果的快速檢查只顯示了 13 個觀察結果，其中包含方法優於模式匹配。 其中，最大的差異是快了大約 4,000 微秒（而模式匹配比快了 0.72 秒）包含）。 將高斯過程模型擬合到數據表明，通過搜索字符串長度交互，搜索字符串長度具有最大的影響，並且頻率很高。 子串長度幾乎沒有影響。 線性模型最初適合數據，但殘差分析表明擬合可能存在問題。 我已經發布了用於此模型和後續模型的 Profiler 的交互式版本到 JMP 公共。
創建了第二個空間填充設計，以更好地表示以下情況包含優於模式匹配。 在這個實驗中，兩個弦長的範圍保持不變，頻率在 0 和 1 之間變化。 結果在空間填充結果 2.JMP . 結果的方向與以前相似，搜索的字符串長度更長，頻率更高，有利於模式匹配。
結合兩組結果並擬合高斯過程模型表明，對於 50,000 個字符的字符串，包含如果頻率低於 2，則更好。 在 100,000 個字符時，頻率下降到 1，在 200,000 個字符時，頻率下降到 0.5 以下。 由此得出的結論是，對於更稀疏的較短字符串，請使用包含並且隨著字符串長度和密度的增加，模式匹配將更加出色。
本月挑戰
使用 Covid19 基因組，找到出現兩次或多次的最長子串。
讓我們討論和分享，而不是等到下一篇文章才回答這個問題。 我相信這是一個非常困難的問題，我們將從公開討論中受益。
請找到另一種在某些情況下可能會更好也可能不會更好的方法（我沒有時間檢查）：
NamesDefaultToHere(1);
str = LoadTextFile("$DOWNLOADS/covid-19.txt");
txt = "AACA";
bStr = CharToBlob(str);
bTxt = CharToBlob(txt);
pos = {};
for(i=1, i<=(Length(str) - Length(txt)), i++,
if(BlobPeek(bStr, i, Length(txt)) == bTxt, InsertInto(pos, i+1));
);
Print(NItems(pos));包括重疊？ 這很容易。 不知道我會怎麼做，但不包括重疊。 此解決方案是提供數據的 5-10 秒答案。 對於具有更長子字符串的更長數據會很痛苦。
這個想法是計算在每個位置 1...文本長度的不同子字符串長度下每個子字符串有多少個出現。 然後問“是否有任何計數> 1？” 對子字符串長度使用二進制搜索來查找產生重複項的最大子字符串長度。 最後應該做更多的工作，看看是否有另一個更有趣的相同長度的子字符串。
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 );
包含方法可以更簡潔地編寫如下。 由於刪除了比較，這還具有運行速度更快的好處，最明顯的是當主字符串中不存在子字符串時。
findSequence = Function({txt, sub},{pos=0,ans={}},
while(pos = contains(txt, sub, pos), insertinto(ans, pos++));
ans;
);
