2016-08-04 6 views
3

次のコードの効率を改善しようとしています。私は与えられた点の前に(Burrows-Wheeler変換を使ってパターンマッチングの一部として)シンボルのすべての出現を数えたいと思う。どのようにシンボルを数えているかにはいくつかの重複があります。しかし、より効率的なコードでなければならないようなものを実装しようとすると、効率が悪いことが判明しました。私はその怠惰な評価とその貧弱な理解が責任であると仮定しています。リスト処理のためのHaskellの最適化Lazy評価で悩む

カウント機能での私の最初の試みは、次のように行った:マッチング機能自体の本体に続いて

count :: Ord a => [a] -> a -> Int -> Int 
count list sym pos = length . filter (== sym) . take pos $ list 

matching str refCol pattern = match 0 (n - 1) (reverse pattern) 
    where n = length str 
     refFstOcc sym = length $ takeWhile (/= sym) refCol 
     match top bottom [] = bottom - top + 1 
     match top bottom (sym : syms) = 
      let topCt = count str sym top 
       bottomCt = count str sym (bottom + 1) 
       middleCt = bottomCt - topCt 
       refCt = refFstOcc sym 
      in if middleCt > 0 
       then match (refCt + topCt) (refCt + bottomCt - 1) syms 
       else 0 

(簡潔にするためにストリップダウン - 私が最初にmemoizingよrefCol内のシンボルの出現、および他のいくつかの詳細も同様です)。

編集:サンプルの使用は次のようになります。1(私は何かを間違えていなかったと仮定した場合)されなければならない

matching "AT$TCTAGT" "$AACGTTTT" "TCG" 

topポインタとbottomの中間にあるすべての情報を再計算しています。これは、文字列に対して4つの可能な選択肢しかない100万文字のDNA文字列を数えたときに加算されます(プロファイリングではこれが大きなボトルネックも、私の時間の48%をbottomCtに、およそ38%をtopCtのために取っています)。参考のために、100万文字列に対してこれを計算し、50パターン(それぞれ1文字と1000文字の間)を一致させようとすると、プログラムの実行には約8.5〜9.5秒かかります。

しかし、私は次の関数を実装しようとした場合:

countBetween :: Ord a => [a] -> a -> Int -> Int -> (Int, Int) 
countBetween list sym top bottom = 
    let (topList, bottomList) = splitAt top list 
     midList = take (bottom - top) bottomList 
     getSyms = length . filter (== sym) 
    in (getSyms topList, getSyms midList) 

(補償するマッチング機能に加えられた変更で)、プログラムを実行するために18〜22秒かかります。

以前の呼び出しを追跡できるマップを渡そうとしましたが、実行に約20秒かかってしまい、メモリ使用量が増加します。

同様に、私はfoldに短絡しましたが、foldrでは20秒、foldlでは14-15と短絡しています。

だから、このコードを書き直して最適化するには、適切なHaskellの方法は何でしょうか? (具体的には、事前計算に関係しないものを探しています。文字列をあま​​り再利用していない可能性があります。なぜこれが起こっているのかを説明します)。

編集:より明確に、私は以下の通りです探しています何:

A)なぜこの動作はHaskellで起こるのでしょうか?怠惰な評価はどのように役割を果たすのですか?コンパイラがcountcountBetweenの関数を書き直すためにどのような最適化をしていますか?

b)リストを複数回トラバースしないようにこの問題に対処する簡単なコードの書き換えは何ですか?私は、それを回避する解決策ではなく、その問題に対処するものを特に探しています。最終的な答えがcountの場合、コードを書くのに最も効率的な方法ですが、それはなぜですか?

+1

サンプルを入力できますか(大幅に縮小しても)? – pdexter

+0

'長さ。フィルター 'をたやすく倍に短縮することができる。 – ThreeFx

+0

'マッチング'の 'str'パラメータの目的は何ですか? 'n'は' length refCol'やそれに類するはずです。 – ErikR

答えて

1

でも利用可能なコード、私は怠惰な評価は、コードのパフォーマンスと大いに関係しているかわかりません。主な問題は、よりパフォーマンスの高い文字列型ではなく、リンクされたリストであるStringの使用だと思います。

let (topList, bottomList) = splitAt top list 

ます 多くの割り当てを意味topListに対応するリンクされたリンクを再作成します。あなたのcountBetween機能で、この呼び出しは

注意。

splitAttake n/drop n を比較する基準ベンチマークは、http://lpaste.net/174526です。 splitAtのバージョンは で約3倍遅く、もちろん、より多くの割り当てがあります。

カウントを「事前計算したくない」場合でも、 は、ByteStringまたはTextのいずれかに切り替えるだけで、大したことになります。

を定義します。

countSyms :: Char -> ByteString -> Int -> Int -> Int 
countSyms sym str lo hi = 
    length [ i | i <- [lo..hi], BS.index str i == sym ] 

、その後を:

countBetween :: ByteString -> Char -> Int -> Int -> (Int,Int) 
countBetween str sym top bottom = (a,b) 
    where a = countSyms sym str 0 (top-1) 
     b = countSyms sym str top (bottom-1) 

をまた、大きなリストにreverseを使用していない - それは リスト全体を再割り当てします。逆の方法でByteString/Textにインデックスを付けるだけです。

メモの数は役に立つかもしれません。それはどのように行われたかによって異なります。

+0

パーフェクト - それはまさに私が探していたものでした。本当にありがとう! – nomicflux

1

matchルーチンの主なポイントは、現在のシンボルsymに基づいて別の間隔 に間隔(bottom,top)を変換する ようです。式は、基本的に です:

ref_fst = index of sym in ref_col 
    -- defined in an outer scope 

match :: Char -> (Int,Int) -> (Int,Int) 
match sym (bottom, top) | bottom > top = (bottom, top) -- if the empty interval 
match sym (bottom, top) = 
    let 
    top_count = count of sym in str from index 0 to top 
    bot_count = count of sym in str from index 0 to bottom 
    mid_count = top_count - bot_count 
    in if mid_count > 0 
     then (ref_fst + bot_count, ref_fst + top_count) 
     else (1,0) -- the empty interval 

そしてmatchingは初期間隔(0, n-1)match を使用してpatternオーバーちょうど倍です。

top_countbot_count両方が事前に計算ルックアップテーブルを用いて効率的に を計算し、そして以下 はそれを行うコードであることができます。

test1を実行すると、間隔 がパターンの各記号を介してどのように変換されるかのトレースが表示されます。

注:オフ・バイ・1の誤差があるかもしれません、と私はハード0に ref_fstをコード化されてきた - 私はこれは 大きなアルゴリズムにどのように適合するかわからないんだけど、基本的な考え方は、健全でなければなりません。

ベクトルが作成された後でも、 は元の文字列にインデックスを付ける必要はありません。 したがって、ここでは (大)DNA配列のByteStringを使用していますが、それは重要ではありません。 の代わりに mkCountsルーチンも同様に動作します。

http://lpaste.net/174288

{-# LANGUAGE OverloadedStrings #-} 

import Data.Vector.Unboxed ((!)) 
import qualified Data.Vector.Unboxed as UV 
import qualified Data.Vector.Unboxed.Mutable as UVM 
import qualified Data.ByteString.Char8 as BS 
import Debug.Trace 
import Text.Printf 
import Data.List 

mkCounts :: BS.ByteString -> UV.Vector (Int,Int,Int,Int) 
mkCounts syms = UV.create $ do 
    let n = BS.length syms 
    v <- UVM.new (n+1) 
    let loop x i | i >= n = return x 
     loop x i = let s = BS.index syms i 
        (a,t,c,g) = x 
        x' = case s of 
          'A' -> (a+1,t,c,g) 
          'T' -> (a,t+1,c,g) 
          'C' -> (a,t,c+1,g) 
          'G' -> (a,t,c,g+1) 
          _ -> x 
       in do UVM.write v i x 
         loop x' (i+1) 
    x <- loop (0,0,0,0) 0 
    UVM.write v n x 
    return v 

data DNA = A | C | T | G 
    deriving (Show) 

getter :: DNA -> (Int,Int,Int,Int) -> Int 
getter A (a,_,_,_) = a 
getter T (_,t,_,_) = t 
getter C (_,_,c,_) = c 
getter G (_,_,_,g) = g 

-- narrow a window 
narrow :: Int -> UV.Vector (Int,Int,Int,Int) -> DNA -> (Int,Int) -> (Int,Int) 

narrow refcol counts sym (lo,hi) | trace msg False = undefined 
    where msg = printf "-- lo: %d hi: %d refcol: %d sym: %s top_cnt: %d bot_count: %d" lo hi refcol (show sym) top_count bot_count 
     top_count = getter sym (counts ! (hi+1)) 
     bot_count = getter sym (counts ! lo) 

narrow refcol counts sym (lo,hi) = 
    let top_count = getter sym (counts ! (hi+1)) 
     bot_count = getter sym (counts ! (lo+0)) 
     mid_count = top_count - bot_count 
    in if mid_count > 0 
     then (refcol + bot_count, refcol + top_count-1) 
     else (lo+1,lo) -- signal an wmpty window 

findFirst :: DNA -> UV.Vector (Int,Int,Int,Int) -> Int 
findFirst sym v = 
    let n = UV.length v 
     loop i | i >= n = n 
     loop i = if getter sym (v ! i) > 0 
       then i 
       else loop (i+1) 
    in loop 0 

toDNA :: String -> [DNA] 
toDNA str = map charToDNA str 

charToDNA :: Char -> DNA 
charToDNA = go 
    where go 'A' = A 
     go 'C' = C 
     go 'T' = T 
     go 'G' = G 

dnaToChar A = 'A' 
dnaToChar C = 'C' 
dnaToChar T = 'T' 
dnaToChar G = 'G' 

first :: DNA -> BS.ByteString -> Int 
first sym str = maybe len id (BS.elemIndex (dnaToChar sym) str) 
    where len = BS.length str 

test2 = do 
-- matching "AT$TCTAGT" "$AACGTTTT" "TCG" 
    let str = "AT$TCTAGT" 
     refcol = "$AACGTTTT" 
     syms = toDNA "TCG" 

     -- hard coded for now 
     -- may be computeed an memoized 
     refcol_G = 4 
     refcol_C = 3 
     refcol_T = 5 

     counts = mkCounts str 
     w0 = (0, BS.length str -1) 

     w1 = narrow refcol_G counts G w0 
     w2 = narrow refcol_C counts C w1 
     w3 = narrow refcol_T counts T w2 

     firsts = (first A refcol, first T refcol, first C refcol, first G refcol) 

    putStrLn $ "firsts: " ++ show firsts 

    putStrLn $ "w0: " ++ show w0 
    putStrLn $ "w1: " ++ show w1 
    putStrLn $ "w2: " ++ show w2 
    putStrLn $ "w3: " ++ show w3 
    let (lo,hi) = w3 
     len = if lo <= hi then hi - lo + 1 else 0 
    putStrLn $ "length: " ++ show len 

matching :: BS.ByteString -> BS.ByteString -> String -> Int 
matching str refcol pattern = 
    let counts = mkCounts str 
     n = BS.length str 
     syms = toDNA (reverse pattern) 
     firsts = (first A refcol, first T refcol, first C refcol, first G refcol) 

     go (lo,hi) sym = narrow refcol counts sym (lo,hi) 
     where refcol = getter sym firsts 

     (lo, hi) = foldl' go (0,n-1) syms 
     len = if lo <= hi then hi - lo + 1 else 0 
    in len 

test3 = matching "AT$TCTAGT" "$AACGTTTT" "TCG" 
+0

コードが更新されました - 私はあなたのアルゴリズムと互換性があると確信しています。 'test2'と' test3'を試してみてください。また、http://lpaste.net/174288でも利用可能です – ErikR

+0

私の前のコメントを取り戻してください - 明らかに、私のArrayの実装には何か問題がありました。上記と同様の方法で配列と 'unfoldr'を使うと、時間の問題をうまく解決できます。しかし、私は今は記憶上の問題がありますので、これをどうやって行うのかはまだ分かりません。 – nomicflux

+0

多くのことは、あなたがどのように配列を構築するかによって異なります。上記のベクトルコードは、約1/10秒で100万の長さの入力に対するカウントを計算することができます。 – ErikR

関連する問題