text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE DeriveGeneric       #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE OverloadedStrings   #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeOperators       #-}
{-# LANGUAGE QuasiQuotes         #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE UndecidableInstances      #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE AllowAmbiguousTypes       #-}
{-# LANGUAGE UndecidableSuperClasses #-}
module Frames.Regression where

import qualified Frames.Misc                  as FU
import qualified Frames.VegaLite.Utils         as FV
import qualified Math.Regression.Regression    as MR
import qualified Math.Regression.LeastSquares  as MR
import qualified Polysemy                      as P
import qualified Knit.Effect.Logger            as Log

import qualified Lucid                         as LH
import qualified Text.Blaze.Html               as BH
import qualified Control.Foldl                 as FL
import qualified Data.List                     as List
import qualified Data.Text                     as T
import qualified Data.Vinyl                    as V
import qualified Data.Vinyl.Functor            as V
import qualified Data.Vinyl.TypeLevel          as V
import           Data.Profunctor               as PF
import qualified Frames                        as F
import qualified Frames.Melt                   as F
import qualified Data.Vector.Storable          as VS

import qualified MachineLearning               as ML
import qualified MachineLearning.Regression    as ML
import qualified Numeric.LinearAlgebra         as LA
import           Numeric.LinearAlgebra.Data     ( R )
import qualified Statistics.Types              as S

import           GHC.TypeLits                   ( Symbol )

type Unweighted = "unweighted" F.:-> Void

type family NonVoidField (f :: (Symbol, Type)) :: Bool where
  NonVoidField '(a,Void) = 'False
  NonVoidField _         = 'True

class BoolVal (w :: Bool) where
  asBool :: Bool

instance BoolVal 'True where
  asBool = True

instance BoolVal 'False where
  asBool = False

data FrameRegressionResult (y :: (Symbol, Type)) (wc :: Bool) (xs :: [(Symbol, Type)]) (w :: (Symbol,Type)) (rs :: [(Symbol, Type)]) where
  FrameUnweightedRegressionResult :: (w ~ Unweighted)
    => MR.RegressionResult R -> FrameRegressionResult y wc xs Unweighted rs
  FrameWeightedRegressionResult :: (V.KnownField w, F.ElemOf rs w)
    => (V.Snd w -> R) -> MR.RegressionResult R -> FrameRegressionResult y wc xs w rs

realRecToList
  :: ( RealFrac b
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => F.Record as
  -> [b]
realRecToList =
  V.recordToList
    . V.rmap (\(V.Compose (V.Dict x)) -> V.Const $ realToFrac x)
    . V.reifyConstraint @Real

instance ( MR.Predictor S.NormalErr a (LA.Vector a) (MR.RegressionResult R)
         , BoolVal wc
         , LA.Element a
         , LA.Numeric a
         , RealFloat a
         , xs F.⊆ rs
         , V.RMap xs
         , V.RecordToList xs
         , V.ReifyConstraint Real F.ElField xs
         , V.NatToInt (V.RLength xs)
         ) => MR.Predictor S.NormalErr a (F.Record rs) (FrameRegressionResult y wc xs w as) where
  predict frr r =
    let rr = regressionResult frr
        addBias l = if asBool @wc then 1 : l else l -- if we regressed with a constant term, we need to put it into the xs for prediction
        va :: LA.Vector a = VS.fromList $ addBias $ realRecToList (F.rcast @xs r)
    in MR.predict rr va

regressionResult :: FrameRegressionResult y wc xs w rs -> MR.RegressionResult R
regressionResult (FrameUnweightedRegressionResult x) = x
regressionResult (FrameWeightedRegressionResult _ x) = x

withConstant
  :: forall y wc xs w rs
   . (BoolVal wc)
  => FrameRegressionResult y wc xs w rs
  -> Bool
withConstant _ = asBool @wc

weightedRegression
  :: forall y wc xs w rs
   . (BoolVal (NonVoidField w))
  => FrameRegressionResult y wc xs w rs
  -> Bool
weightedRegression _ = asBool @(NonVoidField w)

-- make X, y from the data
prepRegression
  :: forall y as f rs
   . ( Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => f (F.Record rs)
  -> (LA.Matrix R, LA.Vector R)
prepRegression dat =
  let
    nCols  = V.natToInt @(V.RLength as)
    yListF = PF.dimap (realToFrac . F.rgetField @y) List.reverse FL.list
    toListOfDoubles :: F.Record as -> [Double] = realRecToList
    mListF = PF.dimap (toListOfDoubles . F.rcast)
                      (List.concat . List.reverse)
                      FL.list
    (yList, mList) = FL.fold ((,) <$> yListF <*> mListF) dat
  in
    (LA.matrix nCols mList, LA.vector yList)

-- make X, y, w from the data
prepWeightedRegression
  :: forall y as w f rs
   . ( Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , F.ElemOf rs w
     , FU.RealField w
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => f (F.Record rs)
  -> (LA.Matrix R, LA.Vector R, LA.Vector R)
prepWeightedRegression dat =
  let
    nCols  = V.natToInt @(V.RLength as)
    yListF = PF.dimap (realToFrac . F.rgetField @y) List.reverse FL.list
    wListF = PF.dimap (realToFrac . F.rgetField @w) List.reverse FL.list
    toListOfDoubles :: F.Record as -> [Double] = realRecToList
    mListF = PF.dimap (toListOfDoubles . F.rcast)
                      (List.concat . List.reverse)
                      FL.list
    (yList, mList, wList) = FL.fold ((,,) <$> yListF <*> mListF <*> wListF) dat
  in
    (LA.matrix nCols mList, LA.vector yList, LA.vector wList)

prettyPrintRegressionResult
  :: forall y wc as w rs
   . ( F.ColumnHeaders '[y]
     , F.ColumnHeaders '[w]
     , F.ColumnHeaders as
     , BoolVal wc
     )
  => (T.Text -> T.Text -> T.Text)
  -> FrameRegressionResult y wc as w rs
  -> S.CL Double
  -> T.Text
prettyPrintRegressionResult headerF res cl =
  let yName   = FV.colName @y
      wName   = FV.colName @w
      xNames' = T.pack <$> F.columnHeaders (Proxy :: Proxy (F.Record as))
      xNames  = if withConstant res then "intercept" : xNames' else xNames'
  in  MR.prettyPrintRegressionResult (headerF yName wName)
                                     xNames
                                     (regressionResult res)
                                     cl

prettyPrintRegressionResultLucid
  :: forall y wc as w rs
   . ( F.ColumnHeaders '[y]
     , F.ColumnHeaders '[w]
     , F.ColumnHeaders as
     , BoolVal wc
     )
  => (T.Text -> T.Text -> T.Text)
  -> FrameRegressionResult y wc as w rs
  -> S.CL Double
  -> LH.Html ()
prettyPrintRegressionResultLucid headerF res cl =
  let yName   = FV.colName @y
      wName   = FV.colName @w
      xNames' = T.pack <$> F.columnHeaders (Proxy :: Proxy (F.Record as))
      xNames  = if withConstant res then "intercept" : xNames' else xNames'
  in  MR.prettyPrintRegressionResultLucid (headerF yName wName)
                                          xNames
                                          (regressionResult res)
                                          cl

prettyPrintRegressionResultBlaze
  :: forall y wc as w rs
   . ( F.ColumnHeaders '[y]
     , F.ColumnHeaders '[w]
     , F.ColumnHeaders as
     , BoolVal wc
     )
  => (T.Text -> T.Text -> T.Text)
  -> FrameRegressionResult y wc as w rs
  -> S.CL Double
  -> BH.Html
prettyPrintRegressionResultBlaze headerF res cl =
  let yName   = FV.colName @y
      wName   = FV.colName @w
      xNames' = T.pack <$> F.columnHeaders (Proxy :: Proxy (F.Record as))
      xNames  = if withConstant res then "intercept" : xNames' else xNames'
  in  MR.prettyPrintRegressionResultBlaze (headerF yName wName)
                                          xNames
                                          (regressionResult res)
                                          cl

keyRecordText
  :: (V.ReifyConstraint Show F.ElField ks, V.RecordToList ks, V.RMap ks)
  => F.Record ks
  -> T.Text
keyRecordText keyRec =
  let keyValuesAsText rks =
        V.recordToList
          . V.rmap (\(V.Compose (V.Dict x)) -> V.Const $ T.pack $ show x)
          $ V.reifyConstraint @Show rks
  in  T.intercalate "; " $ keyValuesAsText keyRec

prettyPrintRegressionResults
  :: forall y wc as w rs k f a
   . ( F.ColumnHeaders '[y]
     , F.ColumnHeaders '[w]
     , F.ColumnHeaders as
     , BoolVal wc
     , BoolVal (NonVoidField w)
     , Foldable f
     , Monoid a
     )
  => (k -> T.Text)
  -> f (k, FrameRegressionResult y wc as w rs)
  -> S.CL R
  -> (  (T.Text -> T.Text -> T.Text)
     -> FrameRegressionResult y wc as w rs
     -> S.CL Double
     -> a
     )
  -> a
  -> a
prettyPrintRegressionResults keyText keyed cl printOne sepEach =
  let headerF res key yName wName = if weightedRegression res
                                    then
                                      "Explaining "
                                        <> yName
                                        <> " ("
                                        <> keyText key
                                        <> "; weights from "
                                        <> wName
                                        <> ")"
                                    else "Explaining " <> yName <> " (" <> keyText key <> ")"
  in  FL.fold
        (FL.Fold (\t (rk, res) -> t <> printOne (headerF res rk) res cl)
                 sepEach
                 id
        )
        keyed

-- explain y in terms of as
leastSquaresByMinimization
  :: forall y as rs f
   . ( Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => Bool
  -> [R]
  -> f (F.Record rs)
  -> [R]
leastSquaresByMinimization wc guess dat =
  let (mX, y)       = prepRegression @y @as dat
      mX1           = if wc then ML.addBiasDimension mX else mX
      (solution, _) = ML.minimize (ML.ConjugateGradientFR 0.1 0.1)
                                  ML.LeastSquares
                                  0.001
                                  20
                                  ML.RegNone
                                  mX1
                                  y
                                  (LA.fromList guess)
  in  LA.toList solution

ordinaryLeastSquares
  :: forall effs y wc as rs f
   . ( Log.LogWithPrefixesLE effs
     , Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , BoolVal wc
     , V.NatToInt (V.RLength as)
     )
  => f (F.Record rs)
  -> P.Sem effs (FrameRegressionResult y wc as Unweighted rs)
ordinaryLeastSquares dat = do
  let (mA, vB)  = prepRegression @y @as dat
      withConst = asBool @wc
  FrameUnweightedRegressionResult <$> MR.ordinaryLS withConst mA vB

weightedLeastSquares
  :: forall y wc as w rs f effs
   . ( Log.LogWithPrefixesLE effs
     , Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , F.ElemOf rs w
     , FU.RealField w
     , BoolVal wc
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => f (F.Record rs)
  -> P.Sem effs (FrameRegressionResult y wc as w rs)
weightedLeastSquares dat = do
  let (mA, vB, vW) = prepWeightedRegression @y @as @w dat
      withConst    = asBool @wc
  FrameWeightedRegressionResult realToFrac <$> MR.weightedLS withConst mA vB vW

-- special case when weights come from observations being population averages of different populations
popWeightedLeastSquares
  :: forall effs y wc as w rs f
   . ( Log.LogWithPrefixesLE effs
     , Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , F.ElemOf rs w
     , FU.RealField w
     , BoolVal wc
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => f (F.Record rs)
  -> P.Sem effs (FrameRegressionResult y wc as w rs)
popWeightedLeastSquares dat = do
  let (mA, vB, vW) = prepWeightedRegression @y @as @w dat
      withConst    = asBool @wc
      vWpop        = LA.cmap sqrt vW -- this is the correct weight for population average, the sqrt of the number averaged in that sample
  FrameWeightedRegressionResult (sqrt . realToFrac)
    <$> MR.weightedLS withConst mA vB vWpop

-- special case when we know residuals are heteroscedastic with variances proportional to given numbers
varWeightedLeastSquares
  :: forall effs y wc as w rs f
   . ( Log.LogWithPrefixesLE effs
     , Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , F.ElemOf rs w
     , FU.RealField w
     , BoolVal wc
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => f (F.Record rs)
  -> P.Sem effs (FrameRegressionResult y wc as w rs)
varWeightedLeastSquares dat = do
  let (mA, vB, vW) = prepWeightedRegression @y @as @w dat
      withConst    = asBool @wc
      vWvar        = LA.cmap (\x -> 1 / sqrt x) vW -- this is the correct weight for given variance
  FrameWeightedRegressionResult (\x -> 1 / sqrt (realToFrac x))
    <$> MR.weightedLS withConst mA vB vWvar

totalLeastSquares
  :: forall effs y wc as rs f
   . ( Log.LogWithPrefixesLE effs
     , Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , BoolVal wc
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => f (F.Record rs)
  -> P.Sem effs (FrameRegressionResult y wc as Unweighted rs)
totalLeastSquares dat = do
  let (mA, vB)  = prepRegression @y @as dat
      withConst = asBool @wc
  FrameUnweightedRegressionResult <$> MR.totalLS withConst mA vB

weightedTLS
  :: forall effs y wc as w rs f
   . ( Log.LogWithPrefixesLE effs
     , Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , BoolVal wc
     , F.ElemOf rs w
     , FU.RealField w
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => f (F.Record rs)
  -> P.Sem effs (FrameRegressionResult y wc as w rs)
weightedTLS dat = do
  let (mA, vB, vW) = prepWeightedRegression @y @as @w dat
      withConst    = asBool @wc
  FrameWeightedRegressionResult realToFrac <$> MR.weightedTLS withConst mA vB vW

popWeightedTLS
  :: forall effs y wc as w rs f
   . ( Log.LogWithPrefixesLE effs
     , Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , BoolVal wc
     , F.ElemOf rs w
     , FU.RealField w
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => f (F.Record rs)
  -> P.Sem effs (FrameRegressionResult y wc as w rs)
popWeightedTLS dat = do
  let (mA, vB, vW) = prepWeightedRegression @y @as @w dat
      withConst    = asBool @wc
      vWpop        = LA.cmap sqrt vW -- this is the correct weight for population average, the sqrt of the number averaged in that sample
  FrameWeightedRegressionResult (sqrt . realToFrac)
    <$> MR.weightedTLS withConst mA vB vWpop

varWeightedTLS
  :: forall effs y wc as w rs f
   . ( Log.LogWithPrefixesLE effs
     , Foldable f
     , as F.⊆ rs
     , F.ElemOf rs y
     , FU.RealField y
     , BoolVal wc
     , F.ElemOf rs w
     , FU.RealField w
     , V.AllConstrained (FU.RealFieldOf rs) as
     , V.RMap as
     , V.RecordToList as
     , V.ReifyConstraint Real F.ElField as
     , V.NatToInt (V.RLength as)
     )
  => f (F.Record rs)
  -> P.Sem effs (FrameRegressionResult y wc as w rs)
varWeightedTLS dat = do
  let (mA, vB, vW) = prepWeightedRegression @y @as @w dat
      withConst    = asBool @wc
      vWvar        = LA.cmap (\x -> 1 / sqrt x) vW -- this is the correct weight for given variance
  FrameWeightedRegressionResult (\x -> 1 / sqrt (realToFrac x))
    <$> MR.weightedTLS withConst mA vB vWvar

{-
-- this is sort of ugly but I think it will work.  And it has a pretty narrow purpose
-- But what will we do when we want logistic regression or whatever?
class (V.AllConstrained RealFrac (V.Unlabeled rs), Real a) => RealFracRecordFromList a rs where
  recordFromList :: [a] -> Maybe (F.Record rs)

instance Real a => RealFracRecordFromList a '[] where
  recordFromList _ = Just $ V.RNil

instance ( V.AllConstrained RealFrac (V.Unlabeled (r : rs))
         , V.AllConstrained RealFrac (V.Unlabeled rs)
         , RealFracRecordFromList a rs
         , V.KnownField r
         , RealFrac (V.Snd r)
         , Real a) => RealFracRecordFromList a (r : rs) where
  recordFromList [] = Nothing
  recordFromList (a : as) = case recordFromList as of
    Nothing -> Nothing
    Just xs -> let x = (realToFrac a) :: V.Snd r in Just $ x &: xs
-}

max_stars_repo_name: teto/Frames-utils

max_issues_repo_name: teto/Frames-utils

max_forks_repo_name: teto/Frames-utils


text:
{-# LANGUAGE AllowAmbiguousTypes, CPP, ConstraintKinds, DataKinds, FlexibleContexts, FlexibleInstances, GADTs, MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedLabels, OverloadedStrings, PartialTypeSignatures, PolyKinds, QuasiQuotes, RankNTypes, ScopedTypeVariables     #-}
{-# LANGUAGE TypeApplications, TypeFamilyDependencies, TypeInType, TypeOperators, UndecidableInstances                               #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -Wno-partial-type-signatures -fconstraint-solver-iterations=10 #-}

{-| Talk to Tensorboard. You will need to run it manually with
'tensorboard --logdir <yourdir>'

It will occasionally update on its own.  We have fairly comprehensive
Tensorboard support and integration with showing images, plots, and graphs.
-}

module Torch.Tensorboard where
import           Control.Monad
import           Control.Monad.Logger
import           Data.Binary.Put
import           Data.Bits
import qualified Data.ByteString              as B
import qualified Data.ByteString              as BS
import qualified Data.ByteString.Internal     as BS
import qualified Data.ByteString.Lazy         as BL
import           Data.Coerce
import           Data.IORef
import qualified Data.Map                     as M
import           Data.ProtoLens
import           Data.ProtoLens.Labels        ()
import           Data.Singletons
import           Data.Singletons.TypeLits
import           Data.Text                    (Text)
import qualified Data.Text                    as T
import qualified Data.Text.Encoding           as T
import           Data.Time.Clock.POSIX
import           Data.Time.Format
import           Data.Time.LocalTime
import           Data.Vector.Storable         (Vector)
import qualified Data.Vector.Storable         as V
import           Foreign.C.Types
import           Foreign.ForeignPtr
import           Foreign.Ptr
import           GHC.Float
import           GHC.Int
import           GHC.Word
import qualified Graphics.Matplotlib          as M
import           Lens.Micro
import           Network.HostName
import qualified Statistics.Quantile          as S
import           System.Directory
import           System.FilePath
import           System.IO
import qualified Tensorboard.Proto.Attributes as PA
import qualified Tensorboard.Proto.Event      as PE
import qualified Tensorboard.Proto.Graph      as PG
import qualified Tensorboard.Proto.Summary    as PS
import qualified Tensorboard.Proto.Tensor     as PT
import qualified Torch.C.Tensor               as C
import qualified Torch.C.Types                as C
import qualified Torch.C.Variable             as C
import           Torch.Images
import           Torch.Internal.CRC32C
import           Torch.Tensor                 as Tensor
import           Torch.Types                  as Tensor
import           Torch.Visualization

-- | `SummaryWriter` keeps track of some basic information so we can talk to
-- Tensorboard. You'll be passing this around everywhere.
data SummaryWriter = SummaryWriter { runName      :: Text
                                   , runDirectory :: Text
                                   , runFilename  :: Text
                                   , lastStep     :: IORef Int
                                   }

-- | Make a new summary writer for a given directory and run name. We will also
-- add some information to the run time like the current host and local timestamp.
summaryWriter :: Text -> Text -> IO SummaryWriter
summaryWriter directory runname = do
  now <- getCurrentTime
  tz <- getCurrentTimeZone
  let time = formatTime defaultTimeLocale "%Y-%m-%d_%H:%M:%S" (utcToZonedTime tz now)
  let timeUnix = formatTime defaultTimeLocale "%s" (utcToZonedTime tz now)
  host <- T.pack <$> getHostName
  let name = runname <> "_" <> host <> "_" <> T.pack time
  let location = T.unpack directory </> T.unpack name
  createDirectoryIfMissing True location
  step <- newIORef 0
  pure $ SummaryWriter { runName = name
                       , runDirectory = T.pack location
                       , runFilename = "events.haskelltorch." <> T.pack timeUnix <> ".tfevents." <> host
                       , lastStep = step
                       }

-- | Advance the writer to the next step. This helps you keep track of the
-- iteration number without needing a manual counter.
nextStep :: SummaryWriter -> IO SummaryWriter
nextStep r = do
  modifyIORef' (lastStep r) (+ 1)
  pure r

-- | Explicitly set the step counter.
setStep :: SummaryWriter -> Int -> IO SummaryWriter
setStep sw s = do
  writeIORef (lastStep sw) s
  pure sw

-- | Add a scalar to Tensorboard with a given name and value.
addScalar :: Real a => SummaryWriter -> Text -> a -> IO ()
addScalar summary name value =
  writeEvent summary (EventSummary name $ SummaryValue $ fromRational $ toRational value)

-- | Add a histogram to Tensorboard, we do some internal work to make sure the
-- number of bins and their sizes make sense.
addHistogram ::
  (TensorConstraints ty 'KCpu sz) =>
  SummaryWriter -> Text -> Tensor ty 'KCpu sz -> IO ()
addHistogram summary name tensor = do
  bins  <- histogramBinsAuto <$> V.map toDoubleC <$> toVector tensor
  withSomeSing (fromIntegral bins) (\ n@(SNat :: SNat m) -> do
    hist_ <- toVector =<< histc @m tensor 0 0
    let hist = V.map toDoubleC hist_
    min   <- toDouble <$> (fromScalar =<< Tensor.min tensor)
    max   <- toDouble <$> (fromScalar =<< Tensor.max tensor)
    let step = (max-min)/2
    let binLimits = V.generate bins (\i -> step * fromIntegral i)
    sum     <- fromScalar =<< Tensor.sum tensor
    sumSq <- fromScalar =<< Tensor.sum =<< Tensor.pow tensor =<< toScalar 2
    writeEvent summary
      (EventSummary name
       (SummaryHistogram
         { shMin = min
         , shMax = max
         , shNum = fromIntegral $ numel tensor
         , shSum = toDouble sum
         , shSumSq = toDouble sumSq
         , shBucketLimit = binLimits
         , shBucket = hist })))

-- | Add an RGB image tensor to the Tensorboard.
addImageRGB :: (KnownNat szw, KnownNat szh, SingI ty)
            => SummaryWriter -> Text -> Tensor ty 'KCpu '[3, szh, szw] -> IO ()
addImageRGB summary name image = do
    (bytes, len) <- V.unsafeToForeignPtr0 <$> rgbTensorToVector image
    (h,w,c,b) <- pure $ case V.toList $ size image of
       [3,h,w] -> (h,w,3,bytes)
       -- TODO Move this to compile time
       _       -> error $ "Don't know how to display a tensor with these dimensions as an RGB image " ++ show (size image)
    writeEvent summary
      (EventSummary name
       (SummaryImage
        { siHeight     = fromIntegral h
        , siWidth      = fromIntegral w
        , siColorspace = c
        , siBytes      = BS.PS b 0 len}))

-- | Add a greyscale image tensor to the Tensorboard.
addImageGrey :: (KnownNat szw, KnownNat szh, SingI ty)
             => SummaryWriter -> Text -> Tensor ty 'KCpu '[1, szh, szw] -> IO ()
addImageGrey summary name image = do
    (bytes, len) <- V.unsafeToForeignPtr0 <$> greyTensorToVector image
    (h,w,c,b) <- pure $ case V.toList $ size image of
       [1,h,w] -> (h,w,1,bytes)
       -- TODO Move this to compile time
       _       -> error $ "Don't know how to display a tensor with these dimensions as an Greyscale image " ++ show (size image)
    writeEvent summary
      (EventSummary name
       (SummaryImage
        { siHeight     = fromIntegral h
        , siWidth      = fromIntegral w
        , siColorspace = c
        , siBytes      = BS.PS b 0 len}))

-- | Add a tensor that contains a grid of RGB images, like a batch of RGB images
-- with the leading dimension being the number of images. Grids are padded with
-- the given Size, constructed with the @Size@ datatype. They are padded with a
-- given constant, usually 0.
addImageRGBGrid ::
  forall (imagesPerRow :: Nat) (padding :: Nat) (rows :: Nat) (szh :: Nat) (szw :: Nat) (nr :: Nat) (ty :: TensorType).
    _
  => SummaryWriter
  -> Text
  -> Size imagesPerRow
  -> Padding padding
  -> TensorTyToHs ty
  -> Tensor ty 'KCpu '[nr, 3, szh, szw]
  -> IO ()
addImageRGBGrid summary name sz padding padValue images = do
  grid <- makeRGBGrid sz padding padValue images
  addImageRGB summary name grid

-- | Add a matplotlib plot to Tensorboard. This uses the Haskell matplotlib
-- bindings from the `matplotlib` package,
-- `https://github.com/abarbu/matplotlib-haskell`.
addPlot :: SummaryWriter -> Text -> M.Matplotlib -> IO ()
addPlot summary name mplot = do
  (fname, handle) <- openTempFile "/tmp/" "plot.png"
  hClose handle
  M.file fname mplot
  (h,w,_) <- fileImageProperties (T.pack fname)
  withSomeSing (fromIntegral h)
    (\ szh@(SNat :: SNat szh) ->
        withSomeSing (fromIntegral w)
          (\ szw@(SNat :: SNat szw) -> do
              tensor <- typed @TByte <$> readRGBTensorFromFile @szh @szw (T.pack fname)
              addImageRGB summary name tensor))
  fileExists <- doesFileExist fname
  when fileExists (removeFile fname)

-- | Add the trace of a computation to Tensorboard as a graph.
-- TODO Support naming nodes, more type information, more device information
addGraph :: SummaryWriter -> Text -> ForeignPtr C.CTracingState -> IO ()
addGraph summary name trace = do
  (inputs, nodes, outputs, ret, block) <-
    withForeignPtr trace C.tracing_state_graph
  inputs' <- mapM (\i -> do
                    b <- cbool <$> C.check_value_tensor i
                    na <- T.pack <$> C.value_name i
                    if b then do
                      sz <- C.value_sizes i
                      ty <- C.value_scalar_type i
                      pure (na, Just sz, Just ty)
                      else
                      pure (na, Nothing, Nothing))
            $ V.toList inputs
  nodes' <- mapM (\(n,nr) -> do
                   ins <- C.node_inputs n
                   outs <- C.node_outputs n
                   kind <- C.node_kind n
                   b <- C.node_has_attribute n "value"
                   attr <- if cbool b then do
                            ak <- C.node_attribute_kind n "value"
                            case ak of
                              C.AttributeKindFloat -> do
                                f <- C.node_get_attribute_float n "value"
                                pure $ Just $ defMessage & #f .~ coerce f
                              C.AttributeKindInt -> do
                                i <- C.node_get_attribute_int n "value"
                                pure $ Just $ defMessage & #i .~ fromIntegral i
                              C.AttributeKindString -> do
                                s <- C.node_get_attribute_string n "value"
                                pure $ Just $ defMessage & #s .~ T.encodeUtf8 (T.pack s)
                              C.AttributeKindTensor -> do
                                tensor <- C.node_get_attribute_tensor n "value"
                                ty <- C.getType tensor
                                tensor' <- C.contiguous_mm tensor (fromIntegral $ fromEnum C.MemoryFormatContiguous)
                                ptr <- newForeignPtr_ =<< castPtr <$> C.data_ptr tensor'
                                len <- fromIntegral <$> C.numel tensor'
                                let get :: V.Storable a => [a]
                                    get = V.toList $ V.unsafeFromForeignPtr0 (castForeignPtr ptr) len
                                pure $ Just $ defMessage & #tensor .~
                                  case ty of
                                    C.ScalarTypeByte ->
                                      defMessage & #dtype .~ PS.DT_UINT8
                                                 & #intVal .~ map (fromIntegral :: Word8 -> Int32) get
                                    C.ScalarTypeChar ->
                                      defMessage & #dtype .~ PS.DT_INT8
                                                 & #intVal .~ map (fromIntegral :: Int8 -> Int32) get
                                    C.ScalarTypeShort ->
                                      defMessage & #dtype .~ PS.DT_INT16
                                                 & #intVal .~ map (fromIntegral :: Int16 -> Int32) get
                                    C.ScalarTypeInt ->
                                      defMessage & #dtype .~ PS.DT_INT32
                                                 & #intVal .~ get
                                    C.ScalarTypeLong ->
                                      defMessage & #dtype .~ PS.DT_INT64
                                                 & #int64Val .~ get
                                    C.ScalarTypeHalf ->
                                      defMessage & #dtype .~ PS.DT_HALF
                                                 & #halfVal .~ get
                                    C.ScalarTypeFloat ->
                                      defMessage & #dtype .~ PS.DT_FLOAT
                                                 & #floatVal .~ get
                                    C.ScalarTypeDouble ->
                                      defMessage & #dtype .~ PS.DT_DOUBLE
                                                 & #doubleVal .~ get
                                    C.ScalarTypeUndefined ->
                                      defMessage & #dtype .~ PS.DT_INVALID
                              _ -> pure Nothing
                     else pure Nothing
                   ins' <- mapM (\v -> do
                                n <- C.value_name v
                                t <- C.value_type v
                                ct <- convertCType t
                                pure $ TempGraphValue v n t ct)
                          $ V.toList ins
                   outs' <- mapM (\v -> do
                                n <- C.value_name v
                                t <- C.value_type v
                                ct <- convertCType t
                                pure $ TempGraphValue v n t ct)
                           $ V.toList outs
                   pure (n, TempGraphNode nr kind ins' outs' attr)) $ zip (V.toList nodes) [0..]
  let nodesMap = M.fromList nodes'
  let insToNodes  = M.fromListWith (++)
                  $ concatMap (\(n, TempGraphNode nr kinds ins outs attr) ->
                             map (\(TempGraphValue v _ _ _) -> (v,[n])) ins) nodes'
  let outsToNodes = M.fromListWith (++)
                  $ concatMap (\(n, TempGraphNode nr kinds ins outs attr) ->
                             map (\(TempGraphValue v _ _ _) -> (v,[n])) outs) nodes'
  let nodeName nr kind = T.replace "prim::" "" $ T.replace "aten::" "" $ T.pack kind <> "_"<> T.pack (show nr)
  writeEvent summary
    (EventGraph name
      (map (\(name, size, ty) ->
               case (size, ty) of
                 (Just sz, Just t) ->
                   defMessage & #name .~ "input_" <> name
                              & #op .~ "Input"
                              & #attr .~ M.fromList [("_output_shapes",
                                                      defMessage & #list .~
                                                       (defMessage & #shape .~
                                                         [defMessage & #dim .~
                                                          map (\x -> defMessage & #size .~ x) (V.toList sz)]))]
                 (_, _) ->
                   defMessage & #name .~ "input_" <> name
                              & #op .~ "Input")
         inputs'
       ++
       map (\(node,(TempGraphNode nr kind ins outs attr)) ->
               defMessage & #name .~ nodeName nr kind
                          & #op .~ T.pack kind
                          & #device .~ T.pack ""
                          & #input .~ map (\(TempGraphValue v n t ct) ->
                                                case M.lookup v outsToNodes of
                                                  Nothing -> T.pack ("input_" ++ n)
                                                  Just [outNode] -> case M.lookup outNode nodesMap of
                                                                   Just (TempGraphNode nr kind _ _ _) -> nodeName nr kind
                                                                   Nothing                            -> error "Can't find node") ins
                          & #attr .~ M.fromList([("_output_shapes",
                                                   defMessage
                                                   & #list .~ (defMessage
                                                               & #shape .~
                                                                  (map (\(TempGraphValue v n t ct) ->
                                                                           case ct of
                                                                             TraceTypeSimple _ _ ->
                                                                               defMessage & #dim .~ [defMessage & #size .~ 1]
                                                                             TraceTypeTensor _ (Just sz) _ _ ->
                                                                               defMessage & #dim .~ (map (\x -> defMessage & #size .~ x)
                                                                                (V.toList sz))
                                                                             _ -> defMessage)
                                                                    outs)))]
                                                ++ (case attr of
                                                     Nothing -> []
                                                     Just a  -> [("value",a)])))
           nodes'))

-- | A convenience function, add the summary statistics and histograms for a
-- tensor.
addTensorSummary :: (Num (TensorTyToHs ty), Num (TensorTyToHsC ty),
                     V.Storable (TensorTyToHs ty), V.Storable (TensorTyToHsC ty),
                     SingI sz, SingI ty) =>
                   SummaryWriter -> Text -> Tensor ty 'KCpu sz -> IO ()
addTensorSummary summary name tensor = do
  mean   <- fromScalar =<< Tensor.mean tensor
  median <- fromScalar =<< Tensor.median tensor
  var    <- fromScalar =<< Tensor.var tensor True
  min    <- fromScalar =<< Tensor.min tensor
  max    <- fromScalar =<< Tensor.max tensor
  addScalar    summary ("summaries/" <> name <> "/mean")   $ toDouble mean
  addScalar    summary ("summaries/" <> name <> "/median") $ toDouble median
  addScalar    summary ("summaries/" <> name <> "/stddev") $ Prelude.sqrt $ toDouble var
  addScalar    summary ("summaries/" <> name <> "/min")    $ toDouble min
  addScalar    summary ("summaries/" <> name <> "/max")    $ toDouble max
  addHistogram summary ("summaries/" <> name <> "/histogram") tensor

-- * Basic heuristics

-- | The numpy auto histogram mode
histogramBinsAuto :: Vector Double -> Int
histogramBinsAuto x | V.length x < 1000 = histogramBinsSturges x
                    | otherwise = histogramBinsFD x

-- | The Freedman-Diaconis histogram bin estimator.
-- bins = 2 \frac{IQR}{n^{1/3}}
histogramBinsFD :: Vector Double -> Int
histogramBinsFD v =
  ceiling $ (V.maximum v - V.minimum v) /
               (2 * (S.weightedAvg 3 4 v - S.weightedAvg 1 4 v)
                  / (fromIntegral (V.length v) ** (1/3)))

-- | Sturges's histogram bin estimator.
histogramBinsSturges :: Vector Double -> Int
histogramBinsSturges v = Prelude.round $ logBase 2 (fromIntegral $ V.length v) + 1

-- * Internal

-- | This is internal, you never need to use it manually. It's what we fill out
-- with the various functions in this module and then serialize for Tensorboard.
data Summary = SummaryValue Double
             | SummaryHistogram { shMin         :: Double
                                , shMax         :: Double
                                , shNum         :: Double
                                , shSum         :: Double
                                , shSumSq       :: Double
                                , shBucketLimit :: Vector Double
                                , shBucket      :: Vector Double }
             | SummaryImage { siHeight     :: Int32
                            , siWidth      :: Int32
                            , siColorspace :: Int32
                            , siBytes      :: BS.ByteString }

-- | This is internal.
data Node = Node { nName   :: Text
                 , nOp     :: Text
                 , nInput  :: [Text]
                 , nDevice :: Text
                 }

-- | Internal. Used to capture the various kinds of entities that Tensorboard
-- can understand.
data Event = EventVersion
           | EventGraph Text [PG.NodeDef]
           | EventSummary Text Summary
           | EventMessage LogLevel Text
           | EventSessionLog
           | EventMetadata
           | EventMetaGraph

-- | Internal
data TempGraphValue = TempGraphValue (Ptr C.CJitValue) String (Ptr C.CType) TraceType
  deriving (Show)
-- | Internal
data TempGraphNode  = TempGraphNode Int String [TempGraphValue] [TempGraphValue] (Maybe PA.AttrValue)
  deriving (Show)

-- | Internal
eventToProto :: Int -> Event -> IO PE.Event
eventToProto step event = do
  (time :: Double) <- realToFrac <$> getPOSIXTime
  let d = case event of
            EventVersion -> defMessage & #fileVersion .~ "brain.Event:2"
            EventMessage lvl txt ->
              defMessage & #logMessage .~ (defMessage & #level .~ (case lvl of
                                                        LevelDebug         -> PE.LogMessage'DEBUG
                                                        LevelInfo          -> PE.LogMessage'INFO
                                                        LevelWarn          -> PE.LogMessage'WARN
                                                        LevelError         -> PE.LogMessage'ERROR
                                                        LevelOther "fatal" -> PE.LogMessage'FATAL
                                                        LevelOther _       -> PE.LogMessage'UNKNOWN)
                                        & #message .~ txt)
            EventSummary tag (SummaryValue val) ->
              defMessage & #summary .~
              (defMessage & #value .~
                [defMessage & #tag .~ tag
                -- TODO Metadata!
                     & #simpleValue .~ double2Float val])
            EventSummary tag sh@SummaryHistogram {} ->
              defMessage & #summary .~
              (defMessage & #value .~
                [defMessage & #tag .~ tag
                -- TODO Metadata!
                     & #histo .~
                     (defMessage
                          & #min .~ shMin sh
                          & #max .~ shMax sh
                          & #num .~ shNum sh
                          & #sum   .~ shSum sh
                          & #sumSquares .~ shSumSq sh
                          & #bucketLimit .~ V.toList (shBucketLimit sh)
                          & #bucket .~ V.toList (shBucket sh))])
            EventSummary tag sh@SummaryImage{} ->
              defMessage & #summary .~
              (defMessage & #value .~
                [defMessage & #tag .~ tag
                -- TODO Metadata!
                     & #image .~
                     (defMessage
                          & #height             .~ siHeight sh
                          & #width              .~ siWidth sh
                          & #colorspace         .~ siColorspace sh
                          & #encodedImageString .~ siBytes sh)])
            EventGraph tag nodes ->
              defMessage & #graphDef .~
                encodeMessage ((defMessage :: PG.GraphDef)
                                & #versions .~ (defMessage & #producer .~ 1
                                                           & #minConsumer .~ 1
                                                           & #badConsumers .~ [])
                                & #version .~ 3
                                & #node .~ nodes)
  pure $ d & #wallTime .~ time
           & #step .~ fromIntegral step

-- | Internal
encodeEvent :: PE.Event -> B.ByteString
encodeEvent ev = BL.toStrict $ runPut $ do
  putByteString len
  putWord32le (maskedCRC32 len)
  putByteString msg
  putWord32le (maskedCRC32 msg)
  where msg = encodeMessage ev
        len = BL.toStrict $ runPut $ putWord64le (fromIntegral $ B.length msg)

-- | Internal
writeEvent :: SummaryWriter -> Event -> IO ()
writeEvent r e = do
  step <- readIORef $ lastStep r
  let fname = T.unpack (runDirectory r) </> (T.unpack $ runFilename r)
  b <- doesFileExist fname
  unless b $ do
    bs <- encodeEvent <$> eventToProto step EventVersion
    B.writeFile fname bs
  bs <- encodeEvent <$> eventToProto step e
  B.appendFile fname bs

-- | Internal
maskedCRC32 :: BS.ByteString -> Word32
maskedCRC32 bs = (crc32c bs `rotateR` 15) + 0xa282ead8

ext: hs

max_stars_repo_name: abarbu/haskell-torch

max_issues_repo_name: abarbu/haskell-torch

max_forks_repo_name: abarbu/haskell-torch


text:
import Data.Complex

-- Cooley-Tukey
fft [] = []
fft [x] = [x]
fft xs = zipWith (+) ys ts ++ zipWith (-) ys ts
    where n = length xs
          ys = fft evens
          zs = fft odds
          (evens, odds) = split xs
          split [] = ([], [])
          split [x] = ([x], [])
          split (x:y:xs) = (x:xt, y:yt) where (xt, yt) = split xs
          ts = zipWith (\z k -> exp' k n * z) zs [0..]
          exp' k n = cis $ -2 * pi * (fromIntegral k) / (fromIntegral n)

main = mapM_ print $ fft [1,1,1,1,0,0,0,0]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
{-# LANGUAGE DeriveGeneric     #-}
{-# LANGUAGE FlexibleContexts  #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE PatternSynonyms   #-}
{-# LANGUAGE RecordWildCards   #-}

module Main where

import           BDSCOD.Aggregation         (AggregatedObservations (..),
                                             AggregationTimes,
                                             aggregateUnscheduledObservations,
                                             maybeAggregationTimes,
                                             pattern AggTimes)
import           BDSCOD.Llhd                (initLlhdState, llhdAndNB,
                                             logPdeStatistics, pdeGF,
                                             pdeStatistics)
import           BDSCOD.Types               (LogLikelihood (..),
                                             MCMCConfiguration (..), MWCSeed,
                                             NegativeBinomial (..), NumLineages,
                                             Observation (..),
                                             ObservedEvent (..),
                                             PDESolution (..), Parameters (..),
                                             packParameters, putLambda, putMu,
                                             putNus, putOmega, putPsi, putRhos,
                                             scheduledTimes, unpackParameters)
import           BDSCOD.Utility             (eventsAsObservations, invLogit)
import           Control.Monad              (liftM2, when)
import           Control.Monad.Except       (ExceptT, runExceptT, throwError)
import           Control.Monad.Reader       (ReaderT, asks, liftIO, runReaderT)
import qualified Data.Aeson                 as Json
import qualified Data.ByteString.Builder    as BBuilder
import qualified Data.ByteString.Lazy       as L
import           Data.ByteString.Lazy.Char8 (pack, singleton)
import qualified Data.Csv                   as Csv
import           Data.List                  (intercalate)
import           Data.Maybe                 (fromJust)
import qualified Data.Vector.Unboxed        as Unboxed
import qualified Epidemic.BDSCOD            as SimBDSCOD
import           Epidemic.Types.Events      (EpidemicEvent (..), asNewickString,
                                             maybeEpidemicTree,
                                             maybeReconstructedTree)
import           Epidemic.Types.Parameter   (AbsoluteTime (..), Probability,
                                             Rate, TimeDelta (..), Timed (..))
import           Epidemic.Types.Population  (Identifier (..), People (..),
                                             Person (..), numPeople)
import           Epidemic.Utility           ()
import qualified Epidemic.Utility           as SimUtil
import           GHC.Generics               (Generic)
import           Numeric.GSL.Minimization   (MinimizeMethod (NMSimplex2),
                                             minimizeV)
import           Numeric.LinearAlgebra      (dot)
import           Numeric.LinearAlgebra.Data (Vector (..), fromList, linspace,
                                             toList)
import           Numeric.MCMC.Metropolis    (Chain (..), asGenIO, chain')
import           System.Environment         (getArgs)
import           System.Random.MWC          (initialize)

-- | Record of the scheduled observation times.
data ScheduledTimes =
  ScheduledTimes
    { stRhoTimes :: [AbsoluteTime]
    , stNuTimes  :: [AbsoluteTime]
    }
  deriving (Show)

-- | These objects define the specifics of the evaluation of LLHD profiles. If a
-- point estimate is given, then that is the central point of the profiles,
-- otherwise the parameters are estimated first. In every case the natural death
-- rate is set to the true values of the simulation. The boolean value is true
-- when the unscheduled observations are to be aggregated, otherwise they are
-- left as they are.
data InferenceConfiguration =
  InferenceConfiguration
    { reconstructedTreeOutputFiles :: (FilePath, FilePath)
    , observationsOutputCsv        :: FilePath
    , llhdOutputCsv                :: FilePath
    , pointEstimatesCsv            :: FilePath
    , maybePointEstimate           :: Maybe Parameters
    , icMaybeTimesForAgg           :: Maybe ([AbsoluteTime], [AbsoluteTime])
    , icMaybeMCMCConfig            :: Maybe MCMCConfiguration
    }
  deriving (Show, Generic)

instance Json.FromJSON InferenceConfiguration

-- | This object configures the whole evaluation of this program and is to be
-- read in from a suitable JSON file.
--
--     * A CSV to write the whole simulation to
--     * The parameters to use for the simulation
--     * The duration of the simulation
--     * The bounds on the size of an acceptable simulation
--     * One inference configuration for the regular data (with and without
--     estimated parameters) and another for the aggregated data with estimated
--     parameters.
--     * A Boolean to toggle printing progress on and off with the
--     @isVerbosePutStrLn@ function.
data Configuration =
  Configuration
    { simulatedEventsOutputCsv :: FilePath
    , simulationParameters :: Parameters
    , simulationDuration :: TimeDelta
    , simulationSizeBounds :: (Int, Int)
    , inferenceConfigurations :: ( InferenceConfiguration
                                 , InferenceConfiguration
                                 , InferenceConfiguration)
    , isVerbose :: Bool
    , configSimulationSeed :: MWCSeed
    }
  deriving (Show, Generic)

instance Json.FromJSON Configuration

-- | Monad for running a simulation, it uses the ExceptT for error messages and
-- the ReaderT for holding program configuration.
type Simulation x = ReaderT Configuration (ExceptT String IO) x

-- | A convenience function for printing output if in verbose mode when in the
-- 'Simulation' monad.
ifVerbosePutStrLn :: String -> Simulation ()
ifVerbosePutStrLn msg =
  do
    beLoud <- asks isVerbose
    when beLoud $ liftIO (putStrLn msg)

-- | This type is used to indicate if parameters are the true ones used in the
-- simulation or estimates parameters and in the case of estimated parameters,
-- what sort of data was used to inform the estimate.
data AnnotatedParameter
  = TrueParameters Parameters
  | EstimatedParametersRegularData Parameters
  | EstimatedParametersAggregatedData Parameters
  deriving (Show, Eq)

-- | A BDSCOD simulation configuration based on the parameters in the
-- environment.
bdscodConfiguration = do
  simParams@(Parameters (pLambda, pMu, pPsi, Timed pRhos, pOmega, Timed pNus)) <-
    asks simulationParameters
  if pLambda > 0 && pMu > 0 && pPsi > 0 && pOmega > 0 && null pRhos && null pNus
    then do
      (TimeDelta simDur) <- asks simulationDuration
      let bdscodConfig =
            SimBDSCOD.configuration (AbsoluteTime simDur) (unpackParameters simParams)
      case bdscodConfig of
        Nothing       -> throwError "Could not construct BDSCOD configuration"
        (Just config) -> return config
    else throwError "Simulation parameters not acceptable for this program..."

-- | Check if the LTT ever returns to 1 after being larger than 1.
pMultipleOrigins :: [EpidemicEvent] -> Bool
pMultipleOrigins = go 1
  where
    numLineage :: People -> NumLineages
    numLineage = fromIntegral . numPeople
    go :: NumLineages -> [EpidemicEvent] -> Bool
    go _ []                         = False
    go n (Removal _ _:ees)          = (n < 3) || go (n-1) ees
    go n (Sampling _ _:ees)         = (n < 3) || go (n-1) ees
    go n (Occurrence _ _:ees)       = (n < 3) || go (n-1) ees
    go n (Catastrophe _ people:ees) = if m > 1 then go m ees else True where m = n - numLineage people
    go n (Disaster _ people:ees)    = if m > 1 then go m ees else True where m = n - numLineage people
    go n (Infection {}:ees)      = go (n+1) ees

-- | Simulate the actual epidemic making sure that the results are acceptable
-- before returning, otherwise try using a different seed.
simulateEpidemic seedInt bdscodConfig = do
  ifVerbosePutStrLn "Running simulateEpidemic..."
  genIO <- liftIO $ prngGen seedInt
  simEvents <-
    liftIO $ SimUtil.simulation' bdscodConfig SimBDSCOD.allEvents genIO
  (sizeLowerBound, sizeUpperBound) <- asks simulationSizeBounds
  if length simEvents > sizeLowerBound && length simEvents < sizeUpperBound && not (pMultipleOrigins simEvents)
    then do
      ifVerbosePutStrLn "simulated an acceptable epidemic..."
      simEventsCsv <- asks simulatedEventsOutputCsv
      ifVerbosePutStrLn $
        "\twriting events to: " ++
        simEventsCsv ++ "\n\tthere are " ++ show (length simEvents) ++ " events."
      liftIO $ L.writeFile simEventsCsv (Csv.encode simEvents)
      return simEvents
    else do
      ifVerbosePutStrLn $ "\tthere where multiple origins: " ++ show (pMultipleOrigins simEvents)
      ifVerbosePutStrLn $ "\trepeating the simulation with seed: " ++ show seedInt
      simulateEpidemic (succ seedInt) bdscodConfig

-- | Take a simulated epidemic and generate the observations, first with full
-- resolution of the event times and the true epidemic parameters, second with
-- the event times and the estimated parameters and third with the sampling
-- times aggregated as defined in the inference configuration. This is done all
-- at the same time because the raw epidemic events are needed to genereate the
-- observations and the aggregated observations.
observeEpidemicThrice ::
     [EpidemicEvent] -- ^ the raw simulated epidemic events
  -> Simulation ( (InferenceConfiguration, [Observation])
                , (InferenceConfiguration, [Observation])
                , (InferenceConfiguration, AggregatedObservations))
observeEpidemicThrice simEvents = do
  (regInfConfig, regInfConfig', aggInfConfig) <- asks inferenceConfigurations
  let maybeRegObs = eventsAsObservations <$> SimBDSCOD.observedEvents simEvents
      maybeAggTimes =
        icMaybeTimesForAgg aggInfConfig >>= uncurry maybeAggregationTimes
      maybeAggObs =
        liftM2 aggregateUnscheduledObservations maybeAggTimes maybeRegObs
      (reconNewickTxt, reconNewickCsv) =
        reconstructedTreeOutputFiles regInfConfig
      maybeNewickData =
        asNewickString (AbsoluteTime 0, Person (Identifier 1)) =<<
        maybeReconstructedTree =<< maybeEpidemicTree simEvents
  case maybeNewickData of
    Just (newickBuilder, newickMetaData) -> do
      ifVerbosePutStrLn $
        "Writing reconstructed tree Newick to " ++ reconNewickTxt
      liftIO $
        L.writeFile reconNewickTxt (BBuilder.toLazyByteString newickBuilder)
      ifVerbosePutStrLn $
        "Writing reconstructed tree tip labels to " ++ reconNewickCsv
      liftIO $ L.writeFile reconNewickCsv (Csv.encode newickMetaData)
    Nothing -> throwError "Could not reconstruct tree..."
  case (maybeRegObs, maybeAggObs) of
    (Just regObs, Just aggObs@(AggregatedObservations _ unboxedAggObs)) -> do
      ifVerbosePutStrLn "Writing the regular observations to CSV"
      liftIO $
        L.writeFile (observationsOutputCsv regInfConfig) (Csv.encode regObs)
      ifVerbosePutStrLn "Writing the aggregated observations to CSV"
      liftIO $
        L.writeFile
          (observationsOutputCsv aggInfConfig)
          (Csv.encode unboxedAggObs)
      return
        ( (regInfConfig, regObs)
        , (regInfConfig', regObs)
        , (aggInfConfig, aggObs))
    (Just _, Nothing) ->
      throwError "Could not evaluate aggregated observations..."
    (Nothing, Just _) -> throwError "Could not evaluate regular observations..."
    (Nothing, Nothing) ->
      throwError "Could not evaluate either set of observations..."

-- | Evaluate the NB posterior approximation of the prevalence for a single
-- point in parameter space.
recordPresentPrevalenceEstimate ::
     InferenceConfiguration
  -> [Observation]
  -> AnnotatedParameter
  -> Simulation ()
recordPresentPrevalenceEstimate InferenceConfiguration {..} obs centerParam =
  let (parametersUsed,singleParams) = case centerParam of
                                        (TrueParameters x) -> ("true_parameters_regular_data" :: L.ByteString,x)
                                        (EstimatedParametersRegularData x) -> ("estimated_parameters_regular_data",x)
                                        (EstimatedParametersAggregatedData x) -> ("estimated_parameters_aggregated_data",x)
      nBVal = pure (parametersUsed, snd $ llhdAndNB obs singleParams initLlhdState)
   in do ifVerbosePutStrLn $ "\twriting NB values to: " ++ pointEstimatesCsv
         liftIO $ L.writeFile pointEstimatesCsv (Csv.encode nBVal)

-- | Run the estimation of the prevalence (at present) on the given set of regular
-- (i.e., disaggregated) observations at the parameters used to simulate the
-- data set and write the result to file.
evaluateLLHD :: InferenceConfiguration
             -> [Observation] -- ^ the regular observations
             -> Simulation ()
evaluateLLHD infConfig obs = do
  ifVerbosePutStrLn "Running evaluateLLHD..."
  simParams <- asks simulationParameters
  ifVerbosePutStrLn "\tUsing non-aggregated data with the true parameters..."
  ifVerbosePutStrLn $ show simParams
  recordPresentPrevalenceEstimate infConfig obs (TrueParameters simParams)

-- | Using regular (i.e., disaggregated) observations, estimate the parameters
-- used in the simulation and run the estimation of the prevalence (at present)
-- using this estimate. Finally generate MCMC samples of the posterior
-- distribution to understand the likelihood surface.
--
-- __NOTE__ This uses the actual simulation parameters as a way to get the
-- scheduled observation times, /they are not used in the inference/. The
-- optimisation starts at a fixed initial condition.
--
estimateLLHD :: InferenceConfiguration
             -> [Observation] -- ^ the regular observations
             -> Simulation ()
estimateLLHD infConfig obs = do
  ifVerbosePutStrLn "Running estimateLLHD..."
  Parameters (_, deathRate, _, _, _, _) <- asks simulationParameters
  let mleParams@(Parameters (estBirthRate, _, _, _, _, _)) = estimateRegularParameters deathRate obs
      annotatedMLE = EstimatedParametersRegularData mleParams
  ifVerbosePutStrLn "\tUsing non-aggregated data the computed MLE is..."
  ifVerbosePutStrLn $ "\tThe estimated birth rates is" ++ show estBirthRate
  recordPresentPrevalenceEstimate infConfig obs annotatedMLE
  runUnscheduledObservationMCMC infConfig deathRate obs annotatedMLE

-- | Using __aggregated__ observations, estimate the parameters of the process
-- and run the estimation of the prevalence (at present) using this estimate.
-- Finally generate MCMC samples of the posterior distribution to understand the
-- likelihood surface.
estimateLLHDAggregated :: InferenceConfiguration
                       -> AggregatedObservations -- ^ the aggregated observations
                       -> Simulation ()
estimateLLHDAggregated infConfig (AggregatedObservations _ obs) = do
  ifVerbosePutStrLn "Running estimateLLHDAggregated..."
  Parameters (_, deathRate, _, _, _, _) <- asks simulationParameters
  let (Just (rhoTimes,nuTimes)) = icMaybeTimesForAgg infConfig
      schedTimes = ScheduledTimes {stRhoTimes = rhoTimes, stNuTimes = nuTimes}
      mleParams = estimateAggregatedParameters deathRate schedTimes obs
      annotatedMLE = EstimatedParametersAggregatedData mleParams
  ifVerbosePutStrLn "\tUsing aggregated data the computed MLE is..."
  ifVerbosePutStrLn $ "\t" ++ show mleParams
  recordPresentPrevalenceEstimate infConfig obs annotatedMLE
  runScheduledObservationMCMC infConfig deathRate schedTimes obs annotatedMLE

-- | Use GSL to estimate the MLE based on the observations given. This uses a
-- simplex method because it seems to be faster and more accurate than the
-- simulated annealing and does not require a derivative.
--
-- __NOTE__ This will fit the model assuming that we have __aggregated__
-- observations, i.e., it will crash if there are unscheduled samples in the
-- data.
--
-- __NOTE__ we fix the death rate to the true value, which is given as an
-- argument, because this is assumed to be known a priori also, since this is
-- being applied to aggregated data we can set the value of psi and omega to
-- zero.
--
-- __NOTE__ The @energyFunc@ uses a regularisation cost to prevent the
-- parameters from wandering off which can occur with smaller data sets.
--
estimateAggregatedParameters :: Rate -- ^ the known death rate
                             -> ScheduledTimes -- ^ the time of scheduled observations
                             -> [Observation] -- ^ the aggregated observations
                             -> Parameters
estimateAggregatedParameters deathRate ScheduledTimes{..} obs =
  let maxIters = 20
      desiredPrec = 1e-1
      initBox = fromList ([0.1, 0.1, 0.1] :: [Rate]) -- lambda, rho and nu (because mu, psi, omega are fixed)
      randInit = fromList $ exp <$> [-1,-2,-2]

      vecAsParams :: Vector Rate -> Maybe Parameters
      vecAsParams x =
        let [lnR1, lnP1, lnP2] = toList x
            p1 = invLogit lnP1
            p2 = invLogit lnP2
            r = exp lnR1
            rts = [(t, p1) | t <- stRhoTimes]
            nts = [(t, p2) | t <- stNuTimes]
            packedParams = packParameters (r, deathRate, 0, rts, 0, nts)
          in if r < 5
             then Just packedParams
             else Nothing

      negLlhd :: Parameters -> LogLikelihood
      negLlhd ps = negate . fst $ llhdAndNB obs ps initLlhdState

      energyFunc :: Vector Rate -> Double
      energyFunc x = maybe 1e6 negLlhd $ vecAsParams x

      (est, _) = minimizeV NMSimplex2 desiredPrec maxIters initBox energyFunc randInit
   in fromJust $ vecAsParams est

-- | This is the main entry point to the actual simulation study, it is its own
-- function to avoid stressing about configuration IO. Since this is within the
-- simulation monad it has access to all the configuration data and can perform
-- IO.
simulationStudy :: Simulation ()
simulationStudy = do
  bdscodConfig <- bdscodConfiguration
  prngSeed <- asks configSimulationSeed
  epiSim <- simulateEpidemic prngSeed bdscodConfig
  (regObs, regObs', aggObs) <- observeEpidemicThrice epiSim
  uncurry evaluateLLHD regObs
  uncurry estimateLLHD regObs'
  uncurry estimateLLHDAggregated aggObs

-- | We use an external file to configure the simulation so that it is easier to
-- try different parameter values and random seeds. The @simulationStudy@
-- function is real driver of this and runs in the @Simulation@ monad.
main :: IO ()
main = do
  configFilePath <- head <$> getArgs
  maybeConfig <- Json.decode <$> L.readFile configFilePath
  case maybeConfig of
    Nothing ->
      putStrLn $ "Could not get configuration from file: " ++ configFilePath
    Just config -> do
      result <- runExceptT (runReaderT simulationStudy config)
      case result of
        Left errMsg -> putStrLn errMsg
        Right _     -> return ()

-- | Use GSL to estimate the MLE based on the observations given. This uses a
-- simplex method because it seems to be faster and more accurate than the
-- simulated annealing.
--
-- __NOTE__ This will fit the model assuming that we have __disaggregated__
-- observations.
--
-- __NOTE__ we fix the death rate to the true value because
-- this is assumed to be known a priori.
--
estimateRegularParameters ::
     Rate -> [Observation] -> Parameters
estimateRegularParameters deathRate obs =
  let maxIters = 50
      desiredPrec = 1e-2
      -- the only parameters are lambda, psi and omega (because mu, rho and nu
      -- are fixed)
      initBox = fromList [0.1, 0.1, 0.1]

      randInit :: Vector Double
      randInit = fromList [-1.5,-2.5,-2.5]

      vecAsParams :: Vector Double -> Maybe Parameters
      vecAsParams x =
        let [lnR1, lnR2, lnR3] = toList x
            r1 = exp lnR1
            r2 = exp lnR2
            r3 = exp lnR3
            packedParams = packParameters (r1, deathRate, r2, [], r3, [])
            rNaught = r1 / (deathRate + r2 + r3)
        in if r1 < 5 && 1 < rNaught && rNaught < 6
              -- constrain the R0 to be less than 5 and a sensible birth rate.
           then Just packedParams
           else Nothing

      negLlhd :: Parameters -> LogLikelihood
      negLlhd ps = negate . fst $ llhdAndNB obs ps initLlhdState

      energyFunc :: Vector Double -> Double
      energyFunc x = maybe 1e6 negLlhd $ vecAsParams x

      est :: Vector Double
      (est, path) =
        minimizeV NMSimplex2 desiredPrec maxIters initBox energyFunc randInit

   in fromJust (vecAsParams est)

-- | Action to run an MCMC analysis based on regular unscheduled observations.
runUnscheduledObservationMCMC :: InferenceConfiguration
                              -> Rate
                              -> [Observation]
                              -> AnnotatedParameter
                              -> Simulation ()
runUnscheduledObservationMCMC InferenceConfiguration {..} deathRate obs (EstimatedParametersRegularData (Parameters (mleR1, _, mleR2, _, mleR3, _))) =
  case icMaybeMCMCConfig of
    (Just mcmcConfig) ->
      let numIters = mcmcNumIters mcmcConfig
          stepSd = mcmcStepSD mcmcConfig
          variableNames = ["lambda", "psi", "omega", "nbSize", "nbProb"]
          x0 = [mleR1, mleR2, mleR3]
          listAsParams [r1, r2, r3] =
            packParameters (r1, deathRate, r2, [], r3, [])
          logPost x = fst $ llhdAndNB obs (listAsParams x) initLlhdState
          prngSeed = mcmcSeed mcmcConfig
          maybeGenQuantityFunc = Just (\x -> snd $ llhdAndNB obs (listAsParams x) initLlhdState)
       in do ifVerbosePutStrLn "Running runUnscheduledObservationMCMC..."
             genIO <- liftIO $ prngGen prngSeed
             chainVals <-
               liftIO $ (asGenIO $ chain' numIters stepSd x0 logPost maybeGenQuantityFunc) genIO
             liftIO $
               L.writeFile
                 (mcmcOutputCSV mcmcConfig)
                 (chainAsByteString variableNames chainVals)
    Nothing -> ifVerbosePutStrLn "No MCMC configuration found!"

runScheduledObservationMCMC ::
     InferenceConfiguration
  -> Rate -- ^ the /known/ death rate
  -> ScheduledTimes
  -> [Observation]
  -> AnnotatedParameter
  -> Simulation ()
runScheduledObservationMCMC InferenceConfiguration {..} deathRate ScheduledTimes {..} obs (EstimatedParametersAggregatedData (Parameters (mleR1, _, _, _, _, _))) =
  case icMaybeMCMCConfig of
    (Just mcmcConfig) ->
      let numIters = mcmcNumIters mcmcConfig
          stepSd = mcmcStepSD mcmcConfig
          variableNames = ["lambda", "rho", "nu", "nbSize", "nbProb"]
          x0 = [mleR1, 0.5, 0.5]
          listAsParams [r1, p1, p2] =
            packParameters
              ( r1
              , deathRate
              , 0
              , [(rt, p1) | rt <- stRhoTimes]
              , 0
              , [(nt, p2) | nt <- stNuTimes])
          logPost :: [Double] -> Double
          logPost x@[r1, p1, p2] =
            if r1 > 0 && 0 < p1 && p1 < 1 && 0 < p2 && p2 < 1
              then fst $ llhdAndNB obs (listAsParams x) initLlhdState
              else (-1 / 0)
          prngSeed = mcmcSeed mcmcConfig
          maybeGenQuantityFunc =
            Just (\x -> snd $ llhdAndNB obs (listAsParams x) initLlhdState)
       in do ifVerbosePutStrLn "Running runScheduledObservationMCMC..."
             genIO <- liftIO $ prngGen prngSeed
             chainVals <-
               liftIO $
               (asGenIO $ chain' numIters stepSd x0 logPost maybeGenQuantityFunc)
                 genIO
             liftIO $
               L.writeFile
                 (mcmcOutputCSV mcmcConfig)
                 (chainAsByteString variableNames chainVals)
             return ()
    Nothing -> ifVerbosePutStrLn "No MCMC configuration found!"

-- | A bytestring representation of the MCMC samples.
chainAsByteString :: [String] -- ^ the names of the elements of the chain
                  -> [Chain [Double] NegativeBinomial] -- ^ the samples in the chain
                  -> L.ByteString
chainAsByteString varNames chainVals =
  let header = pack $ intercalate "," ("llhd" : varNames)
      nbParams cv =
        case chainTunables cv of
          Just (NegBinomSizeProb r p) -> [r, p]
          Just Zero                   -> [0, 0]
          Nothing                     -> [-1, -1]
      records =
        Csv.encode
          [chainScore cv : (chainPosition cv ++ nbParams cv) | cv <- chainVals]
      linebreak = singleton '\n'
   in mconcat [header, linebreak, records]

-- | A generator for random numbers from a seed.
prngGen seed = initialize (Unboxed.fromList [seed])

max_stars_repo_name: aezarebski/timtam

max_issues_repo_name: aezarebski/timtam

max_forks_repo_name: aezarebski/timtam


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE OverloadedStrings #-}
-- | Main entry point for executable, mainly for tests
module Main (main) where

import AI.Network.RNN
import AI.GeneticAlgorithm.Simple
import Control.Monad
import Control.Monad.Random hiding (fromList)
import qualified Data.Text as T
import qualified Data.Text.IO as T
import Data.IORef
import System.Directory
--import Debug.Trace
import System.Environment

import Data.Binary

import qualified Data.ByteString as BS
import qualified Data.ByteString.Lazy as BSL

import Numeric.LinearAlgebra.HMatrix

--main::IO()
--main = do
--    let m :: [Expr Double] = map (\n->Var ("m" ++ show n)) [1..400]
--        v :: [Expr Double] = map (\n->Var ("v" ++ show n)) [1..4]
--        mycost v m = sum $ listMProd m v
--        mygrad = grad (mycost (map auto v)) m
--        --mygrads = gradientDescent (mycost (map auto v))  m
--    print $ mygrad

main :: IO()
main = do
    --txt <- liftM (T.take 100) $ T.readFile "data/tinyshakespeare.txt"
    let txt= "hello world!"
        trainData = train txt
        fn = "output/lstm.one.helloworld"
        generateLength = min 50 (T.length txt)
    print "Text:"
    print txt
    args <- getArgs
    when (length args<2) $ error "rnn gradient <maxgen> | rmsprop <maxgen> | generic <maxgen>"
    let fn2=fn++"."++head args
    ex <- return False -- doesFileExist fn2
    rnn<-case args of
        ("gradient":mg:_) -> do
            let maxGen = read mg
            rd <- readEx fn2 ex trainData
            gradient fn2 rd generateLength maxGen
        ("rmsprop":mg:_) -> do
            let maxGen = read mg
            rd <- readEx fn2 ex trainData
            rmsprop fn2 rd generateLength maxGen
        ("genetic":mg:_) -> do
            let maxGen = read mg
            b <- if ex
               then do
                 rnn1::LSTMNetwork <- fromVector (tdRecSize trainData) <$> decodeFile fn2
                 return $ buildExisting rnn1 trainData
               else return $ buildTD trainData
--            let b = buildTDs trainData
            genetic fn2 b generateLength maxGen
        _ -> error "rnn gradient <maxgen> | rmsprop <maxgen> | generic <maxgen>"
    encodeFile fn2 $ toVector rnn
    print rnn
    where
      readEx fn ex trainData =
          if ex
                then do
                 r::LSTMNetwork <- fromVector (tdRecSize trainData) <$> decodeFile fn
                 return $ RNNData r trainData (fromIntegral $ length $ tdInputs trainData)
                else
                  evalRandIO $ buildTD trainData
      (train, gener) = (textToTrainDataS,generateS)
      genStep maxGen = max 1 $ maxGen `div` 10
      gradient fn (RNNData r td _) generateLength maxGen = learnGradientDescent r td $ test fn generateLength maxGen
      rmsprop fn (RNNData r td _) generateLength maxGen = learnRMSProp r td $ test fn generateLength maxGen
      genetic fn r generateLength maxGen = do
        fitnessList <- newIORef []
        (RNNData rnn _ _) <- runGAIO 64 0.2 r $ stopf2 fn generateLength maxGen fitnessList
        return rnn
      test fn generateLength maxGen rnn td gen= do
        when (mod gen (genStep maxGen) == 0) $ do
            let c= cost rnn td
            print (show gen ++ ":" ++ show c)
            t <- evalRandIO $ gener (tdData td) generateLength (size $ head $ tdInputs td) rnn
            print t
            let fn2 = fn ++ "." ++ (show gen)
            let bs= BSL.toStrict $ encode $ toVector rnn
            BS.writeFile fn2 bs
        return $ gen < maxGen
      stopf2 fn generateLength maxGen fs rd@((RNNData rnn td _),_) gen = do
        _ <- test fn generateLength maxGen rnn td gen
        stopf fs maxGen rd gen

-- | Build LSTM network
buildTD ::(Monad m,RandomGen g) =>  TrainData a -> RandT g m (RNNData LSTMNetwork a)
buildTD td = buildNetworkData td (tdRecSize td) lstmFullSize

buildTDIO ::(Monad m,RandomGen g) =>  TrainData a -> RandT g m (RNNData LSTMIO a)
buildTDIO td =
    let layers = 2
        sz   = 40
     --   tds    = TrainData (tdInputs td) (tdOutputs td) (tdRecSize td,sz,layers,tdRecSize td) (tdData td)
     in buildNetworkData td (tdRecSize td,sz,layers,tdRecSize td) lstmioFullSize

-- | Build LSTM network
buildTDs ::(Monad m,RandomGen g) =>  TrainData a -> RandT g m (RNNData [LSTMNetwork] a)
buildTDs td =
    let layers = 2
     --   tds = TrainData (tdInputs td) (tdOutputs td) (replicate layers $ tdRecSize td) (tdData td)
    in  buildNetworkData td (replicate layers $ tdRecSize td)(sum . map lstmFullSize)

-- | Build from existing data for genetic algorithm, to restart from where we were
buildExisting :: (Monad m,RandomGen g) =>  LSTMNetwork -> TrainData a -> RandT g m (RNNData LSTMNetwork a)
buildExisting rnn1 td@(TrainData is _ _ _) = do
  -- g <- getSplit
  --let rnn2 = evalRand (mutateNetwork rnn1) g
  return $ RNNData rnn1 td (fromIntegral $ length is)

-- | Build from existing data for genetic algorithm, to restart from where we were
buildExistings :: (Monad m,RandomGen g) =>  [LSTMNetwork] -> TrainData a -> RandT g m (RNNData [LSTMNetwork] a)
buildExistings rnns td@(TrainData is _ _ _) = do
  -- g <- getSplit
  --let tds = TrainData (tdInputs td) (tdOutputs td) (replicate (length rnns) $ tdRecSize td) (tdData td)
  --let rnn2 = evalRand (mutateNetwork rnn1) g
  return $ RNNData rnns td (fromIntegral $ length is)

-- | Build from existing data for genetic algorithm, to restart from where we were
buildExistingIO :: (Monad m,RandomGen g) =>  LSTMIO -> TrainData a -> RandT g m (RNNData LSTMIO a)
buildExistingIO rnn1 td@(TrainData is _ _ _) = do
  -- g <- getSplit
  --let tds = TrainData (tdInputs td) (tdOutputs td) (rnnsize rnn1) (tdData td)
  --let rnn2 = evalRand (mutateNetwork rnn1) g
  return $ RNNData rnn1 td (fromIntegral $ length is)

-- DRAGONS
--
--main :: IO()
--main = do
--   -- txt <- T.readFile "data/tinyshakespeare.txt"
--    let txt= "hello world!"
--        (is,os,m) = textToTrainData txt
--        sz = DM.size m
--        fn = "lstm.learn1.helloworld"
--    print sz
--    print $ length is
--    print $ length os
--    ex <- doesFileExist fn
--    (RNNData r _ _ _) <- if ex
--        then do
--          r <- read <$> readFile fn
--          return $ RNNData r [] [] 0
--        else
--          evalRandIO $ build sz is os
--    let ls = toList $ toVector r
--    --print ls
--    --print $ toList $ chead ocs
--    --print ls
--    --print is
--    --print $ lstmList sz ls (toList $ head is)
--    --print $ mapAccumL (lstmList sz) ls (map toList is)
--    --headJet $ tailJet $ jet $ grads
----    let c= cost sz (fromIntegral $ length is) (map (toList) is) (map (toList) os) ls
----    print c
----    let gs= grad (cost sz (auto $ fromIntegral $ length is) (map (map auto . toList) is) (map (map auto . toList) os)) ls
----    --let gs=jacobian (lstmList (lstmFullSize sz) (map auto ls)) (toList $ head os)
----    -- print gs
----    let ls2 = zipWith (\l g->l-g*0.1) ls gs
----    let c2= cost sz (fromIntegral $ length is) (map (toList) is) (map (toList) os) ls2
----    print c2
----    where
----        gFun :: [AD s (Sparse Double)]
----                   -> AD s (Sparse Double)
----        gFun =  map sigmoid
--    let (l,lis,los)=((fromIntegral $ length is),(map toList is),(map toList os))
--    -- lsf <- learn sz ls l lis los 0
--    let myl= length ls
--    print myl
--    -- let adCost = cost sz (auto l) (map (map auto) lis) (map (map auto) los)
--    -- lsf <- learn2 sz ls (replicate myl 0) (replicate myl 0) (replicate myl 0) l lis los 0 (m,(min 50 (T.length txt)))
--    lsf <- learn sz ls l lis los 0 (m,(min 50 (T.length txt)))
--    --print lsf
--    let rnn::LSTMNetwork = fromVector sz (fromList lsf)
--    writeFile fn $ show rnn
--    --print rnn
--    -- print $ length is
--    let s = size $ head is
----    -- print s
----    let tl= T.length txt
----    -- print tl
----    -- print m
----    let c= cost sz l lis los ls
----    print c
----    let c2= cost sz l lis los lsf
----    print c2
--    t <- evalRandIO $ generate m (min 50 (T.length txt)) s rnn
--    print t
--    return ()
--    where
--          test sz ls l is os gen (m,gl) = when ((mod gen 500) == 0) $ do
--                let c= cost sz l is os ls
--                print (show gen ++ ":" ++ show c)
--                let rnn::LSTMNetwork = fromVector sz (fromList ls)
--                t <- evalRandIO $ generate m gl (length $ head is) rnn
--                print t
--          learn sz ls l is os gen testInfo = do
--            test sz ls l is os gen testInfo
--            if gen < 10000
--                then do
--                    let ls2 = learnGradientDescent (cost sz (auto l) (map (map auto) is) (map (map auto) os)) ls
--                    --let gs= grad (cost sz (auto l) (map (map auto) is) (map (map auto) os)) ls
--                    --    ls2 = zipWith (\o g->o-g*0.1) ls gs
--                    learn sz ls2 l is os (gen + 1) testInfo
--                else return ls
--          learn2 sz ls rgs rgs2 ugs l is os gen testInfo= do
--            test sz ls l is os gen testInfo
--            if gen < 10000
--                then do
--                    let gs= force $ grad (cost sz (auto l) (map (map auto) is) (map (map auto) os)) ls
--                        rgup = force $ zipWith (\rg g-> 0.95 * rg + 0.05 * g) rgs ls
--                        rg2up = force $ zipWith (\rg2 g-> 0.95 * rg2 + 0.05 * (g ** 2)) rgs2 ls
--                        ugup = force $ zipWith4 (\ud zg rg rg2 -> 0.9 * ud - 1e-4 * zg / sqrt(rg2 - rg ** 2 + 1e-4)) ugs gs rgup rg2up
--                        -- ls2 = zipWith (\o g->o-g*0.1) ls gs
--                        ls2 = force $ zipWith (+) ls ugup
--                    learn2 sz ls2 rgup rg2up ugup l is os (gen + 1) testInfo
--                else return ls
--          learn3 _ ls _ _ _ _ [] [] _ = return ls
--          learn3 sz ls rgs rgs2 ugs l is os gen = do
--            let batchSize = 100
--                (is1,is2) = splitAt batchSize is
--                (os1,os2) = splitAt batchSize os
--            when ((mod gen 100) == 0) $ do
--                let c= cost sz l is1 os1 ls
--                print (show gen ++ ":" ++ show c)
--            --if gen < 1
--            --    then do
--            let gs= force $ grad (cost sz (auto l) (map (map auto) is1) (map (map auto) os1)) ls
--                rgup = force $ zipWith (\rg g-> 0.95 * rg + 0.05 * g) rgs ls
--                rg2up = force $ zipWith (\rg2 g-> 0.95 * rg2 + 0.05 * (g ** 2)) rgs2 ls
--                ugup = force $ zipWith4 (\ud zg rg rg2 -> 0.9 * ud - 1e-4 * zg / sqrt(rg2 - rg ** 2 + 1e-4)) ugs gs rgup rg2up
--                -- ls2 = zipWith (\o g->o-g*0.1) ls gs
--                ls2 = force $ zipWith (+) ls ugup
--            learn3 sz ls2 rgup rg2up ugup l is2 os2 (gen + 1)
--            --    else return ls

--main1 :: IO ()
--main1 = do
----    n<-createRandomNetwork (RNNDimensions 1 2 3 True)
----    let (n2,out)=evalStep n [1]
----        (n3,out1)=evalStep n2 [3]
----        (n4,out2)=evalSteps n [[1],[3]]
----    print $ n3==n4
----    print $ out1 == (last out2)
----  txt <- liftM (T.take 100) $ T.readFile "data/tinyshakespeare.txt"
--  let -- dim =  (RNNDimensions 4 6 4 True)
--      --is  = [[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,1,0],[0,0,0,1]]
--      --os  = [[0,1,0,0],[0,0,1,0],[0,0,1,0],[0,0,0,1],[0,0,0,0]]
--      txt= "hello world!"
--      (is,os,m) = textToTrainData txt
--      sz = DM.size m
--      fn = "out.lstm.h"
--     --  dim = RNNDimensions sz (sz) sz True
--  print sz
--  fitnessList <- newIORef []
----  rnn1::LSTMNetwork <- read <$> readFile fn
--  ex <- doesFileExist fn
--  b <- if ex
--    then do
--        rnn1::LSTMNetwork <- read <$> readFile fn
--        return $ buildExisting rnn1 sz is os
--    else return $ build sz is os
--  r@(RNNData rnn _ _ _) <- runGAIO 64 0.1 b $ stopf fitnessList 100
----  -- (RNNData rnn _ _) <- evalRandIO $ buildNetworkData dim is os
----  print $ fitness r
--  -- print rnn
--  --let real = snd $ evalSteps rnn is
--  --print $ dataToText m real
--  t <- evalRandIO $ generate m (min 50 (T.length txt)) (size $ head is) rnn
--  print t
--  -- writeFile fn $ show rnn

-- build ::(Monad m,RandomGen g) =>  Int -> [Vector Double] -> [Vector Double] -> RandT g m (RNNData RNNetwork RNNDimensions)
-- build sz is os = buildNetworkData (RNNDimensions sz (sz) sz True) totalDataLength is os

--build ::(Monad m,RandomGen g) =>  TrainData a Int  -> RandT g m (RNNData LSTMNetwork Int)
--build sz td = buildNetworkData sz lstmFullSize is os

max_stars_repo_name: JPMoresmau/rnn

max_issues_repo_name: JPMoresmau/rnn

max_forks_repo_name: JPMoresmau/rnn


text:
import Data.Complex

mandelbrot a = iterate (\z -> z^2 + a) 0 !! 50

main = mapM_ putStrLn [[if magnitude (mandelbrot (x :+ y)) < 2 then '*' else ' '
                           | x <- [-2, -1.9685 .. 0.5]]
                       | y <- [1, 0.95 .. -1]]

max_stars_repo_name: djgoku/RosettaCodeData

max_issues_repo_name: djgoku/RosettaCodeData

max_forks_repo_name: djgoku/RosettaCodeData


text:
module NeuralNet.ExampleSpec (exampleSpec) where

import Test.Hspec
import Numeric.LinearAlgebra
import NeuralNet.Example
import Control.Exception (evaluate)

exampleSpec :: SpecWith ()
exampleSpec =
  describe "NeuralNet.Example" $ do
    describe "createExampleSet" $ do
      it "creates x correctly" $
        let
          examples = [([1, 2, 3], 4), ([5, 6, 7], 8)]
          set = createExampleSet examples
        in exampleSetX set `shouldBe` (fromLists [[1, 5], [2, 6], [3, 7]])

      it "creates y correctly" $
        let
          examples = [([1, 2, 3], 4), ([5, 6, 7], 8)]
          set = createExampleSet examples
        in exampleSetY set `shouldBe` (fromLists [[4, 8]])

      it "throws error for misaligned example sizes" $
        let examples = [([1..10], 4), ([1..2], 8)]
        in evaluate (createExampleSet examples) `shouldThrow` anyErrorCall

      it "throws error for empty example" $
        let examples = [([], 4)]
        in evaluate (createExampleSet examples) `shouldThrow` anyErrorCall

      it "throws error for empty examples" $
        evaluate (createExampleSet []) `shouldThrow` anyErrorCall

max_stars_repo_name: danielholmes/neural-net

max_issues_repo_name: danielholmes/neural-net

max_forks_repo_name: danielholmes/neural-net


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts      #-}
{-|
Module      : Grenade.Core.Pad
Description : Padding layer for 2D and 3D images
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Pad (
    Pad (..)
  ) where

import           Data.Maybe
import           Data.Proxy
import           Data.Serialize
import           Data.Singletons.TypeLits
import           GHC.TypeLits

import           Grenade.Core
import           Grenade.Layers.Internal.Pad

import           Numeric.LinearAlgebra (konst, subMatrix, diagBlock)
import           Numeric.LinearAlgebra.Static (extract, create)

-- | A padding layer for a neural network.
--
--   Pads on the X and Y dimension of an image.
data Pad  :: Nat
          -> Nat
          -> Nat
          -> Nat -> * where
  Pad  :: Pad padLeft padTop padRight padBottom

instance Show (Pad padLeft padTop padRight padBottom) where
  show Pad = "Pad"

instance UpdateLayer (Pad l t r b) where
  type Gradient (Pad l t r b) = ()
  runUpdate _ x _ = x
  createRandom = return Pad

instance Serialize (Pad l t r b) where
  put _ = return ()
  get = return Pad

-- | A two dimentional image can be padded.
instance ( KnownNat padLeft
         , KnownNat padTop
         , KnownNat padRight
         , KnownNat padBottom
         , KnownNat inputRows
         , KnownNat inputColumns
         , KnownNat outputRows
         , KnownNat outputColumns
         , (inputRows + padTop + padBottom) ~ outputRows
         , (inputColumns + padLeft + padRight) ~ outputColumns
         ) => Layer (Pad padLeft padTop padRight padBottom) ('D2 inputRows inputColumns) ('D2 outputRows outputColumns) where
  type Tape (Pad padLeft padTop padRight padBottom) ('D2 inputRows inputColumns) ('D2 outputRows outputColumns)  = ()
  runForwards Pad (S2D input) =
    let padl  = fromIntegral $ natVal (Proxy :: Proxy padLeft)
        padt  = fromIntegral $ natVal (Proxy :: Proxy padTop)
        padr  = fromIntegral $ natVal (Proxy :: Proxy padRight)
        padb  = fromIntegral $ natVal (Proxy :: Proxy padBottom)
        m     = extract input
        r     = diagBlock [konst 0 (padt,padl), m, konst 0 (padb,padr)]
    in  ((), S2D . fromJust . create $ r)
  runBackwards Pad _ (S2D dEdy) =
    let padl  = fromIntegral $ natVal (Proxy :: Proxy padLeft)
        padt  = fromIntegral $ natVal (Proxy :: Proxy padTop)
        nrows = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        ncols = fromIntegral $ natVal (Proxy :: Proxy inputColumns)
        m     = extract dEdy
        vs    = subMatrix (padt, padl) (nrows, ncols) m
    in  ((), S2D . fromJust . create $ vs)

-- | A two dimentional image can be padded.
instance ( KnownNat padLeft
         , KnownNat padTop
         , KnownNat padRight
         , KnownNat padBottom
         , KnownNat inputRows
         , KnownNat inputColumns
         , KnownNat outputRows
         , KnownNat outputColumns
         , KnownNat channels
         , KnownNat (inputRows * channels)
         , KnownNat (outputRows * channels)
         , (inputRows + padTop + padBottom) ~ outputRows
         , (inputColumns + padLeft + padRight) ~ outputColumns
         ) => Layer (Pad padLeft padTop padRight padBottom) ('D3 inputRows inputColumns channels) ('D3 outputRows outputColumns channels) where
  type Tape (Pad padLeft padTop padRight padBottom) ('D3 inputRows inputColumns channels) ('D3 outputRows outputColumns channels)  = ()
  runForwards Pad (S3D input) =
    let padl  = fromIntegral $ natVal (Proxy :: Proxy padLeft)
        padt  = fromIntegral $ natVal (Proxy :: Proxy padTop)
        padr  = fromIntegral $ natVal (Proxy :: Proxy padRight)
        padb  = fromIntegral $ natVal (Proxy :: Proxy padBottom)
        outr  = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        outc  = fromIntegral $ natVal (Proxy :: Proxy outputColumns)
        inr   = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        inc   = fromIntegral $ natVal (Proxy :: Proxy inputColumns)
        ch    = fromIntegral $ natVal (Proxy :: Proxy channels)
        m     = extract input
        padded = pad ch padl padt padr padb inr inc outr outc m
    in  ((), S3D . fromJust . create $ padded)

  runBackwards Pad () (S3D gradient) =
    let padl  = fromIntegral $ natVal (Proxy :: Proxy padLeft)
        padt  = fromIntegral $ natVal (Proxy :: Proxy padTop)
        padr  = fromIntegral $ natVal (Proxy :: Proxy padRight)
        padb  = fromIntegral $ natVal (Proxy :: Proxy padBottom)
        outr  = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        outc  = fromIntegral $ natVal (Proxy :: Proxy outputColumns)
        inr   = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        inc   = fromIntegral $ natVal (Proxy :: Proxy inputColumns)
        ch    = fromIntegral $ natVal (Proxy :: Proxy channels)
        m     = extract gradient
        cropped = crop ch padl padt padr padb inr inc outr outc m
    in  ((), S3D . fromJust . create $ cropped)


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeOperators    #-}
module Filter.Gaussian where

import           Control.Monad       (when)
import           Data.Array.Repa     as R
import           Data.Complex
import           Data.List           as L (length, reverse, (++))
import           Data.Vector.Generic as VG (convert)
import           DFT.Plan
import           Filter.Utils

data Gaussian2DParams = Gaussian2DParams
  { getGaussianFilterSigma :: !Double
  , getGaussianFilterRows  :: !Int
  , getGaussianFilterCols  :: !Int
  } deriving (Show)

{-# INLINE gaussian2D #-}
gaussian2D
  :: (Floating a)
  => a -> Int -> Int -> a
gaussian2D sd i j =
  1 / ((2 * pi) * sd * sd) * exp (-r / (2 * (sd ^ (2 :: Int))))
  where
    r = fromIntegral (i * i + j * j)

gaussian2DFilter ::
     DFTPlan -> Gaussian2DParams -> IO (DFTPlan, Array U DIM2 (Complex Double))
gaussian2DFilter plan (Gaussian2DParams sigma rows cols) = do
  let filterVec =
        VG.convert .
        toUnboxed .
        computeS .
        R.map (:+ 0) . makeFilter2D . R.fromFunction (Z :. cols :. rows) $ \(Z :. i :. j) ->
          gaussian2D sigma (i - (div cols 2)) (j - (div rows 2))
  lock <- getFFTWLock
  (plan1, filterVecF) <- dft1dGPlan lock plan [cols, rows] [0, 1] filterVec
  (newPlan, _) <- idft1dGPlan lock plan1 [cols, rows] [0, 1] filterVecF
  return (newPlan, fromUnboxed (Z :. cols :. rows) . VG.convert $ filterVecF)

{-# INLINE convolveGaussian2D #-}
convolveGaussian2D ::
     (Source r (Complex Double), Shape sh)
  => DFTPlan
  -> Array U DIM2 (Complex Double)
  -> Array r (sh :. Int :. Int) (Complex Double)
  -> IO (Array U (sh :. Int :. Int) (Complex Double))
convolveGaussian2D plan filterF input = do
  let len = rank . extent $ input
  when
    (len < 2)
    (error $ "convolveGaussian2D: rank = " L.++ show len L.++ " is less than 2.")
  inputF <-
    fmap (fromUnboxed (extent input) . VG.convert) .
    dftExecute
      plan
      (DFTPlanID
         DFT1DG
         (L.reverse . listOfShape . extent $ input)
         [len - 2, len - 1]) .
    VG.convert . toUnboxed . computeS . delay $
    input
  fmap (fromUnboxed (extent input) . VG.convert) .
    dftExecute
      plan
      (DFTPlanID
         IDFT1DG
         (L.reverse . listOfShape . extent $ input)
         [len - 2, len - 1]) .
    VG.convert . toUnboxed . computeS . R.traverse2 inputF filterF const $ \fInput fFilter idx@(sh :. i :. j) ->
    fInput idx * fFilter (Z :. i :. j)

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE ScopedTypeVariables #-}

module CSVMatrixParser
  (
    getData
  ) where

import Numeric.LinearAlgebra
import qualified Data.ByteString.Lazy as BL
import Data.Csv
import qualified Data.Vector as V

getCSV :: String -> IO (V.Vector [Double])
getCSV path = do
  csvData <- BL.readFile path
  case decode NoHeader csvData of
    Left err -> return V.empty
    Right v -> return v

splitIO :: [Double] -> ([Double], [Double])
splitIO list = splitAt 5 list

toVectorRow :: ([Double], [Double]) -> (Vector R, Vector R)
toVectorRow (xs,ys) = (vector xs, vector ys)

concatMats :: V.Vector (Vector R, Vector R) -> (Matrix R, Matrix R)
concatMats vec = let list = V.toList vec
                 in (fromRows (map fst list), fromRows (map snd list))

getData :: String -> IO (Matrix R, Matrix R)
getData path = getCSV path >>= return . concatMats . fmap (toVectorRow . splitIO)

concatTupleMats :: (Matrix R, Matrix R) -> (Matrix R, Matrix R) -> (Matrix R, Matrix R)
concatTupleMats (a,x) (b,y) = (a === b, x === y)