varan-0.3: src/RandomSelect.hs

-- Calculates p-values using a permutation test (not currently used)

module RandomSelect where

import System.Random
import Control.Monad.State
import Variants
import AgrestiCoull (confidenceInterval)

pval_max_count :: Int
pval_max_count = 1000

pval_acc :: Double
pval_acc = 0.01

pval :: RandomGen g => g -> ([Counts] -> Double) -> [Counts] -> (Double,Double)
pval g f cs = let thresh = f cs
                  xs = take pval_max_count $ rselect g cs
                  go suc tot (y:ys) = let (a,b) = confidenceInterval 1.0 suc (tot-suc) 
                                      in if b-a <= pval_acc then fromIntegral suc/fromIntegral tot 
                                         else go (suc + if f y >= thresh then 1 else 0) (tot+1) ys
                  go suc tot [] = fromIntegral suc/fromIntegral tot
              in if thresh > 0 then (thresh, go 0 0 xs) else (0,1)

type AlleleSample = (Int,Int,Int,Int)

-- generate an infinite stream of sampled allele distributions
rselect :: RandomGen g => g -> [Counts] -> [[Counts]]
rselect g cs = map fromLists $ fst $ runState sel (g,sum_al,pool_sizes)
  where (_tot,sum_al,pool_sizes) = count cs
        
fromLists :: [AlleleSample] -> [Counts]
fromLists = map f
  where f (a,b,c,d) = C a b c d []
        
sel :: RandomGen g => State (g,AlleleSample,[Int]) [[AlleleSample]]
sel = do
  a <- select1
  as <- sel
  return (a:as)
        
select1 :: RandomGen g => State (g,AlleleSample,[Int]) [AlleleSample]
select1 = do
  -- select randomly for each p_sz
  (g,s_al,p_sz) <- get
  let (g',g'') = split g
  put (g'',s_al,p_sz) -- restore state
  return $ pickNs g' p_sz s_al

-- pick a random sampling of alleles given pool sizes
pickNs :: RandomGen g => g -> [Int] -> AlleleSample -> [AlleleSample]
pickNs _ [] _als = [] -- verify that als is empty
pickNs g (p:ps) als = let
  (g1,g2) = split g
  in pickN' g1 als p : pickNs g2 ps als

-- sample alleles given randomgen, count, and sum alleles (i.e. distribution)
pickN' :: RandomGen g => g -> AlleleSample -> Int -> AlleleSample
pickN' g als = go (0,0,0,0) g
  where
        go acc _ 0 = acc
        go acc g1 cnt = 
          let (i,g2) = randomR (1,sum' als) g1
          in go (add i als acc) g2 (cnt-1)
        sum' (a,b,c,d) = a+b+c+d
        
-- count total alleles, number of each variant, and number in each count
count :: [Counts] -> (Int,AlleleSample,[Int])
count cs = let xs = sumList . map toList $ cs
               ys = map (sum . toList) cs
               tuple [a,b,c,d] = (a,b,c,d)
           in (sum xs,tuple xs,ys)

add :: Int -> AlleleSample -> AlleleSample -> AlleleSample
add i (a,b,c,d) (x,y,z,w)
  | i <= a = (x+1,y,z,w)
  | i <= a+b = (x,y+1,z,w)
  | i <= a+b+c = (x,y,z+1,w)               
  | i <= a+b+c+d = (x,y,z,w+1)               
  | otherwise = error "too high i"

{-
-- remove allele number x from input
pick :: Int -> AlleleSample -> AlleleSample
pick x (a:as) | x <= a = 1 : map (const 0) as
              | otherwise = 0 : pick (x-a) as
pick _ [] = error "pick ran out of input, x too large"
-}