diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -4,6 +4,10 @@ ## Unreleased changes + +## Version 0.5.0.2 + +- Speed up mixture model simulation. - Improve rooting functions. - Improve `Topology` data type (but still a lot to do). - Various additions to the documentation. diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -2,7 +2,7 @@ # The ELynx Suite -Version: 0.5.0.1. +Version: 0.5.0.2. Reproducible evolution made easy.

@@ -73,9 +73,9 @@ slynx --help | head -n -16 - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] [-f|--force] [--no-elynx-file] COMMAND @@ -90,7 +90,7 @@ Specify base name of output file -f,--force Ignore previous analysis and overwrite existing output files. - --no-elynx-file Do not write files for needed for reproducibility. + --no-elynx-file Do not write data required to reproduce an analysis. Available commands: concatenate Concatenate sequences found in input files. @@ -120,9 +120,9 @@ slynx concatenate --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE Concatenate sequences found in input files. @@ -141,9 +141,9 @@ slynx examine --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site] Examine sequences. If data is a multi sequence alignment, additionally analyze columns. @@ -163,9 +163,9 @@ slynx filter-rows --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH] [--shorter-than LENGTH] [--standard-characters] @@ -186,9 +186,9 @@ slynx filter-columns --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE [--standard-chars DOUBLE] @@ -210,9 +210,9 @@ slynx simulate --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL] [-m|--mixture-model MODEL] [-e|--edm-file NAME] @@ -289,9 +289,9 @@ slynx sub-sample --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE (-n|--number-of-sites INT) @@ -319,9 +319,9 @@ slynx translate --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT) (-u|--universal-code CODE) @@ -344,9 +344,9 @@ tlynx --help | head -n -16 - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME] [-f|--force] [--no-elynx-file] COMMAND @@ -361,7 +361,7 @@ Specify base name of output file -f,--force Ignore previous analysis and overwrite existing output files. - --no-elynx-file Do not write files for needed for reproducibility. + --no-elynx-file Do not write data required to reproduce an analysis. Available commands: compare Compare two phylogenetic trees (compute distances and branch-wise differences). @@ -384,9 +384,9 @@ tlynx compare --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect] [-f|--newick-format FORMAT] NAMES @@ -413,9 +413,9 @@ tlynx examine --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT] Compute summary statistics of phylogenetic trees. @@ -437,9 +437,9 @@ tlynx simulate --help - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: tlynx simulate (-t|--nTrees INT) (-n|--nLeaves INT) PROCESS [-u|--sub-sample DOUBLE] [-s|--summary-statistics] @@ -474,9 +474,9 @@ elynx --help | head -n -16 - ELynx Suite version 0.5.0.1. + ELynx Suite version 0.5.0.2. Developed by Dominik Schrempf. - Compiled on December 18, 2020, at 10:14 am, UTC. + Compiled on February 19, 2021, at 13:42 pm, UTC. Usage: elynx COMMAND Validate and redo past ELynx analyses diff --git a/bench/Bench.hs b/bench/Bench.hs new file mode 100644 --- /dev/null +++ b/bench/Bench.hs @@ -0,0 +1,20 @@ +-- | +-- Module : Main +-- Description : Benchmark Markov process along tree +-- Copyright : (c) 2021 Dominik Schrempf +-- License : GPL-3.0-or-later +-- +-- Maintainer : dominik.schrempf@gmail.com +-- Stability : experimental +-- Portability : portable +-- +-- Creation date: Wed Feb 17 17:42:30 2021. +module Main + ( main, + ) +where + +-- Stub. + +main :: IO () +main = return () diff --git a/elynx-markov.cabal b/elynx-markov.cabal --- a/elynx-markov.cabal +++ b/elynx-markov.cabal @@ -1,6 +1,6 @@ cabal-version: 2.2 name: elynx-markov -version: 0.5.0.1 +version: 0.5.0.2 license: GPL-3.0-or-later license-file: LICENSE copyright: Dominik Schrempf (2020) @@ -48,20 +48,19 @@ default-language: Haskell2010 ghc-options: -Wall -Wunused-packages build-depends: - async >=2.2.2 && <2.3, - attoparsec >=0.13.2.4 && <0.14, + async >=2.2.2, + attoparsec >=0.13.2.5, base >=4.7 && <5, - bytestring >=0.10.10.0 && <0.11, - containers >=0.6.2.1 && <0.7, - elynx-seq >=0.5.0.1 && <0.6, - hmatrix >=0.20.0.0 && <0.21, - integration >=0.2.1 && <0.3, - math-functions >=0.3.4.1 && <0.4, - mwc-random >=0.14.0.0 && <0.15, - parallel >=3.2.2.0 && <3.3, - primitive >=0.7.1.0 && <0.8, - statistics >=0.15.2.0 && <0.16, - vector >=0.12.1.2 && <0.13 + bytestring >=0.10.12.0, + containers >=0.6.2.1, + elynx-seq >=0.5.0.2, + hmatrix >=0.20.1, + integration >=0.2.1, + math-functions >=0.3.4.1, + mwc-random >=0.14.0.0, + primitive >=0.7.1.0, + statistics >=0.15.2.0, + vector >=0.12.2.0 test-suite markov-test type: exitcode-stdio-1.0 @@ -80,10 +79,18 @@ ghc-options: -Wall -Wunused-packages build-depends: base >=4.7 && <5, - containers >=0.6.2.1 && <0.7, + containers >=0.6.2.1, elynx-markov -any, - elynx-tools >=0.5.0.1 && <0.6, - hmatrix >=0.20.0.0 && <0.21, - hspec >=2.7.4 && <2.8, - mwc-random >=0.14.0.0 && <0.15, - vector >=0.12.1.2 && <0.13 + elynx-tools >=0.5.0.2, + hmatrix >=0.20.1, + hspec >=2.7.8, + mwc-random >=0.14.0.0, + vector >=0.12.2.0 + +benchmark markov-bench + type: exitcode-stdio-1.0 + main-is: Bench.hs + hs-source-dirs: bench + default-language: Haskell2010 + ghc-options: -Wall -Wunused-packages + build-depends: base >=4.14.1.0 diff --git a/src/ELynx/Data/MarkovProcess/AminoAcid.hs b/src/ELynx/Data/MarkovProcess/AminoAcid.hs --- a/src/ELynx/Data/MarkovProcess/AminoAcid.hs +++ b/src/ELynx/Data/MarkovProcess/AminoAcid.hs @@ -12,18 +12,22 @@ -- -- The order of amino acids is alphabetic. module ELynx.Data.MarkovProcess.AminoAcid - ( lg, + ( -- * Amino acid substitution models + lg, lgCustom, wag, wagCustom, poisson, poissonCustom, gtr20, + + -- * Convenience functions + alphaToPamlVec, + pamlToAlphaVec, ) where import Data.ByteString.Internal (c2w) -import Data.List (elemIndex) import Data.Maybe (fromMaybe) import qualified Data.Vector.Storable as V import Data.Word (Word8) @@ -39,80 +43,84 @@ -- 'pamlToAlphaVec' and 'pamlToAlphaMat'. -- Amno acids in alphabetical order. -aaAlphaOrder :: [Word8] +aaAlphaOrder :: V.Vector Word8 aaAlphaOrder = - map + V.map c2w - [ 'A', - 'C', - 'D', - 'E', - 'F', - 'G', - 'H', - 'I', - 'K', - 'L', - 'M', - 'N', - 'P', - 'Q', - 'R', - 'S', - 'T', - 'V', - 'W', - 'Y' - ] + $ V.fromList + [ 'A', + 'C', + 'D', + 'E', + 'F', + 'G', + 'H', + 'I', + 'K', + 'L', + 'M', + 'N', + 'P', + 'Q', + 'R', + 'S', + 'T', + 'V', + 'W', + 'Y' + ] -- Amino acids in PAML oder. -aaPamlOrder :: [Word8] +aaPamlOrder :: V.Vector Word8 aaPamlOrder = - map + V.map c2w - [ 'A', - 'R', - 'N', - 'D', - 'C', - 'Q', - 'E', - 'G', - 'H', - 'I', - 'L', - 'K', - 'M', - 'F', - 'P', - 'S', - 'T', - 'W', - 'Y', - 'V' - ] - --- -- This is a very slow implementation; since I only convert matrices once it --- -- should not be a problem. A map would be better if performance is an issue. --- pamlIndexToAlphaIndex :: Int -> Int --- pamlIndexToAlphaIndex i = fromMaybe --- (error $ "Could not convert index " ++ show i ++ ".") --- (elemIndex aa aaAlphaOrder) --- where aa = aaPamlOrder !! i + $ V.fromList + [ 'A', + 'R', + 'N', + 'D', + 'C', + 'Q', + 'E', + 'G', + 'H', + 'I', + 'L', + 'K', + 'M', + 'F', + 'P', + 'S', + 'T', + 'W', + 'Y', + 'V' + ] --- This is a very slow implementation; since I only convert matrices once it --- should not be a problem. A map would be better if performance is an issue. alphaIndexToPamlIndex :: Int -> Int alphaIndexToPamlIndex i = fromMaybe (error $ "Could not convert index " ++ show i ++ ".") - (elemIndex aa aaPamlOrder) + (V.elemIndex aa aaPamlOrder) where - aa = aaAlphaOrder !! i + aa = aaAlphaOrder V.! i --- Convert an amino acid vector in PAML order to a vector in alphabetical order. +pamlIndexToAlphaIndex :: Int -> Int +pamlIndexToAlphaIndex i = + fromMaybe + (error $ "Could not convert index " ++ show i ++ ".") + (V.elemIndex aa aaAlphaOrder) + where + aa = aaPamlOrder V.! i + +-- | Convert an amino acid vector in PAML order to a vector in alphabetical order. pamlToAlphaVec :: Vector R -> Vector R pamlToAlphaVec v = build n (\i -> v ! alphaIndexToPamlIndex (round i)) + +-- | Convert an amino acid vector in alphabetical order to a vector in PAML order. +alphaToPamlVec :: Vector R -> Vector R +alphaToPamlVec v = build n (\i -> v ! pamlIndexToAlphaIndex (round i)) -- Convert an amino acid matrix in PAML order to a matrix in alphabetical order. pamlToAlphaMat :: Matrix R -> Matrix R diff --git a/src/ELynx/Data/MarkovProcess/CXXModels.hs b/src/ELynx/Data/MarkovProcess/CXXModels.hs --- a/src/ELynx/Data/MarkovProcess/CXXModels.hs +++ b/src/ELynx/Data/MarkovProcess/CXXModels.hs @@ -19,7 +19,7 @@ ) where -import Data.List.NonEmpty (fromList) +import qualified Data.Vector as V import ELynx.Data.MarkovProcess.AminoAcid import ELynx.Data.MarkovProcess.CXXModelsData import qualified ELynx.Data.MarkovProcess.MixtureModel as M @@ -139,8 +139,8 @@ cxxFromStatDistsAndWeights ws ds = M.fromSubstitutionModels (cxxName n) - (fromList ws) + (V.fromList ws) sms where n = length ds - sms = fromList $ cxxSubstitutionModelsFromStatDists ds + sms = V.fromList $ cxxSubstitutionModelsFromStatDists ds diff --git a/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs b/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs --- a/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs +++ b/src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs @@ -20,9 +20,8 @@ ) where +import qualified Data.Vector as V import qualified Data.ByteString.Lazy.Char8 as BL -import Data.List.NonEmpty (NonEmpty) -import qualified Data.List.NonEmpty as N import qualified ELynx.Data.MarkovProcess.MixtureModel as M import qualified ELynx.Data.MarkovProcess.PhyloModel as P import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S @@ -60,20 +59,20 @@ ++ show alpha splitSubstitutionModel :: - Int -> Double -> S.SubstitutionModel -> NonEmpty S.SubstitutionModel + Int -> Double -> S.SubstitutionModel -> V.Vector S.SubstitutionModel splitSubstitutionModel n alpha sm = renamedSMs where means = getMeans n alpha - scaledSMs = N.map (`S.scale` sm) means - names = N.fromList $ map (("; gamma rate category " ++) . show) [1 :: Int ..] - renamedSMs = N.zipWith S.appendName names scaledSMs + scaledSMs = V.map (`S.scale` sm) means + names = V.fromList $ map (("; gamma rate category " ++) . show) [1 :: Int ..] + renamedSMs = V.zipWith S.appendName names scaledSMs expandSubstitutionModel :: Int -> Double -> S.SubstitutionModel -> M.MixtureModel expandSubstitutionModel n alpha sm = M.fromSubstitutionModels name ws sms where name = S.name sm <> getName n alpha - ws = N.repeat 1.0 + ws = V.replicate n 1.0 sms = splitSubstitutionModel n alpha sm expandMixtureModel :: Int -> Double -> M.MixtureModel -> M.MixtureModel @@ -81,15 +80,17 @@ where name = M.name mm <> getName n alpha means = getMeans n alpha - scaledMMs = N.map (`M.scale` mm) means - names = N.fromList $ map (("; gamma rate category " ++) . show) [1 :: Int ..] - renamedMMs = N.zipWith M.appendNameComponents names scaledMMs + scaledMMs = V.map (`M.scale` mm) means + names = V.fromList $ map (("; gamma rate category " ++) . show) [1 :: Int ..] + renamedMMs = V.zipWith M.appendNameComponents names scaledMMs -- For a given number of rate categories 'n' and a shape parameter 'alpha' (the -- rate or scale is set such that the mean is 1.0), return a list of rates that -- represent the respective categories. Use the mean rate for each category. -getMeans :: Int -> Double -> NonEmpty Double -getMeans n alpha = N.fromList $ means <> pure lastMean +getMeans :: Int -> Double -> V.Vector Double +getMeans n alpha + | n <= 0 = error "getMeans: Number of rate categories is zero or negative." + | otherwise = means <> pure lastMean where gamma = gammaDistr alpha (1.0 / alpha) quantiles = @@ -99,7 +100,7 @@ -- probability mass is 1/n. meanFunc x = fromIntegral n * x * density gamma x -- Only calculate the first (n-1) categories with normal integration. - means = + means = V.fromList [ integralAToB meanFunc (quantiles !! i) (quantiles !! (i + 1)) | i <- [0 .. n - 2] ] diff --git a/src/ELynx/Data/MarkovProcess/MixtureModel.hs b/src/ELynx/Data/MarkovProcess/MixtureModel.hs --- a/src/ELynx/Data/MarkovProcess/MixtureModel.hs +++ b/src/ELynx/Data/MarkovProcess/MixtureModel.hs @@ -32,8 +32,7 @@ ) where -import qualified Data.List.NonEmpty as N -import Data.Semigroup +import qualified Data.Vector as V import ELynx.Data.Alphabet.Alphabet hiding (all) import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S import Prelude @@ -53,40 +52,38 @@ { -- | Name name :: S.Name, alphabet :: Alphabet, - components :: N.NonEmpty Component + components :: V.Vector Component } deriving (Show, Read) -- | Get weights. -getWeights :: MixtureModel -> N.NonEmpty Weight -getWeights = N.map weight . components +getWeights :: MixtureModel -> V.Vector Weight +getWeights = V.map weight . components -- | Get substitution models. -getSubstitutionModels :: MixtureModel -> N.NonEmpty S.SubstitutionModel -getSubstitutionModels = N.map substModel . components +getSubstitutionModels :: MixtureModel -> V.Vector S.SubstitutionModel +getSubstitutionModels = V.map substModel . components -- | Create a mixture model from a list of substitution models. -fromSubstitutionModels :: - S.Name -> N.NonEmpty Weight -> N.NonEmpty S.SubstitutionModel -> MixtureModel -fromSubstitutionModels n ws sms = - if allEqual $ N.toList alphs - then MixtureModel n (N.head alphs) comps - else - error - "fromSubstitutionModels: alphabets of substitution models are not equal." +fromSubstitutionModels :: S.Name -> V.Vector Weight -> V.Vector S.SubstitutionModel -> MixtureModel +fromSubstitutionModels n ws sms + | null ws = error "fromSubstitutionModels: No weights given." + | length ws /= length sms = error "fromSubstitutionModels: Number of weights and substitution models does not match." + | not $ allEqual alphs = error "fromSubstitutionModels: alphabets of substitution models are not equal." + | otherwise = MixtureModel n (V.head alphs) comps where - comps = N.zipWith Component ws sms - alphs = N.map S.alphabet sms - allEqual [] = True - allEqual xs = all (== head xs) $ tail xs + comps = V.zipWith Component ws sms + alphs = V.map S.alphabet sms + allEqual xs | V.null xs = True + | otherwise = V.all (== V.head xs) xs -- | Concatenate mixture models. -concatenate :: S.Name -> N.NonEmpty MixtureModel -> MixtureModel +concatenate :: S.Name -> V.Vector MixtureModel -> MixtureModel concatenate n mms = fromSubstitutionModels n ws sms where - comps = sconcat $ N.map components mms - ws = N.map weight comps - sms = N.map substModel comps + comps = V.concatMap components mms + ws = V.map weight comps + sms = V.map substModel comps scaleComponent :: Double -> Component -> Component scaleComponent s c = c {substModel = s'} where s' = S.scale s $ substModel c @@ -96,16 +93,16 @@ scale s m = m {components = cs'} where cs = components m - cs' = N.map (scaleComponent s) cs + cs' = V.map (scaleComponent s) cs -- | Globally normalize a mixture model so that on average one event happens per -- unit time. normalize :: MixtureModel -> MixtureModel normalize mm = scale (1 / c) mm where - c = sum $ N.zipWith (*) weights scales + c = sum $ V.zipWith (*) weights scales weights = getWeights mm - scales = N.map S.totalRate $ getSubstitutionModels mm + scales = V.map S.totalRate $ getSubstitutionModels mm appendNameComponent :: S.Name -> Component -> Component appendNameComponent n c = c {substModel = s'} @@ -117,4 +114,4 @@ appendNameComponents n m = m {components = cs'} where cs = components m - cs' = N.map (appendNameComponent n) cs + cs' = V.map (appendNameComponent n) cs diff --git a/src/ELynx/Simulate/MarkovProcessAlongTree.hs b/src/ELynx/Simulate/MarkovProcessAlongTree.hs --- a/src/ELynx/Simulate/MarkovProcessAlongTree.hs +++ b/src/ELynx/Simulate/MarkovProcessAlongTree.hs @@ -18,6 +18,7 @@ ( -- * Single rate matrix. simulate, simulateAndFlatten, + simulateAndFlattenPar, -- * Mixture models. simulateMixtureModel, @@ -30,16 +31,26 @@ import Control.Concurrent.Async import Control.Monad import Control.Monad.Primitive -import Control.Parallel.Strategies import Data.Tree -import qualified Data.Vector.Storable as V +import qualified Data.Vector as V import Data.Word (Word32) import ELynx.Data.MarkovProcess.RateMatrix import ELynx.Simulate.MarkovProcess -import Numeric.LinearAlgebra import System.Random.MWC import System.Random.MWC.Distributions (categorical) +-- XXX @performace. The horizontal concatenation might be slow. If so, +-- 'concatenateSeqs' or 'concatenateAlignments' can be used, which directly +-- appends vectors. + +-- A brain f***. As an example, let @xss@ be a list of alignments (i.e., a list +-- of a list of a list of alleles). This function horizontally concatenates the +-- sites. The number of species needs to be same in each alignment. No checks +-- are performed! +horizontalConcat :: [[[a]]] -> [[a]] +horizontalConcat [xs] = xs +horizontalConcat xss = foldr1 (zipWith (++)) xss + toProbTree :: RateMatrix -> Tree Double -> Tree ProbMatrix toProbTree q = fmap (probMatrix q) @@ -54,8 +65,9 @@ -- | Simulate a number of sites for a given substitution model. Only the states -- at the leafs are retained. The states at internal nodes are removed. This has -- a lower memory footprint. --- --- XXX: Improve performance. Use vectors, not lists. + +-- XXX: Improve performance. Use vectors, not lists. I am actually not sure if +-- this improves performance... simulateAndFlatten :: PrimMonad m => Int -> @@ -85,6 +97,31 @@ then return [is'] else concat <$> sequence [simulateAndFlatten' is' t g | t <- f] +-- | See 'simulateAndFlatten', parallel version. +simulateAndFlattenPar :: + Int -> + StationaryDistribution -> + ExchangeabilityMatrix -> + Tree Double -> + GenIO -> + IO [[State]] +simulateAndFlattenPar n d e t g = do + c <- getNumCapabilities + gs <- splitGen c g + let chunks = getChunks c n + q = fromExchangeabilityMatrix e d + pt = toProbTree q t + -- The concurrent map returns a list of [[State]] objects. They have to be + -- concatenated horizontally. + sss <- + mapConcurrently + ( \(num, gen) -> do + is <- getRootStates num d gen + simulateAndFlatten' is pt gen + ) + (zip chunks gs) + return $ horizontalConcat sss + -- | Simulate a number of sites for a given substitution model. Keep states at -- internal nodes. The result is a tree with the list of simulated states as -- node labels. @@ -116,20 +153,21 @@ return $ Node is' f' toProbTreeMixtureModel :: - [RateMatrix] -> Tree Double -> Tree [ProbMatrix] + V.Vector RateMatrix -> Tree Double -> Tree (V.Vector ProbMatrix) toProbTreeMixtureModel qs = - fmap (\a -> [probMatrix q a | q <- qs] `using` parList rpar) + -- XXX: This function is slow. Parallelization? + fmap (\a -> V.map (`probMatrix` a) qs) getComponentsAndRootStates :: PrimMonad m => Int -> - Vector R -> - [StationaryDistribution] -> + V.Vector Double -> + V.Vector StationaryDistribution -> Gen (PrimState m) -> m ([Int], [State]) getComponentsAndRootStates n ws ds g = do cs <- replicateM n $ categorical ws g - is <- sequence [categorical (ds !! c) g | c <- cs] + is <- sequence [categorical (ds V.! c) g | c <- cs] return (cs, is) -- | Simulate a number of sites for a given set of substitution models with @@ -138,27 +176,29 @@ simulateAndFlattenMixtureModel :: PrimMonad m => Int -> - Vector R -> - [StationaryDistribution] -> - [ExchangeabilityMatrix] -> + V.Vector Double -> + V.Vector StationaryDistribution -> + V.Vector ExchangeabilityMatrix -> Tree Double -> Gen (PrimState m) -> - m [[State]] + -- | (IndicesOfComponents, [SimulatedSequenceForEachTip]) + m ([Int], [[State]]) simulateAndFlattenMixtureModel n ws ds es t g = do - let qs = zipWith fromExchangeabilityMatrix es ds + let qs = V.zipWith fromExchangeabilityMatrix es ds pt = toProbTreeMixtureModel qs t (cs, is) <- getComponentsAndRootStates n ws ds g - simulateAndFlattenMixtureModel' is cs pt g + ss <- simulateAndFlattenMixtureModel' is cs pt g + return (cs, ss) simulateAndFlattenMixtureModel' :: (PrimMonad m) => [State] -> [Int] -> - Tree [ProbMatrix] -> + Tree (V.Vector ProbMatrix) -> Gen (PrimState m) -> m [[State]] simulateAndFlattenMixtureModel' is cs (Node ps f) g = do - is' <- sequence [jump i (ps !! c) g | (i, c) <- zip is cs] + is' <- sequence [jump i (ps V.! c) g | (i, c) <- zip is cs] if null f then return [is'] else @@ -189,22 +229,28 @@ -- | See 'simulateAndFlattenMixtureModel', parallel version. simulateAndFlattenMixtureModelPar :: Int -> - Vector R -> - [StationaryDistribution] -> - [ExchangeabilityMatrix] -> + V.Vector Double -> + V.Vector StationaryDistribution -> + V.Vector ExchangeabilityMatrix -> Tree Double -> GenIO -> - IO [[[State]]] + IO ([Int], [[State]]) simulateAndFlattenMixtureModelPar n ws ds es t g = do - let qs = zipWith fromExchangeabilityMatrix es ds + let qs = V.zipWith fromExchangeabilityMatrix es ds pt = toProbTreeMixtureModel qs t - parComp - n - ( \n' g' -> do - (cs, is) <- getComponentsAndRootStates n' ws ds g' - simulateAndFlattenMixtureModel' is cs pt g' - ) - g + -- The concurrent computation returns a list of ([Int], [[State]]) objects. + -- They have to be concatenated horizontally. + csss <- + parComp + n + ( \n' g' -> + do + (cs, is) <- getComponentsAndRootStates n' ws ds g' + ss <- simulateAndFlattenMixtureModel' is cs pt g' + return (cs, ss) + ) + g + return (concatMap fst csss, horizontalConcat $ map snd csss) -- | Simulate a number of sites for a given set of substitution models with -- corresponding weights. Keep states at internal nodes. See also @@ -212,14 +258,14 @@ simulateMixtureModel :: PrimMonad m => Int -> - Vector R -> - [StationaryDistribution] -> - [ExchangeabilityMatrix] -> + V.Vector Double -> + V.Vector StationaryDistribution -> + V.Vector ExchangeabilityMatrix -> Tree Double -> Gen (PrimState m) -> m (Tree [State]) simulateMixtureModel n ws ds es t g = do - let qs = zipWith fromExchangeabilityMatrix es ds + let qs = V.zipWith fromExchangeabilityMatrix es ds pt = toProbTreeMixtureModel qs t (cs, is) <- getComponentsAndRootStates n ws ds g simulateMixtureModel' is cs pt g @@ -230,10 +276,10 @@ (PrimMonad m) => [State] -> [Int] -> - Tree [ProbMatrix] -> + Tree (V.Vector ProbMatrix) -> Gen (PrimState m) -> m (Tree [State]) simulateMixtureModel' is cs (Node ps f) g = do - is' <- sequence [jump i (ps !! c) g | (i, c) <- zip is cs] + is' <- sequence [jump i (ps V.! c) g | (i, c) <- zip is cs] f' <- sequence [simulateMixtureModel' is' cs t g | t <- f] return $ Node is' f'