diff options
| author | Justin Bedo <cu@cua0.org> | 2025-02-26 16:22:42 +1100 |
|---|---|---|
| committer | Justin Bedo <cu@cua0.org> | 2025-02-27 10:41:10 +1100 |
| commit | 75fadaddab4a3e50fa509dcc21d3024fe681ce67 (patch) | |
| tree | 6a16b3f1e1be365c22342b82311178086ad3f617 | |
| parent | 480b6e000afeec19aa65857232f1261d2b21de76 (diff) | |
add single site mh
| -rw-r--r-- | bin/cluster.hs | 17 | ||||
| -rw-r--r-- | cabal.project | 4 |
2 files changed, 13 insertions, 8 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) diff --git a/cabal.project b/cabal.project index 0677c45..38182e1 100644 --- a/cabal.project +++ b/cabal.project @@ -3,5 +3,5 @@ packages: *.cabal source-repository-package type: git location: https://vk3.wtf/cgit/ppl.git - tag: 368ebf9f96e4fded6b69ca609f9db53eb19b4589 - --sha256: 0jm5q15l9advppc5gfjmkpr4zmv53grm46bh1c0bxib0hd8mmpqq + tag: ec8321993c3b5195581e5286b4416b65a9f23fa3 + --sha256: 8ac3ef9672ef43f988c74f28c6ffe62bd0ff958f0b828451bcf8d64f7bf448ec |