自從上一篇文章發布以來，發生了很多事情，包括為JMP 16 ，幾天前才發布，並且對我上一篇詢問讀者遇到的列公式挑戰的帖子沒有回复，我認為我沒有時間或材料來發布值得您花時間閱讀的內容。
最近我有了更多的時間，是時候接受新的挑戰了。 我們將再次列出列公式挑戰，並從新挑戰和新方法開始。 對於這篇文章以及之後的幾篇文章，我將討論基於 Craige Hales 對社區貼文的回應的一般性問題。 克雷格（ @Craige_Hales ），儘管已經從 JMP Development 退休，但仍然在社區中非常活躍，回答用戶問題，保持他的神秘海域博客，並為此博客做出貢獻。 問題是關於在不使用循環的情況下轉置（多維）列表的效率。 我想提出一個我經常想知道的更普遍的問題：“迭代是否有更好的替代方案？如果有，它們是什麼？”在大多數情況下，更好意味著更快，但如果不會過度影響速度，我也會考慮程式碼緊湊性。 我們也將研究不同的迭代方法，看看是否也有差異。 讓我們從搜尋文字開始。
弦樂
生物資訊學中的一個常見任務是搜尋包含遺傳資訊的字串以找到特定的子串模式。 對於 DNA，字串字元僅由四個字母組成：A、T、C 和 G。 文件新冠病毒.TXT （附於本文）包含 Covid19 基因組的遺傳序列，大約長 30,000 個字元（相較之下，人類 1 號染色體最長，大約長 2.9 億個字元）。 搜尋以多快的速度找到長度為 2 到 10 個字元的某些子字串的所有出現位置？
讓我們來看看四種方法，其中三種是迭代的。 無論採用何種技術，請注意子字串匹配可能會重疊。 例如，巴布包含兩次巴巴布，從位置 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 結果完整.JMP附在這篇文章中。 用於產生測試的程式碼位於 Test Covid19 Search.JSL 中。
| 子串長度 |
子字串 |
指數 |
包含 |
帕特匹配 |
| 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 個字元）和頻率（ 0.75 – 每 1,000 個字元 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...length 處，對於不同的子字串長度，每個子字串出現了多少個。 然後問“有沒有大於 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 );
Contains 方法可以寫得更簡潔，如下所示。 由於刪除了比較，這也具有運行速度更快的好處，尤其是當主字串中不存在子字串時。
findSequence = Function({txt, sub},{pos=0,ans={}},
while(pos = contains(txt, sub, pos), insertinto(ans, pos++));
ans;
);
