elynx-markov 0.7.1.0 → 0.7.2.0

raw patch · 12 files changed

+82/−45 lines, 12 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• ChangeLog.md +5/−0
• README.md +4/−4
• elynx-markov.cabal +1/−1
• src/ELynx/MarkovProcess/AminoAcid.hs +14/−14
• src/ELynx/MarkovProcess/CXXModels.hs +2/−0
• src/ELynx/MarkovProcess/GammaRateHeterogeneity.hs +3/−5
• src/ELynx/MarkovProcess/MixtureModel.hs +22/−4
• src/ELynx/MarkovProcess/Nucleotide.hs +8/−8
• src/ELynx/MarkovProcess/SubstitutionModel.hs +12/−2
• test/ELynx/MarkovProcess/AminoAcidSpec.hs +7/−4
• test/ELynx/MarkovProcess/NucleotideSpec.hs +1/−1
• test/ELynx/Simulate/MarkovProcessAlongTreeSpec.hs +3/−2

@@ -5,6 +5,11 @@ ## Unreleased changes  +## Version 0.7.2.0++-   `slynx`: Allow global normalization of mixture models.++ ## Version 0.7.1.0  -   Be less strict with quoted identifiers/names in phylogenetic trees.

@@ -2,7 +2,7 @@  # The ELynx Suite -Version: 0.7.0.0.+Version: 0.7.2.0. Reproducible evolution made easy.  <p align="center"><img src="https://travis-ci.org/dschrempf/elynx.svg?branch=master"/></p>@@ -73,13 +73,13 @@     # OR: stack exec slynx -- --help     # OR: slynx --help -    Up to date-    ELynx Suite version 0.7.0.0.+    ELynx Suite version 0.7.2.0.     Developed by Dominik Schrempf.-    Compiled on July 10, 2022, at 08:48 am, UTC.+    Compiled on May 3, 2023, at 13:18 pm, UTC.          Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME]                  [-f|--force] [--no-elynx-file] COMMAND+           Analyze, and simulate multi sequence alignments.          Available options:

@@ -1,6 +1,6 @@ cabal-version:      3.0 name:               elynx-markov-version:            0.7.1.0+version:            0.7.2.0 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>.

@@ -368,12 +368,12 @@ lgStatDist = pamlToAlphaVec lgStatDistPaml  -- | LG substitution model.-lg :: SubstitutionModel-lg = substitutionModel Protein "LG" [] lgStatDist lgExch+lg :: Normalize -> SubstitutionModel+lg nz = substitutionModel Protein "LG" [] nz lgStatDist lgExch  -- | LG substitution model with maybe a name and a custom stationary distribution.-lgCustom :: Maybe String -> StationaryDistribution -> SubstitutionModel-lgCustom mnm d = substitutionModel Protein nm [] d lgExch+lgCustom :: Maybe String -> Normalize -> StationaryDistribution -> SubstitutionModel+lgCustom mnm nz d = substitutionModel Protein nm [] nz d lgExch   where     nm = fromMaybe "LG-Custom" mnm @@ -608,12 +608,12 @@ wagStatDist = pamlToAlphaVec wagStatDistPaml  -- | LG substitution model.-wag :: SubstitutionModel-wag = substitutionModel Protein "WAG" [] wagStatDist wagExch+wag :: Normalize -> SubstitutionModel+wag nz = substitutionModel Protein "WAG" [] nz wagStatDist wagExch  -- | LG substitution model with maybe a name and a custom stationary distribution.-wagCustom :: Maybe String -> StationaryDistribution -> SubstitutionModel-wagCustom mnm d = substitutionModel Protein nm [] d wagExch+wagCustom :: Maybe String -> Normalize -> StationaryDistribution -> SubstitutionModel+wagCustom mnm nz d = substitutionModel Protein nm [] nz d wagExch   where     nm = fromMaybe "WAG-Custom" mnm @@ -631,17 +631,17 @@ uniformVec = V.replicate n (1 / fromIntegral n)  -- | Poisson substitution model.-poisson :: SubstitutionModel-poisson = substitutionModel Protein "Poisson" [] uniformVec poissonExch+poisson :: Normalize -> SubstitutionModel+poisson nz = substitutionModel Protein "Poisson" [] nz uniformVec poissonExch  -- | Poisson substitution model with maybe a name and a custom stationary distribution.-poissonCustom :: Maybe String -> StationaryDistribution -> SubstitutionModel-poissonCustom mnm d = substitutionModel Protein nm [] d poissonExch+poissonCustom :: Maybe String -> Normalize -> StationaryDistribution -> SubstitutionModel+poissonCustom mnm nz d = substitutionModel Protein nm [] nz d poissonExch   where     nm = fromMaybe "Poisson-Custom" mnm  -- | General time reversible (GTR) substitution model for amino acids.-gtr20 :: [Double] -> StationaryDistribution -> SubstitutionModel-gtr20 es d = substitutionModel Protein "GTR" es d e+gtr20 :: Normalize -> [Double] -> StationaryDistribution -> SubstitutionModel+gtr20 nz es d = substitutionModel Protein "GTR" es nz d e   where     e = exchFromListUpper n es

@@ -117,6 +117,7 @@ cxxSubstitutionModelFromStatDist nComps comp d =   poissonCustom     (Just name)+    S.DoNotNormalize     (normalizeSD d)   where     name = componentName nComps comp@@ -139,6 +140,7 @@ cxxFromStatDistsAndWeights ws ds =   M.fromSubstitutionModels     (cxxName n)+    S.DoNormalize     (V.fromList ws)     sms   where

@@ -58,8 +58,7 @@     ++ "shape parameter "     ++ show alpha -splitSubstitutionModel ::-  Int -> Double -> S.SubstitutionModel -> V.Vector S.SubstitutionModel+splitSubstitutionModel :: Int -> Double -> S.SubstitutionModel -> V.Vector S.SubstitutionModel splitSubstitutionModel n alpha sm = renamedSMs   where     means = getMeans n alpha@@ -67,9 +66,8 @@     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+expandSubstitutionModel :: Int -> Double -> S.SubstitutionModel -> M.MixtureModel+expandSubstitutionModel n alpha sm = M.fromSubstitutionModels name S.DoNotNormalize ws sms   where     name = S.name sm <> getName n alpha     ws = V.replicate n 1.0

@@ -64,15 +64,33 @@ getSubstitutionModels :: MixtureModel -> V.Vector S.SubstitutionModel getSubstitutionModels = V.map substModel . components +normalizeGlobally :: V.Vector Weight -> V.Vector S.SubstitutionModel -> V.Vector S.SubstitutionModel+normalizeGlobally ws ss = V.map (S.scale $ recip c) ss+  where+    cks = V.map S.totalRate ss+    cNoWeights = V.sum $ V.zipWith (*) ws cks+    c = cNoWeights / V.sum ws+ -- | Create a mixture model from a list of substitution models.-fromSubstitutionModels :: S.Name -> V.Vector Weight -> V.Vector S.SubstitutionModel -> MixtureModel-fromSubstitutionModels n ws sms+--+-- If 'S.Normalize' is 'S.DoNormalize', globally normalize the mixture model.+-- Global normalization has no effect if all components are already normalized.+fromSubstitutionModels ::+  S.Name ->+  S.Normalize ->+  V.Vector Weight ->+  V.Vector S.SubstitutionModel ->+  MixtureModel+fromSubstitutionModels n nz 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 = V.zipWith Component ws sms+    smsNormalized = case nz of+      S.DoNormalize -> normalizeGlobally ws sms+      S.DoNotNormalize -> sms+    comps = V.zipWith Component ws smsNormalized     alphs = V.map S.alphabet sms     allEqual xs       | V.null xs = True@@ -80,7 +98,7 @@  -- | Concatenate mixture models. concatenate :: S.Name -> V.Vector MixtureModel -> MixtureModel-concatenate n mms = fromSubstitutionModels n ws sms+concatenate n mms = fromSubstitutionModels n S.DoNotNormalize ws sms   where     comps = V.concatMap components mms     ws = V.map weight comps

@@ -81,12 +81,12 @@ uniformVec = V.replicate n (1 / fromIntegral n)  -- | JC substitution model.-jc :: SubstitutionModel-jc = substitutionModel DNA "JC" [] uniformVec jcExch+jc :: Normalize -> SubstitutionModel+jc nz = substitutionModel DNA "JC" [] nz uniformVec jcExch  -- | F81 substitution model.-f81 :: StationaryDistribution -> SubstitutionModel-f81 d = substitutionModel DNA "F81" [] d jcExch+f81 :: Normalize -> StationaryDistribution -> SubstitutionModel+f81 nz d = substitutionModel DNA "F81" [] nz d jcExch  hkyExch :: Double -> ExchangeabilityMatrix hkyExch k =@@ -94,11 +94,11 @@     [0.0, 1.0, k, 1.0, 1.0, 0.0, 1.0, k, k, 1.0, 0.0, 1.0, 1.0, k, 1.0, 0.0]  -- | HKY substitution model.-hky :: Double -> StationaryDistribution -> SubstitutionModel-hky k d = substitutionModel DNA "HKY" [k] d e where e = hkyExch k+hky :: Normalize -> Double -> StationaryDistribution -> SubstitutionModel+hky nz k d = substitutionModel DNA "HKY" [k] nz d e where e = hkyExch k  -- | HKY substitution model.-gtr4 :: [Double] -> StationaryDistribution -> SubstitutionModel-gtr4 es d = substitutionModel DNA "GTR" es d e+gtr4 :: Normalize -> [Double] -> StationaryDistribution -> SubstitutionModel+gtr4 nz es d = substitutionModel DNA "GTR" es nz d e   where     e = exchFromListUpper n es

@@ -27,6 +27,7 @@     totalRate,      -- * Building substitution models+    Normalize (..),     substitutionModel,      -- * Transformations@@ -92,17 +93,26 @@     s = V.sum v {-# INLINE normalizeSumVec #-} +-- | Should the substitution model be normalized?+data Normalize = DoNormalize | DoNotNormalize+  deriving (Show, Read, Eq)+ -- | Create normalized 'SubstitutionModel'. See 'normalize'. substitutionModel ::   Alphabet ->   Name ->   Params ->+  Normalize ->   R.StationaryDistribution ->   R.ExchangeabilityMatrix ->   SubstitutionModel-substitutionModel c n ps d e =+substitutionModel c n ps nz d e =   if R.isValid d-    then normalize $ SubstitutionModel c n ps (normalizeSumVec d) e+    then+      let sm = SubstitutionModel c n ps (normalizeSumVec d) e+       in case nz of+            DoNormalize -> normalize sm+            DoNotNormalize -> sm     else       error $         "substitionModel: Stationary distribution does not sum to 1.0: "

@@ -16,6 +16,7 @@ import qualified Data.Vector.Storable as V import ELynx.MarkovProcess.AminoAcid import qualified ELynx.MarkovProcess.RateMatrix as R+import ELynx.MarkovProcess.SubstitutionModel (Normalize (DoNormalize)) import qualified ELynx.MarkovProcess.SubstitutionModel as S import ELynx.Tools.Equality import Numeric.LinearAlgebra@@ -482,13 +483,13 @@ statDistUniform = vector $ replicate 20 0.05  statDistLG :: R.StationaryDistribution-statDistLG = S.stationaryDistribution lg+statDistLG = S.stationaryDistribution $ lg DoNormalize  exchLG :: R.ExchangeabilityMatrix-exchLG = S.exchangeabilityMatrix lg+exchLG = S.exchangeabilityMatrix $ lg DoNormalize  rmLG :: R.RateMatrix-rmLG = S.rateMatrix lg+rmLG = S.rateMatrix $ lg DoNormalize  spec :: Spec spec = do@@ -513,11 +514,12 @@               S.rateMatrix $                 lgCustom                   Nothing+                  DoNormalize                   statDistUniform       f `nearlyEqVec` statDistUniform `shouldBe` True   describe "poisson" $     it "stationary distribution is uniform 1/20" $-      R.getStationaryDistribution (S.rateMatrix poisson)+      R.getStationaryDistribution (S.rateMatrix $ poisson DoNormalize)         `nearlyEqVec` statDistUniform         `shouldBe` True   describe "poissonCustom" $@@ -527,5 +529,6 @@               S.rateMatrix $                 poissonCustom                   Nothing+                  DoNormalize                   statDistLGPython       f `nearlyEqVec` statDistLGPython `shouldBe` True

@@ -24,7 +24,7 @@ stationaryDist = fromList [0.2, 0.3, 0.3, 0.2]  hkyModel :: SubstitutionModel-hkyModel = hky 6.0 stationaryDist+hkyModel = hky DoNormalize 6.0 stationaryDist  spec :: Spec spec =

@@ -15,6 +15,7 @@  import Data.Tree import ELynx.MarkovProcess.Nucleotide+import ELynx.MarkovProcess.SubstitutionModel (Normalize (DoNormalize)) import qualified ELynx.MarkovProcess.SubstitutionModel as S import ELynx.Simulate.MarkovProcess import ELynx.Simulate.MarkovProcessAlongTree@@ -50,5 +51,5 @@     tr <- simulate 10 d e testTree gen     (length . rootLabel $ tr) `shouldBe` 10   where-    d = S.stationaryDistribution jc-    e = S.exchangeabilityMatrix jc+    d = S.stationaryDistribution $ jc DoNormalize+    e = S.exchangeabilityMatrix $ jc DoNormalize