前回の投稿直後から、準備を含めて多くのことが行われています。 JMP16 、これはほんの数日前に公開されたもので、読者が遭遇した列数式の課題を尋ねた前回の投稿には返答がなかったため、時間をかけて読む価値のあるものを投稿する時間も資料も自分にはないと思いました。
最近は少し時間ができたので、新たな挑戦をする時期が来ました。 列式の課題については別の機会に取り上げ、新しい課題と新しいアプローチから始めます。 この投稿と、おそらくその後数回にわたって、私は次のことについて話します。コミュニティ投稿に対するクレイジ・ヘイルズの回答に基づく一般的な質問。 クレイジ ( @Craige_Hales ) は、JMP 開発を退職したにもかかわらず、ここコミュニティで非常に活発に活動を続けており、ユーザーの質問に答え、彼の意見を維持しています。未知のブログ、そしてこのブログに貢献しています。 質問は、ループを使用せずに (多次元) リストを転置する効率に関するものでした。 私がよく疑問に思っている、より一般的な質問を提起したいと思います。「反復に代わるより良い方法はありますか? ある場合、それは何ですか?」ほとんどの場合、より良いということはより高速であることを意味しますが、速度をあまり犠牲にしない場合は、コードのコンパクト性も考慮します。 また、反復に対するさまざまなアプローチも見て、違いがあるかどうかを確認します。 まずはテキストの検索から始めましょう。
文字列
バイオインフォマティクスにおける一般的なタスクは、遺伝情報を含む文字列を検索して、特定の部分文字列パターンを見つけることです。 DNA の場合、文字列文字は A、T、C、G の 4 つの文字のみで構成されます。 ファイルCovid19.TXT (この投稿に添付) には、約 30,000 文字の長さの Covid19 ゲノムの遺伝子配列が含まれています (対照的に、最も長いヒト染色体 1 は約 2 億 9,000 万文字の長さです)。 長さ 2 ～ 10 文字の特定の部分文字列のすべての出現場所を見つけるために、どのくらいの速さで検索できるでしょうか?
4 つのアプローチを見てみましょう。そのうち 3 つは反復的です。 手法に関係なく、部分文字列の一致は重複する可能性があることに注意してください。 例えば、バブに2回含まれていますババブ、位置 1 から 3 まで、そして再び 3 から 5 まで。 これら両方の出来事を捉えるものが必要です。 まずはサブ文字列以下のコードを使用します。ここで、文字列内は検索する文字列であり、 inSeq検索する部分文字列です。
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 にあります。
| 部分文字列の長さ |
サブ文字列 |
索引 |
含まれています |
パットマッチ |
| 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 倍高速です。 とは異なりサブ文字列ただし、検索時間は部分文字列の長さに応じて増加します。
最後の反復アプローチでは、含まれています。 このアイデアは、指定された部分文字列を検索し、検索文字列内でその位置をマークし、その位置の 1 文字後ろから開始して部分文字列が見つからなくなるまでこのプロセスを繰り返すというものです。 この方法の魅力は、検索文字列内に部分文字列の出現が少ない場合、検索が非常に高速になることです。 以下のコードは 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;
);
ご覧のとおり、このアプローチは 2 文字の部分文字列のインデックス付け方法よりも開始がわずかに遅くなりますが、部分文字列が長くなるにつれて、つまり、検索する部分文字列が少なくなるにつれて、急速に高速化します。
最後のアプローチでは、私のお気に入りの JSL ツールの 1 つとなっているパターン マッチングを使用します。 コードは速いだけでなく、短いです。
findSequence = Function({inString,inSeq},{i=1,locs={}},
Pat Match(inString,Pat Pos()>>r + inSeq + Pat Test(Insert Into(locs,r+1);0));
locs;
);
2 文字の部分文字列ではわずかに遅くなりますが、3 文字以上の部分文字列の長さの影響はないようです。
明らかに、次を使用した反復的なアプローチは、含まれていますおよびパターン マッチングのアプローチは、長い文字列を検索するための優れた方法です。 しかし、いつ、どの方法で行うのでしょうか? テストはさまざまな部分文字列を使用して行われたため、部分文字列の頻度 (出現回数 * 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 秒も高速でした)。含まれています）。 ガウス過程モデルをデータに当てはめると、検索された文字列の長さの相互作用により、検索された文字列の長さが強い頻度で最大の影響を及ぼしたことが示唆されます。 部分文字列の長さはほとんど影響しませんでした。 最初に線形モデルがデータに適合されましたが、残差分析により、適合に疑問がある可能性があることが示唆されました。 投稿しましたこのモデルとフォローアップ モデルのプロファイラーの対話型バージョンJMPパブリックへ。
2 番目のスペース充填デザインは、以下のケースをより適切に表現するために作成されました。含まれていますパターンマッチングを上回りました。 この実験では、両方の文字列の長さの範囲は以前と変わらず、周波数は 0 と 1 の間で変化しました。 結果は次のとおりですスペース充填結果 2.JMP 。 結果の方向性は以前と同様で、検索文字列の長さが長くなり、頻度が高くなるほどパターン マッチングが有利になりました。
両方の結果セットを組み合わせてガウス過程モデルを当てはめると、50,000 文字の文字列の場合、含まれています周波数が 2 未満の場合に優れています。 100,000 文字では頻度は 1 に低下し、200,000 文字では 0.5 未満になります。 ここからの結論は、よりスパースな短い文字列の場合は、次のようになります。含まれています文字列の長さと密度が増加すると、パターン マッチングがより優れたものになります。
今月のチャレンジ
Covid19 ゲノムを使用して、2 回以上出現する最長の部分文字列を見つけます。
次の投稿まで答えを出さないままにするのではなく、議論して共有しましょう。 これは非常に難しい問題であり、オープンな議論が有益であると私は信じています。
状況によっては、より良い場合もそうでない場合もある別のアプローチを見つけてください (確認する時間がありませんでした)。
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 );
Contains アプローチは、以下のようにより簡潔に書くことができます。 比較が削除されるため、これには、サブ文字列がメイン文字列内に存在しない場合に最も顕著に、実行が少し高速になるという利点もあります。
findSequence = Function({txt, sub},{pos=0,ans={}},
while(pos = contains(txt, sub, pos), insertinto(ans, pos++));
ans;
);
