自从上一篇文章之后，还有很多事情要做，包括准备好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;
);
