text:
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE KindSignatures        #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PolyKinds             #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE StandaloneDeriving    #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}
{-|
Module      : Grenade.Layers.Concat
Description : Concatenation layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental

This module provides the concatenation layer, which runs two chilld layers in parallel and combines their outputs.
-}
module Grenade.Recurrent.Layers.ConcatRecurrent (
    ConcatRecurrent (..)
  ) where

import           Data.Serialize

import           Data.Singletons
import           GHC.TypeLits

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

import           Grenade.Core
import           Grenade.Recurrent.Core

import           Numeric.LinearAlgebra.Static (R, split, ( # ))

-- | A Concatentating Layer.
--
-- This layer shares it's input state between two sublayers, and concatenates their output.
--
-- With Networks able to be Layers, this allows for very expressive composition of complex Networks.
--
-- The Concat layer has a few instances, which allow one to flexibly "bash" together the outputs.
--
-- Two 1D vectors, can go to a 2D shape with 2 rows if their lengths are identical.
-- Any 2 1D vectors can also become a longer 1D Vector.
--
-- 3D images become 3D images with more channels. The sizes must be the same, one can use Pad
-- and Crop layers to ensure this is the case.
data ConcatRecurrent :: Shape -> Type -> Shape -> Type -> Type where
  ConcatRecLeft  :: !x -> !y -> ConcatRecurrent m (Recurrent x) n (FeedForward y)
  ConcatRecRight :: !x -> !y -> ConcatRecurrent m (FeedForward x) n (Recurrent y)
  ConcatRecBoth  :: !x -> !y -> ConcatRecurrent m (Recurrent x) n   (Recurrent y)

instance (Show x, Show y) => Show (ConcatRecurrent m (p x) n (q y)) where
  show (ConcatRecLeft x y)  = "ConcatRecLeft\n" ++ show x ++ "\n" ++ show y
  show (ConcatRecRight x y) = "ConcatRecRight\n" ++ show x ++ "\n" ++ show y
  show (ConcatRecBoth x y)  = "ConcatRecBoth\n" ++ show x ++ "\n" ++ show y

instance (RecurrentUpdateLayer x, UpdateLayer y) => UpdateLayer (ConcatRecurrent m (Recurrent x) n (FeedForward y)) where
  type Gradient (ConcatRecurrent m (Recurrent x) n (FeedForward y)) = (Gradient x, Gradient y)
  runUpdate lr (ConcatRecLeft x y) (x', y')
    = ConcatRecLeft (runUpdate lr x x') (runUpdate lr y y')

instance (RandomLayer x, RandomLayer y) => RandomLayer (ConcatRecurrent m (Recurrent x) n (FeedForward y)) where
  createRandomWith m gen
    = ConcatRecLeft <$> createRandomWith m gen <*> createRandomWith m gen

instance (UpdateLayer x, RecurrentUpdateLayer y) => UpdateLayer (ConcatRecurrent m (FeedForward x) n (Recurrent y)) where
  type Gradient (ConcatRecurrent m (FeedForward x) n (Recurrent y)) = (Gradient x, Gradient y)
  runUpdate lr (ConcatRecRight x y) (x', y')
    = ConcatRecRight (runUpdate lr x x') (runUpdate lr y y')

instance (RandomLayer x, RandomLayer y) => RandomLayer (ConcatRecurrent m (FeedForward x) n (Recurrent y)) where
  createRandomWith m gen
    = ConcatRecRight <$> createRandomWith m gen <*> createRandomWith m gen

instance (RecurrentUpdateLayer x, RecurrentUpdateLayer y) => UpdateLayer (ConcatRecurrent m (Recurrent x) n (Recurrent y)) where
  type Gradient (ConcatRecurrent m (Recurrent x) n (Recurrent y)) = (Gradient x, Gradient y)
  runUpdate lr (ConcatRecBoth x y) (x', y')
    = ConcatRecBoth (runUpdate lr x x') (runUpdate lr y y')

instance (RandomLayer x, RandomLayer y) => RandomLayer (ConcatRecurrent m (Recurrent x) n (Recurrent y)) where
  createRandomWith m gen
    = ConcatRecBoth <$> createRandomWith m gen <*> createRandomWith m gen

instance (RecurrentUpdateLayer x, UpdateLayer y) => RecurrentUpdateLayer (ConcatRecurrent m (Recurrent x) n (FeedForward y)) where
  type RecurrentShape (ConcatRecurrent m (Recurrent x) n (FeedForward y)) = RecurrentShape x

instance (UpdateLayer x, RecurrentUpdateLayer y) => RecurrentUpdateLayer (ConcatRecurrent m (FeedForward x) n (Recurrent y)) where
  type RecurrentShape (ConcatRecurrent m (FeedForward x) n (Recurrent y)) = RecurrentShape y

instance (RecurrentUpdateLayer x, RecurrentUpdateLayer y) => RecurrentUpdateLayer (ConcatRecurrent m (Recurrent x) n (Recurrent y)) where
  type RecurrentShape (ConcatRecurrent m (Recurrent x) n (Recurrent y)) =  RecurrentInputs '[ Recurrent x, Recurrent y ]

instance ( SingI i
         , Layer x i ('D1 m)
         , RecurrentLayer y i ('D1 n)
         , KnownNat o
         , KnownNat m
         , KnownNat n
         , o ~ (m + n)
         , n ~ (o - m)
         , (m <=? o) ~ 'True
         ) => RecurrentLayer (ConcatRecurrent ('D1 m) (FeedForward x) ('D1 n) (Recurrent y)) i ('D1 o) where
  type RecTape (ConcatRecurrent ('D1 m) (FeedForward x) ('D1 n) (Recurrent y)) i ('D1 o) = (Tape x i ('D1 m), RecTape y i ('D1 n))

  runRecurrentForwards (ConcatRecRight x y) s input =
    let (xT, xOut :: S ('D1 m))       = runForwards x input
        (yT, side, yOut :: S ('D1 n)) = runRecurrentForwards y s input
    in case (xOut, yOut) of
        (S1D xOut', S1D yOut') ->
            ((xT, yT), side, S1D (xOut' # yOut'))

  runRecurrentBackwards (ConcatRecRight x y) (xTape, yTape) s (S1D o) =
    let (ox :: R m , oy :: R n) = split o
        (x', xB :: S i)         = runBackwards x xTape (S1D ox)
        (y', side, yB :: S i)   = runRecurrentBackwards y yTape s (S1D oy)
    in  ((x', y'), side, xB + yB)

instance ( SingI i
         , RecurrentLayer x i ('D1 m)
         , Layer y i ('D1 n)
         , KnownNat o
         , KnownNat m
         , KnownNat n
         , o ~ (m + n)
         , n ~ (o - m)
         , (m <=? o) ~ 'True
         ) => RecurrentLayer (ConcatRecurrent ('D1 m) (Recurrent x) ('D1 n) (FeedForward y)) i ('D1 o) where
  type RecTape (ConcatRecurrent ('D1 m) (Recurrent x) ('D1 n) (FeedForward y)) i ('D1 o) = (RecTape x i ('D1 m), Tape y i ('D1 n))

  runRecurrentForwards (ConcatRecLeft x y) s input =
    let (xT, side, xOut :: S ('D1 m)) = runRecurrentForwards x s input
        (yT, yOut :: S ('D1 n))       = runForwards y input
    in case (xOut, yOut) of
        (S1D xOut', S1D yOut') ->
            ((xT, yT), side, S1D (xOut' # yOut'))

  runRecurrentBackwards (ConcatRecLeft x y) (xTape, yTape) s (S1D o) =
    let (ox :: R m , oy :: R n) = split o
        (x', side, xB :: S i)   = runRecurrentBackwards x xTape s (S1D ox)
        (y', yB :: S i)         = runBackwards y yTape (S1D oy)
    in  ((x', y'), side, xB + yB)

instance ( SingI i
         , RecurrentLayer x i ('D1 m)
         , RecurrentLayer y i ('D1 n)
         , Fractional (RecurrentShape x)
         , Fractional (RecurrentShape y)
         , KnownNat o
         , KnownNat m
         , KnownNat n
         , o ~ (m + n)
         , n ~ (o - m)
         , (m <=? o) ~ 'True
         ) => RecurrentLayer (ConcatRecurrent ('D1 m) (Recurrent x) ('D1 n) (Recurrent y)) i ('D1 o) where
  type RecTape (ConcatRecurrent ('D1 m) (Recurrent x) ('D1 n) (Recurrent y)) i ('D1 o) = (RecTape x i ('D1 m), RecTape y i ('D1 n))

  runRecurrentForwards (ConcatRecBoth x y) (sx :~@+> (sy :~@+> RINil)) input =
    let (xT, s'x, xOut :: S ('D1 m)) = runRecurrentForwards x sx input
        (yT, s'y, yOut :: S ('D1 n)) = runRecurrentForwards y sy input
    in case (xOut, yOut) of
        (S1D xOut', S1D yOut') ->
            ((xT, yT), (s'x :~@+> (s'y :~@+> RINil)), S1D (xOut' # yOut'))

  runRecurrentBackwards (ConcatRecBoth x y) (xTape, yTape) (sx :~@+> (sy :~@+> RINil)) (S1D o) =
    let (ox :: R m , oy :: R n) = split o
        (x', s'x, xB :: S i)    = runRecurrentBackwards x xTape sx (S1D ox)
        (y', s'y, yB :: S i)    = runRecurrentBackwards y yTape sy (S1D oy)
    in  ((x', y'), (s'x :~@+> (s'y :~@+> RINil)), xB + yB)

instance (Serialize a, Serialize b) => Serialize (ConcatRecurrent sa (FeedForward a) sb (Recurrent b)) where
  put (ConcatRecRight a b) = put a *> put b
  get = ConcatRecRight <$> get <*> get

instance (Serialize a, Serialize b) => Serialize (ConcatRecurrent sa (Recurrent a) sb (Recurrent b)) where
  put (ConcatRecBoth a b) = put a *> put b
  get = ConcatRecBoth <$> get <*> get

instance (Serialize a, Serialize b) => Serialize (ConcatRecurrent sa (Recurrent a) sb (FeedForward b)) where
  put (ConcatRecLeft a b) = put a *> put b
  get = ConcatRecLeft <$> get <*> get

-------------------- GNum instances --------------------

instance (GNum x, GNum y) => GNum (ConcatRecurrent m (p x) n (q y)) where
  n |* (ConcatRecLeft x y) = ConcatRecLeft (n |* x) (n |* y)
  n |* (ConcatRecRight x y) = ConcatRecRight (n |* x) (n |* y)
  n |* (ConcatRecBoth x y) = ConcatRecBoth (n |* x) (n |* y)
  (ConcatRecLeft x1 y1) |+ (ConcatRecLeft x2 y2) = ConcatRecLeft (x1 |+ x2) (y1 |+ y2)
  (ConcatRecRight x1 y1) |+ (ConcatRecRight x2 y2) = ConcatRecRight (x1 |+ x2) (y1 |+ y2)
  (ConcatRecBoth x1 y1) |+ (ConcatRecBoth x2 y2) = ConcatRecBoth (x1 |+ x2) (y1 |+ y2)

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
{-# LANGUAGE GADTs, TypeOperators, ScopedTypeVariables, Rank2Types #-}
module Math.Probably.GibbsUntyped where

import Math.Probably.Sampler
import Math.Probably.StochFun
import Numeric.LinearAlgebra hiding (flatten)
import Data.Array.ST
import Control.Monad.ST
import Control.Monad
import Data.STRef
import Data.Array
import GHC.Arr

type V = Vector Double

type VArr s = STArray s Int V

type Gibbs s = [VArr s -> Sampler V]

myG :: Gibbs s
myG =  [\_-> return 0.5,
        \_->return 0.2]

{-gibbsSampler :: GTerm a b -> StochFun a b
gibbsSampler gibbs = SF stochfun
    where stochfun (x0, dbls) = undefined

-}

--there is a serious flaw in this approach. we need to replace items
--in x0 when running subsequent samplers. this is considerably more
--difficult.
runGibbs :: Gibbs s -> VArr s -> [Double] -> ST s [VArr s]
runGibbs gbs x0 dbls = do

--what a nightmare. What about the zipper??
  return []

{-nearlyGibbs ::  Gibbs -> StochFun VArr VArr
nearlyGibbs gbs = SF $ \(arr, dbls') -> runST $ do
                    dbls <- newSTRef dbls'
                    melted <- thawSTArray arr
                    forM_ (zip gbs [0..]) $ \(f,idx) -> do
                      rnds <- readSTRef dbls
                      let (x, rnds) = unSF $ f melted
                    ice <- freezeSTArray melted
                    newDbls <- readSTRef dbls
                    return (ice, newDbls) -}

{-
                                                                    (y, dbls') = sf dbls
                                                next = unSF $ nearlyGibbs g
                                                (rest, dbls'') = next (x0, dbls')
                                             in (y :+: rest, dbls'') -}
--nearlyGibbs GNil = sampler $ return TNil

{-

myChain = Mrkv (nearlyGibbs myG2) (2:+:2:+:TNil) id

tst = take 20 `fmap` runMarkovIO myChain

{-iterateG :: GTerm a b -> a -> b
iterateG (f :%: g) x0 = f x0 :+: iterateG g x0
iterateG GNil _ = TNil

myXs :: [Double :+: (Double :+: TNil)]
myXs = iterate (iterateG myG2) (2 :+: (3 :+: TNil))
-}

gibbs2 :: (a-> b) -> (b-> a) -> (a,b) -> (a,b)
gibbs2 f g (x,y) = let y1 = f x
                       x1 = g y1 in (x1,y1)

gibbs3 :: ((a, b)->c) -> ((b,c)-> a) -> ((c,a)-> b) -> (a,b,c) -> (a,b,c)
gibbs3 f g h (x0,y0,z0) = let x1 = g (y0,z0)
                              y1 = h (z0, x1)
                              z1 = f (x1, y1)
                          in (x1,y1,z1)

asD :: Double -> Double
asD = id
-}

ext: hs

max_stars_repo_name: glutamate/probably-base

max_issues_repo_name: glutamate/probably-base

max_forks_repo_name: glutamate/probably-base


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE DeriveFunctor #-}

import Control.Monad ((>=>),forM)
import Control.Applicative (liftA2)
import Data.Complex
import Data.Either
import Data.Fixed
import Data.Functor.Const
import Data.Functor.Identity
import Data.HashMap.Strict
import Data.HashSet
import Data.Hashable
import Data.Int
import Data.IntMap.Strict
import Data.IntSet
import Data.Map
import Data.Monoid
import Data.Ord (Down(..))
import Data.Orphans ()
import Data.Proxy (Proxy(..))
import Data.Ratio
import Data.Semigroup
import Data.Semiring
import Data.Sequence
import Data.Set
import Data.Word
import GHC.Natural
import qualified Data.Foldable as F
import Prelude hiding (Num(..),(^))
import System.IO.Unsafe
import Test.QuickCheck hiding (Fixed)
import Test.QuickCheck.Gen (suchThat)
import Test.QuickCheck.Instances ()
import Test.Tasty (defaultMain, testGroup, TestTree)
import Text.Show.Functions ()
import qualified Data.Map as Map
import qualified Data.Sequence as Seq
import qualified GHC.Num as Num
import qualified Test.QuickCheck.Classes as QCC
import qualified Control.Exception as E

main :: IO ()
main = do
  QCC.lawsCheckMany namedTests
  F.fold pow_prop

forceError :: a -> IO (Either E.ErrorCall a)
forceError = E.try . E.evaluate

pow_prop :: [IO ()]
pow_prop = [quickCheck pow_prop1, quickCheck pow_prop2, quickCheck pow_prop3]

pow_prop1 :: Int -> Negative Int -> Property
pow_prop1 x y = ioProperty $ do
  p <- forceError (x ^ getNegative y)
  pure $ isLeft p

pow_prop2 :: Int -> Positive Int -> Property
pow_prop2 x y = ioProperty $ do
  p <- forceError (x ^ getPositive y)
  pure $ isRight p

pow_prop3 :: Int -> Property
pow_prop3 x = ioProperty $ do
  p <- forceError (x ^ 0)
  case p of
    Left e -> pure False
    Right x -> pure $ x == 1

type Laws = QCC.Laws

semiringLaws :: (Arbitrary a, Show a, Eq a, Semiring a) => Proxy a -> [Laws]
semiringLaws p = [QCC.semiringLaws p]

ringLaws :: (Arbitrary a, Show a, Eq a, Ring a) => Proxy a -> [Laws]
ringLaws p = [QCC.semiringLaws p, QCC.ringLaws p]

namedTests :: [(String, [Laws])]
namedTests =
  [ ("Bool", semiringLaws pBool)
--  , ("Double", ringLaws pDouble) -- needs to be within some epsilon
--  , ("Float", ringLaws pFloat)   -- needs to be within some epsilon
--  , ("Complex", ringLaws pComplex) -- needs to be within some epsilon
  , ("Int", ringLaws pInt)
  , ("Int8", ringLaws pInt8)
  , ("Int16", ringLaws pInt16)
  , ("Int32", ringLaws pInt32)
  , ("Int64", ringLaws pInt64)
  , ("Word", ringLaws pWord)
  , ("Word8", ringLaws pWord8)
  , ("Word16", ringLaws pWord16)
  , ("Word32", ringLaws pWord32)
  , ("Word64", ringLaws pWord64)
  , ("()", ringLaws pUnit)
  , ("Maybe", semiringLaws pMaybe)
  , ("PosRatio", semiringLaws pPosRatio)
  , ("IO", ringLaws pIO)
  , ("Fixed", ringLaws pFixed)
  , ("Identity", ringLaws pIdentity)
  , ("Dual", ringLaws pDual)
  , ("(->)", ringLaws pFunction)
  , ("Down", ringLaws pDown)
  , ("Const", ringLaws pConst)
--  , ("IntMap", semiringLaws pIntMap) -- needs newtypes
  , ("Set", semiringLaws pSet)
--  , ("IntSet", semiringLaws pIntSet) -- needs newtypes
  , ("HashSet", semiringLaws pHashSet)
  , ("HashMap", semiringLaws pHashMap)
  , ("Map", semiringLaws pMap)
  , ("Predicate", semiringLaws pPredicate)
  , ("Equivalence", semiringLaws pEquivalence)
  , ("Op", ringLaws pOp)
  , ("Ap", ringLaws pAp)

  , ("IntSet Sum", semiringLaws pIntSetSum)
  , ("IntSet Product", semiringLaws pIntSetProduct)
  , ("IntSet Min", semiringLaws pIntSetMin)
  , ("IntSet Max", semiringLaws pIntSetMax)

  , ("IntMap Sum", semiringLaws pIntMapSum)
  , ("IntMap Product", semiringLaws pIntMapProduct)
  , ("IntMap Min", semiringLaws pIntMapMin)
  , ("IntMap Max", semiringLaws pIntMapMax)
  ]

#if !(MIN_VERSION_base(4,12,0))
newtype Ap f a = Ap { getAp :: f a }
  deriving (Eq, Functor, Applicative, Show)
instance (Semiring a, Applicative f) => Semiring (Ap f a) where
  zero = pure zero
  one = pure one
  plus = liftA2 plus
  times = liftA2 times
#endif
deriving instance (Arbitrary (f a)) => Arbitrary (Ap f a)

newtype Predicate a = Predicate (a -> Bool)
  deriving (Eq, Show)
deriving instance Semiring (Predicate a)
deriving instance (Arbitrary a, CoArbitrary a) => Arbitrary (Predicate a)

newtype Equivalence a = Equivalence (a -> a -> Bool)
  deriving (Eq, Show)
deriving instance Semiring a => Semiring (Equivalence a)
deriving instance (Arbitrary a, CoArbitrary a) => Arbitrary (Equivalence a)

newtype Op a b = Op (b -> a)
  deriving (Eq, Show)
deriving instance Semiring a => Semiring (Op a b)
deriving instance Ring a => Ring (Op a b)
deriving instance (Arbitrary a, CoArbitrary a, CoArbitrary b) => Arbitrary (Op a b)

pAp = p @(Ap Identity Int)
pPredicate = p @(Predicate Int)
pEquivalence = p @(Equivalence Int)
pOp = p @(Op Int Int)

deriving instance Arbitrary a => Arbitrary (Down a)

-- | Rationals constrained to being positive
newtype PosRatio a = PosRatio { getPosRatio :: Ratio a }
  deriving (Eq, Show, Semiring)

instance (Arbitrary a, Integral a, Ord a) => Arbitrary (PosRatio a) where
  arbitrary =
    let nm = suchThat (arbitrary :: Gen a) (> 0)
        dm = suchThat (arbitrary :: Gen a) (> 0)
    in PosRatio <$> liftA2 (%) nm dm

p :: forall a. Proxy a
p = Proxy

instance (Arbitrary a, Eq b) => Eq (a -> b) where
  (==) f g =
    let x = generate arbitrary :: IO a
        x' = unsafePerformIO x
    in f x' == g x'

instance Show (Endo a) where
  show _ = "<<<Endo>>>"

deriving instance (Arbitrary a, Eq a) => Eq (Endo a)

-- no
instance Show (IO a) where
  show _ = "<<<IO action>>>"

-- no no
instance Eq a => Eq (IO a) where
  x == y = unsafePerformIO $ liftA2 (==) x y

-- this is fine
instance Arbitrary a => Arbitrary (IO a) where
  arbitrary = fmap pure (arbitrary :: Gen a)

deriving instance Hashable a => Hashable (Sum a)

deriving instance Arbitrary (IntSetOf a)
deriving instance Arbitrary v => Arbitrary (IntMapOf k v)

pBool = p @Bool
pDouble = p @Double
pFloat = p @Float
pInt = p @Int
pInt8 = p @Int8
pInt16 = p @Int16
pInt32 = p @Int32
pInt64 = p @Int64
pInteger = p @Integer
pNatural = p @Natural
pWord = p @Word
pWord8 = p @Word8
pWord16 = p @Word16
pWord32 = p @Word32
pWord64 = p @Word64
pUnit = p @()
pMaybe = p @(Maybe Int)
pPosRatio = p @(PosRatio Int)
pIO = p @(IO Int)
pComplex = p @(Complex Double)
pFixed = p @(Fixed E0)
pMin = p @(Min Int)
pMax = p @(Max Int)
pIdentity = p @(Identity Int)
pDual = p @(Dual Int)
pEndo = p @(Endo (Sum Int))
pSum  = p @(Sum Int)
pProduct = p @(Product Int)
pDown = p @(Down Int)
pSeq = p @(Seq Int)
pSet = p @(Set (Sum Int))
pHashSet = p @(HashSet (Sum Int))
pFunction = p @(Int -> Int)
pMap = p @(Map (Sum Int) Int)
pHashMap = p @(HashMap (Sum Int) Int)
pConst = p @(Const Int Int)
pAlt = p @(Alt Maybe Int)

pIntSetSum = p @(IntSetOf (Sum Int))
pIntSetProduct = p @(IntSetOf (Product Int))
pIntSetMin = p @(IntSetOf (Min Int))
pIntSetMax = p @(IntSetOf (Max Int))

pIntMapSum = p @(IntMapOf (Sum Int) Int)
pIntMapProduct = p @(IntMapOf (Product Int) Int)
pIntMapMin = p @(IntMapOf (Min Int) Int)
pIntMapMax = p @(IntMapOf (Max Int) Int)

max_stars_repo_name: erikd/semirings

max_issues_repo_name: erikd/semirings

max_forks_repo_name: erikd/semirings


text:
{-# LANGUAGE ConstraintKinds, FlexibleContexts, FlexibleInstances, GADTs   #-}
{-# LANGUAGE KindSignatures, MultiParamTypeClasses                         #-}
{-# LANGUAGE NoMonomorphismRestriction, TypeFamilies, TypeSynonymInstances #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
module Algebra.Matrix (Matrix(..), delta, companion,
                       gaussReduction, maxNorm, rankWith, det,
                       inverse, inverseWith) where
import           Algebra.Internal
import qualified Algebra.LinkedMatrix as LM
import           Algebra.Prelude.Core hiding (maxNorm, zero)

import           Control.Lens                (both, view, (%~), (&), _3)
import           Control.Monad               (when)
import qualified Data.Matrix                 as DM
import qualified Data.Vector                 as V
import           GHC.Exts                    (Constraint)
import           Numeric.Algebra             (Additive, DecidableZero)
import           Numeric.Algebra             (Monoidal, Multiplicative)
import           Numeric.Algebra             (Unital)
import qualified Numeric.Algebra             as NA
import qualified Numeric.Decidable.Zero      as NA
import qualified Numeric.LinearAlgebra       as LA
import qualified Numeric.LinearAlgebra.Devel as LA
import qualified Prelude                     as P

{-# ANN module "HLint: ignore Redundant bang pattern" #-}

class Matrix mat where
  type Elem mat a :: Constraint
  cmap  :: (Elem mat a, Elem mat b) => (a -> b) -> mat a -> mat b
  empty :: Elem mat b => mat b
  fromLists :: Elem mat a => [[a]] -> mat a
  fromCols  :: Elem mat a => [V.Vector a] -> mat a
  fromCols [] = zero 0 0
  fromCols xs = foldr1 (<||>) $ map colVector xs
  fromRows :: Elem mat a => [V.Vector a] -> mat a
  fromRows [] = zero 0 0
  fromRows xs = foldr1 (<-->) $ map rowVector xs
  toCols :: Elem mat a => mat a -> [V.Vector a]
  toCols m = map (`getCol` m) [1..ncols m]
  toRows :: Elem mat a => mat a -> [V.Vector a]
  toRows m = map (`getRow` m) [1..nrows m]
  ncols :: mat a -> Int
  nrows :: mat a -> Int
  identity :: Elem mat a => Int -> mat a
  diag     :: Elem mat a => V.Vector a -> mat a
  getDiag :: Elem mat a => mat a -> V.Vector a
  trace :: Elem mat a => mat a -> a
  diagProd :: Elem mat a => mat a -> a
  zero :: Elem mat a => Int -> Int -> mat a
  colVector :: Elem mat a => V.Vector a -> mat a
  rowVector :: Elem mat a => V.Vector a -> mat a
  getCol :: Elem mat a => Int -> mat a -> V.Vector a
  getRow :: Elem mat a => Int -> mat a -> V.Vector a
  switchRows :: Elem mat a => Int -> Int -> mat a -> mat a
  scaleRow :: Elem mat a => a -> Int -> mat a -> mat a
  combineRows :: Elem mat a => Int -> a -> Int -> mat a -> mat a
  trans :: Elem mat a => mat a -> mat a
  buildMatrix :: Elem mat a => Int -> Int -> ((Int, Int) -> a) -> mat a
  index :: Elem mat a => Int -> Int -> mat a -> Maybe a
  index i j m = if 1 <= i && i <= nrows m && 1 <= j && j <= ncols m
                then Just $ m ! (i, j)
                else Nothing
  (!) :: Elem mat a => mat a -> (Int, Int) -> a
  (<||>) :: Elem mat a => mat a -> mat a -> mat a
  (<-->) :: Elem mat a => mat a -> mat a -> mat a
  nonZeroRows :: (DecidableZero a, Elem mat a) => mat a -> [Int]
  nonZeroRows = map fst . filter (V.any (not . NA.isZero) . snd) . zip [1..] . toRows
  nonZeroCols :: (DecidableZero a, Elem mat a) => mat a -> [Int]
  nonZeroCols = map fst . filter (V.any (not . NA.isZero) . snd) . zip [1..] . toCols

instance Matrix DM.Matrix where
  type Elem DM.Matrix a = P.Num a
  empty = DM.zero 0 0
  cmap  = fmap
  fromLists = DM.fromLists
  ncols = DM.ncols
  nrows = DM.nrows
  trans = DM.transpose
  identity = DM.identity
  diag v = DM.matrix (V.length v) (V.length v) $ \(i, j) ->
    if i == j then v V.! (i-1) else 0
  getDiag = DM.getDiag
  trace = DM.trace
  diagProd = DM.diagProd
  zero = DM.zero
  colVector = DM.colVector
  rowVector = DM.rowVector
  getCol = DM.getCol
  getRow = DM.getRow
  switchRows = DM.switchRows
  combineRows = DM.combineRows
  scaleRow = DM.scaleRow
  buildMatrix = DM.matrix
  (!) = (DM.!)
  (<||>) = (DM.<|>)
  (<-->) = (DM.<->)

swapIJ :: Eq a => a -> a -> a -> a
swapIJ i j k
  | k == i = j
  | k == j = i
  | otherwise = k

instance Matrix LA.Matrix where
  type Elem LA.Matrix a = (P.Num a, LA.Numeric  a, LA.Element a, LA.Container LA.Vector a)
  empty = LA.fromLists [[]]
  fromLists = LA.fromLists
  cmap = LA.cmap
  ncols = LA.cols
  nrows = LA.rows
  trans = LA.tr
  identity = LA.ident
  fromCols = LA.fromColumns . map (LA.fromList . V.toList)
  diag = LA.diag . LA.fromList . V.toList
  getDiag = V.fromList . LA.toList . LA.takeDiag
  trace = LA.sumElements . LA.takeDiag
  diagProd = LA.prodElements . LA.takeDiag
  zero i j = LA.konst 0 (i, j)
  colVector = LA.asColumn . LA.fromList . V.toList
  rowVector = LA.asRow . LA.fromList . V.toList
  toCols = map (V.fromList . LA.toList) . LA.toColumns
  toRows = map (V.fromList . LA.toList) . LA.toRows
  getCol i = V.fromList . LA.toList . (!! (i - 1)) . LA.toColumns
  getRow i = V.fromList . LA.toList . (!! (i - 1)) . LA.toRows
  switchRows i j m = m LA.? map (swapIJ (i-1) (j-1)) [0.. nrows m - 1]
  combineRows j s i m = LA.mapMatrixWithIndex (\(k,l) v -> if k == j - 1 then s P.* (m ! (i,l+1)) P.+ v else v) m
  buildMatrix w h f = LA.build (w, h) (\i j -> f (toIntLA i+1, toIntLA j+1))
  scaleRow a i = (fst .) $ LA.mutable $ \(k, l) m -> do
    v <- LA.readMatrix m k l
    when (k == i - 1) $
      LA.writeMatrix m k l (a P.* v)
  m ! (i, j) = m `LA.atIndex` (i - 1, j - 1)
  m <||> n = LA.fromColumns $ LA.toColumns m ++ LA.toColumns n
  m <--> n = LA.fromRows $ LA.toRows m ++ LA.toRows n

toIntLA :: LA.Container LA.Matrix e => e -> Int
toIntLA e = fromIntegral $ LA.toZ ((1 LA.>< 1) [e]) `LA.atIndex` (0, 0)

instance Matrix LM.Matrix where
  type Elem LM.Matrix a = (CoeffRing a, Unital a, Monoidal a, Multiplicative a, Additive a)
  cmap = LM.cmap
  (!) m pos = m LM.! (pos & both %~ pred)
  index i j = LM.index (i-1) (j-1)
  empty = LM.empty
  buildMatrix h w f = LM.fromList [((i, j), f (i,j)) | j <- [1..w], i <- [1..h]]
  trans = LM.transpose
  combineRows j s i = LM.combineRows s (i-1) (j-1)
  switchRows i j = LM.switchRows (i-1) (j-1)
  scaleRow a i = LM.scaleRow a (pred i)
  zero = LM.zeroMat
  identity = LM.identity
  diag = LM.diag
  getDiag = LM.getDiag
  diagProd = LM.diagProd
  trace = LM.trace
  getRow = LM.getRow . pred
  getCol = LM.getCol . pred
  nrows = LM.nrows
  ncols = LM.ncols
  fromLists = LM.fromLists
  (<||>) = (LM.<||>)
  (<-->) = (LM.<-->)
  rowVector = LM.rowVector
  colVector = LM.colVector
  nonZeroRows = LM.nonZeroRows
  nonZeroCols = LM.nonZeroCols

delta :: (NA.Monoidal r, NA.Unital r) => Int -> Int -> r
delta i j | i == j = NA.one
          | otherwise = NA.zero

companion :: (KnownNat n, CoeffRing r, Matrix mat,
              Elem mat r, IsMonomialOrder n ord)
          => Ordinal n -> OrderedPolynomial r ord n -> mat r
companion odn poly =
  let deg = totalDegree' poly
      vx  = leadingMonomial (var odn `asTypeOf` poly)
  in buildMatrix deg deg $ \(j, k) ->
  if 1 <= k && k <= deg - 1
  then delta j (k+1)
  else NA.negate $ coeff (NA.pow vx (fromIntegral $ j-1 :: NA.Natural)) poly

-- | @gaussReduction a = (a', p)@ where @a'@ is row echelon form and @p@ is pivoting matrix.
gaussReduction :: (Matrix mat, Elem mat a, Normed a, Eq a, NA.Field a)
               => mat a -> (mat a, mat a)
gaussReduction mat =
  let (a, b, _) = gaussReduction' mat in (a, b)

-- | @gaussReduction a = (a', p)@ where @a'@ is row echelon form and @p@ is pivoting matrix.
gaussReduction' :: (Matrix mat, Elem mat a, Normed a, Eq a, NA.Field a)
               => mat a -> (mat a, mat a, a)
gaussReduction' mat = {-# SCC "gaussRed" #-} go 1 1 mat (identity $ nrows mat) NA.one
  where
    go i j a p acc
      | i > nrows mat || j > ncols mat = (a, p, acc)
      | otherwise =
        let (k, new) = maximumBy (comparing $ norm . snd) [(l, a ! (l, j)) | l <- [i..nrows mat]]
        in if new == NA.zero
           then go i (j + 1) a p NA.zero
           else let prc l a0 p0
                      | l == i = prc (l+1) a0 p0
                      | l > nrows mat = (a0, p0)
                      | otherwise     =
                        let coe = NA.negate (a0 ! (l, j))
                            a'' = combineRows l coe i a0
                            p'' = combineRows l coe i p0
                        in prc (l+1) a'' p''
                    (a', p') = prc 1 (scaleRow (NA.recip new) i $ switchRows i k a)
                                     (scaleRow (NA.recip new) i $ switchRows i k p)
                    offset = if i == k then id else NA.negate
                in go (i+1) (j+1) a' p' (offset $ acc NA.* new)

det :: (Elem mat a, Eq a, NA.Field a, Normed a, Matrix mat)
    => mat a -> a
det = view _3 . gaussReduction'

maxNorm :: (Elem mat a, Normed a, Matrix mat) => mat a -> Norm a
maxNorm = maximum . concat . map (map norm . V.toList) . toRows

rankWith :: (Elem mat r, CoeffRing r, Matrix mat)
         => (mat r -> mat r) -> mat r -> Int
rankWith gauss = length . nonZeroRows . gauss

inverse :: (Elem mat a, Eq a, NA.Field a, Normed a, Matrix mat)
        => mat a -> mat a
inverse = snd . gaussReduction

inverseWith :: (mat a -> (mat a, mat a)) -> mat a -> mat a
inverseWith = (snd .)

max_stars_repo_name: hangingman/computational-algebra

max_issues_repo_name: hangingman/computational-algebra

max_forks_repo_name: hangingman/computational-algebra


text:
module Types where

import System.Random
import Numeric.Vector

data Camera = Camera { c_getPos :: Vec3d
                     , c_getForward :: Vec3d
                     , c_getUp :: Vec3d
                     , c_getRight :: Vec3d
                     , c_getFOV :: Double
                     }
            | OrthoCamera { c_getPos :: Vec3d
                          , c_getForward :: Vec3d
                          , c_getUp :: Vec3d
                          , c_getRight :: Vec3d
                          , c_getSize :: Double
                          } deriving Show

data Scene = Scene { s_getCamera :: Camera
                   , s_getObjects :: [Object]
                   , s_getLights :: [Light]
                   , s_getSkyColor :: Vec3d
                   , s_getWidth :: Int
                   , s_getHeight :: Int
                   , s_getSeed :: Int
                   , s_getSamples :: Int
                   , s_getBounces :: Int
                   } deriving Show

data SceneContext = SceneContext { ss_getScene :: Scene
                                 , ss_getGen :: StdGen
                                 , ss_getHit :: RayHit
                                 } deriving Show

data Ray = Ray Vec3d Vec3d deriving Show
r_getOrigin :: Ray -> Vec3d
r_getOrigin (Ray o _) = o
r_getDir :: Ray -> Vec3d
r_getDir (Ray _ d) = d

data RayHit = RayHit { rh_getPos :: Vec3d
                     , rh_getRay :: Ray
                     , rh_getDistance :: Double
                     , rh_getTri :: Tri
                     , rh_getObj :: Object
                     , rh_getShader :: BRDF
                     , rh_getTexCoord :: Vec2d
                     , rh_getNormal :: Vec3d
                     , rh_getTangent :: Vec3d
                     , rh_getBitangent :: Vec3d
                     , rh_getColor :: Vec3d
                     } deriving Show

data SVOp = VAdd | VSub | VDot | VCross | VMul | VDiv | VAbs | VMod deriving Show
data SOp = Add | Sub | Mul | Div | Abs | Mod | LessThan | GreaterThan deriving Show
data SVector = UV
             | Position
             | Normal
             | VecConst Double Double Double
             | VecMath SVOp SVector SVector
             | CombineXYZ SValue SValue SValue
             deriving Show
data SValue = ValConst Double
            | RayLength
            | ValMath SOp SValue SValue
            | SeparateX SVector | SeparateY SVector | SeparateZ SVector
            | SeparateR SColor | SeparateG SColor | SeparateB SColor
            deriving Show
data SColor = ColVertex
            | ColVector SVector
            | ColConst Double Double Double
            | CombineRGB SValue SValue SValue
            deriving Show
data BRDF = Emission SColor SValue
          | Diffuse SColor
          | Glossy SColor SValue
          | BRDFEmpty
          deriving Show

data Light = Directional { l_getAngle :: Double
                         , l_getStrength :: Double
                         , l_getDirection :: Vec3d
                         , l_getColor :: Vec3d
                         }
           | Point { l_getRadius :: Double
                   , l_getStrength :: Double
                   , l_getColor :: Vec3d
                   } deriving Show

data Object = Object { o_getTris :: [Tri]
                     , o_getShader :: BRDF
                     } deriving Show

data Tri = Tri { t_getP0 :: Vertex
               , t_getP1 :: Vertex
               , t_getP2 :: Vertex
               } deriving Show

data Vertex = Vertex { v_getPos :: Vec3d
                     , v_getColor :: Vec3d
                     , v_getTexCoord :: Vec2d
                     , v_getNormal :: Vec3d
                     , v_getTangent :: Vec3d
                     , v_getBitangent :: Vec3d
                     } deriving Show


text:
-- Kolmogorov-Smirnov tests for distribution
--
-- Note that it's not most powerful test for normality.
module KS (
  tests
  ) where

import qualified Data.Vector.Unboxed as U

import Statistics.Test.KolmogorovSmirnov

import Statistics.Distribution
import Statistics.Distribution.Binomial
import Statistics.Distribution.Exponential
import Statistics.Distribution.Gamma
import Statistics.Distribution.Normal
import Statistics.Distribution.Uniform
import Statistics.Distribution.Beta

import qualified System.Random.MWC               as MWC
import qualified System.Random.MWC.Distributions as MWC

import Test.HUnit hiding (Test)
import Test.Framework
import Test.Framework.Providers.HUnit

tests :: MWC.GenIO -> Test
tests g = testGroup "Kolmogorov-Smirnov"
  [ testCase "standard"           $ testKS standard           MWC.standard    g
  , testCase "normal m=1 s=2"     $ testKS (normalDistr 1 2) (MWC.normal 1 2) g
    -- Gamma distribution
  , testCase "gamma  k=1   θ=1"   $ testKS (gammaDistr  1   1  ) (MWC.gamma  1   1  ) g
  , testCase "gamma  k=0.3 θ=0.4" $ testKS (gammaDistr  0.3 0.4) (MWC.gamma  0.3 0.4) g
  , testCase "gamma  k=0.3 θ=3"   $ testKS (gammaDistr  0.3 3  ) (MWC.gamma  0.3 3  ) g
  , testCase "gamma  k=3   θ=0.4" $ testKS (gammaDistr  3   0.4) (MWC.gamma  3   0.4) g
  , testCase "gamma  k=3   θ=3"   $ testKS (gammaDistr  3   3  ) (MWC.gamma  3   3  ) g
    -- Uniform
  , testCase "uniform -2 .. 3"    $ testKS (uniformDistr (-2) 3) (MWC.uniformR (-2,3)) g
    -- Exponential
  , testCase "exponential l=1"    $ testKS (exponential 1)       (MWC.exponential 1) g
  , testCase "exponential l=3"    $ testKS (exponential 3)       (MWC.exponential 3) g
    -- Beta
  , testCase "beta a=0.3,b=0.5"    $ testKS (betaDistr 0.3 0.5)       (MWC.beta 0.3 0.5) g
  , testCase "beta a=0.1,b=0.8"    $ testKS (betaDistr 0.3 0.5)       (MWC.beta 0.3 0.5) g
  , testCase "beta a=0.8,b=0.1"    $ testKS (betaDistr 0.3 0.5)       (MWC.beta 0.3 0.5) g
  ]

testKS :: (Distribution d) => d -> (MWC.GenIO -> IO Double) -> MWC.GenIO -> IO ()
testKS distr generator g = do
  sample <- U.replicateM 1000 (generator g)
  case kolmogorovSmirnovTest distr 0.01 sample of
    Significant    -> assertFailure "KS test failed"
    NotSignificant -> return ()

max_stars_repo_name: OlivierSohn/mwc-random

max_issues_repo_name: OlivierSohn/mwc-random

max_forks_repo_name: OlivierSohn/mwc-random


text:
{-# LANGUAGE RecordWildCards #-}
module Main where

import Control.Monad.Trans.Except
import System.IO
import Data.Word
import Data.Complex
import Foreign.C.Types
import Data.Monoid

import Control.Error.Util
import Foreign.Storable.Complex
import Options.Applicative
import Pipes as P
import Pipes.Prelude as P
import Data.Vector.Storable as VS hiding ((++))

import SDR.Util
import SDR.RTLSDRStream
import SDR.ArgUtils
import SDR.Serialize as S

data Options = Options {
    fileName   :: FilePath,
    frequency  :: Word32,
    sampleRate :: Word32,
    size       :: Maybe Int
}

optParser :: Parser Options
optParser = Options
          <$> strOption (
                 long "output  "
              <> short 'o'
              <> metavar "FILENAME"
              <> help "Output filename"
              )
          <*> option (fmap fromIntegral parseSize) (
                 long "frequency"
              <> short 'f'
              <> metavar "FREQUENCY"
              <> help "Frequency to tune to"
              )
          <*> option (fmap fromIntegral parseSize) (
                 long "samplerate"
              <> short 'r'
              <> metavar "RATE"
              <> help "Sample rate"
              )
          <*> optional (option (fmap fromIntegral parseSize) (
                 long "size"
              <> short 's'
              <> metavar "SIZE"
              <> help "Size of output file. If omitted, samples will be recorded until the program is killed"
              ))

opt :: ParserInfo Options
opt = info (helper <*> optParser) (fullDesc <> progDesc "Record IQ samples from an RTL2832 based device" <> header "RTLSDR Record")

doIt Options{..} = do
    str <- sdrStream (defaultRTLSDRParams frequency sampleRate) 1 16384
    lift $ withFile fileName WriteMode $ \handle ->
        runEffect $ str >-> maybe P.cat P.take size >-> P.map (interleavedIQUnsigned256ToFloat :: VS.Vector CUChar -> VS.Vector (Complex CFloat)) >-> S.toHandle handle

main = execParser opt >>= exceptT putStrLn return . doIt

max_stars_repo_name: adamwalker/sdr-apps

max_issues_repo_name: adamwalker/sdr-apps

max_forks_repo_name: adamwalker/sdr-apps


text:
module Main where

import           Control.Monad.Par
import           Data.Monoid                              ( (<>) )
import           Types
import qualified Text.Printf                   as Printf
import           Graphics.Gloss
import           Graphics.Gloss.Data.Color
import           Graphics.Gloss.Data.ViewPort
import qualified Numeric.LinearAlgebra         as LinAlg

windowSize :: (Int, Int)
windowSize = (1200, 800)

fps :: Int
fps = 10

data Model = Model
  { earth :: Point
  , earthSpeed :: LinAlg.Vector Float
  , jupiter :: Point
  , jupiterSpeed :: LinAlg.Vector Float
  } deriving (Show, Eq)

earthPosString :: Model -> String
earthPosString m =
  Printf.printf "earth (%.1f, %.1f)" (fst . earth $ m) (snd . earth $ m)

jupiterPosString :: Model -> String
jupiterPosString m =
  Printf.printf "jupiter (%.1f, %.1f)" (fst . jupiter $ m) (snd . jupiter $ m)

initialModel = Model
  { earth        = (-200, 0)
  , earthSpeed   = LinAlg.fromList [0, 20]
  , jupiter      = (0, 0)
  , jupiterSpeed = LinAlg.fromList [0, 0]
  }

main :: IO ()
main = simulate (InWindow "Newtons Planets" windowSize (200, 500))
                black
                fps
                initialModel
                drawModel
                stepFunction

radiusScaleFactor :: Integer -> Float
radiusScaleFactor diameter = fromIntegral diameter / distanceScale

earthRadius :: Float
earthRadius = radiusScaleFactor . diameter $ Earth

jupiterRadius :: Float
jupiterRadius = radiusScaleFactor . diameter $ Jupiter

drawModel :: Model -> Picture
drawModel model =
  let earth'               = color blue (circle earthRadius)
      jupiter'             = color cyan (circle jupiterRadius)
      (earthX  , earthY  ) = earth model
      (jupiterX, jupiterY) = jupiter model
      earthPos             = translate earthX earthY earth'
      jupiterPos           = translate jupiterX jupiterY jupiter'
  in  earthPos <> jupiterPos <> modelInfo (earthPosString model) (-580) (-350) <>
        modelInfo (jupiterPosString model) (-580) (-300)

modelInfo :: String -> Integer -> Integer -> Picture
modelInfo text a b =
  let text'  = Text text
      scaled = scale 0.1 0.1 text'
      moved  = translate (fromIntegral a) (fromIntegral b) scaled
  in  color red moved

stepFunction :: ViewPort -> Float -> Model -> Model
stepFunction vp secs model =
  let earthPoint    = earth model
      jupiterPoint  = jupiter model
      earthSpeed'   = earthSpeed model
      jupiterSpeed' = jupiterSpeed model
      (earthPulled, jupiterPulled)  = runPar $ do
        p1 <- spawnP $ pullOn Earth earthPoint Jupiter jupiterPoint
        p2 <- spawnP $ pullOn Jupiter jupiterPoint Earth earthPoint
        p1' <- get p1
        p2' <- get p2
        return (p1', p2')

  in  Model
        { earth        = newPoint earthPoint earthSpeed'
        , earthSpeed   = alterSpeed earthSpeed' earthPulled
        , jupiter      = newPoint jupiterPoint jupiterSpeed'
        , jupiterSpeed = alterSpeed jupiterSpeed' jupiterPulled
        }

max_stars_repo_name: mariatsji/planets

max_issues_repo_name: mariatsji/planets

max_forks_repo_name: mariatsji/planets


text:
module Shuffle where

import Numeric.LinearAlgebra ((<>))
import qualified System.Random as Random

import qualified Cube as Cube
import Rotation (MatrixRotation)
import qualified Rotation as Rotation

shuffleRotations :: [MatrixRotation]
shuffleRotations = [
    Rotation.rightL
  , Rotation.rightR
  , Rotation.leftL
  , Rotation.leftR
  , Rotation.topL
  , Rotation.topR
  , Rotation.bottomL
  , Rotation.bottomR
  , Rotation.frontL
  , Rotation.frontR
  , Rotation.backL
  , Rotation.backR
  ]

randomRotations :: Int -> [MatrixRotation]
randomRotations seed =
  let gen = Random.mkStdGen seed
      indices = map (\x -> (abs x) `mod` length shuffleRotations) $ Random.randoms gen
  in map (shuffleRotations !!) indices

exampleShuffle =
  let rotations = take 20 $ randomRotations 5
      names = map Rotation.nameMatrixRotation rotations
      rotation = foldl (<>) Rotation.identity rotations
      cube = Rotation.rotate rotation Cube.solvedCube
  in (cube, names)

max_stars_repo_name: runjak/hRubiks

max_issues_repo_name: runjak/hRubiks

max_forks_repo_name: runjak/hRubiks


text:
module AI.GeneticAlgorithm
    ( Genome
    , Individual
    , Population
    , individualGenome
    , individualFitness
    , mutate
    , proportionalSelection
    , getBasePopulation
    , evolve
    , evolves
    , bestAverageWorst ) where

import           AI.NeuralNetwork            as NN
import           Control.Lens
import           Control.Monad.State
import           Control.Parallel.Strategies
import           Data.Function
import           Data.Functor.Identity
import           Data.Maybe
import           Data.Random.Normal
import qualified Data.Vector                 as V
import qualified Data.Vector.Generic         as VG
import           Numeric.LinearAlgebra
import           System.Random

type Genome       = NN.Network
-- An individual with a genome and fitness
data Individual a = Individual Genome a
    deriving (Eq, Show)
type Population a = V.Vector (Individual a)

mut :: Double -> Double -> Double -> State StdGen Double
mut rate amount val = do
    rRate <- state random
    if rRate < rate
       then fmap (\x -> val + x * amount) (state normal)
       else return val

toIndividual :: (Genome -> a) -> Genome -> Individual a
toIndividual fitfunc genome = Individual genome (fitfunc genome)

toPopulation :: (Genome -> a) -> V.Vector Genome -> V.Vector (Individual a)
toPopulation fitfunc genomes = let pop = V.map (toIndividual fitfunc) genomes
                                in pop `using` parTraversable rseq

individualGenome :: Individual a -> Genome
individualGenome (Individual g _) = g

individualFitness :: Individual a -> a
individualFitness (Individual _ f) = f

mutate :: Double -> Double -> Genome -> State StdGen Genome
mutate rate amount genome = do
    let cmut = mut rate amount
    newBiases  <- mapM (VG.mapM    cmut) (genome^.biases)
    newWeights <- mapM (mapMatrixM cmut) (genome^.weights)
    let newNN = execState (biases .= newBiases >> weights .= newWeights) genome
    return newNN

proportionalSelection :: (Random a, Ord a, Num a, NFData a) =>
    Population a -> State StdGen (Population a)
proportionalSelection pop =
    let fitness  = individualFitness <$> pop
        accFit   = V.scanl1 (+) fitness
        choose p = snd $ fromJust $ V.find (\t -> p <= fst t) (V.zip accFit pop)
     in forM pop $ \_ -> do
         p <- state $ randomR (0, V.last accFit)
         return $ choose p

getBasePopulation :: (Genome -> a) -> State StdGen Genome ->
    Int -> State StdGen (Population a)
getBasePopulation fitfunc genomeGen popSize =
    toPopulation fitfunc <$> V.replicateM popSize genomeGen

evolve :: (Genome -> Double) -> (Genome -> State StdGen Genome) ->
    Population Double -> State StdGen (Population Double)
evolve fitfunc mutfunc pop = do
    selected <- proportionalSelection pop
    mutated  <- mapM (\(Individual g _) -> mutfunc g) selected
    return $ toPopulation fitfunc mutated

evolves :: Int -> Int -> State StdGen Genome ->
    (Genome -> Double) -> (Genome -> State StdGen Genome) ->
        State StdGen [Population Double]
evolves generations popSize genomeGen fitfunc mutfunc = do
    basePop <- getBasePopulation fitfunc genomeGen popSize
    let evofunc = evolve fitfunc mutfunc
    evolutions <- V.iterateNM generations evofunc basePop
    return $ V.toList evolutions

bestAverageWorst :: (Fractional a, Ord a) =>
    [Population a] -> [(Individual a, a, Individual a)]
bestAverageWorst = reverse . foldl (\acc pop ->
        let fitness = individualFitness <$> pop
            best = V.maximumBy (compare `on` individualFitness) pop
            avgFit = sum fitness / fromIntegral (V.length fitness)
            worst = V.minimumBy (compare `on` individualFitness) pop
         in (best, avgFit, worst) : acc) []

mapMatrixM :: (Element a, Element b, Monad m) =>
    (a -> m b) -> Matrix a -> m (Matrix b)
mapMatrixM f m = let ll  = toLists m
                     mmM = mapM . mapM
                  in fromLists <$> mmM f ll

max_stars_repo_name: cornelius-sevald/neurocar

max_issues_repo_name: cornelius-sevald/neurocar

max_forks_repo_name: cornelius-sevald/neurocar


text:
module Regression
 ( simpleRegression ) where

import Data.List
import Data.Maybe
import Control.Monad
import Statistics

simpleRegression :: (Floating a) => [a] -> [a] -> (Maybe a, Maybe a)
simpleRegression x y  = (slope, intercept)
                  where
                    slope      = (/) <$> covariance x y <*> variance x
                    intercept  = subtract <$> mean y <*> ((*) <$> slope <*> mean x)


text:
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE BangPatterns #-}

module Main where

import Options.Applicative
import Data.Semigroup ((<>))
import Chart.Heatmap
import Control.Lens
import Control.Monad (forM_)
import Data.List (maximumBy)
import Data.Ord (comparing)
import Data.Colour.SRGB (sRGB, RGB(..))
import Data.Colour.RGBSpace.HSL (hsl)
import Data.Complex (magnitude, realPart, imagPart, Complex(..))
import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Cairo (toFile, FileOptions(..), FileFormat(..))
import Data.Array.CArray.Base (unsafeForeignPtrToCArray, toForeignPtr, CArray(..), ixmapWithIndP)
import Data.Fixed (mod')
import Math.FFT (dftRC)
import Numeric.Signal (analytic_signal, auto_correlation, deriv, Filterable)
import Data.Maybe (catMaybes)
import Sound.MIDI.File
import qualified Data.Vector.Storable as VS
import qualified Data.Vector as V
import qualified Sound.File.Sndfile as SF
import qualified Sound.File.Sndfile.Buffer.Vector as SFV
import qualified Sound.MIDI.File.Save as MS
import qualified Sound.MIDI.File.Event as ME
import qualified Sound.MIDI.Message.Channel as MC
import qualified Sound.MIDI.Message.Channel.Voice as MCV
import qualified Sound.MIDI.Message.Channel.Mode as MCM
import qualified Data.EventList.Relative.TimeBody as E

type Frequency = Double
type Amplitude = Double
type Step = Double

maxFrequency :: Double
maxFrequency = 7000

data NoteLetter = A
        | Ais
        | B
        | C
        | Cis
        | D
        | Dis
        | E
        | F
        | Fis
        | G
        | Gis
        deriving (Show, Eq)

type Octave = Int

data Note = Note NoteLetter Octave deriving (Show, Eq)

getBaseFreq :: NoteLetter -> Frequency
getBaseFreq A   = 27.50
getBaseFreq Ais = 29.14
getBaseFreq B   = 30.87
getBaseFreq C   = 16.35
getBaseFreq Cis = 17.32
getBaseFreq D   = 18.35
getBaseFreq Dis = 19.45
getBaseFreq E   = 20.60
getBaseFreq F   = 21.83
getBaseFreq Fis = 23.12
getBaseFreq G   = 24.50
getBaseFreq Gis = 25.96

getFreq :: Note -> Frequency
getFreq (Note letter n) = getBaseFreq letter * (fromIntegral n + 1)

getNote :: Frequency -> Maybe Note
getNote freq
  | divisor A `mod'` 1 < 0.02 || divisor A `mod'` 1 > 0.98 = Just
    (Note A (findN freq A))
  | divisor Ais `mod'` 1 < 0.02 || divisor Ais `mod'` 1 > 0.98 = Just
    (Note Ais (findN freq Ais))
  | divisor B `mod'` 1 < 0.02 || divisor B `mod'` 1 > 0.98 = Just
    (Note B (findN freq B))
  | divisor C `mod'` 1 < 0.02 || divisor C `mod'` 1 > 0.98 = Just
    (Note C (findN freq C))
  | divisor Cis `mod'` 1 < 0.02 || divisor Cis `mod'` 1 > 0.98 = Just
    (Note Cis (findN freq Cis))
  | divisor D `mod'` 1 < 0.02 || divisor D `mod'` 1 > 0.98 = Just
    (Note D (findN freq D))
  | divisor Dis `mod'` 1 < 0.02 || divisor Dis `mod'` 1 > 0.98 = Just
    (Note Dis (findN freq Dis))
  | divisor E `mod'` 1 < 0.02 || divisor E `mod'` 1 > 0.98 = Just
    (Note E (findN freq E))
  | divisor F `mod'` 1 < 0.02 || divisor F `mod'` 1 > 0.98 = Just
    (Note F (findN freq F))
  | divisor Fis `mod'` 1 < 0.02 || divisor Fis `mod'` 1 > 0.98 = Just
    (Note Fis (findN freq Fis))
  | divisor G `mod'` 1 < 0.02 || divisor G `mod'` 1 > 0.98 = Just
    (Note G (findN freq G))
  | divisor Gis `mod'` 1 < 0.02 || divisor Gis `mod'` 1 > 0.98 = Just
    (Note Gis (findN freq Gis))
  | otherwise = Nothing
 where
  findN f n = round $ logBase 2 (f / getBaseFreq n) :: Int
  divisor note = logBase 2 (freq / getBaseFreq note)

-- samples = N = number of samples
-- n = sample number
hannWindow :: Integer -> Integer -> Double
hannWindow samples n =
  0.5 * (1 - cos (2 * pi * fromIntegral n / (fromIntegral samples - 1)))

sinFunc :: Double -> Double -> Double -> Double -> Double
sinFunc samples cycle amplitude p =
  -amplitude * sin (2 * pi * cycle * p / samples)

cosFunc :: Double -> Double -> Double -> Double -> Double
cosFunc samples cycle amplitude p =
  amplitude * cos (2 * pi * cycle * p / samples)

performFFT :: VS.Vector Double -> Int -> IO (VS.Vector (Complex Double))
performFFT vector intSamples = do
  let (vec, x1, x2) = VS.unsafeToForeignPtr $ VS.take intSamples vector
  carray <- unsafeForeignPtrToCArray vec (x1, x2)
  let fft           = dftRC carray
  let (len, fftptr) = toForeignPtr fft
  return (VS.unsafeFromForeignPtr fftptr 0 len)

-- |Filter an array of maxima by partioning fftabs into bins of 100 elements, then finding the
-- largest value in that bin.
-- peaks are the list of indexes of the maxima.
filterMaxima :: V.Vector Int -> V.Vector Double -> V.Vector Int
filterMaxima peaks fftabs
  | V.null peaks = V.empty
  | otherwise = V.filter (> -1) $ V.map
    (\i ->
      let
        section = V.map
          (\x -> (x, fftabs V.! x))
          (V.filter (\x -> x < V.length fftabs && x >= i && x < i + binSize)
                    peaks
          )
      in  if V.null section
            then -1
            else fst $ V.maximumBy (comparing snd)
            -- list of (index, magnitude)
                                                   section
    )
    (V.fromList [0, binSize .. V.last peaks])
  where binSize = 100

-- | Find the maxima from fftabs by finding all points where the left and right points
-- are less than the current point.
findMaxima :: V.Vector Double -> V.Vector Int
findMaxima fftabs = V.filter
  (\x -> (diffs' V.! x) > 0 && (diffs' V.! (x + 1)) < 0)
  (V.fromList [1 .. (V.length diffs' - 2)])
 where
  diff f x = f x - f (x - 1)
  diffs' = V.map (diff (fftabs V.!)) (V.fromList [1 .. (V.length fftabs - 1)])

-- | Find maxima by finding points where that point is the maximum in it's neighbourhood.
findMaxima2 :: Filterable a => Int -> V.Vector a -> V.Vector (Int, a)
findMaxima2 radius v =
  V.filter (\(i, x) -> x == V.maximum (neighbours i))
    $ V.imap (\i x -> (i, x)) v
  where neighbours x0 = V.ifilter (\i _ -> abs (x0 - i) <= radius) v

findSignificantNotesForSection
  :: Int -> Int -> V.Vector (Complex Double) -> [(Frequency, Double)]
findSignificantNotesForSection sampleRate samples fftSection = fmap
  (\cycle ->
    ( fromIntegral cycle * fromIntegral sampleRate / fromIntegral samples
    , fftabs V.! cycle
    )
  )
  maxima'
 where
  fftabs  = V.map magnitude fftSection :: V.Vector Double
  peaks   = findMaxima fftabs
  peaks'  = filterMaxima peaks fftabs
  maxima  = V.toList peaks' -- tricky: the cycles returned is the cycles+1
  maxima' = filter (\value -> fftabs V.! value > 20) maxima -- Find 'large' maxima

midiPitch :: Note -> Int
midiPitch (Note letter n) =
  round $ 12 * logBase 2 (((2 ^ n) * getBaseFreq letter) / 440) + 69

-- | Find notes that standout in each stft section.
findSignificantNotes
  :: Int -> Int -> [(Step, [Complex Double])] -> String -> IO ()
findSignificantNotes sampleRate samples stft filePrefix = do
  notes <-
    mapM
      (\(step, fft) -> do
        let freqs =
              findSignificantNotesForSection sampleRate samples (V.fromList fft)
        let notes = fmap (\(f, _) -> getNote f)
                         (filter (\(f, _) -> f < maxFrequency) freqs)
        return $ map (\note -> (step, note)) $ catMaybes notes
      )
      stft :: IO [[(Step, Note)]]

  print $ map (map (\(s, note) -> (s, note, midiPitch note))) notes

  let
    sections = map
      (foldl
        (\acc (step, note) -> E.append
          acc
          (E.singleton
            (toElapsedTime 0)
            (ME.MIDIEvent
              (MC.Cons
                (MC.toChannel 1)
                (MC.Voice
                  (MCV.NoteOn (MC.toPitch (midiPitch note)) (MC.toVelocity 20))
                )
              )
            )
          )
        )
        E.empty
      )
      notes

  let
    t =
      E.append
          (foldl (\acc section -> E.append acc (E.delay 100 section))
                 E.empty
                 sections
          )
        $ E.singleton
            (toElapsedTime 100)
            (ME.MIDIEvent (MC.Cons (MC.toChannel 1) (MC.Mode MCM.AllNotesOff)))

  let midi = Cons Serial (Ticks 120) [t] :: T

  MS.toFile (filePrefix ++ ".midi") midi

--    toFile def{_fo_format=SVG} "diff.svg" $ do
--        layout_title .= "First derivative"
--        plot (points "diff" $ zip fftx (VS.toList fftabs))

--    print $ V.map (\x -> (x, fftabs VS.! x)) peaks'
    --toFile def{_fo_format=SVG} "diff2.svg" $ do
    --    layout_title .= "Second derivative"
    --    plot (points "diff" $ zip fftx diffs'')

    --let sinAmplitudes = take 500 $ fmap (\(cycle, _) -> (fromIntegral cycle, fftimg   VS.! cycle / fromIntegral (samples `div` 2))) maxima :: [(Double, Double)]
    --let cosAmplitudes = take 500 $ fmap (\(cycle, _) -> (fromIntegral cycle, fftreals VS.! cycle / fromIntegral (samples `div` 2))) maxima :: [(Double, Double)]

    --print "sin amps"
    --print sinAmplitudes
    --print "cos amps"
    --print cosAmplitudes

    --toFile def{_fo_format=SVG} "sine2.svg" $ do
    --    layout_title .= "curve"
    --    plot (line "data" [zip x y])
    --    plot
    --        (line "fit"
    --        [zip
    --            x
    --            (fmap
    --                (\p -> foldr (\(cycle, amplitude) total -> total + sinFunc (fromIntegral samples) cycle amplitude p) 0 sinAmplitudes +
    --                    foldr (\(cycle, amplitude) total -> total + cosFunc (fromIntegral samples) cycle amplitude p) 0 cosAmplitudes) x)])

    --let signal = VS.take samples $ analytic_signal vector
    --let signalReal = VS.toList $ VS.map realPart signal
    --let signalImag = VS.toList $ VS.map imagPart signal
    ----print $ VS.length signal
    --toFile def{_fo_format=SVG} "signal.svg" $ do
    --    layout_title .= "Hibert transform derivative"
    --    plot (line "Hilbert transform" [zip x (VS.toList $ VS.map imagPart signal)])
    --    plot (line "Original signal" [zip x intList])
    --    plot (line "Analytic signal" [zip x $ VS.toList (VS.map magnitude signal)])

autoCorrelate :: Filterable a => VS.Vector a -> Int -> IO ()
autoCorrelate vector sampleRate = do
  let result = auto_correlation 1000 vector
  let y      = map realToFrac $ VS.toList result :: [Double]
  --let values = [zip (map (\x -> x / fromIntegral sampleRate) [1..]) y] :: [[(Double, Double)]]
  let values = [zip [1 ..] y] :: [[(Double, Double)]]
  toFile def { _fo_format = SVG } "autocorrelation.svg" $ do
    layout_title .= "auto correlation"
    plot (line "" values)

  let localMaxima = findMaxima2 10 (VS.convert result)
  print $ V.map (\(i, x) -> (i, realToFrac x)) localMaxima
  print $ V.imap
    (\j (i, _) -> if j == 0 then j else i - fst (localMaxima V.! (j - 1)))
    localMaxima
  print $ fromIntegral sampleRate / fromIntegral
    (V.maximum
      (V.imap
        (\j (i, _) -> if j == 0 then j else i - fst (localMaxima V.! (j - 1)))
        localMaxima
      )
    )

graphSignal :: Filterable a => VS.Vector a -> Int -> IO ()
graphSignal vector sampleRate = do
  let intList = map realToFrac $ VS.toList vector :: [Double]
  let y       = intList
  let x       = map (\x -> x / fromIntegral sampleRate) [1 .. 200] :: [Double]
  toFile def { _fo_format = SVG } "signal.svg" $ do
    layout_x_axis . laxis_title .= "Elapsed time (s)"
    layout_title .= "signal curve"
    plot (line "" [zip x y])

graphFFTSignal :: VS.Vector (Complex Double) -> Int -> Int -> IO ()
graphFFTSignal fftvec samples sampleRate = do
  let fftreals = VS.toList $ VS.map realPart fftvec :: [Double]
  let fftimg   = VS.toList $ VS.map imagPart fftvec :: [Double]
  let fftabs   = VS.toList $ VS.map magnitude fftvec :: [Double]
  let fftx =
        map (\x -> fromIntegral (x * sampleRate) / fromIntegral samples)
            [0 .. 1000] :: [Double]
  toFile def { _fo_format = SVG } "fft.svg" $ do
    layout_title .= "fft curve"
    layout_x_axis . laxis_title .= "Frequency (Hz)"
    plot (points "abs" $ zip fftx fftabs)
    plot (points "real" $ zip fftx fftreals)
    plot (points "imaginary" $ zip fftx fftimg)

data Options = Options {
    sourceFile :: String,
    stepWidth :: Int
}

args :: ParserInfo Options
args = info
  (Options <$> strOption (long "src" <> metavar "SOURCE") <*> option
    auto
    (long "step-width" <> metavar "STEP_WIDTH" <> value 10000 <> showDefault)
  )
  (fullDesc <> progDesc "Desc" <> header "test")

main :: IO ()
main = do
  opts         <- execParser args
  (i, content) <-
    SF.readFile (sourceFile opts) :: IO (SF.Info, Maybe (SFV.Buffer Double))
  let sampleRate = SF.samplerate i
  let samples    = SF.frames i * SF.channels i
  print samples
  print i
  let v = fmap SFV.fromBuffer content
  case v of
    Nothing     -> print "Nothing!"
    Just vector -> do
      graphSignal vector sampleRate

      -- Run FFT on whole dataset and plot the result
      fftvec <- performFFT vector samples
      graphFFTSignal fftvec samples sampleRate
      findSignificantNotes sampleRate samples [(0, VS.toList fftvec)] "fft"

      -- get auto correlation
      autoCorrelate vector sampleRate

      -- plot histogram
      let (vec, x1, x2) = VS.unsafeToForeignPtr vector
      carray <- unsafeForeignPtrToCArray vec (x1, x2)
      --let stepWidth = fromIntegral samples :: Integer
      let stepWidth = 5000 :: Integer
      print $ "Step width: " ++ show stepWidth
      let cutSamples =
            (samples `div` fromIntegral stepWidth) * fromIntegral stepWidth :: Int

      -- windows of size stepWidth
      let stft =
            fmap (stftSection carray cutSamples stepWidth)
                 [0, stepWidth `div` 2 .. fromIntegral cutSamples] :: [ ( Step
                , [Complex Double]
                )
              ]

      findSignificantNotes sampleRate cutSamples stft "stft"

      -- (window start, fft for window)
      let plots' =
            concatMap
              (spectrogramIntensity (fromIntegral cutSamples) sampleRate)
              stft :: [(Step, Frequency, Amplitude)]

      let stftData  = generateStftData (fromIntegral stepWidth) plots'
      let stftData' = filter (\((_, _), (_, _), z) -> z > 1) stftData
      toFile def { _fo_format = PNG } "stft.png" $ do
        layout_title .= "curve"
        plot
          (liftEC $ do
            plot_rect_title .= "title"
            plot_rect_style .= \amp ->
              let RGB a b c = hsl ((1 - amp) * 240) 1 (amp * 0.6)
              in  solidFillStyle (opaque $ sRGB a b c)
            plot_rect_values .= stftData'
          )

      return ()

-- convert to plottable set of values
generateStftData
  :: Step
  -> [(Step, Frequency, Amplitude)]
  -> [((Double, Double), (Double, Double), Amplitude)]
generateStftData windowSize = map (generateColumn windowSize)

-- column in spectrogram
generateColumn
  :: Step
  -> (Step, Frequency, Amplitude)
  -> ((Double, Double), (Double, Double), Amplitude)
generateColumn windowSize (step, freq, amp) =
  ((step, windowSize + step), (freq, freq + 1), amp)

spectrogramIntensity
  :: Double -> Int -> (Step, [Complex Double]) -> [(Step, Frequency, Amplitude)]
spectrogramIntensity totalSamples sampleRate (step, fft) =
  let rate = fromIntegral sampleRate / totalSamples
  in  filter
        (\(_, f, _) -> f < 4000)
        (map (\(cycle, x) -> (step, cycle * rate, magnitude (x + 0) ** 2))
             (zip [0 ..] fft)
        )

ftSignals :: Double -> [Complex Double] -> [Double] -> [Double]
ftSignals totalSamples fft = fmap
  (\x ->
    foldr
        (\(cycle, amplitude) total ->
          total + sinFunc totalSamples cycle amplitude x
        )
        0
        sinAmplitudes
      + foldr
          (\(cycle, amplitude) total ->
            total + cosFunc totalSamples cycle amplitude x
          )
          0
          cosAmplitudes
  )
 where
  significantSignals = filter (\x -> magnitude x > 5) fft
  cycles             = zip [0 ..] significantSignals
  sinAmplitudes =
    fmap
      (\(cycle, signal) ->
        (fromIntegral cycle, realPart signal / (totalSamples / 2.0))
      )
      cycles :: [(Double, Double)]
  cosAmplitudes =
    fmap
      (\(cycle, signal) ->
        (fromIntegral cycle, imagPart signal / (totalSamples / 2.0))
      )
      cycles :: [(Double, Double)]

-- | Apply window function to the set of samples
windowedSection
  :: CArray Int Double -> Int -> Integer -> Integer -> CArray Int Double
windowedSection data_ totalSamples n stepWidth = ixmapWithIndP
  (0, totalSamples)
  id
  (\_ x i' -> if i' < fromIntegral n || i' >= fromIntegral (n + stepWidth)
    then 0
    else x * hannWindow (fromIntegral stepWidth) (fromIntegral i' - n)
  )
  data_

-- | n = start of step sample number
stftSection
  :: CArray Int Double -> Int -> Integer -> Integer -> (Step, [Complex Double])
stftSection !completeSample totalSamples stepWidth n =
  (fromIntegral n, VS.toList $ VS.unsafeFromForeignPtr fftptr 0 len)
 where
  fft = dftRC (windowedSection completeSample totalSamples n stepWidth)
  (len, fftptr) = toForeignPtr fft

max_stars_repo_name: Deadleg/music-listener

max_issues_repo_name: Deadleg/music-listener

max_forks_repo_name: Deadleg/music-listener


text:
{-# LANGUAGE CPP #-}
{- MatInv.hs
 - 4x4 Matrix invert function.
 - Direct implementation obtained from Maple's CodeGeneration
 -
 - Timothy A. Chagnon
 - CS 636 - Spring 2009
 -}
module MatInv where

#ifdef USING_HMATRIX
import MathHmatrix
import Numeric.LinearAlgebra
inverse4f :: Mat4f -> Mat4f
inverse4f = inv

#else
import MathDirect

-- Analytical solution via Maple (see minv.mw)
inverse4f :: Mat4f -> Mat4f
inverse4f   (Mat4f  (Vec4f m11 m12 m13 m14)
                    (Vec4f m21 m22 m23 m24)
                    (Vec4f m31 m32 m33 m34)
                    (Vec4f m41 m42 m43 m44)) =
    let t1 = m22 * m43 in
    let t3 = m22 * m44 in
    let t5 = m42 * m24 in
    let t7 = m44 * m32 in
    let t9 = m43 * m32 in
    let t11 = m42 * m23 in
    let t14 = m11 * m22 in
    let t15 = m43 * m34 in
    let t17 = m44 * m33 in
    let t19 = m11 * m42 in
    let t20 = m24 * m33 in
    let t22 = m11 * m44 in
    let t23 = m32 * m23 in
    let t25 = m11 * m43 in
    let t26 = m32 * m24 in
    let t28 = m23 * m34 in
    let t30 = m31 * m14 in
    let t32 = m22 * m41 in
    let t33 = m13 * m34 in
    let t35 = m14 * m33 in
    let t37 = m31 * m13 in
    let t39 = m43 * m31 in
    let t40 = m12 * m24 in
    let t42 = m21 * m12 in
    let t44 = -t14 * t15 + t14 * t17 - t19 * t20 - t22 * t23 + t25 * t26 + t19 * t28 + t1 * t30 + t32 * t33 - t32 * t35 - t3 * t37 - t39 * t40 - t17 * t42 in
    let t46 = m42 * m21 in
    let t50 = m41 * m12 in
    let t52 = m21 * m13 in
    let t54 = m41 * m14 in
    let t57 = m44 * m31 in
    let t58 = m12 * m23 in
    let t60 = m41 * m13 in
    let t62 = m21 * m14 in
    let t65 = t5 * t37 - t46 * t33 + t46 * t35 + t15 * t42 + t50 * t20 + t7 * t52 + t54 * t23 - t50 * t28 + t57 * t58 - t60 * t26 - t9 * t62 - t11 * t30 in
    let t67 = 1.0 / (t44 + t65) in
    let t121 = m44 * m21 in
    let t125 = m43 * m21 in
    let t146 = m41 * m32 in
    let t152 = m42 * m31 in
    let t172 = m32 * m11 in
    let t174 = m32 * m21 in
    let t176 = m31 * m12 in
    let r00 = (-t1 * m34 + t3 * m33 - t5 * m33 - t7 * m23 + t9 * m24 + t11 * m34) * t67 in
    let r01 = -(t17 * m12 + m42 * m13 * m34 - m42 * m14 * m33 - t15 * m12 - t7 * m13 + t9 * m14) * t67 in
    let r02 = (t1 * m14 - t3 * m13 - m43 * m12 * m24 + t5 * m13 + m44 * m12 * m23 - t11 * m14) * t67 in
    let r03 = -(-m22 * m13 * m34 + m22 * m14 * m33 - t40 * m33 - m14 * m32 * m23 + t58 * m34 + m13 * m32 * m24) * t67 in
    let r10 = -(t17 * m21 + m41 * m23 * m34 - t15 * m21 - t57 * m23 + t39 * m24 - m41 * m24 * m33) * t67 in
    let r11 = (-t25 * m34 + t22 * m33 + t39 * m14 + t60 * m34 - t54 * m33 - t57 * m13) * t67 in
    let r12 = -(t22 * m23 - t25 * m24 - t121 * m13 - t54 * m23 + t60 * m24 + t125 * m14) * t67 in
    let r13 = (-m11 * m24 * m33 + m11 * m23 * m34 + m24 * m31 * m13 - t52 * m34 + t62 * m33 - m23 * m31 * m14) * t67 in
    let r20 = (t32 * m34 - t3 * m31 + t5 * m31 - t46 * m34 + t7 * m21 - t146 * m24) * t67 in
    let r21 = -(t22 * m32 - t19 * m34 + t152 * m14 + t50 * m34 - t54 * m32 - t57 * m12) * t67 in
    let r22 = (t22 * m22 - t121 * m12 - t54 * m22 - t19 * m24 + t46 * m14 + t50 * m24) * t67 in
    let r23 = -(m34 * m11 * m22 - m34 * m21 * m12 - t30 * m22 - t172 * m24 + t174 * m14 + t176 * m24) * t67 in
    let r30 = -(-t1 * m31 + t32 * m33 - t46 * m33 + t11 * m31 - t146 * m23 + t9 * m21) * t67 in
    let r31 = (-t19 * m33 + t25 * m32 - t39 * m12 + t152 * m13 + t50 * m33 - t60 * m32) * t67 in
    let r32 = -(t25 * m22 - t125 * m12 - t60 * m22 - t19 * m23 + t46 * m13 + t50 * m23) * t67 in
    let r33 = t67 * (m33 * m11 * m22 - m33 * m21 * m12 - t37 * m22 - t172 * m23 + t174 * m13 + t176 * m23) in
    Mat4f   (Vec4f r00 r01 r02 r03)
            (Vec4f r10 r11 r12 r13)
            (Vec4f r20 r21 r22 r23)
            (Vec4f r30 r31 r32 r33)

#endif


text:
module LispParser where

import Control.Monad (liftM)
import Control.Monad.Error.Class (throwError)
import Data.Complex (Complex (..))
import Data.Ratio ((%))
import LispVal
import Numeric (readFloat, readHex, readOct)
import Text.ParserCombinators.Parsec hiding (spaces)

readOrThrow :: Parser a -> String -> ThrowsError a
readOrThrow parser input = case parse parser "lisp" input of
    Left err  -> throwError $ Parser err
    Right val -> return val

readExpr :: String -> ThrowsError LispVal
readExpr = readOrThrow parseExpr

readExprList :: String -> ThrowsError [LispVal]
readExprList = readOrThrow (endBy parseExpr spaces)

parseLisp :: SourceName -> String -> Either ParseError LispVal
parseLisp = parse parseExpr

parseExpr :: Parser LispVal
parseExpr =
  parseAtom
    <|> parseString
    <|> parseCharacter
    <|> try parseFloat
    <|> try parseRational
    <|> try parseComplex
    <|> parseNumber
    <|> parseQuote
    <|> parseQuasiquote
    <|> do
      _ <- char '('
      x <- try parseList <|> parseDottedList
      _ <- char ')'
      return x

parseList :: Parser LispVal
parseList = liftM List $ sepBy parseExpr spaces

parseDottedList :: Parser LispVal
parseDottedList = do
  x <- endBy parseExpr spaces
  xs <- char '.' >> spaces >> parseExpr
  return $ DottedList x xs

parseQuote :: Parser LispVal
parseQuote = do
  _ <- char '\''
  x <- parseExpr
  return $ List [Atom "quote", x]

parseQuasiquote :: Parser LispVal
parseQuasiquote = do
  _ <- char ','
  x <- parseExpr
  return $ List [Atom "quasiquote", x]

parseAtom :: Parser LispVal
parseAtom = do
  first <- letter <|> symbol
  rest <- many (letter <|> digit <|> symbol)
  let atom = first : rest
  return $ case atom of
    "True" -> Bool True
    "False" -> Bool False
    _ -> Atom atom

parseNumber :: Parser LispVal
parseNumber = parseDecimal <|> parseHex <|> parseOct <|> parseDecimalNotation

parseDecimal :: Parser LispVal
parseDecimal = many1 digit >>= (return . Number . read)

parseDecimalNotation :: Parser LispVal
parseDecimalNotation = do
  _ <- try $ string "#d"
  x <- many1 digit
  (return . Number . read) x

parseHex :: Parser LispVal
parseHex = do
  _ <- try $ string "#x"
  x <- many1 hexDigit
  return $ Number (hexToDigit x)

hexToDigit :: (Eq a, Num a) => String -> a
hexToDigit x = fst $ head (readHex x)

parseOct :: Parser LispVal
parseOct = do
  _ <- try $ string "#o"
  x <- many1 octDigit
  return $ Number (octToDigit x)

octToDigit :: (Eq a, Num a) => String -> a
octToDigit x = fst $ head (readOct x)

parseString :: Parser LispVal
parseString = do
  _ <- char '"'
  x <- many $ noneOf "\"\\" <|> escapedChars
  _ <- char '"'
  return $ String x

parseFloat :: Parser LispVal
parseFloat = do
  whole <- many1 digit
  _ <- char '.'
  frac <- many1 digit
  return . Float . fst . head . readFloat $ whole ++ "." ++ frac

parseRational :: Parser LispVal
parseRational = do
  numerator <- many1 digit
  _ <- char '/'
  denominator <- many1 digit
  return $ Rational (read numerator % read denominator)

parseComplex :: Parser LispVal
parseComplex = do
  real <- many1 digit
  _ <- char '+'
  imaginary <- many1 digit
  _ <- char 'i'
  return $ Complex (read real :+ read imaginary)

parseCharacter :: Parser LispVal
parseCharacter = try (string "\\#") >>= (return . String)

escapedChars :: Parser Char
escapedChars = do
  _ <- char '\\'
  x <- oneOf "\\\"\'nrt"
  return $ case x of
    '\\' -> x
    '"' -> x
    '\'' -> x
    'n' -> '\n'
    't' -> '\t'
    'r' -> '\r'
    _ -> x

symbol :: Parser Char
symbol = oneOf "!$%&|*+-/:<=>?@^_~"

spaces :: Parser ()
spaces = skipMany1 space


text:
{-|
Module      : MachineLearning.TIR
Description : TIR expression data structures
Copyright   : (c) Fabricio Olivetti de Franca, 2022
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com
Stability   : experimental
Portability : POSIX

The TIR expression  represents a function of the form:

\[
f(x) = g(\sum_{i}{w_i \cdot t_i(\prod_{j}{x_j^{k_{ij}})}} / (1 + \sum_{i}{w_i \cdot t_i(\prod_{j}{x_j^{k_{ij}})}})
\]

with \(t_i\) being a transformation function, \(g\) an invertible function, \(w_i\) a linear coefficient, \(k_{ij}\)
the interaction strength.

Any given expression can be represented by two lists of terms, with each term
being composed of a transformatioon function and an interaction.
The transformation function is represented by a `Function` sum type.
The interaction is represented as a list of tuple of ints where the key is the
predictor index and the value is the strength of the predictor in this
term. Strengths with a value of zero are omitted.
-}
module MachineLearning.TIR where

import Control.Evolution
import Control.Monad.State.Strict
import Data.List.Split
import Data.SRTree
import System.Random
import Control.DeepSeq                           (NFData, rnf)
import Data.List                                 (delete)
import Data.SRTree.Print                         (showDefault)
import Data.Vector                               (Vector)
import qualified Data.Vector             as V
import qualified Data.Vector.Storable    as VS
import qualified Numeric.LinearAlgebra   as LA
import MachineLearning.Utils.Config

-- | `TIR` is a record type composed of the external
-- function `_funY` of type `Function`, numerator `_p`
-- and denominator `_q` of type `Sigma`
data TIR = TIR { _funY :: Function
               , _p :: Sigma
               , _q :: Sigma
               } deriving Show

instance NFData TIR where
  rnf _ = ()

-- | `Sigma` is just a list of terms `Pi`
type Sigma = [Pi]
-- | `Pi` is a triple composed of a coefficient, a `Function`
-- and a list of tuples where `(ix, k)` represents `x ! ix ^ k`.
type Pi    = (Double, Function, [(Int, Int)])

-- | generates a random integer within the specified range.
randomRng :: (Int, Int) -> Rnd Int
randomRng rng = state $ randomR rng
{-# INLINE randomRng #-}

-- | generates a random integer within the specified range excluding zero.
randomRngNZ :: (Int, Int) -> Rnd Int
randomRngNZ rng = do
  x <- randomRng rng
  if x == 0
    then randomRngNZ rng
    else pure x
{-# INLINE randomRngNZ #-}

-- | picks a random element from a list.
randomFrom :: [a] -> Rnd a
randomFrom xs = do
  ix <- randomRng (0, length xs - 1)
  pure (xs !! ix)
{-# INLINE randomFrom #-}

-- | returns a random index of variables provided by the mutation configuration.
randomVar :: MutationCfg -> Rnd (Maybe Int, MutationCfg)
randomVar params = do
  let vars = _vars params
      n    = length vars
  ix <- randomRng (0, n)
  if ix == n
     then pure (Nothing, params)
     else do let x = vars !! ix
             pure (Just x, params{ _vars=delete x vars })

-- | returns a list of random interactions (tuples of variables indeces and exponentes)
-- with parameters provided by the mutation configuration.
randomVars :: MutationCfg -> Rnd [(Int, Int)]
randomVars params = do
  (v, params') <- randomVar params
  k            <- randomRngNZ $ _kRange params
  case v of
    Nothing  -> pure []
    Just var -> do vs <- randomVars params'
                   pure $ (var, k) : vs

-- | returns a random `Pi`
randomPi :: MutationCfg -> Rnd (Maybe Pi)
randomPi params = do
  pis <- randomVars params
  f   <- randomFrom $ _funs params
  if null pis
    then pure Nothing
    else pure $ Just (1.0, f, pis)

-- | returns a random `Sigma`
randomSigma :: MutationCfg -> Int -> Rnd (Sigma, Int)
randomSigma params budget | budget <= 0 = pure ([], budget)
randomSigma params budget = do
  n <- randomRng (0, budget)
  if n == budget
     then pure ([], budget)
     else do term             <- randomPi params
             (terms, budget') <- randomSigma params (budget - spentBudget term)
             case term of
               Nothing -> pure (terms, budget')
               Just t  -> pure (t:terms, budget')

  where
    spentBudget Nothing           = 0
    spentBudget (Just (_, _, ps)) = 1 -- length ps

-- | returns a random `TIR` expression
randomTIR :: MutationCfg -> Rnd TIR
randomTIR params = do
  yf           <- randomFrom $ _yfuns params
  (p, budget') <- randomSigma params $ _budget params
  (q, _)       <- randomSigma params budget'
  if null p
    then randomTIR params
    else pure (TIR yf p q)

-- | We store thee dataset as a vector of columns.
-- Each vector is stored a `Storable`-based vector.
type Column a   = LA.Vector a

-- | A dataset is a `Vector` of `Column`
type Dataset a  = Vector (Column a)

-- | A constraint is a function that gets a symbolic tree
-- as an input and returns non negative `Double` representing
-- how much a constraint was violated.
type Constraint = SRTree Int Double -> Double

-- | An individual in the population is composed of
-- the chromossome, a vector of fitness, a list of
-- coefficients (for multiclass problems it stores
-- one vector of coefficient per class),
-- the constraint violation, the size of the expression,
-- and the penalty value.
data Individual = Individual { _chromo  :: TIR
                             , _fit     :: [Double]
                             , _weights :: [LA.Vector Double]
                             , _constr  :: Double
                             , _len     :: Int
                             , _penalty :: Double
                             }

-- | creates an unevaluated individual.
createIndividual :: TIR -> Individual
createIndividual tir = Individual tir [] [] 0.0 0 0.0

-- | calculates the penalized fitness.
penalizedFit :: Individual -> Double
penalizedFit t = (head . _fit) t + _penalty t
{-# INLINE penalizedFit #-}

-- | replaces the coefficients of a TIR expression
replaceConsts :: TIR -> V.Vector Double -> TIR
replaceConsts (TIR g p q) ws = TIR g p' q'
  where
    (p', ws1) = runState (traverse replaceWeight p) (V.toList ws)
    (q', ws2) = runState (traverse replaceWeight q) ws1

replaceWeight :: Pi -> State [Double] Pi
replaceWeight (w, g, h) = state $ \ws -> case ws of
                                           (wi:ws') -> ((wi, g, h), ws')
                                           []       -> error $ show h -- ((w, g, h), [])

instance Eq Individual where
    t1 == t2 = penalizedFit t1 == penalizedFit t2
instance Ord Individual where
    t1 <= t2 = penalizedFit t1 <= penalizedFit t2

instance NFData Individual where
  rnf _ = ()

instance Solution Individual where
  _getFitness = head . _fit
  _isFeasible = (<1e-12) . _constr

-- | creates a symbolic tree from a TIR expression.
assembleTree :: Double -> TIR -> SRTree Int Double
assembleTree bias (TIR f p q) = Fun f ((Const bias + assemble p) / (1 + assemble q))
  where
    -- assemble :: Sigma ix val -> SRTree ix val
    assemble []      = 0
    assemble [p']    = mk p'
    assemble (p':ps) = mk p' + assemble ps

    -- mk :: Pi ix val -> SRTree ix val
    mk (v, g, ts) = Const v * Fun g (foldr (\(ix, k) acc -> acc * Pow (Var ix) k) 1 ts)

-- | pretty print a solution.
prettyPrintsolution :: Individual -> String
prettyPrintsolution sol | Prelude.null (_fit sol) = error "unevaluated solution"
prettyPrintsolution sol = concat [ "Expression:\n", (showDefault . assembleTree bias . _chromo) sol, "\n"
                                 , "Fitness: ", (show . head . _fit) sol, "\n"
                                 , "Constraints: ", (show . _constr) sol, "\n"
                                 , "Length: ", (show . _len) sol, "\n"
                                 , "Penalty: ", (show . _penalty) sol, "\n"
                                 ]

  where bias = V.head $ VS.convert $ head $ _weights sol

max_stars_repo_name: folivetti/tir

max_issues_repo_name: folivetti/tir

max_forks_repo_name: folivetti/tir


text:
module Main where

import QuickSpec
import Test.QuickCheck
import Test.QuickCheck.Arbitrary
import Test.QuickCheck.Poly(OrdA(..))
import Data.List hiding ( insert )
import Data.Ord ( comparing )
import Data.Complex

{- Simple laws about addition -}
addSig :: [Sig]
addSig = [
  con "+" ((+) :: Int -> Int -> Int),
  con "0" (0 :: Int)
  ]

{- Simple laws about concat and map -}
concatMapSig :: [Sig]
concatMapSig = [
  con "map" (map :: (A -> B) -> [A] -> [B]),
  con "concat" (concat :: [[A]] -> [A]),
  con "concatMap" (concatMap :: (A -> [B]) -> [A] -> [B]),
  con "[]" ([] :: [A]),
  con ":" ((:) :: A -> [A] -> [A])
  ]

{- Laws about equivalence of two functions -}
rev :: [a] -> [a]
rev [] = []
rev (x:xs) = rev xs ++ [x]

revAcc :: [a] -> [a]
revAcc = aux []
  where
    aux acc [] = acc
    aux acc (x:xs) = aux (x:acc) xs

revSig :: [Sig]
revSig = [
  con "rev" (rev :: [A] -> [A]),
  con "revAcc" (revAcc :: [A] -> [A]),
  con "[]" ([] :: [A])
  ]

{- Laws about sortedness -}
sorted :: Ord a => [a] -> Bool
sorted [] = True
sorted [x] = True
sorted (x:y:xs) = x <= y && sorted (y:xs)

sortedSig :: [Sig]
sortedSig = [
  predicate "sorted" (sorted :: [Int] -> Bool),
  background [
    predicate "<" ((<) :: Int -> Int -> Bool),
    prelude
    ]
  ]

{- Many, many laws about foldr -}
foldrSig :: [Sig]
foldrSig = [
  con "foldr" (foldr :: (A -> B -> B) -> B -> [A] -> B),
  background [
    con "length" (length :: [A] -> Int),
    con "const" (const :: A -> B -> A),
    con "+" ((+) :: Int -> Int -> Int),
    con "1" (1 :: Int),
    con "0" (0 :: Int)
    ]
  ]

{- Custom datatypes similar to https://github.com/nick8325/quickspec/blob/master/examples/Heaps.hs -}
data Tree a = Leaf | Node (Tree a) a (Tree a) deriving Typeable

instance Eq a => Eq (Tree a) where
  t1 == t2 = inorder t1 == inorder t2

instance Ord a => Ord (Tree a) where
  compare = comparing inorder

instance (Ord a, Arbitrary a) => Arbitrary (Tree a) where
  arbitrary = fromList <$> arbitrary

inorder :: Tree a -> [a]
inorder Leaf = []
inorder (Node l v r) = inorder l ++ [v] ++ inorder r

fromList :: Ord a => [a] -> Tree a
fromList = foldr insert Leaf

insert :: Ord a => a -> Tree a -> Tree a
insert v Leaf = Node Leaf v Leaf
insert v (Node l vt r)
  | v <= vt = Node (insert v l) vt r
  | otherwise = Node l vt (insert v r)

mirror :: Tree a -> Tree a
mirror Leaf = Leaf
mirror (Node l v r) = Node (mirror r) v (mirror l)

{- Some laws on trees -}
treeSig1 :: [Sig]
treeSig1 = [
  monoType (Proxy :: Proxy OrdA),
  monoTypeWithVars ["t", "t1", "t2"] (Proxy :: Proxy (Tree OrdA)),
  monoTypeWithVars ["t", "t1", "t2"] (Proxy :: Proxy (Tree Int)),
  con "mirror" (mirror :: Tree OrdA -> Tree OrdA),
  con "rev" (rev :: [OrdA] -> [OrdA]),
  con "fromList" (fromList :: [OrdA] -> Tree OrdA),
  con "insert" (insert :: OrdA -> Tree OrdA -> Tree OrdA)
  ]

{- A real challenge... -}
treeSig2 :: [Sig]
treeSig2 = [
  monoTypeWithVars ["a", "b", "c"] (Proxy :: Proxy OrdA),
  monoTypeWithVars ["t", "t1", "t2"] (Proxy :: Proxy (Tree OrdA)),
  monoTypeWithVars ["t", "t1", "t2"] (Proxy :: Proxy (Tree Int)),
  con "fromList" (fromList :: [OrdA] -> Tree OrdA),
  con "inorder" (inorder :: Tree OrdA -> [OrdA]),
  predicate "sorted" (sorted :: [OrdA] -> Bool)
  ]

{- A "complex" example... observational equivalence -}
complexSig :: [Sig]
complexSig = [
  monoTypeObserve (Proxy :: Proxy Float),
  --monoType (Proxy :: Proxy (Complex Int)),
  monoTypeObserve (Proxy :: Proxy (Complex Float)),
  con "realPart" (realPart :: Complex Float -> Float),
  con "imagPart" (imagPart :: Complex Float -> Float),
  con ":+" ((:+) :: Float -> Float -> Complex Float),
  background [ prelude ]
  ]

main :: IO ()
main = return ()

max_stars_repo_name: phagenlocher/quickspec-examples

max_issues_repo_name: phagenlocher/quickspec-examples

max_forks_repo_name: phagenlocher/quickspec-examples


text:
{-# LANGUAGE BangPatterns #-}

{- |
Module      : Main
Description : Haskell translation of Neural Networks and Deep Learning code sample
Copyright   : (c) FUJII Satohi, 2016
License     : BSD3
Maintainer  : fujiisat@gmaill.com
Stability   : experimental
Portability : portable

My personal Haskell translation of Neural Networks and Deep Learning
code sample
(https://github.com/mnielsen/neural-networks-and-deep-learning,
Copyright (c) 2012-2015 Michael Nielsen)

-}

module Main where

import           Control.Monad         (replicateM, zipWithM)
import           Data.List             (foldl')
import qualified Data.Vector           as V
import           Numeric.LinearAlgebra (( #> ))
import qualified Numeric.LinearAlgebra as H
import           System.Environment    (getArgs)
import           System.Random         (newStdGen, randomIO)
import           System.Random.Shuffle (shuffle')
import           Text.Printf           (printf)

import           MnistData

type Vec = H.Vector H.R
type Mat = H.Matrix H.R
type Data = (Vec, Vec)
data Network = Network
               { numLayers :: !Int
               , sizes     :: V.Vector Int
               , biases    :: V.Vector Vec
               , weights   :: V.Vector Mat
               } deriving (Eq, Show)

main :: IO ()
main = do
  -- ([30, 10, 3.0], [784, 30, 10])
  ([epo, mbsize, e], networkArgs) <-  splitAt 3  <$> getArgs
  putStrLn $ "Epoch: " ++ epo
  putStrLn $ "Minibatch: " ++ mbsize
  putStrLn $ "Eta: " ++ e
  putStrLn $ "Network: " ++ show networkArgs
  (trainingData, validateDate, testData) <- loadDataWrapper
  network <- networkInit (map read networkArgs)
  !result <- sgd network (V.fromList trainingData) (read epo) (read mbsize) (read e) (Just $ V.fromList testData)
  -- print result
  return ()

networkInit :: [Int] -> IO Network
networkInit sizes = do
  b <- map H.flatten <$> mapM (`H.randn` 1) (tail sizes)
  w <- zipWithM H.randn sizes (tail sizes)
  return Network { numLayers = length sizes
                 , sizes     = V.fromList sizes
                 , biases    = V.fromList b
                 , weights   = V.fromList w
                 }

feedforward :: Network -> Vec -> Vec
feedforward (Network n s biases weights) a =
  V.foldl' (\ !a (w, b) -> sigmoidVec ((H.tr w #> a) + b) ) a $ V.zip weights biases

sgd :: Network
    -> V.Vector Data               -- Training data
    -> Int                         -- Epoch
    -> Int                         -- Minibatch size
    -> Double                      -- Eta
    -> Maybe (V.Vector (Vec, Int)) -- Test data
    -> IO Network
sgd network trainingData 0 miniBatchSize eta testData = return network
sgd network trainingData epohcs miniBatchSize eta testData = do
  miniBatches <- divide miniBatchSize <$> shuffleList trainingData
  let network' = V.foldl' (\ !n b -> updateMiniBatch n b eta) network miniBatches
  case testData of
    Just t  -> putStrLn $ printf "Epoch %d: %d / %d" epohcs (evaluate network' t) (length t)
    Nothing -> putStrLn $ printf "Epoch %d complete" epohcs
  sgd network' trainingData (epohcs - 1) miniBatchSize eta testData

updateMiniBatch :: Network -> V.Vector Data -> Double -> Network
updateMiniBatch network miniBatch eta = network { weights=weights', biases=biases'}
  where
    f :: (V.Vector Vec, V.Vector Mat) -> Data -> (V.Vector Vec, V.Vector Mat)
    f nabla mb = nablaNew
      where
        !deltaNabla = backprop network mb
        !nablaNew = ( V.zipWith (+) (fst nabla) (fst deltaNabla)
                    , V.zipWith (+) (snd nabla) (snd deltaNabla))

    (!nablaB, !nablaW) = V.foldl' f zeros miniBatch

    !zeros = ( V.map (\a -> a - a) $ biases network
             , V.map (\a -> a - a) $ weights network )

    !etaN = eta / fromIntegral (length miniBatch)
    !weights' = V.map (\(w, nw) -> w - H.cmap (etaN *) nw) $ V.zip (weights network) nablaW
    !biases'  = V.map (\(b, nb) -> b - H.cmap (etaN *) nb) $ V.zip (biases network) nablaB

backprop :: Network -> Data -> (V.Vector Vec, V.Vector Mat)
backprop network (x, y) = (V.fromList nablaB, V.fromList nablaW)
  where
    (ac0:ac1:activations, z:zs) = V.foldl' forward ([x], []) $ V.zip (biases network) (weights network)
--    !delta = costDerivative ac0 y * sigmoidPrimeVec z -- quadratic
    !delta = costDerivative ac0 y  -- cross entropy
    !nablaIni = ([delta], [H.tr (delta `H.outer` ac1)], delta)
    (nablaB, nablaW, _) = foldl' backward nablaIni
      $ zip3 zs (V.toList $ V.reverse $ weights network) activations

backward :: ([Vec], [Mat], Vec) -> (Vec, Mat, Vec) -> ([Vec], [Mat], Vec)
backward  (nabB, nabW, del) (z, w, ac) = (delta : nabB, H.tr (delta `H.outer` ac) : nabW, delta)
  where
    !sp = sigmoidPrimeVec z
    !delta = (w #> del) * sp

forward :: ([Vec], [Vec]) -> (Vec, Mat) -> ([Vec], [Vec])
forward (a@(ac:acs), zzs) (b, w) = (activation:a, z:zzs)
  where
    !z = (H.tr w #> ac) + b
    !activation = sigmoidVec z

costDerivative :: Vec -> Vec -> Vec
costDerivative outputActivations y = outputActivations - y

sigmoid :: Double -> Double
sigmoid z = 1.0 / (1.0 + exp (-z))

sigmoidVec :: Vec -> Vec
sigmoidVec = H.cmap sigmoid

sigmoidPrime :: Double -> Double
sigmoidPrime z = sigmoid z * (1 - sigmoid z)

sigmoidPrimeVec :: Vec -> Vec
sigmoidPrimeVec = H.cmap sigmoidPrime

evaluate :: Network -> V.Vector (Vec, Int) -> Int
evaluate network testData = length $ V.filter (uncurry (==)) testResults
  where
    testResults = V.map (\(x, y) -> (H.maxIndex (feedforward network x), y)) testData

divide :: Int -> V.Vector a -> V.Vector (V.Vector a)
divide n xs | V.null xs = V.empty
            | otherwise =  V.cons hd (divide n tl)
  where
    (hd, tl) = V.splitAt n xs

shuffleList :: V.Vector a -> IO (V.Vector a)
shuffleList xs = do
  let len = length xs
  gen <- newStdGen
  return $ V.fromList $ shuffle' (V.toList xs) len gen

max_stars_repo_name: fujiisat/nnadl-haskell-study

max_issues_repo_name: fujiisat/nnadl-haskell-study

max_forks_repo_name: fujiisat/nnadl-haskell-study