diff --git a/BioInf/PassiveAggressive.hs b/BioInf/PassiveAggressive.hs --- a/BioInf/PassiveAggressive.hs +++ b/BioInf/PassiveAggressive.hs @@ -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. diff --git a/RNAwolf.cabal b/RNAwolf.cabal --- a/RNAwolf.cabal +++ b/RNAwolf.cabal @@ -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 diff --git a/RNAwolfTrain.hs b/RNAwolfTrain.hs --- a/RNAwolfTrain.hs +++ b/RNAwolfTrain.hs @@ -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" }