{-# LANGUAGE OverloadedStrings #-}
module SLynx.Simulate.Simulate
( simulateCmd,
)
where
import Control.Applicative ((<|>))
import Control.Monad
import Control.Monad.IO.Class
import Control.Monad.Trans.Class
import Control.Monad.Trans.Reader (ask)
import qualified Data.ByteString.Builder as BB
import qualified Data.ByteString.Lazy.Char8 as BL
import Data.List
import Data.Maybe
import qualified Data.Set as Set
import qualified Data.Vector as V
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Unboxed as U
import ELynx.Alphabet.Alphabet as A
import ELynx.Import.MarkovProcess.EDMModelPhylobayes
import ELynx.Import.MarkovProcess.SiteprofilesPhylobayes
import qualified ELynx.MarkovProcess.AminoAcid as MA
import ELynx.MarkovProcess.GammaRateHeterogeneity
import qualified ELynx.MarkovProcess.MixtureModel as MM
import qualified ELynx.MarkovProcess.PhyloModel as MP
import qualified ELynx.MarkovProcess.RateMatrix as MR
import qualified ELynx.MarkovProcess.SubstitutionModel as MS
import ELynx.Sequence.Export.Fasta
import qualified ELynx.Sequence.Sequence as Seq hiding
( name,
)
import ELynx.Simulate.MarkovProcessAlongTree
import ELynx.Tools.ByteString
import ELynx.Tools.Definitions
import ELynx.Tools.ELynx
import ELynx.Tools.Environment
import ELynx.Tools.InputOutput
import ELynx.Tools.Logger
import ELynx.Tools.Options
import ELynx.Tools.Reproduction
import ELynx.Tree
import qualified Numeric.LinearAlgebra as L
import SLynx.Simulate.Options
import SLynx.Simulate.PhyloModel
import System.Random.Stateful
import Text.Printf
dispv :: Int -> VS.Vector L.R -> String
dispv :: Int -> Vector Double -> String
dispv Int
p Vector Double
v = forall a. [a] -> a
head forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
tail forall a b. (a -> b) -> a -> b
$ String -> [String]
lines forall a b. (a -> b) -> a -> b
$ Int -> Matrix Double -> String
L.dispf Int
p (forall a. Storable a => Vector a -> Matrix a
L.asRow Vector Double
v)
dispmi :: Int -> Int -> L.Matrix L.R -> String
dispmi :: Int -> Int -> Matrix Double -> String
dispmi Int
p Int
i Matrix Double
m =
forall a. [a] -> [[a]] -> [a]
intercalate String
"\n" forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a. Int -> a -> [a]
replicate Int
i Char
' ' forall a. [a] -> [a] -> [a]
++) forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
tail forall a b. (a -> b) -> a -> b
$ String -> [String]
lines forall a b. (a -> b) -> a -> b
$ Int -> Matrix Double -> String
L.dispf Int
p Matrix Double
m
getDistLine :: Int -> MR.StationaryDistribution -> BB.Builder
getDistLine :: Int -> Vector Double -> Builder
getDistLine Int
i Vector Double
d =
Int -> Builder
BB.intDec Int
i
forall a. Semigroup a => a -> a -> a
<> Char -> Builder
BB.char8 Char
' '
forall a. Semigroup a => a -> a -> a
<> Builder
s
where
s :: Builder
s = forall a. Monoid a => [a] -> a
mconcat forall a b. (a -> b) -> a -> b
$ forall a. a -> [a] -> [a]
intersperse (Char -> Builder
BB.char8 Char
' ') forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map Double -> Builder
BB.doubleDec forall a b. (a -> b) -> a -> b
$ forall a. Storable a => Vector a -> [a]
VS.toList Vector Double
d
writeSiteDists :: [Int] -> V.Vector MR.StationaryDistribution -> ELynx SimulateArguments ()
writeSiteDists :: [Int] -> Vector (Vector Double) -> ELynx SimulateArguments ()
writeSiteDists [Int]
componentIs Vector (Vector Double)
ds = do
Maybe String
mbn <- GlobalArguments -> Maybe String
outFileBaseName forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Environment a -> GlobalArguments
globalArguments forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall (m :: * -> *) r. Monad m => ReaderT r m r
ask
case Maybe String
mbn of
Maybe String
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return ()
Just String
bn -> forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ String -> ByteString -> IO ()
BL.writeFile (String
bn forall a. Semigroup a => a -> a -> a
<> String
".sitedists") ByteString
output
where
dsPaml :: Vector (Vector Double)
dsPaml = forall a b. (a -> b) -> Vector a -> Vector b
V.map Vector Double -> Vector Double
MA.alphaToPamlVec Vector (Vector Double)
ds
lns :: [Builder]
lns = [Int -> Vector Double -> Builder
getDistLine Int
i Vector Double
d | (Int
i, Int
c) <- forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1 ..] [Int]
componentIs, let d :: Vector Double
d = Vector (Vector Double)
dsPaml forall a. Vector a -> Int -> a
V.! Int
c]
output :: ByteString
output = Builder -> ByteString
BB.toLazyByteString forall a b. (a -> b) -> a -> b
$ forall a. Monoid a => [a] -> a
mconcat forall a b. (a -> b) -> a -> b
$ forall a. a -> [a] -> [a]
intersperse (Char -> Builder
BB.char8 Char
'\n') [Builder]
lns
simulateAlignment ::
(RandomGen g, HasLength e, HasName a) =>
MP.PhyloModel ->
Tree e a ->
Int ->
IOGenM g ->
ELynx SimulateArguments ()
simulateAlignment :: forall g e a.
(RandomGen g, HasLength e, HasName a) =>
PhyloModel
-> Tree e a -> Int -> IOGenM g -> ELynx SimulateArguments ()
simulateAlignment PhyloModel
pm Tree e a
t' Int
n IOGenM g
g = do
let t :: Tree Double
t = Length -> Double
fromLength forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall e. HasLength e => e -> Length
getLength forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall e a. Tree e a -> Tree e
toTreeBranchLabels Tree e a
t'
[[Int]]
leafStates <- case PhyloModel
pm of
MP.SubstitutionModel SubstitutionModel
sm -> forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall g.
RandomGen g =>
Int
-> Vector Double
-> Matrix Double
-> Tree Double
-> IOGenM g
-> IO [[Int]]
simulateAndFlattenPar Int
n Vector Double
d Matrix Double
e Tree Double
t IOGenM g
g
where
d :: Vector Double
d = SubstitutionModel -> Vector Double
MS.stationaryDistribution SubstitutionModel
sm
e :: Matrix Double
e = SubstitutionModel -> Matrix Double
MS.exchangeabilityMatrix SubstitutionModel
sm
MP.MixtureModel MixtureModel
mm -> do
([Int]
cs, [[Int]]
ss) <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall g.
RandomGen g =>
Int
-> Vector Double
-> Vector (Vector Double)
-> Vector (Matrix Double)
-> Tree Double
-> IOGenM g
-> IO ([Int], [[Int]])
simulateAndFlattenMixtureModelPar Int
n Vector Double
ws Vector (Vector Double)
ds Vector (Matrix Double)
es Tree Double
t IOGenM g
g
[Int] -> Vector (Vector Double) -> ELynx SimulateArguments ()
writeSiteDists [Int]
cs Vector (Vector Double)
ds
forall (m :: * -> *) a. Monad m => a -> m a
return [[Int]]
ss
where
ws :: Vector Double
ws = MixtureModel -> Vector Double
MM.getWeights MixtureModel
mm
ds :: Vector (Vector Double)
ds = forall a b. (a -> b) -> Vector a -> Vector b
V.map SubstitutionModel -> Vector Double
MS.stationaryDistribution forall a b. (a -> b) -> a -> b
$ MixtureModel -> Vector SubstitutionModel
MM.getSubstitutionModels MixtureModel
mm
es :: Vector (Matrix Double)
es = forall a b. (a -> b) -> Vector a -> Vector b
V.map SubstitutionModel -> Matrix Double
MS.exchangeabilityMatrix forall a b. (a -> b) -> a -> b
$ MixtureModel -> Vector SubstitutionModel
MM.getSubstitutionModels MixtureModel
mm
let leafNames :: [Name]
leafNames = forall a b. (a -> b) -> [a] -> [b]
map forall a. HasName a => a -> Name
getName forall a b. (a -> b) -> a -> b
$ forall e a. Tree e a -> [a]
leaves Tree e a
t'
code :: Alphabet
code = PhyloModel -> Alphabet
MP.getAlphabet PhyloModel
pm
alph :: Set Character
alph = AlphabetSpec -> Set Character
A.all forall a b. (a -> b) -> a -> b
$ Alphabet -> AlphabetSpec
alphabetSpec Alphabet
code
sequences :: [Sequence]
sequences =
[ ByteString -> ByteString -> Alphabet -> Characters -> Sequence
Seq.Sequence (Name -> ByteString
fromName Name
sName) ByteString
"" Alphabet
code (forall a. Unbox a => [a] -> Vector a
U.fromList forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a. Int -> Set a -> a
`Set.elemAt` Set Character
alph) [Int]
ss)
| (Name
sName, [Int]
ss) <- forall a b. [a] -> [b] -> [(a, b)]
zip [Name]
leafNames [[Int]]
leafStates
]
output :: ByteString
output = [Sequence] -> ByteString
sequencesToFasta [Sequence]
sequences
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS String
""
forall a.
Reproducible a =>
String -> ByteString -> String -> ELynx a ()
out String
"simulated multi sequence alignment" ByteString
output String
".fasta"
summarizeEDMComponents :: [EDMComponent] -> BL.ByteString
summarizeEDMComponents :: [EDMComponent] -> ByteString
summarizeEDMComponents [EDMComponent]
cs =
String -> ByteString
BL.pack forall a b. (a -> b) -> a -> b
$
String
"Empiricial distribution mixture model with "
forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show (forall (t :: * -> *) a. Foldable t => t a -> Int
length [EDMComponent]
cs)
forall a. [a] -> [a] -> [a]
++ String
" components."
reportModel :: MP.PhyloModel -> ELynx SimulateArguments ()
reportModel :: PhyloModel -> ELynx SimulateArguments ()
reportModel PhyloModel
m = do
GlobalArguments
as <- forall a. Environment a -> GlobalArguments
globalArguments forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall (m :: * -> *) r. Monad m => ReaderT r m r
ask
if GlobalArguments -> Bool
writeElynxFile GlobalArguments
as
then
( do
let bn :: Maybe String
bn = GlobalArguments -> Maybe String
outFileBaseName GlobalArguments
as
case Maybe String
bn of
Maybe String
Nothing ->
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS
String
"No output file provided; omit writing machine-readable phylogenetic model."
Just String
_ ->
forall a.
Reproducible a =>
String -> ByteString -> String -> ELynx a ()
out String
"model definition (machine readable)" (String -> ByteString
BL.pack (forall a. Show a => a -> String
show PhyloModel
m) forall a. Semigroup a => a -> a -> a
<> ByteString
"\n") String
".model.gz"
)
else forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS String
"No elynx file required; omit writing machine-readable phylogenetic model."
pretty :: Length -> String
pretty :: Length -> String
pretty = forall r. PrintfType r => String -> r
printf String
"%.5f" forall b c a. (b -> c) -> (a -> b) -> a -> c
. Length -> Double
fromLength
prettyRow :: String -> String -> BL.ByteString
prettyRow :: String -> String -> ByteString
prettyRow String
name String
val = Int -> ByteString -> ByteString
alignLeft Int
33 ByteString
n forall a. Semigroup a => a -> a -> a
<> Int -> ByteString -> ByteString
alignRight Int
8 ByteString
v
where
n :: ByteString
n = String -> ByteString
BL.pack String
name
v :: ByteString
v = String -> ByteString
BL.pack String
val
summarizeLengths :: HasLength e => Tree e a -> BL.ByteString
summarizeLengths :: forall e a. HasLength e => Tree e a -> ByteString
summarizeLengths Tree e a
t =
ByteString -> [ByteString] -> ByteString
BL.intercalate
ByteString
"\n"
[ String -> String -> ByteString
prettyRow String
"Origin height: " forall a b. (a -> b) -> a -> b
$ Length -> String
pretty Length
h,
String -> String -> ByteString
prettyRow String
"Average distance origin to leaves: " forall a b. (a -> b) -> a -> b
$ Length -> String
pretty Length
h',
String -> String -> ByteString
prettyRow String
"Total branch length: " forall a b. (a -> b) -> a -> b
$ Length -> String
pretty Length
b
]
where
n :: Int
n = forall (t :: * -> *) a. Foldable t => t a -> Int
length forall a b. (a -> b) -> a -> b
$ forall e a. Tree e a -> [a]
leaves Tree e a
t
h :: Length
h = forall e a. HasLength e => Tree e a -> Length
height Tree e a
t
h' :: Length
h' = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum (forall e a. HasLength e => Tree e a -> [Length]
distancesOriginLeaves Tree e a
t) forall a. Fractional a => a -> a -> a
/ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
b :: Length
b = forall e a. HasLength e => Tree e a -> Length
totalBranchLength Tree e a
t
roundN :: Int -> Double -> Double
roundN :: Int -> Double -> Double
roundN Int
n Double
v = forall a. Num a => Integer -> a
fromInteger (forall a b. (RealFrac a, Integral b) => a -> b
round forall a b. (a -> b) -> a -> b
$ Double
v forall a. Num a => a -> a -> a
* (Double
10 forall a b. (Num a, Integral b) => a -> b -> a
^ Int
n)) forall a. Fractional a => a -> a -> a
/ (Double
10.0 forall a b. (Fractional a, Integral b) => a -> b -> a
^^ Int
n)
summarizeSM :: MS.SubstitutionModel -> [BL.ByteString]
summarizeSM :: SubstitutionModel -> [ByteString]
summarizeSM SubstitutionModel
sm =
forall a b. (a -> b) -> [a] -> [b]
map String -> ByteString
BL.pack forall a b. (a -> b) -> a -> b
$
(forall a. Show a => a -> String
show (SubstitutionModel -> Alphabet
MS.alphabet SubstitutionModel
sm) forall a. [a] -> [a] -> [a]
++ String
" substitution model: " forall a. [a] -> [a] -> [a]
++ SubstitutionModel -> String
MS.name SubstitutionModel
sm forall a. [a] -> [a] -> [a]
++ String
".")
forall a. a -> [a] -> [a]
: [String
"Parameters: " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show (SubstitutionModel -> Params
MS.params SubstitutionModel
sm) forall a. [a] -> [a] -> [a]
++ String
"." | Bool -> Bool
not (forall (t :: * -> *) a. Foldable t => t a -> Bool
null (SubstitutionModel -> Params
MS.params SubstitutionModel
sm))]
forall a. [a] -> [a] -> [a]
++ case SubstitutionModel -> Alphabet
MS.alphabet SubstitutionModel
sm of
Alphabet
DNA ->
[ String
"Stationary distribution: "
forall a. [a] -> [a] -> [a]
++ Int -> Vector Double -> String
dispv Int
precision (SubstitutionModel -> Vector Double
MS.stationaryDistribution SubstitutionModel
sm)
forall a. [a] -> [a] -> [a]
++ String
".",
String
"Exchangeability matrix:\n"
forall a. [a] -> [a] -> [a]
++ Int -> Int -> Matrix Double -> String
dispmi Int
2 Int
precision (SubstitutionModel -> Matrix Double
MS.exchangeabilityMatrix SubstitutionModel
sm)
forall a. [a] -> [a] -> [a]
++ String
".",
String
"Scale: " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show (Int -> Double -> Double
roundN Int
precision forall a b. (a -> b) -> a -> b
$ SubstitutionModel -> Double
MS.totalRate SubstitutionModel
sm) forall a. [a] -> [a] -> [a]
++ String
"."
]
Alphabet
Protein ->
[ String
"Stationary distribution: "
forall a. [a] -> [a] -> [a]
++ Int -> Vector Double -> String
dispv Int
precision (SubstitutionModel -> Vector Double
MS.stationaryDistribution SubstitutionModel
sm)
forall a. [a] -> [a] -> [a]
++ String
".",
String
"Scale: " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show (Int -> Double -> Double
roundN Int
precision forall a b. (a -> b) -> a -> b
$ SubstitutionModel -> Double
MS.totalRate SubstitutionModel
sm) forall a. [a] -> [a] -> [a]
++ String
"."
]
Alphabet
_ ->
forall a. HasCallStack => String -> a
error
String
"Extended character sets are not supported with substitution models."
summarizeMMComponent :: MM.Component -> [BL.ByteString]
summarizeMMComponent :: Component -> [ByteString]
summarizeMMComponent Component
c =
String -> ByteString
BL.pack String
"Weight: "
forall a. Semigroup a => a -> a -> a
<> (Builder -> ByteString
BB.toLazyByteString forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Builder
BB.doubleDec forall a b. (a -> b) -> a -> b
$ Component -> Double
MM.weight Component
c)
forall a. a -> [a] -> [a]
: SubstitutionModel -> [ByteString]
summarizeSM (Component -> SubstitutionModel
MM.substModel Component
c)
summarizeMM :: MM.MixtureModel -> [BL.ByteString]
summarizeMM :: MixtureModel -> [ByteString]
summarizeMM MixtureModel
m =
[ String -> ByteString
BL.pack forall a b. (a -> b) -> a -> b
$ String
"Mixture model: " forall a. [a] -> [a] -> [a]
++ MixtureModel -> String
MM.name MixtureModel
m forall a. [a] -> [a] -> [a]
++ String
".",
String -> ByteString
BL.pack forall a b. (a -> b) -> a -> b
$ String
"Number of components: " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show Int
n forall a. [a] -> [a] -> [a]
++ String
"."
]
forall a. [a] -> [a] -> [a]
++ [ByteString]
detail
where
n :: Int
n = forall (t :: * -> *) a. Foldable t => t a -> Int
length forall a b. (a -> b) -> a -> b
$ MixtureModel -> Vector Component
MM.components MixtureModel
m
detail :: [ByteString]
detail =
if Int
n forall a. Ord a => a -> a -> Bool
<= Int
100
then
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat
[ String -> ByteString
BL.pack (String
"Component " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show Int
i forall a. [a] -> [a] -> [a]
++ String
":") forall a. a -> [a] -> [a]
: Component -> [ByteString]
summarizeMMComponent Component
c
| (Int
i, Component
c) <- forall a b. [a] -> [b] -> [(a, b)]
zip [Int
1 :: Int ..] (forall a. Vector a -> [a]
V.toList forall a b. (a -> b) -> a -> b
$ MixtureModel -> Vector Component
MM.components MixtureModel
m)
]
else []
summarizePM :: MP.PhyloModel -> [BL.ByteString]
summarizePM :: PhyloModel -> [ByteString]
summarizePM (MP.MixtureModel MixtureModel
mm) = MixtureModel -> [ByteString]
summarizeMM MixtureModel
mm
summarizePM (MP.SubstitutionModel SubstitutionModel
sm) = SubstitutionModel -> [ByteString]
summarizeSM SubstitutionModel
sm
simulateCmd :: ELynx SimulateArguments ()
simulateCmd :: ELynx SimulateArguments ()
simulateCmd = do
SimulateArguments
l <- forall a. Environment a -> a
localArguments forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall (m :: * -> *) r. Monad m => ReaderT r m r
ask
let treeFile :: String
treeFile = SimulateArguments -> String
argsTreeFile SimulateArguments
l
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS String
""
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS forall a b. (a -> b) -> a -> b
$ String
"Read tree from file '" forall a. [a] -> [a] -> [a]
++ String
treeFile forall a. [a] -> [a] -> [a]
++ String
"'."
Tree Phylo Name
tree <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. Parser a -> String -> IO a
parseFileWith (NewickFormat -> Parser (Tree Phylo Name)
newick NewickFormat
Standard) String
treeFile
let t' :: Tree Length Name
t' = forall a c b. (a -> c) -> (b -> c) -> Either a b -> c
either forall a. HasCallStack => String -> a
error forall a. a -> a
id forall a b. (a -> b) -> a -> b
$ forall e a.
HasMaybeLength e =>
Tree e a -> Either String (Tree Length a)
toLengthTree Tree Phylo Name
tree
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS forall a b. (a -> b) -> a -> b
$ String
"Number of leaves: " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show (forall (t :: * -> *) a. Foldable t => t a -> Int
length forall a b. (a -> b) -> a -> b
$ forall e a. Tree e a -> [a]
leaves Tree Length Name
t')
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logInfoB forall a b. (a -> b) -> a -> b
$ forall e a. HasLength e => Tree e a -> ByteString
summarizeLengths Tree Length Name
t'
let edmFile :: Maybe String
edmFile = SimulateArguments -> Maybe String
argsEDMFile SimulateArguments
l
let sProfileFiles :: Maybe [String]
sProfileFiles = SimulateArguments -> Maybe [String]
argsSiteprofilesFiles SimulateArguments
l
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS String
""
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logDebugS String
"Read EDM file or siteprofile files."
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (forall a. Maybe a -> Bool
isJust Maybe String
edmFile Bool -> Bool -> Bool
&& forall a. Maybe a -> Bool
isJust Maybe [String]
sProfileFiles) forall a b. (a -> b) -> a -> b
$
forall a. HasCallStack => String -> a
error String
"Got both: --edm-file and --siteprofile-files."
Maybe [EDMComponent]
edmCs <- case Maybe String
edmFile of
Maybe String
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a. Maybe a
Nothing
Just String
edmF -> do
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS String
"Read EDM file."
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall a. a -> Maybe a
Just forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall a. Parser a -> String -> IO a
parseFileWith Parser [EDMComponent]
phylobayes String
edmF
forall b a. b -> (a -> b) -> Maybe a -> b
maybe
(forall (m :: * -> *) a. Monad m => a -> m a
return ())
(forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logInfoB forall b c a. (b -> c) -> (a -> b) -> a -> c
. [EDMComponent] -> ByteString
summarizeEDMComponents)
Maybe [EDMComponent]
edmCs
Maybe [EDMComponent]
sProfiles <- case Maybe [String]
sProfileFiles of
Maybe [String]
Nothing -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a. Maybe a
Nothing
Just [String]
fns -> do
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS forall a b. (a -> b) -> a -> b
$
String
"Read siteprofiles from "
forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show (forall (t :: * -> *) a. Foldable t => t a -> Int
length [String]
fns)
forall a. [a] -> [a] -> [a]
++ String
" file(s)."
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logDebugS forall a b. (a -> b) -> a -> b
$ String
"The file names are:" forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show [String]
fns
[[EDMComponent]]
xs <- forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (forall a. Parser a -> String -> IO a
parseFileWith Parser [EDMComponent]
siteprofiles) [String]
fns
forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[EDMComponent]]
xs
forall b a. b -> (a -> b) -> Maybe a -> b
maybe
(forall (m :: * -> *) a. Monad m => a -> m a
return ())
(forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logInfoB forall b c a. (b -> c) -> (a -> b) -> a -> c
. [EDMComponent] -> ByteString
summarizeEDMComponents)
Maybe [EDMComponent]
sProfiles
let edmCsOrSiteprofiles :: Maybe [EDMComponent]
edmCsOrSiteprofiles = Maybe [EDMComponent]
edmCs forall (f :: * -> *) a. Alternative f => f a -> f a -> f a
<|> Maybe [EDMComponent]
sProfiles
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS String
"Read model string."
let ms :: Maybe String
ms = SimulateArguments -> Maybe String
argsSubstitutionModelString SimulateArguments
l
mm :: Maybe String
mm = SimulateArguments -> Maybe String
argsMixtureModelString SimulateArguments
l
mws :: Maybe Params
mws = SimulateArguments -> Maybe Params
argsMixtureWeights SimulateArguments
l
eitherPhyloModel' :: Either String PhyloModel
eitherPhyloModel' = Maybe String
-> Maybe String
-> Maybe Params
-> Maybe [EDMComponent]
-> Either String PhyloModel
getPhyloModel Maybe String
ms Maybe String
mm Maybe Params
mws Maybe [EDMComponent]
edmCsOrSiteprofiles
PhyloModel
phyloModel' <- case Either String PhyloModel
eitherPhyloModel' of
Left String
err -> forall (t :: (* -> *) -> * -> *) (m :: * -> *) a.
(MonadTrans t, Monad m) =>
m a -> t m a
lift forall a b. (a -> b) -> a -> b
$ forall a. HasCallStack => String -> a
error String
err
Right PhyloModel
pm -> forall (m :: * -> *) a. Monad m => a -> m a
return PhyloModel
pm
let maybeGammaParams :: Maybe GammaRateHeterogeneityParams
maybeGammaParams = SimulateArguments -> Maybe GammaRateHeterogeneityParams
argsGammaParams SimulateArguments
l
PhyloModel
phyloModel <- case Maybe GammaRateHeterogeneityParams
maybeGammaParams of
Maybe GammaRateHeterogeneityParams
Nothing -> do
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logInfoB forall a b. (a -> b) -> a -> b
$ [ByteString] -> ByteString
BL.unlines forall a b. (a -> b) -> a -> b
$ PhyloModel -> [ByteString]
summarizePM PhyloModel
phyloModel'
forall (m :: * -> *) a. Monad m => a -> m a
return PhyloModel
phyloModel'
Just (Int
n, Double
alpha) -> do
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logInfoB forall a b. (a -> b) -> a -> b
$ ByteString -> [ByteString] -> ByteString
BL.intercalate ByteString
"\n" forall a b. (a -> b) -> a -> b
$ PhyloModel -> [ByteString]
summarizePM PhyloModel
phyloModel'
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS String
""
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
ByteString -> Logger e ()
logInfoB forall a b. (a -> b) -> a -> b
$ ByteString -> [ByteString] -> ByteString
BL.intercalate ByteString
"\n" forall a b. (a -> b) -> a -> b
$ Int -> Double -> [ByteString]
summarizeGammaRateHeterogeneity Int
n Double
alpha
forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Int -> Double -> PhyloModel -> PhyloModel
expand Int
n Double
alpha PhyloModel
phyloModel'
PhyloModel -> ELynx SimulateArguments ()
reportModel PhyloModel
phyloModel
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS String
"Simulate alignment."
let alignmentLength :: Int
alignmentLength = SimulateArguments -> Int
argsLength SimulateArguments
l
forall e.
(HasLock e, HasLogHandles e, HasVerbosity e) =>
String -> Logger e ()
logInfoS forall a b. (a -> b) -> a -> b
$ String
"Length: " forall a. Semigroup a => a -> a -> a
<> forall a. Show a => a -> String
show Int
alignmentLength forall a. Semigroup a => a -> a -> a
<> String
"."
IOGenM StdGen
gen <- forall (m :: * -> *) g. MonadIO m => g -> m (IOGenM g)
newIOGenM forall a b. (a -> b) -> a -> b
$ Int -> StdGen
mkStdGen forall a b. (a -> b) -> a -> b
$ case SimulateArguments -> SeedOpt
argsSeed SimulateArguments
l of
SeedOpt
RandomUnset -> forall a. HasCallStack => String -> a
error String
"simulateCmd: seed not available; please contact maintainer."
RandomSet Int
s -> Int
s
Fixed Int
s -> Int
s
forall g e a.
(RandomGen g, HasLength e, HasName a) =>
PhyloModel
-> Tree e a -> Int -> IOGenM g -> ELynx SimulateArguments ()
simulateAlignment PhyloModel
phyloModel Tree Length Name
t' Int
alignmentLength IOGenM StdGen
gen