diff options
Diffstat (limited to 'bin/cluster.hs')
| -rw-r--r-- | bin/cluster.hs | 17 |
1 files changed, 11 insertions, 6 deletions
diff --git a/bin/cluster.hs b/bin/cluster.hs index 46a1fa9..222834e 100644 --- a/bin/cluster.hs +++ b/bin/cluster.hs @@ -11,7 +11,7 @@ import qualified Data.Map as M import GHC.Compact import Numeric.Log hiding (sum) import Options.Applicative -import PPL hiding (binom) +import PPL hiding (Tree, binom) import qualified Streaming as S import Streaming.Prelude (Of, Stream, each, fold, yield) import qualified Streaming.Prelude as S @@ -110,17 +110,18 @@ main = run =<< execParser opts (parser <**> helper) (fullDesc <> progDesc "Infer a phylogeny from SNV calls in multiple samples" <> header "phylogey - Bayesian phylogeny inference") parser = - Options <$> option auto (long "seed" <> short 's' <> help "random seed" <> showDefault <> value 42 <> metavar "INT") + Options + <$> option auto (long "seed" <> short 's' <> help "random seed" <> showDefault <> value 42 <> metavar "INT") <*> option auto (long "nsamples" <> short 'n' <> help "number of samples from posterior" <> value 100000 <> metavar "INT" <> showDefault) - <*> option probability (long "mhfrac" <> short 'm' <> help "Metropolis-Hastings mutation probability" <> value 0.01 <> metavar "(0,1]" <> showDefault) + <*> option probability (long "mhfrac" <> short 'm' <> help "Metropolis-Hastings mutation probability" <> value 0.01 <> metavar "[0,1]" <> showDefault) <*> argument str (metavar "INPUT") <*> argument str (metavar "OUTPUTDIR") probability = eitherReader $ \arg -> case reads arg of - [(r, "")] -> if r <= 1 && r > 0 then Right r else Left "mhfrac not a valid probability" + [(r, "")] -> if r <= 1 && r >= 0 then Right r else Left "mhfrac not a valid probability" _ -> Left "mhfrac not a valid probability" -takeWithProgress :: S.MonadIO m => Int -> Stream (Of a) m r -> Stream (Of a) m () +takeWithProgress :: (S.MonadIO m) => Int -> Stream (Of a) m r -> Stream (Of a) m () takeWithProgress n str = do pb <- S.liftIO $ newProgressBar defStyle 10 (Progress 0 n ()) S.mapM (update pb) $ S.take n str @@ -135,7 +136,11 @@ run opts = do g = mkStdGen $ seed opts parsed <- compact $ map (map dbl . tail . words) lines hSetBuffering stdout NoBuffering - m@((ps, cl), _) <- S.fold_ (\l r -> if mml l < mml r then l else r) (([[]], []), -1 / 0) id . takeWithProgress (nsamples opts) $ mh g (mhfrac opts) (model $ getCompact parsed) + m@((ps, cl), _) <- + S.fold_ (\l r -> if mml l < mml r then l else r) (([[]], []), -1 / 0) id . takeWithProgress (nsamples opts) $ + if mhfrac opts <= 0 + then ssmh g (model $ getCompact parsed) + else mh g (mhfrac opts) (model $ getCompact parsed) writeFile (outputPath opts <> "/props") . unlines $ map (intercalate "," . map show) ps writeFile (outputPath opts <> "/tree") . unlines $ map show cl writeFile (outputPath opts <> "/mml") $ show (mml m) |