packages feed

elynx-tree 0.6.1.1 → 0.7.0.1

raw patch · 51 files changed

+469/−341 lines, 51 filesdep +randomdep −mwc-randomdep −primitivePVP ok

version bump matches the API change (PVP)

Dependencies added: random

Dependencies removed: mwc-random, primitive

API changes (from Hackage documentation)

+ ELynx.Tree.Mrca: findNode :: (Ord a, Show a) => a -> Tree e a -> Either String (Set a)
+ ELynx.Tree.Mrca: getPathToMrca :: (Ord a, Show a) => Set a -> Tree e a -> Either String Path
+ ELynx.Tree.Mrca: getTreeAtMrca :: (Ord a, Show a) => Set a -> Tree e a -> Either String (Tree e a)
+ ELynx.Tree.Mrca: isAncestor :: Ord a => Set a -> Tree e a -> Bool
+ ELynx.Tree.Mrca: isMrca :: Ord a => Set a -> Tree e a -> Bool
+ ELynx.Tree.Rooted: duplicateLabels :: Ord a => Tree e a -> Bool
- ELynx.Tree.Simulate.Coalescent: simulate :: PrimMonad m => Int -> Gen (PrimState m) -> m (Tree Length Int)
+ ELynx.Tree.Simulate.Coalescent: simulate :: StatefulGen g m => Int -> g -> m (Tree Length Int)
- ELynx.Tree.Simulate.PointProcess: simulate :: PrimMonad m => Int -> TimeSpec -> Rate -> Rate -> Gen (PrimState m) -> m (PointProcess Int Double)
+ ELynx.Tree.Simulate.PointProcess: simulate :: StatefulGen g m => Int -> TimeSpec -> Rate -> Rate -> g -> m (PointProcess Int Double)
- ELynx.Tree.Simulate.PointProcess: simulateNReconstructedTrees :: PrimMonad m => Int -> Int -> TimeSpec -> Rate -> Rate -> Gen (PrimState m) -> m (Forest Length Int)
+ ELynx.Tree.Simulate.PointProcess: simulateNReconstructedTrees :: StatefulGen g m => Int -> Int -> TimeSpec -> Rate -> Rate -> g -> m (Forest Length Int)
- ELynx.Tree.Simulate.PointProcess: simulateReconstructedTree :: PrimMonad m => Int -> TimeSpec -> Rate -> Rate -> Gen (PrimState m) -> m (Tree Length Int)
+ ELynx.Tree.Simulate.PointProcess: simulateReconstructedTree :: StatefulGen g m => Int -> TimeSpec -> Rate -> Rate -> g -> m (Tree Length Int)

Files

ChangeLog.md view
@@ -5,9 +5,16 @@ ## Unreleased changes  +## Version 0.7.0.1++-   Random 1.2: Parallel functions now require an \`IOGenM\` random number+    generator.+-   Fix splitting of the random number generator.++ ## Version 0.6.1.1 --   Remove plotting functionality (gnuplot incompatible with ghc921).+-   Remove plotting functionality (gnuplot incompatible with ghc922). -   Read files strictly. -   Refactor; flatten model hierarchy.
README.md view
@@ -2,7 +2,7 @@  # The ELynx Suite -Version: 0.6.1.1.+Version: 0.7.0.0. Reproducible evolution made easy.  <p align="center"><img src="https://travis-ci.org/dschrempf/elynx.svg?branch=master"/></p>@@ -74,9 +74,9 @@     # OR: slynx --help      Up to date-    ELynx Suite version 0.6.1.1.+    ELynx Suite version 0.7.0.0.     Developed by Dominik Schrempf.-    Compiled on February 22, 2022, at 15:10 pm, UTC.+    Compiled on July 10, 2022, at 08:48 am, UTC.          Usage: slynx [-v|--verbosity VALUE] [-o|--output-file-basename NAME]                  [-f|--force] [--no-elynx-file] COMMAND@@ -144,16 +144,16 @@     # OR: slynx simulate --help      Up to date-    ELynx Suite version 0.6.1.1.+    ELynx Suite version 0.7.0.0.     Developed by Dominik Schrempf.-    Compiled on February 22, 2022, at 15:10 pm, UTC.+    Compiled on July 10, 2022, at 08:48 am, UTC.          Usage: slynx simulate (-t|--tree-file Name) [-s|--substitution-model MODEL]                           [-m|--mixture-model MODEL] [-e|--edm-file NAME]                           [-p|--siteprofile-files NAMES]                           [-w|--mixture-model-weights "[DOUBLE,DOUBLE,...]"]                           [-g|--gamma-rate-heterogeneity "(NCAT,SHAPE)"]-                          (-l|--length NUMBER) [-S|--seed [INT]]+                          (-l|--length NUMBER) [-S|--seed INT]       Simulate multi sequence alignments.          Available options:@@ -175,8 +175,7 @@       -g,--gamma-rate-heterogeneity "(NCAT,SHAPE)"                                Number of gamma rate categories and shape parameter       -l,--length NUMBER       Set alignment length to NUMBER-      -S,--seed [INT]          Seed for random number generator; list of 32 bit-                               integers with up to 256 elements (default: random)+      -S,--seed INT            Seed for random number generator (default: random)       -h,--help                Show this help text          Substitution models:
bench/Bench.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  Bench -- Description :  Various benchmarks--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -22,7 +22,7 @@ import ELynx.Tree.Simulate.PointProcess import Length import Lens-import System.Random.MWC+import System.Random.Stateful import Tree  treeFileMany :: FilePath@@ -32,7 +32,7 @@ getManyTrees = parseFileWith (someNewick Standard) treeFileMany  hugeTree :: IO (Tree Length Int)-hugeTree = create >>= simulateReconstructedTree 50000 Random 1.0 0.9+hugeTree = newIOGenM (mkStdGen 0) >>= simulateReconstructedTree 50000 Random 1.0 0.9  sinN :: Int -> Length -> Length sinN n x = iterate sin x !! n
bench/Length.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  Length -- Description :  Benchmark length newtype--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
bench/Lens.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  Lens -- Description :  Benchmark lens operations--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
bench/Tree.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  Tree -- Description :  Benchmark some tree algorithms--- Copyright   :  (c) 2021 Dominik Schrempf+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
elynx-tree.cabal view
@@ -1,81 +1,82 @@-cabal-version:  3.0-name:           elynx-tree-version:        0.6.1.1-synopsis:       Handle phylogenetic trees-description:    Examine, compare, and simulate phylogenetic trees in a reproducible way. Please see the README on GitHub at <https://github.com/dschrempf/elynx>.-category:       Bioinformatics-homepage:       https://github.com/dschrempf/elynx#readme-bug-reports:    https://github.com/dschrempf/elynx/issues-author:         Dominik Schrempf-maintainer:     dominik.schrempf@gmail.com-copyright:      Dominik Schrempf (2021)-license:        GPL-3.0-or-later-license-file:   LICENSE-build-type:     Simple+cabal-version:      3.0+name:               elynx-tree+version:            0.7.0.1+synopsis:           Handle phylogenetic trees+description:+  Examine, compare, and simulate phylogenetic trees in a reproducible way. Please see the README on GitHub at <https://github.com/dschrempf/elynx>. +category:           Bioinformatics+homepage:           https://github.com/dschrempf/elynx#readme+bug-reports:        https://github.com/dschrempf/elynx/issues+author:             Dominik Schrempf+maintainer:         dominik.schrempf@gmail.com+copyright:          2021 Dominik Schrempf+license:            GPL-3.0-or-later+license-file:       LICENSE+build-type:         Simple extra-source-files:-    README.md-    ChangeLog.md-    data/ConnectA.tree-    data/ConnectB.tree-    data/ConnectConstraints.tree-    data/Many.trees-    data/MicrospoEDM32.tree-    data/MicrospoEDM64.tree-    data/MicrospoPoisson.tree-    data/Multifurcating.tree-    data/NewickConstraint.tree-    data/NewickRevBayes.tree-    data/NewickRooted.tree-    data/Newick.tree-    data/RevBayes.trees-    data/SimpleTree.nex-    data/TreeDist.trees-    data/UltraMetric.tree+  ChangeLog.md+  data/ConnectA.tree+  data/ConnectB.tree+  data/ConnectConstraints.tree+  data/Many.trees+  data/MicrospoEDM32.tree+  data/MicrospoEDM64.tree+  data/MicrospoPoisson.tree+  data/Multifurcating.tree+  data/Newick.tree+  data/NewickConstraint.tree+  data/NewickRevBayes.tree+  data/NewickRooted.tree+  data/RevBayes.trees+  data/SimpleTree.nex+  data/TreeDist.trees+  data/UltraMetric.tree+  README.md  source-repository head-  type: git+  type:     git   location: https://github.com/dschrempf/elynx  library   exposed-modules:-      ELynx.Topology-      ELynx.Topology.Phylogeny-      ELynx.Topology.Rooted-      ELynx.Tree-      ELynx.Tree.Bipartition-      ELynx.Tree.Distance-      ELynx.Tree.Length-      ELynx.Tree.Name-      ELynx.Tree.Parallel-      ELynx.Tree.Partition-      ELynx.Tree.Phylogeny-      ELynx.Tree.Rooted-      ELynx.Tree.Splittable-      ELynx.Tree.Support-      ELynx.Tree.Zipper-      ELynx.Tree.Distribution.BirthDeath-      ELynx.Tree.Distribution.BirthDeathCritical-      ELynx.Tree.Distribution.BirthDeathCriticalNoTime-      ELynx.Tree.Distribution.BirthDeathNearlyCritical-      ELynx.Tree.Distribution.CoalescentContinuous-      ELynx.Tree.Distribution.TimeOfOrigin-      ELynx.Tree.Distribution.TimeOfOriginNearCritical-      ELynx.Tree.Distribution.Types-      ELynx.Tree.Export.Newick-      ELynx.Tree.Export.Nexus-      ELynx.Tree.Import.Newick-      ELynx.Tree.Import.Nexus-      ELynx.Tree.Simulate.Coalescent-      ELynx.Tree.Simulate.PointProcess-  other-modules:-      Paths_elynx_tree-  autogen-modules:-      Paths_elynx_tree-  hs-source-dirs: src-  ghc-options: -Wall -Wunused-packages+    ELynx.Topology+    ELynx.Topology.Phylogeny+    ELynx.Topology.Rooted+    ELynx.Tree+    ELynx.Tree.Bipartition+    ELynx.Tree.Distance+    ELynx.Tree.Distribution.BirthDeath+    ELynx.Tree.Distribution.BirthDeathCritical+    ELynx.Tree.Distribution.BirthDeathCriticalNoTime+    ELynx.Tree.Distribution.BirthDeathNearlyCritical+    ELynx.Tree.Distribution.CoalescentContinuous+    ELynx.Tree.Distribution.TimeOfOrigin+    ELynx.Tree.Distribution.TimeOfOriginNearCritical+    ELynx.Tree.Distribution.Types+    ELynx.Tree.Export.Newick+    ELynx.Tree.Export.Nexus+    ELynx.Tree.Import.Newick+    ELynx.Tree.Import.Nexus+    ELynx.Tree.Length+    ELynx.Tree.Mrca+    ELynx.Tree.Name+    ELynx.Tree.Parallel+    ELynx.Tree.Partition+    ELynx.Tree.Phylogeny+    ELynx.Tree.Rooted+    ELynx.Tree.Simulate.Coalescent+    ELynx.Tree.Simulate.PointProcess+    ELynx.Tree.Splittable+    ELynx.Tree.Support+    ELynx.Tree.Zipper++  other-modules:    Paths_elynx_tree+  autogen-modules:  Paths_elynx_tree+  hs-source-dirs:   src+  ghc-options:      -Wall -Wunused-packages   build-depends:-      aeson+    , aeson     , attoparsec     , base     , bytestring@@ -86,37 +87,39 @@     , double-conversion     , elynx-nexus     , math-functions-    , mwc-random     , parallel-    , primitive+    , random     , statistics+   default-language: Haskell2010  test-suite tree-test-  type: exitcode-stdio-1.0-  main-is: Spec.hs+  type:             exitcode-stdio-1.0+  main-is:          Spec.hs   other-modules:-      ELynx.ClassLaws-      ELynx.Topology.Arbitrary-      ELynx.Topology.RootedSpec-      ELynx.Tree.Arbitrary-      ELynx.Tree.BipartitionSpec-      ELynx.Tree.DistanceSpec-      ELynx.Tree.PartitionSpec-      ELynx.Tree.PhylogenySpec-      ELynx.Tree.RootedSpec-      ELynx.Tree.SupportSpec-      ELynx.Tree.Export.NewickSpec-      ELynx.Tree.Export.NexusSpec-      ELynx.Tree.Import.NewickSpec-      ELynx.Tree.Import.NexusSpec-      Paths_elynx_tree-  autogen-modules:-      Paths_elynx_tree-  hs-source-dirs: test-  ghc-options: -Wall -Wunused-packages -threaded -rtsopts -with-rtsopts=-N+    ELynx.ClassLaws+    ELynx.Topology.Arbitrary+    ELynx.Topology.RootedSpec+    ELynx.Tree.Arbitrary+    ELynx.Tree.BipartitionSpec+    ELynx.Tree.DistanceSpec+    ELynx.Tree.Export.NewickSpec+    ELynx.Tree.Export.NexusSpec+    ELynx.Tree.Import.NewickSpec+    ELynx.Tree.Import.NexusSpec+    ELynx.Tree.PartitionSpec+    ELynx.Tree.PhylogenySpec+    ELynx.Tree.RootedSpec+    ELynx.Tree.SupportSpec+    Paths_elynx_tree++  autogen-modules:  Paths_elynx_tree+  hs-source-dirs:   test+  ghc-options:+    -Wall -Wunused-packages -threaded -rtsopts -with-rtsopts=-N+   build-depends:-      attoparsec+    , attoparsec     , base     , bytestring     , containers@@ -126,26 +129,30 @@     , hspec     , QuickCheck     , quickcheck-classes+   default-language: Haskell2010  benchmark tree-bench-  type: exitcode-stdio-1.0-  main-is: Bench.hs+  type:             exitcode-stdio-1.0+  main-is:          Bench.hs   other-modules:-      Length-      Lens-      Tree-      Paths_elynx_tree-  autogen-modules:-      Paths_elynx_tree-  hs-source-dirs: bench-  ghc-options: -Wall -Wunused-packages -threaded -rtsopts -with-rtsopts=-N+    Length+    Lens+    Paths_elynx_tree+    Tree++  autogen-modules:  Paths_elynx_tree+  hs-source-dirs:   bench+  ghc-options:+    -Wall -Wunused-packages -threaded -rtsopts -with-rtsopts=-N+   build-depends:-      base >=4.7 && <5+    , base         >=4.7 && <5     , criterion     , elynx-tools     , elynx-tree     , microlens-    , mwc-random     , parallel+    , random+   default-language: Haskell2010
src/ELynx/Topology.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Topology -- Description :  Rooted topologies--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
src/ELynx/Topology/Phylogeny.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Topology.Phylogeny -- Description :  Phylogenetic topologies--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -66,8 +66,8 @@ intersect ts   | S.null lvsCommon = Left "intersect: Intersection of leaves is empty."   | otherwise = case sequence [dropLeavesWith (predicate ls) t | (ls, t) <- zip leavesToDrop ts] of-    Nothing -> Left "intersect: A topology is empty."-    Just ts' -> Right ts'+      Nothing -> Left "intersect: A topology is empty."+      Just ts' -> Right ts'   where     -- Leaf sets.     lvss = map (S.fromList . leaves) ts
src/ELynx/Topology/Rooted.hs view
@@ -4,7 +4,7 @@ -- | -- Module      :  ELynx.Topology.Rooted -- Description :  Topologies--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
src/ELynx/Tree.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Tree -- Description :  Rooted trees--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -32,6 +32,9 @@     -- * Phylogenies     module ELynx.Tree.Phylogeny, +    -- * Most recent common ancestors+    module ELynx.Tree.Mrca,+     -- * Partitions and distances     module ELynx.Tree.Bipartition,     module ELynx.Tree.Partition,@@ -52,6 +55,7 @@ import ELynx.Tree.Import.Newick import ELynx.Tree.Import.Nexus import ELynx.Tree.Length+import ELynx.Tree.Mrca import ELynx.Tree.Name import ELynx.Tree.Parallel import ELynx.Tree.Partition
src/ELynx/Tree/Bipartition.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Tree.Bipartition -- Description :  Bipartitions on trees--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -146,7 +146,7 @@   [Set a] getComplementaryLeaves p (Node _ _ ts) =   [ S.unions $ p : take i lvsChildren ++ drop (i + 1) lvsChildren-    | i <- [0 .. (n -1)]+    | i <- [0 .. (n - 1)]   ]   where     n = length ts@@ -158,8 +158,8 @@ bipartitions' p (Node _ p' []) = either (const S.empty) S.singleton $ bp p p' bipartitions' p t@(Node _ p' ts) =   S.unions $-    either (const S.empty) S.singleton (bp p p') :-      [bipartitions' c s | (c, s) <- zip cs ts]+    either (const S.empty) S.singleton (bp p p')+      : [bipartitions' c s | (c, s) <- zip cs ts]   where     cs = getComplementaryLeaves p t @@ -200,7 +200,7 @@   Map (Bipartition a) e bipartitionToBranch' p t@(Node b p' ts) =   M.unionsWith (<>) $-    either (const M.empty) (`M.singleton` b) (bp p p') :-      [bipartitionToBranch' c s | (c, s) <- zip cs ts]+    either (const M.empty) (`M.singleton` b) (bp p p')+      : [bipartitionToBranch' c s | (c, s) <- zip cs ts]   where     cs = getComplementaryLeaves p t
src/ELynx/Tree/Distance.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Tree.Distance -- Description :  Compute distances between trees--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -55,11 +55,11 @@ symmetric :: Ord a => Tree e1 a -> Tree e2 a -> Either String Int symmetric t1 t2   | S.fromList (leaves t1) /= S.fromList (leaves t2) =-    Left "symmetric: Trees contain different leaves."+      Left "symmetric: Trees contain different leaves."   | otherwise = do-    bps1 <- bipartitions t1-    bps2 <- bipartitions t2-    return $ length $ symmetricDifference bps1 bps2+      bps1 <- bipartitions t1+      bps2 <- bipartitions t2+      return $ length $ symmetricDifference bps1 bps2  countIncompatibilities :: (Show a, Ord a) => Set (Bipartition a) -> Set (Partition a) -> Int countIncompatibilities bs ms =@@ -102,22 +102,22 @@ incompatibleSplits :: (Show a, Ord a) => Tree e1 a -> Tree e2 a -> Either String Int incompatibleSplits t1 t2   | S.fromList (leaves t1) /= S.fromList (leaves t2) =-    Left "incompatibleSplits: Trees contain different leaves."+      Left "incompatibleSplits: Trees contain different leaves."   | otherwise = do-    -- Bipartitions.-    bs1 <- bipartitions t1-    bs2 <- bipartitions t2-    -- traceShowM $ "bs1" ++ show (S.map bpHuman bs1)-    -- traceShowM $ "bs2" ++ show (S.map bpHuman bs2)-    let -- Putative incompatible bipartitions of trees one and two, respectively.-        putIncBs1 = bs1 S.\\ bs2-        putIncBs2 = bs2 S.\\ bs1-    -- Partitions.-    ms1 <- partitions t1-    ms2 <- partitions t2-    -- traceShowM $ "putIncBs1 " ++ show (S.map bpHuman putIncBs1)-    -- traceShowM $ "putIncBs2 " ++ show (S.map bpHuman putIncBs2)-    return $ countIncompatibilities putIncBs1 ms2 + countIncompatibilities putIncBs2 ms1+      -- Bipartitions.+      bs1 <- bipartitions t1+      bs2 <- bipartitions t2+      -- traceShowM $ "bs1" ++ show (S.map bpHuman bs1)+      -- traceShowM $ "bs2" ++ show (S.map bpHuman bs2)+      let -- Putative incompatible bipartitions of trees one and two, respectively.+          putIncBs1 = bs1 S.\\ bs2+          putIncBs2 = bs2 S.\\ bs1+      -- Partitions.+      ms1 <- partitions t1+      ms2 <- partitions t2+      -- traceShowM $ "putIncBs1 " ++ show (S.map bpHuman putIncBs1)+      -- traceShowM $ "putIncBs2 " ++ show (S.map bpHuman putIncBs2)+      return $ countIncompatibilities putIncBs1 ms2 + countIncompatibilities putIncBs2 ms1  -- | Compute branch score distance between two trees. --@@ -128,10 +128,10 @@ branchScore :: (HasLength e1, HasLength e2, Ord a) => Tree e1 a -> Tree e2 a -> Either String Double branchScore t1 t2   | S.fromList (leaves t1) /= S.fromList (leaves t2) =-    Left "branchScoreWith: Trees contain different leaves."+      Left "branchScoreWith: Trees contain different leaves."   | otherwise = do-    bpToBr1 <- bipartitionToBranch $ first (Sum . getLength) t1-    bpToBr2 <- bipartitionToBranch $ first (Sum . getLength) t2-    let dBs = M.unionWith (-) bpToBr1 bpToBr2-        dsSquared = foldl' (\acc e -> acc + e * e) 0 dBs-    return $ sqrt $ fromLength $ getSum dsSquared+      bpToBr1 <- bipartitionToBranch $ first (Sum . getLength) t1+      bpToBr2 <- bipartitionToBranch $ first (Sum . getLength) t2+      let dBs = M.unionWith (-) bpToBr1 bpToBr2+          dsSquared = foldl' (\acc e -> acc + e * e) 0 dBs+      return $ sqrt $ fromLength $ getSum dsSquared
src/ELynx/Tree/Distribution/BirthDeath.hs view
@@ -4,7 +4,7 @@ -- | --   Module      :  ELynx.Tree.Distribution.BirthDeath --   Description :  Birth and death distribution---   Copyright   :  (c) Dominik Schrempf 2021+--   Copyright   :  2021 Dominik Schrempf --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -57,8 +57,8 @@   | otherwise = t1 * t2   where     d = l - m-    t1 = (1.0 - exp (- d * x)) / (l - m * exp (- d * x))-    t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))+    t1 = (1.0 - exp (-d * x)) / (l - m * exp (-d * x))+    t2 = (l - m * exp (-d * t)) / (1.0 - exp (-d * t))  instance D.ContDistr BirthDeathDistribution where   density = density@@ -72,23 +72,23 @@   | otherwise = d ** 2 * t1 * t2   where     d = l - m-    t1 = exp (- d * x) / ((l - m * exp (- d * x)) ** 2)-    t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))+    t1 = exp (-d * x) / ((l - m * exp (-d * x)) ** 2)+    t2 = (l - m * exp (-d * t)) / (1.0 - exp (-d * t))  -- | Inverted cumulative probability distribution 'cumulative'. See also -- 'D.ContDistr'. quantile :: BirthDeathDistribution -> Double -> Time quantile (BDD t l m) p   | p >= 0 && p <= 1 =-    res+      res   | otherwise =-    error $-      "BirthDeath.quantile: p must be in range [0,1] but got "-        ++ show p-        ++ "."+      error $+        "BirthDeath.quantile: p must be in range [0,1] but got "+          ++ show p+          ++ "."   where     d = l - m-    t2 = (l - m * exp (- d * t)) / (1.0 - exp (- d * t))+    t2 = (l - m * exp (-d * t)) / (1.0 - exp (-d * t))     res = (-1.0 / d) * log ((1.0 - p * l / t2) / (1.0 - p * m / t2))  instance D.ContGen BirthDeathDistribution where
src/ELynx/Tree/Distribution/BirthDeathCritical.hs view
@@ -4,7 +4,7 @@ -- | --   Module      :  ELynx.Tree.Distribution.BirthDeathCritical --   Description :  Birth and death distribution---   Copyright   :  (c) Dominik Schrempf 2021+--   Copyright   :  2021 Dominik Schrempf --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -71,12 +71,12 @@ quantile :: BirthDeathCriticalDistribution -> Double -> Time quantile (BDCD t l) p   | p >= 0 && p <= 1 =-    res+      res   | otherwise =-    error $-      "PointProcess.quantile: p must be in [0,1] range. Got: "-        ++ show p-        ++ "."+      error $+        "PointProcess.quantile: p must be in [0,1] range. Got: "+          ++ show p+          ++ "."   where     res = p * t / (1 + l * t - l * p * t)
src/ELynx/Tree/Distribution/BirthDeathCriticalNoTime.hs view
@@ -4,7 +4,7 @@ -- | --   Module      :  ELynx.Tree.Distribution.BirthDeathCriticalNoTime --   Description :  Birth and death distribution---   Copyright   :  (c) Dominik Schrempf 2021+--   Copyright   :  2021 Dominik Schrempf --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -68,9 +68,9 @@ quantile :: BirthDeathCriticalNoTimeDistribution -> Double -> Time quantile (BDCNTD l) p   | p >= 0 && p <= 1 =-    p / (l - l * p)+      p / (l - l * p)   | otherwise =-    error $ "PointProcess.quantile: p must be in [0,1]. Got: " ++ show p ++ "."+      error $ "PointProcess.quantile: p must be in [0,1]. Got: " ++ show p ++ "."  instance D.ContGen BirthDeathCriticalNoTimeDistribution where   genContVar = D.genContinuous
src/ELynx/Tree/Distribution/BirthDeathNearlyCritical.hs view
@@ -4,7 +4,7 @@ -- | --   Module      :  ELynx.Tree.Distribution.BirthDeathNearlyCritical --   Description :  Birth and death distribution---   Copyright   :  (c) Dominik Schrempf 2021+--   Copyright   :  2021 Dominik Schrempf --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -60,7 +60,7 @@   | otherwise = o0 + o1   where     o0 = s * (1.0 + t * l) / t / (1.0 + s * l)-    o1 = (- s * s + s * t) * (m - l) / (2.0 * t * (1.0 + s * l) ** 2)+    o1 = (-s * s + s * t) * (m - l) / (2.0 * t * (1.0 + s * l) ** 2)  instance D.ContDistr BirthDeathNearlyCriticalDistribution where   density = density@@ -81,27 +81,27 @@ quantile :: BirthDeathNearlyCriticalDistribution -> Double -> Time quantile (BDNCD t l m) p   | p >= 0 && p <= 1 =-    res+      res   | otherwise =-    error $-      "PointProcess.quantile: p must be in [0,1] range. Got: "-        ++ show p-        ++ "."+      error $+        "PointProcess.quantile: p must be in [0,1] range. Got: "+          ++ show p+          ++ "."   where     den = l * (-3.0 + 2.0 * t * (-1.0 + p) * l) + m     t1 = (2.0 + t * (l - 4.0 * p * l + m)) / den     t2Nom =       4.0         + t-        * ( l-              * (4.0 + t * l + 8.0 * p * (1.0 + t * l))-              + 2.0-              * (2.0 + t * l - 4.0 * p * (1.0 + t * l))-              * m-              + t-              * m-              * m-          )+          * ( l+                * (4.0 + t * l + 8.0 * p * (1.0 + t * l))+                + 2.0+                  * (2.0 + t * l - 4.0 * p * (1.0 + t * l))+                  * m+                + t+                  * m+                  * m+            )     t2 = t2Nom / (den ** 2)     res = 0.5 * (t1 + sqrt t2)
src/ELynx/Tree/Distribution/CoalescentContinuous.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Tree.Distribution.CoalescentContinuous -- Description :  Distribution of coalescent times--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
src/ELynx/Tree/Distribution/TimeOfOrigin.hs view
@@ -4,7 +4,7 @@ -- | --   Module      :  ELynx.Tree.Distribution.TimeOfOrigin --   Description :  Distribution of time of origin for birth and death trees---   Copyright   :  (c) Dominik Schrempf 2021+--   Copyright   :  2021 Dominik Schrempf --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -56,7 +56,7 @@   | otherwise = te ** fromIntegral n   where     d = l - m-    te = l * (1.0 - exp (- d * x)) / (l - m * exp (- d * x))+    te = l * (1.0 - exp (-d * x)) / (l - m * exp (-d * x))  instance D.ContDistr TimeOfOriginDistribution where   density = density@@ -70,7 +70,7 @@   where     d = l - m     n = fromIntegral nn-    ex = exp (- d * x)+    ex = exp (-d * x)     t1 = 1.0 - ex     t2 = l - m * ex @@ -79,12 +79,12 @@ quantile :: TimeOfOriginDistribution -> Double -> Time quantile (TOD n' l m) p   | p >= 0 && p <= 1 =-    -1.0 / d * log (t1 / t2)+      -1.0 / d * log (t1 / t2)   | otherwise =-    error $-      "PointProcess.quantile: p must be in [0,1] range. Got: "-        ++ show p-        ++ "."+      error $+        "PointProcess.quantile: p must be in [0,1] range. Got: "+          ++ show p+          ++ "."   where     d = l - m     n = fromIntegral n'
src/ELynx/Tree/Distribution/TimeOfOriginNearCritical.hs view
@@ -4,7 +4,7 @@ -- | --   Module      :  ELynx.Tree.Distribution.TimeOfOriginNearCritical --   Description :  Distribution of time of origin for birth and death trees---   Copyright   :  (c) Dominik Schrempf 2021+--   Copyright   :  2021 Dominik Schrempf --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -80,15 +80,15 @@ quantile :: TimeOfOriginNearCriticalDistribution -> Double -> Time quantile (TONCD n' l m) p   | p >= 0 && p <= 1 =-    t1 + t2nom / t2den+      t1 + t2nom / t2den   | otherwise =-    error $-      "PointProcess.quantile: p must be in [0,1] range. Got: "-        ++ show p-        ++ "."+      error $+        "PointProcess.quantile: p must be in [0,1] range. Got: "+          ++ show p+          ++ "."   where     n = fromIntegral n'-    t1 = - p ** (1 / n) / ((-1 + p ** (1 / n)) * l)+    t1 = -p ** (1 / n) / ((-1 + p ** (1 / n)) * l)     t2nom = p ** (2 / n) * (m - l)     t2den = 2 * (-1 + p ** (1 / n)) ** 2 * l ** 2
src/ELynx/Tree/Distribution/Types.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Tree.Distribution.Types -- Description :  Data types for distributions on trees--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
src/ELynx/Tree/Export/Newick.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Tree.Export.Newick -- Description :  Export tree objects to Newick format--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
src/ELynx/Tree/Export/Nexus.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Tree.Export.Nexus -- Description :  Export trees to Nexus files--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3 -- -- Maintainer  :  dominik.schrempf@gmail.com
src/ELynx/Tree/Import/Newick.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Tree.Import.Newick -- Description :  Import Newick trees--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -225,7 +225,8 @@     *> optional brackets     *> treeRevBayes     <* char ';'-    <* skipWhile isSpace <?> "newickRevBayes"+    <* skipWhile isSpace+    <?> "newickRevBayes"  treeRevBayes :: Parser (Tree Phylo Name) treeRevBayes = branchedRevBayes <|> leafRevBayes <?> "treeRevBayes"
src/ELynx/Tree/Import/Nexus.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Tree.Import.Nexus -- Description :  Import trees from Nexus files--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3 -- -- Maintainer  :  dominik.schrempf@gmail.com
src/ELynx/Tree/Length.hs view
@@ -5,7 +5,7 @@ -- | -- Module      :  ELynx.Tree.Length -- Description :  Labels having a length--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
+ src/ELynx/Tree/Mrca.hs view
@@ -0,0 +1,112 @@+-- |+-- Module      :  ELynx.Tree.Mrca+-- Description :  Most recent common ancestors+-- Copyright   :  2022 Dominik Schrempf+-- License     :  GPL-3.0-or-later+--+-- Maintainer  :  dominik.schrempf@gmail.com+-- Stability   :  experimental+-- Portability :  portable+--+-- Creation date: Wed Jun 29 15:57:09 2022.+--+-- Specify nodes using most recent common ancestors (MRCA).+module ELynx.Tree.Mrca+  ( isAncestor,+    isMrca,+    getPathToMrca,+    getTreeAtMrca,+    findNode,+  )+where++import Data.Either+import Data.List+import Data.Maybe+import qualified Data.Set as S+import ELynx.Tree.Rooted+import ELynx.Tree.Zipper++-- | Test if the root node of the given tree is an ancestor of the given leaves.+isAncestor :: Ord a => S.Set a -> Tree e a -> Bool+--                      True if an x of xs is not in the collection of leaves.+--                False      if an x of xs is not in the collection of leaves. -> OK.+isAncestor xs t = not $ any (`S.notMember` lvs) xs+  where+    lvs = S.fromList $ leaves t++-- | Test if the root node of the given tree is the MRCA of the given leaves.+isMrca :: Ord a => S.Set a -> Tree e a -> Bool+--                                    True if any daughter forest is an ancestor.+--                               False     if any daughter forest is an ancestor. -> OK.+isMrca xs t = isAncestor xs t && not (any (isAncestor xs) (forest t))++-- | Get the path to the MRCA of the given leaves on the given tree.+--+-- Return 'Left' if:+--+-- - The tree has duplicate leaves.+--+-- - The MRCA cannot be found.+getPathToMrca :: (Ord a, Show a) => S.Set a -> Tree e a -> Either String Path+getPathToMrca ss tr+  | duplicateLeaves tr = Left "getPathToMrca: Tree contains duplicate leaves."+  | otherwise = tail <$> go 0 tr+  where+    go i t+      | isMrca ss t = Right [i]+      --                                    One path will be (Right p).+      | isAncestor ss t =+          case find isRight [go j t' | (j, t') <- zip [0 ..] (forest t)] of+            -- Use 'error' because this should never happen since one subtree has+            -- to contain the MRCA.+            Nothing -> error $ "getPathToMrca: BUG; I am ancestor but no MRCA found for: " <> show ss+            Just xs -> fmap (i :) xs+      | otherwise = Left $ "getPathToMrca: Could not get MRCA for: " <> show ss <> "."++getTreeAtMrca :: (Ord a, Show a) => S.Set a -> Tree e a -> Either String (Tree e a)+getTreeAtMrca ss tr+  | duplicateLeaves tr = Left "getTreeAtMrca: Tree contains duplicate leaves."+  | otherwise = go tr+  where+    go t+      | isMrca ss t = Right t+      --                                One tree will be (Right t).+      | isAncestor ss t = case find isRight $ map go (forest t) of+          -- Use 'error' because this should never happen since one subtree has+          -- to contain the MRCA.+          Nothing -> error $ "getTreeAtMrca: BUG; I am ancestor but no MRCA found for: " <> show ss+          Just x -> x+      | otherwise = Left $ "getTreeAtMrca: Could not get MRCA for: " <> show ss <> "."++-- | Find a node on a tree.+--+-- If the node is found, the node is specified by the MRCA of the returned set+-- of leaves. The set will contain one element if the node is a leaf, or two+-- elements, if the node is internal.+--+-- Return 'Left' if:+--+-- - The tree has duplicate labels.+--+-- - The node cannot be found.+--+-- - The node is a degree two node.+findNode :: (Ord a, Show a) => a -> Tree e a -> Either String (S.Set a)+findNode n t+  | duplicateLabels t = Left "findNode: tree contains duplicate labels"+  | otherwise = case go n t of+      Nothing -> Left $ "findNode: node not found: " <> show n+      Just x -> x+  where+    go x (Node _ y ts)+      | x == y = case ts of+          [] -> Just $ Right $ S.singleton x+          (l : r : _) -> Just $ Right $ S.fromList [head $ leaves l, head $ leaves r]+          _ -> Just $ Left $ "findNode: degree two node found" <> show n+      | otherwise = case catMaybes $ map (go x) ts of+          [] -> Nothing+          [z] -> Just z+          -- Use 'error' because this should never happen since we check for+          -- duplicate labels.+          _ -> error $ "findNode: BUG; more internal nodes found" <> show n
src/ELynx/Tree/Name.hs view
@@ -5,7 +5,7 @@ -- | -- Module      :  ELynx.Tree.Name -- Description :  Trees with named nodes--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
src/ELynx/Tree/Parallel.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Tree.Parallel -- Description :  Evaluation strategies for trees--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -40,11 +40,11 @@ parTree :: (NFData e, NFData a) => Int -> Strategy (Tree e a) parTree n t@(Node br lb ts)   | n == 1 = do-    ts' <- myParList rdeepseq ts-    return $ Node br lb ts'+      ts' <- myParList rdeepseq ts+      return $ Node br lb ts'   | n >= 2 = do-    ts' <- myParList (parTree (n -1)) ts-    return $ Node br lb ts'+      ts' <- myParList (parTree (n - 1)) ts+      return $ Node br lb ts'   | otherwise = rdeepseq t  branchFoldMap :: (e -> f) -> (f -> f -> f) -> Tree e a -> f
src/ELynx/Tree/Partition.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Tree.Partition -- Description :  Partitions on rose trees--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -109,8 +109,8 @@ partitions' _ (Node _ _ []) = S.empty partitions' p t@(Node _ _ ts) =   S.unions $-    either (const S.empty) S.singleton (pt (p : map label ts)) :-    zipWith partitions' cs ts+    either (const S.empty) S.singleton (pt (p : map label ts))+      : zipWith partitions' cs ts   where     cs = getComplementaryLeaves p t
src/ELynx/Tree/Phylogeny.hs view
@@ -5,7 +5,7 @@ -- | -- Module      :  ELynx.Tree.Phylogeny -- Description :  Phylogenetic trees--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -120,8 +120,8 @@ intersect ts   | S.null lvsCommon = Left "intersect: Intersection of leaves is empty."   | otherwise = case sequence [dropLeavesWith (predicate ls) t | (ls, t) <- zip leavesToDrop ts] of-    Nothing -> Left "intersect: A tree is empty."-    Just ts' -> Right ts'+      Nothing -> Left "intersect: A tree is empty."+      Just ts' -> Right ts'   where     -- Leaf sets.     lvss = map (S.fromList . leaves) ts@@ -186,10 +186,10 @@   | length lvLst /= S.size lvSet = Left "rootAt: Tree has duplicate leaves."   | toSet b /= lvSet = Left "rootAt: Bipartition does not match leaves of tree."   | otherwise = do-    ts <- roots t-    case find (\x -> bipartition x == Right b) ts of-      Nothing -> Left "rootAt': Bipartition not found on tree."-      Just t' -> Right t'+      ts <- roots t+      case find (\x -> bipartition x == Right b) ts of+        Nothing -> Left "rootAt': Bipartition not found on tree."+        Just t' -> Right t'   where     lvLst = leaves t     lvSet = S.fromList $ leaves t@@ -289,7 +289,7 @@     tR = Node bL' c $ tail ts  complementaryForests :: Tree e a -> Forest e a -> [Forest e a]-complementaryForests t ts = [t : take i ts ++ drop (i + 1) ts | i <- [0 .. (n -1)]]+complementaryForests t ts = [t : take i ts ++ drop (i + 1) ts | i <- [0 .. (n - 1)]]   where     n = length ts
src/ELynx/Tree/Rooted.hs view
@@ -6,7 +6,7 @@ -- | -- Module      :  ELynx.Tree.Rooted -- Description :  Rooted trees with labeled branches--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -77,6 +77,7 @@     setLabel,     modifyLabel,     labels,+    duplicateLabels,     setLabels,     identify, @@ -295,6 +296,10 @@   where     squish (Node _ lb ts) xs = lb : foldr squish xs ts +-- | Check if a tree has duplicate labels.+duplicateLabels :: Ord a => Tree e a -> Bool+duplicateLabels = duplicates . labels+ -- | Set node labels in pre-order. -- -- Return 'Nothing' if the provided list of node labels is too short.@@ -353,9 +358,9 @@ dropNodesWith p (Node br lb ts)   | p lb = Nothing   | otherwise =-    if null ts'-      then Nothing-      else Just $ Node br lb ts'+      if null ts'+        then Nothing+        else Just $ Node br lb ts'   where     ts' = mapMaybe (dropNodesWith p) ts
src/ELynx/Tree/Simulate/Coalescent.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Tree.Simulate.Coalescent -- Description :  Generate coalescent trees--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -14,46 +14,45 @@   ) where -import Control.Monad.Primitive import ELynx.Tree.Distribution.CoalescentContinuous import ELynx.Tree.Length import ELynx.Tree.Rooted import Statistics.Distribution-import System.Random.MWC+import System.Random.Stateful  -- | Simulate a coalescent tree with @n@ leaves. The branch lengths are in units -- of effective population size. simulate ::-  (PrimMonad m) =>+  StatefulGen g m =>   -- | Number of leaves.   Int ->-  Gen (PrimState m) ->+  g ->   m (Tree Length Int) simulate n = simulate' n 0 trs   where     trs = [Node 0 i [] | i <- [0 .. n - 1]]  simulate' ::-  (PrimMonad m) =>+  StatefulGen g m =>   Int ->   Int ->   Forest Length Int ->-  Gen (PrimState m) ->+  g ->   m (Tree Length Int) simulate' n a trs g   | n <= 0 = error "Cannot construct trees without leaves."   | n == 1 && length trs /= 1 = error "Too many trees provided."   | n == 1 && length trs == 1 = return $ head trs   | otherwise = do-    -- Indices of the leaves to join will be i-1 and i.-    i <- uniformR (1, n - 1) g-    -- The time of the coalescent event.-    t <- toLengthUnsafe <$> genContVar (coalescentDistributionCont n) g-    let trs' = map (modifyStem (+ t)) trs -- Move time 't' up on the tree.-        tl = trs' !! (i - 1)-        tr = trs' !! i-        -- Join the two chosen trees.-        tm = Node 0 a [tl, tr]-        -- Take the trees on the left, the merged tree, and the trees on the right.-        trs'' = take (i - 1) trs' ++ [tm] ++ drop (i + 1) trs'-    simulate' (n - 1) a trs'' g+      -- Indices of the leaves to join will be i-1 and i.+      i <- uniformRM (1, n - 1) g+      -- The time of the coalescent event.+      t <- toLengthUnsafe <$> genContVar (coalescentDistributionCont n) g+      let trs' = map (modifyStem (+ t)) trs -- Move time 't' up on the tree.+          tl = trs' !! (i - 1)+          tr = trs' !! i+          -- Join the two chosen trees.+          tm = Node 0 a [tl, tr]+          -- Take the trees on the left, the merged tree, and the trees on the right.+          trs'' = take (i - 1) trs' ++ [tm] ++ drop (i + 1) trs'+      simulate' (n - 1) a trs'' g
src/ELynx/Tree/Simulate/PointProcess.hs view
@@ -3,7 +3,7 @@ -- | --   Module      :  ELynx.Tree.Simulate.PointProcess --   Description :  Point process and functions---   Copyright   :  (c) Dominik Schrempf 2021+--   Copyright   :  2021 Dominik Schrempf --   License     :  GPL-3.0-or-later -- --   Maintainer  :  dominik.schrempf@gmail.com@@ -28,7 +28,6 @@ where  import Control.Monad-import Control.Monad.Primitive import Data.Function import Data.List import Data.Sequence (Seq)@@ -43,9 +42,7 @@ import ELynx.Tree.Length import ELynx.Tree.Rooted import qualified Statistics.Distribution as D-  ( genContVar,-  )-import System.Random.MWC+import System.Random.Stateful  -- Require near critical process if birth and death rates are closer than this value. epsNearCriticalPointProcess :: Double@@ -67,10 +64,10 @@ sortListWithIndices xs = sortBy (compare `on` fst) $ zip xs ([0 ..] :: [Int])  -- Insert element into random position of list.-randomInsertList :: PrimMonad m => a -> [a] -> Gen (PrimState m) -> m [a]+randomInsertList :: StatefulGen g m => a -> [a] -> g -> m [a] randomInsertList e v g = do   let l = length v-  i <- uniformR (0, l) g+  i <- uniformRM (0, l) g   return $ take i v ++ [e] ++ drop i v  -- | A __point process__ for \(n\) points and of age \(t_{or}\) is defined as@@ -98,7 +95,7 @@ -- | Sample a point process using the 'BirthDeathDistribution'. The names of the -- points will be integers. simulate ::-  (PrimMonad m) =>+  StatefulGen g m =>   -- | Number of points (samples).   Int ->   -- | Time of origin or MRCA.@@ -107,58 +104,57 @@   Rate ->   -- | Death rate.   Rate ->-  -- | Generator.-  Gen (PrimState m) ->+  g ->   m (PointProcess Int Double) simulate n ts l m g   | n < 1 = error "Number of samples needs to be one or larger."   | l < 0.0 = error "Birth rate needs to be positive."   | otherwise = case ts of-    Random -> simulateRandom n l m g-    Origin t -> simulateOrigin n t l m g-    Mrca t -> simulateMrca n t l m g+      Random -> simulateRandom n l m g+      Origin t -> simulateOrigin n t l m g+      Mrca t -> simulateMrca n t l m g  -- No time of origin given. We also don't need to take care of the conditioning -- (origin or MRCA). simulateRandom ::-  PrimMonad m =>+  StatefulGen g m =>   Int ->   Double ->   Double ->-  Gen (PrimState m) ->+  g ->   m (PointProcess Int Double) simulateRandom n l m g   | -- There is no formula for the over-critical process.     m > l =-    error-      "simulateRandom: Please specify height if mu > lambda."+      error+        "simulateRandom: Please specify height if mu > lambda."   | -- For the critical process, we have no idea about the time of origin, but can     -- use a specially derived distribution.     m =~= l =-    do-      !vs <- replicateM (n - 1) (D.genContVar (BDCNTD l) g)-      -- The length of the root branch will be 0.-      let t = maximum vs-      return $ PointProcess [0 .. (n - 1)] vs t+      do+        !vs <- replicateM (n - 1) (D.genContVar (BDCNTD l) g)+        -- The length of the root branch will be 0.+        let t = maximum vs+        return $ PointProcess [0 .. (n - 1)] vs t   | -- For the near critical process, we use a special distribution.     abs (m - l) <= epsNearCriticalTimeOfOrigin =-    do-      t <- D.genContVar (TONCD n l m) g-      simulateOrigin n t l m g+      do+        t <- D.genContVar (TONCD n l m) g+        simulateOrigin n t l m g   | -- For a sub-critical branching process, we can use the formula from Tanja Stadler.     otherwise =-    do-      t <- D.genContVar (TOD n l m) g-      simulateOrigin n t l m g+      do+        t <- D.genContVar (TOD n l m) g+        simulateOrigin n t l m g  -- Time of origin is given. simulateOrigin ::-  PrimMonad m =>+  StatefulGen g m =>   Int ->   Time ->   Double ->   Double ->-  Gen (PrimState m) ->+  g ->   m (PointProcess Int Double) simulateOrigin n t l m g   | t < 0.0 = error "simulateOrigin: Time of origin needs to be positive."@@ -171,38 +167,38 @@     -- 2. The near critical branching process.     -- 3. Normal values :).     m =~= l = do-    !vs <- replicateM (n - 1) (D.genContVar (BDCD t l) g)-    return $ PointProcess [0 .. (n - 1)] vs t+      !vs <- replicateM (n - 1) (D.genContVar (BDCD t l) g)+      return $ PointProcess [0 .. (n - 1)] vs t   | abs (m - l) <= epsNearCriticalPointProcess = do-    !vs <- replicateM (n - 1) (D.genContVar (BDNCD t l m) g)-    return $ PointProcess [0 .. (n - 1)] vs t+      !vs <- replicateM (n - 1) (D.genContVar (BDNCD t l m) g)+      return $ PointProcess [0 .. (n - 1)] vs t   | otherwise = do-    !vs <- replicateM (n - 1) (D.genContVar (BDD t l m) g)-    return $ PointProcess [0 .. (n - 1)] vs t+      !vs <- replicateM (n - 1) (D.genContVar (BDD t l m) g)+      return $ PointProcess [0 .. (n - 1)] vs t  -- Time of Mrca is given. simulateMrca ::-  PrimMonad m =>+  StatefulGen g m =>   Int ->   Time ->   Double ->   Double ->-  Gen (PrimState m) ->+  g ->   m (PointProcess Int Double) simulateMrca n t l m g   | t < 0.0 = error "simulateMrca: Time of MRCA needs to be positive."   | m =~= l = do-    !vs <- replicateM (n - 2) (D.genContVar (BDCD t l) g)-    vs' <- randomInsertList t vs g-    return $ PointProcess [0 .. (n - 1)] vs' t+      !vs <- replicateM (n - 2) (D.genContVar (BDCD t l) g)+      vs' <- randomInsertList t vs g+      return $ PointProcess [0 .. (n - 1)] vs' t   | abs (m - l) <= epsNearCriticalPointProcess = do-    !vs <- replicateM (n - 2) (D.genContVar (BDNCD t l m) g)-    vs' <- randomInsertList t vs g-    return $ PointProcess [0 .. (n - 1)] vs' t+      !vs <- replicateM (n - 2) (D.genContVar (BDNCD t l m) g)+      vs' <- randomInsertList t vs g+      return $ PointProcess [0 .. (n - 1)] vs' t   | otherwise = do-    !vs <- replicateM (n - 2) (D.genContVar (BDD t l m) g)-    vs' <- randomInsertList t vs g-    return $ PointProcess [0 .. (n - 1)] vs' t+      !vs <- replicateM (n - 2) (D.genContVar (BDD t l m) g)+      vs' <- randomInsertList t vs g+      return $ PointProcess [0 .. (n - 1)] vs' t  -- Sort the values of a point process and their indices to be (the indices -- that they will have while creating the tree).@@ -226,7 +222,7 @@  -- | See 'simulateReconstructedTree', but n times. simulateNReconstructedTrees ::-  (PrimMonad m) =>+  (StatefulGen g m) =>   -- | Number of trees   Int ->   -- | Number of points (samples)@@ -237,8 +233,7 @@   Rate ->   -- | Death rate   Rate ->-  -- | Generator (see 'System.Random.MWC')-  Gen (PrimState m) ->+  g ->   m (Forest Length Int) simulateNReconstructedTrees nT nP t l m g   | nT <= 0 = return []@@ -248,7 +243,7 @@ -- 'toReconstructedTree') possibly with specific height and a fixed number of -- leaves according to the birth and death process. simulateReconstructedTree ::-  (PrimMonad m) =>+  StatefulGen g m =>   -- | Number of points (samples)   Int ->   -- | Time of origin or MRCA@@ -257,8 +252,7 @@   Rate ->   -- | Death rate   Rate ->-  -- | Generator (see 'System.Random.MWC')-  Gen (PrimState m) ->+  g ->   m (Tree Length Int) simulateReconstructedTree n t l m g = do   PointProcess ns vs o <- simulate n t l m g@@ -284,7 +278,7 @@   | -- -- Test is deactivated.     -- -- | otherwise = if isReconstructed treeOrigin then treeOrigin else error "Error in algorithm."     otherwise =-    treeOrigin+      treeOrigin   where     (vsSorted, isSorted) = sortPP pp     !lvs = S.fromList [Node 0 p [] | p <- ps]
src/ELynx/Tree/Splittable.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Tree.Splittable -- Description :  Splittable branch labels--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
src/ELynx/Tree/Support.hs view
@@ -5,7 +5,7 @@ -- | -- Module      :  ELynx.Tree.Support -- Description :  Labels with support values--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
src/ELynx/Tree/Zipper.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Tree.Zipper -- Description :  Zippers on rooted rose trees with branch labels--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com@@ -133,13 +133,13 @@     | null ts -> Nothing     | length ts <= n -> Nothing     | otherwise ->-      Just $-        Pos-          { current = head rs',-            before = reverse ls',-            after = tail rs',-            parents = (before pos, br, lb, after pos) : parents pos-          }+        Just $+          Pos+            { current = head rs',+              before = reverse ls',+              after = tail rs',+              parents = (before pos, br, lb, after pos) : parents pos+            }     where       (ls', rs') = splitAt n ts @@ -151,12 +151,12 @@     | null ts -> error "goChildUnsafe: Forest is empty."     | length ts <= n -> error "goChildUnsafe: Forest is too short."     | otherwise ->-      Pos-        { current = head rs',-          before = reverse ls',-          after = tail rs',-          parents = (before pos, br, lb, after pos) : parents pos-        }+        Pos+          { current = head rs',+            before = reverse ls',+            after = tail rs',+            parents = (before pos, br, lb, after pos) : parents pos+          }     where       (ls', rs') = splitAt n ts
test/ELynx/ClassLaws.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.ClassLaws -- Description :  Unit tests for ELynx.ClassLaws--- Copyright   :  (c) 2021 Dominik Schrempf+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Topology/Arbitrary.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Topology.Arbitrary -- Description :  Arbitrary instances for topologies--- Copyright   :  (c) 2021 Dominik Schrempf+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Topology/RootedSpec.hs view
@@ -1,7 +1,7 @@ -- | -- Module      :  ELynx.Topology.RootedSpec -- Description :  Unit tests for ELynx.Topology.Rooted--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/Arbitrary.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Tree.Arbitrary -- Description :  Arbitrary instances for trees--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/BipartitionSpec.hs view
@@ -2,7 +2,7 @@  -- | -- Module      :  ELynx.Tree.BipartitionSpec--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/DistanceSpec.hs view
@@ -3,7 +3,7 @@  -- | -- Module      :  ELynx.Tree.DistanceSpec--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/Export/NewickSpec.hs view
@@ -2,7 +2,7 @@  -- | -- Module      :  ELynx.Tree.Export.NewickSpec--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/Export/NexusSpec.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Tree.Export.NexusSpec -- Description :  Test export of trees in Nexus files--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3 -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/Import/NewickSpec.hs view
@@ -2,7 +2,7 @@  -- | -- Module      :  ELynx.Tree.Import.NewickSpec--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/Import/NexusSpec.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Tree.Import.NexusSpec -- Description :  Test import of trees in Nexus files--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3 -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/PartitionSpec.hs view
@@ -2,7 +2,7 @@  -- | -- Module      :  ELynx.Tree.PartitionSpec--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/PhylogenySpec.hs view
@@ -4,7 +4,7 @@ -- | -- Module      :  ELynx.Tree.PhylogenySpec -- Description :  Unit tests for ELynx.Tree.Phylogeny--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/RootedSpec.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Tree.RootedSpec -- Description :  Unit tests for ELynx.Tree.Rooted--- Copyright   :  (c) Dominik Schrempf 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com
test/ELynx/Tree/SupportSpec.hs view
@@ -3,7 +3,7 @@ -- | -- Module      :  ELynx.Tree.SupportSpec -- Description :  Unit tests for ELynx.Tree.SupportSpec--- Copyright   :  (c) Dominik Schrempf, 2021+-- Copyright   :  2021 Dominik Schrempf -- License     :  GPL-3.0-or-later -- -- Maintainer  :  dominik.schrempf@gmail.com