RNAwolf 0.3.1.0 → 0.3.2.0

raw patch · 3 files changed

+163/−108 lines, 3 filesdep +splitPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: split

API changes (from Hackage documentation)

Files

• BioInf/PassiveAggressive.hs +65/−34
• RNAwolf.cabal +12/−3
• RNAwolfTrain.hs +86/−71

@@ -16,42 +16,64 @@  module BioInf.PassiveAggressive where -import qualified Data.Vector.Unboxed as VU-import Data.List as L-import Data.Set as S import Control.Arrow+import Control.DeepSeq+import Control.Parallel (pseq)+import Data.List as L import Data.Map as M+import Data.Set as S+import qualified Data.Vector.Unboxed as VU import Text.Printf  import Biobase.TrainingData import BioInf.Keys  import qualified BioInf.Params as P-import qualified BioInf.Params.Import as P import qualified BioInf.Params.Export as P+import qualified BioInf.Params.Import as P  import Statistics.ConfusionMatrix import Statistics.PerformanceMetrics -import Data.PrimitiveArray as PA-import Data.PrimitiveArray.Ix  ---- | Default implementation of P/A.+-- | Default implementation of P/A. We return a data structure that contains+-- all 'changes' required from this run, the 'enerDif' or energy difference+-- between the known and the predicted structure, and a structural difference+-- score. Furthermore, some errors are being reported in 'errors'.+--+-- The energy difference can be (i) in that case, a wrongly predicted structure+-- has better (lower) energy than the known one. (ii) It can be zero, then we+-- have either found a co-optimal structural or the correct structure. (iii) In+-- some cases, it can be positive, this is a formal error, but will not abort+-- the program. (The calling program may opt to abort on (not . null $ errors).+--+-- The structural difference is [0..1] with "0" for structurally identical+-- predictions and known structures and otherwise growing toward "1" for bad+-- predictions where nothing is correct. -defaultPA :: Double -> P.Params -> TrainingData -> (P.Params,Double,Double,[(Int,Double)])+defaultPA :: Double -> P.Params -> TrainingData -> PA defaultPA aggressiveness params td@TrainingData{..}-  | L.null $ pOnly++kOnly = (params,0,1,[])-  | sty >= 0.999 = (params,0,1,[])-  | otherwise = ( new-                , tau-                , sty-                , changes-                )+  | L.null $ pOnly++kOnly = PA+      { changes = []+      , enerDif = edif+      , accMeas = struc+      , errors  = []+      }+  | struc >= 0.999 = PA+      { changes = []+      , enerDif = edif+      , accMeas = struc+      , errors  = []+      }+  | otherwise = PA+      { changes = changes+      , enerDif = edif+      , accMeas = struc+      , errors  = eError+      }   where-    -- create new vector-    new = P.fromList . VU.toList $ VU.accum (\v pm -> v+pm) cur changes+    -- calculate changes     pFeatures = featureVector primary predicted     kFeatures = featureVector primary secondary     pOnly = pFeatures L.\\ kFeatures@@ -61,24 +83,16 @@     cur = VU.fromList . P.toList $ params     pScore = sum . L.map (cur VU.!) $ pFeatures     kScore = sum . L.map (cur VU.!) $ kFeatures-    -- weight calculation-    tau-      | kScore + epsilon < pScore-          = error $ "S(known) < S(predicted)\n" ++ errorKnownTooGood td cur kFeatures pFeatures-      |  sty >  0.999-      && kScore+epsilon < pScore-          = error $ "S(known) < S(predicted)\n" ++ errorKnownTooGood td cur kFeatures pFeatures-      | sty >= 0.999-          = 0-      | otherwise-          = val-      where-        val = min aggressiveness $ (kScore - pScore + sqrt (1-sty)) / (numChanges ^ 2)-    sty = case fmeasure (mkConfusionMatrix td) of -- currently optimizing using F_1+    edif = kScore - pScore+    eError = if edif <= 0+               then []+               else ["S(known) < S(predicted)\n" ++ errorKnownTooGood td cur kFeatures pFeatures]+    struc = case fmeasure (mkConfusionMatrix td) of -- currently optimizing using F_1             Left  _ -> 1             Right v -> v-    -- special constants-    epsilon = 0.1+    -- weight calculation+    tau = min aggressiveness $ ( (min 0 $ kScore - pScore) + sqrt (1-struc)+                               ) / (numChanges ^ 2)  -- | In case that the known structure has a score 'epsilon' better than the -- predicted, we have an error condition, as this should never be the case.@@ -89,6 +103,23 @@     ++ printf "%s\n%s\n" primary (concat $ intersperse "\n" comments)   kScore = sum . L.map (curPs VU.!) $ kFeatures   pScore = sum . L.map (curPs VU.!) $ pFeatures++-- | Return a lot of information from each P/A call. We do not return the new+-- 'Params' anymore, only a list of changes. This allows us to do some things.+-- If the implementation of 'Params' is switched, we can update in place; or we+-- can perform calculations in parallel.++data PA = PA+  { changes :: [(Int,Double)] -- (index of change, change)+  , enerDif :: Double         -- how much pressure from a wrong energy difference+  , accMeas :: Double         -- accuracy measure+  , errors  :: [String]       -- if something strange happens+  } deriving (Show)++instance NFData PA where+  rnf PA{..} = rnf changes `pseq` rnf enerDif `pseq` rnf accMeas `pseq` rnf errors++-- * Instances  -- | Pull in the statistical interface. From the confusion matrix, we -- automagically get everything we need.

@@ -1,5 +1,5 @@ name:           RNAwolf-version:        0.3.1.0+version:        0.3.2.0 author:         Christian Hoener zu Siederdissen, Stephan H Bernhart, Peter F Stadler, Ivo L Hofacker copyright:      Christian Hoener zu Siederdissen, 2010-2011 homepage:       http://www.tbi.univie.ac.at/software/rnawolf/@@ -48,8 +48,16 @@                 changes. Please send a mail, if you encounter strange behaviour                 or bugs.                 .-                Last Changes:+                Changes in 0.3.2.0                 .+                * simpler training wrapper+                .+                * added parallelism option for multi-core systems (reduce+                  iteration time for the cost of a possible reduction in+                  training efficiency; but should be worth it)+                .+                Changes in 0.3.1.0+                .                 * fixed bugs introduced by bulge/interior/multi-loops  Flag llvm@@ -95,11 +103,12 @@  executable RNAwolfTrain   build-depends:+    split,     cmdargs == 0.7.*   main-is:     RNAwolfTrain.hs   ghc-options:-    -O2 -rtsopts+    -O2 -rtsopts -threaded   if flag(llvm)     ghc-options:       -fllvm

@@ -26,31 +26,34 @@  module Main where +import Control.Applicative+import Control.Arrow import Control.Monad-import System.Console.CmdArgs-import Text.Printf-import Data.List+import Control.Parallel (pseq)+import Control.Parallel.Strategies import Data.Function (on)-import System.Random-import Control.Applicative+import Data.List+import Data.List.Split (splitEvery) import Data.Ord-import Control.Arrow-import qualified Data.Vector.Unboxed as VU import qualified Data.Map as M+import qualified Data.Vector.Unboxed as VU+import System.Console.CmdArgs+import System.Random+import Text.Printf  import Biobase.Primary+import Biobase.Secondary.Diagrams import Biobase.TrainingData import Biobase.TrainingData.Import import Statistics.ConfusionMatrix import Statistics.PerformanceMetrics-import Biobase.Secondary.Diagrams +import BioInf.Keys import BioInf.Params as P import BioInf.Params.Export as P import BioInf.Params.Import as P-import BioInf.RNAwolf import BioInf.PassiveAggressive-import BioInf.Keys+import BioInf.RNAwolf   @@ -89,86 +92,85 @@  doIteration :: Options -> [TrainingData] -> (P.Params,[Double]) -> Int -> IO (P.Params,[Double]) doIteration o@Options{..} xs' (!p,rhos) !k = do-  xs <- shuffle xs'+  xs <- fmap (splitEvery parallelism) $ shuffle xs'   when (Iteration `elem` verbose) $ do     putStrLn "\n======================================"     printf "# INFO iteration: %4d / %4d starting\n"             k             numIterations+    printf "# INFO folding %d elements, maximal length: %d\n"+            (length xs')+            (maximum $ map (length . primary) xs')     putStrLn "======================================\n"-  (newp,totalchange,rhosum,cooptimality) <- foldM (foldTD o $ length xs) (p,0,0,0) $ zip xs [1..]+  let indices = mapAccumL (\acc x -> (acc+x,(acc+1,acc+x))) 0 $ map length xs+  (newp,rs) <- foldM (foldTD o $ length xs) (p,[]) . zip xs . snd $ indices   let drctch = sum $ zipWith (\x y -> abs $ x-y) (P.toList p) (P.toList newp)-  let rho = rhosum / genericLength xs+  let rhosum = sum $ map accMeas rs+  let rho = rhosum / (sum . map genericLength $ xs)   when (Iteration `elem` verbose) $ do     putStrLn "\n======================================"-    printf "# INFO iteration: %4d / %4d ended\n"+    printf "# INFO iteration: %4d / %4d ended, rho: %4.2f (%4.2f, %4.2f)\n"             k             numIterations-    printf "# INFO sum tau: %7.2f, total change: %7.2f, avg.rho: %4.2f, avg.coopt: %5d\n"-            totalchange-            drctch             rho-            (cooptimality `div` length xs)-    putStr "# INFO history:"+            (minimum $ map accMeas rs)+            (maximum $ map accMeas rs)+    putStr "# INFO rho history:"     zipWithM_ (printf " %4d %4.2f") [1::Int ..] $ rhos++[rho]     putStrLn ""-    {--    print $ sum $ map abs $ P.toList newp-    print $ minimum $ P.toList newp-    print $ maximum $ P.toList newp-    -}     putStrLn "======================================\n"   writeFile (printf "%04d.db" k) . show $ newp   return (newp,rhos++[rho]) --- | Fold one 'TrainingData', print some info and stuff+-- | Fold one 'TrainingData' element and return the suggested changes and+-- additional information. -foldTD :: Options -> Int -> (P.Params,Double,Double,Int) -> (TrainingData,Int) -> IO (P.Params,Double,Double,Int)-foldTD o@Options{..} total (!p,accChange,rhosum,cooptimality) (td@TrainingData{},k) = do-  let pri = mkPrimary $ primary td-  let tables = rnaWolf p pri-  let bs' = let f x = td{predicted = x} in-            map (first f) -            . take (maybe 1 id maxLoss)-            $ rnaWolfBacktrack p pri 0.001 tables-  printf "co-opts: %d\n" $ length bs'---  mapM_ print bs'---  print $ rnaWolfOptimal tables-  let bs = pure $ minimumBy (comparing (fmeasure . mkConfusionMatrix . fst)) bs'-  case bs of-    [(x,ddd)] -> do-      let fV = featureVector (primary x) (predicted x)-      let pVU = VU.fromList . P.toList $ p-      let sss = map (pVU VU.!) fV-      when (abs (ddd - sum sss) > 0.001) $ do-        printf "SCORE DIFFERENCE, backtracking score: %f, sum features: %f\n"  ddd   (sum sss) -- , " ", map (vks M.!) fV, " ", sss)-        mapM_ print $ zip (map (vks M.!) fV) sss-        print "You have found a bug, now write choener to have him fix it!"-      let (newp,tau,rho,votes) = defaultPA aggressiveness p-                                $ x { comments = [ show ddd-                                                 , simpleViewer (primary x) $ secondary x-                                                 , simpleViewer (primary x) $ predicted x-                                                 , show $ predicted x-                                                 ]-                                    }-      when (Single `elem` verbose) $ do-        printf "# INFO currently at: %4d / %4d (s-tau: %7.4f, rho: %5.2f, changes: %4d)\n"-                k-                total-                (tau * genericLength votes)-                rho-                (length votes)-      when (Detailed `elem` verbose) $ do-        putStrLn $ take (length $ primary x) . concatMap show . concat . repeat $ [0..9]-        putStrLn $ primary x-        putStrLn $ simpleViewer (primary x) $ secondary x-        putStrLn $ simpleViewer (primary x) $ predicted x-        when (AllPairs `elem` verbose) $ do-          mapM_ print $ predicted x-      when (errorOnError && abs (ddd - sum sss) > 0.0001) $ error "error-ing out"-      return (newp,accChange + tau * genericLength votes, rhosum+rho, cooptimality + length bs')-    _       -> error $ "no prediction for: " ++ show td+foldOne :: Options -> P.Params -> TrainingData -> (TrainingData,PA)+foldOne o@Options{..} p td+  | null bs   = (td   , PA [] 0 0 ["no prediction for: " ++ primary td])+  | otherwise = (fst worst, ret)+  where+    pri = mkPrimary $ primary td+    tables = rnaWolf p pri+    bs = let f x = td{predicted = x} in+         map (first f) . take (maybe 1 id maxLoss) $ rnaWolfBacktrack p pri 0.001 tables+    worst = minimumBy (comparing (fmeasure . mkConfusionMatrix . fst)) bs+    runPA (x,score) = defaultPA aggressiveness p+            $ x { comments =+                    [ show score+                    , simpleViewer (primary x) $ secondary x+                    , simpleViewer (primary x) $ predicted x+                    , show $ predicted x+                    ]+                }+    ret = runPA worst +-- | Folding of 'TrainingData' elements.++foldTD :: Options -> Int -> (P.Params,[PA]) -> ([TrainingData],(Int,Int)) -> IO (P.Params,[PA])+foldTD o@Options{..} total (!p,oldresults) (ts,(f,t)) = do+  -- At this point, we trade most efficient optimization with increased parallelism, if that option is >1+  let parfolds = map (foldOne o p) ts+  let !results = let xs = map snd parfolds in xs `using` (parList rdeepseq)+  let cs = concatMap changes results+  let cur = VU.fromList . P.toList $ p+  let new = P.fromList . VU.toList $ VU.accum (\v pm -> v+pm) cur cs+  let rhosum = sum $ map accMeas results+  let rho = rhosum / genericLength ts+  let rhosumR = sum . map accMeas $ oldresults ++ results+  let rhoR = rhosumR / genericLength (oldresults ++ results)+  when (Single `elem` verbose) $ do+    printf "# INFO parallel: %4d - %4d, avg.rho: %4.2f, running rho: %4.2f\n"+            f t+            rho+            rhoR+  -- detailed information on each folded structure+  mapM_ (printDetailed o . fst) parfolds+  return $ pseq (rdeepseq results)+         ( new+         , oldresults ++ results+         )+ -- | simple viewer...  simpleViewer s xs = foldl f (replicate (length s) '.') xs where@@ -187,8 +189,19 @@       l = head $ drop k str       post = drop (k+1) str +-- | print out detailed information on a folded candidate +printDetailed :: Options -> TrainingData -> IO ()+printDetailed Options{..} x = do+  when (Detailed `elem` verbose) $ do+    putStrLn $ take (length $ primary x) . concatMap show . concat . repeat $ [0..9]+    putStrLn $ primary x+    putStrLn $ simpleViewer (primary x) $ secondary x+    putStrLn $ simpleViewer (primary x) $ predicted x+    when (AllPairs `elem` verbose) $ do+      mapM_ print $ predicted x + -- ** program options  data Options = Options@@ -201,6 +214,7 @@   , maxLoss :: Maybe Int   , aggressiveness :: Double   , errorOnError :: Bool+  , parallelism :: Int   } deriving (Show,Data,Typeable)  data Verbose@@ -216,10 +230,11 @@   , trainingData = [] &= help "training data elements to read"   , maxLength = Nothing &= help "[dev] only train using elements of length or less"   , numIterations = 50 &= help "how many optimizer iterations"-  , verbose = [Iteration] &= help "select verbosity options: single, iteration, detailed (all switch on different verbosity options)"+  , verbose = [] &= help "select verbosity options: single, iteration, detailed (all switch on different verbosity options)"   , maxLoss = Nothing &= help "use maxLoss optimization instead of prediction-based, requires maximal number of instances to search for maxLoss (default: not used)"   , aggressiveness = 1 &= help "maximal tau for each round"   , errorOnError = False &= help "error out if an error is detected (default: false)"+  , parallelism = 1 &= help "perform more than one prediction concurrently. Will probably reduce the effectiveness of the algorithm but allow to use more than one core; call with +RTS -N -RTS"   }