text:
{-# LANGUAGE AllowAmbiguousTypes        #-}
{-# LANGUAGE DeriveGeneric              #-}
{-# LANGUAGE FlexibleContexts           #-}
{-# LANGUAGE GADTs                      #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE InstanceSigs               #-}
{-# LANGUAGE LambdaCase                 #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE StandaloneDeriving         #-}
{-# LANGUAGE TypeApplications           #-}
{-# LANGUAGE TypeFamilies               #-}
{-# LANGUAGE TypeFamilyDependencies     #-}
{-# LANGUAGE TypeInType                 #-}
{-# LANGUAGE TypeOperators              #-}
{-# LANGUAGE UndecidableInstances       #-}
{-# LANGUAGE ViewPatterns               #-}

module Numeric.BLAS.HMatrix (
    HM(..)
  , HM'
  , flatten
  , unflatten
  ) where

import           Control.Applicative
import           Control.DeepSeq
import           Control.Monad
import           Control.Monad.Trans.State
import           Data.Finite.Internal
import           Data.Foldable
import           Data.Function
import           Data.Kind
import           Data.Maybe
import           Data.MonoTraversable
import           Data.Monoid                  (Endo(..))
import           Data.Singletons
import           Data.Singletons.Prelude
import           Data.Singletons.TypeLits
import           Data.Type.Product
import           GHC.Generics                 (Generic)
import           GHC.TypeLits
import           Numeric.BLAS hiding          (outer)
import           Numeric.LinearAlgebra.Static
import           Type.Family.List hiding      (Reverse)
import qualified Data.Vector                  as UV
import qualified Data.Vector.Generic.Sized    as VG
import qualified Data.Vector.Sized            as V
import qualified Data.Vector.Storable         as UVS
import qualified Data.Vector.Storable.Sized   as VS
import qualified Numeric.LinearAlgebra        as LA

type family HM' (s :: [Nat]) = (h :: Type) | h -> s where
    HM' '[]            = Double
    HM' '[n]           = R n
    HM' '[m,n]         = L m n
    HM' (n ': m ': ms) = V.Vector n (HM' (m ': ms))

newtype HM :: [Nat] -> Type where
    HM  :: { getHM :: HM' b }
        -> HM b
  deriving (Generic)

type instance Element (HM s) = Double

instance BLAS HM where

    transp (HM x) = HM (tr x)

    scal α (HM x)        = HM (konst α * x)
    axpy α (HM x) (HM y) = HM (konst α * x + y)
    dot    (HM x) (HM y) = x <.> y
    norm2  (HM x)        = norm_2 x
    asum   (HM x)        = norm_1 x

    iamax
        :: forall n. KnownNat n
        => HM '[n + 1]
        -> Finite (n + 1)
    iamax  (HM x)        = withKnownNat (SNat @n %:+ SNat @1) $
        Finite . fromIntegral
      . LA.maxIndex . extract
      $ abs x

    gemv α (HM a) (HM x) = \case
        Just (β, HM y) -> HM (konst α * (a #> x) + konst β * y)
        Nothing        -> HM (konst α * (a #> x))
    ger  α (HM x) (HM y) = \case
        Just (HM a) -> HM (konst α * outer x y + a)
        Nothing     -> HM (konst α * outer x y)
    gemm α (HM a) (HM b) = \case
        Just (β, HM c) -> HM (konst α * (a <> b) + konst β * c)
        Nothing        -> HM (konst α * (a <> b))

instance Tensor HM where
    type Scalar  HM = Double

    gen = \case
      SNil -> \f -> HM $ f Ø
      n@SNat `SCons` SNil -> \f -> HM . fromJust . create $
        LA.build (fromIntegral (fromSing n))
          (f . only . Finite . round)
      m@SNat `SCons` (n@SNat `SCons` SNil) -> \f -> HM . fromJust . create $
        LA.build (fromIntegral (fromSing m), fromIntegral (fromSing n))
          (\i j -> f (Finite (round i) :< Finite (round j) :< Ø))
      SNat `SCons` ns@(_ `SCons` (_ `SCons` _)) -> \f -> HM $
        V.generate_ $ \i -> getHM $ gen ns (\js -> f (i :< js))

    genA = \case
      SNil -> \f -> HM <$> f Ø
      n@SNat `SCons` SNil -> \f ->
        fmap (HM . vector) . traverse (f . only . Finite) $
          [0 .. fromSing n - 1]
      m@SNat `SCons` (n@SNat `SCons` SNil) -> \f ->
        fmap (HM . matrix) . traverse f $
          [ Finite i :< Finite j :< Ø | j <- [0 .. fromSing m - 1]
                                      , i <- [0 .. fromSing n - 1]
          ]
      SNat `SCons` ns@(_ `SCons` (_ `SCons` _)) -> \f -> fmap HM $
        sequenceA . V.generate_ $ \i -> getHM <$> genA ns (\js -> f (i :< js))

    tkonst s = HM . hkonst s

    tsum :: forall s. SingI s => HM s -> Double
    tsum = go sing . getHM
      where
        go :: Sing ns -> HM' ns -> Double
        go = \case
          SNil                                      -> id
          SNat `SCons` SNil                         -> LA.sumElements . extract
          SNat `SCons` (SNat `SCons` SNil)          -> LA.sumElements . extract
          SNat `SCons` ns@(_ `SCons` (_ `SCons` _)) -> sum . fmap (go ns)

    tmap :: forall s. SingI s => (Double -> Double) -> HM s -> HM s
    tmap = omap

    tzip
        :: forall s. SingI s
        => (Double -> Double -> Double)
        -> HM s -> HM s -> HM s
    tzip f (HM x0) (HM y0) = HM $ hzip f sing x0 y0

    tindex
        :: forall s. SingI s
        => Prod Finite s
        -> HM s
        -> Double
    tindex ix0 = go sing ix0 . getHM
      where
        go :: Sing ns -> Prod Finite ns -> HM' ns -> Double
        go = \case
          SNil -> \case
            Ø -> id
          SNat `SCons` SNil -> \case
            i :< Ø -> (`LA.atIndex` fromIntegral i) . extract
          SNat `SCons` (SNat `SCons` SNil) -> \case
            i :< j :< Ø -> (`LA.atIndex` (fromIntegral i, fromIntegral j)) . extract
          SNat `SCons` ns@(_ `SCons` (_ `SCons` _)) -> \case
            i :< js -> go ns js . (`V.index` i)

    tslice
        :: forall n m. ()
        => ProdMap Slice n m
        -> HM n
        -> HM m
    tslice sl0 (HM x0) = HM $ go sl0 x0
      where
        go  :: forall ns ms. ()
            => ProdMap Slice ns ms
            -> HM' ns
            -> HM' ms
        go = \case
          PMZ -> \x -> x
          PMS (Slice sL sC@SNat sR) PMZ -> \xs -> fromJust . create $
            let l = fromIntegral $ fromSing sL
                c = fromIntegral $ fromSing sC
            in  withKnownNat (sL %:+ sC %:+ sR) $
                  LA.subVector l c (extract xs)
          PMS (Slice sLy sCy@SNat sRy) (PMS (Slice sLx sCx@SNat sRx) PMZ) ->
              \xs -> fromJust . create $
            let lx = fromIntegral $ fromSing sLx
                ly = fromIntegral $ fromSing sLy
                cx = fromIntegral $ fromSing sCx
                cy = fromIntegral $ fromSing sCy
            in  withKnownNat (sLy %:+ sCy %:+ sRy) $
                withKnownNat (sLx %:+ sCx %:+ sRx) $
                  LA.subMatrix (ly, lx) (cy, cx) (extract xs)
          PMS (Slice sL sC@SNat _) pm@(PMS _ (PMS _ _)) -> \xs ->
            let l = fromIntegral $ fromSing sL
                c = fromIntegral $ fromSing sC
            in  fmap (go pm)
                  . fromJust . V.toSized
                  . UV.take c
                  . UV.drop l
                  . V.fromSized
                  $ xs

    tconv
        :: forall n m s. ()
        => Sing n
        -> ProdMap Conv m s
        -> HM (m >: n)      -- ^ mask
        -> HM s
        -> HM (s >: n)
    tconv sn@SNat pm0 (HM m0) (HM x0) = HM $ go pm0 m0 x0
      where
        go  :: forall ms ss. ()
            => ProdMap Conv ms ss
            -> HM' (ms >: n)
            -> HM' ss
            -> HM' (ss >: n)
        go = \case
          PMZ -> \m x -> konst x * m
          PMS (Conv sM sS str off) PMZ -> \m x ->
            undefined

    treshape
        :: (SingI s1, Product s1 ~ Product s2)
        => Sing s2
        -> HM s1
        -> HM s2
    treshape s = unflatten s . flatten

    tload
        :: Sing s
        -> V.Vector (Product s) Double
        -> HM s
    tload s = unflatten s . VG.convert

    textract
        :: SingI s
        => HM s
        -> V.Vector (Product s) Double
    textract = VG.convert . flatten

-- hconv
--     :: forall o ms ns. (KnownNat o)
--     => DoubleProd Sing ms ns
--     -> HM' (o ': ms)
--     -> HM' ns
--     -> HM' (o ': ns)
-- hconv = \case
--     DPZ -> \ms x -> konst x * ms
--     DPS SNat sn@SNat DPZ -> \m x -> fromJust . create $
--       let c = LA.conv2 (extract m) (LA.asRow (extract x))
--           o = fromInteger (natVal (Proxy @o))
--           left = o `div` 2
--       in  LA.subMatrix (0,left) (o, fromInteger (fromSing sn)) c
--     DPS smx@SNat snx@SNat (DPS smy@SNat sny@SNat DPZ) -> \m x ->
--       -- todo: vectorize with im2colV
--       flip fmap m $ \msk -> fromJust . create $
--         let c = LA.conv2 (extract msk) (extract x)
--             left = fromInteger (fromSing smx) `div` 2
--             top  = fromInteger (fromSing smy) `div` 2
--         in  LA.subMatrix (left, top) (fromInteger (fromSing snx), fromInteger (fromSing sny)) c
--     dp0@(DPS (SNat :: Sing m0) (SNat :: Sing n0) (DPS _ _ (DPS _ _ _) :: DoubleProd Sing ms0 ns0)) ->
--               \(ms :: V.Vector o (V.Vector m0 (HM' ms0))) (xs :: V.Vector n0 (HM' ns0)) ->
--       flip fmap ms $ \(msk :: V.Vector m0 (HM' ms0)) -> hconv1 dp0 msk xs

-- hconv1
--     :: forall m s. ()
--     => DoubleProd Sing m s
--     -> HM' m
--     -> HM' s
--     -> HM' s
-- hconv1 = \case
--     DPZ -> (*)
--     DPS sm@SNat sn@SNat DPZ -> \m x -> fromJust . create $
--       let c = LA.conv (extract m) (extract x)
--           left = fromInteger (fromSing sm) `div` 2
--       in  UVS.slice left (fromInteger (fromSing sn)) c
--     DPS smx@SNat snx@SNat (DPS smy@SNat sny@SNat DPZ) -> \m x -> fromJust . create $
--       let c = LA.conv2 (extract m) (extract x)
--           left = fromInteger (fromSing smx) `div` 2
--           top  = fromInteger (fromSing smy) `div` 2
--       in  LA.subMatrix (left, top) (fromInteger (fromSing snx), fromInteger (fromSing sny)) c
--     DPS (SNat :: Sing m0) (SNat :: Sing n0) dps@(DPS _ _ (DPS _ _ _) :: DoubleProd Sing ms0 ns0) ->
--               \(ms :: V.Vector m0 (HM' ms0)) (xs :: V.Vector n0 (HM' ns0)) ->
--       let s   :: Sing ns0
--           s   = prodSing $ secondDP dps
--           cl :: V.Vector n0 (V.Vector m0 (HM' ns0))
--           cl = im2colV (hkonst s 0) xs
--       in  fmap (hsum s . V.zipWith (hconv1 dps) ms) cl

flatten :: SingI s => HM s -> VS.Vector (Product s) Double
flatten = hflatten sing . getHM

unflatten :: Sing s -> VS.Vector (Product s) Double -> HM s
unflatten s = HM . hunflatten s

hflatten
    :: Sing s
    -> HM' s
    -> VS.Vector (Product s) Double
hflatten = \case
    SNil -> VS.singleton
    SNat `SCons` SNil -> fromJust . VS.toSized . extract
    sn@SNat `SCons` (sm@SNat `SCons` SNil) -> case sn %:* sm of
      SNat -> fromJust . VS.toSized . LA.flatten . extract
    SNat `SCons` ns@(_ `SCons` (_ `SCons` _)) ->
      VG.convert . V.concatMap (VG.convert . hflatten ns)

hunflatten
    :: Sing s
    -> VS.Vector (Product s) Double
    -> HM' s
hunflatten = \case
    SNil -> VS.head
    SNat `SCons` SNil -> fromJust . create . VS.fromSized
    SNat `SCons` (sm@SNat `SCons` SNil) -> fromJust . create . LA.reshape (fromInteger (fromSing sm)) . VS.fromSized
    sn@SNat `SCons` ns@(_ `SCons` (_ `SCons` _)) -> case sProduct ns of
      sp@SNat -> fromJust
            . V.fromList
            . evalState (replicateM (fromInteger (fromSing sn)) $
                           hunflatten ns . fromJust . VS.toSized <$> state (UVS.splitAt (fromInteger (fromSing sp)))
                        )
            . VS.fromSized

im2col
    :: forall m n o. (KnownNat m, KnownNat n, KnownNat o)
    => L n m
    -> L n o
im2col = undefined

im2colV
    :: forall m n a. (KnownNat m, KnownNat n)
    => a
    -> V.Vector n a
    -> V.Vector n (V.Vector m a)
im2colV pad (V.fromSized->v) = V.generate $ \i ->
      fromJust . V.toSized $ UV.slice i m padded
  where
    padded = UV.concat [UV.replicate left pad, v, UV.replicate right pad]
    m :: Int
    m  = fromIntegral $ natVal (Proxy @m)
    left  = m `div` 2
    right = m - left

hadd :: Sing s -> HM' s -> HM' s -> HM' s
hadd = \case
    SNil                                     -> (+)
    SNat `SCons` SNil                        -> (+)
    SNat `SCons` (SNat `SCons` SNil)         -> (+)
    SNat `SCons` s@(_ `SCons` (_ `SCons` _)) -> liftA2 (hadd s)

hsum :: Foldable f => Sing s -> f (HM' s) -> HM' s
hsum s = foldl' (hadd s) (hkonst s 0)

hkonst :: Sing s -> Double -> HM' s
hkonst = \case
    SNil                                     -> id
    SNat `SCons` SNil                        -> konst
    SNat `SCons` (SNat `SCons` SNil)         -> konst
    SNat `SCons` s@(_ `SCons` (_ `SCons` _)) -> pure . hkonst s

hzip
    :: (Double -> Double -> Double)
    -> Sing s
    -> HM' s
    -> HM' s
    -> HM' s
hzip f = go
  where
    go :: Sing t -> HM' t -> HM' t -> HM' t
    go = \case
      SNil                             -> f
      SNat `SCons` SNil                -> zipWithVector f
      SNat `SCons` (SNat `SCons` SNil) ->
        (\xs ys -> matrix (zipWith f xs ys))
           `on` (concat . LA.toLists . extract)
      SNat `SCons` ns@(_ `SCons` (_ `SCons` _)) ->
         V.zipWith (go ns)

-- firstDP
--     :: DoubleProd f as bs
--     -> Prod f as
-- firstDP = \case
--     DPZ        -> Ø
--     DPS x _ xs -> x :< firstDP xs

-- secondDP
--     :: DoubleProd f as bs
--     -> Prod f bs
-- secondDP = \case
--     DPZ        -> Ø
--     DPS _ x xs -> x :< secondDP xs

prodSing
    :: Prod Sing as
    -> Sing as
prodSing = \case
    Ø       -> SNil
    x :< xs -> x `SCons` prodSing xs

instance SingI s => MonoFunctor (HM s) where
    omap f = HM . go sing . getHM
      where
        go :: Sing ns -> HM' ns -> HM' ns
        go = \case
          SNil                                      -> f
          SNat `SCons` SNil                         -> dvmap f
          SNat `SCons` (SNat `SCons` SNil)          -> dmmap f
          SNat `SCons` ns@(_ `SCons` (_ `SCons` _)) -> fmap (go ns)

helems :: forall s. SingI s => HM s -> [Double]
helems = flip appEndo [] . go sing . getHM
  where
    go :: Sing ns -> HM' ns -> Endo [Double]
    go = \case
      SNil                                      -> \x -> Endo (x:)
      SNat `SCons` SNil                         -> Endo . (++) . LA.toList . extract
      SNat `SCons` (SNat `SCons` SNil)          -> foldMap (Endo . (++)) . LA.toLists . extract
      SNat `SCons` ns@(_ `SCons` (_ `SCons` _)) -> foldMap (go ns)

instance SingI s => MonoFoldable (HM s) where
    ofoldMap f     = foldMap f . helems
    ofoldr f z     = foldr f z . helems
    ofoldl' f z    = foldl' f z . helems
    otoList        = helems
    oall f         = all f . helems
    oany f         = any f . helems
    onull          = (== 0) . olength
    olength _      = fromIntegral (product (fromSing (sing @_ @s)))
    olength64      = fromIntegral . olength
    ocompareLength = ocompareLength . helems
    otraverse_ f   = traverse_ f . helems
    ofor_ x        = for_ (helems x)
    omapM_ f       = traverse_ f . helems
    oforM_ x       = for_ (helems x)
    ofoldlM f x    = foldlM f x . helems
    ofoldMap1Ex f  = ofoldMap1Ex f . helems
    ofoldr1Ex f    = ofoldr1Ex f . helems
    ofoldl1Ex' f   = ofoldl1Ex' f . helems
    headEx         = headEx . helems
    lastEx         = lastEx . helems
    maximumByEx f  = maximumByEx f . helems
    minimumByEx f  = minimumByEx f . helems

deriving instance NFData (HM' s)     => NFData (HM s)
deriving instance Show (HM' s)       => Show (HM s)
deriving instance Num (HM' s)        => Num (HM s)
deriving instance Fractional (HM' s) => Fractional (HM s)
deriving instance Floating (HM' s)   => Floating (HM s)

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE ScopedTypeVariables, BangPatterns #-}
module ECC.Tester (eccMain, eccPrinter, eccMerger) where

import ECC.Types
import System.Random.MWC
import Control.Monad
--import Control.Monad.Primitive (PrimMonad, PrimState)
import Control.Monad.Trans (MonadIO(liftIO))
import Control.Concurrent
import Data.Word
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed  as U
import System.Random.MWC.Distributions
import Data.Monoid
import Statistics.Resampling.Bootstrap
import System.IO
import Data.Char
import System.Environment
import ECC.Estimate
import GHC.Conc
import Control.Concurrent.ParallelIO
import Numeric
import Data.Time.Clock
import qualified Data.UUID as UUID
import System.Directory
import qualified Data.Csv as CSV
import qualified Data.ByteString.Lazy as LBS
import Statistics.Types (Estimate (..), ConfInt(..), CL, confidenceLevel)
import qualified Data.Map.Strict as Map
import Data.Map.Strict (Map)
import Data.List as List
import Data.Ratio
import Text.Read
import Data.Word (Word32)
import Data.Maybe (isJust)

data Options = Options
        { codenames :: [String]
        , ebN0s     :: [EbN0]
        , verbose   :: Int
        , enough    :: Enough
        , logDir    :: FilePath -- directory or file to put the log
        , cmp       :: Bool
        , replay    :: Maybe FilePath -- file to run replays
        , seed      :: Maybe Word32
        } deriving Show

data Enough = BitErrorCount Int
            | MessageCount Int
        deriving Show

data TestRun = TestRun
       { test_encode_time    :: !Double  -- in seconds
       , test_decode_time    :: !Double -- in seconds
       , test_ber            :: !BEs
       }
      deriving Show

instance Monoid TestRun where
    mempty = TestRun 0 0 mempty
    TestRun en1 de1 ber1 `mappend` TestRun en2 de2 ber2
        = TestRun (en1 + en2) (de1 + de2) (ber1 `mappend` ber2)

defaultLogDir :: String
defaultLogDir = "log/"

-- | Give a 'Code' (set of possible Error Correcting Codes) and a printer, run the tests.
eccMain :: Code -> (Options -> [String] -> IO (ECC IO -> EbN0 -> TestRun -> IO Bool)) -> IO ()
eccMain code k = do
        args <- getArgs
        if null args
         then error $ "usage: <name> [-v<n>] [-b<n>] [-m<n>] [-l<file or dir/>] [-c] <EbN0_1> <EbN0_2> ... <Code Name> <Code Name>"
                   ++ "\ncodes: " ++ show code
         else eccTester (parseOptions args $ defaultOptions defaultLogDir) code k

parseOptions :: [String] -> Options -> Options
parseOptions (('-':'v':ns):rest) o
        | all isDigit ns = (parseOptions rest o) { verbose = read ns }
parseOptions (('-':'b':ns):rest) o
        | all isDigit ns = (parseOptions rest o) { enough = BitErrorCount (read ns) }
parseOptions (('-':'m':ns):rest) o
        | all isDigit ns = (parseOptions rest o) { enough = MessageCount (read ns) }
parseOptions (('-':'l':lg):rest) o
        | otherwise = (parseOptions rest o) { logDir = lg }
parseOptions (('-':'r':r):rest) o
        | otherwise = (parseOptions rest o) { replay = Just r }
parseOptions (('-':'s':s):rest) o
        | otherwise = (parseOptions rest o) { seed = Just (read s) }
parseOptions ("-c":rest) o
        | otherwise       = (parseOptions rest o) { cmp = True }
parseOptions (arg:rest) o =
        case reads arg of
          [(ebN0::EbN0,"")] -> opts { ebN0s = ebN0 : ebN0s opts }
          _                 -> opts { codenames = arg : codenames opts }
  where
     opts = parseOptions rest o
parseOptions [] o = o

defaultOptions :: String -> Options
defaultOptions d = Options {
    codenames = [],
    ebN0s = [],
    verbose = 0,
    enough = BitErrorCount 1000,
    logDir = d,
    cmp = False,
    replay = Nothing,
    seed = Nothing
  }

-- | A basic printer for our tests. Currently, we report on powers of two,
-- and accept a value if there are at least 1000 bit errors (say).
-- We also output a log file into a diretory, using a UUID style file.
eccPrinter :: Options -> [String] -> IO (ECC f -> EbN0 -> TestRun -> IO Bool)
eccPrinter opts names = do

   let tab1 = maximum (map length names)

   let rjust n xs = take (n - length xs) (cycle " ") ++ xs

   -- This generator is uses for the bootstrapping, and generating the UUID
   gen :: GenIO <- createSystemRandom
   start <- getCurrentTime

   (w1,w2,w3,w4) <- uniform gen
   let uuid = UUID.fromWords w1 w2 w3 w4

   logFileName <- if last (logDir opts) == '/' then do
       createDirectoryIfMissing True $ logDir opts
       return $ logDir opts ++ UUID.toString uuid
     else
       return $ logDir opts

     -- Don't write to the default log directory if a seed was manually
     -- provided (to avoid skewing the statistics, if the same seed is
     -- manually used for a number of runs).
   unless (logDir opts == defaultLogDir && isJust (seed opts))
     $ writeHeader logFileName

   putStrLn $ "#" ++
              rjust 7    "Time" ++ " " ++
              rjust tab1 "ECC" ++ " " ++
              rjust 5    "EbN0" ++ " " ++
              rjust 8    "Packets" ++ " " ++
              rjust 10   "Encode/s" ++ " " ++
              rjust 10   "Decode/s" ++ " " ++
              rjust 8    "Errors" ++ " " ++
              rjust 8    "BER" ++ " " ++
              ""

   return $ \  ecc ebN0 (TestRun tEn tDe bes) -> do
           est <- if sizeBEs bes <= 2
                  then return Nothing
                  else estimate gen 0.95 (message_length ecc) bes
           let accept = case enough opts of
                          BitErrorCount n -> sumBEs bes >= n
                          MessageCount n  -> sizeBEs bes >= n

           now <- getCurrentTime
           let diff = diffUTCTime now start

           putStrLn $
                    rjust 8 (showFFloat (Just 2) (realToFrac diff) "") ++ " " ++
                    rjust tab1 (name ecc) ++ " " ++
                    rjust 5 (showFFloat (Just 2) ebN0 "") ++ " " ++
                    rjust 8 (show (sizeBEs bes)) ++ " " ++
                    rjust 10 (showFFloat (Just 4) (fromIntegral (sizeBEs bes)/tEn) "") ++ " " ++
                    rjust 10 (showFFloat (Just 4) (fromIntegral (sizeBEs bes)/tDe) "") ++ " " ++
                    rjust 8 (show (sumBEs bes)) ++ " " ++
                    (case est of
                      Just e -> " " ++ rjust 20 (showEstimate e)
                      Nothing -> " 0.00e-0") ++
                    if accept then "." else ","
           hFlush stdout

           -- Log file, with raw data

           unless (logDir opts == defaultLogDir && isJust (seed opts))
             $ LBS.appendFile logFileName $ CSV.encode [ Result
                 ( realToFrac diff :: Double )
                 ( name ecc        :: String )
                 ( ebN0            :: Double )
                 ( sizeBEs bes     :: Int )
                 ( tEn        )
                 ( tDe        )
                 ( sumBEs bes      :: Int )
                 ( estPoint <$> est )
                 ( (confIntLDX . estError) <$> est )
                 ( (confIntUDX . estError) <$> est )
                 ( (confidenceLevel . confIntCL . estError) <$> est )
                 ( message_length ecc )
                 ( codeword_length ecc )
                 ( bes )
                 ]
           return accept

-- eccTester runs testECC at different ebN0 levels.
eccTester :: Options -> Code -> (Options -> [String] -> IO (ECC IO -> EbN0 -> TestRun -> IO Bool)) -> IO ()
eccTester opts (Code _ initializeC finalizeC f) k = do
   print opts
   let debug n msg | n <= verbose opts  = putStrLn msg
                   | otherwise  = return ()

   -- This generator is uses for the generation of bits
   gen :: GenIO <- case seed opts of
                     Just s -> initialize (V.singleton s)
                     _      -> createSystemRandom

   let mg :: (Show a, Show b, Eq b, Monad m) => (a -> m b) -> (a -> m b) -> a -> m b
       mg f1 f2 a = do
           r1 <- f1 a
           r2 <- f2 a
           if r1 == r2
            then return r1
            else fail $ "merge of two codes failed at runtime" ++ show (a,r1,r2)

   let jn :: Monad m => ECC m -> ECC m -> ECC m
       jn ecc1 ecc2
            | message_length ecc1 /= message_length ecc2
             = error "trying to combine two codes with different message lengths"
            | codeword_length ecc1 /= codeword_length ecc2
            = error "trying to combine two codes with different codeword lengths"
            | otherwise
            = ECC { name = name ecc1 ++ "+" ++ name ecc2
                  , encode = mg (encode ecc1) (encode ecc2)
                  , decode = mg (decode ecc1) (decode ecc2)
                  , message_length = message_length ecc1
                  , codeword_length = codeword_length ecc2
                  }

   vars <- initializeC

   eccs <- liftM concat
            $ sequence
            $ map (f vars)
            $ map splitCodename
            $ codenames opts

   let eccs' = if cmp opts then [foldr1 jn eccs] else eccs

   k2 <- k opts (map name eccs)

   -- First, try replay (which have fixed EbN0s)
   replay_log <- case replay opts of
                   Nothing -> pure []
                   Just fileName -> do
                        str <- readFile fileName
                        return [ r :: Log
                               | ln <- lines str
                               , r <- case readMaybe ln of
                                        Nothing -> []
                                        Just x -> [x]
                               ]
   let replayMe [] = return ()
       replayMe (Message m:TxCodeword tx:RxCodeword ebN0 rx:rest) = do
           sequence_
               [ do bes <- rxECC (verbose opts) ecc m ebN0 rx
                    k2 ecc ebN0 bes
               | ecc <- eccs'
               ]
           replayMe rest
       replayMe (_:rest) = replayMe rest

   replayMe replay_log

   -- Then, loop over all given EbN0s.
   sequence_
          [ sequence_
                [ do testECC (verbose opts) gen ebN0 ecc k2
                | ebN0 <- ebN0s opts
                ]
          | ecc <- eccs'
          ]

   finalizeC vars

   stopGlobalPool -- TODO: use a local pool instead?

-- Running a *multi* run of an ECC, giving a compounted TestRun result to the
-- continuation, and stopping when the continuation returns True.
testECC :: Int -> GenIO -> EbN0 -> ECC IO -> (ECC IO -> EbN0 -> TestRun -> IO Bool) -> IO ()
testECC verb gen ebN0 ecc k = do
   let debug n msg | n <= verb  = putStrLn msg
                   | otherwise  = return ()

   cap <- getNumCapabilities

   let ourpar = sequence -- parallel
   let cap = 1           -- for now

   let loop (n:ns) !bes = do
        debug 1 $ "trying " ++ show n ++ " messages"
        let real_par = min (cap * 16) n
        bess <- ourpar   [ do foldM (\ !a !_ -> do bes0 <- runECC verb gen ebN0 ecc
                                                   return $ a <> bes0) mempty [1..(n `div` real_par)]
                         | _ <- [1..real_par]
                         ]
        let bes1 = mconcat (bes:bess)
        let errs = sumBEs $ test_ber bes1
        debug 1 $ "found " ++ show errs ++ " so far"
        okay <- k ecc ebN0 bes1
        if okay then return () else loop ns bes1

   -- We run with the cap twice, then double each time
   loop (cap : iterate (*2) cap) mempty

splitCodename :: String -> [String]
splitCodename = words . map (\ c -> if c == '/' then ' ' else c)

-- Running a *single* run of an ECC, getting a single bit error count
runECC :: Int -> GenIO -> EbN0 -> ECC IO -> IO TestRun
runECC verb gen ebN0 !ecc = do

  let debug n msg | n <= verb  = putStrLn msg
                  | otherwise  = return ()

  debug 3 $ "starting message"
  !mess0  <- liftM U.fromList $ sequence [ uniform gen | _ <- [1..message_length ecc]]
  debug 3 $ "generated message"
  debug 4 $ show (Message mess0)

  start_encoding <- getCurrentTime

  !code0  <- encode ecc mess0
  debug 3 $ "encoded message"
  debug 4 $ show (TxCodeword code0)

  end_encoding <- getCurrentTime

  !rx <- txRx_EbN0 ebN0 (rateOf ecc) gen code0
  debug 3 $ "tx/rx'd message"
  debug 4 $ show (RxCodeword ebN0 rx)

  rxECC verb ecc mess0 (realToFrac $ diffUTCTime end_encoding start_encoding) rx

rxECC :: Int -> ECC IO -> U.Vector Bool -> Double -> U.Vector Double -> IO TestRun
rxECC verb ecc mess0 encoding_time rx = do

  let debug n msg | n <= verb  = putStrLn msg
                  | otherwise  = return ()

  start_decoding <- getCurrentTime

  !(!mess1,parity) <- decode ecc rx

  debug 3 $ "decoded message"
  debug 4 $ "parity: " ++ show parity
  debug 4 $ show (Decoded (name ecc) mess1)

  when (U.length mess0 /= U.length mess1) $ do
    error $ "before and after codeword different lengths" ++ show (U.length mess0,U.length mess1)

  let !bitErrorCount = length [ () | (b,a) <- U.toList $ U.zip mess0 mess1, a /= b ]
  debug 2 $ show bitErrorCount ++ " bit errors in message"

  end_decoding <- getCurrentTime

  return $ TestRun encoding_time
                   (realToFrac $ diffUTCTime end_decoding start_decoding)
         $ eventBEs bitErrorCount

{-

    Sigma2 = max(( 1/10^(EbNoDB(jj)/10) )/2, 1e-9);
    Lc    = 2*sqrt(Ec)/Sigma2;

    for ii=1:NumBlocks(jj)

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Bits = 0.5 >= rand(1,k);
        CodedBits = mod(Bits*G,2);

        % Modulate the CodedBits
        s     = sqrt(Ec)*(2*CodedBits - 1);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % AWGN channel
        noise  = sqrt(Sigma2)*randn(size(s));
        r(1:n) = s + noise;
-}

-- Adding randomness
txRx_EbN0 :: EbN0 -> Rate -> GenIO -> U.Vector Bool -> IO (U.Vector Double)
txRx_EbN0 ebnoDB rate gen xs
        | isNaN ebnoDB = return $ U.map soft
                                $ xs
txRx_EbN0 ebnoDB rate gen xs = do
        rs :: U.Vector Double  <- U.fromList <$> sequence [ standard gen | _ <- U.toList xs ]
        return $ U.map (* lc)
               $ U.zipWith (+) (U.map (* sqrt sigma2) rs)
                               (U.map (* sqrt ec)
                                  $ U.map soft
                                  $ xs)
     where
         sigma2 = ((1/10) ** (ebnoDB/10)) / 2
         ec     = fromIntegral (numerator rate) / fromIntegral (denominator rate)
         lc     = 2 * sqrt ec / sigma2

-- eccMerge merges all the input files.
eccMerger :: IO ()
eccMerger = do
    args <- getArgs
    if null args
     then error $ "usage: <name> [-v<n>] [-b<n>] [-m<n>] [-l<file or dir/>] [-c] log/<file> log/<file>"
     else do
       let opts = parseOptions args $ defaultOptions "merge/"
--       print opts
       -- now we try load the files.
       logs <- sequence
         [ do bs <- LBS.readFile nm
              case CSV.decode CSV.HasHeader bs of
                Left msg -> fail msg
                Right t  -> return $ V.toList $ (t :: V.Vector Result)
         | nm <- codenames opts
         ]

       -- first, pick only the final results. We can not merge partual with final results.
       let sumMe :: Result -> Map String (Map Double Result) -> Map String (Map Double Result)
           sumMe r m = flip (Map.insert (res_name r)) m $
                   case Map.lookup (res_name r) m of
                     Nothing -> Map.singleton (res_ebN0 r) r
                     Just m' -> flip (Map.insert (res_ebN0 r)) m' $
                       case Map.lookup (res_ebN0 r) m' of
                         Nothing -> r
                         Just r' | res_size r >= res_size r' -> r
                                 | otherwise                 -> r'

       let logs' = map (foldr sumMe Map.empty) logs
       -- next, merge the different runs
       let mergeRes :: Result -> Result -> Result
           mergeRes r1 r2 = r1
              { timestamp = 0.0
              , res_size = res_size r1 + res_size r2
              , res_tEn  = res_tEn r1 + res_tEn r2
              , res_tDe  = res_tDe r1 + res_tDe r2
              , res_sum  = res_sum r1 + res_sum r2
              , res_ber  = Nothing
              , res_ldx  = Nothing
              , res_udx  = Nothing
              , res_conf = Nothing
              , res_bes  = res_bes r1 <> res_bes r2
              }
       let sumLog = Map.unionsWith (Map.unionWith mergeRes) logs'
--       print sumLog

       -- setup the printer

       -- Options -> [String] -> IO (ECC f -> EbN0 -> TestRun -> IO Bool)
       p <- eccPrinter opts (Map.keys sumLog)

       sequence_ [ do
         sequence_ [ p ecc db testrun
                   | db <- List.sort $ Map.keys m
                   , let r = m Map.! db
                     -- These are both fake, for the purposes of the log printing
                   , let ecc = ECC (res_name r) undefined undefined (res_mesg r) (res_code r)
                   , let testrun = TestRun (res_tEn r) (res_tDe r) (res_bes r)
                   ]
                 | (nm,m) <- Map.toList sumLog
                 ]

data Result = Result
 { timestamp  :: Double
 , res_name   :: String
 , res_ebN0   :: Double
 , res_size   :: Int     -- number of packets tested
 , res_tEn    :: Double
 , res_tDe    :: Double
 , res_sum    :: Int     -- total bit errors
 , res_ber    :: Maybe Double
 , res_ldx    :: Maybe Double
 , res_udx    :: Maybe Double
 , res_conf   :: Maybe Double
 , res_mesg   :: Int        -- length of message
 , res_code   :: Int        -- length of codeword
 , res_bes    :: BEs
 } deriving Show

instance CSV.ToRecord Result where
    toRecord (Result a b c d e f g h i j k l m n) = V.fromList [
        CSV.toField a, CSV.toField b, CSV.toField c, CSV.toField d, CSV.toField e, CSV.toField f,
        CSV.toField g, CSV.toField h, CSV.toField i, CSV.toField j, CSV.toField k, CSV.toField l,
        CSV.toField m, CSV.toField (show n)]

instance CSV.FromRecord Result where
    parseRecord v
        | n == 14    = Result        <$> CSV.unsafeIndex v 0
                                     <*> CSV.unsafeIndex v 1
                                     <*> CSV.unsafeIndex v 2
                                     <*> CSV.unsafeIndex v 3
                                     <*> CSV.unsafeIndex v 4
                                     <*> CSV.unsafeIndex v 5
                                     <*> CSV.unsafeIndex v 6
                                     <*> CSV.unsafeIndex v 7
                                     <*> CSV.unsafeIndex v 8
                                     <*> CSV.unsafeIndex v 9
                                     <*> CSV.unsafeIndex v 10
                                     <*> CSV.unsafeIndex v 11
                                     <*> CSV.unsafeIndex v 12
                                     <*> (read <$> CSV.unsafeIndex v 13)
        | otherwise = fail "Fail to read Result"
          where
            n = V.length v

writeHeader :: String -> IO ()
writeHeader logFileName =
   LBS.writeFile logFileName $ CSV.encode [
               [ "Time"
               , "ECC"
               , "EbN0"
               , "Packets"
               , "Encodes"
               , "Decodes"
               , "Errors"
               , "BER"
               , "LDX"
               , "UDX"
               , "Confidence"
               , "Message"
               , "Codeword"
               , "Internal"
               ]]

ext: hs

max_stars_repo_name: ku-fpg/ecc-manifold

max_issues_repo_name: ku-fpg/ecc-manifold

max_forks_repo_name: ku-fpg/ecc-manifold


text:
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE UndecidableInstances #-}
module FEECa.PolynomialTest where

import Control.Monad    ( liftM, liftM2 )

import FEECa.Internal.Vector
import FEECa.Internal.Simplex
import FEECa.Internal.Spaces
import FEECa.Polynomial

import FEECa.Utility.Utility
import FEECa.Utility.Combinatorics
import FEECa.Utility.Print

import Properties
import FEECa.Utility.Test
import qualified Test.QuickCheck as Q

-- =======
-- import qualified FEECa.Internal.MultiIndex as MI

--import qualified Numeric.LinearAlgebra.HMatrix as M
-- >>>>>>> efa221f... Cleaned up further tests (definitely shouldn't use `fromDouble`).

------------------------------------------------------------------------------
-- Dimension of the space to be tested. Must be greater than zero.
------------------------------------------------------------------------------

n :: Int
n = 2

------------------------------------------------------------------------------
-- Generate a random polynomial of dimension n as fixed by the above parameter.
-- The Polynomial may be of type Constant or Polynomial. If it is of type
-- Polynomial, degree the degree is chosen between 0 and 10 and the constants
-- are chosen randomly.
------------------------------------------------------------------------------

instance (Ring r, Q.Arbitrary r) => Q.Arbitrary (Polynomial r) where
  arbitrary = Q.frequency [(4,arbitraryPolynomial n), (1,arbitraryConstant)]

arbitraryPolynomial :: (Ring r, Q.Arbitrary r) => Int -> Q.Gen (Polynomial r)
arbitraryPolynomial n = do
  r <- Q.choose (0,10)
  s <- Q.choose (0,15)
  let mis = Q.vectorOf s (arbitraryMI n r)
      cs  = Q.vectorOf s Q.arbitrary
  liftM polynomial (liftM2 zip cs mis)

arbitraryConstant :: (Ring r, Q.Arbitrary r) => Q.Gen (Polynomial r)
arbitraryConstant = liftM constant Q.arbitrary

------------------------------------------------------------------------------
-- Generate random vectors of dimension n, so that they can be used to evaluate
-- the randomly generated polynomials.
------------------------------------------------------------------------------

instance (Field f, Q.Arbitrary f) => Q.Arbitrary (Vector f) where
  arbitrary = arbitraryVector n

------------------------------------------------------------------------------
-- Abstract propreties of arithmetic on polynomials. Addition, Subtraction and
-- multiplication of polynomials must commute with the evaluate operator.
------------------------------------------------------------------------------

propArithmetic :: (EuclideanSpace v, Function f v, VectorSpace f, Ring f,
                    r ~ Scalar v, r ~ Scalar f)
                => (r -> r -> Bool)
                -> f -> f -> v -> r
                -> Bool
propArithmetic eq x y v c =
     homomorphicEv addV add && homomorphicEv mul mul && homomorphicEv sub sub
    && prop_operator_commutativity eq (sclV c) (mul c) atV x
  where homomorphicEv o1 o2 = prop_homomorphism eq o1 o2 atV x y
        atV                 = evaluate v

------------------------------------------------------------------------------
-- Concrete arithmetic properties for polynomials defined over rationals.
------------------------------------------------------------------------------

prop_arithmetic_rf :: Polynomial Rational -> Polynomial Rational
                   -> Vector Rational -> Rational
                   -> Bool
prop_arithmetic_rf = propArithmetic (==)

--------------------------------------------------------------------------------
-- Derivation of Polynomials
--------------------------------------------------------------------------------

-- Linearity
propDerivation_linear :: (EuclideanSpace v, Function f v, VectorSpace f,
                          r ~ Scalar v, r ~ Scalar f)
                       => v -> v -> r -> f -> f
                       -> Bool
propDerivation_linear v1 v2 = prop_linearity (==) (evaluate v2 . derive v2)

-- Product rule
propDerivationProduct :: (EuclideanSpace v, Function f v, Ring f)
                        => v -> v -> f -> f
                        -> Bool
propDerivationProduct v1 v2 f g =
    ev (add (mul g (d f)) (mul f (d g))) == (ev . d) (mul f g)
  where ev = evaluate v1
        d  = derive v2

prop_derivation_product :: Vector Rational -> Vector Rational
                        -> Polynomial Rational -> Polynomial Rational
                        -> Bool
prop_derivation_product = propDerivationProduct

prop_derivation_productD :: Vector Double -> Vector Double
                         -> Polynomial Double -> Polynomial Double
                         -> Bool
prop_derivation_productD v1 v2 f g =
    ev (add (mul g (d f)) (mul f (d g))) `eqNum` (ev . d) (mul f g)
  where ev = evaluate v1
        d  = derive v2

-- Test for polynomials using exact arithmetic.
prop_arithmetic_rational :: Polynomial Rational -> Polynomial Rational
                         -> Vector Rational -> Rational
                         -> Bool
prop_arithmetic_rational = propArithmetic (==)

prop_derivation_linear_rational :: Vector Rational -> Vector Rational
                                -> Rational
                                -> Polynomial Rational -> Polynomial Rational
                                -> Bool
prop_derivation_linear_rational = propDerivation_linear

prop_derivation_product_rational :: Vector Rational -> Vector Rational
                                 -> Rational
                                 -> Polynomial Rational -> Polynomial Rational
                                 -> Bool
prop_derivation_product_rational = propDerivation_linear

--------------------------------------------------------------------------------
-- Barycentric Coordinates
--------------------------------------------------------------------------------

-- | Generate random simplex of dimesion 1 <= n <= 10.
instance (EuclideanSpace v, Q.Arbitrary (Scalar v)) => Q.Arbitrary (Simplex v) where
  arbitrary = do n <- Q.choose (1, 3)
                 k <- Q.choose (1, n)
                 arbitrarySubsimplex k n

-- TODO: Fails for n > 3 apparently due to numerical instability. To investigate
-- further.
prop_barycentric :: Simplex (Vector Double) -> Bool
prop_barycentric t =
    allEq [[evaluate v b | v <- vs] | b <- bs] oneLists
  where allEq l1 l2 = and $ zipWith (\l3 l4 -> (and (zipWith eqNum l3 l4))) l1 l2
        bs          = barycentricCoordinates t
        vs          = vertices t
        k           = topologicalDimension t
        oneLists    = map (map fromInt) (sumRLists (k+1) 1)

--------------------------------------------------------------------------------
-- Gradients of Barycentric Coordinates
--------------------------------------------------------------------------------

-- Helper function to create list of vector with 1.0 as first component
d0vectors :: Field a => Int -> [Vector a]
d0vectors n = [fromList $ [mulId]
                 ++ (replicate (i-1) addId)
                 ++ [addInv mulId]
                 ++ (replicate (n - i) addId) | i <- [1..n]]

-- Local gradients of barycentric coordinates, i.e. taken w.r.t the barycentric
-- coordinates themselves. Here we ensure that for a random simplex, if we move
-- from one corner a to corner b (in barycentric coordinates) the multiplication
-- of the  different vector with jacobian yields the corner b
-- (also in barycentric coordinates). This is equivalent to the vectors
-- [x_0, ..., x_n] with x_0 = 1 and x_i = -1 for any i > 0 being eigenvectors of
-- the gradient matrix.
prop_local_gradients :: Simplex (Vector Double) -> Bool
prop_local_gradients t = all (\v -> and $ zipWith eqNum (toList ((mult grads v)::Vector Double)) (toList v)) vs
  where grads     = localBarycentricGradients t
        vs        = d0vectors n :: [Vector Double]
        mult vs v = fromList [dot w v | w <- vs]
        n         = topologicalDimension t

return []
testPolynomial = $quickCheckWithAll

max_stars_repo_name: Airini/FEECa

max_issues_repo_name: Airini/FEECa

max_forks_repo_name: Airini/FEECa


text:
{- Assignment 1
 - Name: Omar Alkersh
 - Date: TODO add of completion
 -}
module Assign_1_ExtraCredit where

import           Data.Complex
-- see https://www.stackage.org/haddock/lts-8.24/base-4.9.1.0/Data-Complex.html

macid = "alkersho"

{- -----------------------------------------------------------------
 - cubicQ
 - -----------------------------------------------------------------
 - Description: Finds Q, refer to website
 -}
cubicQ :: Double -> Double -> Double -> Double
cubicQ a b c = ((3*a*c)-b^2)/(9*a^2)

{- -----------------------------------------------------------------
 - cubicR
 - -----------------------------------------------------------------
 - Description: Finds R using
 -}
cubicR :: Double -> Double -> Double -> Double -> Double
cubicR a b c d = ((9*a*b*c)-(27*a^2*d)-(2*b^3))/(54*a^3)

{- -----------------------------------------------------------------
 - cubicDisc
 - -----------------------------------------------------------------
 - Description: Finds the discriminant
 - when the disc > 0, there is one real solution
 - when the disc = 0, there are two real solutions
 - when the disc < 0, there are three real solutions
 -}
cubicDisc :: Double -> Double -> Double
cubicDisc q r = q^3 + r^2

--Returns the cube root of a number
cubeRoot :: Double -> Double
cubeRoot x = signum x * abs x ** (1/3)

{- -----------------------------------------------------------------
 - cubicS
 - -----------------------------------------------------------------
 - Description: Finds S using Q and R
 -}
cubicS :: Double -> Double -> Double
cubicS q r = cubeRoot(r+sqrt(q^3 + r^2))

cubicComplexS :: Double -> Double -> Complex Double
cubicComplexS q r = if cubicDisc q r >= 0
 then (cubicS q r :+ 0)
 else ((r + sqrt(-(q^3 + r^2))):+1)**(1/3)

{- -----------------------------------------------------------------
 - cubicT
 - -----------------------------------------------------------------
 - Description: Finds T using Q and R
 -}
cubicT :: Double -> Double -> Double
cubicT q r = cubeRoot(r-sqrt(q^3 + r^2))

cubicComplexT :: Double -> Double -> Complex Double
cubicComplexT q r = if cubicDisc q r >= 0
 then (cubicT q r :+ 0)
 else ((r - sqrt(-(q^3 + r^2))):+1)**(1/3)

{- -----------------------------------------------------------------
 - cubicRealSolutions
 - -----------------------------------------------------------------
 - Description: Finds the roots/solutions of the cubic equaition
 - when the disc > 0, there is one real solution
 - when the disc = 0, there are two real solutions
 - when the disc < 0, there are three real solutions
 -}
cubicRealSolutions :: Double -> Double -> Double -> Double -> [Double]
cubicRealSolutions a b c d =
  let
    q = cubicQ a b c
    r = cubicR a b c d
    s = cubicS q r
    t = cubicT q r
    x1 = s + t - (b / (3*a))
    x2 = -((s+t)/2)-(b/3*a)+((sqrt 3 / 2)*(s-t))
    tol = 1e-5
    in if cubicDisc q r > 0
     then [x1]
     else if abs (cubicDisc q r) == 0.0
       then [x1, x2, x2]
       else []

cubicComplexSolutions :: Double -> Double -> Double -> Double -> [Complex Double]
cubicComplexSolutions a b c d =
  let
    q = cubicQ a b c
    r = cubicR a b c d
    s = cubicComplexS q r
    t = cubicComplexT q r
    x1 = (realPart s + realPart t - (b / (3*a))):+ (imagPart s + imagPart t)
    fstTrm = (realPart s + realPart t)/2 :+ (imagPart s + imagPart t)/2 --returns a complex number
    sndTrm = b/(3*a) -- returns a real number
    sDt = (realPart s - realPart t) :+ (imagPart s - imagPart t)
    trdTrm = (-(((sqrt 3)/2)*imagPart sDt)) :+ ((sqrt 3)/2)*realPart sDt-- returns a complex number
    x2 = (-(realPart fstTrm) - sndTrm + (realPart trdTrm)):+( -(imagPart fstTrm) + (imagPart trdTrm))
    x3 = (-realPart fstTrm - sndTrm - realPart trdTrm):+( - imagPart fstTrm - imagPart trdTrm)
    in [x1, x2, x3]


text:
{-# LANGUAGE DeriveGeneric #-}

import qualified Network as N
import Network.HTTP.Client
import Network.Socket

import System.IO
import System.Environment

import Control.Concurrent
import Control.Monad

import qualified Data.ByteString.Lazy.Char8 as BL
import Data.Aeson (FromJSON, ToJSON, decode, encode)
import Data.String

import GHC.Generics (Generic)

import qualified ArgumentsParser as AP
import qualified Data.Text as DT
import qualified GetWordsFile as GWF
import qualified PorterStemmer as PS
import qualified StatisticsSlave as Stat
import qualified KeyWordOccurrence as KWO

main :: IO ()
main = do
    ip <- getArgs
    let ip' = ip !! 0
    print $ "Connecting to Master @" ++ ip' ++ ":4444"
    handle <- N.connectTo ip' (N.PortNumber 4444)
    putStrLn "Connected"
    incomingOrder <- newChan
    outgoingOrder <- newChan
    socket <- N.listenOn (N.PortNumber 5000)
    forkIO $ receiveOrder incomingOrder socket
    forkIO $ executeOrder incomingOrder outgoingOrder ip'
    forkIO $ forever $ sendWork outgoingOrder ip'
    forever $ do
        putStrLn "Slave is active..."
        threadDelay 500000

-- |Receives and broadcasts next order on incomingOrder channel
receiveOrder    :: Chan [[String]]  -- ^ Channel used to broadcast incoming order
                -> N.Socket         -- ^ Socket used to accept orders from Master
                -> IO ()
receiveOrder incomingOrder socket = do
    (handle, host, port) <- N.accept socket
    putStrLn "Receiving orders..."
    code <- hGetLine handle
    let order = read code :: [[String]]
    print $ "Order received : " ++ show order
    writeChan incomingOrder order
    hClose handle
    receiveOrder incomingOrder socket

-- | Reads and executes next order
executeOrder    :: Chan [[String]]  -- ^ Channel used to broadcast incoming order
                -> Chan String      -- ^ Channel used to transfer completed order to socket used for sending work back to Master
                -> String           -- ^ IP address of Master
                -> IO ()
executeOrder incomingOrder outgoingOrder ip = do
    order <- readChan incomingOrder
    let parsedArguments = map (AP.parseArguments ip) order
    let urls = map extractURLs parsedArguments
    putStrLn "Attempting to download all contents"
    contents <- getAllContent urls
    putStrLn "All descriptions downloaded"
    let work = process parsedArguments contents
    writeChan outgoingOrder work
    executeOrder incomingOrder outgoingOrder ip

-- |Sends the work done to the master
sendWork    :: Chan String  -- ^ Channel where completed work is found
            -> String       -- ^ IP address of Master
            -> IO ()
sendWork outgoingOrder ip = do
    work <- readChan outgoingOrder
    handle <- N.connectTo ip (N.PortNumber 4446)
    putStrLn "Attempting to send work back"
    hPutStrLn handle work
    putStrLn "Work sent back"
    hClose handle

-- |Gets all forth elements of a list of tuples of 4 elements
extractURLs :: [(String, Int, Int, String)] -- ^ A list of tuples like [(keyword, season, episode, url)]
            -> [String]                     -- ^ A list of String : [url]
extractURLs [] = []
extractURLs (x:xs) = (getURL x) : (extractURLs xs)

-- |Gets the forth element of a tuple of 4 elements
getURL  :: (String, Int, Int, String)   -- ^ A tuple like (keyword, season, episode, url)
        -> String                       -- ^ A String : url
getURL (_,_,_,u) = u

-- |Executes GET request
get :: String       -- ^ The url to be downloaded
    -> IO String    -- ^ The content of the url
get url = do
    manager <- newManager defaultManagerSettings
    request <- parseRequest url
    response <- httpLbs request manager
    return $ BL.unpack $ responseBody response

-- |Executes a list of lists of GET requests
getAllContent   :: [[String]]       -- ^ A list of lists of url to be downloaded
                -> IO [[String]]    -- ^ A list of lists of contents of url
getAllContent [] = return []
getAllContent (s:sx) = do
    c <- getContent s
    cx <- getAllContent sx
    return (c:cx)

-- |Executes a list of GET requests
getContent  :: [String]     -- ^ A list of url to be downloaded
            -> IO [String]  -- ^ A list of contents of url
getContent [] = return []
getContent (e:ex) = do
    c <- get e
    cx <- getContent ex
    return (c:cx)

-- |Exchanges the forth element of a list of lists of tuples with an element of another list of lists
updateAllContent    :: [[(String, Int, Int, String)]]   -- ^ A list of lists of tuples [[(keyword, season, episode, url)]]
                    -> [[String]]                       -- ^ A list of lists of String [[body]]
                    -> [[(String, Int, Int, String)]]   -- ^ A list of lists which is exchanged the forth element with the other list of lists [[(keyword, season, episode, body)]]
updateAllContent [] _ = []
updateAllContent (x:xs) (y:ys) = (updateContent x y) : (updateAllContent xs ys)

-- |Exchanges the forth element of a list of tuples with an element of another list
updateContent   ::  [(String, Int, Int, String)]    -- ^ A list of tuples [(keyword, season, episode, url)]
                -> [String]                         -- ^ A list of String [body]
                -> [(String, Int, Int, String)]     -- ^ A list which is exchanged the forth element with the other list [(keyword, season, episode, body)]
updateContent [] [] = []
updateContent ((a,b,c,url):xs) (body:ys) = (a, b, c, body) : (updateContent xs ys)

-- |Returns the first element of the head of a list of lists of tuples
getKeyword  :: [[(String, Int, Int, String)]]   -- ^ A list of lists of tuples like [[(keyword, season, episode, url)]]
            -> String                           -- ^ The first element of the tuples : keyword
getKeyword [] = ""
getKeyword (((keyword,_,_,_):_):_) = keyword

-- |Returns a process depending on a keyword ("N/A" returns the number of occurrences of every words, a keyword returns the number of occurrences of this word)
process :: [[(String, Int, Int, String)]]   -- ^ A list of lists of tuples like [[(keyword, season, episode, url)]]
        -> [[String]]                       -- ^ A list of lists with the text to stem [[a JSON content]]
        -> String                           -- ^ a JSON content containing the results of the process depending on the keyword
process parsedArguments descriptions = do
    let updatedContent = updateAllContent parsedArguments descriptions
    if getKeyword updatedContent == "N/A"
        then show $ statisticProcess updatedContent
        else show $ map processEachSeason updatedContent

-- |Returns a list of list of String with the season number, the episode number, and a list of pairs ("word", occurrence number)
statisticProcess    :: [[(String, Int, Int, String)]]   -- ^ A list of lists of tuples with the text to stem [[("N/A", season, episode, a JSON content)]]
                    -> [[String]]                       -- ^ A list like [["1,1,[("a",1),("b",2)]","1,2,[("d",1)]"],["2,1,[("f",1),("g",2)]","2,2,[("i",1)]"]]
statisticProcess l = Stat.statisticsEachSlave $ statisticProcessAllSeason l

-- |On a list of lists, deletes the stopwords of every text, applies the Stemming algorithm and returns all stemmed words
statisticProcessAllSeason   :: [[(String, Int, Int, String)]]   -- ^ A list of lists of tuples with the text to stem [[("N/A", season, episode, a JSON content)]]
                            -> [[(String, Int, Int, String)]]   -- ^ A list of lists of tuples with all stemmed words [[("N/A", season, episode, text composed of stemmed words)]]
statisticProcessAllSeason l = map statisticProcessEachSeason l

-- |On a list, deletes the stopwords of every text, applies the Stemming algorithm and returns all stemmed words
statisticProcessEachSeason  :: [(String, Int, Int, String)] -- ^ A list of tuples with the text to stem [("N/A", season, episode, a JSON content)]
                            -> [(String, Int, Int, String)] -- ^ A list of tuples with all stemmed words [("N/A", season, episode, text composed of stemmed words)]
statisticProcessEachSeason l = map statisticProcessEachEpisode l

-- |Deletes the stopwords of a text, applies the Stemming algorithm and returns all stemmed words
statisticProcessEachEpisode :: (String, Int, Int, String)   -- ^ A tuple with the text to stem ("N/A", season, episode, a JSON content) like ("N/A", 1, 1, text)
                            -> (String, Int, Int, String)   -- ^ A tuple with all stemmed words ("N/A", season, episode, text composed of stemmed words) like ("N/A", 1, 1, text)
statisticProcessEachEpisode (a,b,c,text) =
    let gwf = GWF.getWordsFile text in
    let ps = DT.unpack . DT.toLower . DT.pack $ PS.porterStemmer gwf in -- "jon\nhand\nking\n...\narmi\n"
    (a,b,c,ps)

data Serie = Serie { season :: Int, name :: String, episode :: Int, plot :: String }
    deriving (Show, Generic)

instance FromJSON Serie
instance ToJSON Serie

-- |On a list, extracts a data, takes a part of a JSON content and returns a list with the number of occurrence of a keyword in a text with the Stemming algorithm
processEachSeason   :: [(String, Int, Int, String)] -- ^ A list of tuples (keyword, season, episode, a JSON content) like ("jon", 1, 1, text)
                    -> [[Int]]                      -- ^ A list of lists composed of 3 Int [season, episode, occurrence] like [1,1,7]
processEachSeason [] = []
processEachSeason (x:xs) = (processEachEpisode x) : (processEachSeason xs) -- [("jon",1,1,"jon"),("jon",1,1,"jon jon"),("jon",1,1,"jon jon jon")]

-- |Converts body to a data and applies the Stemming algorithm on the plot
processEachEpisode  :: (String, Int, Int, String)   -- ^ A tuple (keyword, season, episode, a JSON content) like ("jon", 1, 1, text)
                    -> [Int]                        -- ^ A list composed of 3 Int [season, episode, occurrence] like [1,1,7]
processEachEpisode (keyword, season, episode, body) = do
    let req = decode $ BL.pack body :: Maybe Serie
    case req :: Maybe Serie of
        Nothing  -> []
        Just req -> stemmingProcess (keyword, season, episode, (show.plot) req)

-- |Suppresses stopwords, stems the keyword and the text, and returns a triplet with the two Int and the number of occurrence of the first argument in the text of the forth argument
stemmingProcess :: (String, Int, Int, String)   -- ^ A tuple (keyword, season, episode, text) like ("jon", 1, 1, "Jon the Hand of the King ... an army")
                -> [Int]                        -- ^ A list composed of 3 Int [season, episode, occurrence] like [1,1,7]
stemmingProcess (keyword, season, episode, text) =
    let gwf = GWF.getWordsFile text in
    let kw = (DT.unpack . DT.toLower . DT.pack) (lines (PS.porterStemmer keyword) !! 0) in
    let ps = DT.unpack . DT.toLower . DT.pack $ PS.porterStemmer gwf in
    [season,episode,KWO.keyWordOccurrence kw ps]

max_stars_repo_name: thomaslecointre/Distributed-Haskell

max_issues_repo_name: thomaslecointre/Distributed-Haskell

max_forks_repo_name: thomaslecointre/Distributed-Haskell


text:
-- | Data structure of compiler. https://en.wikibooks.org/wiki/Write_Yourself_a_Scheme_in_48_Hours/

module LispVal
       ( LispVal(..)
       , Env
       , unwordsList
       ) where

import           Data.Complex
import           Data.IORef
import           Data.Ratio                 (denominator, numerator)
import           ErrorCheckingAndExceptions
import           System.IO                  (Handle)

data LispVal = Atom String
             | List [LispVal]
             | DottedList [LispVal] LispVal
             | Number Integer
             | Float Double
             | Ratio Rational
             | Complex (Complex Double)
             | Character Char
             | String String
             | Bool Bool
             | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
             | Func { params :: [String], vararg :: Maybe String, body :: [LispVal], closure :: Env}
             | IOFunc ([LispVal] -> IOThrowsError LispVal)
             | Port Handle

showVal :: LispVal -> String
showVal (Atom x) = x
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (Ratio x) = show (numerator x) ++ ('/' : show (denominator x))
showVal (Complex x) = show (realPart x) ++ ('+' : show (imagPart x) ++ "i")
showVal (Float x) = show x
showVal (Number x) = show x
showVal (Character x) = "#\\" ++ show x
showVal (String x) = "\"" ++ x ++ "\""
showVal (List x) = "(" ++ unwordsList x ++ ")"
showVal (DottedList first rest)= "(" ++ unwordsList first ++ " . " ++ showVal rest ++ ")"
showVal (PrimitiveFunc _) = "<primitive>"
showVal (Func {params = args, vararg = varargs, body = body, closure = env}) =
  "(lambda (" ++ unwords (map show args) ++
  (case varargs of
     Nothing -> ""
     Just arg -> " . " ++ arg) ++ ") ...)"
showVal (Port _) = "<IO port>"
showVal (IOFunc _) = "<IO primitive>"

unwordsList :: [LispVal] -> String
unwordsList = unwords . map showVal

instance Show LispVal where
  show = showVal

type Env = IORef [(String, IORef LispVal)]

max_stars_repo_name: zeqing-guo/Haskell-Exercises

max_issues_repo_name: zeqing-guo/Haskell-Exercises

max_forks_repo_name: zeqing-guo/Haskell-Exercises


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

module Main where
import Debug.Trace (traceShow)
import Data.Either.Combinators (fromRight, maybeToRight)
import BDSCOD.Llhd (llhdAndNB,initLlhdState, intervalLlhd)
import BDSCOD.Types
import BDSCOD.Utility
import Control.Monad (zipWithM)
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 qualified Data.Csv as Csv
import Data.List (intercalate, intersperse,nub)
import Data.Maybe (fromMaybe)
import qualified Data.Vector.Unboxed as Unboxed
import qualified Epidemic.BDSCOD as SimBDSCOD
import Epidemic.Types.Events
  ( EpidemicEvent(..)
  , asNewickString
  , eventTime
  , maybeEpidemicTree
  , maybeReconstructedTree
  )
import Epidemic.Types.Parameter
import Epidemic.Types.Population (Person(..), Identifier(..))
import qualified Epidemic.Utility as SimUtil
import GHC.Generics
import Numeric.Minimisation.Powell (minimise)
import Numeric.LinearAlgebra.Data (linspace, toList)
import Numeric.LinearAlgebra.HMatrix
import System.Environment (getArgs)
import System.Random.MWC
import Data.Word (Word32)
import System.IO.Unsafe (unsafePerformIO)

-- | Values of this type are used to specify an analysis of the data.
data InferenceConfiguration =
  InferenceConfiguration
    {
    -- | The time in the simulation at which the inference is carried out.
    inferenceTime :: AbsoluteTime
    ,
    -- | A couple of files to write the Newick representation to and the node labels.
    reconstructedTreeOutputFiles :: (FilePath, FilePath)
    ,
    -- | Where to write the sequence of observed events.
    observationsOutputCsv :: FilePath
    ,
    -- | Where to write the likelihood evaluations.
    llhdOutputCsv :: FilePath
    ,
    -- | Where to write the point estimate of prevalence.
    pointEstimatesCsv :: FilePath
    }
  deriving (Show, Generic)

-- | These objects describe the evaluation mesh when looking at the likelihood
-- profiles. This is useful because it allows us to provide these details at
-- runtime rather than hardcoding them here.
data LlhdProfileMesh =
  LlhdProfileMesh
  { lpmMeshSize :: Int
  , lpmLambdaBounds :: (Rate,Rate)
  , lpmMuBounds :: (Rate,Rate)
  , lpmPsiBounds :: (Rate,Rate)
  , lpmRhoBounds :: (Probability,Probability)
  , lpmOmegaBounds :: (Rate,Rate)
  , lpmNuBounds :: (Probability,Probability)
  } deriving (Show, Generic)

-- | This object configures the whole evaluation of this program and is to be
-- read in from a suitable JSON file.
data AppConfiguration =
  AppConfiguration
    { simulatedEventsOutputCsv :: FilePath
    , simulationParameters :: Parameters
    , simulationDuration :: TimeDelta
    , simulationSizeBounds :: (Int,Int)
    , simulationSeed :: Word32
    , inferenceConfigurations :: [InferenceConfiguration]
    , acLlhdProfileMesh :: LlhdProfileMesh
    }
  deriving (Show, Generic)

instance Json.FromJSON AppConfiguration

instance Json.FromJSON InferenceConfiguration

instance Json.FromJSON LlhdProfileMesh

type Simulation x = ReaderT AppConfiguration (ExceptT String IO) x

-- | This type is used to indicate if parameters are the true ones used in the
-- simulation or estimates parameters.
data ParameterKind
  = SimulationParameters
  | EstimatedParameters
  deriving (Show, Eq)

-- | A BDSCOD simulation configuration based on the parameters in the
-- environment.
bdscodConfiguration = do
  simParams <- asks simulationParameters
  (TimeDelta simDurDouble) <- asks simulationDuration
  let bdscodConfig =
        SimBDSCOD.configuration (AbsoluteTime simDurDouble) (unpackParameters simParams)
  case bdscodConfig of
    Just x -> do liftIO $ putStrLn (show simParams)
                 return x
    Nothing -> throwError "Could not construct BDSCOD configuration"

-- | Simulate the transmission process part of the epidemic making sure that the
-- results are acceptable in terms of the number of observed events before
-- returning a filtration of the events, i.e., the data that was availble at
-- several points in time.
--
-- __NOTE__ the filteration must happen before these are processed into
-- observations since the observations do not accumulate chronologically due to
-- birth events which can occur in the past due to new observations in the
-- present.
partialSimulatedEpidemic seedWord bdscodConfig =
  do
    gen <- liftIO $ initialize (Unboxed.fromList [seedWord])
    simEvents <- liftIO $ SimUtil.simulation' bdscodConfig SimBDSCOD.allEvents gen
    (sizeLowerBound,sizeUpperBound) <- asks simulationSizeBounds
    if length simEvents > sizeLowerBound && length simEvents < sizeUpperBound
      then do infTimes <- map inferenceTime <$> asks inferenceConfigurations
              simEventsCsv <- asks simulatedEventsOutputCsv
              liftIO $ L.writeFile simEventsCsv (Csv.encode simEvents)
              return [filter (\e -> eventTime e <= infTime) simEvents | infTime <- infTimes]
      else do liftIO $ putStrLn "Repeating epidemic simulation..."
              partialSimulatedEpidemic (seedWord + 1) bdscodConfig

-- | Run the actual observation of the simulation and record the results before
-- returning the dataset of observations generated by this epidemic.
simulatedObservations :: InferenceConfiguration
                      -> [EpidemicEvent]
                      -> Simulation (InferenceConfiguration,[Observation])
simulatedObservations infConfig@InferenceConfiguration{..} simEvents = do
  let Just (newickBuilder,newickMetaData) =
        do eTree <- maybeEpidemicTree simEvents
           rTree <- maybeReconstructedTree eTree
           asNewickString (AbsoluteTime 0, Person (Identifier 1)) rTree
      maybeObs = eventsAsObservations <$> SimBDSCOD.observedEvents simEvents
      (reconNewickTxt,reconNewickCsv) = reconstructedTreeOutputFiles
  case maybeObs of
    (Just obs) ->
      do
        liftIO $ L.writeFile reconNewickTxt (BBuilder.toLazyByteString newickBuilder)
        liftIO $ L.writeFile reconNewickCsv (Csv.encode newickMetaData)
        liftIO $ L.writeFile observationsOutputCsv (Csv.encode obs)
        return (infConfig,obs)
    _ -> throwError "Failed to simulate observations."

-- | If there is a unique timed value return that. This is used to make it
-- easier to extract the values of timed parameters so that you don't need to
-- store all of them.
uniqueTimedValue :: Eq x => Timed x -> Maybe x
uniqueTimedValue (Timed txs) = case txs of
  [] -> Nothing
  txs' -> let xs = nub [snd tx | tx <- txs']
              isUnique = 1 == length xs
            in if isUnique then Just (head xs) else Nothing

-- | Evaluate the NB posterior approximation of the prevalence for a single
-- point in parameter space and the LLHD over a list of points and write all of
-- the results to CSV.
recordLlhdCrossSections :: InferenceConfiguration
                          -> [Observation]
                          -> (Parameters,ParameterKind,[Parameters])
                          -> Simulation ()
recordLlhdCrossSections InferenceConfiguration {..} obs (singleParams, paramKind, evalParams) =
  let comma = BBuilder.charUtf8 ','
      parametersUsed = show paramKind
      parametersUsed' = BBuilder.stringUtf8 parametersUsed
      llhdVals = safeLlhdFromInit obs <$> evalParams
      nBValAndParams =
        Csv.encode . (:[]) $
        ( parametersUsed
        , show $ length obs
        , fromRight Zero $ snd <$> llhdAndNB obs singleParams initLlhdState
        , show $ getLambda singleParams
        , show $ getMu singleParams
        , show $ getPsi singleParams
        , show <$> uniqueTimedValue $ getRhos singleParams
        , show $ getOmega singleParams
        , show <$> uniqueTimedValue $ getNus singleParams )
      doublesAsString =
        BBuilder.toLazyByteString .
        mconcat .
        intersperse comma . (parametersUsed' :) . map BBuilder.doubleDec
   in liftIO $
      do
        L.appendFile llhdOutputCsv (doublesAsString llhdVals)
        L.appendFile pointEstimatesCsv nBValAndParams

-- | Evaluate the LLHD function on cross-sections centred at the true values
evaluateLLHD :: InferenceConfiguration -> [Observation] -> Simulation ()
evaluateLLHD infConfig obs = do
  simParams <- asks simulationParameters -- get the actual parameters used to simulate the observations
  llhdProfMesh <- asks acLlhdProfileMesh -- this specifies the cross section to use.
  let evalParams = crossSectionParameters llhdProfMesh simParams
  recordLlhdCrossSections infConfig obs (simParams,SimulationParameters,evalParams)

-- | Estimate the parameters of the of the model and then evaluate the LLHD
-- profiles and prevalence and append this to the file. The first value of the
-- CSV output now describes which parameters where used to to evaluate these
-- things.
estimateLLHD :: InferenceConfiguration -> [Observation] -> Simulation ()
estimateLLHD infConfig obs = do
  simParams@(Parameters (_,deathRate,_,_,_,_)) <- asks simulationParameters
  llhdProfMesh <- asks acLlhdProfileMesh
  liftIO $ putStrLn $ "\t\tEstimating parameters using " ++ show (length obs) ++ " observations..."
  let schedTimes = scheduledTimes simParams
      mleParams = estimateParameters deathRate schedTimes obs -- get the MLE estimate of the parameters
      evalParams = crossSectionParameters llhdProfMesh mleParams
  recordLlhdCrossSections infConfig obs (mleParams,EstimatedParameters,evalParams)

-- | A safe version of the log-likelihood assuming that you want to start from a
-- simple initial condition.
--
-- __NOTE__ This conditions the process against extinction.
--
safeLlhdFromInit :: [Observation] -> Parameters -> LogLikelihood
safeLlhdFromInit obs p@(Parameters (λ, μ, ψ, _, ω, _)) =
    let
      -- log-likelihood
      llhd = fromRight (-1e6) $ fst <$> llhdAndNB obs p initLlhdState
      -- log of the probability of not going extinct.
      lpne = log (1 - (μ + ψ + ω) / λ)
    in llhd - lpne

-- | Estimate of the MLE. This uses a simplex method.
--
-- __NOTE__ the ugly case statement means that this should work both with and
-- without scheduled observations.
--
-- __NOTE__ we fix the death rate to the true value because
-- this is assumed to be known a priori.
--
estimateParameters :: Rate -> ([AbsoluteTime],[AbsoluteTime]) -> [Observation] -> Parameters
estimateParameters deathRate sched obs =
  let energyFunc v2p x =
        negate $ safeLlhdFromInit obs (v2p x)
      mini rI v2p = minimise (energyFunc v2p) rI
  in case sched of
       -- there are no scheduled observations
       ([],[]) ->
         let -- randInit = [-1.5,-2.2,-2.2] -- initial point to start
           randInit = log <$> [0.228,4.8e-2,2.6e-2] -- _GROSS_ initial point as true value
           vec2Param vec =
             let [lnR1, lnR2, lnR3] = vec
             in packParameters ( exp lnR1
                               , deathRate
                               , exp lnR2
                               , []
                               , exp lnR3
                               , [])
         in case mini randInit vec2Param of
              Right (est,_,grad) -> traceShow grad $ vec2Param est
              Left msg -> error msg
       -- there is at least one scheduled observation
       (rhoTs, nuTs) ->
         let randInit = [-1.5,-2.5,-1,-2.5,-1]
             vec2Param vec =
               let [lnR1, lnR2, logitP1, lnR3, logitP2] = vec
                   rhoVal = invLogit logitP1
                   nuVal = invLogit logitP2
                   timed v ts = zip ts (repeat v)
               in packParameters ( exp lnR1
                                 , deathRate
                                 , exp lnR2
                                 , timed rhoVal rhoTs
                                 , exp lnR3
                                 , timed nuVal nuTs)
         in case mini randInit vec2Param of
              Right (est,_,grad) -> traceShow grad $ vec2Param est
              Left msg -> error msg

-- | List of parameters required to plot the cross sections.
crossSectionParameters :: LlhdProfileMesh -> Parameters -> [Parameters]
crossSectionParameters LlhdProfileMesh{..} ps =
  let mesh = toList . linspace lpmMeshSize
      lambdaMesh = mesh lpmLambdaBounds
      muMesh = mesh lpmMuBounds
      psiMesh = mesh lpmPsiBounds
      rhoProbMesh = mesh lpmRhoBounds
      omegaMesh = mesh lpmOmegaBounds
      nuProbMesh = mesh lpmNuBounds
      (rhoTimes,nuTimes) = scheduledTimes ps
      rhoMesh = [Timed [(t,r) | t <- rhoTimes] | r <- rhoProbMesh]
      nuMesh = [Timed [(t,n) | t <- nuTimes] | n <- nuProbMesh]
      -- for each dimension, construct a list of parameter values with the
      -- values from the mesh
      apply f = map (f ps)
      [lPs,mPs,pPs,oPs] =
        zipWith apply [putLambda,putMu,putPsi,putOmega] [lambdaMesh,muMesh,psiMesh,omegaMesh]
      [rPs,nPs] = zipWith apply [putRhos,putNus] [rhoMesh,nuMesh]
  in concat $
     if (rhoTimes,nuTimes) /= ([],[])
     then [lPs,mPs,pPs,rPs,oPs,nPs]
     else [lPs,mPs,pPs,oPs]

-- | Definition of the complete simulation study.
simulationStudy :: Simulation ()
simulationStudy = do
  bdscodConfig <- bdscodConfiguration
  liftIO $ putStrLn "\tRunning epidemic simulation"
  simSeed <- asks simulationSeed
  pEpi <- partialSimulatedEpidemic simSeed bdscodConfig
  infConfigs <- asks inferenceConfigurations
  liftIO $ putStrLn "\tExtracting observations from full simulation"
  pObs <- zipWithM simulatedObservations infConfigs pEpi
  liftIO $ putStrLn "\tEvaluating LLHD on cross-sections about true parameters"
  mapM_ (uncurry evaluateLLHD) pObs
  liftIO $ putStrLn "\tEvaluating LLHD on cross-sections about estimated parameters"
  mapM_ (uncurry estimateLLHD) pObs

main' :: IO ()
main' = do
  let configFilePath = "./examples/simulation-study-time-series/ts-config.json"
  maybeConfig <- getConfiguration configFilePath
  case maybeConfig of
    Nothing ->
      putStrLn $ "Could not get configuration from file: " ++ configFilePath
    Just config -> do
      putStrLn $ "Succeeded in reading configuration from file: " ++ configFilePath
      result <- runExceptT (runReaderT simulationStudy config)
      case result of
        Left errMsg -> putStrLn errMsg
        Right _ -> return ()

main :: IO ()
main = do
  configFilePath <- head <$> getArgs
  maybeConfig <- getConfiguration configFilePath
  case maybeConfig of
    Nothing ->
      putStrLn $ "Could not get configuration from file: " ++ configFilePath
    Just config -> do
      putStrLn $ "Succeeded in reading configuration from file: " ++ configFilePath
      result <- runExceptT (runReaderT simulationStudy config)
      case result of
        Left errMsg -> putStrLn errMsg
        Right _ -> return ()

-- | Attempt to read a configuration object from the given filepath.
getConfiguration :: FilePath -> IO (Maybe AppConfiguration)
getConfiguration fp = Json.decode <$> L.readFile fp

max_stars_repo_name: aezarebski/timtam

max_issues_repo_name: aezarebski/timtam

max_forks_repo_name: aezarebski/timtam


text:
module Dna.Search (
    dnaFromSepll, composingKMer,
    initDeBrujinGraphFromKMers,
    initDeBrujinGraphFromDna, fibrabbit) where

import Dna.Graph
import Dna.Nucleotides
import Data.List
import Control.Applicative
import Data.Maybe
import Control.Monad
import Numeric.LinearAlgebra.Devel

composingKMer ::  Int -> Dna -> DnaList
composingKMer n dna
    | lengthDna dna >= n = takeDna n dna : composingKMer n (tailDna dna)
    | otherwise = []

same one two = one `isSubsequenceOf` two && two `isSubsequenceOf` one

mergeHeadTail [] r = r
mergeHeadTail (lx:lxs) r@(rx:rxs)
    | lx == rx = r
    | otherwise = lx: mergeHeadTail lxs r

isConnected [] _ = False
isConnected _ [] = False
isConnected (_:lxs) r = lxs `isPrefixOf` r

dnaFromSepll :: DnaList -> Dna
dnaFromSepll [] = nullDna
dnaFromSepll [dna] = dna
dnaFromSepll (x:xs) = dnaSpellStep 1 x xs
    where
        dnaSpellStep dropStep growDna dnas = case find (isDnaPrefixOf $ dropDna dropStep growDna) dnas of
            Just newDna -> dnaSpellStep (dropStep+1) (takeDna dropStep growDna +:+ newDna) dnas
            Nothing -> growDna

initDeBrujinGraphFromDna :: Int -> Dna -> Graph Dna
initDeBrujinGraphFromDna n dnas = initOverlaGraph nodes $ \ m findIndex ->
        forM_ adjacentNodes $ \ (from, to) ->
            modifyMatrix m (findIndex from) (findIndex to) (1+)
    where
        nodes = composingKMer n dnas
        adjacentNodes = zip (init nodes) (tail nodes)

initDeBrujinGraphFromKMers :: DnaList -> Graph Dna
initDeBrujinGraphFromKMers kmers = initOverlaGraph nodes $ \m findIdxNode ->
    forM_ nodes $ \from ->
        forM_ (map tailDna (filter (from `isDnaPrefixOf`) kmers)) $ \to ->
            modifyMatrix m (findIdxNode from) (findIdxNode to) (1+)
        where nodes = nub $ map (liftDna init) kmers ++ map tailDna kmers

fibrabbit 1 _ = 1
fibrabbit 2 _ = 1
fibrabbit n k = fibrabbit (n-1) k + (k * fibrabbit (n-2) k)

max_stars_repo_name: mathk/NautilusDna

max_issues_repo_name: mathk/NautilusDna

max_forks_repo_name: mathk/NautilusDna


text:
{-# OPTIONS -fglasgow-exts #-}

import GHC.Prim;
import Data.Complex;
import Data.Array;
import GHC.Base;

type Complex_type = Complex Double;
type Array_type b = Array Int b;
type Assoc_type a = (Int,a);
type Descr_type = (Int,Int);

w2i x = word2Int# x;
i2w x = int2Word# x;

abortstr str = error ("abort:" ++ str);

delay x = abortstr "delay not implemented";

fix :: (x -> x) -> x;
fix f = fix_f where {fix_f = f fix_f};

force x = x; -- error "force not implemented"

iff :: Bool -> x -> x -> x;
iff b x y = if b then x else y;

iffrev :: x -> x -> Bool -> x;
iffrev y x b = if b then x else y;

miraseq :: x -> y -> y;
miraseq x y = seq_const y x;  -- x should be marked #STRICT
seq_const x y = x;

pair :: [x] -> Bool;
pair [] = False;
pair x = True;

entier :: Double -> Double;
entier x = fromIntegral (floor x);

land_i :: Int -> Int -> Int;
land_i (I# x) (I# y) = I# (w2i (and# (i2w x) (i2w y)));

lnot_i :: Int -> Int;
lnot_i (I# x) = I# (w2i (not# (i2w x)));

lor_i :: Int -> Int -> Int;
lor_i (I# x) (I# y) = I# (w2i (or# (i2w x) (i2w y)));

lshift_i :: Int -> Int -> Int;
lshift_i (I# x) (I# y) = I# (w2i (shiftL# (i2w x) y));

rshift_i :: Int -> Int -> Int;
rshift_i (I# x) (I# y) = I# (w2i (shiftRL# (i2w x) y));

write x = abortstr "write not implemented";

descr :: Int -> Int -> Descr_type;
descr l u = (l,u);

destr_update :: Array_type x -> Int -> x -> Array_type x;
destr_update ar i x = ar // [(i,x)];

indassoc :: Assoc_type x -> Int;
indassoc (i,v) = i;

lowbound :: Descr_type -> Int;
lowbound (l,u) = l;

tabulate :: (Int -> x) -> Descr_type -> Array_type x;
tabulate f (l,u) = array (l,u) [(i, f i) | i <- [l..u]];

upbound :: Descr_type -> Int;
upbound (l,u) = u;

update :: Array_type x -> Int -> x -> Array_type x;
update ar i x = ar // [(i,x)];

valassoc :: Assoc_type x -> x;
valassoc (i,v) = v;


text:
module Statistics.Information.Mixed.MutualInfo where

import Data.List
import Data.Matrix
import qualified Data.Vector as Vector
import qualified Statistics.Information.Continuous.Entropy as CE
import Statistics.Information.Utils.List
import Statistics.Information.Utils.Random
import System.Random

micd :: (RandomGen g, Eq a) => g -> Int -> Int -> Matrix Double -> Matrix a ->
        Double
micd g k base xs ys | nrows xs == nrows ys =
  abs $ hx - sum [ycontrib ypg | ypg <- yvals]
  where
    (g1, g2) = split g
    hx = CE.entropy g1 k base xs
    ys' = toLists ys
    yvals = probs ys' `zip` splitN (length (nub ys')) g2
    ycontrib ((y, p), g') =
      let xgiveny = [getRow i xs | i <- [1..nrows xs],
                     (Vector.toList $ getRow i ys) == y] in
        if k <= length xgiveny - 1 then
          p * CE.entropy g' k base (fromLists (map Vector.toList xgiveny))
        else
          p * hx

midc :: (RandomGen g, Eq a) => g -> Int -> Int -> Matrix a -> Matrix Double ->
        Double
midc g k base xs ys = micd g k base ys xs

max_stars_repo_name: eligottlieb/Shannon

max_issues_repo_name: eligottlieb/Shannon

max_forks_repo_name: eligottlieb/Shannon


text:
{-# OPTIONS_GHC -fno-warn-dodgy-exports #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE NoImplicitPrelude #-}

module Sound
(
    -- * Reexports
    -- $reexports
    module Control.Monad
    , module Data.Complex
    , module Data.Int
    -- * Input, output, and conversion
    , module Sound.Io.Snd
    , module Sound.Io
    , module Sound.IoPtr
    , module Sound.Sox
    , module Sound.IoSox
    , module Sound.Portaudio
    -- * Abstraction ladder
    , module Sound.Class
    , module Sound.Pair
    , module Sound.Compile
    , module Sound.InfList
    , module Sound.List
    -- * Buffer
    , module Sound.Buffer
    , module Sound.Ptr
    -- * Time, tempo
    , module Sound.Time
    , module Sound.Tempo
    -- * Sampling
    , module Sound.Sample
    -- * Decibels, contours, envelopes, limiting, clamping
    , module Sound.Amplitude
    , module Sound.Ramp
    -- * Wavetable
    , module Sound.Table
    -- * White noise
    , module Sound.Random
    -- * Frequency modulation
    , module Sound.Fm
    -- * Integration
    , module Sound.Int
    -- * Instrument
    , module Sound.Perform
    -- * MIDI input
    , module Sound.Midi
    -- * Error handling, type hints
    , module Sound.Hint
    -- * Filters: zeros and poles
    , module Sound.Filter
    -- * Unclear
    , module Sound.Function
    -- * Failed experiments
    , module Sound.Abstract
    , module Sound.GenBody
    , module Sound.IoFail
    , module Sound.StreamVector
)
where

import Data.Complex
import Data.Int

import Control.Applicative
import Control.Monad

import Sound.Abstract
import Sound.Amplitude
import Sound.Buffer
import Sound.Class
import Sound.Compile
import Sound.Filter
import Sound.Fm
import Sound.Function
import Sound.GenBody
import Sound.Hint
import Sound.InfList
import Sound.Int
import Sound.Io
import Sound.Io.Snd ()
import Sound.IoFail
import Sound.IoPtr
import Sound.IoSox
import Sound.List ()
import Sound.Midi
import Sound.Pair
import Sound.Perform
import Sound.Portaudio ()
import Sound.Ptr
import Sound.Ramp
import Sound.Random
import Sound.Sample
import Sound.Sox hiding (Rate)
import Sound.StreamVector
import Sound.Table
import Sound.Tempo
import Sound.Time

import qualified Prelude as P

{- $reexports
If you use this module unqualified, you have to hide Prelude imports to avoid name clashes.
The portable way to do that is by adding this import to your code:

@
import Prelude ()
@

This module replaces 'P..' and 'P.id' from "Prelude"
with '.' and 'id' from "Control.Category".
This module hides 'P.seq' from "Prelude".

This module reexports everything else from all modules listed here.
-}

max_stars_repo_name: edom/work

max_issues_repo_name: edom/work

max_forks_repo_name: edom/work


text:
{-# LANGUAGE BangPatterns,
             ScopedTypeVariables,
             RecordWildCards,
             FlexibleContexts,
             TypeFamilies,
             GeneralizedNewtypeDeriving #-}

module AI.HsANN.ComplexNetwork
    (
    -- * Types
      ComplexNetwork(..)
    , ActivationFunction
    , ActivationFunctionDerivative
    , Sample
    , Samples
    , (-->)

    -- * Creating a neural network
    , createComplexNetwork
    , fromWeightMatrices

    -- * Computing a neural network's output
    , output
    , sigmoid
    , sigmoid'
    , complexSigmoid
    , complexSigmoid'
    , randComplex

    -- * Training a neural network
    , trainUntil
    , trainNTimes
    , trainUntilErrorBelow
    , quadError

    -- * Loading and saving a neural network
    , loadComplexNetwork
    , saveComplexNetwork
    ) where

import Codec.Compression.Zlib     (compress, decompress)
import Data.Binary                (Binary(..), encode, decode)
import Data.List                  (foldl')
import Foreign.Storable           (Storable)
import qualified Data.ByteString.Lazy  as B
import qualified Data.Vector           as V
import qualified Data.Complex as C
import Debug.Trace

import System.Random.MWC
import Numeric.LinearAlgebra
import Data.Functor ((<$>))

-- | Our feed-forward neural network type. Note the 'Binary' instance, which means you can use
--   'encode' and 'decode' in case you need to serialize your neural nets somewhere else than
--   in a file (e.g over the network)
newtype ComplexNetwork a = ComplexNetwork
                 { matrices   :: V.Vector (Matrix a) -- ^ the weight matrices
                 } deriving Show

instance (Element a, Binary a) => Binary (ComplexNetwork a) where
  put (ComplexNetwork ms) = put . V.toList $ ms
  get = (ComplexNetwork . V.fromList) `fmap` get

--instance (Complex a, Variate a) => Variate (Complex a) where
--	uniform g = uniform g :+ uniform g
--	uniformR (x,y) g = uniformR (x,y) g :+ uniformR (x,y) g

-- | The type of an activation function, mostly used for clarity in signatures
type ActivationFunction a = a -> a

-- | The type of an activation function's derivative, mostly used for clarity in signatures
type ActivationFunctionDerivative a = a -> a

--

randComplex :: (Variate a, Num a) => IO(Complex a)
randComplex = withSystemRandom . asGenST $ \gen -> do
  r <- uniformR (-1,1) gen
  i <- uniformR (-1,1) gen
  return (r :+ i)

randComplexList :: (Variate a, Num a) => Int -> IO([Complex a])
randComplexList 0 = return []
randComplexList n = do
  c <- randComplex
  lst <- randComplexList (n-1)
  return (c : lst)

randComplexMatrix :: (Variate a, Storable a, Num a) => (Int,Int) -> IO(Matrix (Complex a))
randComplexMatrix (rows,cols) = do
  lst <- randComplexList (rows*cols)
  return (reshape cols $ Numeric.LinearAlgebra.fromList lst)

-- | The following creates a neural network with 'n' inputs and if 'l' is [n1, n2, ...]
--   the net will have n1 neurons on the first layer, n2 neurons on the second, and so on
--   ending with k neurons on the output layer, with random weight matrices as a courtesy of
-- 'System.Random.MWC.uniformR'.
-- > createComplexNetwork n l k
createComplexNetwork :: (Variate a, Storable a, Num a) => Int -> [Int] -> Int -> IO (ComplexNetwork (Complex a))
createComplexNetwork nInputs hiddens nOutputs =
  fmap ComplexNetwork $ go dimensions V.empty
  where
        go [] !ms         = return ms
        go ((!n,!m):ds) ms = do
          !mat <- randComplexMatrix (n,m)
          go ds (ms `V.snoc` mat)
        dimensions      = zip (hiddens ++ [nOutputs]) $
                              (nInputs : hiddens)
{-# INLINE createComplexNetwork #-}

-- | Creates a neural network with exactly the weight matrices given as input here.
--   We don't check that the numbers of rows/columns are compatible, etc.
fromWeightMatrices :: Storable a => V.Vector (Matrix (Complex a)) -> ComplexNetwork (Complex a)
fromWeightMatrices ws = ComplexNetwork ws
{-# INLINE fromWeightMatrices #-}

-- The `join [input, 1]' trick  below is a courtesy of Alberto Ruiz
-- <http://dis.um.es/~alberto/>. Per his words:
--
-- "The idea is that the constant input in the first layer can be automatically transferred to the following layers
-- by the learning algorithm (by setting the weights of a neuron to 1,0,0,0,...). This allows for a simpler
-- implementation and in my experiments those networks are able to easily solve non linearly separable problems."

--without the trick (of constant input)

output :: (Floating (Vector a), Product a, Storable a, Num (Vector a)) => ComplexNetwork a -> ActivationFunction a -> Vector a -> Vector a
output (ComplexNetwork{..}) act input = V.foldl' f input matrices
  where f !inp m = mapVector act $ m <> inp
{-# INLINE output #-}

-- | Computes and keeps the output of all the layers of the neural network with the given activation function
outputs :: (Floating (Vector a), Product a, Storable a, Num (Vector a)) => ComplexNetwork a -> ActivationFunction a -> Vector a -> V.Vector (Vector a)
outputs (ComplexNetwork{..}) act input = V.scanl f input matrices
  where f !inp m = mapVector act $ m <> inp
{-# INLINE outputs #-}

conju (x:+y) = x:+(-y)

deltas :: (Floating b, Floating (Vector a), Floating a, Product a, Storable a, Num (Vector a), Container Vector a, a ~ Complex b) => ComplexNetwork a -> ActivationFunctionDerivative a -> V.Vector (Vector a) -> Vector a -> V.Vector (Matrix a)
deltas (ComplexNetwork{..}) act' os expected = V.zipWith outer (V.tail ds) (V.init (V.map (mapVector conju) os))
  where !dl = (V.last os - expected) * (deriv $ mapVector conju (V.last os)) -- = (out_last - target)*(f'(netin_last)) = dev_last * deriv(out_last) = partial_last
        !ds = V.scanr f dl (V.zip os matrices) -- generates dev_i
        f (!o, m) !del = deriv o * ((ctrans m) <> del) -- dev_k = deriv(out_k)*dev_k = deriv(out_k)*(WEIGHTS <> partial_k+1)
        deriv = mapVector act'
{-# INLINE deltas #-}

updateComplexNetwork :: (Floating b, Floating (Vector a), Floating a, Product a, Storable a, Num (Vector a), Container Vector a, a ~ Complex b) => a -> ActivationFunction a -> ActivationFunctionDerivative a -> ComplexNetwork a -> Sample a -> ComplexNetwork a
updateComplexNetwork alpha act act' n@(ComplexNetwork{..}) (input, expectedOutput) = ComplexNetwork $ V.zipWith (+) matrices corr
    where !xs = outputs n act input
          !ds = deltas n act' xs expectedOutput
          !corr = V.map (scale (-alpha)) ds
{-# INLINE updateComplexNetwork #-}

-- | Input vector and expected output vector
type Sample a = (Vector a, Vector a)

-- | List of 'Sample's
type Samples a = [Sample a]

-- | Handy operator to describe your learning set, avoiding unnecessary parentheses. It's just a synonym for '(,)'.
--   Generally you'll load your learning set from a file, a database or something like that, but it can be nice for
--   quickly playing with hnn or for simple problems where you manually specify your learning set.
--   That is, instead of writing:
--
-- > samples :: Samples Double
-- > samples = [ (fromList [0, 0], fromList [0])
-- >           , (fromList [0, 1], fromList [1])
-- >           , (fromList [1, 0], fromList [1])
-- >           , (fromList [1, 1], fromList [0])
-- >           ]
--
--   You can write:
--
-- > samples :: Samples Double
-- > samples = [ fromList [0, 0] --> fromList [0]
-- >           , fromList [0, 1] --> fromList [1]
-- >           , fromList [1, 0] --> fromList [1]
-- >           , fromList [1, 1] --> fromList [0]
-- >           ]
(-->) :: Vector a -> Vector a -> Sample a
(-->) = (,)

backpropOnce :: (Floating (Vector a), Floating b, Floating a, Product a, Num (Vector a), Container Vector a, a ~ Complex b) => a -> ActivationFunction a -> ActivationFunctionDerivative a -> ComplexNetwork a -> Samples a -> ComplexNetwork a
backpropOnce rate act act' n samples = foldl' (updateComplexNetwork rate act act') n samples
{-# INLINE backpropOnce #-}

-- | Generic training function.
--
-- The first argument is a predicate that will tell the backpropagation algorithm when to stop.
-- The first argument to the predicate is the epoch, i.e the number of times the backprop has been
-- executed on the samples. The second argument is /the current network/, and the third is the list of samples.
-- You can thus combine these arguments to create your own criterion.
--
-- For example, if you want to stop learning either when the network's quadratic error on the samples,
-- using the tanh function, is below 0.01, or after 1000 epochs, whichever comes first, you could
-- use the following predicate:
--
-- > pred epochs net samples = if epochs == 1000 then True else quadError tanh net samples < 0.01
--
-- You could even use 'Debug.Trace.trace' to print the error, to see how the error evolves while it's learning,
-- or redirect this to a file from your shell in order to generate a pretty graphics and what not.
--
-- The second argument (after the predicate) is the learning rate. Then come the activation function you want,
-- its derivative, the initial neural network, and your training set.
-- Note that we provide 'trainNTimes' and 'trainUntilErrorBelow' for common use cases.
trainUntil :: (Floating (Vector a), Floating b, Floating a, Product a, Num (Vector a), Container Vector a, a ~ Complex b) => (Int -> ComplexNetwork a -> Samples a -> Bool) -> a -> ActivationFunction a -> ActivationFunctionDerivative a -> ComplexNetwork a -> Samples a -> ComplexNetwork a
trainUntil pr learningRate act act' net samples = go net 0
  where go n !k | pr k n samples = n
                | otherwise      = case backpropOnce learningRate act act' n samples of
                                    n' -> go n' (k+1)
{-# INLINE trainUntil #-}

-- | Trains the neural network with backpropagation the number of times specified by the 'Int' argument,
-- using the given learning rate (second argument).
trainNTimes :: (Floating (Vector a), Floating b, Floating a, Product a, Num (Vector a), Container Vector a, a ~ Complex b) => Int -> a -> ActivationFunction a -> ActivationFunctionDerivative a -> ComplexNetwork a -> Samples a -> ComplexNetwork a
trainNTimes n = trainUntil (\k _ _ -> k > n)
{-# INLINE trainNTimes #-}

-- | Quadratic error on the given training set using the given activation function. Useful to create
-- your own predicates for 'trainUntil'.
quadError :: (Floating (Vector a), Floating b, Floating a, Num (Vector a), Num (RealOf a), Product a, a ~ Complex b) => ActivationFunction a -> ComplexNetwork a -> Samples a -> RealOf a
quadError act net samples = foldl' (\err (inp, out) -> err + (norm2 $ output net act inp - out)) 0 samples
{-# INLINE quadError #-}

-- | Trains the neural network until the quadratic error ('quadError') comes below the given value (first argument),
-- using the given learning rate (second argument).
--
-- /Note/: this can loop pretty much forever when you're using a bad architecture for the problem, or unappropriate activation
-- functions.
trainUntilErrorBelow :: (Floating (Vector a), Floating b, Floating a, Product a, Num (Vector a), Ord a, Container Vector a, Num (RealOf a), a ~ RealOf a, Show a, a ~ Complex b) => a -> a -> ActivationFunction a -> ActivationFunctionDerivative a -> ComplexNetwork a -> Samples a -> ComplexNetwork a
trainUntilErrorBelow x rate act = trainUntil (\_ n s -> quadError act n s < x) rate act
{-# INLINE trainUntilErrorBelow #-}

-- | The sigmoid function:  1 / (1 + exp (-x))
sigmoid :: Floating a => a -> a
sigmoid !x = 1 / (1 + exp (-x))
{-# INLINE sigmoid #-}

-- | Derivative of the sigmoid function: sigmoid x * (1 - sigmoid x)
sigmoid' :: Floating a => a -> a
sigmoid' !x = case sigmoid x of
  s -> s * (1 - s)
{-# INLINE sigmoid' #-}

complexSigmoid :: Floating a => Complex a -> Complex a
complexSigmoid (x :+ y) = (sigmoid x) :+ (sigmoid y)

complexSigmoid' :: Floating a => Complex a -> Complex a
complexSigmoid' (x :+ y) = (sigmoid' x) :+ (sigmoid' y)

-- | Loading a neural network from a file (uses zlib compression on top of serialization using the binary package).
--   Will throw an exception if the file isn't there.
loadComplexNetwork :: (Storable a, Element a, Binary a) => FilePath -> IO (ComplexNetwork a)
loadComplexNetwork fp = decode . decompress <$> B.readFile fp
{-# INLINE loadComplexNetwork #-}

-- | Saving a neural network to a file (uses zlib compression on top of serialization using the binary package).
saveComplexNetwork :: (Storable a, Element a, Binary a) => FilePath -> ComplexNetwork a -> IO ()
saveComplexNetwork fp = B.writeFile fp . compress . encode
{-# INLINE saveComplexNetwork #-}

max_stars_repo_name: EditResearch/HsANN

max_issues_repo_name: EditResearch/HsANN

max_forks_repo_name: EditResearch/HsANN


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Packed.ST
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Mutable packed matrices.
--
module Numeric.LinearAlgebra.Packed.ST (
    -- * Mutable packed matrices
    STPacked,
    IOPacked,
    create,

    -- * Read-only packed matrices
    RPacked(..),

    -- * Conversions between mutable and immutable packed matrices
    freeze,

    -- * Creating new packed matrices
    new_,

    -- * Copying matrices
    newCopy,

    -- * Vector views of packed matrices
    withVectorM,

    -- * Packed matrix views of vectors
    withFromVector,
    withFromVectorM,

    ) where

import Numeric.LinearAlgebra.Packed.Base

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas