text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE LambdaCase          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving  #-}
{-# LANGUAGE TypeOperators       #-}

import Control.Monad
import Control.Monad.Random
import Data.List
import Data.Maybe
import Data.Singletons
import Data.Singletons.Prelude
import Data.Singletons.TypeLits
import Numeric.LinearAlgebra.Static
import System.Environment
import Text.Read

data Weights i o = W { wBiases :: !(R o)
                     , wNodes  :: !(L o i)
                     }                      -- an "o x i" layer

data Network :: Nat -> [Nat] -> Nat -> * where
    O     :: !(Weights i o)
          -> Network i '[] o
    (:&~) :: KnownNat h
          => !(Weights i h)
          -> !(Network h hs o)
          -> Network i (h ': hs) o
infixr 5 :&~

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 :: (KnownNat i, KnownNat o)
         => Weights i o
         -> R i
         -> R o
runLayer (W wB wN) v = wB + wN #> v

runNet :: (KnownNat i, KnownNat o)
       => Network i hs o
       -> R i
       -> R o
runNet = \case
   O w -> \(!v) -> logistic (runLayer w v)
   (w :&~ n') -> \(!v) -> let v' = logistic (runLayer w v)
                          in runNet n' v'

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

randomNet :: forall m i hs o. (MonadRandom m, KnownNat i, SingI hs, KnownNat o)
          => m (Network i hs o)
randomNet = go sing
  where
    go :: forall h hs'. KnownNat h
       => Sing hs'
       -> m (Network h hs' o)
    go = \case
        SNil            ->     O <$> randomWeights
        SNat `SCons` ss -> (:&~) <$> randomWeights <*> go ss

train :: forall i hs o. (KnownNat i, KnownNat o)
      => Double           -- ^ learning rate
      -> R i              -- ^ input vector
      -> R o              -- ^ target vector
      -> Network i hs o   -- ^ network to train
      -> Network i hs o
train rate x0 target = fst . go x0
  where
    go  :: forall j js. KnownNat j
        => R j              -- ^ input vector
        -> Network j js o   -- ^ network to train
        -> (Network j js o, R j)
    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 - konst rate * dEdy
              wN'  = wN - konst 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 - konst rate * dEdy
              wN'  = wN - konst 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 - 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 :: Network 2 '[16, 8] 1 <- randomNet
    let trained = foldl' trainEach net0 (zip inps outs)
          where
            trainEach :: (KnownNat i, SingI hs, KnownNat o)
                      => Network i hs o
                      -> (R i, R o)
                      -> Network i hs o
            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 :: KnownNat n => R n -> (R n, 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)

-- * Exercises
--

popLayer :: Network i (h ': hs) o -> (Weights i h, Network h hs o)
popLayer (w :&~ n) = (w, n)

addNetWeights :: forall i hs o. (KnownNat i, KnownNat o)
              => Sing hs
              -> Network i hs o
              -> Network i hs o
              -> Network i hs o
addNetWeights = \case
    SNil -> \case
      O w1 -> \case
        O w2      -> O (addW w1 w2)
    SNat `SCons` ss -> \case
      w1 :&~ n1 -> \case
        w2 :&~ n2 -> addW w1 w2 :&~ addNetWeights ss n1 n2
  where
    addW :: (KnownNat n, KnownNat m)
         => Weights n m
         -> Weights n m
         -> Weights n m
    addW (W b1 w1) (W b2 w2) = W (b1 + b2) (w1 + w2)

hiddenSing :: Network i hs o -> Sing hs
hiddenSing = \case
    O _      -> SNil
    _ :&~ n' -> SNat `SCons` hiddenSing n'

max_stars_repo_name: mstksg/inCode

max_issues_repo_name: mstksg/inCode

max_forks_repo_name: mstksg/inCode


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE EmptyDataDecls #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE ViewPatterns #-}
-- |
-- Module     : Data.Matrix.Symmetric.Mutable
-- Copyright  : Copyright (c) 2012 Aleksey Khudyakov <alexey.skladnoy@gmail.com>
-- License    : BSD3
-- Maintainer : Aleksey Khudyakov <alexey.skladnoy@gmail.com>
-- Stability  : experimental
--
-- Symmetric and hermitian matrices
module Data.Matrix.Symmetric.Mutable (
    -- * Data types
    MSymmetric
  , MHermitian
    -- ** Implementation
  , MSymmetricRaw(..)
  , IsSymmetric
  , IsHermitian
    -- * Function
  , new
  , symmetricAsDense
  , hermitianAsDense
  , symmIndex
    -- * Complex number
  , NumberType
  , IsReal
  , IsComplex
  , castSymmetric
  , Conjugate(..)
  ) where

import Control.Monad
import Control.Monad.Primitive

import Data.Complex  (Complex,conjugate)
import Data.Typeable (Typeable)
import           Data.Vector.Storable.Internal
import qualified Data.Vector.Generic.Mutable as M

-- import Foreign.Ptr
import Foreign.Marshal.Array ( advancePtr )
import Foreign.ForeignPtr
import Foreign.Storable

import Data.Internal
import Data.Vector.Storable.Strided.Mutable
import Data.Matrix.Generic.Mutable
import Data.Matrix.Dense.Mutable   (MMatrix(..))

import Unsafe.Coerce

----------------------------------------------------------------
-- Data types
----------------------------------------------------------------

-- | Symmetric/hermitian matrix. Whether it's symmetric of hermitian
--   is determined by type tag. See 'IsSymmetric' and 'IsHermitian'.
--
--   Storage takes n² elements and data is stored in column major
--   order. Fields are
--
-- * Order of matrix
--
-- * Leading dimension size
--
-- * Pointer to data
data MSymmetricRaw tag s a = MSymmetricRaw
                             {-# UNPACK #-} !Int -- Order of matrix
                             {-# UNPACK #-} !Int -- Leading dimension size
                             {-# UNPACK #-} !(ForeignPtr a)
                      deriving (Typeable)

-- | Type tag for symmetric matrices.
data IsSymmetric

-- | Type tag for hermitian matrices.
data IsHermitian

-- | Mutable symmetric matrix
type MSymmetric = MSymmetricRaw IsSymmetric

-- | Mutable hermitian matrix
type MHermitian = MSymmetricRaw IsHermitian

instance Storable a => IsMMatrix (MSymmetricRaw IsSymmetric) a where
  basicRows (MSymmetricRaw n _ _) = n
  {-# INLINE basicRows #-}
  basicCols (MSymmetricRaw _ n _) = n
  {-# INLINE basicCols #-}
  basicIsIndexMutable _ _     = True
  {-# INLINE basicIsIndexMutable #-}
  basicUnsafeRead  (MSymmetricRaw _ lda fp) (symmIndex -> (!i,!j))
    = unsafePrimToPrim
    $ withForeignPtr fp (`peekElemOff` (i + j*lda))
  {-# INLINE basicUnsafeRead #-}
  basicUnsafeWrite (MSymmetricRaw _ lda fp) (symmIndex -> (!i,!j)) x
    = unsafePrimToPrim
    $ withForeignPtr fp $ \p -> pokeElemOff p (i + j*lda) x
  {-# INLINE basicUnsafeWrite #-}
  basicCloneShape = new . rows
  {-# INLINE basicCloneShape #-}
  basicClone = cloneSym
  {-# INLINE basicClone #-}

instance (Conjugate a, Storable a) => IsMMatrix (MSymmetricRaw IsHermitian) a where
  basicRows (MSymmetricRaw n _ _) = n
  {-# INLINE basicRows #-}
  basicCols (MSymmetricRaw _ n _) = n
  {-# INLINE basicCols #-}
  basicIsIndexMutable _ _     = True
  {-# INLINE basicIsIndexMutable #-}
  basicUnsafeRead  (MSymmetricRaw _ lda fp) (!i,!j)
    = unsafePrimToPrim
    $ case () of
       _| i > j     -> conjugateNum `liftM`  withForeignPtr fp (`peekElemOff` (j + i*lda))
        | otherwise ->                       withForeignPtr fp (`peekElemOff` (i + j*lda))
  {-# INLINE basicUnsafeRead #-}
  basicUnsafeWrite (MSymmetricRaw _ lda fp) (!i,!j) x
    = unsafePrimToPrim
    $ case () of
       _| i > j     -> withForeignPtr fp $ \p -> pokeElemOff p (j + i*lda) (conjugateNum x)
        | otherwise -> withForeignPtr fp $ \p -> pokeElemOff p (i + j*lda) x
  {-# INLINE basicUnsafeWrite #-}
  basicCloneShape = new . rows
  {-# INLINE basicCloneShape #-}
  basicClone = cloneSym
  {-# INLINE basicClone #-}

-- | Choose index so upper part of matrix is accessed
symmIndex :: (Int,Int) -> (Int,Int)
{-# INLINE symmIndex #-}
symmIndex (i,j)
  | i > j     = (j,i)
  | otherwise = (i,j)

-- | Allocate new matrix. It works for both symmetric and hermitian
--   matrices.
new :: (PrimMonad m, Storable a)
    => Int                      -- ^ Matrix order
    -> m (MSymmetricRaw tag (PrimState m) a)
{-# INLINE new #-}
new n = do
  fp <- unsafePrimToPrim $ mallocVector $ n * n
  return $ MSymmetricRaw n n fp

-- | Convert to dense matrix. Dense matrix will use same buffer as
--   symmetric
symmetricAsDense
  :: (PrimMonad m, Storable a)
  => MSymmetricRaw IsSymmetric (PrimState m) a
  -> m (MMatrix (PrimState m) a)
symmetricAsDense (MSymmetricRaw n lda fp) = do
  let m = MMatrix n n lda fp
  forM_ [0 .. n-2] $ \i -> do
    let row = MVector (n-i-1) lda $ updPtr (`advancePtr` (  i + (i+1) * lda)) fp
        col = MVector (n-i-1) 1   $ updPtr (`advancePtr` (1+i +    i  * lda)) fp
    M.move col row
  return m

-- | Convert to dense matrix. Dense matrix will use same buffer as
--   symmetric
hermitianAsDense
  :: (PrimMonad m, Storable a, Conjugate a)
  => MSymmetricRaw IsSymmetric (PrimState m) a
  -> m (MMatrix (PrimState m) a)
hermitianAsDense (MSymmetricRaw n lda fp) = do
  let m = MMatrix n n lda fp
  forM_ [0 .. n-2] $ \i -> do
    let len = n - i - 1
        row = MVector (n-i-1) lda $ updPtr (`advancePtr` (  i + (i+1) * lda)) fp
        col = MVector (n-i-1) 1   $ updPtr (`advancePtr` (1+i +    i  * lda)) fp
    forM_ [0 .. len - 1] $ \j -> do
      M.unsafeWrite col j . conjugateNum =<< M.unsafeRead row j
  return m

----------------------------------------------------------------
-- Real/Complex distinction
----------------------------------------------------------------

-- | Type tag for real numbers
data IsReal

-- | Type tag for complex numbers
data IsComplex

type family NumberType a :: *

type instance NumberType Float       = IsReal
type instance NumberType Double      = IsReal
type instance NumberType (Complex a) = IsComplex

-- | Cast between symmetric and hermitian matrices is data parameter
--   is real.
castSymmetric :: (NumberType a ~ IsReal)
              => MSymmetricRaw tag s a -> MSymmetricRaw tag' s a
{-# INLINE castSymmetric #-}
castSymmetric = unsafeCoerce

-- | Conjugate which works for both real (noop) and complex values.
class Conjugate a where
  conjugateNum :: a -> a

instance Conjugate Float where
  conjugateNum = id
  {-# INLINE conjugateNum #-}

instance Conjugate Double where
  conjugateNum = id
  {-# INLINE conjugateNum #-}

instance RealFloat a => Conjugate (Complex a) where
  conjugateNum = conjugate
  {-# INLINE conjugateNum #-}

----------------------------------------------------------------
-- Helpers
----------------------------------------------------------------

-- Get n'th column of matrix as mutable vector. Internal since part of
-- vector contain junk
unsafeGetCol :: Storable a => MSymmetricRaw tag s a -> Int -> MVector s a
{-# INLINE unsafeGetCol #-}
unsafeGetCol (MSymmetricRaw n lda fp) i
  = MVector n 1 $ updPtr (`advancePtr` (i*lda)) fp

cloneSym :: (Storable a, PrimMonad m)
         => MSymmetricRaw tag (PrimState m) a
         -> m (MSymmetricRaw tag (PrimState m) a)
{-# INLINE cloneSym #-}
cloneSym m@(MSymmetricRaw n _ _) = do
  q <- new n
  forM_ [0 .. n - 1] $ \i ->
    M.unsafeCopy (unsafeGetCol q i) (unsafeGetCol m i)
  return q

ext: hs

max_stars_repo_name: Shimuuar/blas-lapack

max_issues_repo_name: Shimuuar/blas-lapack

max_forks_repo_name: Shimuuar/blas-lapack


text:
module Duck.Analysis where

import Duck.Types
import Duck.Plot
import Duck.Utils
import Duck.Statistics
import Duck.Hypothesis

import qualified Statistics.Resampling.Bootstrap as B

import Criterion
import Criterion.Main.Options (defaultConfig)
import Criterion.Internal (runAndAnalyseOne)
import Criterion.Monad (withConfig)
import Control.DeepSeq (NFData)
import qualified Criterion.Types as CT

-- |The worker function. Takes a set of parameters and
-- a function to evaluate and produces a full report of
-- results. It is required that the function to be evaluated
-- is of type (a -> b), where 'a' is sizable and willing to
-- have test cases generated, and 'b' needs to have a normal
-- form.
runSingleTest :: (Cased a, Sized a, NFData b) => Parameters -> (a -> b) -> IO FullReport
runSingleTest pr src = do
  putStrIf ((verbosity pr) >= Moderate) "Starting benchmarks..."
  testInstances <- genInstances $ genSizeList (range pr) (iterations pr)
  let sizes = [fromIntegral $ dimSize vl | vl <- testInstances]
  (expTimes, expStds, otlr) <- runSource pr (src) testInstances sizes
  putStrIf ((verbosity pr) >= Moderate) "Crunching results..."
  return (FullReport {experiment = expTimes,
                      expStd = expStds,
                      sizes = sizes,
                      outlierEffect = otlr / (fromIntegral $ length sizes)
                     })

-- |Takes an hypothesis complexity function and a full
-- report and fits the experimental data to the analytical
-- data to obtain a confidence value and a proportionality
-- constant (in seconds).
testHypothesis :: Hypothesis -> FullReport -> Report
testHypothesis hp fr = generateReport ana fr
  where ana = [hp (fromIntegral vl) | vl <- (sizes fr)]

-- |Performs a grouped hypothesis test, takes a full report
-- and a set of named hypothesis, as well as a grouping function
-- (some examples defined in Duck.Statistics) and produces a list
-- of reports.
testGroup :: FullReport -> [(String, Hypothesis)] ->
             ([(Double, (String, Report))] -> [(Double, (String, Report))])
             -> GroupReport
testGroup fr xs flt = map snd fltrd
  where ls = map (\(name, f) -> (name, testHypothesis f fr)) xs
        ind = zip (map (confidence . snd) ls) ls
        fltrd = flt ind

-- |Takes experimental data and an hypothesis and plots both.
plotHypothesis :: FullReport -> Hypothesis -> Report -> String -> IO ()
plotHypothesis fr hp r ttype = plotReportHypothesis sz exp ana ttype
  where sz = (sizes fr)
        exp = (experiment fr)
        ana = [(propConstant r) * hp (fromIntegral vl) | vl <- (sizes fr)]

-- |Runs a source function in a set of inputs with a given size
-- and produces a tuple of: ([Mean], [Standard deviation], Outlier effect),
-- wrapped in the IO monad.
-- A quiet verbosity produces no output, a moderate one produces a summary
-- and a full one details every step.
runSource :: (NFData b) => Parameters -> (a -> b) -> [a] -> [Int] -> IO ([Double], [Double], Double)
runSource _ _ [] [] = return ([], [], 0)
runSource pr src (inst:is) (sz:ss) = do
  let res = runAndAnalyseOne 1 "" (nf src inst)
  CT.Analysed rep <- withConfig (defaultConfig { CT.verbosity = CT.Quiet,
                                                 CT.timeLimit = (timePerTest pr) }) res
  let sampPoint = CT.anMean $ CT.reportAnalysis rep
      sampVar = CT.anStdDev $ CT.reportAnalysis rep
      sampOutlier = CT.anOutlierVar $ CT.reportAnalysis rep

  putStrIf ((verbosity pr) == Full) $ show sampPoint
  putStrIf ((verbosity pr) == Full) $ show sz
  putStrIf ((verbosity pr) == Full) $ show sampOutlier
  putStrIf ((verbosity pr) == Full) $ show sampVar
  putStrIf ((verbosity pr) == Full) ""
  putStrIf ((verbosity pr) == Moderate && (1 + length ss) `mod` 5 == 0) $
    (show $ 1 + length ss) ++ " iterations left (about " ++
    (show $ (round $ timePerTest pr) * 4 * (1 + length ss) `div` 3) ++ " seconds left)"

  (avgs, stds, otlr) <- runSource pr src is ss
  return ((B.estPoint sampPoint):avgs,
          (B.estPoint sampVar):stds,
          otlr + CT.ovFraction sampOutlier)

-- |Generates a report from a full report and a list
-- of analytical results.
generateReport :: [Double] -> FullReport -> Report
generateReport ana fr = Report {confidence = conf,
                                propConstant = const}
  where (const, conf) = regress ana (experiment fr)

max_stars_repo_name: gangsterveggies/duck-analyzer

max_issues_repo_name: gangsterveggies/duck-analyzer

max_forks_repo_name: gangsterveggies/duck-analyzer


text:
{-# LANGUAGE ScopedTypeVariables, TypeFamilies #-}

module Fuml.Supervised where

import qualified Fuml.Base.LinearRegression as LR
import qualified Fuml.Base.Logistic as LG
import Fuml.Core
import Numeric.LinearAlgebra
import Control.Monad.Identity

ols :: Supervisor Identity Double (Vector Double) Double
ols = Supervisor $ \_ theData ->
              let beta = LR.ols theData
              in return $ linPredict beta

olsWeighted :: Vector Double -> Supervisor Identity Double (Vector Double) Double
olsWeighted wvs
   = Supervisor $ \_ theData ->
              let beta = LR.wols wvs theData
              in return $ linPredict beta

ridge :: Double -> Supervisor Identity Double (Vector Double) Double
ridge gammac = Supervisor $ \_ xys ->
  let npars = size $ fst $ head xys
      gamma = diag $ konst gammac npars
      beta = LR.ridge gamma xys
  in return $ linPredict beta

linPredict :: Vector Double -> Predict (Vector Double) Double
linPredict beta = Predict beta (`dot` beta)

max_stars_repo_name: avctrh/filopodia

max_issues_repo_name: avctrh/filopodia

max_forks_repo_name: avctrh/filopodia


text:
{-# LANGUAGE BangPatterns #-}
module Main where
import           Control.Monad
import           Network
import           Numeric.LinearAlgebra
import           System.Directory      (getCurrentDirectory)
import           System.IO

netConfig :: IO NeuralNet
netConfig = initWB NN
    { weights   = []
    , biases    = []
    , eta       = 3.0
    , epochs    = 30
    , layers    = 2
    , layerSize = 30
    , batchSize = 10  -- has to be a divisor of 50000
    , trainData = []
    , testData  = []
    }

main :: IO ()
main = do
    dir <- getCurrentDirectory

    trainingData <- loadData
    testingData <- loadTestData
    net <- netConfig
    model <- train net {trainData = trainingData, testData = testingData}

    forM_ (enum $ weights model) $ \(i, m) -> do
        let path = (dir ++ "/data/result/w" ++ (show i) ++ ".txt") :: FilePath
        saveMatrix path "%g" m
        return ()
    forM_ (enum $ biases model) $ \(i, b) -> do
        let path = (dir ++ "/data/result/b" ++ (show i) ++ ".txt") :: FilePath
        saveMatrix path "%g" $ asColumn b
        return ()
    where
      enum :: [a] -> [(Int, a)]
      enum a = zip (enumFrom 0) a

max_stars_repo_name: WhateverHappns/haskell-mnist

max_issues_repo_name: WhateverHappns/haskell-mnist

max_forks_repo_name: WhateverHappns/haskell-mnist


text:
import qualified Statistics.Distribution.Normal as Stats
import qualified Amby as Am

main :: IO ()
main = do
  z <- Am.random Stats.standard 10000
  Am.save $ Am.rugPlot' z

max_stars_repo_name: jsermeno/amby

max_issues_repo_name: jsermeno/amby

max_forks_repo_name: jsermeno/amby


text:
module Mnist
    ( Image
    , Label
    , DataSet
    , loadMnist
    ) where

import Control.Monad
import Numeric.LinearAlgebra
import Network.HTTP.Simple (parseRequest, httpLBS, getResponseBody)
import System.Directory (doesFileExist)
import qualified Data.ByteString.Lazy as BL
import qualified Codec.Compression.GZip as GZ (compress, decompress)
import Data.Binary (encode, decode)

type Image = Matrix R
type Label = Vector R
type DataSet = (Image, Label)

baseUrl = "http://yann.lecun.com/exdb/mnist"
keyFiles = [
            ("train_img", "train-images-idx3-ubyte.gz"),
            ("train_label", "train-labels-idx1-ubyte.gz"),
            ("test_img", "t10k-images-idx3-ubyte.gz"),
            ("test_label", "t10k-labels-idx1-ubyte.gz")
          ]

assetsDir = "assets"
pickleFile = "mnist.dat"
imgSize = 784

generatePath :: String -> String
generatePath p = assetsDir ++ "/" ++ p

download :: String -> IO ()
download fn = do
    let p = generatePath fn

    e <- doesFileExist p
    unless e $ do
      putStrLn $ "Downloading " ++ fn ++ " ..."
      res <- httpLBS =<< parseRequest (baseUrl ++ "/" ++ fn)
      BL.writeFile p (getResponseBody res)
      putStrLn "Done"

downloadMnist :: [(String, String)] -> IO ()
downloadMnist [] = return ()
downloadMnist (x:xs) = do
    download $ snd x
    downloadMnist xs

toDoubleList :: BL.ByteString -> [Double]
toDoubleList = fmap (read . show . fromEnum) . BL.unpack

loadLabel :: String -> IO Label
loadLabel fn = do
    c <- fmap GZ.decompress (BL.readFile $ generatePath fn)
    return . vector . toDoubleList $ BL.drop 8 c

loadImg :: String -> IO Image
loadImg fn = do
    c <- fmap GZ.decompress (BL.readFile $ generatePath fn)
    return . matrix imgSize . toDoubleList $ BL.drop 16 c

convertDataset :: IO [DataSet]
convertDataset = do
    tri <- loadImg . snd . head $ keyFiles
    trl <- loadLabel . snd . (!!1) $ keyFiles
    ti <- loadImg . snd . (!!2) $ keyFiles
    tl <- loadLabel . snd . (!!3) $ keyFiles

    return [(tri, trl), (ti, tl)]

createPickle :: String -> [DataSet] -> IO ()
createPickle p ds = BL.writeFile p $ (GZ.compress . encode) ds

loadPickle :: String -> IO [DataSet]
loadPickle p = do
    eds <- BL.readFile p
    return $ (decode . GZ.decompress) eds

initMnist :: IO ()
initMnist = do
    downloadMnist keyFiles
    putStrLn "Creating binary DataSet file ..."
    createPickle (generatePath pickleFile) =<< convertDataset
    putStrLn "Done"

normalizeImg :: Bool -> [DataSet] -> IO [DataSet]
normalizeImg f ds@[tr, t]
    | f         = return [ ((/255) $ fst tr, snd tr), ((/255) $ fst t, snd t) ]
    | otherwise = return ds

loadMnist :: Bool -> IO [DataSet]
loadMnist nml = do
    let p = generatePath pickleFile
    e <- doesFileExist p
    unless e initMnist

    normalizeImg nml =<< loadPickle p


text:
module JuliaSet where

import Data.Complex
import Graphics.Image

-- This will work for any function with one complex number z as input.
-- Constants and coefficients must be set through partial application

-- in the output, x is the real portion of the result and y is the imaginary
-- portion, with these values laid out as one would expect with a cartesian
-- grid
visualizeJuliaSet :: (Complex Double -> Complex Double) -- f
                -> Double -- escapeRadius
                -> Int -- width in pixels
                -> Int -- height in pixels
                -> Int -- maxIter
                -> String -- output file name
                -> String -- array representation
                -> IO ()
visualizeJuliaSet f escapeRadius width height maxIter outputFilename arrType
  | arrType == "VU" = writeImage outputFilename
      (makeImageR VU (width, height) g)
  | arrType == "VS" = writeImage outputFilename
      (makeImageR VS (width, height) g)
  | arrType == "RSU" = writeImage outputFilename
      (makeImageR RSU (width, height) g)
  | arrType == "RPU" = writeImage outputFilename
      (makeImageR RPU (width, height) g)
  | arrType == "RSS" = writeImage outputFilename
      (makeImageR RSS (width, height) g)
  | arrType == "RPS" = writeImage outputFilename
      (makeImageR RPS (width, height) g)
  | otherwise = writeImage outputFilename (makeImageR RPU (width, height) g)
    where
      g = pixelToJuliaSetValue f escapeRadius width height maxIter

pixelToJuliaSetValue :: (Complex Double -> Complex Double) -- f
                     -> Double -- escapeRadius
                     -> Int -- width
                     -> Int -- height
                     -> Int -- maxIter
                     -> (Int, Int) -- i, j
                     -> Pixel RGB Double -- result: pixel value
pixelToJuliaSetValue f escapeRadius width height maxIter (i, j)
    = toPixelRGB (PixelHSI pixelValue 0.5 0.5)
      where
        pixelValue :: Double
        pixelValue = (fromIntegral (julia f (x :+ y) escapeRadius maxIter))/
                    (fromIntegral maxIter+1)
        x = -escapeRadius + (2*escapeRadius/(fromIntegral width)) *
          (fromIntegral j)
        y = -(-escapeRadius + (2*escapeRadius/(fromIntegral height)) *
          (fromIntegral i))

-- complex and real parts of z must be between -escapeRadius and escapeRadius
julia :: (Complex Double -> Complex Double) -- f
      -> Complex Double -- z
      -> Double -- escapeRadius
      -> Int -- maxIter
      -> Int -- result
julia f z escapeRadius maxIter = juliaRecursive f z escapeRadius 1 maxIter

-- not meant to be called by user; use julia instead
juliaRecursive :: (Complex Double -> Complex Double)
      -> Complex Double
      -> Double -- escapeRadius
      -> Int -- i
      -> Int -- maxIter
      -> Int -- result
juliaRecursive f z escapeRadius i maxIter
  | i > maxIter = maxIter
  | Data.Complex.magnitude(z') > escapeRadius = i
  | otherwise = juliaRecursive f z' escapeRadius (i+1) maxIter
    where
      z' = f z


text:
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE ConstraintKinds       #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedLabels      #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}

{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-|
Module      : Grenade.Layers.Convolution
Maintainter : Theo Charalambous
Description : Convolution layer with support for padding and bias

Convolutions with and without biases are supported, allowing for faster
implementation when we know there is no bias.
-}

module Grenade.Layers.Convolution where

import           Data.Function                       ((&))
import           Data.Kind                           (Type)
import           Data.List                           (foldl1')
import           Data.Maybe
import           Data.Proxy
import           Data.Serialize
import           Data.Singletons.TypeLits            hiding (natVal)
import           GHC.TypeLits
import           Control.DeepSeq
import           Lens.Micro                          ((^.))

import           Numeric.LinearAlgebra               hiding (R, konst,
                                                      uniformSample)
import qualified Numeric.LinearAlgebra               as LA
import           Numeric.LinearAlgebra.Static        hiding (build, toRows, (&),
                                                      (|||), size)

import           Grenade.Core
import           Grenade.Layers.Internal.Convolution
import           Grenade.Layers.Internal.Update
import           Grenade.Onnx
import           Grenade.Utils.LinearAlgebra
import           Grenade.Utils.ListStore

type OutputShapeIsOkay (input :: Nat) (pad :: Nat) (kernel :: Nat) (strides :: Nat) (output :: Nat)
  = ( strides * (output - 1) <= (input - kernel + pad)
    , (input - kernel + pad) <= (output * strides ) - 1 )

data HasBias = WithBias | WithoutBias

data ConvPadding = NoPadding | SameUpper | SameLower | Padding Nat Nat Nat Nat

data Convolution :: HasBias
                 -> ConvPadding
                 -> Nat -- Number of channels, for the first layer this could be RGB for instance.
                 -> Nat -- Number of filters, this is the number of channels output by the layer.
                 -> Nat -- The number of rows in the kernel filter
                 -> Nat -- The number of column in the kernel filter
                 -> Nat -- The row stride of the convolution filter
                 -> Nat -- The columns stride of the convolution filter
                 -> Type where
  Convolution :: ( KnownNat channels
                 , KnownNat filters
                 , KnownNat kernelRows
                 , KnownNat kernelColumns
                 , KnownNat strideRows
                 , KnownNat strideColumns
                 , KnownNat kernelFlattened
                 , kernelFlattened ~ (kernelRows * kernelColumns * channels))
              => !(L kernelFlattened filters) -- The kernel filter weights
              -> !(ListStore (L kernelFlattened filters)) -- The last kernel update (or momentum)
              -> Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns

  BiasConvolution :: ( KnownNat channels
                     , KnownNat filters
                     , KnownNat kernelRows
                     , KnownNat kernelColumns
                     , KnownNat strideRows
                     , KnownNat strideColumns
                     , KnownNat kernelFlattened
                     , kernelFlattened ~ (kernelRows * kernelColumns * channels))
                  => !(L kernelFlattened filters) -- The kernel filter weights
                  -> !(R filters) -- The bias weights
                  -> !(ListStore (L kernelFlattened filters)) -- The last kernel update (or momentum)
                  -> Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns

data Convolution' :: HasBias
                  -> ConvPadding
                  -> Nat -- Number of channels, for the first layer this could be RGB for instance.
                  -> Nat -- Number of filters, this is the number of channels output by the layer.
                  -> Nat -- The number of rows in the kernel filter
                  -> Nat -- The number of column in the kernel filter
                  -> Nat -- The row stride of the convolution filter
                  -> Nat -- The columns stride of the convolution filter
                  -> Type where
  Convolution' :: ( KnownNat channels
                  , KnownNat filters
                  , KnownNat kernelRows
                  , KnownNat kernelColumns
                  , KnownNat strideRows
                  , KnownNat strideColumns
                  , KnownNat kernelFlattened
                  , kernelFlattened ~ (kernelRows * kernelColumns * channels))
               => !(L kernelFlattened filters) -- The kernel filter weights
               -> Convolution' 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns

  BiasConvolution' :: ( KnownNat channels
                      , KnownNat filters
                      , KnownNat kernelRows
                      , KnownNat kernelColumns
                      , KnownNat strideRows
                      , KnownNat strideColumns
                      , KnownNat kernelFlattened
                      , kernelFlattened ~ (kernelRows * kernelColumns * channels))
                   => !(L kernelFlattened filters) -- The kernel filter weights
                   -> !(R filters) -- The bias weights
                   -> Convolution' 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns

{---------------------------}
{--    Layer instances    --}
{---------------------------}

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat ((kernelRows * kernelColumns) * channels)
         , KnownNat (filters * ((kernelRows * kernelColumns) * channels))
         ) => RandomLayer (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  createRandomWith m gen = do
    wN <- getRandomMatrix i i m gen
    wB <- getRandomVector i i m gen
    return $ BiasConvolution wN wB mkListStore
    where
      i = natVal (Proxy :: Proxy ((kernelRows * kernelColumns) * channels))

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat ((kernelRows * kernelColumns) * channels)
         , KnownNat (filters * ((kernelRows * kernelColumns) * channels))
         ) => RandomLayer (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  createRandomWith m gen = do
    wN <- getRandomMatrix i i m gen
    return $ Convolution wN mkListStore
    where
      i = natVal (Proxy :: Proxy ((kernelRows * kernelColumns) * channels))

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * channels)
         , kernelFlattened ~ (kernelRows * kernelColumns * channels)
         ) => UpdateLayer (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  type Gradient (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) = (Convolution' 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns)

  type MomentumStore (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns)  = ListStore (L (kernelRows * kernelColumns * channels) filters)

  runUpdate opt@OptSGD{} x@(Convolution oldKernel store) (Convolution' kernelGradient) =
    let  momentum = getData opt x store
         result = descendMatrix opt (MatrixValuesSGD oldKernel kernelGradient momentum)
         newStore = setData opt x store (matrixMomentum result)
    in Convolution (matrixActivations result) newStore
  runUpdate opt@OptAdam{} x@(Convolution oldKernel store) (Convolution' kernelGradient) =
    let (m, v) = toTuple $ getData opt x store
        result = descendMatrix opt (MatrixValuesAdam (getStep store) oldKernel kernelGradient m v)
        newStore = setData opt x store [matrixM result, matrixV result]
    in Convolution (matrixActivations result) newStore
    where toTuple [m ,v] = (m, v)
          toTuple xs = error $ "unexpected input of length " ++ show (length xs) ++ "in toTuple in Convolution.hs"

  reduceGradient grads = Convolution' $ dmmap (/ (fromIntegral $ length grads)) (foldl1' add (map (\(Convolution' x) -> x) grads))

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * channels)
         , kernelFlattened ~ (kernelRows * kernelColumns * channels)
         ) => UpdateLayer (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  type Gradient (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) = (Convolution' 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns)

  type MomentumStore (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns)  = ListStore (L (kernelRows * kernelColumns * channels) filters)

  runUpdate      = undefined
  reduceGradient = undefined

{-----------------------------------}
{--    No Bias Layer instances    --}
{-----------------------------------}

-- | A three dimensional image (or 2d with many channels) can have
--   an appropriately sized convolution filter run across it.
--   Case without bias vector.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , KnownNat channels
         , OutputShapeIsOkay inputRows 0 kernelRows strideRows outputRows
         , OutputShapeIsOkay inputCols 0 kernelCols strideCols outputCols
         ) => Layer (Convolution 'WithoutBias 'NoPadding channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) where

  type Tape (Convolution 'WithoutBias 'NoPadding channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) = S ('D3 inputRows inputCols channels)

  runForwards (Convolution kernel _) (S3D input) =
    let ex = extract input
        ek = extract kernel
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters)
        cs = fromIntegral $ natVal (Proxy :: Proxy channels)

        out  = forwardConv2d ex cs ix iy ek fs kx ky sx sy ox oy 0 0
    in  (S3D input, S3D . fromJust . create $ out)

  runBackwards (Convolution kernel _) (S3D input) (S3D dEdy) =
    let ex = extract input
        ek = extract kernel
        ey = extract dEdy
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters)
        cs = fromIntegral $ natVal (Proxy :: Proxy channels)

        (dx, dw)  = backwardConv2d ex cs ix iy ek fs kx ky sx sy 0 0 0 0 ey ox oy
    in  (Convolution' . fromJust . create $ dw, S3D . fromJust . create $ dx)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , KnownNat channels
         , KnownNat padt
         , KnownNat padl
         , KnownNat padb
         , KnownNat padr
         , OutputShapeIsOkay inputRows (padt + padb) kernelRows strideRows outputRows
         , OutputShapeIsOkay inputCols (padl + padr) kernelCols strideCols outputCols
         ) => Layer (Convolution 'WithoutBias ('Padding padl padt padr padb) channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) where

  type Tape (Convolution 'WithoutBias ('Padding padl padt padr padb) channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) = S ('D3 inputRows inputCols channels)

  runForwards (Convolution kernel _) (S3D input) =
    let ex = extract input
        ek = extract kernel
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters)
        cs = fromIntegral $ natVal (Proxy :: Proxy channels)
        pl = fromIntegral $ natVal (Proxy :: Proxy padl)
        pt = fromIntegral $ natVal (Proxy :: Proxy padt)

        out  = fromJust . create $ forwardConv2d ex cs ix iy ek fs kx ky sx sy ox oy pl pt
    in  (S3D input, S3D out)

  runBackwards (Convolution kernel _) (S3D input) (S3D dEdy) =
    let ex = extract input
        ek = extract kernel
        ey = extract dEdy
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows )
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols )
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters   )
        cs = fromIntegral $ natVal (Proxy :: Proxy channels  )
        pl = fromIntegral $ natVal (Proxy :: Proxy padl      )
        pt = fromIntegral $ natVal (Proxy :: Proxy padt      )
        pr = fromIntegral $ natVal (Proxy :: Proxy padr      )
        pb = fromIntegral $ natVal (Proxy :: Proxy padb      )

        (dx, dw) = backwardConv2d ex cs ix iy ek fs kx ky sx sy pl pt pr pb ey ox oy
    in  (Convolution' . fromJust . create $ dw, S3D . fromJust . create $ dx)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat channels
         ) => Layer (Convolution 'WithoutBias 'SameUpper channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 inputRows inputCols filters) where

  type Tape (Convolution 'WithoutBias 'SameUpper channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 inputRows inputCols filters) = S ('D3 inputRows inputCols channels)

  runForwards (Convolution kernel _) (S3D input) =
    let ex = extract input
        ek = extract kernel
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows )
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols )
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters   )
        cs = fromIntegral $ natVal (Proxy :: Proxy channels  )

        pady = (ix - 1) * sx + kx - ix
        padx = (iy - 1) * sy + ky - iy

        padt = div pady 2
        padl = div padx 2

        out  = fromJust . create $ forwardConv2d ex cs ix iy ek fs kx ky sx sy ix iy padl padt
    in  (S3D input, S3D out)

  runBackwards (Convolution kernel _) (S3D input) (S3D dEdy) =
    let ex = extract input
        ek = extract kernel
        ey = extract dEdy
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows )
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols )
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters   )
        cs = fromIntegral $ natVal (Proxy :: Proxy channels  )

        pady = (ix - 1) * sx + kx - ix
        padx = (iy - 1) * sy + ky - iy

        padl = div padx 2
        padt = div pady 2
        padr = padx - padl
        padb = pady - padt

        (dx, dw) = backwardConv2d ex cs ix iy ek fs kx ky sx sy padl padt padr padb ey ix iy
    in  (Convolution' . fromJust . create $ dw, S3D . fromJust . create $ dx)

-- | A three dimensional image (or 2d with many channels) can have
--   an appropriately sized convolution filter run across it.
--   Case without bias vector.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat channels
         ) => Layer (Convolution 'WithoutBias 'SameLower channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 inputRows inputCols filters) where

  type Tape (Convolution 'WithoutBias 'SameLower channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 inputRows inputCols filters) = S ('D3 inputRows inputCols channels)

  runForwards (Convolution kernel _) (S3D input) =
    let ex = extract input
        ek = extract kernel
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters)
        cs = fromIntegral $ natVal (Proxy :: Proxy channels)

        pady = (ix - 1) * sx + kx - ix
        padx = (iy - 1) * sy + ky - iy

        padt = pady - div pady 2
        padl = padx - div padx 2
        out  = fromJust . create $ forwardConv2d ex cs ix iy ek fs kx ky sx sy ix iy padl padt

    in  (S3D input, S3D out)

  runBackwards (Convolution kernel _) (S3D input) (S3D dEdy) =
    let ex = extract input
        ek = extract kernel
        ey = extract dEdy
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows )
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols )
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters   )
        cs = fromIntegral $ natVal (Proxy :: Proxy channels  )

        pady = (ix - 1) * sx + kx - ix
        padx = (iy - 1) * sy + ky - iy

        padt = pady - padb
        padl = padx - padr
        padr = div padx 2
        padb = div pady 2

        (dx, dw) = backwardConv2d ex cs ix iy ek fs kx ky sx sy padl padt padr padb ey ix iy
    in  (Convolution' . fromJust . create $ dw, S3D . fromJust . create $ dx)

{--------------------------------}
{--    Bias Layer instances    --}
{--------------------------------}

-- | A three dimensional image (or 2d with many channels) can have
--   an appropriately sized convolution filter run across it.
--   Case with bias vector
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , KnownNat channels
         , OutputShapeIsOkay inputRows 0 kernelRows strideRows outputRows
         , OutputShapeIsOkay inputCols 0 kernelCols strideCols outputCols
         ) => Layer (Convolution 'WithBias 'NoPadding channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) where

  type Tape (Convolution 'WithBias 'NoPadding channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) = S ('D3 inputRows inputCols channels)

  runForwards (BiasConvolution kernel biases _) (S3D input) =
    let ex = extract input
        ek = extract kernel
        eb = extract biases
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters)
        cs = fromIntegral $ natVal (Proxy :: Proxy channels)

        out  = fromJust . create $ forwardBiasConv2d ex cs ix iy eb ek fs kx ky sx sy ox oy 0 0

    in  (S3D input, S3D out)

  runBackwards (BiasConvolution kernel _ _) (S3D input) (S3D dEdy) =
    let ex = extract input
        ek = extract kernel
        ey = extract dEdy
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters)
        cs = fromIntegral $ natVal (Proxy :: Proxy channels)

        (dx', dw', db') = backwardBiasConv2d ex cs ix iy ek fs kx ky sx sy 0 0 0 0 ey ox oy
        dx              = fromJust . create $ dx'
        dw              = fromJust . create $ dw'
        db              = fromJust . create $ db'
    in  (BiasConvolution' dw db, S3D dx)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , KnownNat channels
         , KnownNat padt
         , KnownNat padl
         , KnownNat padb
         , KnownNat padr
         , OutputShapeIsOkay inputRows (padt + padb) kernelRows strideRows outputRows
         , OutputShapeIsOkay inputCols (padl + padr) kernelCols strideCols outputCols
         ) => Layer (Convolution 'WithBias ('Padding padl padt padr padb) channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) where

  type Tape (Convolution 'WithBias ('Padding padl padt padr padb) channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) = S ('D3 inputRows inputCols channels)

  runForwards (BiasConvolution kernel biases _) (S3D input) =
    let ex = extract input
        ek = extract kernel
        eb = extract biases
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows )
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols )
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters   )
        cs = fromIntegral $ natVal (Proxy :: Proxy channels  )
        pl = fromIntegral $ natVal (Proxy :: Proxy padl      )
        pt = fromIntegral $ natVal (Proxy :: Proxy padt      )

        out  = fromJust . create $ forwardBiasConv2d ex cs ix iy eb ek fs kx ky sx sy ox oy pl pt

    in  (S3D input, S3D out)

  runBackwards (BiasConvolution kernel _ _) (S3D input) (S3D dEdy) =
    let ex = extract input
        ek = extract kernel
        ey = extract dEdy
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows )
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols )
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters   )
        cs = fromIntegral $ natVal (Proxy :: Proxy channels  )
        pl = fromIntegral $ natVal (Proxy :: Proxy padl      )
        pt = fromIntegral $ natVal (Proxy :: Proxy padt      )
        pr = fromIntegral $ natVal (Proxy :: Proxy padr      )
        pb = fromIntegral $ natVal (Proxy :: Proxy padb      )

        (dx', dw', db') = backwardBiasConv2d ex cs ix iy ek fs kx ky sx sy pl pt pr pb ey ox oy
        dx              = fromJust . create $ dx'
        dw              = fromJust . create $ dw'
        db              = fromJust . create $ db'
    in  (BiasConvolution' dw db, S3D dx)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat channels
         ) => Layer (Convolution 'WithBias 'SameUpper channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 inputRows inputCols filters) where

  type Tape (Convolution 'WithBias 'SameUpper channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 inputRows inputCols filters) = S ('D3 inputRows inputCols channels)

  runForwards (BiasConvolution kernel bias _) (S3D input) =
    let ex = extract input
        ek = extract kernel
        eb = extract bias
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters)
        cs = fromIntegral $ natVal (Proxy :: Proxy channels)

        pady = (ix - 1) * sx + kx - ix
        padx = (iy - 1) * sy + ky - iy

        padt = div pady 2
        padl = div padx 2

        out  = fromJust . create $ forwardBiasConv2d ex cs ix iy eb ek fs kx ky sx sy ix iy padl padt

    in  (S3D input, S3D out)

  runBackwards (BiasConvolution kernel _ _) (S3D input) (S3D dEdy) =
    let ex = extract input
        ek = extract kernel
        ey = extract dEdy
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows )
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols )
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters   )
        cs = fromIntegral $ natVal (Proxy :: Proxy channels  )

        pady = (ix - 1) * sx + kx - ix
        padx = (iy - 1) * sy + ky - iy

        padl = div padx 2
        padt = div pady 2
        padr = padx - padl
        padb = pady - padt

        (dx', dw', db') = backwardBiasConv2d ex cs ix iy ek fs kx ky sx sy padl padt padr padb ey ix iy
        dx              = fromJust . create $ dx'
        dw              = fromJust . create $ dw'
        db              = fromJust . create $ db'
    in  (BiasConvolution' dw db, S3D dx)

-- | A three dimensional image (or 2d with many channels) can have
--   an appropriately sized convolution filter run across it.
--   Case without bias vector.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat channels
         ) => Layer (Convolution 'WithBias 'SameLower channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 inputRows inputCols filters) where

  type Tape (Convolution 'WithBias 'SameLower channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 inputRows inputCols filters) = S ('D3 inputRows inputCols channels)

  runForwards (BiasConvolution kernel bias _) (S3D input) =
    let ex = extract input
        ek = extract kernel
        eb = extract bias
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters)
        cs = fromIntegral $ natVal (Proxy :: Proxy channels)

        pady = (ix - 1) * sx + kx - ix
        padx = (iy - 1) * sy + ky - iy

        padt = pady - div pady 2
        padl = padx - div padx 2

        out  = fromJust . create $ forwardBiasConv2d ex cs ix iy eb ek fs kx ky sx sy ix iy padl padt
    in  (S3D input, S3D out)

  runBackwards (BiasConvolution kernel _ _) (S3D input) (S3D dEdy) =
    let ex = extract input
        ek = extract kernel
        ey = extract dEdy
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows )
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols )
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        fs = fromIntegral $ natVal (Proxy :: Proxy filters   )
        cs = fromIntegral $ natVal (Proxy :: Proxy channels  )

        pady = (ix - 1) * sx + kx - ix
        padx = (iy - 1) * sy + ky - iy

        padr = div padx 2
        padb = div pady 2
        padl = padx - padr
        padt = pady - padb

        (dx', dw', db') = backwardBiasConv2d ex cs ix iy ek fs kx ky sx sy padl padt padr padb ey ix iy
        dx              = fromJust . create $ dx'
        dw              = fromJust . create $ dw'
        db              = fromJust . create $ db'
    in  (BiasConvolution' dw db, S3D dx)

{---------------------------------}
{--    Other Layer instances    --}
{---------------------------------}

-- | A two dimentional image may have a convolution filter applied to it.
--   Case without bias vector
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , OutputShapeIsOkay inputRows 0 kernelRows strideRows outputRows
         , OutputShapeIsOkay inputCols 0 kernelCols strideCols outputCols
         ) => Layer (Convolution 'WithoutBias 'NoPadding 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 outputRows outputCols filters) where
  type Tape (Convolution 'WithoutBias 'NoPadding 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 outputRows outputCols filters) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    runForwards c (S3D input :: S ('D3 inputRows inputCols 1))

  runBackwards c tape grads =
    case runBackwards c tape grads of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) ->  (c', S2D back)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         ) => Layer (Convolution 'WithoutBias 'SameUpper 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 inputRows inputCols filters) where
  type Tape (Convolution 'WithoutBias 'SameUpper 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 inputRows inputCols filters) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    runForwards c (S3D input :: S ('D3 inputRows inputCols 1))

  runBackwards c tape grads =
    case runBackwards c tape grads of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) ->  (c', S2D back)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         ) => Layer (Convolution 'WithoutBias 'SameLower 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 inputRows inputCols filters) where
  type Tape (Convolution 'WithoutBias 'SameLower 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 inputRows inputCols filters) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    runForwards c (S3D input :: S ('D3 inputRows inputCols 1))

  runBackwards c tape grads =
    case runBackwards c tape grads of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) ->  (c', S2D back)

-- | A two dimentional image may have a convolution filter applied to it.
--   Case with bias vector.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , OutputShapeIsOkay inputRows 0 kernelRows strideRows outputRows
         , OutputShapeIsOkay inputCols 0 kernelCols strideCols outputCols
         ) => Layer (Convolution 'WithBias 'NoPadding 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 outputRows outputCols filters) where
  type Tape (Convolution 'WithBias 'NoPadding 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 outputRows outputCols filters) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    runForwards c (S3D input :: S ('D3 inputRows inputCols 1))

  runBackwards c tape grads =
    case runBackwards c tape grads of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) ->  (c', S2D back)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         ) => Layer (Convolution 'WithBias 'SameUpper 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 inputRows inputCols filters) where
  type Tape (Convolution 'WithBias 'SameUpper 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 inputRows inputCols filters) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    runForwards c (S3D input :: S ('D3 inputRows inputCols 1))

  runBackwards c tape grads =
    case runBackwards c tape grads of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) ->  (c', S2D back)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         ) => Layer (Convolution 'WithBias 'SameLower 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 inputRows inputCols filters) where
  type Tape (Convolution 'WithBias 'SameLower 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 inputRows inputCols filters) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    runForwards c (S3D input :: S ('D3 inputRows inputCols 1))

  runBackwards c tape grads =
    case runBackwards c tape grads of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) ->  (c', S2D back)

-- | A two dimensional image may have a convolution filter applied to it producing
--   a two dimensional image if both channels and filters is 1.
--   Case without bias vector.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputCols
         , KnownNat outputRows
         , OutputShapeIsOkay inputRows 0 kernelRows strideRows outputRows
         , OutputShapeIsOkay inputCols 0 kernelCols strideCols outputCols
         ) => Layer (Convolution 'WithoutBias 'NoPadding 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 outputRows outputCols) where
  type Tape (Convolution 'WithoutBias 'NoPadding 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 outputRows outputCols) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    case runForwards c (S3D input :: S ('D3 inputRows inputCols 1)) of
      (tps, S3D back :: S ('D3 outputRows outputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    case runBackwards c tape (S3D grads :: S ('D3 outputRows outputCols 1)) of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) -> (c', S2D back)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         ) => Layer (Convolution 'WithoutBias 'SameUpper 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 inputRows inputCols) where
  type Tape (Convolution 'WithoutBias 'SameUpper 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 inputRows inputCols) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    case runForwards c (S3D input :: S ('D3 inputRows inputCols 1)) of
      (tps, S3D back :: S ('D3 inputRows inputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    case runBackwards c tape (S3D grads :: S ('D3 inputRows inputCols 1)) of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) -> (c', S2D back)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         ) => Layer (Convolution 'WithoutBias 'SameLower 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 inputRows inputCols) where
  type Tape (Convolution 'WithoutBias 'SameLower 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 inputRows inputCols) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    case runForwards c (S3D input :: S ('D3 inputRows inputCols 1)) of
      (tps, S3D back :: S ('D3 inputRows inputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    case runBackwards c tape (S3D grads :: S ('D3 inputRows inputCols 1)) of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) -> (c', S2D back)

-- | A two dimensional image may have a convolution filter applied to it producing
--   a two dimensional image if both channels and filters is 1.
--   Case with bias vector.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputCols
         , KnownNat outputRows
         , OutputShapeIsOkay inputRows 0 kernelRows strideRows outputRows
         , OutputShapeIsOkay inputCols 0 kernelCols strideCols outputCols
         ) => Layer (Convolution 'WithBias 'NoPadding 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 outputRows outputCols) where
  type Tape (Convolution 'WithBias 'NoPadding 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 outputRows outputCols) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    case runForwards c (S3D input :: S ('D3 inputRows inputCols 1)) of
      (tps, S3D back :: S ('D3 outputRows outputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    case runBackwards c tape (S3D grads :: S ('D3 outputRows outputCols 1)) of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) -> (c', S2D back)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         ) => Layer (Convolution 'WithBias 'SameUpper 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 inputRows inputCols) where
  type Tape (Convolution 'WithBias 'SameUpper 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 inputRows inputCols) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    case runForwards c (S3D input :: S ('D3 inputRows inputCols 1)) of
      (tps, S3D back :: S ('D3 inputRows inputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    case runBackwards c tape (S3D grads :: S ('D3 inputRows inputCols 1)) of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) -> (c', S2D back)

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         ) => Layer (Convolution 'WithBias 'SameLower 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 inputRows inputCols) where
  type Tape (Convolution 'WithBias 'SameLower 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 inputRows inputCols) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    case runForwards c (S3D input :: S ('D3 inputRows inputCols 1)) of
      (tps, S3D back :: S ('D3 inputRows inputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    case runBackwards c tape (S3D grads :: S ('D3 inputRows inputCols 1)) of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) -> (c', S2D back)

-- | A three dimensional image can produce a 2D image from a convolution with 1 filter
--   Case without bias vector.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputCols
         , KnownNat outputRows
         , KnownNat channels
         , OutputShapeIsOkay inputRows 0 kernelRows strideRows outputRows
         , OutputShapeIsOkay inputCols 0 kernelCols strideCols outputCols
         ) => Layer (Convolution 'WithoutBias 'NoPadding channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 outputRows outputCols) where
  type Tape (Convolution 'WithoutBias 'NoPadding channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 outputRows outputCols) = S ('D3 inputRows inputCols channels)
  runForwards c input =
    case runForwards c input of
      (tps, S3D back :: S ('D3 outputRows outputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    runBackwards c tape (S3D grads :: S ('D3 outputRows outputCols 1))

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat channels
         ) => Layer (Convolution 'WithoutBias 'SameUpper channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 inputRows inputCols) where
  type Tape (Convolution 'WithoutBias 'SameUpper channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 inputRows inputCols) = S ('D3 inputRows inputCols channels)
  runForwards c input =
    case runForwards c input of
      (tps, S3D back :: S ('D3 inputRows inputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    runBackwards c tape (S3D grads :: S ('D3 inputRows inputCols 1))

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat channels
         ) => Layer (Convolution 'WithoutBias 'SameLower channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 inputRows inputCols) where
  type Tape (Convolution 'WithoutBias 'SameLower channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 inputRows inputCols) = S ('D3 inputRows inputCols channels)
  runForwards c input =
    case runForwards c input of
      (tps, S3D back :: S ('D3 inputRows inputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    runBackwards c tape (S3D grads :: S ('D3 inputRows inputCols 1))

-- | A three dimensional image can produce a 2D image from a convolution with 1 filter
--   Case with bias vector.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputCols
         , KnownNat outputRows
         , KnownNat channels
         , OutputShapeIsOkay inputRows 0 kernelRows strideRows outputRows
         , OutputShapeIsOkay inputCols 0 kernelCols strideCols outputCols
         ) => Layer (Convolution 'WithBias 'NoPadding channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 outputRows outputCols) where
  type Tape (Convolution 'WithBias 'NoPadding channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 outputRows outputCols) = S ('D3 inputRows inputCols channels)
  runForwards c input =
    case runForwards c input of
      (tps, S3D back :: S ('D3 outputRows outputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    runBackwards c tape (S3D grads :: S ('D3 outputRows outputCols 1))

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat channels
         ) => Layer (Convolution 'WithBias 'SameUpper channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 inputRows inputCols) where
  type Tape (Convolution 'WithBias 'SameUpper channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 inputRows inputCols) = S ('D3 inputRows inputCols channels)
  runForwards c input =
    case runForwards c input of
      (tps, S3D back :: S ('D3 inputRows inputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    runBackwards c tape (S3D grads :: S ('D3 inputRows inputCols 1))

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat channels
         ) => Layer (Convolution 'WithBias 'SameLower channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 inputRows inputCols) where
  type Tape (Convolution 'WithBias 'SameLower channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 inputRows inputCols) = S ('D3 inputRows inputCols channels)
  runForwards c input =
    case runForwards c input of
      (tps, S3D back :: S ('D3 inputRows inputCols 1)) ->  (tps, S2D back)

  runBackwards c tape (S2D grads) =
    runBackwards c tape (S3D grads :: S ('D3 inputRows inputCols 1))

{--------------------}
{-- ONNX Instances --}
{--------------------}

instance OnnxOperator (Convolution 'WithBias padding channels filters kernelRows kernelCols strideRows strideCols) where
  onnxOpTypeNames _ = ["Conv"]

instance OnnxOperator (Convolution 'WithoutBias padding channels filters kernelRows kernelCols strideRows strideCols) where
  onnxOpTypeNames _ = ["Conv"]

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat channels
         , KnownNat (kernelRows * kernelCols * channels)
         ) => OnnxLoadable (Convolution 'WithoutBias 'NoPadding channels filters kernelRows kernelCols strideRows strideCols) where
  loadOnnxNode inits node = do
    node `doesNotHaveAttribute` "auto_pad"

    node & hasSupportedDilations
    node & hasSupportedGroup

    (node `hasMatchingShape` "kernel_shape") kernelShape
    (node `hasMatchingShape` "strides"     ) strideShape
    hasCorrectPadding node (Proxy :: Proxy 0) (Proxy :: Proxy 0) (Proxy :: Proxy 0) (Proxy :: Proxy 0)

    case node ^. #input of
      [_, w] -> do
        filterWeights <- tr <$> readInitializerTensorIntoMatrix inits w
        return (Convolution filterWeights mkListStore)
      _ -> onnxIncorrectNumberOfInputs
      where
        kernelShape = [natVal (Proxy :: Proxy kernelRows), natVal (Proxy :: Proxy kernelCols)]
        strideShape = [natVal (Proxy :: Proxy strideRows), natVal (Proxy :: Proxy strideCols)]

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat channels
         , KnownNat (kernelRows * kernelCols * channels)
         ) => OnnxLoadable (Convolution 'WithoutBias 'SameUpper channels filters kernelRows kernelCols strideRows strideCols) where
  loadOnnxNode inits node = do
    node `doesNotHaveAttribute` "auto_pad"

    node & hasSupportedDilations
    node & hasSupportedGroup

    (node `hasMatchingShape` "kernel_shape") kernelShape
    (node `hasMatchingShape` "strides"     ) strideShape

    case node ^. #input of
      [_, w] -> do
        filterWeights <- tr <$> readInitializerTensorIntoMatrix inits w
        return (Convolution filterWeights mkListStore)
      _ -> onnxIncorrectNumberOfInputs
      where
        kernelShape = [natVal (Proxy :: Proxy kernelRows), natVal (Proxy :: Proxy kernelCols)]
        strideShape = [natVal (Proxy :: Proxy strideRows), natVal (Proxy :: Proxy strideCols)]

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat channels
         , KnownNat padl, KnownNat padt, KnownNat padr, KnownNat padb
         ) => OnnxLoadable (Convolution 'WithoutBias ('Padding padl padt padr padb) channels filters kernelRows kernelCols strideRows strideCols) where
  loadOnnxNode inits node = do
    node `doesNotHaveAttribute` "auto_pad"

    node & hasSupportedDilations
    node & hasSupportedGroup

    (node `hasMatchingShape` "kernel_shape") kernelShape
    (node `hasMatchingShape` "strides"     ) strideShape

    hasCorrectPadding node (Proxy :: Proxy padl) (Proxy :: Proxy padr) (Proxy :: Proxy padt) (Proxy :: Proxy padb)

    case node ^. #input of
      [_, w] -> do
        filterWeights <- tr <$> readInitializerTensorIntoMatrix inits w
        return (Convolution filterWeights mkListStore)
      _ -> onnxIncorrectNumberOfInputs
      where
        kernelShape = [natVal (Proxy :: Proxy kernelRows), natVal (Proxy :: Proxy kernelCols)]
        strideShape = [natVal (Proxy :: Proxy strideRows), natVal (Proxy :: Proxy strideCols)]

instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat channels
         , KnownNat (kernelRows * kernelCols * channels)
         ) => OnnxLoadable (Convolution 'WithBias padding channels filters kernelRows kernelCols strideRows strideCols) where
  loadOnnxNode inits node = do
    node & hasSupportedDilations
    node & hasSupportedGroup

    (node `hasMatchingShape` "kernel_shape") kernelShape
    (node `hasMatchingShape` "strides"     ) strideShape

    -- todo: proper checking to see if auto_pad attribute is valid

    case node ^. #input of
      [_, w, b] -> do
        filterWeights <- tr <$> readInitializerTensorIntoMatrix inits w
        filterBias    <- readInitializerVector inits b
        return (BiasConvolution filterWeights filterBias mkListStore)
      _ -> onnxIncorrectNumberOfInputs
      where
        kernelShape = [natVal (Proxy :: Proxy kernelRows), natVal (Proxy :: Proxy kernelCols)]
        strideShape = [natVal (Proxy :: Proxy strideRows), natVal (Proxy :: Proxy strideCols)]

{--------------------}
{-- Misc Instances --}
{--------------------}

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         ) =>
         LayerOptimizerData (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'SGD) where
  type MomentumExpOptResult (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'SGD) = L (kernelRows * kernelColumns * channels) filters
  type MomentumDataType (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'SGD) = L (kernelRows * kernelColumns * channels) filters
  getData opt x store = head $ getListStore opt x store
  setData opt x store = setListStore opt x store . return
  newData _ _ = konst 0

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         ) =>
         LayerOptimizerData (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'SGD) where
  type MomentumExpOptResult (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'SGD) = L (kernelRows * kernelColumns * channels) filters
  type MomentumDataType (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'SGD) = L (kernelRows * kernelColumns * channels) filters
  getData = undefined
  setData = undefined
  newData = undefined

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         ) => LayerOptimizerData (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'Adam) where
  type MomentumExpOptResult (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'Adam)  = [L (kernelRows * kernelColumns * channels) filters]
  type MomentumDataType (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'Adam) = L (kernelRows * kernelColumns * channels) filters
  getData = getListStore
  setData = setListStore
  newData _ _ = konst 0

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * channels)
         ) => LayerOptimizerData (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'Adam) where
  type MomentumExpOptResult (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'Adam)  = [L (kernelRows * kernelColumns * channels) filters]
  type MomentumDataType (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) (Optimizer 'Adam) = L (kernelRows * kernelColumns * channels) filters
  getData = undefined
  setData = undefined
  newData = undefined

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns) => FoldableGradient (Convolution' hasBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  mapGradient f (Convolution' kernelGradient) = Convolution' (dmmap f kernelGradient)
  mapGradient f (BiasConvolution' kernelGradient biasGradient) = BiasConvolution' (dmmap f kernelGradient) (dvmap f biasGradient)
  squaredSums (Convolution' kernelGradient) = [sumM . squareM $ kernelGradient]
  squaredSums (BiasConvolution' kernelGradient biasGradient) = [sumM . squareM $ kernelGradient, sumV . squareV $ biasGradient ]

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * channels)
         ) => Serialize (Convolution 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  put (Convolution w store) = do
    putListOf put . toList . flatten . extract $ w
    put (fmap (toList . flatten . extract) store)
  get = do
      let f  = fromIntegral $ natVal (Proxy :: Proxy filters)
      wN    <- maybe (fail "Vector of incorrect size") return . create . reshape f . LA.fromList =<< getListOf get
      store <- fmap (fromMaybe (error "Vector of incorrect size") . create . reshape f . LA.fromList)  <$> get
      return $ Convolution wN store

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * channels)
         ) => Serialize (Convolution 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  put (BiasConvolution w b store) = do
    putListOf put . toList . flatten . extract $ w
    putListOf put . toList . extract $ b
    put (fmap (toList . flatten . extract) store)
  get = do
      let f  = fromIntegral $ natVal (Proxy :: Proxy filters)
      wN    <- maybe (fail "Vector of incorrect size") return . create . reshape f . LA.fromList =<< getListOf get
      wB    <- maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get
      store <- fmap (fromMaybe (error "Vector of incorrect size") . create . reshape f . LA.fromList)  <$> get
      return $ BiasConvolution wN wB store

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * channels)
         ) => Serialize (Convolution' 'WithoutBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  put (Convolution' w) = do
    putListOf put . toList . flatten . extract $ w
  get = do
      let f  = fromIntegral $ natVal (Proxy :: Proxy filters)
      wN    <- maybe (fail "Vector of incorrect size") return . create . reshape f . LA.fromList =<< getListOf get
      return $ Convolution' wN

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * channels)
         ) => Serialize (Convolution' 'WithBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  put (BiasConvolution' w b) = do
    putListOf put . toList . flatten . extract $ w
    putListOf put . toList . extract $ b
  get = do
      let f  = fromIntegral $ natVal (Proxy :: Proxy filters)
      wN    <- maybe (fail "Vector of incorrect size") return . create . reshape f . LA.fromList =<< getListOf get
      wB    <- maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get
      return $ BiasConvolution' wN wB

instance NFData (Convolution hasBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  rnf (BiasConvolution a b c) = rnf a `seq` rnf b `seq` rnf c
  rnf (Convolution a b)       = rnf a `seq` rnf b

instance NFData (Convolution' hasBias padding channels filters kernelRows kernelColumns strideRows strideColumns) where
  rnf (BiasConvolution' a b) = rnf a `seq` rnf b
  rnf (Convolution' a)       = rnf a

instance Show (Convolution hasBias padding c f k k' s s') where
  show (BiasConvolution _ _ _) = "Bias Convolution"
  show (Convolution _ _)       = "Convolution"

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
module Parser where
import Data.Binary
import Data.Binary.Get
import qualified Data.ByteString.Lazy as BS
import Codec.Compression.GZip (decompress)
import Control.Monad
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Devel
import qualified Data.Vector as V
import Control.Monad.ST

type Pixel = Word8
type Image = V.Vector (Matrix Float)
type Label = Vector Float

decodeImages :: Get [Image]
decodeImages = do
    mc <- getWord32be
    guard (mc == 0x00000803)
    [d1,d2,d3] <- many 3 getWord32be
    guard (d2 == 28 && d3 == 28)
    many d1 pic
  where
    pic :: Get Image
    pic = do
      bs <- getByteString (28*28)
      return . V.singleton . reshape 28 . toVecDouble . fromByteString $ bs
    toVecDouble :: Vector Pixel -> Vector Float
    -- mapVectorM requires the monad be strict
    -- so Identity monad shall not be used
    toVecDouble v = runST $ mapVectorM (return . (/255) . fromIntegral) v

decodeLabels :: Get [Label]
decodeLabels = do
    mc <- getWord32be
    guard (mc == 0x00000801)
    d1 <- getWord32be
    many d1 lbl
  where
    lbl :: Get Label
    lbl = do
      v <- fromIntegral <$> (get :: Get Word8)
      return $ fromList (replicate v 0 ++ [1] ++ replicate (9-v) 0)

many :: (Integral n, Monad m) => n -> m a -> m [a]
many cnt dec = sequence (replicate (fromIntegral cnt) dec)

trainingData :: IO ([Image], [Label])
trainingData = do
    s <- decompress <$> BS.readFile "tdata/train-images-idx3-ubyte.gz"
    t <- decompress <$> BS.readFile "tdata/train-labels-idx1-ubyte.gz"
    return (runGet decodeImages s, runGet decodeLabels t)

testData :: IO ([Image], [Label])
testData = do
    s <- decompress <$> BS.readFile "tdata/t10k-images-idx3-ubyte.gz"
    t <- decompress <$> BS.readFile "tdata/t10k-labels-idx1-ubyte.gz"
    return (runGet decodeImages s, runGet decodeLabels t)

max_stars_repo_name: pierric/neural-network

max_issues_repo_name: pierric/neural-network

max_forks_repo_name: pierric/neural-network


text:
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE FlexibleContexts #-}
module Data.Array.Accelerate.Utility.Lift.Exp (
   Unlift,
   Unlifted,
   Tuple,
   unlift,
   modify,
   modify2,
   modify3,
   modify4,
   Exp(Exp), expr, atom,
   unliftPair,
   unliftTriple,
   unliftQuadruple,
   asExp,
   mapFst,
   mapSnd,
   fst3,
   snd3,
   thd3,
   indexCons,
   ) where

import qualified Data.Array.Accelerate.Data.Complex as Complex
import qualified Data.Array.Accelerate as A
import Data.Complex (Complex((:+)))
import Data.Array.Accelerate ((:.)((:.)))

import qualified Data.Tuple.HT as Tuple
import Data.Tuple.HT (mapTriple)

{- |
This class simplifies untupling of expressions.
If you have a function

> g :: ((Exp a, Exp b), Exp (c,d)) -> (Exp e, Exp f)

you cannot apply it to an array @arr :: Array sh ((a,b),(c,d))@ using 'A.map'.
Here, the 'modify' function helps:

> modify ((expr,expr),expr) g :: Exp ((a,b),(c,d)) -> Exp (e,f)

The 'expr'-pattern tells, how deep the tuple shall be unlifted.
This way you can write:

> A.map
>    (Exp.modify ((expr,expr),expr) $ \((a,b), cd) -> g ((a,b), cd))
>    arr

'modify' is based on 'unlift'.
In contrast to 'A.unlift' it does not only unlift one level of tupels,
but is guided by an 'expr'-pattern.
In the example I have demonstrated,
how the pair @(a,b)@ is unlifted, but the pair @(c,d)@ is not.
For the result tuple, 'modify' simply calls 'A.lift'.
In contrast to 'A.unlift',
'A.lift' lifts over all tupel levels until it obtains a single 'Exp'.
-}
class
   (A.Elt (Tuple pattern), A.Plain (Unlifted pattern) ~ Tuple pattern) =>
      Unlift pattern where
   type Unlifted pattern
   type Tuple pattern
   unlift :: pattern -> A.Exp (Tuple pattern) -> Unlifted pattern

modify ::
   (A.Lift A.Exp a, Unlift pattern) =>
   pattern ->
   (Unlifted pattern -> a) ->
   A.Exp (Tuple pattern) -> A.Exp (A.Plain a)
modify p f = A.lift . f . unlift p

modify2 ::
   (A.Lift A.Exp a, Unlift patternA, Unlift patternB) =>
   patternA ->
   patternB ->
   (Unlifted patternA -> Unlifted patternB -> a) ->
   A.Exp (Tuple patternA) -> A.Exp (Tuple patternB) -> A.Exp (A.Plain a)
modify2 pa pb f a b = A.lift $ f (unlift pa a) (unlift pb b)

modify3 ::
   (A.Lift A.Exp a, Unlift patternA, Unlift patternB, Unlift patternC) =>
   patternA ->
   patternB ->
   patternC ->
   (Unlifted patternA -> Unlifted patternB -> Unlifted patternC -> a) ->
   A.Exp (Tuple patternA) -> A.Exp (Tuple patternB) ->
   A.Exp (Tuple patternC) -> A.Exp (A.Plain a)
modify3 pa pb pc f a b c =
   A.lift $ f (unlift pa a) (unlift pb b) (unlift pc c)

modify4 ::
   (A.Lift A.Exp a,
    Unlift patternA, Unlift patternB, Unlift patternC, Unlift patternD) =>
   patternA ->
   patternB ->
   patternC ->
   patternD ->
   (Unlifted patternA -> Unlifted patternB ->
    Unlifted patternC -> Unlifted patternD -> a) ->
   A.Exp (Tuple patternA) -> A.Exp (Tuple patternB) ->
   A.Exp (Tuple patternC) -> A.Exp (Tuple patternD) -> A.Exp (A.Plain a)
modify4 pa pb pc pd f a b c d =
   A.lift $ f (unlift pa a) (unlift pb b) (unlift pc c) (unlift pd d)

instance (A.Elt a) => Unlift (Exp a) where
   type Unlifted (Exp a) = A.Exp a
   type Tuple (Exp a) = a
   unlift _ = id

data Exp e = Exp

expr :: Exp e
expr = Exp

{-# DEPRECATED atom "use expr instead" #-}
-- | for compatibility with accelerate-utility-0.0
atom :: Exp e
atom = expr

instance (Unlift pa, Unlift pb) => Unlift (pa,pb) where
   type Unlifted (pa,pb) = (Unlifted pa, Unlifted pb)
   type Tuple (pa,pb) = (Tuple pa, Tuple pb)
   unlift (pa,pb) ab =
      (unlift pa $ A.fst ab, unlift pb $ A.snd ab)

instance
   (Unlift pa, Unlift pb, Unlift pc) =>
      Unlift (pa,pb,pc) where
   type Unlifted (pa,pb,pc) = (Unlifted pa, Unlifted pb, Unlifted pc)
   type Tuple (pa,pb,pc) = (Tuple pa, Tuple pb, Tuple pc)
   unlift (pa,pb,pc) =
      mapTriple (unlift pa, unlift pb, unlift pc) . A.unlift

instance (Unlift pa, A.Slice (Tuple pa), int ~ Exp Int) => Unlift (pa :. int) where
   type Unlifted (pa :. int) = Unlifted pa :. A.Exp Int
   type Tuple (pa :. int) = Tuple pa :. Int
   unlift (pa:.pb) ab =
      (unlift pa $ A.indexTail ab) :. (unlift pb $ A.indexHead ab)

instance (Unlift p) => Unlift (Complex p) where
   type Unlifted (Complex p) = Complex (Unlifted p)
   type Tuple (Complex p) = Complex (Tuple p)
   unlift (preal:+pimag) z =
      unlift preal (Complex.real z)
      :+
      unlift pimag (Complex.imag z)

unliftPair :: (A.Elt a, A.Elt b) => A.Exp (a,b) -> (A.Exp a, A.Exp b)
unliftPair = A.unlift

unliftTriple ::
   (A.Elt a, A.Elt b, A.Elt c) => A.Exp (a,b,c) -> (A.Exp a, A.Exp b, A.Exp c)
unliftTriple = A.unlift

unliftQuadruple ::
   (A.Elt a, A.Elt b, A.Elt c, A.Elt d) =>
   A.Exp (a,b,c,d) -> (A.Exp a, A.Exp b, A.Exp c, A.Exp d)
unliftQuadruple = A.unlift

asExp :: A.Exp a -> A.Exp a
asExp = id

mapFst ::
   (A.Elt a, A.Elt b, A.Elt c) =>
   (A.Exp a -> A.Exp b) -> A.Exp (a,c) -> A.Exp (b,c)
mapFst f = modify (expr,expr) $ \(a,c) -> (f a, c)

mapSnd ::
   (A.Elt a, A.Elt b, A.Elt c) =>
   (A.Exp b -> A.Exp c) -> A.Exp (a,b) -> A.Exp (a,c)
mapSnd f = modify (expr,expr) $ \(a,b) -> (a, f b)

fst3 ::
   (A.Elt a, A.Elt b, A.Elt c) =>
   A.Exp (a,b,c) -> A.Exp a
fst3 = modify (expr,expr,expr) Tuple.fst3

snd3 ::
   (A.Elt a, A.Elt b, A.Elt c) =>
   A.Exp (a,b,c) -> A.Exp b
snd3 = modify (expr,expr,expr) Tuple.snd3

thd3 ::
   (A.Elt a, A.Elt b, A.Elt c) =>
   A.Exp (a,b,c) -> A.Exp c
thd3 = modify (expr,expr,expr) Tuple.thd3

indexCons ::
   (A.Slice ix) => A.Exp ix -> A.Exp Int -> A.Exp (ix :. Int)
indexCons ix n = A.lift $ ix:.n