bugfixesnondiploid

1 files changed, 8 insertions, 4 deletions

diff --git a/bin/cluster.hs b/bin/cluster.hs
index e78e261..1fe4d2c 100644
--- a/bin/cluster.hs
+++ b/bin/cluster.hs
@@ -1,4 +1,5 @@
 {-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TemplateHaskell #-}
 {-# LANGUAGE ViewPatterns #-}
 
 module Main where
@@ -9,6 +10,7 @@ import Data.Foldable (toList)
 import Data.List hiding (group)
 import qualified Data.Map as M
 import GHC.Compact
+import qualified Language.Haskell.TH.Syntax as TH
 import Numeric.Log hiding (sum)
 import Options.Applicative
 import PPL hiding (Tree, binom)
@@ -48,7 +50,7 @@ instance Foldable Tree where
       bftrav ((Tree a l r) : ts) = f a <> bftrav (ts <> [l, r])
 
 {-# INLINE group #-}
-group (a : b : c: d: rs) = (a, b, c, d) : group rs
+group (a : b : c : d : rs) = (a, b, c, d) : group rs
 group [] = []
 group s = error $ "unexpected number of columns, expecting 4:" <> show s
 
@@ -59,7 +61,7 @@ group s = error $ "unexpected number of columns, expecting 4:" <> show s
 treeFromList (x : xs) = Tree x (treeFromList lpart) (treeFromList rpart)
   where
     (lpart, rpart) = unzip $ pairs xs
-    pairs (a:b:rs) = (a, b) : pairs rs
+    pairs (a : b : rs) = (a, b) : pairs rs
 
 -- Constrain trees so leaves sum to node value
 normTree :: Tree Double -> Tree Double
@@ -85,12 +87,14 @@ model xs = do
   mapM_ scoreLog $ zipWith likelihood params xs
   let cls = take (length xs) clusters
       k = maximum cls + 1
-      n = length (head xs) `div` 2
+      n = length (head xs) `div` 4
   pure (map (take k) $ take n ps, cls)
   where
     likelihood ps cnts = product $ zipWith go ps (group cnts)
       where
-        go p (c, d, e, f) = binom d c (min 1 (p * fromIntegral (1+e-f) / fromIntegral (1+f)))
+        go p (c, d, e, f) = binom d c (min (1 - 1e-10) (p * fromIntegral (1 + e - f) / fromIntegral (1 + f)))
+
+    eps = $(TH.lift (2 * until ((== 1) . (1 +)) (/ 2) (1 :: Double)))
 
 -- Tabulate list
 tabulate xs = M.elems $ M.fromListWith (+) [(c, 1) | c <- xs]