自上一篇文章发布以来，发生了很多事情，包括为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;
);
