7

Haskellを使用してバイオインフォマティクスの領域からモチーフ検索アルゴリズムを実装しています。分枝限定法による中央値文字列検索以外のアルゴリズムの詳細については説明しません。マルチコアの速度を上げるために、並行アプローチ(および後でSTMアプローチ)を実装することによって実装をより面白くすることを計画していましたが、フォローフラグを使用してコンパイルした後

$ ghc -prof -auto-all -O2 -fllvm -threaded -rtsopts --make main 

プロファイルを印刷すると、何か面白い(そしておそらく明らかな)ものが見つかりました:

COST CENTRE      entries  %time %alloc  
hammingDistance  34677951  47.6   14.7  
motifs           4835446   43.8   71.1  

マルチコアプログラミングに近づくことなく、驚くべきスピードアップが得られることは明らかです(ただし、それは実行されており、いくつかの優れたテストデータを見つけて、そのための基準を整理する必要があります)。

とにかく、これらの機能は両方とも純粋に機能的であり、決して同時ではありません。とてもシンプルなこともしているので、時間がかかってびっくりしました。それらのコードは次のとおりです。

data NukeTide = A | T | C | G deriving (Read, Show, Eq, Ord, Enum)

type Motif = [NukeTide] 

hammingDistance :: Motif -> Motif -> Int
hammingDistance [] [] = 0
hammingDistance xs [] = 0 -- optimistic
hammingDistance [] ys = 0 -- optimistic
hammingDistance (x:xs) (y:ys) = case (x == y) of
    True  -> hammingDistance xs ys
    False -> 1 + hammingDistance xs ys

motifs :: Int -> [a] -> [[a]]
motifs n nukeTides = [ take n $ drop k nukeTides | k <- [0..length nukeTides - n] ]

hammingDistanceに対する2つの引数のうち、改善の余地がある場合は、xsがxの長さになり、ysがそれ以下になると実際に想定できることに注意してください。

ご覧のとおり、hammingDistanceは、ヌクレオチドのリストである2つのモチーフ間のハミング距離を計算します。モチーフ関数は数値とリストを受け取り、その長さのすべてのサブ文字列を返します。例:

> motifs 3 "hello world"
["hel","ell","llo","lo ","o w"," wo","wor","orl","rld"]

関連するアルゴリズムプロセスは非常に単純なので、これをさらに最適化する方法を考えることはできません。しかし、私はどこに向かうべきかについて2つの推測があります。

  1. HammingDistance:私が使用しているデータ型(NukeTidesと[])は遅い/不器用です。私はそれらの実装に精通していないので、これは単なる推測ですが、自分のデータ型を定義することは、より読みやすくなりますが、おそらく私が意図するよりも多くのオーバーヘッドを伴うと思います。また、パターンマッチングは私にとって異質であり、それが些細なことなのか、費用がかかるのかはわかりません。
  2. モチーフ:これを正しく読んでいると、すべてのメモリ割り当ての70%がモチーフによって行われ、いつかガベージコレクションが必要になると思います。繰り返しになりますが、万能リストを使用すると、コストが非常に不明確になるため、速度が低下したり、リストの理解が遅くなったりする可能性があります。

ここでの通常の手順について誰かアドバイスはありますか?データ型が問題である場合、配列は正しい答えでしょうか?(箱に入っていると聞きました)

助けてくれてありがとう。

編集:これらの2つの関数が呼び出される方法を説明すると便利かもしれないと思いました。

totalDistance :: Motif -> Int
totalDistance motif = sum $ map (minimum . map (hammingDistance motif) . motifs l) dna

この関数は別の関数の結果であり、ツリー内のノードに渡されます。ツリーの各ノードで、totalDistanceを使用してノードをスコアリングし、ヌクレオチドの評価(長さ<= n、つまり== nの場合はリーフノード)が実行されます。それ以降は、典型的な分枝限定アルゴリズムになります。

編集:ジョンは、モチーフのコストを実質的に排除するために行った変更を印刷するように依頼しました。

scoreFunction :: DNA -> Int -> (Motif -> Int)
scoreFunction dna l = totalDistance
    where
        -- The sum of the minimum hamming distance in each line of dna
        -- is given by totalDistance motif
        totalDistance motif = sum $ map (minimum . map (hammingDistance motif)) possibleMotifs
        possibleMotifs = map (motifs l) dna -- Previously this was computed in the line above

元の投稿では明確にしませんでしたが、scoreFunctionは1回だけ呼び出され、結果はツリートラバーサル/分枝限定法で渡され、バインドされてノードの評価に使用されます。振り返ってみると、あらゆる段階でモチーフを再計算することは、私が行った中で最も明るいことの1つではありません。

4

3 に答える 3

7

の各アプリケーションは最初からリストをトラバースする必要があるため、の定義はmotifs必要以上にトラバースを実行しているように見えます。代わりに次dropを使用して実装します。Data.List.tails

motifs2 :: Int -> [a] -> [[a]]
motifs2 n nukeTides = map (take n) $ take count $ tails nukeTides
  where count = length nukeTides - n + 1

GHCiで簡単に比較すると、違いがわかります(sum . map length評価を強制するために使用)。

*Main> let xs = concat (replicate 10000 [A, T, C, G])
(0.06 secs, 17914912 bytes)
*Main> sum . map length $ motifs 5 xs
199980
(3.47 secs, 56561208 bytes)
*Main> sum . map length $ motifs2 5 xs
199980
(0.15 secs, 47978952 bytes)
于 2011-10-09T15:36:54.053 に答える
6

の定義hammingDistanceは、おそらくそれよりもはるかに効率的ではありません。

hammingDistance (x:xs) (y:ys) = case (x == y) of
    True  -> hammingDistance xs ys
    False -> 1 + hammingDistance xs ys

haskellの怠惰のため、これは(最悪の場合)次のように拡張されます。

(1 + (1 + (1 + ...)))

これはスタック上にサンクとして存在し、使用された場合にのみ減少します。これが実際に問題であるかどうかは、呼び出しサイト、コンパイラオプションなどに依存するため、この問題を完全に回避する形式でコードを記述することをお勧めします。

一般的な解決策は、厳密なアキュムレータを使用して末尾再帰形式を作成することですが、この場合、次のような高階関数を使用できます。

hammingDistance :: Motif -> Motif -> Int
hammingDistance xs ys = length . filter (uncurry (==)) $ zip xs ys

比較のために、ここに末尾再帰の実装があります

hammingDistance :: Motif -> Motif -> Int
hammingDistance xs ys = go 0 xs ys
  where
    go !acc [] [] = acc
    go !acc xs [] = acc -- optimistic
    go !acc [] ys = acc -- optimistic
    go !acc (x:xs) (y:ys) = case (x == y) of
      True  -> go acc xs ys
      False -> go (acc+1) xs ys

これはBangPatterns拡張機能を使用してアキュムレータを厳密に評価するように強制します。そうでない場合、現在の定義と同じ問題が発生します。

他の質問のいくつかに直接答えるには:

  1. パターンマッチングは簡単です
  2. リストと配列のどちらを使用するかは、主にデータの作成方法と消費方法によって異なります。この場合、リストが最適なタイプである可能性があります。特に、リストが作成されたときにすべて消費され、リスト全体をメモリに保存する必要がない場合は、問題ないはずです。ただし、リストをメモリに保持すると、スペースのオーバーヘッドが大きくなります。

使用パターン

これらの関数を使用する方法は、いくつかの追加の作業も行うと思います。

(minimum . map (hammingDistance motif) . motifs l

必要なのは最小値だけなので、不要hammingDistanceな余分な値をたくさん計算している可能性があります。私はこれに対する2つの解決策を考えることができます:

オプション1.hammingDistanceThresh :: Motif -> Int -> Motif -> Intしきい値を超えると停止する新しい関数を定義します。少し奇妙なタイプの順序は、次のように折りたたんで使用しやすくすることです。

let motifs' = motifs l
in foldl' (hammingDistanceThresh motif) (hammingDistance motif $ head motifs') (tail motifs')

オプション2。怠惰な自然数タイプを定義する場合Int、の結果にsの代わりにそれを使用できますhammingDistance。次に、必要なだけのハミング距離が計算されます。

最後の注意:を使用-auto-allすると、他のプロファイリングオプションよりもはるかに遅いコードが頻繁に生成されます。-auto最初に使用してから、必要に応じて手動のSCC注釈を追加することをお勧めします。

于 2011-10-09T19:41:57.320 に答える
2

そうです...私は限界に達することに抵抗できず、プレーンメタルっぽいパックビットの実装を書きました:

{-# language TypeSynonymInstances #-}
{-# language BangPatterns #-}

import Data.Bits
import Data.Word


data NukeTide = A | T | C | G deriving (Read, Show, Eq, Ord, Enum)

type UnpackedMotif = [NukeTide] 

type PackageType = Word32
nukesInPackage = 16 :: Int
allSetMask = complement 0 :: PackageType


-- Be careful to have length of motif == nukesInPackage here!
packNukesToWord :: UnpackedMotif -> PackageType
packNukesToWord = packAt 0
    where packAt _ [] = 0
          packAt i (m:ml) =    (b0 m .&. bit i)
                           .|. (b1 m .&. bit (i+1))
                           .|. packAt (i+2) ml
          b0 A = 0
          b0 T = allSetMask
          b0 C = 0
          b0 G = allSetMask
          b1 A = 0
          b1 T = 0
          b1 C = allSetMask
          b1 G = allSetMask

unpackNukesWord :: PackageType -> UnpackedMotif
unpackNukesWord = unpackNNukesFromWord nukesInPackage

unpackNNukesFromWord :: Int -> PackageType -> UnpackedMotif
unpackNNukesFromWord = unpackN
    where unpackN 0 _ = []
          unpackN i w = (nukeOf $ w .&. r2Mask):(unpackN (i-1) $ w`shiftR`2)
          nukeOf bs
           | bs == 0      = A
           | bs == bit 0  = T
           | bs == bit 1  = C
           | otherwise    = G
          r2Mask = (bit 1 .|. bit 0) :: PackageType


data PackedMotif = PackedMotif { motifPackets::[PackageType]
                               , nukesInLastPack::Int        }
 -- note nukesInLastPack will never be zero; motifPackets must be [] to represent empty motifs.
packNukes :: UnpackedMotif -> PackedMotif
packNukes m = case remain of
               [] -> PackedMotif [packNukesToWord takeN] (length takeN)
               r  -> prAppend (packNukesToWord takeN) (packNukes r)
    where (takeN, remain) = splitAt nukesInPackage m
          prAppend w (PackedMotif l i) = PackedMotif (w:l) i

unpackNukes :: PackedMotif -> UnpackedMotif
unpackNukes (PackedMotif l i) = unpack l i
  where unpack [l] i = unpackNNukesFromWord i l
        unpack (l:ls) i = unpackNukesWord l ++ unpack ls i
        unpack [] _ = []

instance Show PackedMotif where
  show = show . unpackNukes



class Nukes a where
  pLength :: a -> Int
  shiftLN1 :: a -> a
  hammingDistance :: a -> a -> Int
  motifs :: Int -> a -> [a]

instance Nukes PackageType where
  pLength _ = nukesInPackage
  shiftLN1 = (`shiftR`2)
  hammingDistance !x !y = fromIntegral $ abt (x `xor` y)
      where abt !b = bbt(b.&.a0Mask .|. ((b.&.a1Mask) `shiftR` 1))
            bbt !b = sbt $ (b.&.r16Mask) + (b `shiftR` nukesInPackage)
            sbt !b = (r2Mask .&. b)             + (r2Mask .&. (b`shiftR`2))
                   + (r2Mask .&. (b`shiftR`4))  + (r2Mask .&. (b`shiftR`6))
                   + (r2Mask .&. (b`shiftR`8))  + (r2Mask .&. (b`shiftR`10))
                   + (r2Mask .&. (b`shiftR`12)) + (r2Mask .&. (b`shiftR`14))
            a0Mask = 0x55555555 :: PackageType
            a1Mask = 0xAAAAAAAA :: PackageType
            r16Mask = 0xFFFF :: PackageType
            r2Mask = 0x3 :: PackageType
  motifs 0 _ = []
  motifs l x = x : motifs (l-1) (shiftLN1 x)


maskNukesBut :: Int -> PackageType -> PackageType
maskNukesBut i = ( ( allSetMask `shiftR` (2*(nukesInPackage - i)) ) .&.)

instance Nukes PackedMotif where
  pLength (PackedMotif (x:xs) ix) = nukesInPackage * (length xs) + ix
  pLength _ = 0
  shiftLN1 ξ@(PackedMotif [] _) = ξ
  shiftLN1 (PackedMotif [x] ix) | ix>1       = PackedMotif [x`shiftR`2] (ix-1)
                                | otherwise  = PackedMotif [] nukesInPackage
  shiftLN1 (PackedMotif (x:x':xs) ix)
        = PackedMotif (( shiftLN1 x .|. pnext ):sxs) resLMod
      where sxs = motifPackets $ shiftLN1 (PackedMotif (x':xs) ix)
            pnext = shiftL (x'.&.0x3) 30
            resLMod = if ix > 1 then (ix-1) else nukesInPackage
  hammingDistance xs ys = go 0 xs ys
    where
      go :: Int -> PackedMotif -> PackedMotif -> Int
      go !acc (PackedMotif [x] ix) (PackedMotif [y] iy)
       | ix > iy    = acc + (hammingDistance y $ maskNukesBut iy x)
       | otherwise  = acc + (hammingDistance x $ maskNukesBut ix y)
      go !acc (PackedMotif [x] ix) (PackedMotif (y:ys) iy)
        = acc + (hammingDistance x $ maskNukesBut ix y)
      go !acc (PackedMotif (x:xs) ix) (PackedMotif [y] iy)
        = acc + (hammingDistance y $ maskNukesBut iy x)
      go !acc (PackedMotif (x:xs) ix) (PackedMotif (y:ys) iy)
        = go (acc + hammingDistance x y) (PackedMotif xs ix) (PackedMotif ys iy)
      go !acc _ _ = acc
  motifs l ξ
     | l>0        = fShfts (min nukesInPackage $ pLength ξ + 1 - l) ξ >>= ct
     | otherwise  = []
    where fShfts k χ | k > 0      = χ : fShfts (k-1) (shiftLN1 χ)
                     | otherwise  = []
          ct (PackedMotif ys iy) = case remain of
                [] -> if (length takeN - 1) * nukesInPackage + iy >= l
                       then [PackedMotif takeN lMod] else []
                _  -> PackedMotif takeN lMod : ct(PackedMotif (tail ys) iy)
            where (takeN, remain) = splitAt lQuot ys
          (lQuot,lMod) = case l `quotRem` nukesInPackage of
                   (i,0) -> (i, nukesInPackage)
                   (i,m) -> (i+1, m)

これは、関数を使用してUnpackedMotif = [NukeTide]sから使用できます。packNukes

*BioNuke0> motifs 23 $ packNukes $ take 27 $ cycle [A,T,G,C,A]
[[A,T,G,C,A,A,T,G,C,A,A,T,G,C,A,A,T,G,C,A,A,T,G],[T,G,C,A,A,T,G,C,A,A,T,G,C,A,A,T,G,C,A,A,T,G,C],[G,C,A,A,T,G,C,A,A,T,G,C,A,A,T,G,C,A,A,T,G,C,A],[C,A,A,T,G,C,A,A,T,G,C,A,A,T,G,C,A,A,T,G,C,A,A],[A,A,T,G,C,A,A,T,G,C,A,A,T,G,C,A,A,T,G,C,A,A,T]]

*BioNuke0> hammingDistance (packNukes [A,T,G,C,A,A,T,G]) (packNukes [A,T,C,C,A,T,G])
3

*BioNuke0> map (hammingDistance (packNukes $ take 52 $ cycle [A,T,C,C,A,T,G])) (motifs 52 $ packNukes $ take 523 $ cycle [A,T,C,C,A,T,G])
[0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44,38,52,0,52,36,45,44]

パフォーマンスを元のバージョンとまだ比較していませんが、代数的データ型の実装よりもかなり高速であるはずです。さらに、スペース効率の高いストレージ形式を簡単に提供します。

于 2011-10-10T02:12:44.340 に答える