text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DefaultSignatures     #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE InstanceSigs          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PolyKinds             #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-|
Module      : Grenade.Core.Network
Description : Core definition of a Neural Network
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental

This module defines the core data types and functions
for non-recurrent neural networks.
-}

module Grenade.Core.Network (
    Network (..)
  , CreatableNetwork (..)
  , Gradients (..)
  , Tapes (..)
  , GNum (..)
  , FoldableGradient (..)

  , l2Norm
  , clipByGlobalNorm
  , clipByValue
  , runNetwork
  , runGradient
  , applyUpdate
  , randomNetwork
  , randomNetworkInitWith
  ) where

import           Control.DeepSeq
import           Control.Monad.IO.Class
import           Control.Monad.Primitive          (PrimBase, PrimState)
import           Control.Parallel.Strategies
import           Data.Default
import           Data.Serialize
import           Data.Singletons
import           Data.Singletons.Prelude
import qualified Data.Vector.Storable             as V
import           GHC.TypeLits                     (KnownNat)
import           Numeric.LinearAlgebra.Static
import           System.Random.MWC

#if MIN_VERSION_base(4,9,0)
import           Data.Kind                        (Type)
#endif

import           Grenade.Core.Layer
import           Grenade.Core.NetworkInitSettings
import           Grenade.Core.NetworkSettings
import           Grenade.Core.Optimizer
import           Grenade.Core.Shape
import           Grenade.Layers.Internal.CUDA     (setCudaTriggerSize)
import           Grenade.Types

-- | Type of a network.
--
--   The @[*]@ type specifies the types of the layers.
--
--   The @[Shape]@ type specifies the shapes of data passed between the layers.
--
--   Can be considered to be a heterogeneous list of layers which are able to
--   transform the data shapes of the network.
data Network :: [Type] -> [Shape] -> Type where
    NNil  :: SingI i
          => Network '[] '[i]

    (:~>) :: (SingI i, SingI h, Layer x i h)
          => !x
          -> !(Network xs (h ': hs))
          -> Network (x ': xs) (i ': h ': hs)
infixr 5 :~>

instance Show (Network '[] '[i]) where
  show NNil = "NNil"
instance (Show x, Show (Network xs rs)) => Show (Network (x ': xs) (i ': rs)) where
  show (x :~> xs) = show x ++ " ~> " ++ show xs

instance NFData (Network '[] '[ i]) where
  rnf NNil = ()
instance (NFData x, NFData (Network xs rs)) => NFData (Network (x ': xs) (i ': rs)) where
  rnf ((!x) :~> (!xs)) = rnf x `seq` rnf xs

-- | Gradient of a network.
--
--   Parameterised on the layers of the network.
data Gradients :: [Type] -> Type where
   GNil  :: Gradients '[]

   (:/>) :: UpdateLayer x
         => !(Gradient x)
         -> !(Gradients xs)
         -> Gradients (x ': xs)

instance NFData (Gradients '[]) where
  rnf GNil = ()
instance (NFData (Gradient x), NFData (Gradients xs)) => NFData (Gradients (x ': xs)) where
  rnf (g :/> gs) = rnf g `seq` rnf gs

instance Serialize (Gradients '[]) where
  put GNil = put ()
  get = return GNil
instance (UpdateLayer x, Serialize (Gradient x), Serialize (Gradients xs)) => Serialize (Gradients (x ': xs)) where
  put (g :/> gs) = put g >> put gs
  get = (:/>) <$> get <*> get

-- | Wegnert Tape of a network.
--
--   Parameterised on the layers and shapes of the network.
data Tapes :: [Type] -> [Shape] -> Type where
   TNil  :: SingI i
         => Tapes '[] '[i]

   (:\>) :: (SingI i, SingI h, Layer x i h)
         => !(Tape x i h)
         -> !(Tapes xs (h ': hs))
         -> Tapes (x ': xs) (i ': h ': hs)

instance NFData (Tapes '[] '[i]) where
  rnf TNil       = ()

instance (NFData (Tape x i h), NFData (Tapes xs (h ': hs))) => NFData (Tapes (x ': xs) (i ': h ': hs)) where
  rnf (t :\> ts) = rnf t `seq` rnf ts

-- | Running a network forwards with some input data.
--
--   This gives the output, and the Wengert tape required for back
--   propagation.
runNetwork :: forall layers shapes.
              Network layers shapes
           -> S (Head shapes)
           -> (Tapes layers shapes, S (Last shapes))
runNetwork = go
    where
  go  :: forall js ss. (Last js ~ Last shapes)
      => Network ss js
      -> S (Head js)
      -> (Tapes ss js, S (Last js))
  go (layer :~> n) !x =
    let (tape, forward) =

          runForwards layer x
        (tapes, answer) = go n forward
    in  (tape :\> tapes, answer)

  go NNil !x
      = (TNil, x)

-- | Running a loss gradient back through the network.
--
--   This requires a Wengert tape, generated with the appropriate input
--   for the loss.
--
--   Gives the gradients for the layer, and the gradient across the
--   input (which may not be required).
runGradient :: forall layers shapes.
               Network layers shapes
            -> Tapes layers shapes
            -> S (Last shapes)
            -> (Gradients layers, S (Head shapes))
runGradient net tapes o =
  go net tapes
    where
  go  :: forall js ss. (Last js ~ Last shapes)
      => Network ss js
      -> Tapes ss js
      -> (Gradients ss, S (Head js))
  go (layer :~> n) (tape :\> nt) =
    let (gradients, feed)  = go n nt
        (layer', backGrad) = runBackwards layer tape feed
    in  (layer' :/> gradients, backGrad)

  go NNil TNil
      = (GNil, o)

-- | Apply one step of stochastic gradient descent across the network.
applyUpdate :: Optimizer opt -> Network layers shapes -> Gradients layers -> Network layers shapes
applyUpdate rate (layer :~> rest) (gradient :/> grest) =
  let layer' = runUpdate rate layer gradient
      rest' = applyUpdate rate rest grest `using` rpar
   in layer' :~> rest'
applyUpdate _ NNil GNil = NNil

-- | Apply network settings across the network.
applySettingsUpdate :: NetworkSettings -> Network layers shapes -> Network layers shapes
applySettingsUpdate settings (layer :~> rest) =
  let layer' = runSettingsUpdate settings layer
      layers' = applySettingsUpdate settings rest `using` rpar
   in layer' :~> layers'
applySettingsUpdate _ NNil = NNil

-- | A network can easily be created by hand with (:~>), but an easy way to
--   initialise a random network is with the @randomNetworkWith@ function.
class CreatableNetwork (xs :: [Type]) (ss :: [Shape])
  -- | Create a network with randomly initialised weights.
  --
  --   Calls to this function will not compile if the type of the neural
  --   network is not sound.
  where
  randomNetworkWith :: PrimBase m => NetworkInitSettings -> Gen (PrimState m) -> m (Network xs ss)

-- | Create a random network using uniform distribution.
randomNetwork :: (MonadIO m, CreatableNetwork xs ss) => m (Network xs ss)
randomNetwork = randomNetworkInitWith def

-- | Create a random network using the specified weight initialization method.
randomNetworkInitWith :: (MonadIO m, CreatableNetwork xs ss) => NetworkInitSettings -> m (Network xs ss)
randomNetworkInitWith m = liftIO $ withSystemRandom . asGenST $ \gen -> randomNetworkWith m gen

instance SingI i => CreatableNetwork '[] '[i] where
  randomNetworkWith initCfg  _ = setCudaTriggerSize (gpuTriggerSize initCfg) >> return NNil

instance (SingI i, SingI o, Layer x i o, RandomLayer x, CreatableNetwork xs (o ': rs)) => CreatableNetwork (x ': xs) (i ': o ': rs) where
  randomNetworkWith m gen = (:~>) <$> createRandomWith m gen <*> randomNetworkWith m gen

-- | Add very simple serialisation to the network
instance SingI i => Serialize (Network '[] '[i]) where
  put NNil = pure ()
  get = return NNil

instance (SingI i, SingI o, Layer x i o, Serialize x, Serialize (Network xs (o ': rs))) => Serialize (Network (x ': xs) (i ': o ': rs)) where
  put (x :~> r) = put x >> put r
  get = (:~>) <$> get <*> get

-- | Ultimate composition.
--
--   This allows a complete network to be treated as a layer in a larger network.
instance UpdateLayer (Network sublayers subshapes) where
  type Gradient (Network sublayers subshapes) = Gradients sublayers
  runUpdate = applyUpdate
  runSettingsUpdate = applySettingsUpdate

instance FoldableGradient (Gradients '[]) where
  mapGradient _ GNil = GNil
  squaredSums GNil = []

instance (NFData (Gradient x), NFData (Gradients xs), FoldableGradient (Gradient x), FoldableGradient (Gradients xs)) => FoldableGradient (Gradients (x ': xs)) where
  mapGradient f (x :/> xs) =
    let x' = mapGradient f x `using` rparWith rdeepseq
        xs' = mapGradient f xs `using` rparWith rdeepseq
     in x' :/> xs'
  squaredSums (x :/> xs) = (squaredSums x `using` rparWith rdeepseq) ++ (squaredSums xs `using` rparWith rdeepseq)

-- | Get the L2 Norm of a Foldable Gradient.
l2Norm :: (FoldableGradient x) => x -> RealNum
l2Norm grad = sqrt . sum $ squaredSums grad

-- | Clip the gradients by the global norm.
clipByGlobalNorm :: (FoldableGradient (Gradients xs)) => RealNum -> Gradients xs -> Gradients xs
clipByGlobalNorm c grads =
  let divisor = sqrt . sum $ squaredSums grads
   in if divisor > c
        then mapGradient (* (c / divisor)) grads
        else grads

-- | Clip the gradients by the value v. Preforms @min v (max(-v) x)@ on all gradients.
clipByValue :: (FoldableGradient (Gradients xs)) => RealNum -> Gradients xs -> Gradients xs
clipByValue v grads = mapGradient (min v . (max (-v))) grads

instance CreatableNetwork sublayers subshapes => RandomLayer (Network sublayers subshapes) where
  createRandomWith = randomNetworkWith

-- | Ultimate composition.
--
--   This allows a complete network to be treated as a layer in a larger network.
instance (i ~ (Head subshapes), o ~ (Last subshapes)) => Layer (Network sublayers subshapes) i o where
  type Tape (Network sublayers subshapes) i o = Tapes sublayers subshapes
  runForwards  = runNetwork
  runBackwards = runGradient

--------------------------------------------------

-- | Grenade Num class.
--
-- This allows for instance scalar multiplication of the weights, which is useful for slowly adapting networks, e.g. NN'
-- <- \tau * NN' + (1-\tau) * NN. Or one could sum up some gradients in parallel and apply them at once: @applyUpdate lp
-- net $ foldl1 (|+) ...@aq.
class GNum a where
  (|*) :: Rational -> a -> a
  (|+) :: a -> a -> a
  sumG :: [a] -> a
  default sumG :: [a] -> a
  sumG = foldl1 (|+)

infixl 7 |*
infixr 5 |+

instance (SingI i) => GNum (Network '[] '[ i]) where
  _ |* NNil = NNil
  _ |+ NNil = NNil
  sumG _ = NNil

instance (SingI i, SingI o, Layer x i o, NFData x, NFData (Network xs (o ': rs)), GNum x, GNum (Network xs (o ': rs))) => GNum (Network (x ': xs) (i ': o ': rs)) where
  s |* (x :~> xs) =
    let x' = (s |* x)
        xs' = (s |* xs) `using` rparWith rdeepseq
     in x' :~> xs'
  (x :~> xs) |+ (y :~> ys) =
    let x' = (x |+ y)
        xs' = (xs |+ ys) `using` rparWith rdeepseq
     in x' :~> xs'
  sumG xs = sumG (map (\(l :~> _) -> l) xs) :~> (sumG (map (\(_ :~> ls) -> ls) xs) `using` rparWith rdeepseq)

instance GNum (Gradients '[]) where
  _ |* GNil = GNil
  _ |+ GNil = GNil
  sumG _ = GNil

instance (GNum a) => GNum [a] where
  r |* xs = fmap (r |*) xs
  xs |+ ys = zipWith (|+) xs ys
  sumG = map sumG

instance (UpdateLayer x, GNum (Gradient x), GNum (Gradients xs), NFData (Gradient x), NFData (Gradients xs)) => GNum (Gradients (x ': xs)) where
  s |* (x :/> xs) =
    let x' = (s |* x)
        xs' = (s |* xs) `using` rparWith rdeepseq
     in x' :/> xs'
  (x :/> xs) |+ (y :/> ys) =
    let x' = (x |+ y)
        xs' = (xs |+ ys) `using` rparWith rdeepseq
     in x' :/> xs'
  sumG xs = sumG (map (\(x :/> _) -> x) xs) :/> (sumG (map (\(_ :/> xs') -> xs') xs) `using` rparWith rdeepseq)

instance GNum () where
  _ |* () = ()
  _ |+ () = ()
  sumG _ = ()

instance (GNum a, GNum b) => GNum (a, b) where
  s |* (a, b) = (s |* a, s |* b)
  (a1, b1) |+ (a2, b2) = (a1 |+ a2, b1 |+ b2)
  sumG xs = (sumG (map fst xs), sumG (map snd xs))

instance (KnownNat m) => GNum (R m) where
  s |* vec = dvmap (fromRational s *) vec
  (|+) = (+)
  sumG xs = sum xs

instance (KnownNat m, KnownNat n) => GNum (L m n) where
  s |* mat = dmmap (fromRational s *) mat
  (|+) = (+)
  sumG xs = sum xs

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
{-# LANGUAGE CPP, UndecidableInstances, TypeSynonymInstances, FlexibleInstances, OverloadedStrings, DefaultSignatures, PatternSynonyms, ViewPatterns, TypeOperators, FlexibleContexts #-}
module Pure.Data.Cond where

import Pure.Data.Default
import Pure.Data.JSON
import Pure.Data.Lifted
import Pure.Data.Try
import Pure.Data.Txt as Txt
import Pure.Data.View

import qualified Data.IntMap as IntMap
import qualified Data.IntSet as IntSet
import qualified Data.Map as Map
import qualified Data.Sequence as Seq
import qualified Data.Set as Set

import qualified Data.HashMap.Strict as HM

import qualified Data.Vector as Vector

import Control.Applicative
import Control.Exception
import Data.Complex
import Data.Int
import Data.Monoid
import Data.Ratio
import Data.Word

import GHC.Generics

-- | Cond, for conditional, is a class used for representing, and checking for,
-- nil values.
--
-- Cond is lawless, but instances should satisfy the ideal that a nil value
-- represents either failure, or as little information as possible.
--
-- Caveat emptor: Cond usage /will/ introduce bugs.
--
-- For `Any` and `All`, Cond's nil represents failure:
--
-- > (mempty :: Any) == Any False == (nil :: Any)
--
-- > (mempty :: All) == All True /= All False == (nil :: All)
--
-- For numbers, Cond's nil is zero.
--
-- For containers, Cond's nil is the empty container.
--
-- For text, Cond's nil is the empty text.
--
-- For Product and Sum, Cond's nil is the lited nil.
--
-- > (nil :: Cond a => Product a) == Product nil
--
-- > (nil :: Cond a => Sum a) == Sum nil
--
-- For Ordering:
--
-- > (nil :: Ordering) == LT
--
-- For Either:
--
-- > (nil :: Cond a => Either a b) == Left nil
--
-- For SomeException:
--
-- > (nil :: SomeException) == toException GenericException
-- > (isNil :: SomeException -> Bool) == True
-- > (notNil :: SomeException -> Bool) == False
--
class Cond a where
  nil :: a
  default nil :: (Generic a, GCond (Rep a)) => a
  nil = to gnil

  isNil :: a -> Bool
  default isNil :: (Generic a, GCond (Rep a)) => a -> Bool
  isNil = gisNil . from

  notNil :: a -> Bool
  notNil = not . isNil

infixr 1 ?
(?) :: (Cond x) => x -> a -> a -> a
(?) x t e = if notNil x then t else e

infixr 1 !?
(!?) :: (Cond x) => x -> a -> a -> a
(!?) x t e = if isNil x then t else e

may :: Cond a => (b -> a) -> Maybe b -> a
may = maybe nil

infixr 6 #
(#) :: (Cond x, Default a) => x -> a -> a
(#) b t = b Pure.Data.Cond.? t $ def

cond :: (Cond a, Default a) => a -> a
cond a = a # a

instance Cond [a] where
  nil = []
  isNil [] = True
  isNil _ = False

#ifdef __GHCJS__
instance Cond JSV where
  nil = nullJSV
  isNil = isNull
#endif

instance Cond a => Cond (Either a b) where
  nil = Left nil
  isNil (Left a) = isNil a
  isNil _ = False

data Nil = Nil deriving (Show)
instance Exception Nil

instance Cond SomeException where
  nil = toException Nil
  isNil _ = True
  notNil _ = False

instance Cond (Maybe a) where
  nil = Nothing
  isNil Nothing = True
  isNil _ = False

instance Cond (Try a) where
  nil = Failed
  isNil Failed = True
  isNil _ = False

instance Cond () where
  nil = ()
  isNil _ = True

instance Cond Any where
  nil = Any False
  isNil (Any a) = not a

instance Cond All where
  nil = All False
  isNil (All a) = not a

instance Cond Ordering where
  nil = LT
  isNil LT = True
  isNil _ = False

instance Cond (Last a) where
  nil = Last Nothing
  isNil (Last Nothing) = True
  isNil _ = False

instance Cond (First a) where
  nil = First Nothing
  isNil (First Nothing) = True
  isNil _ = False

instance (Cond a) => Cond (Sum a) where
  nil = Sum nil
  isNil (Sum a) = isNil a

instance (Cond a) => Cond (Product a) where
  nil = Product nil
  isNil (Product a) = isNil a

instance (Cond a) => Cond (Const a b) where
  nil = Const nil
  isNil (Const a) = isNil a

instance Cond (ZipList a) where
  nil = ZipList []
  isNil (ZipList []) = True
  isNil _ = False

instance (Cond a) => Cond (Dual a) where
  nil = Dual nil
  isNil (Dual d) = isNil d

instance Cond Bool where
  nil = False
  isNil = not

instance Cond Txt where
  nil = ""
  isNil = Txt.null

instance {-# OVERLAPPABLE #-} (Num n, Eq n) => Cond n where
  nil = fromInteger 0
  isNil = (== fromInteger 0)

instance (Cond a, Cond b) => Cond (a,b) where
  nil = (nil,nil)
  isNil (a,b) = isNil a && isNil b

instance (Cond a, Cond b, Cond c) => Cond (a,b,c) where
  nil = (nil,nil,nil)
  isNil (a,b,c) = isNil a && isNil b && isNil c

instance (Cond a, Cond b, Cond c, Cond d) => Cond (a,b,c,d) where
  nil = (nil,nil,nil,nil)
  isNil (a,b,c,d) = isNil a && isNil b && isNil c && isNil d

instance (Cond a, Cond b, Cond c, Cond d, Cond e) => Cond (a,b,c,d,e) where
  nil = (nil,nil,nil,nil,nil)
  isNil (a,b,c,d,e) = isNil a && isNil b && isNil c && isNil d && isNil e

instance (Cond a, Cond b, Cond c, Cond d, Cond e, Cond f) => Cond (a,b,c,d,e,f) where
  nil = (nil,nil,nil,nil,nil,nil)
  isNil (a,b,c,d,e,f) = isNil a && isNil b && isNil c && isNil d && isNil e && isNil f

instance (Cond a, Cond b, Cond c, Cond d, Cond e, Cond f, Cond g) => Cond (a,b,c,d,e,f,g) where
  nil = (nil,nil,nil,nil,nil,nil,nil)
  isNil (a,b,c,d,e,f,g) = isNil a && isNil b && isNil c && isNil d && isNil e && isNil f && isNil g

instance Cond (Set.Set a) where
  nil = Set.empty
  isNil = Set.null

instance Cond (Seq.Seq a) where
  nil = Seq.empty
  isNil = Seq.null

instance Cond (Map.Map a b) where
  nil = Map.empty
  isNil = Map.null

instance Cond (IntMap.IntMap a) where
  nil = IntMap.empty
  isNil = IntMap.null

instance Cond IntSet.IntSet where
  nil = IntSet.empty
  isNil = IntSet.null

instance {-# OVERLAPPABLE #-} Cond (HM.HashMap a b) where
  nil = HM.empty
  isNil = HM.null

instance Cond (Vector.Vector a) where
  nil = Vector.empty
  isNil = Vector.null

instance {-# OVERLAPPING #-} Cond Obj where
  nil = mempty
  isNil o =
#ifdef __GHCJS__
    Prelude.null $ objectAssocs o
#else
    o == mempty
#endif

instance {-# OVERLAPPING #-} Cond Value where
  nil =
#ifdef __GHCJS__
    nullValue
#else
    Null
#endif
  isNil =
#ifdef __GHCJS__
    (== nullValue)
#else
    (== Null)
#endif

instance Cond View where
  nil = NullView Nothing
  isNil (NullView _) = True
  isNil _ = False

class GCond f where
  gnil :: f a
  gisNil :: f a -> Bool

instance GCond V1 where
  gnil = undefined
  gisNil _ = True

instance GCond U1 where
  gnil = U1
  gisNil U1 = True

instance (Cond a) => GCond (K1 i a) where
  gnil = K1 nil
  gisNil (K1 a) = isNil a

instance (GCond a, GCond b) => GCond (a :*: b) where
  gnil = gnil :*: gnil
  gisNil (a :*: b) = gisNil a && gisNil b

instance (GCond a, GCond b) => GCond (a :+: b) where
  gnil = L1 gnil
  gisNil (L1 n) = gisNil n
  gisNil _ = False

instance (GCond a) => GCond (M1 i c a) where
  gnil = M1 gnil
  gisNil (M1 a) = gisNil a

ext: hs

max_stars_repo_name: grumply/pure-cond

max_issues_repo_name: grumply/pure-cond

max_forks_repo_name: grumply/pure-cond


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

{-# LANGUAGE AllowAmbiguousTypes   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedLabels      #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-|
Module      : Grenade.Layers.Add
Description : Layer to add some constant, with shape broadcasting
Maintainer  : Theo Charalambous
License     : BSD2
Stability   : experimental
-}

module Grenade.Layers.Add
  (
  -- * Layer instances
    Add (..)

  -- * Helper functions
  , initAdd
  ) where

import           Control.DeepSeq              (NFData (..))

import           Data.Kind                    (Type)
import           Data.Maybe                   (fromJust)
import           Data.Proxy
import           Data.Serialize
import           GHC.TypeLits

import qualified Numeric.LinearAlgebra        as LA
import           Numeric.LinearAlgebra.Static (R)
import qualified Numeric.LinearAlgebra.Static as H

import           Lens.Micro                   ((^.))

import           Grenade.Core
import           Grenade.Layers.Internal.Add
import           Grenade.Onnx

data Add :: Nat -- The number of channels of the bias
         -> Nat -- The number of rows of the bias
         -> Nat -- The number of columns of the bias
         -> Type where
  Add  :: ( KnownNat channels
          , KnownNat rows
          , KnownNat columns )
          => R (channels * rows * columns)
          -> Add channels rows columns

instance UpdateLayer (Add c h w) where
  type Gradient (Add c h w) = ()
  runUpdate _ x _  = x
  reduceGradient _ = ()

instance (KnownNat c, KnownNat h, KnownNat w ) => RandomLayer (Add c h w) where
  createRandomWith _ _ = pure initAdd

-- | Initialize an Add layer with 0 bias
initAdd :: forall c h w. ( KnownNat c, KnownNat h, KnownNat w )
        => Add c h w
initAdd =
  let c'      = fromIntegral $ natVal (Proxy :: Proxy c)
      h'      = fromIntegral $ natVal (Proxy :: Proxy h)
      w'      = fromIntegral $ natVal (Proxy :: Proxy w)
      zeroes  = replicate (c' * h' * w') 0
      bias    = H.fromList zeroes   :: R (c * h * w)
  in Add bias

instance ( KnownNat c, KnownNat h, KnownNat w ) => Serialize (Add c h w) where
  put (Add bias) = putListOf put . LA.toList . H.extract $ bias
  get            = do
    bias <- maybe (fail "Vector of incorrect size") return . H.create . LA.fromList =<< getListOf get
    return $ Add bias

-- | A three dimensions image can have a vector added to it with size equal to the number of channels
--   of the image, this corresponds to adding the nth element of the bias to every pixel of the nth
--   channel of the input image
instance ( KnownNat i, KnownNat j, KnownNat k ) => Layer (Add k 1 1) ('D3 i j k) ('D3 i j k) where
  type Tape (Add k 1 1) ('D3 i j k) ('D3 i j k) = ()

  runForwards (Add b) (S3D m)
    = let c  = fromIntegral $ natVal (Proxy :: Proxy k)
          h  = fromIntegral $ natVal (Proxy :: Proxy i)
          w  = fromIntegral $ natVal (Proxy :: Proxy j)
          m' = H.extract m
          b' = H.extract b
          r  = addPerChannel c h w m' b'
      in  ((), S3D . fromJust . H.create $ r)

  runBackwards = undefined

instance OnnxOperator (Add c h w) where
  onnxOpTypeNames _ = ["Add"]

instance (KnownNat c, KnownNat h, KnownNat w) => OnnxLoadable (Add c h w) where
  loadOnnxNode inits node = case node ^. #input of
    [bias, _] -> do
      loadedBias <- readInitializerTensorIntoVector inits bias
      return $ Add loadedBias
    _ -> onnxIncorrectNumberOfInputs

instance NFData (Add c h w) where
  rnf (Add bias) = rnf bias

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE OverloadedStrings #-}

module GameData
    ( ResultSimple (ResultSimple, evSimpleCol, evSimpleRow, sdSimpleCol, sdSimpleRow, weightsColsSimple, weightsRowsSimple)
    , GameSimple (GameSimple, gameSName, gameSCols, gameSRows, gameSMatrixCol, gameSMatrixRow, gameSOutcomes)
    , Game (Game, gameCName, attCName, defCName, gameData, outcomesC)
    , Result (Result, evc, evr, sdc, sdr, weightsAtts, weightsDefs)
    , Opt
    , OptOutcome
    ) where

import Data.Function
import Data.List
import Data.Maybe
import Data.Text (Text, unpack)
import Numeric.LinearAlgebra hiding (build)
import Numeric.LinearAlgebra.Data hiding (build)
import Formatting
import Formatting.Formatters hiding (build)
import qualified Formatting.Formatters as F
import Formatting.Buildable

data ResultSimple = ResultSimple {evSimpleCol::Double, evSimpleRow::Double, sdSimpleCol::Double, sdSimpleRow::Double, weightsColsSimple::[Double], weightsRowsSimple::[Double]}
    deriving (Show)

data GameSimple = GameSimple {gameSName::Text, gameSCols::[Text], gameSRows::[Text], gameSMatrixCol::(Matrix Double), gameSMatrixRow::(Matrix Double), gameSOutcomes::ResultSimple}
instance Buildable GameSimple where
    build g = let
                  (ResultSimple evc evr sdc sdr wc wr) = gameSOutcomes g
                  n = gameSName g
                  m1 = gameSMatrixCol g
                  m2 = gameSMatrixRow g
                  c = zip (gameSCols g) wc
                  r = zip (gameSRows g) wr
              in
                  bformat (stext % ":\n Column's Game: " % shown % "\n Row's Game: " % shown % "\n Column player: " % shown % "\n Row player: " % shown % "\n EVs: " % shown % "\n SDs: " % shown) n m1 m2 c r (evc, evr) (sdc, sdr)
instance Show GameSimple where
    show g = formatToString F.build g

type Opt = (Text, Maybe Double)
type OptOutcome = (Opt, Opt, Double, Double)

data Result = Result {evc::Double, evr::Double, sdc::Double, sdr::Double, weightsAtts::[(Text,Double)], weightsDefs::[(Text,Double)]}
    deriving (Eq, Ord)
instance Show Result where
    show r = "\n EVs: " ++ (show (evc r, evr r)) ++ "\n SDs: " ++ (show (sdc r, sdr r)) ++ "\n Attacker Options: " ++ (show . weightsAtts $ r) ++ "\n Defender Options: " ++ (show . weightsDefs $ r)

data Game = Game {gameCName::Text, attCName::Text, defCName::Text, gameData::[OptOutcome], outcomesC::Result}
instance Ord Game where
    compare = compare `on` gameCName
instance Eq Game where
    a == b = and [gameCName a == gameCName b, attCName a == attCName b, defCName a == defCName b]
instance Show Game where
    show (Game "" _ _ gdata gout) = show gout
    show (Game gname "" "" gdata gout) = (unpack gname) ++ "\n" ++ (show gout)
    show (Game gname attname defname gdata gout) = (unpack gname) ++ " (" ++ (unpack attname) ++ " vs " ++ (unpack defname) ++ ")" ++ (show gout)

max_stars_repo_name: StaccatoSemibreve/Mixup-Analyser

max_issues_repo_name: StaccatoSemibreve/Mixup-Analyser

max_forks_repo_name: StaccatoSemibreve/Mixup-Analyser


text:
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE FlexibleInstances                        #-}
{-# LANGUAGE MultiParamTypeClasses                    #-}
{-# LANGUAGE PatternSynonyms                          #-}
{-# LANGUAGE RankNTypes                               #-}
{-# LANGUAGE TypeFamilies                             #-}
{-# LANGUAGE TypeOperators                            #-}
{-# LANGUAGE UndecidableInstances                     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

module Backprop.Learn.Model.Neural (
  -- * Fully connected
  -- ** Feed-forward
    FC, pattern FC, FCp, fcBias, fcWeights
  -- *** With activation function
  , FCA, pattern FCA, _fcaActivation
  -- ** Recurrent
  , FCR, pattern FCR, FCRp, fcrBias, fcrInputWeights, fcrStateWeights
  -- *** With activation function
  , FCRA, pattern FCRA, _fcraStore, _fcraActivation
  ) where

import           Backprop.Learn.Model.Combinator
import           Backprop.Learn.Model.Regression
import           Backprop.Learn.Model.State
import           Data.Tuple
import           GHC.TypeNats
import           Lens.Micro
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop
import qualified Numeric.LinearAlgebra.Static          as H

-- | Fully connected feed-forward layer with bias.  Parameterized by its
-- initialization distribution.
--
-- Note that this has no activation function; to use as a model with
-- activation function, chain it with an activation function using 'RMap',
-- ':.~', etc.; see 'FCA' for a convenient type synonym and constructor.
--
-- Without any activation function, this is essentially a multivariate
-- linear regression.
--
-- With the logistic function as an activation function, this is
-- essentially multivariate logistic regression. (See 'Logistic')
type FC i o = LinReg i o

pattern FC :: FC i o
pattern FC = LinReg

-- | Convenient synonym for an 'FC' post-composed with a simple
-- parameterless activation function.
type FCA i o = RMap (R o) (R o) (FC i o)

-- | Construct an 'FCA' using a generating function and activation
-- function.
--
-- Some common ones include 'logistic' and @'vmap' 'reLU'@.
pattern FCA
    :: (forall s. Reifies s W => BVar s (R o) -> BVar s (R o))          -- ^ '_fcaActivation'
    -> FCA i o
pattern FCA { _fcaActivation } = RM _fcaActivation FC

type FCp = LRp

fcWeights :: Lens (FCp i o) (FCp i' o) (L o i) (L o i')
fcWeights = lrBeta

fcBias :: Lens' (FCp i o) (R o)
fcBias = lrAlpha

-- | Fully connected recurrent layer with bias.
--
-- Parameterized by its initialization distributions, and also the function
-- to compute the new state from previous input.
--
-- @
-- instance 'Learn' ('R' i) (R o) ('FCR' h i o) where
--     type 'LParamMaybe' (FCR h i o) = ''Just' ('FCRp' h i o)
--     type 'LStateMaybe' (FCR h i o) = 'Just (R h)
-- @
type FCR h i o = Recurrent (R (i + h)) (R i) (R h) (R o) (FC (i + h) o)

-- | Construct an 'FCR'
pattern FCR
    :: (KnownNat h, KnownNat i)
    => (forall s. Reifies s W => BVar s (R o) -> BVar s (R h))          -- ^ '_fcrSTore'
    -> FCR h i o
pattern FCR { _fcrStore } <-
    Rec { _recLoop  = _fcrStore
        , _recLearn = FC
        }
  where
    FCR s = Rec { _recSplit = H.split
                , _recJoin  = (H.#)
                , _recLoop  = s
                , _recLearn = FC
                }
{-# COMPLETE FCR #-}

-- | Convenient synonym for an 'FCR' post-composed with a simple
-- parameterless activation function.
type FCRA h i o = RMap (R o) (R o) (FCR h i o)

-- | Construct an 'FCRA' using a generating function and activation
-- function.
--
-- Some common ones include 'logistic' and @'vmap' 'reLU'@.
pattern FCRA
    :: (KnownNat h, KnownNat i)
    => (forall s. Reifies s W => BVar s (R o) -> BVar s (R h))          -- ^ '_fcraStore'
    -> (forall s. Reifies s W => BVar s (R o) -> BVar s (R o))          -- ^ '_fcraActivation'
    -> FCRA h i o
pattern FCRA { _fcraStore, _fcraActivation }
    = RM _fcraActivation (FCR _fcraStore)

type FCRp h i o = FCp (i + h) o

lensIso :: (s -> (a, x)) -> ((b, x) -> t) -> Lens s t a b
lensIso f g h x = g <$> _1 h (f x)

fcrInputWeights
    :: (KnownNat h, KnownNat i, KnownNat i', KnownNat o)
    => Lens (FCRp h i o) (FCRp h i' o) (L o i) (L o i')
fcrInputWeights = fcWeights . lensIso H.splitCols (uncurry (H.|||))

fcrStateWeights
    :: (KnownNat h, KnownNat h', KnownNat i, KnownNat o)
    => Lens (FCRp h i o) (FCRp h' i o) (L o h) (L o h')
fcrStateWeights = fcWeights . lensIso (swap . H.splitCols) (uncurry (H.|||) . swap)

fcrBias :: Lens' (FCRp h i o) (R o)
fcrBias = fcBias

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

module Test.Grenade.Layers.PadCrop where

import           Grenade

import           Hedgehog

import           Numeric.LinearAlgebra.Static ( norm_Inf )

import           Test.Hedgehog.Hmatrix

import           Data.Serialize
import           Data.Either
import           System.Random.MWC             (create)

type PadCropNet3D = Network '[Pad 2 3 4 6, Crop 2 3 4 6] '[ 'D3 7 9 5, 'D3 16 15 5, 'D3 7 9 5 ]
type PadCropNet2D = Network '[Pad 2 3 4 6, Crop 2 3 4 6] '[ 'D2 7 9, 'D2 16 15, 'D2 7 9 ]

prop_pad_crop :: Property
prop_pad_crop =
  let net :: PadCropNet3D
      net = Pad :~> Crop :~> NNil
  in  property $
    forAll genOfShape >>= \(d :: S ('D3 7 9 5)) ->
      let (tapes, res)  = runForwards  net d
          (_    , grad) = runBackwards net tapes d
      in  do assert $ d ~~~ res
             assert $ grad ~~~ d

prop_pad_crop_2d :: Property
prop_pad_crop_2d =
  let net :: Network '[Pad 2 3 4 6, Crop 2 3 4 6] '[ 'D2 7 9, 'D2 16 15, 'D2 7 9 ]
      net = Pad :~> Crop :~> NNil
  in  property $
    forAll genOfShape >>= \(d :: S ('D2 7 9)) ->
      let (tapes, res)  = runForwards  net d
          (_    , grad) = runBackwards net tapes d
      in  do assert $ d ~~~ res
             assert $ grad ~~~ d

prop_pad_crop_is_serializable :: Property
prop_pad_crop_is_serializable = withTests 1 $ property $ do
  gen <- evalIO create
  pad  :: Pad 2 3 4 6  <- evalIO $ createRandomWith UniformInit gen
  crop :: Crop 2 3 4 6 <- evalIO $ createRandomWith UniformInit gen
  let net :: PadCropNet3D
      net = pad :~> crop :~> NNil
      bs = encode net
      dec = decode bs :: Either String PadCropNet3D
  assert $ isRight dec

prop_can_show_pad_crop :: Property
prop_can_show_pad_crop =
  let net :: PadCropNet3D
      net = Pad :~> Crop :~> NNil
  in (withTests 1 . property) $ show net `seq` success

prop_can_update_pad_crop_and_use_in_batches :: Property
prop_can_update_pad_crop_and_use_in_batches =
  let net :: PadCropNet3D
      net = Pad :~> Crop :~> NNil
  in  property $
    forAll genOfShape >>= \(d :: S ('D3 7 9 5)) ->
      let (tapes, _) = runForwards  net d
          (v    , _)   = runBackwards net tapes d
      in do
        runUpdate defSGD net v `seq` success
        runUpdate defAdam net v `seq` success
        reduceGradient @PadCropNet3D [v] `seq` success

(~~~) :: S x -> S x -> Bool
(S1D x) ~~~ (S1D y) = norm_Inf (x - y) < 0.00001
(S2D x) ~~~ (S2D y) = norm_Inf (x - y) < 0.00001
(S3D x) ~~~ (S3D y) = norm_Inf (x - y) < 0.00001

tests :: IO Bool
tests = checkParallel $$(discover)


text:
module Seed
  ( Seed
  , createSeed
  , arrowsAt
  , drawSeed
  )
where

import           Arrow                          ( Arrow
                                                , arrange
                                                , endpoint
                                                , updateArrow
                                                )
import qualified Fourier
import           Data.Complex
import           Graphics.Gloss

-- A seed represents an initial state of arrow frequencies, directions, and magnitudes
-- Any complex function can be traced out by changing this initial state.
-- It also contains some drawing and scaling information.
data Seed = Seed { _arrows :: [(Arrow, Integer)], _spaceScale :: Float, _timeScale :: Float, _trailLength :: Float }

createSeed :: (Float -> Complex Float) -> Int -> Float -> Float -> Float -> Seed
createSeed f n spaceScale timeScale trailLength = Seed
  { _arrows      = arrows
  , _spaceScale  = spaceScale
  , _timeScale   = timeScale
  -- measure trail length by percentage of a time-scaled second
  , _trailLength = trailLength
  }
 where
  -- order arrow frequencies 0, 1, -1, 2, -2, 3, -3, etc.
  freqs             = take n (0 : concat [ [x, -x] | x <- [1 ..] ])
  startingPositions = map (Fourier.coefficient f) freqs
  arrows            = zip startingPositions freqs

arrowsAt :: Seed -> Float -> [Arrow]
arrowsAt seed t = map (\(arr, freq) -> updateArrow arr freq t) (_arrows seed)

drawSeed :: Seed -> Float -> Picture
drawSeed seed seconds = scale s s $ Pictures [vectors, tip, trail]
 where
  s         = _spaceScale seed
  t         = _timeScale seed * seconds
  positions = arrange $ arrowsAt seed t
  (x, y)    = last positions
  -- draw stuff
  vectors   = color (light aquamarine) . line $ positions
  tip       = color red . translate x y $ circleSolid (5 / _spaceScale seed)
  trail     = color white . line $ map (endpoint . arrowsAt seed) ts
  -- the trail draws through points with a step size of 1/1000 of a (time-scaled) second
  tf        = truncateTo 3 t
  ti        = max 0 (tf - _trailLength seed)
  ts        = [ti, ti + 0.001 .. t]

-- truncate a real number to the given number of decimal places
truncateTo :: (RealFrac a) => Integer -> a -> a
truncateTo n = (/ 10 ^ n) . fromIntegral . truncate . (* 10 ^ n)

max_stars_repo_name: juliagracedefoor/fourier-visualizer

max_issues_repo_name: juliagracedefoor/fourier-visualizer

max_forks_repo_name: juliagracedefoor/fourier-visualizer


text:
{-# LANGUAGE ConstraintKinds #-}

module Element
( Element
, StructElem(..)
, getElementNodes
, getNumNodes
, getConnectivity
, getElementNumber
, computeJacobian
, computeJacobianDet
, dndx
) where

import qualified Numeric.LinearAlgebra         as L
import qualified Numeric.LinearAlgebra.HMatrix as HMat

import           Basis
import           Node
import           Quadrature
import           ShapeFcns

-- Utility functions
allCombinations :: [Int] -> [[Int]]
allCombinations x = mapM (const x) [1..(length x)]

-- Synonym for types required by HMatrix
type FracElNum a = (Fractional a, L.Element a, L.Numeric a)

-- Element types
data StructElem a = StructElem [Node a] Int deriving (Show,Eq)
data Simplex a    = Simplex    [Node a] Int deriving (Show,Eq)

-- This typeclass defines which functions operate for any element.  Different
-- element types are defined as instances of this typeclass.
class Element e where
  getElementNodes    :: (Fractional a) => e a -> [Node a]
  getNumNodes        :: (Fractional a) => e a -> Int
  getConnectivity    :: (Fractional a) => e a -> [Int]
  getElementNumber   :: (Fractional a) => e a -> Int
  faceNormals        :: (Fractional a) => e a -> Int -> [a]
  computeJacobian    :: (Basis b,ShapeFcn s,FracElNum a) => e a -> s b -> [a] -> L.Matrix a
  computeJacobianDet :: (Basis b,ShapeFcn s,FracElNum a,L.Field a) => e a -> s b -> [a] -> a
  dndx               :: (Basis b,ShapeFcn s,FracElNum a,L.Field a) => e a -> s b -> [a] -> Int -> L.Vector a

instance Element StructElem where

  getElementNodes  (StructElem nodes _)   = nodes
  getNumNodes      (StructElem nodes _)   = length nodes
  getConnectivity  (StructElem nodes _)   = map nodeNumber nodes
  getElementNumber (StructElem _ elemNum) = elemNum

  computeJacobian  (StructElem nodes elemNum) shpFcn coords = HMat.mul matA matB
    where
      dim  = getDimension shpFcn
      matA = L.tr $ L.fromLists $ map (dndXi shpFcn coords) [0..(getShapeFcnOrder shpFcn + 1)^dim - 1]
      matB = L.fromLists $ map nodeCoordinates nodes

  computeJacobianDet (StructElem nodes elemNum) shpFcn coords = L.det $ computeJacobian (StructElem nodes elemNum) shpFcn coords
  dndx elem shpFcn coords idx = HMat.app (HMat.inv (computeJacobian elem shpFcn coords)) (L.fromList $ dndXi shpFcn coords idx)

  faceNormals (StructElem nodes elemNum) 1 = [-1.0, 1.0]  -- Specialized for linear elements
  faceNormals _ _ = error "Face normals can only be computed for 1D elements for now."


text:
import Numeric.LinearAlgebra
import HVX

n = 4
a = EConst $ (n><n) ([0, 2, 0, 0,
                      -3, 0, 0, 0,
                       0, 0, 1, 0,
                       0, 0, 0, 4] :: [Double])
b = EConst $ (n><n) ([4, 8, -1, 0,
                      0, -2, 1, 0,
                      9, -3, 0, 3,
                      0, 0,  2, 0] :: [Double])
c = EConst $ (n><1) [1..]

x = EVar "x"
y = EVar "y"
constNeg3 = EConst $ scale (-3.0) $ konst 1.0 (n, 1)
constTwos = EConst $ scale 2.0 $ konst 1.0 (n, 1)

ans = subgradMinimize
  (norm 2 (a *~ x) +~ norm 2 (b *~ y) +~ norm 2 (c +~ x +~ neg y))
  [y <=~ constTwos, x <=~ constNeg3]
  (decNonSumStep 10.0) 5000
  [("x", konst 0.0 (n, 1)), ("y", konst 0.0 (n, 1))]

main :: IO ()
main = print ans

max_stars_repo_name: chrisnc/hvx

max_issues_repo_name: chrisnc/hvx

max_forks_repo_name: chrisnc/hvx


text:
#!/usr/bin/env stack
-- stack --resolver lts-5.15 --install-ghc runghc --package hmatrix --package MonadRandom

-- FROM: https://github.com/mstksg/inCode/blob/master/code-samples/dependent-haskell/NetworkUntyped.hs

{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE ScopedTypeVariables #-}

import Control.Monad
import Control.Monad.Random
import Data.List
import Data.Maybe
import Numeric.LinearAlgebra
import System.Environment
import Text.Read

data Weights = W { wBiases :: !(Vector Double)  -- n
                 , wNodes  :: !(Matrix Double)  -- n x m
                 }                              -- "m to n" layer

data Network = Last !Weights
             | Link !Weights Network

logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))

logistic' :: Floating a => a -> a
logistic' x = logix * (1 - logix)
  where
    logix = logistic x

runLayer :: Weights -> Vector Double -> Vector Double
runLayer (W wB wN) v = wB + wN #> v

runNet :: Network -> Vector Double -> Vector Double
runNet (O w)      !v = logistic (runLayer w v)
runNet (w :&~ n') !v = let v' = logistic (runLayer w v)
                       in  runNet n' v'

randomWeights :: MonadRandom m => Int -> Int -> m Weights
randomWeights i o = do
    seed1 :: Int <- getRandom
    seed2 :: Int <- getRandom
    let wB = randomVector  seed1 Uniform o * 2 - 1
        wN = uniformSample seed2 o (replicate i (-1, 1))
    return $ W wB wN

randomNet :: MonadRandom m => Int -> [Int] -> Int -> m Network
randomNet i []     o =     O <$> randomWeights i o
randomNet i (h:hs) o = (:&~) <$> randomWeights i h <*> randomNet h hs o

train :: Double           -- ^ learning rate
      -> Vector Double    -- ^ input vector
      -> Vector Double    -- ^ target vector
      -> Network          -- ^ network to train
      -> Network
train rate x0 target = fst . go x0
  where
    go :: Vector Double    -- ^ input vector
       -> Network          -- ^ network to train
       -> (Network, Vector Double)
    -- handle the output layer
    go !x (O w@(W wB wN))
        = let y    = runLayer w x
              o    = logistic y
              -- the gradient (how much y affects the error)
              --   (logistic' is the derivative of logistic)
              dEdy = logistic' y * (o - target)
              -- new bias weights and node weights
              wB'  = wB - scale rate dEdy
              wN'  = wN - scale rate (dEdy `outer` x)
              w'   = W wB' wN'
              -- bundle of derivatives for next step
              dWs  = tr wN #> dEdy
          in  (O w', dWs)
    -- handle the inner layers
    go !x (w@(W wB wN) :&~ n)
        = let y          = runLayer w x
              o          = logistic y
              -- get dWs', bundle of derivatives from rest of the net
              (n', dWs') = go o n
              -- the gradient (how much y affects the error)
              dEdy       = logistic' y * dWs'
              -- new bias weights and node weights
              wB'  = wB - scale rate dEdy
              wN'  = wN - scale rate (dEdy `outer` x)
              w'   = W wB' wN'
              -- bundle of derivatives for next step
              dWs  = tr wN #> dEdy
          in  (w' :&~ n', dWs)

netTest :: MonadRandom m => Double -> Int -> m String
netTest rate n = do
    inps <- replicateM n $ do
      s <- getRandom
      return $ randomVector s Uniform 2 * 2 - 1
    let outs = flip map inps $ \v ->
                 if v `inCircle` (fromRational 0.33, 0.33)
                      || v `inCircle` (fromRational (-0.33), 0.33)
                   then fromRational 1
                   else fromRational 0
    net0 <- randomNet 2 [16,8] 1
    let trained = foldl' trainEach net0 (zip inps outs)
          where
            trainEach :: Network -> (Vector Double, Vector Double) -> Network
            trainEach nt (i, o) = train rate i o nt

        outMat = [ [ render (norm_2 (runNet trained (vector [x / 25 - 1,y / 10 - 1])))
                   | x <- [0..50] ]
                 | y <- [0..20] ]
        render r | r <= 0.2  = ' '
                 | r <= 0.4  = '.'
                 | r <= 0.6  = '-'
                 | r <= 0.8  = '='
                 | otherwise = '#'

    return $ unlines outMat
  where
    inCircle :: Vector Double -> (Vector Double, Double) -> Bool
    v `inCircle` (o, r) = norm_2 (v - o) <= r

main :: IO ()
main = do
    args <- getArgs
    let n    = readMaybe =<< (args !!? 0)
        rate = readMaybe =<< (args !!? 1)
    putStrLn "Training network..."
    putStrLn =<< evalRandIO (netTest (fromMaybe 0.25   rate)
                                     (fromMaybe 500000 n   )
                            )

(!!?) :: [a] -> Int -> Maybe a
xs !!? i = listToMaybe (drop i xs)

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
import Control.Monad (replicateM)
import System.Random.MWC (create, GenIO)
import System.Random.MWC.Distributions (beta, bernoulli, dirichlet)
import Numeric.LinearAlgebra.HMatrix
       (Matrix, matrix, Vector, vector, konst, fromRows, (!), toList, fromList, size)
import qualified Data.Vector.Storable as V

import Utils (multinomial)
import VB (lQ, calculateEQZ, calculateQZ, calculateEQmu, calculateQmu, expectationBeta, em)

generateSampleMatrix :: Double  -- ^ alpha parameter for Beta distribution
                     -> Double  -- ^ beta parameter for Beta distribution
                     -> Vector Double  -- ^ alpha parameter for Dirichlet distr
                     -> Int  -- ^ number of samples
                     -> Int  -- ^ vector size
                     -> Int  -- ^ Number of mixins
                     -> IO (Matrix Double, Vector Double)
generateSampleMatrix a b alpha n d k = do
    gen <- create
    pivT <- dirichlet (toList alpha) gen
    let piv = fromList pivT
    muK <-
        replicateM
            k
            (replicateM
                 d
                 (beta a b gen) >>=
             \m ->
                  return (vector m))
    let mu = fromRows muK
    genData <-
        replicateM
            n
            (generateVectors piv mu gen)
    return (fromRows genData, piv)
    where generateVectors :: Vector Double
                          -> Matrix Double
                          -> GenIO
                          -> IO (Vector Double)
          generateVectors piv mu gen = do
              zn <- multinomial piv gen
              let Just j = V.findIndex (== 1) zn
              xn <-
                  mapM
                      (\i ->
                            bernoulli
                                ((mu ! j) !
                                 i)
                                gen >>=
                            \b ->
                                 if b
                                     then return 1.0
                                     else return 0.0)
                      [0 .. d - 1]
              return (vector xn)

-- | L(q) function must increase
testLQ = do
  let k = 6
      n = 1000
      d = 3
      a = 0.6
      b = 1
      alpha = konst (1 / 6) k
  (mX, piV) <- generateSampleMatrix a b alpha n d k
  let test_k = 10
      test_alpha = konst 0.1 test_k
      test_a = matrix d [0.5 | _ <- [1..test_k*d]]
      test_b = matrix d [0.5 | _ <- [1..test_k*d]]
  gen <- create
  pivT <- dirichlet (toList test_alpha) gen
  let piv = fromList pivT
  let qZ = fromRows (replicate n piv)
  let eQz = calculateEQZ qZ
  let qMu = (test_a, test_b)
  let eQmu = calculateEQmu qMu
  print (lQ mX eQz eQmu)
  let qZ' = calculateQZ mX eQmu eQz
      eQz' = calculateEQZ qZ'
      qMu' = calculateQmu mX eQz' qMu
      eQmu' = calculateEQmu qMu'
  print (lQ mX eQz' eQmu')
  let qZ'' = calculateQZ mX eQmu' eQz'
      eQz'' = calculateEQZ qZ''
      qMu'' = calculateQmu mX eQz'' qMu'
      eQmu'' = calculateEQmu qMu''
  print (lQ mX eQz'' eQmu'')

-- | EM first step. Can't imagine.
testEM = do
  let k = 6
      n = 1000
      d = 3
      a = 0.6
      b = 1
      alpha = konst (1 / 6) k
  (mX, piV) <- generateSampleMatrix a b alpha n d k
  let test_k = 10
      test_alpha = konst 0.1 test_k
      test_a = 0.5
      test_b = 0.5
  gen <- create
  pivT <- dirichlet (toList test_alpha) gen
  let piv = fromList pivT
  print "Original pi vector"
  print piV
  print "Calculated pi vector"
  print (em mX piv test_alpha test_a test_b)

main :: IO ()
main = do
  putStrLn "Test suite not yet implemented"
  let k = 3
      n = 10
      d = 3
      a = 1
      b = 1
      alpha = konst 0.01 k
  d <- generateSampleMatrix a b alpha n d k
  print d
  putStrLn "Test L(q). Must increase."
  testLQ
  putStrLn "Test EM work or not :)"
  testEM
  return ()

max_stars_repo_name: DbIHbKA/mnistvb

max_issues_repo_name: DbIHbKA/mnistvb

max_forks_repo_name: DbIHbKA/mnistvb


text:
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE TypeOperators       #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
module Test.Grenade.Layers.Internal.Convolution where

import           Grenade.Layers.Internal.Convolution

import           Numeric.LinearAlgebra hiding (uniformSample, konst, (===))

import           Hedgehog
import qualified Hedgehog.Gen as Gen
import qualified Hedgehog.Range as Range

import qualified Test.Grenade.Layers.Internal.Reference as Reference
import           Test.Hedgehog.Compat

prop_im2col_col2im_symmetrical_with_kernel_stride =
  let factors n = [x | x <- [1..n], n `mod` x == 0]
  in  property $ do
        height   <- forAll $ choose 2 100
        width    <- forAll $ choose 2 100
        kernel_h <- forAll $ (height `div`)    <$> Gen.element (factors height)
        kernel_w <- forAll $ (width `div`)     <$> Gen.element (factors width)
        input    <- forAll $ (height >< width) <$> Gen.list (Range.singleton $ height * width) (Gen.realFloat $ Range.linearFracFrom 0 (-100) 100)

        let stride_h = kernel_h
        let stride_w = kernel_w
        let out      = col2im kernel_h kernel_w stride_h stride_w height width . im2col kernel_h kernel_w stride_h stride_w $ input
        input === out

prop_im2col_col2im_behaves_as_reference =
  let ok extent kernel = [stride | stride <- [1..extent], (extent - kernel) `mod` stride == 0]
  in  property $ do
        height   <- forAll (choose 2 100)
        width    <- forAll (choose 2 100)
        kernel_h <- forAll (choose 2 (height - 1))
        kernel_w <- forAll (choose 2 (width - 1))
        stride_h <- forAll (Gen.element (ok height kernel_h))
        stride_w <- forAll (Gen.element (ok width kernel_w))
        input    <- forAll ((height >< width) <$> Gen.list (Range.singleton $ height * width) (Gen.realFloat $ Range.linearFracFrom 0 (-100) 100))

        let outFast       = im2col kernel_h kernel_w stride_h stride_w input
        let retFast       = col2im kernel_h kernel_w stride_h stride_w height width outFast

        let outReference  = Reference.im2col kernel_h kernel_w stride_h stride_w  input
        let retReference  = Reference.col2im kernel_h kernel_w stride_h stride_w height width outReference

        outFast === outReference
        retFast === retReference

tests :: IO Bool
tests = checkParallel $$(discover)

max_issues_repo_name: Alien-Inc/grenade

max_forks_repo_name: Alien-Inc/grenade


text:
module Main where

import Numeric.LinearAlgebra
import Data.Function (on)
import Data.List (minimumBy)

import Lib

type MarcovChain = Vector R -> Vector R

type State = (Int,Int)

type Window = Int

type Steps = Int

type Data t = t State

type States t = t State

classify :: Steps -> [Double] -> [Int]
classify steps ds = map (roundToNearest . (* 100) . (/ maxValue)) ds
  where
    maxValue = maximum ds

    increments = statePercents steps

    roundToNearest :: Double -> Int
    roundToNearest d = case dropWhile (\(_,y) -> fromIntegral y <= d) $ zip increments (tail increments) of
                          ((x,y):_) -> if d - fromIntegral x >= fromIntegral y - d then y else x
                          _         -> last increments

makeMarcovChainFromData :: Foldable t => Window -> Steps -> Data t -> MarcovChain
makeMarcovChainFromData = undefined

stateProduct :: Window -> Steps -> States []
stateProduct winSize steps = [(i,j) | i <- statePercents steps, j <- [1..winSize]]

statePercents :: Steps -> [Int]
statePercents steps = [0,stepSize..200]
  where
    stepSize = round $ 200 / fromIntegral steps

predict :: MarcovChain -> Vector R -> Vector R
predict = ($)

main :: IO ()
main = print "Crypto"

max_stars_repo_name: FelixHolzwarth/MarkovPrediction

max_issues_repo_name: FelixHolzwarth/MarkovPrediction

max_forks_repo_name: FelixHolzwarth/MarkovPrediction


text:
{-# OPTIONS_GHC -Wall #-}

module Numeric.FFT.HMatrix.Matrix (
  -- * Complex-to-complex transforms
  dft
  , idft
    -- * Real-to-complex transforms
  , dftR2C
  , dftC2R
    -- * Real-to-real transforms
    -- ** Discrete cosine transforms
  , dct1
  , dct2
  , dct3
  , dct4
    -- ** Discrete sine transforms
  , dst1
  , dst2
  , dst3
  , dst4
  ) where

import Foreign.Storable ( Storable )
import Numeric.LinearAlgebra.Data ( Matrix, Complex )
import Numeric.LinearAlgebra ( Element )
import qualified Numeric.LinearAlgebra.Data as M
import Numeric.FFT.Vector.Unnormalized ( Transform )
import qualified Numeric.FFT.Vector.Unnormalized as FFT

transResize :: (Storable a, Element b) =>
               (Matrix b -> Int)
            -> ((Int, Int) -> Transform b a)
            -> Matrix b
            -> Matrix a
transResize calcCols f m = M.reshape newcols . FFT.run (f (r, c)) $ M.flatten m
  where
    c = M.cols m
    r = M.rows m
    newcols = calcCols m

trans :: (Storable a, Element b) =>
         ((Int, Int) -> Transform b a)
      -> Matrix b
      -> Matrix a
trans = transResize (M.cols)

-- |
-- >>> M.toLists . dft $ M.fromLists [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
-- [[45.0 :+ 0.0,(-4.5) :+ 2.598076211353316,(-4.5) :+ (-2.598076211353316)],[(-13.5) :+ 7.794228634059947,0.0 :+ 0.0,0.0 :+ 0.0],[(-13.5) :+ (-7.794228634059947),0.0 :+ 0.0,0.0 :+ 0.0]]
dft :: Matrix (Complex Double) -> Matrix (Complex Double)
dft = trans FFT.dft2d

-- |
-- >>> M.toLists . idft $ M.fromLists [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
-- [[45.0 :+ 0.0,(-4.5) :+ (-2.598076211353316),(-4.5) :+ 2.598076211353316],[(-13.5) :+ (-7.794228634059947),0.0 :+ 0.0,0.0 :+ 0.0],[(-13.5) :+ 7.794228634059947,0.0 :+ 0.0,0.0 :+ 0.0]]
idft :: Matrix (Complex Double) -> Matrix (Complex Double)
idft = trans FFT.idft2d

-- |
-- >>> M.toLists . dftR2C $ M.fromLists [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
-- [[45.0 :+ 0.0,(-4.5) :+ 2.598076211353316],[(-13.5) :+ 7.794228634059947,0.0 :+ 0.0],[(-13.5) :+ (-7.794228634059947),0.0 :+ 0.0]]
dftR2C :: Matrix Double -> Matrix (Complex Double)
dftR2C = transResize (\m -> M.cols m `div` 2 + 1) FFT.dftR2C2D

-- |
--
-- >>> let m = M.fromLists [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
-- >>> let mtest = M.fromLists [[33,-3,-3,-9],[3,-3,-9,-3],[3,0,0,0]]
-- >>> M.maxElement (mtest - dftC2R m) < 1e-12
-- True
dftC2R :: Matrix (Complex Double) -> Matrix Double
dftC2R = transResize (\m -> 2 * (M.cols m - 1)) FFT.dftC2R2D

dct1 :: Matrix Double -> Matrix Double
dct1 = trans FFT.dct1_2D

dct2 :: Matrix Double -> Matrix Double
dct2 = trans FFT.dct2_2D

dct3 :: Matrix Double -> Matrix Double
dct3 = trans FFT.dct3_2D

dct4 :: Matrix Double -> Matrix Double
dct4 = trans FFT.dct4_2D

dst1 :: Matrix Double -> Matrix Double
dst1 = trans FFT.dst1_2D

dst2 :: Matrix Double -> Matrix Double
dst2 = trans FFT.dst2_2D

dst3 :: Matrix Double -> Matrix Double
dst3 = trans FFT.dst3_2D

dst4 :: Matrix Double -> Matrix Double
dst4 = trans FFT.dst4_2D

max_stars_repo_name: peddie/hmatrix-fftw

max_issues_repo_name: peddie/hmatrix-fftw

max_forks_repo_name: peddie/hmatrix-fftw


text:
import Numeric.LinearAlgebra
import Graphics.Plot
import Numeric.GSL.Special(erf_Z, erf)

sombrero n = f x y where
    (x,y) = meshdom range range
    range = linspace n (-2,2)
    f x y = exp (-r2) * cos (2*r2) where
        r2 = x*x+y*y

f x =  sin x + 0.5 * sin (5*x)

gaussianPDF = erf_Z
cumdist x = 0.5 * (1+ erf (x/sqrt 2))

main = do
    let x = linspace 1000 (-4,4)
    mplot [f x]
    mplot [x, mapVector cumdist x, mapVector gaussianPDF x]
    mesh (sombrero 40)


text:
import Data.Complex
import System.Environment
import System.Process

import Lib

func c z = z^2 + c

phi = (1 + sqrt 5) / 2

main = do
    args <- getArgs

    let real:imag:_ = if null args then [1 - phi, 0] else read <$> args

        points = (:+) <$> [-1.5,-1.495..1.5] <*>  [-1.2,-1.195..1.2]
        jset = c2t <$> filter (not . diverges 50 (func (real :+ imag))) points

        fileName = "julia" ++ show real ++ "-" ++ show imag ++ ".svg"
        title = "Julia Set for z^2 + " ++ show real ++ " + " ++ show imag ++ "i"

    plotPoints fileName title 1.7 jset
    callCommand $ "inkview " ++ fileName

max_stars_repo_name: ahaym/dynamics

max_issues_repo_name: ahaym/dynamics

max_forks_repo_name: ahaym/dynamics


text:
module MachineLearningTest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Plus

import MachineLearning

import qualified Numeric.LinearAlgebra as LA

x = LA.matrix 3 [17, 19, 29]
x2 = LA.matrix 9 [17, 19, 29, 289, 323, 493, 361, 551, 841]
x3 = LA.matrix 19 [17, 19, 29, 289, 323, 493, 361, 551, 841, 4913, 5491, 8381, 6137, 9367, 14297, 6859, 10469, 15979, 24389]
eps = 1e-5

tests = [testGroup "mapFeatures" [
            testCase "mapfeatures 1" $ assertMatrix "" eps x (mapFeatures 1 x)
            , testCase "mapfeatures 2" $ assertMatrix "" eps x2 (mapFeatures 2 x)
            , testCase "mapfeatures 3" $ assertMatrix "" eps x3 (mapFeatures 3 x)
            ]
        ]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE Strict #-}
module STC.Multiplication where

import           Control.Concurrent.Async
import           Control.Monad            as M
import           Data.Array.Repa          as R
import           Data.Complex
import           Data.List                as L
import           Data.Vector.Storable     as VS
import           DFT.Plan
import           Filter.Utils
import           STC.DFTArray
import           Utils.List
import           Utils.Parallel

{-# INLINE multiplyPinwheelBasis #-}
multiplyPinwheelBasis :: DFTPlan -> DFTArray -> DFTArray -> IO DFTArray
multiplyPinwheelBasis plan (DFTArray rows cols thetaFreqs rFreqs vecs1) (DFTArray _ _ _ _ vecs2) = do
  dftVecs1 <-
    dftExecuteBatchP plan (DFTPlanID DFT1DG [cols, rows] [0, 1]) .
    L.map
      (VS.convert .
       toUnboxed .
       computeS . makeFilter2D . fromUnboxed (Z :. cols :. rows) . VS.convert) $
    vecs1
  dftVecs2 <- dftExecuteBatchP plan (DFTPlanID DFT1DG [cols, rows] [0, 1]) vecs2
  let dftVecs3 = L.zipWith (VS.zipWith (*)) dftVecs1 dftVecs2
  DFTArray rows cols thetaFreqs rFreqs <$>
    dftExecuteBatchP plan (DFTPlanID IDFT1DG [cols, rows] [0, 1]) dftVecs3

{-# INLINE multiplyPinwheelBasisBatch #-}
multiplyPinwheelBasisBatch ::
     DFTPlan -> Int -> DFTArray -> DFTArray -> IO DFTArray
multiplyPinwheelBasisBatch plan numBatch (DFTArray rows cols thetaFreqs rFreqs vecs1) (DFTArray _ _ _ _ vecs2) = do
  let forwardPlanID = DFTPlanID DFT1DG [cols, rows] [0, 1]
      backwardPlanID = DFTPlanID IDFT1DG [cols, rows] [0, 1]
      forwardPlan = getDFTPlan plan forwardPlanID
      backwardPlan = getDFTPlan plan backwardPlanID
  fmap (DFTArray rows cols thetaFreqs rFreqs . L.concat) .
    mapConcurrently
      (M.mapM
         (\(vec1, vec2) -> do
            dftVec1 <-
              dftExecuteWithPlan forwardPlanID forwardPlan .
              VS.convert .
              toUnboxed .
              computeS .
              makeFilter2D . fromUnboxed (Z :. cols :. rows) . VS.convert $
              vec1
            dftVec2 <- dftExecuteWithPlan forwardPlanID forwardPlan vec2
            dftExecuteWithPlan backwardPlanID backwardPlan $
              VS.zipWith (*) dftVec1 dftVec2)) .
    divideListN numBatch $
    L.zip vecs1 vecs2

{-# INLINE multiplyPinwheelBasisBatch1 #-}
multiplyPinwheelBasisBatch1 ::
     DFTPlan -> Int -> VS.Vector (Complex Double) -> DFTArray -> IO DFTArray
multiplyPinwheelBasisBatch1 plan numBatch dftVec1 (DFTArray rows cols thetaFreqs rFreqs vecs2) = do
  let forwardPlanID = DFTPlanID DFT1DG [cols, rows] [0, 1]
      backwardPlanID = DFTPlanID IDFT1DG [cols, rows] [0, 1]
      forwardPlan = getDFTPlan plan forwardPlanID
      backwardPlan = getDFTPlan plan backwardPlanID
  fmap (DFTArray rows cols thetaFreqs rFreqs . L.concat) .
    mapConcurrently
      (M.mapM
         (\vec2 -> do
            dftVec2 <- dftExecuteWithPlan forwardPlanID forwardPlan vec2
            dftExecuteWithPlan backwardPlanID backwardPlan $
              VS.zipWith (*) dftVec1 dftVec2)) .
    divideListN numBatch $
    vecs2

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
module Types (
  LispVal(..),
  LispError(..),
  ThrowsError,
  trapError,
  extractValue,
  Env,
  nullEnv,
  IOThrowsError,
  liftThrows,
  runIOThrows,
  showVal
) where

import           Control.Monad.Error
import           Data.Complex
import           Data.IORef
import           System.IO
import           Text.ParserCombinators.Parsec

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

showVal :: LispVal -> String
showVal (String contents) = "\"" ++ contents ++ "\""
showVal (Atom name)       = name
showVal (Number contents) = show contents
showVal (Bool True)       = "#t"
showVal (Bool False)      = "#f"
showVal (List contents)   = "(" ++ unwordsList contents ++ ")"
showVal (DottedList head tail) = "(" ++ unwordsList head ++ " . " ++ showVal tail ++ ")"
showVal (Character content) = case content of
  ' '  -> "#\\space"
  '\n' -> "#\\newline"
  _    -> "#\\" ++ [content]
showVal (Float x) = show x
showVal (Ratio x) = show x
showVal (Complex x) = show x
showVal (PrimitiveFunc _) = "<primitive>"
showVal Func { params = args, vararg = varargs, body = body, closure = env } =
  "(lambda (" ++ unwords (map show args) ++ (
    case varargs of
      Nothing  -> ""
      Just arg -> " . " ++arg
  ) ++ ") ...)"
showVal (Port _) = "<IO port>"
showVal (IOFunc _) = "<IO primitive>"

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

instance Show LispVal where show = showVal
instance Eq LispVal where
  (==) a b = show a == show b

data LispError = NumArgs Integer [LispVal]
               | TypeMismatch String LispVal
               | Parser ParseError
               | BadSpecialForm String LispVal
               | NotFunction String String
               | UnboundVar String String
               | Default String

showError :: LispError -> String
showError (UnboundVar message varname) = message ++ ": " ++ varname
showError (BadSpecialForm message form) = message ++ ": " ++ show form
showError (NotFunction message func) = message ++ ": " ++ show func
showError (NumArgs expected found) =
  "Expected " ++ show expected
  ++ " args; found values " ++ unwordsList found
showError (TypeMismatch expected found) =
  "Invalid type: expected " ++ expected
  ++ ", found " ++ show found
showError (Parser parseErr) = "Parse error at " ++ show parseErr
showError (Default err) = "LispError: " ++ show err

instance Show LispError where show = showError
instance Error LispError where
  noMsg = Default "An error has occurred"
  strMsg = Default

type ThrowsError = Either LispError

trapError action = catchError action (return . show)

extractValue :: ThrowsError a -> a
extractValue (Right val) = val

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

nullEnv :: IO Env
nullEnv = newIORef []

type IOThrowsError = ErrorT LispError IO

liftThrows :: ThrowsError a -> IOThrowsError a
liftThrows (Left err)  = throwError err
liftThrows (Right val) = return val

runIOThrows :: IOThrowsError String -> IO String
runIOThrows action = fmap extractValue (runErrorT (trapError action))


text:
module Data.SQP
       ( Problem (..)
       , optimize
       ) where

import Numeric.LinearAlgebra.HMatrix
import Numeric.Minimization.QuadProgPP

import Debug.Trace

data Problem = Problem
               { _cost :: Vector Double -> Double
               , _approxCost :: Vector Double
                             -> (Matrix Double, Vector Double, Double)
               , _trueIneqs :: Vector Double -> Vector Double
               , _approxAffineIneqs :: Vector Double
                                    -> (Matrix Double, Vector Double)
               , _eqs :: Vector Double -> Vector Double
               , _approxEqs :: Vector Double -> (Matrix Double, Vector Double)
               , _numVariables :: Int
               , _numIneqs :: Int
               , _numEqs :: Int
               }

stepAcceptanceThreshold :: Double
stepAcceptanceThreshold = 0.25

trustShrinkFactor :: Double
trustShrinkFactor = 0.1

trustExpandFactor :: Double
trustExpandFactor = 1.5

constraintPenaltyScalingFactor :: Double
constraintPenaltyScalingFactor = 10

minTrustSize :: Double
minTrustSize = 1e-4

minModelImprove :: Double
minModelImprove = 1e-4

minModelImproveRatio :: Double
minModelImproveRatio = negate $ read "Infinity"

constraintSatisfactionThreshold :: Double
constraintSatisfactionThreshold = 1e-4

initTrustSize :: Double
initTrustSize = 0.1

initPenalty :: Double
initPenalty = 10

evalMerit :: Problem -> Double -> Vector Double -> Double
evalMerit problem penaltyParam x =
  let cost = _cost problem x
      penalty =
        penaltyParam * sum (map (max 0.0) $ toList $ _trueIneqs problem x) +
        penaltyParam * sum (map abs $ toList $ _eqs problem x)
  in  cost + penalty

optimize :: Problem
         -> Vector Double -- xInitial
         -> (Vector Double, Double)
optimize problem xInit =
  findSuitableConstraintPenalty problem xInit initTrustSize initPenalty

findSuitableConstraintPenalty :: Problem
                              -> Vector Double -- x
                              -> Double -- trust region size
                              -> Double -- constraint penalty
                              -- xNew, trueMerit
                              -> (Vector Double, Double)
findSuitableConstraintPenalty problem x trustSize penaltyParam =
  let (xNew, trueMerit, newTrustSize) =
        reconvexify problem x trustSize penaltyParam
      ineqsSatisfied = all (<= constraintSatisfactionThreshold) $
                         toList $ _trueIneqs problem xNew
      eqsSatisfied = all ((<= constraintSatisfactionThreshold) . abs) $
                       toList $ _eqs problem xNew
      constraintsSatisfied = ineqsSatisfied && eqsSatisfied
  in  if constraintsSatisfied
      then (xNew, trueMerit)
      else let penalty' = penaltyParam * constraintPenaltyScalingFactor
           in  findSuitableConstraintPenalty
                 problem xNew newTrustSize penalty'

evalQuadratic :: (Matrix Double, Vector Double, Double)
              -> Vector Double
              -> Double
evalQuadratic (m, b, c) x =
  (x `dot` ((0.5 * (m #> x)) + b)) + c

evalLinear :: (Matrix Double, Vector Double)
           -> Vector Double
           -> Vector Double
evalLinear (m, b) x = (m #> x) + b

reconvexify :: Problem
            -> Vector Double -- x
            -> Double -- trust region size
            -> Double -- constraint penalty
            -- xNew, newTrueMerit, newTrustSize
            -> (Vector Double, Double, Double)
reconvexify problem x trustSize penaltyParam =
  let convexCost = _approxCost problem x
      convexIneqs = _approxAffineIneqs problem x
      convexEqs = _approxEqs problem x
      trueMerit = evalMerit problem penaltyParam x
      trustResult = findSuitableTrustStep problem x convexCost
                      convexIneqs convexEqs trueMerit trustSize
                      penaltyParam
  in  case trustResult of
        Reconvexify xNew newMerit newTrustSize ->
          reconvexify problem xNew newTrustSize penaltyParam
        Finished xNew newMerit newTrustSize ->
          (xNew, newMerit, newTrustSize)

data TrustStepResult = -- x merit trustSize
                       Reconvexify (Vector Double) Double Double
                     | Finished (Vector Double) Double Double
  deriving (Show)

findSuitableTrustStep :: Problem
                      -> Vector Double -- x
                      -- convexified cost
                      -> (Matrix Double, Vector Double, Double)
                      -- convexified inequality constraints
                      -> (Matrix Double, Vector Double)
                      -- convexified equality constraints
                      -> (Matrix Double, Vector Double)
                      -> Double -- old true merit
                      -> Double -- trust region size
                      -> Double -- constraint penalty parameter
                      -- xNew, newTrueMerit, newTrustSize
                      -> TrustStepResult
findSuitableTrustStep
  problem x convexCost convexIneq convexEq oldTrueMerit trustSize penaltyParam =
  let trustStep =
        trustRegionStep problem x convexCost convexIneq convexEq oldTrueMerit
          trustSize penaltyParam
  in  case trace (show trustStep) trustStep of
        Reject -> let newTrustSize = trustSize * trustShrinkFactor
                  in  findSuitableTrustStep problem x convexCost convexIneq
                        convexEq oldTrueMerit newTrustSize penaltyParam
        Accept xNew newTrueMerit ->
          let newTrustSize = trustSize * trustExpandFactor
          in  Reconvexify xNew newTrueMerit newTrustSize
        Converged xNew newTrueMerit ->
          let newTrustSize =
                max trustSize $ (minTrustSize / trustShrinkFactor) * 1.5
          in  Finished xNew newTrueMerit newTrustSize

data IterationResult = Reject
                     -- xNew newMerit
                     | Accept (Vector Double) Double
                     | Converged (Vector Double) Double
  deriving (Show)

trustRegionStep :: Problem
                -> Vector Double -- x
                -- convexified cost
                -> (Matrix Double, Vector Double, Double)
                -- convexified inequality constraints
                -> (Matrix Double, Vector Double)
                -- convexified equality constraints
                -> (Matrix Double, Vector Double)
                -> Double -- old true merit
                -> Double -- trust region size
                -> Double -- constraint penalty parameter
                -> IterationResult
trustRegionStep
  problem x convexCost convexIneq convexEq oldTrueMerit trustSize penaltyParam =
    let (xNew, modelMerit) = solveQuadraticSubproblem problem x convexCost
                               convexIneq convexEq trustSize penaltyParam
        trueMerit =  evalMerit problem penaltyParam xNew
        trueImprove = oldTrueMerit - trueMerit
        modelImprove = oldTrueMerit - modelMerit
    in  if modelImprove < -1e-2 -- why so high? our convexification?
        then error $ "Model improvement got worse: " ++
               show (modelImprove, modelMerit, oldTrueMerit, trueMerit)
        else
          if trustSize < minTrustSize ||
             modelImprove < minModelImprove ||
             modelImprove / oldTrueMerit < minModelImproveRatio
          then Converged x oldTrueMerit
          else
            if trueImprove > 0.0 &&
               trueImprove / modelImprove > stepAcceptanceThreshold
            then Accept xNew trueMerit
            else Reject

solveQuadraticSubproblem :: Problem
                         -> Vector Double
                         -- convex cost
                         -> (Matrix Double, Vector Double, Double)
                         -- convexified inequality constraints
                         -> (Matrix Double, Vector Double)
                         -- convexified equality constraints
                         -> (Matrix Double, Vector Double)
                         -> Double -- trust region size
                         -> Double -- constraint penalty parameter
                         -> (Vector Double, Double) -- new x, model merit
solveQuadraticSubproblem
  problem x (costMatrix, costVector, costConstant)
  (approxIneqMatrix, approxIneqVector) (approxEqMatrix, approxEqVector)
  trustSize penaltyParam =
  -- We approximate each nonlinear inequality as |ax + b|^+. For each
  -- of these, we introduce a new optimization variable t (i.e., a
  -- slack variable) that comes with two inequalities:
  --
  -- 0 <= t
  --
  -- ax + b <= t
  --
  -- Each nonlinear equality as |ax + b|. For each of these, two new
  -- optimization variables s and t with following properties:
  --
  -- -s <= 0
  -- -t <= 0
  -- s - t = ax + b
  let numIneqs = _numIneqs problem
      numEqs = _numEqs problem
      numVariables = _numVariables problem

      costMatrixWithSlacks =
        diagBlock [ costMatrix
                  , konst 0.0 (numIneqs, numIneqs)
                  , konst 0.0 (2 * numEqs, 2 * numEqs)]

      -- LET'S MAKE COST SPD
      (l, v) = eigSH costMatrixWithSlacks
      positiveEigVals = cmap (max 1e-12) l
      augmentedCostMatrix = v <> diag positiveEigVals <> tr v

      augmentedCostVector = vjoin [costVector
                                  , konst 1.0 numIneqs
                                  , konst 1.0 $ 2 * numEqs]

      augmentedIneqMatrix =
        ((scalar penaltyParam * approxIneqMatrix) |||
         negate (ident numIneqs) |||
         konst 0.0 (numIneqs, 2 * numEqs)) ===
         (konst 0.0 (numIneqs + 2 * numEqs, numVariables) |||
          negate (ident $ numIneqs + 2 * numEqs))
      augmentedIneqVector =
        vjoin [ scalar penaltyParam * approxIneqVector
              , konst 0.0 $ numIneqs + 2 * numEqs ]

      ineqMatrixWithTrustConstraints =
        augmentedIneqMatrix ===
        (ident numVariables ||| konst 0.0 (numVariables, numIneqs + 2 * numEqs)) ===
        (negate (ident numVariables) ||| konst 0.0 (numVariables, numIneqs + 2 * numEqs))
      ineqVectorWithTrustConstraints =
        vjoin [ augmentedIneqVector
              , negate (x + scalar trustSize)
              , x - scalar trustSize ]

      augmentedEqMatrix =
        scalar penaltyParam * approxEqMatrix |||
        konst 0.0 (numEqs, numIneqs) |||
        negate (ident numEqs) |||
        ident numEqs
      augmentedEqVector = scalar penaltyParam * approxEqVector

      -- inequalities are given as ax + b <= 0, but quadprog++ wants
      -- a'x + b' >= 0
      result = solveQuadProg
        (augmentedCostMatrix, augmentedCostVector)
        (Just (augmentedEqMatrix, augmentedEqVector)) $
        Just ((-ineqMatrixWithTrustConstraints),
              (-ineqVectorWithTrustConstraints))

  in  case result of
        Left e -> error $ show e
        Right (xNewAugmented, newMeritHomogenous) ->
          let xNew = subVector 0 numVariables xNewAugmented
              newMerit = newMeritHomogenous + costConstant

              -- Compute the merit of the approximation without regularization
              cost = evalQuadratic (costMatrix, costVector, costConstant) xNew
              ineqViolations =
                evalLinear (approxIneqMatrix, approxIneqVector) xNew
              ineqPenalty =
                penaltyParam * (sum $ map (max 0.0) $ toList ineqViolations)
              eqViolations = evalLinear (approxEqMatrix, approxEqVector) xNew
              eqPenalty =
                penaltyParam * (sum $ map abs $ toList eqViolations)
              meritNoReg = cost + ineqPenalty + eqPenalty

              -- DEBUG
              meritError = abs $ newMerit - meritNoReg
          in  if meritError / (abs meritNoReg) > 0.5
              then error $ "merits don't match!" ++ show (meritNoReg, newMerit)
              else (xNew, meritNoReg)

max_stars_repo_name: giogadi/hs-trajopt

max_issues_repo_name: giogadi/hs-trajopt

max_forks_repo_name: giogadi/hs-trajopt


text:
{-# LANGUAGE OverloadedStrings #-}

module Main where

import qualified    Data.Text.IO as T
import qualified    Data.Text as T
import              Data.Text (Text)
import              Data.Text.Encoding (encodeUtf8)

import qualified    Data.ByteString.Lazy as B
import              Data.ByteString.Lazy (ByteString)
import qualified    Data.Map as Map
import              Data.Map (Map)
import qualified    Data.Vector as V
import              Data.Vector (Vector)
import              Data.Csv
import              Data.Monoid

import              Parser
import              Type

import              Control.Monad
import              Control.Monad.Primitive
import              System.Random.MWC
import              Statistics.Distribution
import              Statistics.Distribution.Laplace

type County a = Map Text a
type District a = Map Text (County a)
type FlatDistrict a = Map Text a

iteration :: Int
iteration = 1000

main :: IO ()
main = do
    raw <- B.readFile "asset/sim_2.csv"
    case parse raw of
        Left err -> putStrLn err
        Right rows -> do
            let nation = buildUp rows
            calculateNation nation
            calculateCounty nation
            calculateDistrict nation

calculateNation :: District (Vector Row) -> IO ()
calculateNation nation = do
    putStrLn "全國平均就醫距離"
    let header = V.fromList [encodeUtf8 "全國"] :: Vector Name
    let result = [NationResult (nationAvg nation)]
    B.writeFile "result/original/nationwide.csv" (encodeByName header result)

    putStrLn "全國平均就醫距離 (added noise)"
    forM_ [-9 .. 5] $ \level -> do
        results <- replicateM iteration (nationAddNoise level nation)
        B.writeFile ("result/pertubated/nationwide/2^" <> show level <> ".csv") (encodeByName header results)
--
calculateCounty :: District (Vector Row) -> IO ()
calculateCounty nation = do
    putStrLn "各縣市平均就醫距離"
    let header = V.fromList (map encodeUtf8 (Map.keys nation)) :: Vector Name
    let avgs = fmap countyAvg nation
    let result = [CountyResult avgs]
    B.writeFile "result/original/countywide.csv" (encodeByName header result)

    withSystemRandom . asGenIO $ \gen -> do
        putStrLn "各縣市平均就醫距離 (added noise)"
        forM_ [-4 .. 10] $ \level -> do
            results <- replicateM iteration (countyAddNoise level nation gen)
            B.writeFile ("result/pertubated/countywide/2^" <> show level <> ".csv") (encodeByName header results)

calculateDistrict :: District (Vector Row) -> IO ()
calculateDistrict nation = do
    let flattened = flatten nation
    putStrLn "各鄉鎮市區的平均就醫距離"
    let header = V.fromList $ map encodeUtf8 $ Map.keys flattened
    let result = [DistrictResult (flatten $ fmap (fmap districtAvg) nation)]
    B.writeFile "result/original/districtwide.csv" (encodeByName header result)

    putStrLn "各鄉鎮市區的平均就醫距離 (added noise)"

    withSystemRandom . asGenIO $ \gen -> do
        forM_ [5 .. 20] $ \level -> do
            results <- replicateM iteration (districtAddNoise level flattened gen)
            B.writeFile ("result/pertubated/districtwide/2^" <> show level <> ".csv") (encodeByName header results)

nationAddNoise :: Int -> District (Vector Row) -> IO NationResult
nationAddNoise level nation = do
    let avg = nationAvg nation
    pertubated <- addNoise scale avg
    return (NationResult pertubated)
    where
        size = Map.foldl (\s elem -> s + (Map.foldl (\s vec -> s + V.length vec) 0 elem)) 0 nation
        eps = 2.0 ^^ level
        scale = 400000.0 / (fromIntegral size * eps)
            -- pertubated <- addNoise' scale (replicate iteration (nationwideAvg mapping))

countyAddNoise :: Int -> District (Vector Row) -> Gen (PrimState IO) -> IO CountyResult
countyAddNoise level nation gen = do
    pertubated <- mapM (\county -> do
        let avg = countyAvg county
        noise <- genContVar (laplace 0 (scale county)) gen
        return (avg + noise)) nation
    return (CountyResult pertubated)

    where
        countySize = Map.foldl (\s vec -> s + V.length vec) 0
        eps = 2.0 ^^ level
        scale county = 400000.0 / (fromIntegral (countySize county) * eps)

districtAddNoise :: Int -> FlatDistrict (Vector Row) -> Gen (PrimState IO) -> IO DistrictResult
districtAddNoise level flattened gen = do
    pertubated <- mapM (\district -> do
        let avg = districtAvg district
        noise <- genContVar (laplace 0 (scale district)) gen
        return (avg + noise)) flattened
    return (DistrictResult pertubated)

    where
        districtSize = V.length -- Map.foldl (\s vec -> s + V.length vec) 0
        eps = 2.0 ^^ level
        scale district = 400000.0 / (fromIntegral (districtSize district) * eps)

data NationResult = NationResult Double deriving (Show)
data CountyResult = CountyResult (County Double) deriving (Show)
data DistrictResult = DistrictResult (FlatDistrict Double) deriving (Show)

instance ToNamedRecord NationResult where
    toNamedRecord (NationResult value) = namedRecord
        [ (encodeUtf8 "全國") .= value ]

instance ToNamedRecord CountyResult where
    toNamedRecord (CountyResult value)
        = namedRecord (map (\(key, val) -> (encodeUtf8 key) .= val) (Map.assocs value))

flatten :: District a -> FlatDistrict a
flatten = Map.fromList . concat . map fuseName . Map.assocs
    where
        fuseName :: (Text, County a) -> [(Text, a)]
        fuseName (countyName, county) = map (\(districtName, distrct) -> (countyName <> districtName, distrct)) (Map.assocs county)

instance ToNamedRecord DistrictResult where
    toNamedRecord (DistrictResult nation)
        = namedRecord
            $ Map.elems
            $ Map.mapWithKey (\name value -> encodeUtf8 name .= value)
            $ nation

districtAvg :: Vector Row -> Double
districtAvg = average . V.toList . fmap distance

countyAvg :: County (Vector Row) -> Double
countyAvg = average . concat . Map.elems . fmap (V.toList . fmap distance)
    -- average . Map.elems . fmap districtAvg
    -- Map.foldl (\ acc distrct -> acc + districtAvg distrct) 0.0

nationAvg :: District (Vector Row) -> Double
nationAvg = average . concat . Map.elems . fmap (concat . Map.elems . fmap (V.toList . fmap distance))

buildUp :: Vector Row -> District (Vector Row)
buildUp = foldl insertCounty Map.empty
    where

        districtSingleton :: Row -> County (Vector Row)
        districtSingleton row = Map.singleton (snd (key row)) (V.singleton row)

        insertCounty :: District (Vector Row) -> Row -> District (Vector Row)
        insertCounty m row = Map.insertWith
            (insertDistrict row)
            (fst (key row))         -- key for TierCounty
            (districtSingleton row)    -- key for TierCounty (which is TierDistrict)
            m

        insertDistrict :: Row -> County (Vector Row) -> County (Vector Row) -> County (Vector Row)
        insertDistrict row _ old = Map.insertWith
            (V.++)
            (snd (key row))         -- key for TierDistrict
            (V.singleton row)       -- value for TierDistrict
            old

        key :: Row -> (Text, Text)
        key = T.splitAt 3 . townshipName . township . from

        townshipName :: Township -> Text
        townshipName (Township name _) = name

--------------------------------------------------------------------------------
-- Queries
--------------------------------------------------------------------------------

average :: [Double] -> Double
average vec = sum vec / fromIntegral (length vec)

distance :: Row -> Double
distance row = sqrt ((fromX - toX) * (fromX - toX) + (fromY - toY) * (fromY - toY))
    where
        (fromX, fromY) = coord (from row)
        (toX  , toY  ) = coord (to row)

--------------------------------------------------------------------------------
-- Noise
--------------------------------------------------------------------------------

sampleLaplace :: Double -> IO Double
sampleLaplace scale = withSystemRandom . asGenST $ genContVar (laplace 0 scale)

addNoise :: Double -> Double -> IO Double
addNoise scale x = do
    noise <- sampleLaplace scale
    return (x + noise)

-- -- 研究問題
--
-- -- 1. 各鄉鎮市區的平均就醫距離
-- -- 2. 各縣市的平均就醫距離
-- -- 3. 全國的平均就醫距離
--
-- -- 距離上限: 400 公里

max_stars_repo_name: banacorn/dp-trial

max_issues_repo_name: banacorn/dp-trial

max_forks_repo_name: banacorn/dp-trial


text:
-- DTMF Analyzer
-- Haskell School Project (Faculty of Mathematics and Physics, Charles University in Prague)
-- Author: Lukas Jendele
import System.Environment
import qualified Data.ByteString.Lazy as B
import qualified Data.Complex as C
import Data.Binary.Get
import Data.Word
import qualified Numeric as Numeric
import qualified Data.List as List
import qualified Data.Set as Set
import qualified Data.Map as Map
import qualified Data.Function as DataFunction

nthPrimitiveSquareRoot :: Int -> C.Complex Float
nthPrimitiveSquareRoot 0 = 1
nthPrimitiveSquareRoot n = C.conjugate $ C.mkPolar 1 (2*pi/(fromIntegral n))

-- returns filters odd-indexed elements
takeOdd :: [C.Complex a] -> [C.Complex a]
takeOdd lst = map snd $ filter (\(x,_) -> x `mod` 2 == 1) $ zip [0..] lst
-- returns filters even-indexed elements
takeEven :: [C.Complex a] -> [C.Complex a]
takeEven lst = map snd $ filter (\(x,_) -> x `mod` 2 == 0) $ zip [0..] lst

--FFT see: http://mj.ucw.cz/vyuka/ads/44-fft.pdf
fftHelper1 :: Int -> Int -> C.Complex Float -> [C.Complex Float] -> [C.Complex Float] -> [C.Complex Float]
fftHelper1 _ _ _ [] [] = []
fftHelper1 i n omega (e:evenLst) (o:oddLst) | i == n    = []
                                            | otherwise = (e + ((omega^i) * o)):(fftHelper1 (i+1) n omega evenLst oddLst)
--FFT see: http://mj.ucw.cz/vyuka/ads/44-fft.pdf
fftHelper2 :: Int -> Int -> C.Complex Float -> [C.Complex Float] -> [C.Complex Float] -> [C.Complex Float]
fftHelper2 _ _ _ [] [] = []
fftHelper2 i n omega (e:evenLst) (o:oddLst) | i == n    = []
                                            | otherwise = (e - ((omega^i) * o)):(fftHelper2 (i+1) n omega evenLst oddLst)

--FFT see: http://mj.ucw.cz/vyuka/ads/44-fft.pdf
--fft :: (Integral a) => a -> C.Complex Float -> [C.Complex Float] -> [C.Complex Float]
fft :: Int -> C.Complex Float -> [C.Complex Float] -> [C.Complex Float]
fft n omega lst = if n == 1 then (take 1 lst) else
                    let evenLst = (fft (n `div` 2) (omega^2) (takeEven lst))
                        oddLst  = (fft (n `div` 2) (omega^2) (takeOdd lst))
                    in ((fftHelper1 0 (n `div` 2) omega evenLst oddLst) ++ (fftHelper2 0 (n `div` 2) omega evenLst oddLst))

--fftForward :: (Integral a, RealFloat a) => a -> [C.Complex Float] -> [C.Complex Float]
fftForward :: Int -> [C.Complex Float] -> [C.Complex Float]
fftForward n lst = fft n (nthPrimitiveSquareRoot n) lst

-- divides complex number by components
div' :: C.Complex Float -> Int -> C.Complex Float
div' c r = ((C.realPart c) / (fromIntegral r)) C.:+ ((C.imagPart c) / (fromIntegral r))

--fftBackward :: (Integral a, RealFloat a) => a -> [C.Complex a] -> [C.Complex a]
fftBackward :: Int -> [C.Complex Float] -> [C.Complex Float]
fftBackward n lst = map (`div'` n) $ fft n ((C.conjugate (nthPrimitiveSquareRoot n))) lst

main = do
    (filePath:_) <- getArgs
    analyzeFile filePath

--loads the WaveFile from a file given a path
parseFile filePath = do
    fileContents <- B.readFile filePath
    let result = runGet getWaveFile fileContents
    putStrLn (filePath ++ " loaded!")
    return (result)

data WaveFile = WaveFile {
    chunkID ::      Int,
    fileSize ::     Int,
    riffType ::     Int,
    fmtId    ::     Int,
    fmtSize ::      Int,
    fmtCode ::      Int,
    channels ::     Int,
    sampleRate ::   Int,
    fmtAvgBPS ::    Int,
    fmtBlockAlign ::Int,
    bitDepth ::     Int,
    dataID ::       Int,
    dataSize ::     Int,
    fileData ::     [Float]
} -- deriving (Show)

instance Show WaveFile where
    show wav = "FileSize: " ++ (show $ fileSize wav) ++ "\nDataSize: " ++ (show $ dataSize wav) ++ "\nSampleRate: " ++ (show $ sampleRate wav) ++ "\nBitDepth: " ++ (show $ bitDepth wav) ++ "\nChannels: " ++ (show $ channels wav) ++ "\n"

-- parses the header of wav file and reads its contents
getWaveFile :: Get WaveFile
getWaveFile = do
    w_chunkID <- getWord32le
    let chunkID = (fromIntegral w_chunkID)
    w_fileSize <- getWord32le
    let fileSize = fromIntegral w_fileSize
    w_riffType <- getWord32le
    let riffType = fromIntegral w_riffType
    w_fmtID <- getWord32le
    let fmtID = fromIntegral w_fmtID
    w_fmtSize <- getWord32le
    let fmtSize = fromIntegral w_fmtSize
    w_fmtCode <- getWord16le
    let fmtCode = fromIntegral w_fmtCode
    w_channels <- getWord16le
    let channels = fromIntegral w_channels
    w_sampleRate <- getWord32le
    let sampleRate = fromIntegral w_sampleRate
    w_fmtAvgBPS <- getWord32le
    let fmtAvgBPS = fromIntegral w_fmtAvgBPS
    w_fmtBlockAlign <- getWord16le
    let fmtBlockAlign = fromIntegral w_fmtBlockAlign
    w_bitDepth <- getWord16le
    let bitDepth = fromIntegral w_bitDepth
    w_dataID <- getWord32le
    let dataID = fromIntegral w_dataID
    w_dataSize <- getWord32le
    let dataSize = fromIntegral w_dataSize
    fileData <- getWaveContents 0 dataSize
    return $! WaveFile chunkID fileSize riffType fmtID fmtSize fmtCode channels sampleRate fmtAvgBPS fmtBlockAlign bitDepth dataID dataSize fileData

--getWaveContents :: (Integral a) => a -> a-> Get [Float]
getWaveContents i size = do
    if i == size
       then return []
       else do
           fl <- getFloat
           rest <- getWaveContents (i+2) size
           return (fl:rest)

-- turns UInt16 to Int16
wordToInt :: Word16 -> Int
wordToInt word | word < notWord = (fromIntegral word :: Int)
               | otherwise      = (-1) * (fromIntegral notWord :: Int)
               where notWord = (-word)

-- read 16bit Float from the stream
getFloat :: Get Float
getFloat = do
    word <- getWord16le
    let signedWord = wordToInt word
    return ((fromIntegral signedWord :: Float)/(2^15))

-- Duration of the file
wavDuration :: WaveFile -> Float
wavDuration wav = ((fromIntegral $ dataSize wav) * (fromIntegral $ channels wav) * 8) / ((fromIntegral $ sampleRate wav) * (fromIntegral $ bitDepth wav) )

-- turns the list into a list of small intervals
intervalLength = 100

intervals_internal [] = []
intervals_internal lst@(head:tail) = (take (intervalLength) lst) : (intervals_internal tail)

intervals lst = filter (\l -> (length l) == (intervalLength)) (intervals_internal lst)

-- calculates the average of the absolute values in the list
my_average :: [Float] -> Float
my_average lst = (sum $ map (abs) lst) / (fromIntegral $ length lst)

-- get the interval of values from to in msec
getInterval :: Float -> Float -> WaveFile -> [Float]
getInterval from to wav = let rate = ((fromIntegral $ sampleRate wav) :: Float)
                              fromIndex = floor $ (rate * from)
                              toIndex = floor $ (rate * to)
                          in take (toIndex - fromIndex) (drop fromIndex $ fileData wav)

-- helper function for getBeepIntervals
getBeeps _ [] _ = []
getBeeps i (interval:tailOfIntervals) globalAverage = (i, if globalAverage < intervalAverage then True else False, head interval):(getBeeps (i+1) tailOfIntervals globalAverage)
                                                   where intervalAverage = my_average interval
-- handy function for debugging
indexToTime index wav = (fromIntegral index) / (fromIntegral $ sampleRate wav)

filterBeepValues [] = []
filterBeepValues ((index, isQuiet, value):tail)     |   (not isQuiet)     = []
                                                    |   otherwise         = (index, value):(filterBeepValues tail)

-- splits the wav file into into intervals and returns only non-quiet ones
-- it takes little windows and if the abs average value of the window > abs average of the file, the window is regarded as non-quite
getBeepIntervals wav = let wavData = fileData wav
                           globalAverage = my_average wavData
                           listOfIntervals = intervals wavData
                           beeps = getBeeps 0 listOfIntervals globalAverage
                           listOfBeeps = List.groupBy (\(_, isQuiet1, _) (_, isQuiet2, _) -> isQuiet1 == isQuiet2) beeps
                       in filter (\lst -> length lst > 9) $ filter (not . null) $ map (filterBeepValues) listOfBeeps
                           -- firstLast = map (\lst -> (head lst, last lst)) listOfBeeps
                       -- in  filter (\t -> (snd $ fst t)) firstLast

{-
DTMF keypad frequencies (with sound clips)
            5/1209 Hz 	6/1336 Hz 	7/1477 Hz 	8/1633 Hz
1/697 Hz  	 1 / 5   	 2 / 6    	 3 / 7   	 A / 8
2/770 Hz 	 4 / 10  	 5 / 12   	 6 / 14  	 B / 16
3/852 Hz 	 7 / 15 	 8 / 18   	 9 / 21 	 C / 24
4/941 Hz 	 * / 20 	 0 / 24 	 # / 28 	 D / 32

Each frequency button and frequency has its index. If you multiply the indices of two DTMF frequencies, you get index of the button
-}

-- constant - range around the precise DTMF freq
freqLimit = 15

dtmfFrequencies = [697,770,852,941,1209,1336,1477,1633]

dtmfFrequencyRanges_internal :: [Int] -> [Set.Set Int]
dtmfFrequencyRanges_internal [] = []
dtmfFrequencyRanges_internal (freq:rest) = (Set.fromList [(freq - freqLimit)..(freq + freqLimit)]):(dtmfFrequencyRanges_internal rest)

-- set for looking up if frequency is regarded as DTMF frequency or not
dtmfFrequencyRanges :: [Set.Set Int]
dtmfFrequencyRanges = dtmfFrequencyRanges_internal dtmfFrequencies

-- dictionary to look up chars representing buttons, we look up by the index of the button
dtmfTonesDict = Map.fromList $ zip ([x*y | x<- [1..4], y<-[5..8]]) ['1','2','3','A','4','5','6','B','7','8','9','C','*','0','#','D']

-- returns 0 if no otherwise index of frequency from table above
isDTMFFreq :: Int -> Int
isDTMFFreq freq = snd $ foldl f (1,0) dtmfFrequencyRanges
                where f acc@(index, result) set = if result /= 0 then acc
                                                  else
                                                    if freq `Set.member` set then (index+1, index)
                                                    else (index+1, 0)

-- given FFT Output, it returns recognized DTMF frequencies
analyzeFFTOutput :: Int -> Int -> [C.Complex Float] -> Int -> Float -> [(Int, Float)]
analyzeFFTOutput i n (x:xs) smplRate avrgValue   | i == n          = []
                                                 | i > (n `div` 2) = []
                                                 | otherwise       = let freq = (fromIntegral i) / ((fromIntegral n) / (fromIntegral smplRate))  -- calculate the frequency
                                                                         dfmtFreq = isDTMFFreq $ floor freq
                                                                         currentVal = (C.realPart $ abs x)
                                                                     in if dfmtFreq == 0 then restOfRecurrsion
                                                                        else if currentVal > avrgValue then (dfmtFreq, currentVal):(restOfRecurrsion)
                                                                             else restOfRecurrsion
                                                                     where restOfRecurrsion = analyzeFFTOutput (i+1) n xs smplRate avrgValue

maybeCharToChar :: Maybe Char -> Char
maybeCharToChar Nothing = '?'
maybeCharToChar (Just c) = c

-- analyzes interval of function and returns which button it is dialing
analyzeInterval :: [(Int, Float)] -> Int -> Char
analyzeInterval interval smplRate = let len = length interval
                                        floorPow2 = if len <= 1 then 0 else 2 ^ (floor $ logBase 2 (fromIntegral len)) -- calculate the n for FFT
                                        --startIndex = fst $ head interval
                                        --endIndex = fst $ last interval
                                        cmplxValues = map ((C.:+0).snd) interval -- make the values complex and trim the index
                                        fftImage = fftForward floorPow2 cmplxValues
                                        avrgValue = (sum $ map (C.realPart . abs) fftImage) / (fromIntegral floorPow2)
                                        frequencies = reverse $ List.sortBy (compare `DataFunction.on` snd) $ analyzeFFTOutput 0 floorPow2 fftImage smplRate avrgValue
                                        tone = if (length frequencies) < 2 then 0 else (fst $ head frequencies) * (fst $ head (tail frequencies))
                                    in  maybeCharToChar $ Map.lookup tone dtmfTonesDict

analyzeIntervals :: [[(Int, Float)]] -> Int -> String
analyzeIntervals [] _ = []
analyzeIntervals (interval:rest) smplRate = (analyzeInterval interval smplRate):(analyzeIntervals rest smplRate)

analyzeFile filePath = do
    wav <- parseFile filePath
    let beeps = getBeepIntervals wav
    putStrLn (analyzeIntervals beeps (sampleRate wav))

max_stars_repo_name: jendelel/DTMFAnalyzerHaskell

max_issues_repo_name: jendelel/DTMFAnalyzerHaskell

max_forks_repo_name: jendelel/DTMFAnalyzerHaskell


text:
module HopfieldMat ( initialWeights
                   , energy
                   , feed
                   , train
                   ) where

import Numeric.LinearAlgebra

initialWeights :: Int -> Matrix R
initialWeights l = (l><l) $ repeat 0

train :: Vector R -> Matrix R -> Matrix R
train p w = (outer p p) * zeroDiagonal + w
  where
    zeroDiagonal = (sp><sp)(repeat 1) - ident sp
    sp = size p

feed :: Matrix R -> Vector R -> Vector R
feed w = cmap signum . (w #>)

energy :: Matrix R -> Vector R -> Double
energy w p = -0.5 * (sumElements $ p * (w #> p))

max_stars_repo_name: mfandl/hml

max_issues_repo_name: mfandl/hml

max_forks_repo_name: mfandl/hml