store latex tables

cspollard · Apr 4, 2019 · 909daa6 · 909daa6
1 parent bee4edf
commit 909daa6
Show file tree

Hide file tree

Showing 2 changed files with 125 additions and 18 deletions.
diff --git a/app/Main.hs b/app/Main.hs
@@ -4,7 +4,7 @@
 
 module Main where
 
-import           Control.Monad        (void)
+import           Control.Lens
 import           Data.Aeson
 import           Data.Aeson.Types     (Parser, parseEither)
 import qualified Data.ByteString.Lazy as BS
@@ -16,13 +16,15 @@ import           Model
 import           Options.Applicative  hiding (Parser, auto)
 import qualified Options.Applicative  as OA
 import           RunModel
-import           System.IO            (BufferMode (..), hSetBuffering, stdout)
+import           System.IO            (BufferMode (..), IOMode (..), hPutStrLn,
+                                       hSetBuffering, stdout, withFile)
 
 data InArgs =
   InArgs
-    { nsamps  :: Int
-    , outfile :: String
-    , infile  :: String
+    { nsamps     :: Int
+    , outfile    :: String
+    , tablesfile :: String
+    , infile     :: String
     }
 
 inArgs :: OA.Parser InArgs
@@ -34,7 +36,11 @@ inArgs =
     )
   <*> strOption
     ( long "outfile"
-    <> help "text file to record to"
+    <> help "text file to record samples to"
+    )
+  <*> strOption
+    ( long "tablesfile"
+    <> help "text file to write tables to"
     )
   <*> strOption
     ( long "infile"
@@ -44,6 +50,7 @@ inArgs =
 opts :: ParserInfo InArgs
 opts = info (helper <*> inArgs) fullDesc
 
+
 main :: IO ()
 main = do
   InArgs {..} <- execParser opts
@@ -56,10 +63,26 @@ main = do
 
   -- then try to parsse to our data, Model, and ModelParams
   -- NB: need to give explicit types here so the parser knows what to look for.
-  case parseEither parseModel =<< values of
-    Left err -> error err
-    Right (dataH, model, modelparams)
-      -> void $ runModel nsamps outfile dataH model modelparams
+  (params, covariances) <-
+    case parseEither parseModel =<< values of
+      Left err -> error err
+      Right (dataH, model, modelparams)
+        -> runModel nsamps outfile dataH model modelparams
+
+  let uncerts = posteriorMatrices params covariances
+
+  withFile tablesfile WriteMode $ \h -> do
+      hPutStrLn h "absolute uncertainties:"
+      hPutStrLn h . latextable $ view _1 <$> uncerts
+
+      hPutStrLn h ""
+      hPutStrLn h "relative uncertainties:"
+      hPutStrLn h . latextable $ view _2 <$> uncerts
+
+      hPutStrLn h ""
+      hPutStrLn h "correlations:"
+      hPutStrLn h . latextable $ view _3 <$> uncerts
+
 
 
 parseModel

diff --git a/src/RunModel.hs b/src/RunModel.hs
@@ -10,7 +10,7 @@
 {-# LANGUAGE TypeFamilies              #-}
 
 module RunModel
-  ( runModel, X.TDigest
+  ( runModel, X.TDigest, latextable, posteriorMatrices
   ) where
 
 import           Control.Foldl                 (FoldM (..))
@@ -19,7 +19,9 @@ import           Control.Lens
 import           Control.Monad                 (when)
 import           Data.Foldable                 (fold)
 import qualified Data.HashMap.Strict           as M
-import           Data.List                     (intersperse)
+import           Data.List                     (intercalate, intersperse, nub,
+                                                sortBy)
+import           Data.Maybe                    (fromMaybe)
 import           Data.Monoid                   ((<>))
 import           Data.Reflection               (Reifies)
 import           Data.TDigest                  as X
@@ -42,7 +44,7 @@ import           System.IO                     (BufferMode (..), IOMode (..),
                                                 hFlush, hPutStrLn,
                                                 hSetBuffering, stdout, withFile)
 import           System.Random.MWC.Probability
-import           Text.Printf                   (printf)
+import           Text.Printf                   (PrintfArg, printf)
 
 toError :: Either String c -> c
 toError = either error id
@@ -139,7 +141,7 @@ runModel nsamps outfile dataH model' modelparams = do
       reluncerts = zipWith (\x v -> (/x) <$> v) (V.toList start') absuncerts
 
       mpnameslist = V.toList mpnames
-      matToTable = latextable mpnameslist
+      matToTable = latextable' mpnameslist
 
   putStrLn "prediction given starting params:"
   print predstart
@@ -233,8 +235,8 @@ runModel nsamps outfile dataH model' modelparams = do
         tdigestf = F.Fold (flip insert) mempty id
 
         pairwise :: [a] -> [(a, a)]
-        pairwise xs@(x:ys) = fmap (x,) xs ++ pairwise ys
-        pairwise _         = []
+        pairwise ys@(y:ys') = fmap (y,) ys ++ pairwise ys'
+        pairwise _          = []
 
         vcovf n =
           F.premap (V.fromList . pairwise . V.toList)
@@ -273,6 +275,88 @@ runModel nsamps outfile dataH model' modelparams = do
     return (hmval, hmcov)
 
 
+
+
+posteriorMatrices
+  :: Ord a
+  => M.HashMap T.Text (Maybe a, TDigest comp)
+  -> M.HashMap (T.Text, T.Text) Double
+  -> M.HashMap (T.Text, T.Text) (Double, Double, Double)
+posteriorMatrices params covariances =
+  let params' =
+        fromMaybe (error "error getting mode or quantiles") . sequence
+        $ quant <$> params
+
+      quant (mx, y) = do
+        x <- mx
+        q16 <- quantile 0.16 y
+        q50 <- quantile 0.50 y
+        q84 <- quantile 0.84 y
+        return (x, (q16, q50, q84))
+
+      symmetrize hm =
+        let insert1 h (x, y) = M.insert (y, x) (h M.! (x, y)) h
+        in foldl insert1 hm $ M.keys hm
+
+
+      uncerts =
+        flip M.mapMaybeWithKey (symmetrize covariances) $ \(name, name') cov ->
+          let var =
+                fromMaybe (error "missing variance")
+                $ M.lookup (name, name) covariances
+
+              var' =
+                fromMaybe (error "missing variance")
+                $ M.lookup (name', name') covariances
+
+              mean1 =
+                fromMaybe (error "missing best fit value") $ do
+                  (_, (_, q50, _)) <- M.lookup name params'
+                  return q50
+
+              corr = cov / sqrt var / sqrt var'
+              absuncert = abs $ cov / sqrt var'
+              reluncert = absuncert / mean1
+
+
+          in
+            if T.isPrefixOf "normtruthbin" name
+                && not (T.isPrefixOf "truthbin" name')
+                && not (T.isPrefixOf "recobin" name')
+                && name' /= "llh"
+              then Just (absuncert, reluncert, corr)
+              else Nothing
+
+  in uncerts
+
+
+latextable :: PrintfArg a => M.HashMap (T.Text, T.Text) a -> String
+latextable m =
+  let ks = M.keys m
+      srt s s' = if T.isPrefixOf "normtruthbin" s then LT else s `compare` s'
+
+      poinames = sortBy srt . nub $ fst <$> ks
+      npnames = sortBy srt . nub $ snd <$> ks
+      fmtLine npname =
+        T.unpack (paramToName npname)
+        ++ " & "
+        ++ intercalate " & "
+            ( printf "%.3f" . (M.!) m . (,npname)
+              <$> poinames
+            )
+        ++ " \\\\"
+
+  in unlines $
+    [ "\\begin{tabular}{ l " ++ fold (replicate (length poinames) "| r ") ++ "}"
+    , " & " ++ intercalate " & " (T.unpack . paramToName <$> poinames) ++ " \\\\"
+    , "\\hline"
+    ]
+    ++ (fmtLine <$> npnames)
+    ++ ["\\end{tabular}"]
+
+
+
+
 -- taken directly from the `ad` package, but now testing for NaNs.
 
 -- | The 'gradientDescent' function performs a multivariate
@@ -312,8 +396,8 @@ gradientAscent'
 gradientAscent' f = gradientDescent' (negate . f)
 
 
-latextable :: [T.Text] -> [[Double]] -> T.Text
-latextable names mat =
+latextable' :: [T.Text] -> [[Double]] -> T.Text
+latextable' names mat =
   let fmtLine npname vals =
         paramToName npname
         <> " & "