packages feed

elynx-markov 0.3.0 → 0.3.1

raw patch · 16 files changed

+206/−230 lines, 16 filesdep +asyncdep +attoparsecdep −elynx-treedep −hspec-megaparsecdep −megaparsec

Dependencies added: async, attoparsec

Dependencies removed: elynx-tree, hspec-megaparsec, megaparsec

Files

ChangeLog.md view
@@ -2,7 +2,14 @@ # Changelog for ELynx  -## Version 0.3.1; in development+## Version 0.3.2; in development+++## Version 0.3.1++-   Use Attoparsec.+-   Use ByteString consistently.+-   Remove elynx-tools dependency from libaries.   ## Version 0.3.0
README.md view
@@ -2,7 +2,7 @@  # The ELynx Suite -Version: 0.3.0.+Version: 0.3.1. Reproducible evolution made easy.  <p align="center"><img src="https://travis-ci.org/dschrempf/elynx.svg?branch=master"/></p>@@ -90,9 +90,9 @@      slynx --help | head -n -16 -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME]                   [-f|--force] COMMAND@@ -136,9 +136,9 @@      slynx concatenate --help -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: slynx concatenate (-a|--alphabet NAME) INPUT-FILE       Concatenate sequences found in input files.@@ -157,9 +157,9 @@      slynx examine --help -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: slynx examine (-a|--alphabet NAME) INPUT-FILE [--per-site]       Examine sequences. If data is a multi sequence alignment, additionally analyze columns.@@ -179,9 +179,9 @@      slynx filter-rows --help -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: slynx filter-rows (-a|--alphabet NAME) INPUT-FILE [--longer-than LENGTH]                               [--shorter-than LENGTH] [--standard-characters]@@ -202,9 +202,9 @@      slynx filter-columns --help -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: slynx filter-columns (-a|--alphabet NAME) INPUT-FILE                                  [--standard-chars DOUBLE]@@ -226,9 +226,9 @@      slynx simulate --help -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL]                            [-m|--mixture-model MODEL] [-e|--edm-file NAME] @@ -305,9 +305,9 @@      slynx sub-sample --help -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: slynx sub-sample (-a|--alphabet NAME) INPUT-FILE                             (-n|--number-of-sites INT)@@ -335,9 +335,9 @@      slynx translate --help -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: slynx translate (-a|--alphabet NAME) INPUT-FILE (-r|--reading-frame INT)                            (-u|--universal-code CODE)@@ -360,9 +360,9 @@      tlynx --help | head -n -16 -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: tlynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME]                   [-f|--force] COMMAND@@ -399,9 +399,9 @@      tlynx compare --help -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: tlynx compare [-n|--normalize] [-b|--bipartitions] [-t|--intersect]                           [-f|--newick-format FORMAT] NAMES@@ -428,9 +428,9 @@      tlynx examine --help -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: tlynx examine INPUT-FILE [-f|--newick-format FORMAT]       Compute summary statistics of phylogenetic trees.@@ -452,9 +452,9 @@      tlynx simulate --help -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: tlynx simulate [-t|--nTrees INT] [-n|--nLeaves INT] [-H|--height DOUBLE]                            [-M|--condition-on-mrca] [-l|--lambda DOUBLE] @@ -491,9 +491,9 @@      elynx --help | head -n -16 -    ELynx Suite version 0.3.0.+    ELynx Suite version 0.3.1.     Developed by Dominik Schrempf.-    Compiled on July 30, 2020, at 12:40 pm, UTC.+    Compiled on August 3, 2020, at 12:12 pm, UTC.          Usage: elynx COMMAND       Validate and redo past ELynx analyses
elynx-markov.cabal view
@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: d6c9194049e2cfd023ea390aefa6b761f0297f5484539c466a7f5b646e156877+-- hash: c05f62a676b8e0fa2f4741dbd0fc99542f112c2421694e8e1b3a5e59ab5f235b  name:           elynx-markov-version:        0.3.0+version:        0.3.1 synopsis:       Simulate molecular sequences along trees description:    Examine, modify, and simulate molecular sequences in a reproducible way. Please see the README on GitHub at <https://github.com/dschrempf/elynx>. category:       Bioinformatics@@ -50,16 +50,15 @@       src   ghc-options: -Wall -fllvm   build-depends:-      base >=4.7 && <5+      async+    , attoparsec+    , base >=4.7 && <5     , bytestring     , containers     , elynx-seq-    , elynx-tools-    , elynx-tree     , hmatrix     , integration     , math-functions-    , megaparsec     , mwc-random     , parallel     , primitive@@ -86,10 +85,8 @@     , containers     , elynx-markov     , elynx-tools-    , elynx-tree     , hmatrix     , hspec-    , hspec-megaparsec     , mwc-random     , vector   default-language: Haskell2010
src/ELynx/Data/MarkovProcess/AminoAcid.hs view
@@ -22,13 +22,14 @@   ) 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) import ELynx.Data.Alphabet.Alphabet import ELynx.Data.MarkovProcess.RateMatrix import ELynx.Data.MarkovProcess.SubstitutionModel-import ELynx.Tools import Numeric.LinearAlgebra  n :: Int@@ -321,10 +322,16 @@ lgExch :: ExchangeabilityMatrix lgExch = pamlToAlphaMat $ exchFromListLower n lgExchRawPaml +normalizeSumVec :: Vector Double -> Vector Double+normalizeSumVec v = V.map (/ s) v+  where+    s = V.sum v+{-# INLINE normalizeSumVec #-}+ -- Stationary distribution in PAML order. lgStatDistPaml :: StationaryDistribution lgStatDistPaml =-  normalizeSumVec 1.0 $+  normalizeSumVec $     fromList       [ 0.079066,         0.055941,@@ -564,7 +571,7 @@ -- WAG stationary distribution in PAML order. wagStatDistPaml :: StationaryDistribution wagStatDistPaml =-  normalizeSumVec 1.0 $+  normalizeSumVec $     fromList       [ 0.0866279,         0.043972,@@ -602,15 +609,22 @@   where     nm = fromMaybe "WAG-Custom" mnm +matrixSetDiagToZero :: Matrix R -> Matrix R+matrixSetDiagToZero m = m - diag (takeDiag m)+{-# INLINE matrixSetDiagToZero #-}+ uniformExch :: ExchangeabilityMatrix uniformExch = matrixSetDiagToZero $ matrix n $ replicate (n * n) 1.0  poissonExch :: ExchangeabilityMatrix poissonExch = uniformExch +uniformVec :: Vector Double+uniformVec = V.replicate n (1 / fromIntegral n)+ -- | Poisson substitution model. poisson :: SubstitutionModel-poisson = substitutionModel Protein "Poisson" [] (uniformVec n) poissonExch+poisson = substitutionModel Protein "Poisson" [] uniformVec poissonExch  -- | Poisson substitution model with maybe a name and a custom stationary distribution. poissonCustom :: Maybe String -> StationaryDistribution -> SubstitutionModel
src/ELynx/Data/MarkovProcess/CXXModels.hs view
@@ -10,7 +10,7 @@ -- -- Creation date: Tue Feb 26 16:44:33 2019. ----- Quang, L. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for+-- Quang, BL. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for -- phylogenetic reconstruction. Bioinformatics 24, 2317–2323 (2008). -- -- XXX: For now, I only provide Poisson exchangeabilities.
src/ELynx/Data/MarkovProcess/CXXModelsData.hs view
@@ -10,7 +10,7 @@ -- -- Creation date: Tue Feb 26 17:17:35 2019. ----- Quang, L. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for+-- Quang, BL. S., Gascuel, O. & Lartillot, N. Empirical profile mixture models for -- phylogenetic reconstruction. Bioinformatics 24, 2317–2323 (2008). module ELynx.Data.MarkovProcess.CXXModelsData   ( c10StatDists,
src/ELynx/Data/MarkovProcess/GammaRateHeterogeneity.hs view
@@ -20,7 +20,7 @@   ) where -import qualified Data.ByteString.Lazy.Char8 as L+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@@ -32,10 +32,10 @@ import Prelude hiding (repeat)  -- | Short summary of gamma rate heterogeneity parameters.-summarizeGammaRateHeterogeneity :: Int -> Double -> [L.ByteString]+summarizeGammaRateHeterogeneity :: Int -> Double -> [BL.ByteString] summarizeGammaRateHeterogeneity n alpha =   map-    L.pack+    BL.pack     [ "Discrete gamma rate heterogeneity.",       "Number of categories: " ++ show n,       "Shape parameter of gamma distribution: " ++ show alpha,@@ -106,9 +106,8 @@     -- The last category has to be calculated with an improper integration.     lastMean = integralAToInf meanFunc (quantiles !! (n - 1)) --- The error of integration. eps :: Double-eps = 1e-6+eps = 1e-8  -- The integration method to use method :: (Double -> Double) -> Double -> Double -> [Result]
src/ELynx/Data/MarkovProcess/MixtureModel.hs view
@@ -14,11 +14,10 @@ module ELynx.Data.MarkovProcess.MixtureModel   ( -- * Types     Weight,-    Component,-    MixtureModel (name),+    Component (weight, substModel),+    MixtureModel (name, alphabet, components),      -- * Getters-    getAlphabet,     getWeights,     getSubstitutionModels, @@ -30,26 +29,14 @@     scale,     normalize,     appendNameComponents,--    -- * Output-    summarizeComponent,-    summarize,   ) where -import qualified Data.ByteString.Builder as L-import qualified Data.ByteString.Lazy.Char8 as L-import Data.List.NonEmpty hiding (zip)+import qualified Data.List.NonEmpty as N import Data.Semigroup-import ELynx.Data.Alphabet.Alphabet+import ELynx.Data.Alphabet.Alphabet hiding (all) import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S-import ELynx.Tools-import Prelude hiding-  ( head,-    length,-    map,-    zipWith,-  )+import Prelude  -- | Mixture model component weight. type Weight = Double@@ -66,43 +53,40 @@   { -- | Name     name :: S.Name,     alphabet :: Alphabet,-    components :: NonEmpty Component+    components :: N.NonEmpty Component   }   deriving (Show, Read) --- | Get alphabet used with mixture model. Throws error if components use--- different 'Alphabet's.-getAlphabet :: MixtureModel -> Alphabet-getAlphabet = alphabet- -- | Get weights.-getWeights :: MixtureModel -> NonEmpty Weight-getWeights = map weight . components+getWeights :: MixtureModel -> N.NonEmpty Weight+getWeights = N.map weight . components  -- | Get substitution models.-getSubstitutionModels :: MixtureModel -> NonEmpty S.SubstitutionModel-getSubstitutionModels = map substModel . components+getSubstitutionModels :: MixtureModel -> N.NonEmpty S.SubstitutionModel+getSubstitutionModels = N.map substModel . components  -- | Create a mixture model from a list of substitution models. fromSubstitutionModels ::-  S.Name -> NonEmpty Weight -> NonEmpty S.SubstitutionModel -> MixtureModel+  S.Name -> N.NonEmpty Weight -> N.NonEmpty S.SubstitutionModel -> MixtureModel fromSubstitutionModels n ws sms =-  if allEqual $ toList alphs-    then MixtureModel n (head alphs) comps+  if allEqual $ N.toList alphs+    then MixtureModel n (N.head alphs) comps     else       error         "fromSubstitutionModels: alphabets of substitution models are not equal."   where-    comps = zipWith Component ws sms-    alphs = map S.alphabet sms+    comps = N.zipWith Component ws sms+    alphs = N.map S.alphabet sms+    allEqual [] = True+    allEqual xs = all (== head xs) $ tail xs  -- | Concatenate mixture models.-concatenate :: S.Name -> NonEmpty MixtureModel -> MixtureModel+concatenate :: S.Name -> N.NonEmpty MixtureModel -> MixtureModel concatenate n mms = fromSubstitutionModels n ws sms   where-    comps = sconcat $ map components mms-    ws = map weight comps-    sms = map substModel comps+    comps = sconcat $ N.map components mms+    ws = N.map weight comps+    sms = N.map substModel comps  scaleComponent :: Double -> Component -> Component scaleComponent s c = c {substModel = s'} where s' = S.scale s $ substModel c@@ -112,16 +96,16 @@ scale s m = m {components = cs'}   where     cs = components m-    cs' = map (scaleComponent s) cs+    cs' = N.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 $ zipWith (*) weights scales+    c = sum $ N.zipWith (*) weights scales     weights = getWeights mm-    scales = map S.totalRate $ getSubstitutionModels mm+    scales = N.map S.totalRate $ getSubstitutionModels mm  appendNameComponent :: S.Name -> Component -> Component appendNameComponent n c = c {substModel = s'}@@ -133,29 +117,4 @@ appendNameComponents n m = m {components = cs'}   where     cs = components m-    cs' = map (appendNameComponent n) cs---- | Summarize a mixture model component; lines to be printed to screen or log.-summarizeComponent :: Component -> [L.ByteString]-summarizeComponent c =-  L.pack "Weight: "-    <> (L.toLazyByteString . L.doubleDec $ weight c) :-  S.summarize (substModel c)---- | Summarize a mixture model; lines to be printed to screen or log.-summarize :: MixtureModel -> [L.ByteString]-summarize m =-  [ L.pack $ "Mixture model: " ++ name m ++ ".",-    L.pack $ "Number of components: " ++ show n ++ "."-  ]-    ++ detail-  where-    n = length $ components m-    detail =-      if n <= 100-        then-          concat-            [ L.pack ("Component " ++ show i ++ ":") : summarizeComponent c-              | (i, c) <- zip [1 :: Int ..] (toList $ components m)-            ]-        else []+    cs' = N.map (appendNameComponent n) cs
src/ELynx/Data/MarkovProcess/Nucleotide.hs view
@@ -25,10 +25,10 @@   ) where +import qualified Data.Vector.Storable as V import ELynx.Data.Alphabet.Alphabet import ELynx.Data.MarkovProcess.RateMatrix import ELynx.Data.MarkovProcess.SubstitutionModel-import ELynx.Tools import Numeric.LinearAlgebra hiding (normalize)  -- XXX: Another idea of structuring the code. This would probably be cleaner in@@ -77,9 +77,12 @@       0.0     ] +uniformVec :: Vector Double+uniformVec = V.replicate n (1 / fromIntegral n)+ -- | JC substitution model. jc :: SubstitutionModel-jc = substitutionModel DNA "JC" [] d jcExch where d = uniformVec n+jc = substitutionModel DNA "JC" [] uniformVec jcExch  -- | F81 substitution model. f81 :: StationaryDistribution -> SubstitutionModel
src/ELynx/Data/MarkovProcess/PhyloModel.hs view
@@ -18,11 +18,9 @@ module ELynx.Data.MarkovProcess.PhyloModel   ( PhyloModel (..),     getAlphabet,-    summarize,   ) where -import qualified Data.ByteString.Lazy.Char8 as L import ELynx.Data.Alphabet.Alphabet import qualified ELynx.Data.MarkovProcess.MixtureModel as M import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S@@ -34,10 +32,5 @@  -- | Extract code from phylogenetic model. getAlphabet :: PhyloModel -> Alphabet-getAlphabet (MixtureModel mm) = M.getAlphabet mm+getAlphabet (MixtureModel mm) = M.alphabet mm getAlphabet (SubstitutionModel sm) = S.alphabet sm---- | Summarize a phylogenetic model; lines to be printed to screen or log.-summarize :: PhyloModel -> [L.ByteString]-summarize (MixtureModel mm) = M.summarize mm-summarize (SubstitutionModel sm) = S.summarize sm
src/ELynx/Data/MarkovProcess/RateMatrix.hs view
@@ -34,7 +34,7 @@   ) where -import ELynx.Tools+import qualified Data.Vector.Storable as V import Numeric.LinearAlgebra hiding (normalize) import Numeric.SpecFunctions import Prelude hiding ((<>))@@ -50,14 +50,21 @@ -- | Stationary distribution of a rate matrix. type StationaryDistribution = Vector R +epsRelaxed :: Double+epsRelaxed = 1e-5+ -- | True if distribution sums to 1.0. isValid :: StationaryDistribution -> Bool-isValid d = nearlyEqWith eps' (norm_1 d) 1.0+isValid d = epsRelaxed > abs (norm_1 d - 1.0)  -- | Normalize a stationary distribution so that the elements sum to 1.0. normalizeSD :: StationaryDistribution -> StationaryDistribution normalizeSD d = d / scalar (norm_1 d) +matrixSetDiagToZero :: Matrix R -> Matrix R+matrixSetDiagToZero m = m - diag (takeDiag m)+{-# INLINE matrixSetDiagToZero #-}+ -- | Get average number of substitutions per unit time. totalRateWith :: StationaryDistribution -> RateMatrix -> Double totalRateWith d m = norm_1 $ d <# matrixSetDiagToZero m@@ -95,6 +102,15 @@   ExchangeabilityMatrix -> StationaryDistribution -> RateMatrix fromExchangeabilityMatrix em d = setDiagonal $ em <> diag d +eps :: Double+eps = 1e-12++normalizeSumVec :: V.Vector Double -> V.Vector Double+normalizeSumVec v = V.map (/ s) v+  where+    s = V.sum v+{-# INLINE normalizeSumVec #-}+ -- | Get stationary distribution from 'RateMatrix'. Involves eigendecomposition. -- If the given matrix does not satisfy the required properties of transition -- rate matrices and no eigenvector with an eigenvalue nearly equal to 0 is@@ -104,8 +120,8 @@ -- function)? getStationaryDistribution :: RateMatrix -> StationaryDistribution getStationaryDistribution m =-  if magnitude (eVals ! i) `nearlyEq` 0-    then normalizeSumVec 1.0 distReal+  if eps > abs (magnitude (eVals ! i))+    then normalizeSumVec distReal     else       error         "getStationaryDistribution: Could not retrieve stationary distribution."
src/ELynx/Data/MarkovProcess/SubstitutionModel.hs view
@@ -17,9 +17,10 @@     Params,     SubstitutionModel, -    -- * Lenses and other accessors+    -- * Accessors     alphabet,     name,+    params,     stationaryDistribution,     exchangeabilityMatrix,     rateMatrix,@@ -32,16 +33,12 @@     scale,     normalize,     appendName,--    -- * Output-    summarize,   ) where -import qualified Data.ByteString.Lazy.Char8 as L+import qualified Data.Vector.Storable as V import ELynx.Data.Alphabet.Alphabet import qualified ELynx.Data.MarkovProcess.RateMatrix as R-import ELynx.Tools import qualified Numeric.LinearAlgebra as LinAlg  -- | Name of substitution model; abstracted and subject to change.@@ -89,6 +86,12 @@ totalRate :: SubstitutionModel -> Double totalRate sm = R.totalRate (rateMatrix sm) +normalizeSumVec :: V.Vector Double -> V.Vector Double+normalizeSumVec v = V.map (/ s) v+  where+    s = V.sum v+{-# INLINE normalizeSumVec #-}+ -- | Create normalized 'SubstitutionModel'. See 'normalize'. substitutionModel ::   Alphabet ->@@ -99,13 +102,11 @@   SubstitutionModel substitutionModel c n ps d e =   if R.isValid d-    then normalize $ SubstitutionModel c n ps d' e+    then normalize $ SubstitutionModel c n ps (normalizeSumVec d) e     else       error $         "substitionModel: Stationary distribution does not sum to 1.0: "           ++ show d-  where-    d' = normalizeSumVec 1.0 d  -- | Scale the rate of a substitution model by given factor. scale :: Double -> SubstitutionModel -> SubstitutionModel@@ -121,29 +122,3 @@ -- | Abbend to name. appendName :: Name -> SubstitutionModel -> SubstitutionModel appendName n sm = sm {name = n'} where n' = name sm <> n---- | Summarize a substitution model; lines to be printed to screen or log.-summarize :: SubstitutionModel -> [L.ByteString]-summarize sm =-  map L.pack $-    (show (alphabet sm) ++ " substitution model: " ++ name sm ++ ".") :-    ["Parameters: " ++ show (params sm) ++ "." | not (null (params sm))]-      ++ case alphabet sm of-        DNA ->-          [ "Stationary distribution: "-              ++ dispv precision (stationaryDistribution sm)-              ++ ".",-            "Exchangeability matrix:\n"-              ++ dispmi 2 precision (exchangeabilityMatrix sm)-              ++ ".",-            "Scale: " ++ show (roundN precision $ totalRate sm) ++ "."-          ]-        Protein ->-          [ "Stationary distribution: "-              ++ dispv precision (stationaryDistribution sm)-              ++ ".",-            "Scale: " ++ show (roundN precision $ totalRate sm) ++ "."-          ]-        _ ->-          error-            "Extended character sets are not supported with substitution models."
src/ELynx/Import/MarkovProcess/EDMModelPhylobayes.hs view
@@ -1,3 +1,5 @@+{-# LANGUAGE OverloadedStrings #-}+ -- | -- Module      :  ELynx.Import.MarkovProcess.EDMModelPhylobayes -- Description :  Import stationary distributions from Phylobayes format@@ -10,63 +12,53 @@ -- -- Creation date: Tue Jan 29 12:12:55 2019. module ELynx.Import.MarkovProcess.EDMModelPhylobayes-  ( Parser,-    EDMComponent,+  ( EDMComponent,     phylobayes,   ) where +import Control.Applicative import Control.Monad-import qualified Data.ByteString.Lazy.Char8 as L+import qualified Data.Attoparsec.ByteString as AS+import qualified Data.Attoparsec.ByteString.Char8 as AC+import qualified Data.ByteString.Lazy.Char8 as BL import qualified Data.Vector.Storable as V-import Data.Void import ELynx.Data.MarkovProcess.MixtureModel-  ( Weight,-  )-import ELynx.Tools-import Text.Megaparsec-import Text.Megaparsec.Byte-import Text.Megaparsec.Byte.Lexer hiding (space) --- | Shortcut.-type Parser = Parsec Void L.ByteString- -- | An empirical mixture model component has a weight and a stationary -- distribution. type EDMComponent = (Weight, V.Vector Double)  -- | Parse stationary distributions from Phylobayes format.-phylobayes :: Parser [EDMComponent]-phylobayes = do+phylobayes :: AS.Parser [EDMComponent]+phylobayes = (AS.<?> "phylobayes") $ do   n <- headerLine   k <- kComponentsLine-  cs <- count k $ dataLine n-  _ <- many newline *> eof <?> "phylobayes"+  cs <- AS.count k $ dataLine n+  _ <- AC.skipWhile AC.isSpace+  _ <- AS.endOfInput   return cs -horizontalSpace :: Parser ()-horizontalSpace = skipMany $ char (c2w ' ') <|> tab--headerLine :: Parser Int+headerLine :: AS.Parser Int headerLine = do-  n <- decimal-  _ <- horizontalSpace+  n <- AC.decimal+  _ <- AS.skipWhile AC.isHorizontalSpace   -- XXX: This should be more general, but then we also want to ensure that the   -- order of states is correct.   _ <--    chunk (L.pack "A C D E F G H I K L M N P Q R S T V W Y")-      <|> chunk (L.pack "A C G T")-  _ <- many newline <?> "headerLine"+    AS.string (BL.toStrict "A C D E F G H I K L M N P Q R S T V W Y")+      <|> AS.string (BL.toStrict "A C G T")+  _ <- AC.skipWhile AC.isSpace AS.<?> "headerLine"   return n -kComponentsLine :: Parser Int-kComponentsLine = decimal <* newline <?> "kComponentsLine"+kComponentsLine :: AS.Parser Int+kComponentsLine = AC.decimal <* AC.skipWhile AC.isSpace AS.<?> "kComponentsLine" -dataLine :: Int -> Parser EDMComponent-dataLine n = do-  weight <- float-  _ <- horizontalSpace-  vals <- float `sepEndBy1` horizontalSpace+dataLine :: Int -> AS.Parser EDMComponent+dataLine n = (AS.<?> "dataLine") $ do+  wght <- AC.double+  _ <- AS.skipWhile AC.isHorizontalSpace+  vals <- AC.double `AC.sepBy1` AS.skipWhile AC.isHorizontalSpace+  _ <- AC.skipWhile AC.isSpace   when (length vals /= n) (error "Did not find correct number of entries.")-  _ <- space <?> "dataLine"-  return (weight, V.fromList vals)+  return (wght, V.fromList vals)
src/ELynx/Import/MarkovProcess/SiteprofilesPhylobayes.hs view
@@ -21,24 +21,21 @@   ) where +import Control.Applicative import Control.Monad+import qualified Data.Attoparsec.ByteString as AS+import qualified Data.Attoparsec.ByteString.Char8 as AC import Data.Containers.ListUtils (nubInt) import qualified Data.Vector.Storable as V import ELynx.Import.MarkovProcess.EDMModelPhylobayes-  ( EDMComponent,-    Parser,-  )-import ELynx.Tools-import Text.Megaparsec-import Text.Megaparsec.Byte-import Text.Megaparsec.Byte.Lexer  -- | Parse stationary distributions from Phylobayes format.-siteprofiles :: Parser [EDMComponent]-siteprofiles = do+siteprofiles :: AS.Parser [EDMComponent]+siteprofiles = (AS.<?> "phylobayes siteprofiles") $ do   _ <- headerLines   cs <- many dataLine-  _ <- many newline *> eof <?> "phylobayes siteprofiles"+  _ <- AC.skipWhile AC.isSpace+  _ <- AS.endOfInput   let ls = map length cs       nLs = length $ nubInt ls   when@@ -46,28 +43,20 @@     (error "The site profiles have a different number of entries.")   return cs -horizontalSpace :: Parser ()-horizontalSpace = skipMany $ char (c2w ' ') <|> tab--line :: Parser ()-line = do-  _ <- many $ noneOf [c2w '\n']-  pure ()+line :: AS.Parser ()+line = AS.skipWhile (not . AC.isEndOfLine)  -- For now, just ignore the header.-headerLines :: Parser ()-headerLines = do-  _ <- line-  _ <- many newline <?> "headerLine"-  pure ()+headerLines :: AS.Parser ()+headerLines = line *> AC.skipWhile AC.isSpace AS.<?> "headerLine" -dataLine :: Parser EDMComponent-dataLine = do+dataLine :: AS.Parser EDMComponent+dataLine = (AS.<?> "dataLine") $ do   -- Ignore site number.-  _ <- decimal :: Parser Integer-  _ <- horizontalSpace+  _ <- AC.decimal :: AS.Parser Int+  _ <- AS.skipWhile AC.isHorizontalSpace   -- Also ignore additional white space on line.-  vals <- float `sepEndBy1` horizontalSpace-  _ <- many newline <?> "dataLine"+  vals <- AC.double `AS.sepBy1` AS.skipWhile AC.isHorizontalSpace+  _ <- AC.skipWhile AC.isSpace   -- Set the weight to 1.0 for all sites.   return (1.0, V.fromList vals)
src/ELynx/Simulate/MarkovProcessAlongTree.hs view
@@ -1,3 +1,5 @@+{-# LANGUAGE ScopedTypeVariables #-}+ -- | --   Description :  Work with transition probability matrices on rooted trees --   Copyright   :  (c) Dominik Schrempf 2017@@ -24,15 +26,18 @@   ) where +import Control.Concurrent+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 Data.Word (Word32) import ELynx.Data.MarkovProcess.RateMatrix import ELynx.Simulate.MarkovProcess-import ELynx.Tools import Numeric.LinearAlgebra-import System.Random.MWC (Gen, GenIO)+import System.Random.MWC import System.Random.MWC.Distributions (categorical)  toProbTree :: RateMatrix -> Tree Double -> Tree ProbMatrix@@ -159,8 +164,28 @@       concat         <$> sequence [simulateAndFlattenMixtureModel' is' cs t g | t <- f] --- | See 'simulateAndFlattenMixtureModel', parallel version;--- needs to be run in IO monad.+getChunks :: Int -> Int -> [Int]+getChunks c n = ns+  where+    n' = n `div` c+    r = n `mod` c+    ns = replicate r (n' + 1) ++ replicate (c - r) n'++splitGen :: PrimMonad m => Int -> Gen (PrimState m) -> m [Gen (PrimState m)]+splitGen n gen+  | n <= 0 = return []+  | otherwise = do+    seeds :: [V.Vector Word32] <- replicateM (n -1) $ uniformVector gen 256+    fmap (gen :) (mapM initialize seeds)++parComp :: Int -> (Int -> GenIO -> IO b) -> GenIO -> IO [b]+parComp num fun gen = do+  ncap <- getNumCapabilities+  let chunks = getChunks ncap num+  gs <- splitGen ncap gen+  mapConcurrently (uncurry fun) (zip chunks gs)++-- | See 'simulateAndFlattenMixtureModel', parallel version. simulateAndFlattenMixtureModelPar ::   Int ->   Vector R ->
test/ELynx/Data/MarkovProcess/AminoAcidSpec.hs view
@@ -13,6 +13,7 @@   ) where +import qualified Data.Vector.Storable as V import ELynx.Data.MarkovProcess.AminoAcid import qualified ELynx.Data.MarkovProcess.RateMatrix as R import qualified ELynx.Data.MarkovProcess.SubstitutionModel as S@@ -20,9 +21,15 @@ import Numeric.LinearAlgebra import Test.Hspec +normalizeSumVec :: Vector Double -> Vector Double+normalizeSumVec v = V.map (/ s) v+  where+    s = V.sum v+{-# INLINE normalizeSumVec #-}+ statDistLGPython :: R.StationaryDistribution statDistLGPython =-  normalizeSumVec 1.0 $+  normalizeSumVec $     fromList       [ 0.079066,         0.012937,