{- Extremely simple program to introduce duplicates -}

module Duplicator where

import System.IO
import System.Environment (getArgs)
import Control.Monad (when)
import Bio.Sequence.Fasta

import Statistics

main :: IO ()
main = do
  args <- getArgs
  case args of 
       (p:rest) -> do
           let p' = read p 
               i = case rest of [] -> hReadFasta stdin 
                                [f] -> readFasta f
           when (p'<=0 || p'>=1) $ error "Probability should be between 0 and 1"
           i >>= evalRandIO . dup p' >>= hWriteFasta stdout
       _ -> error "Usage: duplicator pr\n   where pr is the (recursive) probability of duplicating a sequence"

dup :: RandomGen g => Double -> [Sequence] -> Rand g [Sequence]
dup p (x:xs) = do
  r <- getRandomR (0,1)
  ys <- if r < p then dup p (x:xs)
        else dup p xs
  return (x:ys)
dup _ [] = return []