elynx-markov 0.5.0.1 → 0.5.0.2

raw patch · 9 files changed

+272/−189 lines, 9 filesdep −paralleldep ~asyncdep ~attoparsecdep ~basePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies removed: parallel

Dependency ranges changed: async, attoparsec, base, bytestring, containers, elynx-seq, elynx-tools, hmatrix, hspec, integration, math-functions, mwc-random, primitive, statistics, vector

API changes (from Hackage documentation)

Files

• ChangeLog.md +4/−0
• README.md +29/−29
• bench/Bench.hs +20/−0
• elynx-markov.cabal +27/−20
• src/ELynx/Data/MarkovProcess/AminoAcid.hs +69/−61
• src/ELynx/Data/MarkovProcess/CXXModels.hs +3/−3
• src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs +14/−13
• src/ELynx/Data/MarkovProcess/MixtureModel.hs +24/−27
• src/ELynx/Simulate/MarkovProcessAlongTree.hs +82/−36

@@ -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.

@@ -2,7 +2,7 @@  # The ELynx Suite -Version: 0.5.0.1.+Version: 0.5.0.2. Reproducible evolution made easy.  <p align="center"><img src="https://travis-ci.org/dschrempf/elynx.svg?branch=master"/></p>@@ -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

@@ -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 ()

@@ -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

@@ -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

@@ -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

@@ -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]       ]

@@ -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

@@ -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'