Nussinov78 0.0.1.3 → 0.1.0.0

raw patch · 5 files changed

+209/−208 lines, 5 filesdep +ghc-primdep −BiobaseXNAdep ~ADPfusiondep ~PrimitiveArraydep ~primitivePVP ok

version bump matches the API change (PVP)

Dependencies added: ghc-prim

Dependencies removed: BiobaseXNA

Dependency ranges changed: ADPfusion, PrimitiveArray, primitive, vector

API changes (from Hackage documentation)

Files

• BioInf/GAPlike.hs +164/−0
• BioInf/Nussinov78.hs +0/−167
• C/nussinov.c +9/−6
• Nussinov78.cabal +15/−32
• Nussinov78.hs +21/−3

@@ -0,0 +1,164 @@+{-# LANGUAGE NoMonomorphismRestriction #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE UndecidableInstances #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE PackageImports #-}+{-# LANGUAGE ScopedTypeVariables #-}++module BioInf.GAPlike where++import Control.Monad+import Control.Monad.Primitive+import Control.Monad.ST+import Control.Monad.ST+import Data.Char (toUpper, ord)+import Data.Primitive+import Data.Vector.Fusion.Stream as S+import Data.Vector.Fusion.Stream.Monadic as SM+import Data.Vector.Fusion.Stream.Size+import Data.Vector.Fusion.Util+import Prelude as P+import "PrimitiveArray" Data.Array.Repa.Index+import qualified Data.Vector.Unboxed as VU++import ADP.Fusion.GAPlike+import Data.PrimitiveArray as PA+import Data.PrimitiveArray.Zero.Unboxed as PA++import Debug.Trace+import Control.Arrow (second)++++-- The signature++type Signature m a r =+  ( ()   -> a+  , Char -> a    -> a+  , a    -> Char -> a+  , Char -> a    -> Char -> a+  , a    -> a    -> a+  , SM.Stream m a -> m r+  )++-- the grammar++gNussinov (empty,left,right,pair,split,h) s b e =+  ( s, (  empty <<< e         |||+          left  <<< b % s     |||+          right <<<     s % b |||+          pair  <<< b % s % b |||+          split <<<  s' % s'  ... h)+  ) where s' = transToN s+{-# INLINE gNussinov #-}++-- pairmax algebra++aPairmax :: (Monad m) => Signature m Int Int+aPairmax = (empty,left,right,pair,split,h) where+  empty   _   = 0+  left    b s = s+  right s b   = s+  pair  l s r = if basepair l r then 1+s else -999999+  {-# INLINE [0] pair #-}+  split  l r  = l+r+  h = SM.foldl1' max+  basepair l r = f l r where+    f 'C' 'G' = True+    f 'G' 'C' = True+    f 'A' 'U' = True+    f 'U' 'A' = True+    f 'G' 'U' = True+    f 'U' 'G' = True+    f _   _   = False+  {-# INLINE basepair #-}+{-# INLINE aPairmax #-}++aPretty :: (Monad m) => Signature m String (SM.Stream m String)+aPretty = (empty,left,right,pair,split,h) where+  empty _     = ""+  left  b s   = "." P.++ s+  right   s b = s P.++ "."+  pair  l s r = "(" P.++ s P.++ ")"+  split l   r = l P.++ r+  h = return . id+{-# INLINE aPretty #-}++type CombSignature m e b =+  ( () -> (e, m (SM.Stream m b))+  , Char -> (e, m (SM.Stream m b)) -> (e, m (SM.Stream m b))+  , (e, m (SM.Stream m b)) -> Char -> (e, m (SM.Stream m b))+  , Char -> (e, m (SM.Stream m b)) -> Char -> (e, m (SM.Stream m b))+  , (e, m (SM.Stream m b)) -> (e, m (SM.Stream m b)) -> (e, m (SM.Stream m b))+  , SM.Stream m (e, m (SM.Stream m b)) -> m (SM.Stream m b)+  )++instance Show (Id [String]) where+  show xs = show $ unId xs++(<**)+  :: (Monad m, Eq b, Eq e, Show e, Show (m [b]))+  => Signature m e e+  -> Signature m b (SM.Stream m b)+  -> CombSignature m e b+(<**) f s = (empty,left,right,pair,split,h) where+  (emptyF,leftF,rightF,pairF,splitF,hF) = f+  (emptyS,leftS,rightS,pairS,splitS,hS) = s++  empty e         = (emptyF e   , return $ SM.singleton (emptyS e))+  left b (x,ys)   = (leftF b x  , ys >>= return . SM.map (\y -> leftS b y  ))+  right  (x,ys) b = (rightF x b , ys >>= return . SM.map (\y -> rightS  y b))+  pair l (x,ys) r = (pairF l x r, ys >>= return . SM.map (\y -> pairS l y r))+  split (x,ys) (s,ts) = (splitF x s, ys >>= \ys' -> ts >>= \ts' -> return $ SM.concatMap (\y -> SM.map (\t -> splitS y t) ts') ys')+  h xs = do+    hfs <- hF $ SM.map fst xs+    let phfs = SM.concatMapM snd . SM.filter ((hfs==) . fst) $ xs+    -- trace (">>>" P.++ show (hfs, SM.toList phfs)) $ hS phfs+    hS phfs+++-- * Boilerplate and driver, will be moved to library++nussinov78 inp = (arr ! (Z:.0:.n),bt) where+  (_,Z:._:.n) = bounds arr+  len  = P.length inp+  vinp = VU.fromList . P.map toUpper $ inp+  arr  = runST (nussinov78Fill $ vinp)+  bt   = backtrack vinp arr+{-# NOINLINE nussinov78 #-}++-- type TBL s = Tbl E (PA.MArr0 s DIM2 Int)++nussinov78Fill :: forall s . VU.Vector Char -> ST s (Arr0 DIM2 Int)+nussinov78Fill inp = do+  let n = VU.length inp+  t' <- fromAssocsM (Z:.0:.0) (Z:.n:.n) 0 []+  let t = mtblE t'+  let b = Chr inp+  let e = Empty+  fillTable $ gNussinov aPairmax t b e+  freeze t'+{-# NOINLINE nussinov78Fill #-}++fillTable :: PrimMonad m => (MTbl E (MArr0 (PrimState m) DIM2 Int), ((Int,Int) -> m Int)) -> m ()+fillTable (MTbl tbl, f) = do+  let (_,Z:.n:._) = boundsM tbl+  forM_ [n,n-1..0] $ \i -> forM_ [i..n] $ \j -> do+    v <- f (i,j)+    v `seq` writeM tbl (Z:.i:.j) v+{-# INLINE fillTable #-}++-- * backtracking++backtrack (inp :: VU.Vector Char) (tbl :: PA.Arr0 DIM2 Int) = unId . SM.toList . unId $ g (0,n) where+  n = VU.length inp+  c = Chr inp+  e = Empty+  t = bttblE tbl (g :: BTfun Id String)+  (_,g) = gNussinov (aPairmax <** aPretty) t c e+{-# INLINE backtrack #-}+

@@ -1,167 +0,0 @@-{-# LANGUAGE PackageImports #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE DoAndIfThenElse #-}-{-# LANGUAGE RecordWildCards #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE RankNTypes #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE NoMonomorphismRestriction #-}---- | Strict, scalar nussinov78 algorithm.--module BioInf.Nussinov78 where--import Control.Arrow (first,second,(***))-import Control.Monad-import Control.Monad.Primitive-import Control.Monad.ST-import Control.Monad.State.Lazy-import "PrimitiveArray" Data.Array.Repa.Index-import qualified Data.Vector.Fusion.Stream as P-import qualified Data.Vector.Fusion.Stream.Monadic as S-import qualified Data.Vector.Unboxed as VU--import Biobase.Primary-import Biobase.Secondary.Vienna--import Data.PrimitiveArray-import Data.PrimitiveArray.Unboxed.Zero--import ADP.Fusion-import ADP.Fusion.Monadic-import ADP.Fusion.Monadic.Internal------ | Simple RNA folding with basepair maximization.--nussinov78 inp = arr `seq` bt where-  (_,Z:._:.n) = bounds arr-  arr = runST (nussinov78Fill . mkPrimary $ inp)-  bt  = nussinov78BT (mkPrimary inp) arr-{-# NOINLINE nussinov78 #-}---- | The actual Nussinov78 folding algorithm.--nussinov78Fill :: Primary -> ST s (Arr0 DIM2 Int)-nussinov78Fill inp = do-  let base = base' inp-      {-# INLINE base #-}-  let n = let (_,Z:.l) = bounds inp in l+1-  s <- fromAssocsM (Z:.0:.0) (Z:.n:.n) 0 []--  fillTable s (-                nil   <<< empty               |||-                left  <<< base -~~ s          |||-                right <<<          s ~~- base |||-                pair  <<< base -~~ s ~~- base |||-                split <<<       s +~+ s       ... h-              )-  freeze s-{-# INLINE nussinov78Fill #-}---- | Fill the single table with values in an orderly fashion. The order in--- which we fill depends on the algorithm.--fillTable :: PrimMonad m => MArr0 (PrimState m) DIM2 Int -> (DIM2 -> m Int) -> m ()-fillTable tbl f = do-  let (_,Z:.n:._) = boundsM tbl-  forM_ [n,n-1..0] $ \i -> forM_ [i..n] $ \j -> do-    v <- f (Z:.i:.j)-    writeM tbl (Z:.i:.j) v-    return ()-{-# INLINE fillTable #-}---- | Request the single character enclosed by (i,i+1), with i+1==j--base' :: Primary -> DIM2 -> (Scalar Nuc)-base' inp (Z:.i:.j) = Scalar $ inp ! (Z:.i)-{-# INLINE base' #-}---- | True, if the subword at ij is empty.--empty :: DIM2 -> Scalar Bool-empty (Z:.i:.j) = Scalar $ i==j---- | The base case of our recursion.--nil :: Bool -> Int-nil b = if b then 0 else -999999---- | A single nucleotide to the left. Note that "x" is monadic. In 'nussinov'--- we are in the ST monad, here we just know that we are in a monad.--left :: Nuc -> Int -> Int-left l x = x-{-# INLINE left #-}---- | A single nucleotide to the right.--right :: Int -> Nuc -> Int-right x r = x-{-# INLINE right #-}---- | Pair function--pair :: Nuc -> Int -> Nuc -> Int-pair l x r-  | basepair l r = x+1-  | otherwise    = -999999-{-# INLINE pair #-}---- | Combine the partition of x next-to y.--split :: Int -> Int -> Int-split = (+)-{-# INLINE split #-}---- | Determine if two characters form a legal basepair.--basepair l r-  | mkViennaPair (l,r) /= vpNS = True-basepair _   _   = False-{-# INLINE basepair #-}---- | the grammar makes sure that we at least have "nil #<< empty" in the stream-h = S.foldl1' max-{-# INLINE h #-}------ * backtrace secondary structures--type Backtrace = (Int,String)--nussinov78BT :: Primary -> Arr0 DIM2 Int -> [Backtrace]-nussinov78BT inp s = P.toList $ grammar (Z:.0:.n) where-  base = base' inp-  n = let (_,(Z:._:.l)) = bounds s in l-  s' :: DIM2 -> Scalar (DIM2,Int)-  s' ij = Scalar $ (ij,s!ij)--  grammar :: DIM2 -> P.Stream Backtrace-  grammar = (-      nilBT   <<< empty                |||-      leftBT  <<< base -~~ s'          |||-      rightBT <<<          s' ~~- base |||-      pairBT  <<< base -~~ s' ~~- base |||-      splitBT <<<      s' +~+ s'       ..@ hBT)--  nilBT :: Bool -> (Int, P.Stream Backtrace)-  nilBT b      = if b then (0, P.singleton (0,"")) else (0, P.empty)--  leftBT :: Nuc -> (DIM2,Int) -> (Int, P.Stream Backtrace)-  leftBT _ (ij,x)   = (x, P.map (second ("."++)) $ grammar ij)--  rightBT :: (DIM2,Int) -> Nuc -> (Int, P.Stream Backtrace)-  rightBT (ij,x) _  = (x, P.map (second (++".")) $ grammar ij)--  pairBT :: Nuc -> (DIM2,Int) -> Nuc -> (Int, P.Stream Backtrace)-  pairBT l (ij,x) r = if basepair l r then (x+1, P.map (second (\a -> "("++a++")")) $ grammar ij) else (0, P.empty)--  splitBT :: (DIM2,Int) -> (DIM2,Int) -> (Int, P.Stream Backtrace)-  splitBT (ij,x) (kl,y)  = (x+y, P.concatMap (\(s1,bts1) -> P.map (\(s2,bts2) -> (s1+s2,bts1++bts2)) $ grammar kl) $ grammar ij)--  hBT ij = P.concatMap (\(score,bts) -> P.map (first (const score)) bts) . P.filter (\(score,bts) -> score == s!ij)-

@@ -1,16 +1,19 @@-//#include <stdlib.h>-//#include <stdio.h>+#include <stdlib.h>+#include <stdio.h>+#include <string.h> +int nussinov (int, char *);+ int main () {-  char *p;+  char p[10000];   int n;   int e;-  while (1==scanf ("%as", &p)) { // only GNU C+  while (1==scanf ("%9999s", &p)) { // only GNU C     n = strlen(p);     e = nussinov (n, p);     printf ("%s %d\n", p, e);-    free(p);-    p=0;+    //free(p);+    //p=0;   };   return 0; }

@@ -1,5 +1,5 @@ name:           Nussinov78-version:        0.0.1.3+version:        0.1.0.0 author:         Christian Hoener zu Siederdissen, 2011-2012 copyright:      Christian Hoener zu Siederdissen, 2011-2012 homepage:       http://www.tbi.univie.ac.at/~choener/adpfusion@@ -18,15 +18,14 @@                 .                 This algorithm is simple enough to be used as a tutorial-type                 example. It also shows that efficient code is possible. The-                ADPfusion code compared to C is slower by a factor of only 1.8.+                ADPfusion code compared to C is slower by a factor of only 1.2.                 We plan to improve upon this.                 .                 A number of helper functions currently present in                 BioInf.Nussinov78 will later move in their own library.                 .-                If possible, build using the GHC llvm backend, and GHC-7.2.2.-                GHC-7.4.x produces very bad code on my system, please benchmark-                using 7.2.2.+                Build using GHC-7.6.1, the new code generator and llvm for best+                performance.                 .                 For comparison, a version of the algorithm written in C is                 available under C/nussinov.c. Use at least@@ -37,50 +36,34 @@ Extra-Source-Files:   C/nussinov.c -Flag llvm-  description: build using llvm backend-  default: True  - library   build-depends:     base >= 4 && < 5,     mtl            >= 2,-    primitive      == 0.4.*   ,-    vector         == 0.9.*   ,-    PrimitiveArray == 0.2.2.0 ,-    BiobaseXNA     == 0.6.2.5 ,-    ADPfusion      == 0.0.1.2+    primitive      == 0.5.*   ,+    vector         == 0.10.*  ,+    PrimitiveArray == 0.4.0.0 ,+--    BiobaseXNA     == 0.6.2.5 ,+    ADPfusion      == 0.1.* ,+    ghc-prim   exposed-modules:-    BioInf.Nussinov78+--    BioInf.Nussinov78+    BioInf.GAPlike   ghc-options:-    -O2-    -funbox-strict-fields-    -funfolding-use-threshold100-    -funfolding-keeness-factor100-  if flag (llvm)-    ghc-options:-      -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts+    -O2 -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts    executable Nussinov78   build-depends:---    from above   main-is:     Nussinov78.hs   other-modules:-    BioInf.Nussinov78+--    BioInf.Nussinov78   ghc-options:-    -rtsopts-    -O2-    -funbox-strict-fields-    -funfolding-use-threshold100-    -funfolding-keeness-factor100-  if flag (llvm)-    ghc-options:-      -fllvm -optlo-O3 -optlo-inline -optlo-std-compile-opts+    -fnew-codegen -fllvm -O2 -funbox-strict-fields -optlo-O3 -optlo-std-compile-opts   

@@ -6,16 +6,34 @@ module Main where  import Text.Printf+import System.Environment -import BioInf.Nussinov78+-- import BioInf.Nussinov78+import qualified BioInf.GAPlike as G    main = do-  xs <- fmap lines getContents-  mapM_ doNussinov78 xs+  as <- getArgs+  print as+  case as of+  {-+    [] -> do xs <- fmap lines getContents+             mapM_ doNussinov78 xs -}+    ["gaplike"] -> do xs <- fmap lines getContents+                      mapM_ (doGAPlike 0) xs+    ["gaplike",k] -> do xs <- fmap lines getContents+                        mapM_ (doGAPlike (read k)) xs +{- doNussinov78 inp = do   putStrLn inp   let rs = nussinov78 inp   mapM_ (\(e,bt) -> putStr bt >> printf " %5d\n" e) $ take 10 rs+-}++doGAPlike :: Int -> String -> IO ()+doGAPlike k inp = do+  let (n,bt) = G.nussinov78 inp+  n `seq` printf "%s %d\n" inp n+  mapM_ putStrLn $ take k bt