text:
{-# LANGUAGE AllowAmbiguousTypes    #-}
{-# LANGUAGE ConstraintKinds        #-}
{-# LANGUAGE DataKinds              #-}
{-# LANGUAGE DeriveDataTypeable     #-}
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE GADTs                  #-}
{-# LANGUAGE InstanceSigs           #-}
{-# LANGUAGE KindSignatures         #-}
{-# LANGUAGE LambdaCase             #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE PatternSynonyms        #-}
{-# LANGUAGE RankNTypes             #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE TupleSections          #-}
{-# LANGUAGE TypeApplications       #-}
{-# LANGUAGE TypeFamilies           #-}
{-# LANGUAGE TypeFamilyDependencies #-}
{-# LANGUAGE TypeInType             #-}
{-# LANGUAGE TypeOperators          #-}
{-# LANGUAGE UndecidableInstances   #-}
{-# LANGUAGE ViewPatterns           #-}

module Backprop.Learn.Model.Combinator (
    Chain(..)
  , (~++)
  , chainParamLength
  , chainStateLength
  , LearnFunc(..), learnFunc
  , LMap(..), RMap(..)
  , Dimap, pattern DM, _dmPre, _dmPost, _dmLearn
  , (.~)
  , nilLF, onlyLF
  , (:.~)(..)
  , (:&&&)(..)
  , KAutoencoder, kAutoencoder, kaeEncoder, kaeDecoder
  , Feedback(..), feedbackId
  , FeedbackTrace(..), feedbackTraceId
  ) where

import           Backprop.Learn.Model.Class
import           Backprop.Learn.Model.Function
import           Control.Applicative
import           Control.Category
import           Control.Monad
import           Control.Monad.Primitive
import           Control.Monad.Trans.State
import           Data.Bifunctor
import           Data.Kind
import           Data.Type.Equality
import           Data.Type.Length
import           Data.Type.Mayb                        as Mayb
import           Data.Type.NonEmpty
import           Data.Type.Tuple hiding                (T2(..), T3(..))
import           Data.Typeable
import           GHC.TypeNats
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop
import           Prelude hiding                        ((.), id)
import           Type.Class.Known
import           Type.Class.Witness
import           Type.Family.List                      as List
import qualified Data.Type.Tuple                       as T
import qualified Data.Vector.Sized                     as SV
import qualified System.Random.MWC                     as MWC

-- | Chain components linearly, retaining the ability to deconstruct at
-- a later time.
data Chain :: [Type] -> Type -> Type -> Type where
    CNil  :: Chain '[] a a
    (:~>) :: (Learn a b l, KnownMayb (LParamMaybe l), KnownMayb (LStateMaybe l))
          => l
          -> Chain ls        b c
          -> Chain (l ': ls) a c
  deriving (Typeable)
infixr 5 :~>

instance ( ListC (Backprop List.<$> LParams ls)
         , ListC (Backprop List.<$> LStates ls)
         )
      => Learn a b (Chain ls a b) where
    type LParamMaybe (Chain ls a b) = NETup Mayb.<$> ToNonEmpty (LParams ls)
    type LStateMaybe (Chain ls a b) = NETup Mayb.<$> ToNonEmpty (LStates ls)

    runLearn      = runChainLearn
    runLearnStoch = runChainLearnStoch

runChainLearn
    :: (Reifies s W, ListC (Backprop List.<$> LParams ls), ListC (Backprop List.<$> LStates ls))
    => Chain ls a b
    -> Mayb (BVar s) (NETup Mayb.<$> ToNonEmpty (LParams ls))
    -> BVar s a
    -> Mayb (BVar s) (NETup Mayb.<$> ToNonEmpty (LStates ls))
    -> (BVar s b, Mayb (BVar s) (NETup Mayb.<$> ToNonEmpty (LStates ls)))
runChainLearn = \case
  CNil -> \_ x _ -> (x, N_)
  (l :: l) :~> ls ->
    let lenPs = chainParamLength ls
        lenSs = chainStateLength ls
    in case knownMayb @(LParamMaybe l) of
      N_ -> \ps x -> case knownMayb @(LStateMaybe l) of
        N_ -> \ss -> flip (runChainLearn ls ps) ss
                   . runLearnStateless l N_
                   $ x
        J_ _ -> case lenSs of
          LZ -> \case
            J_ (isoVar (tOnly . netT) (NETT . onlyT)->s) ->
              let (y, J_ s') = runLearn      l  N_ x (J_ s)
                  (z, _    ) = runChainLearn ls ps y N_
              in  (z, J_ $ isoVar (NETT . onlyT) (tOnly . netT) s')
          LS _ -> \case
            J_ ss -> lenSs //
              let (y, J_ s' ) = runLearn      l  N_ x (J_ (ss ^^. netHead))
                  ssTail      = isoVar NETT netT $ ss ^^. netTail
                  (z, J_ ss') = runChainLearn ls ps y (J_ ssTail)
              in  (z, J_ $ isoVar2 NET unNet s' $ isoVar netT NETT ss')
      J_ _ -> case lenPs of
        LZ -> \case
          J_ (isoVar (tOnly . netT) (NETT . onlyT)->p) -> \x -> case knownMayb @(LStateMaybe l) of
            N_ -> \ss -> flip (runChainLearn ls N_) ss
                       . runLearnStateless l (J_ p)
                       $ x
            J_ _ -> case lenSs of
              LZ -> \case
                J_ (isoVar (tOnly . netT) (NETT . onlyT)->s) ->
                  let (y, J_ s') = runLearn      l  (J_ p)  x (J_ s)
                      (z, _    ) = runChainLearn ls N_      y N_
                  in  (z, J_ $ isoVar (NETT . onlyT) (tOnly . netT) s')
              LS _ -> \case
                J_ ss -> lenSs //
                  let (y, J_ s' ) = runLearn      l  (J_ p) x (J_ (ss ^^. netHead))
                      ssTail      = isoVar NETT netT $ ss ^^. netTail
                      (z, J_ ss') = runChainLearn ls N_     y (J_ ssTail)
                  in  (z, J_ $ isoVar2 NET unNet s' $ isoVar netT NETT ss')
        LS _ -> \case
          J_ ps -> \x -> lenPs //
            let psHead = ps ^^. netHead
                psTail = isoVar NETT netT $ ps ^^. netTail
            in  case knownMayb @(LStateMaybe l) of
                  N_ -> \ss -> flip (runChainLearn ls (J_ psTail)) ss
                             . runLearnStateless l (J_ psHead)
                             $ x
                  J_ _ -> case lenSs of
                    LZ -> \case
                      J_ (isoVar (tOnly . netT) (NETT . onlyT)->s) ->
                        let (y, J_ s') = runLearn      l  (J_ psHead) x (J_ s)
                            (z, _    ) = runChainLearn ls (J_ psTail) y N_
                        in  (z, J_ $ isoVar (NETT . onlyT) (tOnly . netT) s')
                    LS _ -> \case
                      J_ ss -> lenSs //
                        let (y, J_ s' ) = runLearn      l  (J_ psHead) x (J_ (ss ^^. netHead))
                            ssTail      = isoVar NETT netT $ ss ^^. netTail
                            (z, J_ ss') = runChainLearn ls (J_ psTail) y (J_ ssTail)
                        in  (z, J_ $ isoVar2 NET unNet s' $ isoVar netT NETT ss')

runChainLearnStoch
    :: ( Reifies s W
       , ListC (Backprop List.<$> LParams ls)
       , ListC (Backprop List.<$> LStates ls)
       , PrimMonad m
       )
    => Chain ls a b
    -> MWC.Gen (PrimState m)
    -> Mayb (BVar s) (NETup Mayb.<$> ToNonEmpty (LParams ls))
    -> BVar s a
    -> Mayb (BVar s) (NETup Mayb.<$> ToNonEmpty (LStates ls))
    -> m (BVar s b, Mayb (BVar s) (NETup Mayb.<$> ToNonEmpty (LStates ls)))
runChainLearnStoch = \case
  CNil -> \_ _ x _ -> pure (x, N_)
  (l :: l) :~> ls -> \g ->
    let lenPs = chainParamLength ls
        lenSs = chainStateLength ls
    in case knownMayb @(LParamMaybe l) of
      N_ -> \ps x -> case knownMayb @(LStateMaybe l) of
        N_ -> \ss -> flip (runChainLearnStoch ls g ps) ss
                 <=< runLearnStochStateless l g N_
                   $ x
        J_ _ -> case lenSs of
          LZ -> \case
            J_ (isoVar (tOnly . netT) (NETT . onlyT)->s) -> do
              (y, s') <- second fromJ_
                     <$> runLearnStoch      l  g N_ x (J_ s)
              (z, _ ) <- runChainLearnStoch ls g ps y N_
              pure (z, J_ $ isoVar (NETT . onlyT) (tOnly . netT) s')
          LS _ -> \case
            J_ ss -> lenSs // do
              (y, s' ) <- second fromJ_
                      <$> runLearnStoch      l  g N_ x (J_ (ss ^^. netHead))
              let ssTail = isoVar NETT netT $ ss ^^. netTail
              (z, ss') <- second fromJ_
                      <$> runChainLearnStoch ls g ps y (J_ ssTail)
              pure  (z, J_ $ isoVar2 NET unNet s' $ isoVar netT NETT ss')
      J_ _ -> case lenPs of
        LZ -> \case
          J_ (isoVar (tOnly . netT) (NETT . onlyT)->p) -> \x -> case knownMayb @(LStateMaybe l) of
            N_ -> \ss -> flip (runChainLearnStoch ls g N_) ss
                     <=< runLearnStochStateless l g (J_ p)
                       $ x
            J_ _ -> case lenSs of
              LZ -> \case
                J_ (isoVar (tOnly . netT) (NETT . onlyT)->s) -> do
                  (y, s') <- second fromJ_
                         <$> runLearnStoch      l  g (J_ p)  x (J_ s)
                  (z, _ ) <- runChainLearnStoch ls g N_      y N_
                  pure (z, J_ $ isoVar (NETT . onlyT) (tOnly . netT) s')
              LS _ -> \case
                J_ ss -> lenSs // do
                  (y, s' ) <- second fromJ_
                          <$> runLearnStoch      l  g (J_ p) x (J_ (ss ^^. netHead))
                  let ssTail = isoVar NETT netT $ ss ^^. netTail
                  (z, ss') <- second fromJ_
                          <$> runChainLearnStoch ls g N_     y (J_ ssTail)
                  pure (z, J_ $ isoVar2 NET unNet s' $ isoVar netT NETT ss')
        LS _ -> \case
          J_ ps -> \x -> lenPs //
            let psHead = ps ^^. netHead
                psTail = isoVar NETT netT $ ps ^^. netTail
            in  case knownMayb @(LStateMaybe l) of
                  N_ -> \ss -> flip (runChainLearnStoch ls g (J_ psTail)) ss
                           <=< runLearnStochStateless l g (J_ psHead)
                             $ x
                  J_ _ -> case lenSs of
                    LZ -> \case
                      J_ (isoVar (tOnly . netT) (NETT . onlyT)->s) -> do
                        (y, s') <- second fromJ_
                               <$> runLearnStoch      l  g (J_ psHead) x (J_ s)
                        (z, _ ) <- runChainLearnStoch ls g (J_ psTail) y N_
                        pure (z, J_ $ isoVar (NETT . onlyT) (tOnly . netT) s')
                    LS _ -> \case
                      J_ ss -> lenSs // do
                        (y, s' ) <- second fromJ_
                                <$> runLearnStoch      l  g (J_ psHead) x (J_ (ss ^^. netHead))
                        let ssTail = isoVar NETT netT $ ss ^^. netTail
                        (z, ss') <- second fromJ_
                                <$> runChainLearnStoch ls g (J_ psTail) y (J_ ssTail)
                        pure (z, J_ $ isoVar2 NET unNet s' $ isoVar netT NETT ss')

-- | Appending 'Chain'
(~++)
    :: forall ls ms a b c. ()
    => Chain ls a b
    -> Chain ms b c
    -> Chain (ls ++ ms) a c
(~++) = \case
    CNil     -> id
    (l :: l) :~> (ls :: Chain ls' a' b) ->
      case assocMaybAppend @(LParamMaybe l) @(LParams ls') @(LParams ms) known of
        Refl -> case assocMaybAppend @(LStateMaybe l) @(LStates ls') @(LStates ls) known of
          Refl -> \ms -> (l :~> (ls ~++ ms))
            \\ appendLength (chainParamLength ls) (chainParamLength ms)
            \\ appendLength (chainStateLength ls) (chainStateLength ms)

chainParamLength
    :: Chain ls a b
    -> Length (LParams ls)
chainParamLength = \case
    CNil -> LZ
    (_ :: l) :~> ls -> case knownMayb @(LParamMaybe l) of
      N_   -> chainParamLength ls
      J_ _ -> LS $ chainParamLength ls

chainStateLength
    :: Chain ls a b
    -> Length (LStates ls)
chainStateLength = \case
    CNil -> LZ
    (_ :: l) :~> ls -> case knownMayb @(LStateMaybe l) of
      N_   -> chainStateLength ls
      J_ _ -> LS $ chainStateLength ls

appendLength
    :: forall as bs. ()
    => Length as
    -> Length bs
    -> Length (as ++ bs)
appendLength LZ     = id
appendLength (LS l) = LS . appendLength l

assocMaybAppend
    :: forall a bs cs. ()
    => Mayb P a
    -> (MaybeToList a ++ (bs ++ cs)) :~: ((MaybeToList a ++ bs) ++ cs)
assocMaybAppend = \case
    N_   -> Refl
    J_ _ -> Refl

-- | Data type representing trainable models.
--
-- Useful for performant composition, but you lose the ability to decompose
-- parts.
data LearnFunc :: Maybe Type -> Maybe Type -> Type -> Type -> Type where
    LF :: { _lfRunLearn
               :: forall q. Reifies q W
               => Mayb (BVar q) p
               -> BVar q a
               -> Mayb (BVar q) s
               -> (BVar q b, Mayb (BVar q) s)
          , _lfRunLearnStoch
               :: forall m q. (PrimMonad m, Reifies q W)
               => MWC.Gen (PrimState m)
               -> Mayb (BVar q) p
               -> BVar q a
               -> Mayb (BVar q) s
               -> m (BVar q b, Mayb (BVar q) s)
          }
       -> LearnFunc p s a b
  deriving (Typeable)

learnFunc
    :: Learn a b l
    => l
    -> LearnFunc (LParamMaybe l) (LStateMaybe l) a b
learnFunc l = LF { _lfRunLearn      = runLearn l
                 , _lfRunLearnStoch = runLearnStoch l
                 }

instance Learn a b (LearnFunc p s a b) where
    type LParamMaybe (LearnFunc p s a b) = p
    type LStateMaybe (LearnFunc p s a b) = s

    runLearn      = _lfRunLearn
    runLearnStoch = _lfRunLearnStoch

instance Category (LearnFunc p s) where
    id = LF { _lfRunLearn      = \_ -> (,)
            , _lfRunLearnStoch = \_ _ x -> pure . (x,)
            }
    f . g = LF { _lfRunLearn  = \p x s0 ->
                    let (y, s1) = _lfRunLearn g p x s0
                    in  _lfRunLearn f p y s1
               , _lfRunLearnStoch = \gen p x s0 -> do
                    (y, s1) <- _lfRunLearnStoch g gen p x s0
                    _lfRunLearnStoch f gen p y s1
               }

-- | Compose two 'LearnFunc' on lists.
(.~)
    :: forall ps qs ss ts a b c.
     ( ListC (Backprop List.<$> ps)
     , ListC (Backprop List.<$> qs)
     , ListC (Backprop List.<$> ss)
     , ListC (Backprop List.<$> ts)
     , ListC (Backprop List.<$> (ss ++ ts))
     , Known Length ps
     , Known Length ss
     , Known Length ts
     )
    => LearnFunc ('Just (T ps        )) ('Just (T ss         )) b c
    -> LearnFunc ('Just (T qs        )) ('Just (T ts         )) a b
    -> LearnFunc ('Just (T (ps ++ qs))) ('Just (T (ss ++ ts ))) a c
f .~ g = LF { _lfRunLearn  = \(J_ psqs) x (J_ ssts) -> appendLength @ss @ts known known //
                let (y, J_ ts) = _lfRunLearn g (J_ (psqs ^^. tDrop @ps @qs known))
                                               x
                                               (J_ (ssts ^^. tDrop @ss @ts known))
                    (z, J_ ss) = _lfRunLearn f (J_ (psqs ^^. tTake @ps @qs known))
                                               y
                                               (J_ (ssts ^^. tTake @ss @ts known))
                in  (z, J_ $ isoVar2 (tAppend @ss @ts) (tSplit @ss @ts known)
                                     ss ts
                    )
            , _lfRunLearnStoch = \gen (J_ psqs) x (J_ ssts) -> appendLength @ss @ts known known // do
                (y, ts) <- second fromJ_
                       <$> _lfRunLearnStoch g gen (J_ (psqs ^^. tDrop @ps @qs known))
                                                  x
                                                  (J_ (ssts ^^. tDrop @ss @ts known))
                (z, ss) <- second fromJ_
                       <$> _lfRunLearnStoch f gen (J_ (psqs ^^. tTake @ps @qs known))
                                                   y
                                                   (J_ (ssts ^^. tTake @ss @ts known))
                pure  (z, J_ $ isoVar2 (tAppend @ss @ts) (tSplit @ss @ts known)
                                       ss ts
                      )
            }

-- | Identity of '.~'
nilLF :: LearnFunc ('Just (T '[])) ('Just (T '[])) a a
nilLF = id

-- | 'LearnFunc' with singleton lists; meant to be used with '.~'
onlyLF
    :: forall p s a b. (KnownMayb p, MaybeC Backprop p, KnownMayb s, MaybeC Backprop s)
    => LearnFunc p s a b
    -> LearnFunc ('Just (T (MaybeToList p))) ('Just (T (MaybeToList s))) a b
onlyLF f = LF
    { _lfRunLearn = \(J_ ps) x ssM@(J_ ss) -> case knownMayb @p of
        N_ -> case knownMayb @s of
          N_ -> (second . const) ssM
              $ _lfRunLearn f N_ x N_
          J_ _ -> second (J_ . isoVar onlyT tOnly . fromJ_)
                $ _lfRunLearn f N_ x (J_ (isoVar tOnly onlyT ss))
        J_ _ ->
          let p = isoVar tOnly onlyT ps
          in  case knownMayb @s of
                N_ -> (second . const) ssM
                    $ _lfRunLearn f (J_ p) x N_
                J_ _ -> second (J_ . isoVar onlyT tOnly . fromJ_)
                      $ _lfRunLearn f (J_ p) x (J_ (isoVar tOnly onlyT ss))
    , _lfRunLearnStoch = \g (J_ ps) x ssM@(J_ ss) -> case knownMayb @p of
        N_ -> case knownMayb @s of
          N_ -> (fmap . second . const) ssM
              $ _lfRunLearnStoch f g N_ x N_
          J_ _ -> (fmap . second) (J_ . isoVar onlyT tOnly . fromJ_)
                $ _lfRunLearnStoch f g N_ x (J_ (isoVar tOnly onlyT ss))
        J_ _ ->
          let p = isoVar tOnly onlyT ps
          in  case knownMayb @s of
                N_ -> (fmap . second . const) ssM
                    $ _lfRunLearnStoch f g (J_ p) x N_
                J_ _ -> (fmap . second) (J_ . isoVar onlyT tOnly . fromJ_)
                      $ _lfRunLearnStoch f g (J_ p) x (J_ (isoVar tOnly onlyT ss))
    }

-- | Compose two layers sequentially
--
-- Note that this composes in the opposite order of '.' and ':.:', for
-- consistency with the rest of the library.
--
-- The basic autoencoder is simply @l ':.~' m@, where @'Learn' a b l@ and
-- @'Learn' b a m@.  However, for sparse autoencoder, look at
-- 'Autoencoder'.
data (:.~) :: Type -> Type -> Type where
    (:.~) :: l -> m -> l :.~ m

infixr 5 :.~

instance ( Learn a b l
         , Learn b c m
         , KnownMayb (LParamMaybe l)
         , KnownMayb (LParamMaybe m)
         , KnownMayb (LStateMaybe l)
         , KnownMayb (LStateMaybe m)
         , MaybeC Backprop (LParamMaybe l)
         , MaybeC Backprop (LParamMaybe m)
         , MaybeC Backprop (LStateMaybe l)
         , MaybeC Backprop (LStateMaybe m)
         )
      => Learn a c (l :.~ m) where
    type LParamMaybe (l :.~ m) = TupMaybe (LParamMaybe l) (LParamMaybe m)
    type LStateMaybe (l :.~ m) = TupMaybe (LStateMaybe l) (LStateMaybe m)

    runLearn (l :.~ m) pq x st = (z, tupMaybe T2B s' t')
      where
        (p, q) = splitTupMaybe @_ @(LParamMaybe l) @(LParamMaybe m)
                  (\(T2B v u) -> (v, u))
                  pq
        (s, t) = splitTupMaybe @_ @(LStateMaybe l) @(LStateMaybe m)
                  (\(T2B v u) -> (v, u))
                  st
        (y, s') = runLearn l p x s
        (z, t') = runLearn m q y t

    runLearnStoch (l :.~ m) g pq x st = do
        (y, s') <- runLearnStoch l g p x s
        (z, t') <- runLearnStoch m g q y t
        pure (z, tupMaybe T2B s' t')
      where
        (p, q) = splitTupMaybe @_ @(LParamMaybe l) @(LParamMaybe m)
                  (\(T2B v u) -> (v, u))
                  pq
        (s, t) = splitTupMaybe @_ @(LStateMaybe l) @(LStateMaybe m)
                  (\(T2B v u) -> (v, u))
                  st

-- | Pre-compose a pure parameterless function to a model.
--
-- An @'LMap' b a@ takes a model from @b@ and turns it into a model from
-- @a@.
data LMap :: Type -> Type -> Type -> Type where
    LM :: { _lmPre   :: forall s. Reifies s W => BVar s a -> BVar s b
          , _lmLearn :: l
          }
       -> LMap b a l

-- | Post-compose a pure parameterless function to a model.
--
-- An @'Rmap' b c@ takes a model returning @b@ and turns it into
-- a model returning @c@.
data RMap :: Type -> Type -> Type -> Type where
    RM :: { _rmPost  :: forall s. Reifies s W => BVar s b -> BVar s c
          , _rmLearn :: l
          }
       -> RMap b c l

-- | Pre- and post-compose pure parameterless functions to a model.
--
-- A @'Dimap' b c a d@ takes a model from @b@ to @c@ and turns it into
-- a model from @a@ to @d@.
--
-- @
-- instance 'Learn' b c => Learn a d ('Dimap' b c a d l) where
--     type 'LParamMaybe' (Dimap b c a d l) = LParamMaybe l
--     type 'LStateMaybe' (Dimap b c a d l) = LStateMaybe l
-- @
type Dimap b c a d l = RMap c d (LMap b a l)

-- | Constructor for 'Dimap'
pattern DM
    :: (forall s. Reifies s W => BVar s a -> BVar s b)
    -> (forall s. Reifies s W => BVar s c -> BVar s d)
    -> l
    -> Dimap b c a d l
pattern DM { _dmPre, _dmPost, _dmLearn } = RM _dmPost (LM _dmPre _dmLearn)

instance Learn b c l => Learn a c (LMap b a l) where
    type LParamMaybe (LMap b a l) = LParamMaybe l
    type LStateMaybe (LMap b a l) = LStateMaybe l

    runLearn (LM f l) p x = runLearn l p (f x)
    runLearnStoch (LM f l) g p x = runLearnStoch l g p (f x)

instance Learn a b l => Learn a c (RMap b c l) where
    type LParamMaybe (RMap b c l) = LParamMaybe l
    type LStateMaybe (RMap b c l) = LStateMaybe l

    runLearn (RM f l) p x = first f . runLearn l p x
    runLearnStoch (RM f l) g p x = (fmap . first) f . runLearnStoch l g p x

-- | Take a model and turn it into a model that runs its output into itself
-- several times, returning the final result.  Parameterized by the number
-- of repeats, and the function to process the output to become the next
-- input.
--
-- I don't know why anyone would ever want this.
--
-- See 'FeedbackTrace' if you want to observe all of the intermediate
-- outputs.
data Feedback :: Type -> Type -> Type -> Type where
    FB :: { _fbTimes :: Int
          , _fbFeed  :: forall s. Reifies s W => BVar s b -> BVar s a
          , _fbLearn :: l
          }
       -> Feedback a b l
  deriving (Typeable)

-- | Construct a 'Feedback' from an endofunction (a function that returns
-- a value fo the same type as its input) by simply providing the output
-- directly as the next intput.
feedbackId :: Int -> l -> Feedback a a l
feedbackId n = FB n id

instance Learn a b l => Learn a b (Feedback a b l) where
    type LParamMaybe (Feedback a b l) = LParamMaybe l
    type LStateMaybe (Feedback a b l) = LStateMaybe l

    runLearn (FB n f l) p = runState
                          . foldr (>=>) go (replicate (n - 1) (fmap f . go))
      where
        go = state . runLearn l p

    runLearnStoch (FB n f l) g p = runStateT
                                 . foldr (>=>) go (replicate (n - 1) (fmap f . go))
      where
        go = StateT . runLearnStoch l g p

-- | Take a model and turn it into a model that runs its output into itself
-- several times, and returns all of the intermediate outputs.
-- Parameterized by the function to process the output to become the next
-- input.
--
-- See 'Feedback' if you only need the final result.
--
-- Compare also to 'Unroll'.
data FeedbackTrace :: Nat -> Type -> Type -> Type -> Type where
    FBT :: { _fbtFeed  :: forall s. Reifies s W => BVar s b -> BVar s a
           , _fbtLearn :: l
           }
        -> FeedbackTrace n a b l
  deriving (Typeable)

-- | Construct a 'FeedbackTrace' from an endofunction (a function that
-- returns a value fo the same type as its input) by simply providing the
-- output directly as the next intput.
feedbackTraceId :: l -> FeedbackTrace n a a l
feedbackTraceId = FBT id

instance (Learn a b l, KnownNat n, Backprop b)
      => Learn a (ABP (SV.Vector n) b) (FeedbackTrace n a b l) where
    type LParamMaybe (FeedbackTrace n a b l) = LParamMaybe l
    type LStateMaybe (FeedbackTrace n a b l) = LStateMaybe l

    runLearn (FBT f l) p x0 =
          second snd
        . runState (collectVar . ABP <$> SV.replicateM (state go))
        . (x0,)
      where
        go (x, s) = (y, (f y, s'))
          where
            (y, s') = runLearn l p x s

    runLearnStoch (FBT f l) g p x0 =
          (fmap . second) snd
        . runStateT (collectVar . ABP <$> SV.replicateM (StateT go))
        . (x0,)
      where
        go (x, s) = do
          (y, s') <- runLearnStoch l g p x s
          pure (y, (f y, s'))

-- | "Fork"/"Fan out" two models from the same input.  Useful for
data (:&&&) :: Type -> Type -> Type where
    (:&&&) :: l -> m -> l :&&& m

infixr 3 :&&&

instance ( Learn a b l
         , Learn a c m
         , KnownMayb (LParamMaybe l)
         , KnownMayb (LParamMaybe m)
         , KnownMayb (LStateMaybe l)
         , KnownMayb (LStateMaybe m)
         , MaybeC Backprop (LParamMaybe l)
         , MaybeC Backprop (LParamMaybe m)
         , MaybeC Backprop (LStateMaybe l)
         , MaybeC Backprop (LStateMaybe m)
         , Backprop b
         , Backprop c
         )
      => Learn a (T.T2 b c) (l :&&& m) where
    type LParamMaybe (l :&&& m) = TupMaybe (LParamMaybe l) (LParamMaybe m)
    type LStateMaybe (l :&&& m) = TupMaybe (LStateMaybe l) (LStateMaybe m)

    runLearn (l :&&& m) pq x st = ( T2B y z
                                  , tupMaybe T2B s' t'
                                  )
      where
        (p, q) = splitTupMaybe @_ @(LParamMaybe l) @(LParamMaybe m)
                  (\(T2B v u) -> (v, u))
                  pq
        (s, t) = splitTupMaybe @_ @(LStateMaybe l) @(LStateMaybe m)
                  (\(T2B v u) -> (v, u))
                  st
        (y, s') = runLearn l p x s
        (z, t') = runLearn m q x t

    runLearnStoch (l :&&& m) g pq x st = do
        (y, s') <- runLearnStoch l g p x s
        (z, t') <- runLearnStoch m g q x t
        pure ( T2B y z
             , tupMaybe T2B s' t'
             )
      where
        (p, q) = splitTupMaybe @_ @(LParamMaybe l) @(LParamMaybe m)
                  (\(T2B v u) -> (v, u))
                  pq
        (s, t) = splitTupMaybe @_ @(LStateMaybe l) @(LStateMaybe m)
                  (\(T2B v u) -> (v, u))
                  st

-- | K-sparse autoencoder.  A normal autoencoder is simply @l ':.~' m@;
-- however a k-sparse autoencoder attempts to ensure that the encoding has
-- about @k@ active components for every input.
--
-- <http://www.ericlwilkinson.com/blog/2014/11/19/deep-learning-sparse-autoencoders>
--
-- @
-- instance ('Learn' a ('R' n) l, Learn (R n) a m) => Learn a a ('KAutoencoder' n l m) where
--     type 'LParamMaybe' ('KAutoencoder' n l m) = 'TupMaybe' (LParamMaybe l) (LParamMaybe m)
--     type 'LStateMaybe' ('KAutoencoder' n l m) = 'TupMaybe' (LStateMaybe l) (LStateMaybe m)
-- @
--
-- To "encode" after training this, get the encoder with 'kaeEncoder'.
type KAutoencoder n l m = RMap (R n) (R n) l :.~ m

-- | Construct a 'KAutoencoder'.
--
-- Note that this only has a 'Learn' instance of @l@ outputs @'R' n@ and
-- @m@ takes @'R' n@.  Also, for this to be an actual autoencoder, the
-- input of @l@ has to be the same as the output of @m@.
kAutoencoder
    :: KnownNat n
    => l
    -> m
    -> Int
    -> KAutoencoder n l m
kAutoencoder l m k = RM (kSparse k) l
                 :.~ m

kaeEncoder
    :: KAutoencoder n l m
    -> l
kaeEncoder (RM _ l :.~ _) = l

kaeDecoder
    :: KAutoencoder n l m
    -> m
kaeDecoder (_ :.~ m) = m

-- TODO: KL-divergence autoencoders, a la
-- <http://ufldl.stanford.edu/tutorial/unsupervised/Autoencoders/>.

-- -- | Simple /sparse/ autoencoder that outputs the average activation of
-- -- a vector encoding as well as the result.
-- --
-- -- The traditional autoencoder is simply @l ':.~' m@.  However, the enforce
-- -- sparsity, you have to also be able to observe the average activation of
-- -- your encoding for your loss function (typically 'klDivergence' for
-- -- Kullback-Leibler divergence)
-- --
-- -- See <http://ufldl.stanford.edu/tutorial/unsupervised/Autoencoders/>.
-- --
-- -- @
-- -- instance ('Learn' a ('R' n) l, Learn (R n) a l, 1 '<=' n) => Learn a ('T2' 'Double' a) where
-- --     type 'LParamMaybe' ('Autoencoder' n l m) = 'TupMaybe' (LParamMaybe l) (LParamMaybe m)
-- --     type 'LStateMaybe' ('Autoencoder' n l m) = 'TupMaybe' (LStateMaybe l) (LStateMaybe m)
-- -- @
-- --
-- -- To /use/ the autoencoder after training it, just pattern match on 'AE'
-- -- or use '_aeEncode' (or '_aeDecode')
-- type Autoencoder n l m = l :.~ (FixedFunc (R n) Double :&&& m)
--
-- -- | Construct an 'Autoencoder' by giving the encoder and decoder.
-- pattern AE
--     :: (Learn a (R n) l, Learn (R n) a m, 1 <= n, KnownNat n)
--     => l                    -- ^ '_aeEncode'
--     -> m                    -- ^ '_aeDecode'
--     -> Autoencoder n l m
-- pattern AE { _aeEncode, _aeDecode } <- _aeEncode :.~ (_ :&&& _aeDecode)
--   where
--     AE e d = e :.~ (FF mean :&&& d)

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
-- |
-- Module      : BattleHack.Window
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created September 18 2015

-- TODO | -
--        -

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module BattleHack.Window where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Control.Monad      (forM_, when)
import Control.Concurrent (threadDelay, forkIO)
import Text.Printf
import Data.IORef
import Data.Complex
import Data.StateVar
import qualified Data.Map  as M

import qualified Graphics.Rendering.Cairo as Cairo
import           Graphics.UI.Gtk          as Gtk
import           Graphics.UI.Gtk          (AttrOp(..), on)

-- Internal module imports
import           BattleHack.Types
import qualified BattleHack.Render  as Render
import qualified BattleHack.Events  as Events
import qualified BattleHack.Audio   as Audio

--------------------------------------------------------------------------------------------------------------------------------------------
-- Functions
--------------------------------------------------------------------------------------------------------------------------------------------
-- |
-- TODO: Return canvas and frame too (?)
-- TODO: Allow for configurations
create :: Vector -> IO (Window, DrawingArea)
create (winx:+winy) = do
  initGUI
  window <- windowNew
  frame  <- frameNew

  set window [windowTitle := "Chordially"]

  canvas <- drawingAreaNew
  containerAdd frame canvas
  set window [ containerChild := frame ]
  windowSetDefaultSize window (round winx) (round winy)
  windowSetIconFromFile window "assets/images/gclef.png"

  widgetAddEvents canvas [PointerMotionMask] -- MouseButton1Mask
  widgetShowAll window

  return (window, canvas)

-- | Register event listeners
-- TODO: Rename (?)
bindevents :: Window -> DrawingArea -> IORef AppState -> IO ()
bindevents window canvas stateref = do
  window `on` deleteEvent        $ Events.ondelete      stateref
  window `on` motionNotifyEvent  $ Events.onmousemotion stateref

  window `on` buttonPressEvent   $ Events.onmousedown stateref
  window `on` buttonReleaseEvent $ Events.onmouseup   stateref

  window `on` keyPressEvent      $ Events.onkeydown stateref
  window `on` keyReleaseEvent    $ Events.onkeyup   stateref

  canvas `on` scrollEvent        $ Events.onwheelscrool stateref
  canvas `on` draw               $ Events.ondraw        stateref

  return ()

ext: hs


text:
{-# LANGUAGE FlexibleContexts, NamedFieldPuns #-}

module School.Unit.MultiNoulli
( multiNoulli ) where

import Conduit ((.|), ConduitM, MonadResource, ZipSource(..),
                getZipSource, mapC, mapM_C, mapMC)
import Control.Monad.Except (MonadError, throwError)
import Numeric.LinearAlgebra (Container, Element, Matrix, Vector, assoc, cols, fromRows,
                              rows, takeColumns, toColumns, toList, toLists)
import School.FileIO.FileHeader (FileHeader)
import School.FileIO.Source (source)
import School.Types.Error (Error)
import School.Types.LiftResult (LiftResult)
import School.Types.Slinky (Slinky(..), slinkyAppend, slinkySingleton)
import School.Unit.CostParams (CostParams(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitForward (ForwardStack)
import School.Unit.UnitGradient (UnitGradient(..))
import School.Unit.CostFunction (SetupCost, CostFunction(..), defSetupCost)

toTarget :: (Element a, RealFrac a)
         => Vector a
         -> CostParams
toTarget =
  BatchClassTarget . map round . toList

compute :: (Element a, Fractional a, Num a)
        => UnitActivation a
        -> Slinky CostParams
        -> Either Error a
compute (BatchActivation input)
        (SNode (BatchClassTarget target) _) =
  let factor = (-1) / (fromIntegral . rows $ input)
  in Right
   . (*factor)
   . sum
   . zipWith (flip (!!)) target
   . toLists
   $ input
compute _ _ = Left "MultiNoulli expects batch activation and batch class target"

deriv :: (Container Vector a, Fractional a, Num a)
      => UnitActivation a
      -> Slinky CostParams
      -> Either Error (UnitGradient a)
deriv (BatchActivation input)
      (SNode (BatchClassTarget target) _) =
  let factor = (-1) / (fromIntegral . rows $ input)
      c = cols input
  in Right
   . BatchGradient
   . fromRows
   . map (\idx -> assoc c 0 [(idx, factor)])
   $ target
deriv _ _ = Left "MultiNoulli expects batch activation and batch class target"

prepare :: (Element a, RealFrac a, MonadError Error m)
        => Maybe FilePath
        -> ConduitM (ForwardStack a)
                    (ForwardStack a)
                    m
                    ()
prepare Nothing = mapMC $ \(activations, cParams) ->
  case activations of
    [BatchActivation input] -> do
      let c = cols input
      let activation = takeColumns (c - 1) input
      let newParams = toTarget . last . toColumns $ input
      return ([BatchActivation activation], slinkyAppend newParams cParams)
    _ -> throwError "MultiNoulli setup expects single batch activation"
prepare _ = mapC id

setupError :: (MonadError Error m) => ConduitM () (ForwardStack a) m ()
setupError =
  mapM_C . const $ throwError "Need both filepath and file header"

appendTarget :: CostParams -> Matrix a -> ForwardStack a
appendTarget costParams =
  flip (,) (slinkySingleton costParams) . pure . BatchActivation

setup :: (LiftResult m, MonadError Error m, MonadResource m)
      => Maybe FilePath
      -> Maybe FileHeader
      -> SetupCost a m
setup Nothing Nothing byteSource = defSetupCost byteSource
setup Nothing (Just _) _ = setupError
setup (Just _) Nothing _ = setupError
setup (Just path) (Just header) matrixSource =
  getZipSource $ targetZipSource <*> ZipSource matrixSource
  where targetSource = source header path
                    .| mapC BatchClassTarget
                    .| mapMC (return . appendTarget)
        targetZipSource = ZipSource targetSource

multiNoulli :: ( Container Vector a, LiftResult m, RealFrac a
               , MonadError Error m, MonadResource m )
            => Maybe FilePath
            -> Maybe FileHeader
            -> CostFunction a m
multiNoulli path header = CostFunction { computeCost = compute
                                       , derivCost = deriv
                                       , prepareCost = prepare path
                                       , setupCost = setup path header
                                       }

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
-- Principal component analysis

import Numeric.LinearAlgebra
import System.Directory(doesFileExist)
import System.Process(system)
import Control.Monad(when)

type Vec = Vector Double
type Mat = Matrix Double

-- Vector with the mean value of the columns of a matrix
mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a

-- covariance matrix of a list of observations stored as rows
cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
    where xc = x - asRow (mean x)

-- creates the compression and decompression functions from the desired number of components
pca :: Int -> Mat -> (Vec -> Vec , Vec -> Vec)
pca n dataSet = (encode,decode)
  where
    encode x = vp <> (x - m)
    decode x = x <> vp + m
    m = mean dataSet
    c = cov dataSet
    (_,v) = eigSH' c
    vp = takeRows n (trans v)

norm = pnorm PNorm2

main = do
    ok <- doesFileExist ("mnist.txt")
    when (not ok)  $ do
        putStrLn "\nTrying to download test datafile..."
        system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
        system("gunzip mnist.txt.gz")
        return ()
    m <- loadMatrix "mnist.txt" -- fromFile "mnist.txt" (5000,785)
    let xs = takeColumns (cols m -1) m -- the last column is the digit type (class label)
    let x = toRows xs !! 4  -- an arbitrary test Vec
    let (pe,pd) = pca 10 xs
    let y = pe x
    print y  -- compressed version
    print $ norm (x - pd y) / norm x --reconstruction quality

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# LANGUAGE AllowAmbiguousTypes                      #-}
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE DefaultSignatures                        #-}
{-# LANGUAGE DeriveDataTypeable                       #-}
{-# LANGUAGE DeriveFoldable                           #-}
{-# LANGUAGE DeriveFunctor                            #-}
{-# LANGUAGE DeriveGeneric                            #-}
{-# LANGUAGE DeriveTraversable                        #-}
{-# LANGUAGE DerivingVia                              #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE FunctionalDependencies                   #-}
{-# LANGUAGE GeneralizedNewtypeDeriving               #-}
{-# LANGUAGE MultiParamTypeClasses                    #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE StandaloneDeriving                       #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeFamilies                             #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# OPTIONS_GHC -Wno-redundant-constraints            #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

-- |
-- Module      : Numeric.Opto.Update
-- Copyright   : (c) Justin Le 2019
-- License     : BSD3
--
-- Maintainer  : justin@jle.im
-- Stability   : experimental
-- Portability : non-portable
--
-- A unified interface for values in vector spaces that can be added and
-- scaled (purely and also in-place), and measured.
module Numeric.Opto.Update (
    Linear(..), sumLinear, gAdd, gZeroL, gScale
  , Metric(..), gDot, gNorm_inf, gNorm_0, gNorm_1, gNorm_2, gQuadrance
  , LinearInPlace(..), sumLinearInPlace
  , linearWit
  ) where

import           Control.DeepSeq
import           Control.Monad.Primitive
import           Data.Coerce
import           Data.Complex
import           Data.Data
import           Data.Finite
import           Data.Foldable
import           Data.Function
import           Data.Maybe
import           Data.Mutable
import           Data.Semigroup
import           Data.Vinyl hiding                   ((:~:))
import           GHC.Generics                        (Generic)
import           GHC.TypeLits
import           Generics.OneLiner
import           Numeric.Opto.Ref
import           Unsafe.Coerce
import qualified Data.Vector.Generic                 as VG
import qualified Data.Vector.Generic.Mutable.Sized   as SVGM
import qualified Data.Vector.Generic.Sized           as SVG
import qualified Numeric.LinearAlgebra               as UH
import qualified Numeric.LinearAlgebra.Static        as H
import qualified Numeric.LinearAlgebra.Static.Vector as H

-- | If @a@ is an instance of @'Linear' c@, you can /add/ together values
-- of @a@, and /scale/ them using @c@s.
--
-- For example, if you have a vector of doubles, you can add them together
-- component-wise, and scale them by multiplying every item by the scalar.
--
-- Mathematically, this means that @a@ forms something like a module or
-- vector space over @c@, where @c@ can be any 'Num' instance.
class Num c => Linear c a | a -> c where
    -- | Add together @a@s.  Should be associative.
    --
    -- @
    -- x .+. (y .+. z) == (x .+. y) .+. z
    -- @
    --
    -- If @a@ is an instance of 'Num', this can be just 'Prelude.+'.
    (.+.) :: a -> a -> a

    -- | The "zero" @a@, meant to form an identity with '.+.'.
    --
    -- @
    -- x .+. zeroL == x
    -- zeroL .+. y == y
    -- @
    --
    -- If @a@ is an instance of 'Num', this can be just 0.
    zeroL :: a

    -- | Scale an @a@ by a factor @c@.  Should distribute over '.+.'.
    --
    -- @
    -- a .* (x .+. y) == (a .* x) .+. (a .* y)
    -- a .* (b .* c)  == (a * b) .* c
    -- @
    (.*)  :: c -> a -> a

    infixl 6 .+.
    infixl 7 .*

    default (.+.) :: (ADTRecord a, Constraints a (Linear c)) => a -> a -> a
    (.+.) = gAdd @c

    default zeroL :: (ADTRecord a, Constraints a (Linear c)) => a
    zeroL = gZeroL @c

    default (.*) :: (ADTRecord a, Constraints a (Linear c)) => c -> a -> a
    (.*)  = gScale

-- | Sum over a 'Foldable' container of @'Linear' c a@
sumLinear :: (Linear c a, Foldable t) => t a -> a
sumLinear = foldl' (.+.) zeroL

-- | An implementation of '.+.' that works for records where every field is
-- an instance of @'Linear' c@ (that is, every field is additive and can be
-- scaled by the same @c@).
gAdd :: forall c a. (ADTRecord a, Constraints a (Linear c)) => a -> a -> a
gAdd = binaryOp @(Linear c) (.+.)

-- | An implementation of 'zeroL' that works for records where every field
-- is an instance of @'Linear' c@ (that is, every field is additive and can
-- be scaled by the same @c@).
gZeroL :: forall c a. (ADTRecord a, Constraints a (Linear c)) => a
gZeroL = nullaryOp @(Linear c) zeroL

-- | An implementation of '.*' that works for records where every field
-- is an instance of @'Linear' c@ (that is, every field is additive and can
-- be scaled by the same @c@).
gScale :: forall c a. (ADTRecord a, Constraints a (Linear c)) => c -> a -> a
gScale c = unaryOp @(Linear c) (c .*)

-- | Class for values supporting an inner product and various norms.
class Linear c a => Metric c a where
    infixl 7 <.>
    -- | Sum of component-wise product
    (<.>)    :: a -> a -> c
    -- | Maximum absolute component.  Is undefined if no components exist.
    norm_inf :: a -> c
    -- | Number of non-zero components
    norm_0   :: a -> c
    -- | Sum of absolute components
    norm_1   :: a -> c
    -- | Square root of sum of squared components
    norm_2    :: a -> c
    -- | Sum of squared components
    quadrance :: a -> c

    default (<.>) :: (ADT a, Constraints a (Metric c)) => a -> a -> c
    (<.>) = gDot
    default norm_inf :: (ADT a, Constraints a (Metric c), Ord c) => a -> c
    norm_inf = gNorm_inf
    default norm_0 :: (ADT a, Constraints a (Metric c)) => a -> c
    norm_0 = gNorm_0
    default norm_1 :: (ADT a, Constraints a (Metric c)) => a -> c
    norm_1 = gNorm_1
    default norm_2 :: Floating c => a -> c
    norm_2 = sqrt . quadrance
    default quadrance :: (ADT a, Constraints a (Metric c)) => a -> c
    quadrance = gQuadrance

-- | An implementation of 'gDot' that works for records where every field
-- is an instance of @'Metric' c@.
gDot :: forall c a. (ADT a, Constraints a (Metric c), Num c) => a -> a -> c
gDot x = getSum . mzipWith @(Metric c) (\x' -> Sum . (x' <.>)) x

-- | An implementation of 'norm_inf' that works for records where every
-- field is an instance of @'Metric' c@.
gNorm_inf :: forall c a. (ADT a, Constraints a (Metric c), Ord c) => a -> c
gNorm_inf = getMax
          . fromMaybe (error "norm_inf: Divergent infinity norm")
          . getOption
          . gfoldMap @(Metric c) (Option . Just . Max . abs . norm_inf)

-- | An implementation of 'norm_0' that works for records where every field
-- is an instance of @'Metric' c@.
gNorm_0 :: forall c a. (ADT a, Constraints a (Metric c), Num c) => a -> c
gNorm_0 = getSum . gfoldMap @(Metric c) (Sum . norm_0)

-- | An implementation of 'norm_1' that works for records where every field
-- is an instance of @'Metric' c@.
gNorm_1 :: forall c a. (ADT a, Constraints a (Metric c), Num c) => a -> c
gNorm_1 = getSum . gfoldMap @(Metric c) (Sum . norm_1)

-- | An implementation of 'norm_2' that works for records where every field
-- is an instance of @'Metric' c@.
gNorm_2 :: forall c a. (ADT a, Constraints a (Metric c), Floating c) => a -> c
gNorm_2 = sqrt . gQuadrance

-- | An implementation of 'quadrance' that works for records where every
-- field is an instance of @'Metric' c@.
gQuadrance :: forall c a. (ADT a, Constraints a (Metric c), Num c) => a -> c
gQuadrance = getSum . gfoldMap @(Metric c) (Sum . quadrance)

-- | Instaces of 'Linear' that support certain in-place mutations.
-- Inspired by the BLAS Level 1 API.  A @'LinearInPlace' m v c a@ means
-- that @v@ is a mutable reference to an @a@ that can be updated as an
-- action in monad @m@.
class (Mutable s a, Linear c a) => LinearInPlace s c a where
    -- | Add a value in-place.
    (.+.=) :: (PrimMonad m, PrimState m ~ s) => Ref s a -> a -> m ()

    -- | Scale a value in-place.
    (.*=)  :: (PrimMonad m, PrimState m ~ s) => Ref s a -> c -> m ()

    -- | Add a scaled value in-place.
    (.*+=) :: (PrimMonad m, PrimState m ~ s) => Ref s a -> (c, a) -> m ()

    r .+.= x      = modifyRef' r (.+. x)
    r  .*= c      = modifyRef' r (c .*)
    r .*+= (c, x) = modifyRef' r ((c .* x) .+.)

    infix 4 .+.=
    infix 4 .*=
    infix 4 .*+=

-- | Given some starting reference @v@, add every item in a foldable
-- container into that reference in-place.
sumLinearInPlace :: (LinearInPlace s c a, Foldable t, PrimMonad m, PrimState m ~ s) => Ref s a -> t a -> m ()
sumLinearInPlace v = mapM_ (v .+.=)

-- | Newtype wrapper that gives "simple" 'Linear', 'Metric', and
-- 'LinearInPlace' instances for instances of 'Num'.
--
-- This can be used with the /-XDerivingVia/ extension:
--
-- @
-- deriving via (LinearNum Double) instance Linear Double Double
-- deriving via (LinearNum Double) instance Metric Double Double
-- instance Ref s Double v => LinearInPlace s v Double Double
-- @
newtype LinearNum a = LinearNum { getLinearNum :: a }
  deriving ( Show, Eq, Ord
           , Functor, Foldable, Traversable
           , Enum, Bounded
           , Num, Fractional, Floating, Real, Integral, RealFrac, RealFloat
           , Generic, Typeable, Data
           )
instance NFData a => NFData (LinearNum a)

instance Num a => Linear a (LinearNum a) where
    (.+.) = (+)
    zeroL = 0
    (.*)  = coerce ((*) :: a -> a -> a)
instance Num a => Metric a (LinearNum a) where
    (<.>)     = coerce ((*) :: a -> a -> a)
    norm_inf  = coerce (abs :: a -> a)
    norm_0    = coerce (abs . signum :: a -> a)
    norm_1    = coerce (abs :: a -> a)
    norm_2    = coerce (abs :: a -> a)
    quadrance = coerce ((^ (2 :: Int)) :: a -> a)

deriving via (LinearNum Int)         instance Linear Int Int
deriving via (LinearNum Integer)     instance Linear Integer Integer
deriving via (LinearNum Rational)    instance Linear Rational Rational
deriving via (LinearNum Float)       instance Linear Float Float
deriving via (LinearNum Double)      instance Linear Double Double
deriving via (LinearNum (Complex a)) instance RealFloat a => Linear (Complex a) (Complex a)

deriving via (LinearNum Int)         instance Metric Int Int
deriving via (LinearNum Integer)     instance Metric Integer Integer
deriving via (LinearNum Rational)    instance Metric Rational Rational
deriving via (LinearNum Float)       instance Metric Float Float
deriving via (LinearNum Double)      instance Metric Double Double
deriving via (LinearNum (Complex a)) instance RealFloat a => Metric (Complex a) (Complex a)

instance Mutable s Int                        => LinearInPlace s Int Int
instance Mutable s Integer                    => LinearInPlace s Integer Integer
instance Mutable s Rational                   => LinearInPlace s Rational Rational
instance Mutable s Float                      => LinearInPlace s Float Float
instance Mutable s Double                     => LinearInPlace s Double Double
instance (Mutable s (Complex a), RealFloat a) => LinearInPlace s (Complex a) (Complex a)

instance (Num a, VG.Vector v a, KnownNat n) => Linear a (SVG.Vector v n a) where
    (.+.)    = (+)
    zeroL    = 0
    c .* xs  = SVG.map (c *) xs

instance (Floating a, Ord a, VG.Vector v a, KnownNat n) => Metric a (SVG.Vector v n a) where
    xs <.> ys = SVG.sum (xs * ys)
    norm_inf  = SVG.foldl' (\x y -> max (abs x) y) 0
    norm_0    = fromIntegral . SVG.length
    norm_1    = SVG.sum . abs
    quadrance = SVG.sum . (^ (2 :: Int))

instance (Num a, mv ~ VG.Mutable v, VG.Vector v a, KnownNat n)
      => LinearInPlace s a (SVG.Vector v n a) where
    r .+.= xs = flip SVG.imapM_ xs $ \i x ->
      SVGM.modify r (+ x) i
    r .*= c = forM_ finites $ \i ->
      SVGM.modify r (c *) i
    r .*+= (c, xs) = flip SVG.imapM_ xs $ \i x ->
      SVGM.modify r (+ (c * x)) i

instance KnownNat n => Linear Double (H.R n) where
    (.+.)   = (+)
    zeroL   = 0
    c .* xs = H.konst c * xs
instance KnownNat n => Metric Double (H.R n) where
    (<.>)     = (H.<.>)
    norm_inf  = H.norm_Inf
    norm_0    = H.norm_0
    norm_1    = H.norm_1
    norm_2    = H.norm_2
    quadrance = (**2) . H.norm_2
instance KnownNat n => LinearInPlace s Double (H.R n) where
    MR v .+.= x = v .+.= H.rVec x
    MR v  .*= c = v  .*= c
    MR v .*+= (c, x) = v .*+= (c, H.rVec x)

instance (KnownNat n, KnownNat m) => Linear Double (H.L n m) where
    (.+.)   = (+)
    zeroL   = 0
    c .* xs = H.konst c * xs
instance (KnownNat n, KnownNat m) => Metric Double (H.L n m) where
    (<.>)     = (UH.<.>) `on` UH.flatten . H.extract
    norm_inf  = UH.maxElement . H.extract . abs
    norm_0    = sum . map norm_0 . H.toRows
    norm_1    = UH.sumElements . H.extract
    norm_2    = UH.norm_2 . UH.flatten . H.extract
    quadrance = (**2) . norm_2
instance (KnownNat n, KnownNat k) => LinearInPlace s Double (H.L n k) where
    ML v .+.= x = v .+.= H.lVec x
    ML v  .*= c = v  .*= c
    ML v .*+= (c, x) = v .*+= (c, H.lVec x)

instance (Linear c a, Linear c b) => Linear c (a, b) where
instance (Linear c a, Linear c b, Linear c d) => Linear c (a, b, d) where
instance (Linear c a, Linear c b, Linear c d, Linear c e) => Linear c (a, b, d, e) where
instance (Linear c a, Linear c b, Linear c d, Linear c e, Linear c f) => Linear c (a, b, d, e, f) where

instance (Metric c a, Metric c b, Ord c, Floating c) => Metric c (a, b)
instance (Metric c a, Metric c b, Metric c d, Ord c, Floating c) => Metric c (a, b, d)
instance (Metric c a, Metric c b, Metric c d, Metric c e, Ord c, Floating c) => Metric c (a, b, d, e)
instance (Metric c a, Metric c b, Metric c d, Metric c e, Metric c f, Ord c, Floating c) => Metric c (a, b, d, e, f)

instance (Mutable s (a, b), Linear c a, Linear c b) => LinearInPlace s c (a, b)
instance (Mutable s (a, b, d), Linear c a, Linear c b, Linear c d) => LinearInPlace s c (a, b, d)
instance (Mutable s (a, b, d, e), Linear c a, Linear c b, Linear c d, Linear c e) => LinearInPlace s c (a, b, d, e)
instance (Mutable s (a, b, d, e, f), Linear c a, Linear c b, Linear c d, Linear c e, Linear c f) => LinearInPlace s c (a, b, d, e, f)

instance Linear c (f a) => Linear c (Rec f '[a]) where
    x :& RNil .+. y :& RNil = (x .+. y) :& RNil
    zeroL = zeroL :& RNil
    c .* (x :& RNil) = (c .* x) :& RNil

instance (Linear c (f a), Linear c (Rec f (b ': bs))) => Linear c (Rec f (a ': b ': bs)) where
    x :& xs .+. y :& ys = (x .+. y) :& (xs .+. ys)
    zeroL = zeroL :& zeroL
    c .* (x :& xs) = (c .* x) :& (c .* xs)

instance Metric c (f a) => Metric c (Rec f '[a]) where
    (x :& RNil) <.> (y :& RNil) = x <.> y
    norm_inf (x :& RNil) = norm_inf x
    norm_0 (x :& RNil) = norm_0 x
    norm_1 (x :& RNil) = norm_1 x
    norm_2 (x :& RNil) = norm_2 x
    quadrance (x :& RNil) = quadrance x

instance (Floating c, Ord c, Metric c (f a), Metric c (Rec f (b ': bs))) => Metric c (Rec f (a ': b ': bs)) where
    (<.>) = undefined
    norm_inf (x :& xs) = norm_0 x `max` norm_0 xs
    norm_0 (x :& xs) = norm_0 x + norm_0 xs
    norm_1 (x :& xs) = norm_1 x + norm_1 xs
    norm_2 = sqrt . quadrance
    quadrance (x :& xs) = quadrance x + quadrance xs

instance LinearInPlace s c (f a) => LinearInPlace s c (Rec f '[a]) where
    (RecRef v :& RNil) .+.= (x :& RNil)    = v .+.= x
    (RecRef v :& RNil) .*=  c              = v .*=  c
    (RecRef v :& RNil) .*+= (c, x :& RNil) = v .*+= (c, x)

instance ( LinearInPlace s c (f a)
         , LinearInPlace s c (Rec f (b ': bs))
         )
        => LinearInPlace s c (Rec f (a ': b ': bs)) where
    (RecRef v :& vs) .+.= (x :& xs)    = do v .+.= x; vs .+.= xs
    (RecRef v :& vs) .*=  c            = do v .*=  c; vs .*=  c
    (RecRef v :& vs) .*+= (c, x :& xs) = do v .*+= (c, x); vs .*+= (c, xs)

-- | If @a@ and @b@ are both 'Linear' instances, then if @a@ is equal to
-- @b@, their scalars @c@ and @d@ must also be equal.  This is necessary
-- because GHC isn't happy with the functional dependency for some reason.
linearWit
    :: forall a c d. (Linear c a, Linear d a)
    => (c :~: d)
linearWit = unsafeCoerce Refl

max_stars_repo_name: mstksg/opto

max_issues_repo_name: mstksg/opto

max_forks_repo_name: mstksg/opto


text:
{-# LANGUAGE FlexibleContexts #-}

module Geometry where

import Graphics.Gloss as G
import Graphics.Gloss.Data.ViewPort as G
import Graphics.Gloss.Data.ViewState as G
import Data.Array.Accelerate as A
import Data.Array.Accelerate.Data.Complex as A
import qualified Prelude as P
import qualified GHC.Float as P
import qualified Data.Complex as C

type FractalPoint = Complex Double

-- TODO: phantom types

width :: Int
width = 400

height :: Int
height = 400

blWorldByWidthHeight :: (Int, Int) -> Point
blWorldByWidthHeight (w, h) = (-P.fromIntegral w P./ 2.0, -P.fromIntegral h P./ 2.0)

urWorldByWidthHeight :: (Int, Int) -> Point
urWorldByWidthHeight (w, h) = (P.fromIntegral w P./ 2.0, P.fromIntegral h P./ 2.0)

blWorldInit :: Point
blWorldInit = blWorldByWidthHeight (width, height)

urWorldInit :: Point
urWorldInit = urWorldByWidthHeight (width, height)

blComplexInit :: FractalPoint
blComplexInit = (-1) :+ (-1)

urComplexInit :: FractalPoint
urComplexInit = (1) :+ (1)

convertWorldToComplex :: Point -> FractalPoint -> Point -> FractalPoint
convertWorldToComplex w_ur c_ur (x, y) = xComplex C.:+ yComplex
                                where xRel = x P./ P.fst w_ur
                                      yRel = y P./ P.snd w_ur
                                      xComplex = P.float2Double xRel P.* C.realPart c_ur
                                      yComplex = P.float2Double yRel P.* C.imagPart c_ur

makeComplexGrid :: Int -> Int -> FractalPoint -> FractalPoint -> Acc (A.Matrix FractalPoint)
makeComplexGrid w h bl ur = use $ A.fromList (Z:.h :. w) $ do
    imag <- P.reverse $ make1DComplexGrid h (C.imagPart bl) (C.imagPart ur)
    real <- make1DComplexGrid w (C.realPart bl) (C.realPart ur)
    P.return (real :+ imag)
        where
            make1DComplexGrid n l r = (P.+l) . (P./ P.fromIntegral n) . (P.*(r P.- l)) . P.fromIntegral P.<$> [0..n - 1]

max_stars_repo_name: psurkov/Fractal-generator

max_issues_repo_name: psurkov/Fractal-generator

max_forks_repo_name: psurkov/Fractal-generator


text:
module Lib
  ( lab3
  ) where

import Data.Complex
import Data.List
import Generator
import Graphics.Gnuplot.Simple
import qualified Graphics.Gnuplot.Terminal.PNG as PNG

dft :: [Float] -> Int -> Complex Float
dft xs k =
  sum $ zipWith (\element index -> toComplex element * tc index) xs [0.0 ..]
  where
    tc n = toComplex (cos $ arg n) - toComplex (sin $ arg n) * imOne
    arg n = 2.0 * pi * fromIntegral k * n / fromIntegral (length xs)

fft :: [Float] -> Int -> Bool -> Complex Float
fft xs k lower
  | lower = evenPart + expPart * oddPart
  | otherwise = evenPart - expPart * oddPart
  where
    evenPart = dft (evens xs) k
    oddPart = dft (odds xs) k
    expPart = exp $ -2 * pi * fromIntegral k * imOne / fromIntegral (length xs)

evens :: [a] -> [a]
evens (x:xs) = x : odds xs
evens _ = []

odds :: [a] -> [a]
odds (_:xs) = evens xs
odds _ = []

toComplex :: Num a => a -> Complex a
toComplex x = x :+ 0

imOne :: Complex Float
imOne = 0.0 :+ 1.0

printTwoLists :: Show a => [a] -> [a] -> IO ()
printTwoLists xs ys =
  putStrLn . unlines . map (intercalate "      \t\t" . map show) . transpose $
  xs : [ys]

lab3 :: IO ()
lab3 = do
  let sHarmonics = 8
  let sAmplitude = 100
  let sFrequency = 1200
  let sPhase = 100
  let sTime = 60
  let signal =
        [ randomSignal sHarmonics t sAmplitude sFrequency sPhase
        | t <- [0 .. sTime - 1]
        ]
  let spectral_dft = [dft signal k | k <- [0 .. sTime - 1]]
  let spectral_fft =
        [fft signal k True | k <- [0 .. n1 - 1]] ++
        [fft signal k False | k <- [0 .. n2 - 1]]
        where
          n1 = (ceiling :: Float -> Int) $ fromIntegral (length signal) / 2
          n2 = (floor :: Float -> Int) $ fromIntegral (length signal) / 2
  putStrLn . unlines . map show $ signal
  printTwoLists spectral_dft spectral_fft
  plotLists
    [ Key Nothing
    , XLabel "Time"
    , YLabel "Value"
    , Title "Random signal"
    , terminal (PNG.cons "signal.png")
    ]
    [signal]
  plotLists
    [ Key Nothing
    , XLabel "Time"
    , YLabel "Value"
    , Title "DFT"
    , terminal (PNG.cons "dft.png")
    ]
    [map realPart spectral_dft, map imagPart spectral_dft]
  plotLists
    [ Key Nothing
    , XLabel "Time"
    , YLabel "Value"
    , Title "FFT"
    , terminal (PNG.cons "fft.png")
    ]
    [map realPart spectral_fft, map imagPart spectral_fft]

max_stars_repo_name: arjaz/Fourier-transformation

max_issues_repo_name: arjaz/Fourier-transformation

max_forks_repo_name: arjaz/Fourier-transformation


text:
{-# LANGUAGE BangPatterns , RankNTypes, GADTs, DataKinds #-}

{- | The 'Numerical.HBLAS.BLAS' module provides a fully general
yet type safe BLAS API.

When in doubt about the semantics of an operation,
consult your system's BLAS api documentation, or just read the documentation
for
<https://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm the Intel MKL BLAS distribution>

A few basic notes about how to invoke BLAS routines.

Many BLAS operations take one or more arguments of type 'Transpose'.
'Tranpose' has the following different constructors, which tell BLAS
routines what transformation to implicitly apply to an input matrix @mat@ with dimension @n x m@.

*  'NoTranspose' leaves the matrix @mat@ as is.

* 'Transpose' treats the @mat@ as being implicitly transposed, with dimension
    @m x n@. Entry @mat(i,j)@ being treated as actually being the entry
    @mat(j,i)@. For Real matrices this is also the matrix adjoint operation.
    ie @Tranpose(mat)(i,j)=mat(j,i)@

*  'ConjNoTranspose' will implicitly conjugate @mat@, which is a no op for Real ('Float' or 'Double') matrices, but for
'Complex Float' and 'Complex Double' matrices, a given matrix entry @mat(i,j)==x':+'y@
will be treated as actually being  @conjugate(mat)(i,j)=y':+'x@.

* 'ConjTranpose' will implicitly transpose and conjugate the input matrix.
ConjugateTranpose acts as matrix adjoint for both real and complex matrices.

The *gemm operations  work as follows (using 'sgemm' as an example):

* @'sgemm trLeft trRight alpha beta left right result'@, where @trLeft@ and @trRight@
are values of type 'Transpose' that respectively act on the matrices @left@ and @right@.

* the generalized matrix computation thusly formed can be viewed as being
@result = alpha * trLeft(left) * trRight(right) + beta * result@

the *gemv operations are akin to the *gemm operations, but with @right@ and @result@
being vectors rather than matrices.

the *trsv operations solve for @x@ in the equation @A x = y@ given @A@ and @y@.
The 'MatUpLo' argument determines if the matrix should be treated as upper or
lower triangular and 'MatDiag' determines if the triangular solver should treat
the diagonal of the matrix as being all 1's or not.  A general pattern of invocation
would be @'strsv' matuplo  tranposeMatA  matdiag  matrixA  xVector@.
A key detail to note is that the input vector is ALSO the result vector,
ie 'strsv' and friends updates the vector place.

-}

module Numerical.HBLAS.BLAS(
        GemvFun
        ,GemmFun
        ,TrsvFun

        ,dgemm
        ,sgemm
        ,cgemm
        ,zgemm

        ,sgemv
        ,dgemv
        ,cgemv
        ,zgemv

        ,sger
        ,dger

        ,strsv
        ,dtrsv
        ,ctrsv
        ,ztrsv
            ) where

import Numerical.HBLAS.UtilsFFI
import Numerical.HBLAS.BLAS.FFI
import Numerical.HBLAS.BLAS.Internal
import Control.Monad.Primitive
import Data.Complex

sgemm :: PrimMonad m=>  GemmFun Float  orient  (PrimState m) m
sgemm =  gemmAbstraction "sgemm"  cblas_sgemm_safe cblas_sgemm_unsafe (\x f -> f x )

dgemm :: PrimMonad m=>  GemmFun  Double orient  (PrimState m) m
dgemm = gemmAbstraction "dgemm"  cblas_dgemm_safe cblas_dgemm_unsafe  (\x f -> f x )

cgemm :: PrimMonad m=>  GemmFun (Complex Float) orient  (PrimState m) m
cgemm = gemmAbstraction "cgemm" cblas_cgemm_safe cblas_cgemm_unsafe  withRStorable_

zgemm :: PrimMonad m=>  GemmFun (Complex Double) orient  (PrimState m) m
zgemm = gemmAbstraction "zgemm"  cblas_zgemm_safe cblas_zgemm_unsafe withRStorable_

sgemv :: PrimMonad m => GemvFun Float orient (PrimState m) m
sgemv = gemvAbstraction "sgemv" cblas_sgemv_safe cblas_sgemv_unsafe (flip ($))

dgemv :: PrimMonad m => GemvFun Double orient (PrimState m) m
dgemv = gemvAbstraction "dgemv" cblas_dgemv_safe cblas_dgemv_unsafe (flip ($))

cgemv :: PrimMonad m => GemvFun (Complex Float) orient (PrimState m) m
cgemv = gemvAbstraction "cgemv" cblas_cgemv_safe cblas_cgemv_unsafe withRStorable_

zgemv :: PrimMonad m => GemvFun (Complex Double) orient (PrimState m) m
zgemv = gemvAbstraction "zgemv" cblas_zgemv_safe cblas_zgemv_unsafe withRStorable_
strsv :: PrimMonad m => TrsvFun Float orient (PrimState m) m
strsv = trsvAbstraction "strsv" cblas_strsv_safe cblas_strsv_unsafe

dtrsv :: PrimMonad m => TrsvFun Double orient (PrimState m) m
dtrsv = trsvAbstraction "dtrsv" cblas_dtrsv_safe cblas_dtrsv_unsafe

ctrsv :: PrimMonad m => TrsvFun (Complex Float) orient (PrimState m) m
ctrsv = trsvAbstraction "ctrsv" cblas_ctrsv_safe cblas_ctrsv_unsafe

ztrsv :: PrimMonad m => TrsvFun (Complex Double) orient (PrimState m) m
ztrsv = trsvAbstraction "ztrsv" cblas_ztrsv_safe cblas_ztrsv_unsafe

sger :: PrimMonad m => GerFun Float orient (PrimState m) m
sger = gerAbstraction "sger" cblas_sger_safe cblas_sger_unsafe

dger :: PrimMonad m => GerFun Double orient (PrimState m) m
dger = gerAbstraction "dger" cblas_dger_safe cblas_dger_unsafe

max_stars_repo_name: archblob/hblas

max_issues_repo_name: archblob/hblas

max_forks_repo_name: archblob/hblas


text:
-- ch03 example 3.7

import Statistics.Distribution
import Statistics.Distribution.Binomial

b :: BinomialDistribution
b = binomial 5 0.05

main = do
  putStrLn $ "If 5 samples is tested, the probabiliy of observing at least one defective: "
          ++ show (1 - probability b 0)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE NamedFieldPuns #-}
{-# LANGUAGE TypeFamilies #-}

module Math.Markov where

import Control.Monad.Random (getRandoms, Rand, RandomGen)
import Data.Counter (count)
import Data.Default (Default(..))
import Data.List (nub, sort)
import Data.List.Split (splitOn)
import qualified Data.Map as M
import Data.Maybe (fromMaybe)
import qualified Data.Set as S
import qualified Data.Vector as VB
import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA

import Misc.Utils (ngrams)
import Math.Dist (DiscreteDist(..), IndexMap, JointDiscreteDistribution(..), Prob, Simulable(..), bisect, getProb, normalizeVec, trainDiscrete, valsToIndices)
import Math.Multiarray ((^!^), (^:^), MultiArray(..), multiarray)

-- MARKOV MODEL --

data MarkovModel v a = Markov { pi0 :: DiscreteDist v a, trans :: LA.Matrix Prob }
    deriving (Eq, Show)

rowStochasticMatrix :: [V.Vector Double] -> LA.Matrix Prob
rowStochasticMatrix rows = LA.fromRows $ normalizeVec <$> rows

-- construct a MarkovModel from an initial distribution and a transition matrix (entry i, j is probability of transition from i to j, so that the matrix is row-stochastic)
markov :: DiscreteDist v a -> [V.Vector Prob] -> MarkovModel v a
markov dist rows
    | length rows /= n               = error $ "transition matrix must have " ++ show n ++ " rows"
    | any (/= n) $ V.length <$> rows = error $ "all rows of transition matrix must have length " ++ show n
    | otherwise                      = Markov {pi0 = dist, trans = mat}
    where
        n = length $ vals dist
        mat = rowStochasticMatrix rows

instance (Enum v, Num v, Ord v, Ord a, Show a) => JointDiscreteDistribution MarkovModel v a where
    -- vars are always [0, 1, 2, ...]
    getVars :: MarkovModel v a -> [v]
    getVars _ = [0..]
    getJointProbEvent :: MarkovModel v a -> [[a]] -> Prob
    getJointProbEvent dist valss = exp $ getJointLogProbEvent dist valss
    getJointLogProb :: MarkovModel v a -> [a] -> Double
    getJointLogProb _ []                       = 0.0
    getJointLogProb (Markov {pi0, trans}) vals = sum logProbs
        where
            indices = valsToIndices pi0 vals
            transitionIndices = zip indices (tail indices)
            logProbs = log <$> ((probs pi0 V.! (head indices)) : (LA.atIndex trans <$> transitionIndices))
    -- marginal distribution is no longer a simple Markov chain
    marginalizeIndices = undefined
    -- infinite list of marginal distributions
    marginals :: MarkovModel v a -> [DiscreteDist v a]
    marginals (Markov {pi0, trans}) = dists
        where
            pis = normalizeVec <$> iterate (LA.<# trans) (probs pi0)
            vals' = vals pi0
            cdfs = V.scanl' (+) 0.0 <$> pis
            dists = [Discrete {var = i, vals = vals', probs = pmf, cdf = cdf'} | (i, pmf, cdf') <- zip3 [0..] pis cdfs]
    -- conditional distribution is no longer a simple Markov chain
    conditionOn = undefined
    -- distribution is infinite so can't be represented as a DiscreteDist
    forgetJoint = undefined

-- computes conditional distribution on a single variable given the nearest preceding and succeeding variable (if existent)
-- input is t and a pair ((k, val), (l, val')), and we compute P(X_t | X_{t - k} = val, X_{t + l} = val')
-- each subpair can also be Nothing, indicating the absence of such a condition
markovConditionOn :: (Integral v, Ord a) => MarkovModel v a -> v -> (Maybe (Int, a), Maybe (Int, a)) -> DiscreteDist v a
markovConditionOn (Markov {pi0, trans}) t pair = Discrete {var = t, vals = vals', probs = pmf', cdf = cdf'}
    where
        stepRow :: V.Vector Prob -> Int -> V.Vector Prob
        stepRow vec k = (iterate (LA.<# trans) vec) !! k
        stepCol :: V.Vector Prob -> Int -> V.Vector Prob
        stepCol vec k = (iterate (trans LA.#>) vec) !! k
        pis = iterate (LA.<# trans) (probs pi0)
        vals' = vals pi0
        n = VB.length vals'
        indVec :: Int -> V.Vector Prob
        indVec i = V.generate n (\j -> if (j == i) then 1.0 else 0.0)
        vecMul = V.zipWith (*)
        t' = (fromInteger $ toInteger t) :: Int
        valIdx = M.fromList $ zip (VB.toList (vals pi0)) [0..]
        val2Index = (valIdx M.!)
        pmf = case pair of
            (Nothing, Nothing)              -> pis !! t'
            (Just (k, val), Nothing)        -> stepRow (indVec $ val2Index val) k
            (Nothing, Just (l, val'))       -> stepCol (indVec $ val2Index val') l `vecMul` (pis !! t')
            (Just (k, val), Just (l, val')) -> stepCol (indVec $ val2Index val') l `vecMul` stepRow (indVec $ val2Index val) k
        pmf' = normalizeVec pmf
        cdf' = V.scanl' (+) 0.0 pmf'

-- trains a Markov model from data
-- smoothing constant is how much to add to each entry of transition matrix, and to each entry of initial vector
trainMarkovModel :: (Ord a, Show a) => Maybe [a] -> Double -> [[a]] -> MarkovModel Integer a
trainMarkovModel vals smooth chains = markov pi0 (LA.toRows trans)
    where
        allVals = fromMaybe (sort $ nub $ concat chains) vals
        chains' = filter (not . null) chains
        pi0 = trainDiscrete 0 (Just allVals) smooth (head <$> chains)
        indices = valsToIndices pi0 <$> chains
        transitionPairs = concat [zip inds (tail inds) | inds <- indices]
        n = length allVals
        addSmooth (key, val) = (key, val + smooth)
        trans = LA.assoc (n, n) smooth (addSmooth <$> (M.toList $ count transitionPairs))

-- samples a Markov model (infinitely)
markovSample :: (RandomGen g) => MarkovModel v a -> Rand g [a]
markovSample (Markov {pi0, trans}) = do
    rs <- getRandoms
    let indices = bisect (head rs) cdf0 : [bisect r (transCdfs VB.! i) | (i, r) <- zip indices (tail rs)]
    return $ (vs VB.!) <$> indices
    where
        vs = vals pi0
        cdf0 = cdf pi0
        transCdfs = VB.fromList $ (V.scanl' (+) 0.0) <$> (LA.toRows trans)

instance (b ~ [a]) => Simulable (MarkovModel v) a b where
    sample :: (RandomGen g) => MarkovModel v a -> Rand g [a]
    sample = markovSample

-- given a distribution on integers steps, creates a MarkovModel representing the random walk on the i.i.d. steps
-- must be bounded within a finite range
boundedIntegerRandomWalk :: DiscreteDist v Int -> (Int, Int) -> MarkovModel Integer Int
boundedIntegerRandomWalk dist (low, high) = markov pi0 rows
    where
        n = high - low + 1
        pi0 = trainDiscrete 0 (Just [low..high]) 0.0 [0]
        rows = [V.generate n (\j -> getProb dist (j + low - i)) | i <- [low..high]]

-- gets backward transition probabilities P(X_n = i | X_{n+1} = j)
backwardTransitions :: MarkovModel v a -> [V.Vector Prob]
backwardTransitions (Markov {trans}) = normalizeVec <$> (LA.toRows $ LA.tr trans)

-- TERMINATING MARKOV MODEL --

-- type for values in a terminating Markov chain
data Token a = Token a | Stop
    deriving (Eq, Ord, Show)

isToken :: Token a -> Bool
isToken (Token _) = True
isToken _         = False

fromToken :: Token a -> a
fromToken (Token t) = t
fromToken _         = error "invalid token"

-- converts elements to Tokens and appends a final Stop element
-- this Nothing value can be used as a terminating value for Markov chain generation
terminate :: [a] -> [Token a]
terminate = (++ [Stop]) . map Token

trainTerminatingMarkovModel :: (Ord a, Show a) => Maybe [a] -> Double -> [[a]] -> MarkovModel Integer (Token a)
trainTerminatingMarkovModel vals smooth chains = trainMarkovModel tokens smooth (terminate <$> chains)
    where tokens = (Stop :) . map Token <$> vals

instance {-# OVERLAPPING #-} (b ~ [a]) => Simulable (MarkovModel v) (Token a) b where
    -- override behavior so sampling terminates upon reaching Stop
    sample = fmap (map fromToken . takeWhile isToken) . markovSample

-- N-GRAM MODELS --

-- stores n, alphabet, mapping from alphabet to indices, and transition probability matrices for each Markovity level
data NgramModel v a = Ngrams Int (VB.Vector a) (IndexMap a) [MultiArray Prob]
    deriving (Eq, Show)

instance (Enum v, Num v, Ord v, Ord a, Show a) => JointDiscreteDistribution NgramModel v a where
    getVars :: NgramModel v a -> [v]
    getVars _ = [0..]
    getJointProbEvent :: NgramModel v a -> [[a]] -> Prob
    getJointProbEvent dist valss = exp $ getJointLogProbEvent dist valss
    getJointLogProb :: NgramModel v a -> [a] -> Double
    getJointLogProb (Ngrams n _ valIdx mats) vals = sum logprobs
        where
            indices = (valIdx M.!) <$> vals
            idxGrams = [take i indices | i <- [1..(n - 1)]] ++ ngrams n indices
            matSeq = init mats ++ repeat (last mats)
            logprobs = log <$> zipWith (^!^) matSeq idxGrams
    marginalizeIndices = undefined
    marginals = undefined
    conditionOn = undefined
    forgetJoint = undefined

trainNgramModel :: (Ord a, Show a) => Int -> Maybe [a] -> Double -> [[a]] -> NgramModel Integer a
trainNgramModel n alphabet smooth chains = Ngrams n (VB.fromList alphabet') valIdx transArrs
    where
        alphabet' = fromMaybe (sort $ nub $ concat chains) alphabet
        m = length alphabet'
        valIdx = M.fromList $ zip alphabet' [0..]  -- mapping from alphabet values to indices
        idxChains = map (valIdx M.!) <$> chains
        getTransitions :: Int -> MultiArray Prob
        getTransitions k = multiarray (replicate k m) (concat rows)
            where
                pairs = concat $ (map (splitAt (k - 1)) . ngrams k) <$> idxChains
                cumProds = (m ^) <$> reverse [0..(k - 2)]
                transitionPairs = [(sum $ zipWith (*) cumProds fromInds, head toInds) | (fromInds, toInds) <- pairs]
                addSmooth (key, val) = (key, val + smooth)
                trans = LA.assoc (m ^ (k - 1), m) smooth (addSmooth <$> (M.toList $ count transitionPairs))
                rows = V.toList . normalizeVec <$> LA.toRows trans
        transArrs = getTransitions <$> [1..n]

-- samples an n-gram model (infinitely)
ngramSample :: (RandomGen g) => NgramModel v a -> Rand g [a]
ngramSample (Ngrams n alphabet _ transArrs) = do
    rs <- getRandoms
    let indices = [] : [rotateInds arr prev r | (prev, arr, r) <- zip3 indices arrs rs]
    return [alphabet VB.! (last inds) | inds <- tail indices]
    where
        getTail :: [Int] -> [Int]
        getTail inds = if (length inds < n) then inds else tail inds
        rotateInds :: MultiArray Prob -> [Int] -> Double -> [Int]
        rotateInds arr inds r = indTail ++ [nextInd]
            where
                indTail = getTail inds
                row = entries $ arr ^:^ (((Just . pure) <$> indTail) ++ [Nothing])
                cdf = V.scanl' (+) 0.0 row
                nextInd = bisect r cdf
        arrs = transArrs ++ (repeat $ last transArrs)

instance (b ~ [a]) => Simulable (NgramModel v) a b where
    sample :: (RandomGen g) => NgramModel v a -> Rand g [a]
    sample = ngramSample

-- condSlice :: (Ord a) => NgramModel v a -> [a] -> [(Int, a, Prob)]
-- condSlice (Ngrams n vals valIdx transArrs) vals' = zip3 [0..] (VB.toList vals) (V.toList pdf)
--     where
--         lastN = reverse $ take (n - 1) $ reverse vals'
--         inds = (valIdx M.!) <$> lastN
--         arr = transArrs !! length inds
--         pdf = entries $ arr ^:^ (((Just . pure) <$> inds) ++ [Nothing])

-- TERMINATING N-GRAM MODELS --

trainTerminatingNgramModel :: (Ord a, Show a) => Int -> Maybe [a] -> Double -> [[a]] -> NgramModel Integer (Token a)
trainTerminatingNgramModel n alphabet smooth chains = trainNgramModel n toks smooth (terminate <$> chains)
    where toks = (Stop :) . map Token <$> alphabet

instance {-# OVERLAPPING #-} (b ~ [a]) => Simulable (NgramModel v) (Token a) b where
    -- override behavior so sampling terminates upon reaching Stop
    sample = fmap (map fromToken . takeWhile isToken) . ngramSample

-- TEXT N-GRAM MODELS --

data CharNgramParams = CharNgramParams {
      charAlphabet :: Maybe String,   -- alphabet to use
      charSmooth :: Double,           -- smoothing constant
      charLineSep :: String,          -- line separator string
      charTransform :: Char -> Char,  -- function to transform each character
      charFilter :: Bool              -- if True, filter out missing characters
}

instance Default CharNgramParams where
    def = CharNgramParams { charAlphabet = Nothing, charSmooth = 1e-8, charLineSep = "\n", charTransform = id, charFilter = True }

trainCharNgramModel :: Int -> CharNgramParams -> String -> NgramModel Integer (Token Char)
trainCharNgramModel n params text = trainTerminatingNgramModel n (charAlphabet params) (charSmooth params) lines''
    where
        lines = splitOn (charLineSep params) text
        lines' = map (charTransform params) <$> lines
        lines'' = case (charAlphabet params) of
            Just alphabet -> if (charFilter params)
                                then filter (`S.member` alphabetSet) <$> lines'
                                else lines'
                where alphabetSet = S.fromList alphabet
            Nothing       -> lines'

max_stars_repo_name: jeremander/Parnassus

max_issues_repo_name: jeremander/Parnassus

max_forks_repo_name: jeremander/Parnassus


text:
module SimpleNet where

import System.Random
import Numeric.LinearAlgebra

import NeuralNetwork

newtype SimpleNet = SimpleNet { unSimpleNet :: Matrix R }
  deriving (Eq, Show)

makeRandomNet :: Int -> Int -> IO SimpleNet
makeRandomNet rows cols = SimpleNet <$> rand rows cols

-- |
-- >>> net = SimpleNet $ (2 >< 3) [0.47355232, 0.99773930, 0.84668094, 0.85557411, 0.03563661, 0.69422093]
-- >>> x = vector [0.6, 0.9]
-- >>> predict net x
-- [1.054148091,0.630716529,1.132807401]
predict :: SimpleNet -> Vector R -> Vector R
predict (SimpleNet w) x = x <# w

-- |
-- >>> net = SimpleNet $ (2 >< 3) [0.47355232, 0.99773930, 0.84668094, 0.85557411, 0.03563661, 0.69422093]
-- >>> x = vector [0.6, 0.9]
-- >>> t = vector [0, 0, 1]
-- >>> loss net x t
-- 9.746346258628767
loss :: SimpleNet -> Vector R -> Vector R -> Double
loss net@(SimpleNet w) x t = crossEntropyError y t
  where
    z = predict net x
    y = softmax z

max_stars_repo_name: kayhide/deep-learning-from-scratch

max_issues_repo_name: kayhide/zero-to-deep-learning

max_forks_repo_name: kayhide/zero-to-deep-learning


text:
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE OverloadedStrings         #-}
{-# LANGUAGE PolyKinds                 #-}
{-# LANGUAGE RankNTypes                #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# OPTIONS_GHC -fwarn-incomplete-patterns #-}

module Math.Regression.Regression where

import qualified Math.HMatrixUtils             as HU

import qualified Colonnade                     as C

import qualified Data.List                     as List
import qualified Data.Text                     as T
import qualified Data.Vector.Storable          as V
import qualified Lucid                         as LH
import qualified Lucid.Colonnade               as LC
import qualified Statistics.Distribution       as S
import qualified Statistics.Distribution.FDistribution
                                               as S
import qualified Statistics.Distribution.StudentT
                                               as S
import qualified Statistics.Types              as S
import qualified Text.Blaze.Colonnade          as BC
import qualified Text.Blaze.Html5              as BH
import qualified Text.Blaze.Html5.Attributes   as BHA
import qualified Text.Printf                   as TP

import           Numeric.LinearAlgebra          ( (#>)
                                                , (<.>)
                                                )
import qualified Numeric.LinearAlgebra         as LA
import           Numeric.LinearAlgebra.Data     ( Matrix
                                                , R
                                                , Vector
                                                )
import qualified Say
-- NB: for normal errors, we will request the ci and get it.  For interval errors, we may request it but we can't produce
-- what we request, we can only produce what we have.  So we make an interface for both.  We take CL as input to predict
-- *and* return it as output, in the prediction.

class Predictor e a b p where
  predict :: p -> b -> S.Estimate e a -- do we need a predict :: (b, b) -> (a, a) for uncertain inputs?

predictFromEstimateAtConfidence
  :: (RealFrac a, Predictor S.NormalErr a b p)
  => Double
  -> p
  -> S.CL Double
  -> b
  -> (a, a)
predictFromEstimateAtConfidence dof p cl b =
  let
    S.Estimate pt (S.NormalErr sigma) = predict p b
    prob =
      S.confidenceLevel $ S.mkCLFromSignificance (S.significanceLevel cl / 2)
    predCI =
      2 * S.quantile (S.studentTUnstandardized dof 0 (realToFrac sigma)) prob
  in
    (pt, realToFrac predCI)

data RegressionResult a = RegressionResult
                          {
                            parameterEstimates :: [S.Estimate S.NormalErr a]
                          , degreesOfFreedom   :: Double -- since N may be an effective N
                          , meanSquaredError   :: a
                          , rSquared           :: Double
                          , adjRSquared        :: Double
                          , fStatistic         :: Maybe Double -- we can't compute this for fits with one parameter
                          , covariances        :: Matrix Double
                          } deriving (Show)

instance (LA.Element a, LA.Numeric a, RealFloat a) => Predictor S.NormalErr a (Vector a) (RegressionResult a) where
  predict rr va =
    let mse = meanSquaredError rr
        sigmaPred = sqrt (mse + va <.> (LA.cmap realToFrac (covariances rr) #> va))
        prediction = va LA.<.> V.fromList (S.estPoint <$> parameterEstimates rr)
    in S.estimateNormErr prediction sigmaPred

data NamedEstimate a = NamedEstimate { regressorName :: T.Text, regressorEstimate :: a, regressorCI :: a, regressorPValue :: a }

namedEstimates :: [T.Text] -> RegressionResult R -> S.CL R -> [NamedEstimate R]
namedEstimates pNames res cl =
  let prob = S.confidenceLevel
        $ S.mkCLFromSignificance (S.significanceLevel cl / 2)
      dof = degreesOfFreedom res
      sigma e = S.normalError $ S.estError e
      dist e = S.studentTUnstandardized dof 0 (sigma e)
      ci e = 2 * S.quantile (dist e) prob
      pValue e = S.complCumulative (dist e) (abs $ S.estPoint e)
  in  List.zipWith (\n e -> NamedEstimate n (S.estPoint e) (ci e) (pValue e))
                   pNames
                   (parameterEstimates res)

namedEstimatesColonnade
  :: S.CL R -> C.Colonnade C.Headed (NamedEstimate R) T.Text
namedEstimatesColonnade cl =
  C.headed "parameter" regressorName
    <> C.headed "estimate" (T.pack . TP.printf "%4.3f" . regressorEstimate)
    <> C.headed
         (  T.pack (TP.printf "%.0f" (100 * S.confidenceLevel cl))
         <> "% confidence"
         )
         (T.pack . TP.printf "%4.3f" . regressorCI)
    <> C.headed "p-value" (T.pack . TP.printf "%4.3f" . regressorPValue)

namedSummaryStats :: RegressionResult R -> [(T.Text, T.Text)]
namedSummaryStats r =
  let
    p             = List.length (parameterEstimates r)
    effN          = degreesOfFreedom r + realToFrac p
    d1            = realToFrac $ p - 1
    d2            = effN - realToFrac p
    pValM = fmap (S.complCumulative (S.fDistributionReal d1 d2)) (fStatistic r)
    printNum      = T.pack . TP.printf "%4.3f"
    printMaybeNum = maybe "N/A" printNum
  in
    [ ("R-Squared"         , printNum (rSquared r))
    , ("Adj. R-squared"    , printNum (adjRSquared r))
    , ("F-stat"            , printMaybeNum (fStatistic r))
    , ("p-value"           , printMaybeNum pValM)
    , ("Mean Squared Error", printNum (meanSquaredError r))
    ]

namedSummaryStatsColonnade :: C.Colonnade C.Headed (T.Text, T.Text) T.Text
namedSummaryStatsColonnade =
  C.headed "Summary Stat." fst <> C.headed "Value" snd

prettyPrintRegressionResult
  :: T.Text -> [T.Text] -> RegressionResult R -> S.CL R -> T.Text
prettyPrintRegressionResult header xNames r cl =
  let nEsts = namedEstimates xNames r cl
      nSS   = namedSummaryStats r
  in  header
      <> "\n"
      <> T.pack (C.ascii (T.unpack <$> namedEstimatesColonnade cl) nEsts)
      <> T.pack (C.ascii (fmap T.unpack namedSummaryStatsColonnade) nSS)

prettyPrintRegressionResultLucid
  :: T.Text -> [T.Text] -> RegressionResult R -> S.CL R -> LH.Html ()
prettyPrintRegressionResultLucid header xNames r cl = do
  let nEsts = namedEstimates xNames r cl
      nSS   = namedSummaryStats r
      toCell t = LC.Cell [LH.style_ "border: 1px solid black"] (LH.toHtmlRaw t)
  LH.div_
      [ LH.style_
          "display: inline-block; padding: 7px; border-collapse: collapse"
      ]
    $ do
        LH.span_ (LH.toHtmlRaw header)
        LC.encodeCellTable
          [LH.style_ "border: 1px solid black; border-collapse: collapse"]
          (toCell <$> namedEstimatesColonnade cl)
          nEsts
        LC.encodeCellTable
          [LH.style_ "border: 1px solid black; border-collapse: collapse"]
          (fmap toCell namedSummaryStatsColonnade)
          nSS

prettyPrintRegressionResultBlaze
  :: T.Text -> [T.Text] -> RegressionResult R -> S.CL R -> BH.Html
prettyPrintRegressionResultBlaze header xNames r cl = do
  let nEsts = namedEstimates xNames r cl
      nSS   = namedSummaryStats r
      toCell t = BC.Cell (BHA.style "border: 1px solid black") (BH.toHtml t)
  BH.div
    BH.! BHA.style
           "display: inline-block; padding: 7px; border-collapse: collapse"
    $    do
           BH.span (BH.toHtml header)
           BC.encodeCellTable
             (BHA.style "border: 1px solid black; border-collapse: collapse")
             (toCell <$> namedEstimatesColonnade cl)
             nEsts
           BC.encodeCellTable
             (BHA.style "border: 1px solid black; border-collapse: collapse")
             (fmap toCell namedSummaryStatsColonnade)
             nSS

data FitStatistics a = FitStatistics { fsRSquared :: a, fsAdjRSquared :: a, fsFStatistic :: Maybe a}

goodnessOfFit
  :: MonadIO m
  => Int
  -> Vector R
  -> Maybe (Vector R)
  -> Vector R
  -> m (FitStatistics R)
goodnessOfFit pInt vB vWM vU = do
  let
    n  = LA.size vB
    p  = realToFrac pInt
    vW = maybe (V.fromList $ List.replicate n (1.0 / realToFrac n))
               (\v -> LA.scale (1 / realToFrac (LA.sumElements v)) v)
               vWM
    mW     = LA.diag vW
    vWB    = mW #> vB
--      vWU = mW #> vU
    meanWB = LA.sumElements vWB
    effN   = 1 / (vW <.> vW)
    ssTot  = let x = LA.cmap (\y -> y - meanWB) vB in x <.> (mW #> x)
    ssRes  = vU <.> (mW #> vU) -- weighted ssq of residuals
    rSq    = 1 - (ssRes / ssTot)
    arSq =
      1 - (1 - rSq) * realToFrac (effN - 1) / realToFrac (effN - p - 1)
    fStatM = if pInt > 1
      then Just (((ssTot - ssRes) / (p - 1.0)) / (ssRes / (effN - p)))
      else Nothing
  Say.say $ "n=" <> show n
  Say.say $ "p=" <> show p
--  Log.log Log.Diagnostic $ "vW=" <> show vW
  Say.say $ "effN=" <> show effN
  Say.say $ "ssTot=" <> show ssTot
  Say.say $ "ssRes=" <> show ssRes
  return $ FitStatistics rSq arSq fStatM

estimates :: Matrix R -> Vector R -> [S.Estimate S.NormalErr R]
estimates cov means =
  let sigmas = LA.cmap sqrt (LA.takeDiag cov)
  in  List.zipWith S.estimateNormErr
                   (LA.toList means)
                   (LA.toList sigmas)

eickerHeteroscedasticityEstimator
  :: MonadIO m
  => Matrix R
  -> Vector R
  -> Vector R
  -> m (Matrix R)
eickerHeteroscedasticityEstimator mA vB vB' = do
  HU.checkVectorMatrix "b" "A" vB mA
  HU.checkVectorMatrix "b'" "A" vB' mA
  let u     = vB - vB'
      diagU = LA.diag $ V.zipWith (*) u u
      mA'   = LA.tr mA
      mC    = LA.inv (mA' LA.<> mA)
  return $ mC LA.<> (mA' LA.<> diagU LA.<> mA) LA.<> mC

max_stars_repo_name: adamConnerSax/Frames-utils

max_issues_repo_name: adamConnerSax/Frames-utils

max_forks_repo_name: adamConnerSax/Frames-utils


text:
{-# LANGUAGE BangPatterns #-}
module ListStats (mean, var, stDiv, nCr,  varM, qf,unStdNorm,unStdNormSplits, -- cov, corr,
                  stDivM, stNorm, erf, cn, erfInv, binomial, cdf, cdF, meanBy, stDivBy,
                  normal, gammaln,  adjustedMatrix, adjustedMatrixN, adjustedMatrixPN, normalP,
                  module ListAlgebra,
                 ) where

import Numeric
import Numeric.LinearAlgebra
import qualified Data.Vector as V (map, fromList)
import Data.List
import Control.DeepSeq (NFData)
import LINreg2 (foldl'Rnf )
import Data.Ord (comparing)
import ListAlgebra

----------------------------------------------------------------------------------------
--- mean, variance and correlations
----------------------------------------------------------------------------------------

insertZeros :: Num a =>  [a] -> [Int] -> [a]
insertZeros xs   []   = xs
insertZeros xs (a:as) =
        let   (front,back) = splitAt (a-1) xs  in
            (1 :front) ++ (0 : insertZeros back  as)

-- map a function then calculate the standard deviation
stDivBy :: (a -> Double) -> [a] -> Double
stDivBy f xs = (sumBy (sqr . (\a -> a - ns) . f) xs) / k
    where
       (sm,k) = foldl' (\(n,m) a -> (f a + n, m + 1)) (0,0) xs
       ns     = sm / k
       sqr x  = x * x

-- map a function then calculate the mean
meanBy :: (a -> Double) -> [a] -> Double
meanBy g xs = sm / k
   where (sm,k) = foldl' (\(n,m) a -> (g a + n, m + 1)) (0,0) xs

----
mean  :: (NFData a, Fractional a) => [a] -> a
mean     xs                    = fst suml / (snd suml - 1)
                    where suml = foldl'Rnf  kF (0, 1) xs
                          kF (!n , !r) b =  (n + b, 1 + r)
---
var  :: (NFData a, Floating a) => [a] -> a
var     xs                    =  (fst suml / snd suml)
                   where suml = foldl'Rnf  kF (0, 0) xs
                         kF (!n , !r) b = (n + (b - k)^2 , 1 + r)
                         k    = mean xs
-- variance of a list of values with the mean given
varM  :: (NFData a , Floating a) => a -> [a] -> a
varM     k  xs              =  fst suml / snd suml
                 where suml = foldl'Rnf  kF (0, 0) xs
                       kF (!n , !r) b = (n + (b - k)^2 , 1 + r)
--
--- the standard deviation without the Mean given
stDiv   :: (NFData a , Floating a) =>  [a] -> a
stDiv                       = sqrt . var

---- standard deviation with the mean given
stDivM        k             = sqrt . (varM k)

--- standardise by subtracting the mean and dividing by the standard deviation
--stNorm :: [Double] -> [Double]
stNorm       xs             = map (\ a -> (a - m) / std) xs
                 where m    = mean xs
                       std  = stDivM m xs
---
unStdNorm  xs ys  = map (\a -> std * a + m)  ys
                where  m    = mean xs
                       std  = stDivM m xs

unStdNormSplits  xs zs = unStdNorm (insertZeros xs zs)

---   selecting the 'best' sublist of lenght n of a list of values
bestN :: (Ord a , Ord b)  => (a -> b) -> Int -> [a] -> [a]
bestN  comp   lim  xs
    | length xs <= lim   =  xs
    | length back <= lim =  myBest lim comp xs
    | otherwise          =  min (sortBy (comparing comp) front) (bestN comp lim back)
    where
        (front, back) = splitAt lim xs
        myBest n  f   = myTake n .  sortBy (comparing f)
        mn            = lim `div` 2

--- computed the weighted average of the residuals of values, which represent
--- points on the y-axis plotted against x-values [1 .. n], where n is the
-- length of the list. The weight biased towards is determind as the
wgtAvgRd ::  (NFData a, Enum a, Fractional a) => Int -> [[a]]  -> a
wgtAvgRd    n   xs     = weightsSum / len
    where
        (weightsSum , len) = foldl'Rnf (getWeightedAvg (2 * n)) (0,0) xs
        getWeightedAvg k (wa,ln) ys  =   (wa + wVal , ln  + 1)
            where
                f     = fromIntegral
                len   = f $ length ys
                ks    = [1 .. len]
                sumSqrs  as bs = foldl'Rnf (\a (b,c) -> a + b * c) 0 (zip as bs)
                wVal  =  (sumSqrs ks ys) / (sumSqrs ks ks)   -- (f k / len)  *

--------------------------------------------------------------------------------------------
--  Defining the X matrix for the MML function
--------------------------------------------------------------------------------------------
-- We combine adjustMatrix, fillZeros and allPieces to
-- give the final adjusted matrix. This matrix is standardised, then transposed,
-- but the transposeition can be removed if we require that to occur later on
adjustedMatrix  xs   =  [( tS . adjustMatrix . fillZeroes_aux n 0 1) yss | yss <- allPieces xs]
                         where n  = length xs
                               tS =  fromLists . map (1:) . transpose . map stNorm

adjustedMatrixN m xs =  [tS $ adjustMatrix $ fillZeroes_aux n 0 1 yss | yss <- nPieces m xs]
                         where n  = length xs
                               tS = map (1:). transpose . map stNorm
adjustedMatrixPN k m xs = [(tS $ adjustMatrix $ fillZeroes_aux n 0 1 yss, cuts) |
                                (yss,cuts) <- nPiecesP k m Nothing xs, not $ null cuts]
                         where n  = length xs
                               tS = map (1:). transpose . map stNorm
--
-- test function
joins :: Eq a => [a] -> [[a]]
joins  xs | null xs = [] -- && (length (concat a) > 5)
          | otherwise = joins_aux xs ++ joins (tail xs)
    where
        joins_aux (x : y : ys) = ([x,y] :  (map (x :) $ joins_aux (y : ys)))
        joins_aux  _           = []
--------------------------------------------------------------------------------------------
-- calculating the binomial coefficient, and normal prob density
-------------------------------------------------------------------------------------------
nCr ::  (Integral a, Ord a) => a -> a -> a
nCr   n    r   | n < 0 || r < 0  =  error "can find combinations for negative numbers"
nCr   n   0                      =  1
nCr   0   _                      =  0
nCr   !n   !k                    =  nCr (n-1) (k-1) * n `div` k

--- binomial distribution
binomial :: (Integral a, Ord b, Fractional b, Ord a) => a -> a -> b -> b
binomial   n   r   p  | p > 1 || p < 0   = error "Invalid probability value" ---n nCr r *
binomial   n   k   p  | k >= 0 && k <= n = let  nck = fromIntegral (nCr n k) in nck * (p^k) * (1 - p)^(n - k)
                      | otherwise        = 0

-- value of x from a normal distribution with mean 'mn' and standard deviation "std"
normal x mn std = (1 / (std * sqrt (2 * pi))) * exp(- 0.5 * ((x - mn)^2) / (std ^2))
--normal x mn std = cdf (x + eps) mn std  - cdf (x - eps) mn std
    where
        eps = 0.0101095
        --z   = (x - mn) / std

-- calculating the normal with a given precision
normalP eps x mn std = cdf (x + eps) mn std  - cdf (x - eps) mn std
------------------------------------------------------------------------------------------------
--- approximating the error functions and its inverse
----------------------------------------------------------------------------------------------
erf	x	| x  <  0    =  -1 *  erf (- 1 * x)
		| otherwise  = 1 - t * exp k
    where
        t  = 1 / (1 + 0.5 * abs x)
	t3 = t * (- 0.82215223 +  t *  0.17087277)
        t2 = t * ( 0.27886807 +  t * (-1.13520398 +  t *  1.48851587 + t3))
        t4 = t * ( 0.37409196 +  t * ( 0.09678418 +  t *  (- 0.18628806) + t2))
	k  = - x * x -  1.26551223 + t *  1.00002368 +  t4
--

--cn  :: Int -> Double
cn	 0	=	1
cn	 n	=	foldr (\k  b -> b + (cn k  * cn (n - 1 - k)) / ((k + 1) * (2 * k + 1))) 0 [0.. n - 1]

-- inverse of the error function
--erfInv :: (Fractional a, Floating a, Enum a) => a -> a
erfInv	y | y <  -1 || y > 1  = error "erfInv argument out of range"
          | abs y == 1        = pol (-y* log (0.0))
          | y < - 0.7         = pol (- 1 * folD z (head c) (tail c) / folD z 1 d)
          | y <   0.7         = pol ( y *  folD (y^2) (head a) (tail a) / folD (y^2) 1 b)
          | otherwise         = pol (folD z1 (head c) (tail c)  / folD z1 1 d)
                  where a     = [0.886226899, -1.645349621,  0.914624893, -0.140543331]
                        b     = [-2.118377725,  1.442710462, -0.329097515,  0.012229801]
                        c     = [-1.970840454, -1.624906493,  3.429567803,  1.641345311]
                        d     = [3.543889200,  1.637067800]
                        z     = sqrt(- log ((1.0+y)/2.0))
                        z1    = sqrt(- log ((1.0+y)/2.0))
                        pol x = x - (erf(x) - y) / (2.0 /(sqrt pi)  * exp(-x*x))
                        folD  zz a ls = foldl'Rnf  (\a b -> a + zz * (fst b + snd b * zz)) a (zipp ls)

-- commulative density function
cdf   x mean std | std < 0      = error"Negative variance to CDF"
cdf   x mean std | otherwise    =  0.5 * (1 + erf ((x - mean) / (std * sqrt 2)))

-- the CDF for the standard normal distribution
cdF  n  =  cdf  n 0 1

-- the q function, which is 1 - CDF (x) (here given for the standard normal dist)
qf  x  =  1 - cdF x
--------------------------------------------------------------------------------------
--- The gamm function takel  from wikipedia
--------------------------------------------------------------------------------------
gammaln :: Double -> Double
gammaln xx = -tmp' + log(2.5066282746310005 * ser' / xx)
			where
			tmp' = (xx+5.5) - (xx+0.5)*log(xx+5.5)
			ser' = ser + sum (zipWith (/) cof [xx+1..])
			cof  = [76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.001208650973866179,-0.000005395239384953]
			ser  = 1.000000000190015

max_stars_repo_name: rawlep/MML

max_issues_repo_name: rawlep/MML

max_forks_repo_name: rawlep/MML


text:
module Numeric.Loss where

import Dependent.Size
import Numeric.Vector.Sized

class Loss (size :: Size) fn where
  lossDerivative :: fn -> SizedArray size -> SizedArray size -> SizedArray size

max_stars_repo_name: mixed-signals/mixed-signals

max_issues_repo_name: mixed-signals/mixed-signals

max_forks_repo_name: mixed-signals/mixed-signals