text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE DeriveGeneric       #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE LambdaCase          #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving  #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeOperators       #-}

module Lib where

import           Control.Monad.Random
import           Data.Binary                  as B
import           Data.Singletons
import           Data.Singletons.Prelude
import           Data.Singletons.TypeLits
import           GHC.Generics                 (Generic)
import           GHC.Natural
import           Numeric.LinearAlgebra.Static
import           Prelude

{-# ANN module ("HLint: ignore Eta reduce" :: Prelude.String) #-}
{-# ANN module ("HLint: ignore Redundant lambda" :: Prelude.String) #-}

------------------------------------------------------------------------------
--import           Test.Hspec
------------------------------------------------------------------------------
{-
https://blog.jle.im/entry/practical-dependent-types-in-haskell-2.html

Practical Dependent Types in Haskell 2: Existential Neural Networks and Types at Runtime
by Justin Le ♦ Thursday June 30, 2016

using dependent types to write type-safe neural networks!

Part 1 dealt with types that are fixed at compile-time.

Part 2 types that depend runtime info
-}
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

runLayer
  :: (KnownNat i, KnownNat o)
  => Weights i o
  -> R i
  -> R o
runLayer (W wB wN) v = wB + wN #> v

-- * Network
--

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 :&~

deriving instance (KnownNat i, KnownNat o) => Show (Network i hs o)

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

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

type for neural networks:

:k Network
Network :: Nat -> [Nat] -> Nat -> *
           ^      ^        ^
           |      |        size of the output vector
           |      list of hidden layer sizes
           size of input vector

above requires knowing entire structure of the network at compile-time

------------------------------------------------------------------------------
An Existential Crisis

two main ways to solve this issue : actually just two styles of doing the same thing.

------------------------------------------------------------------------------
Types hiding behind constructors

a user only needs to know input and output vectors, so Neural Network is R i -> R o

downside of having the structure in the type
- can’t store them in the same list or data structure
  Network 10 '[5,3] 1  and Network 10 '[5,2] 1

If written without the internal structure in the type then could:

data OpaqueNet i o

OpaqueNet a useful type for users
- only exposes types relevant to usage/API
- can store in list or MVar : [OpaqueNet 10 3] and MVar (OpaqueNet 10 3)
- serialize/deserialize

OpaqueNet as an "existential" wrapper over Network ('hs' is existential)
-}
data OpaqueNet :: Nat -> Nat -> * where
  ONet :: Network i hs o -> OpaqueNet i o
{-
Given        Network 6 '[10,6,3] 2
can create OpaqueNet 6           2
When using ONet constructor, the structure of the hidden layers disappears from the type.

Use the Network inside by pattern matching on ONet:
-}
runOpaqueNet
  :: (KnownNat i, KnownNat o)
  => OpaqueNet i o
  -> R i
  -> R o
runOpaqueNet (ONet n) x = runNet n x

numHiddens :: OpaqueNet i o -> Int
numHiddens (ONet n) = go n
 where
  go :: Network i hs o -> Int
  go  = \case
    O     _  -> 0
    _ :&~ n' -> 1 + go n'
{-
ScopedTypeVariables
- enables bringing hs back into scope:

case oN of
  ONet (n :: Network i hs o) -> ...

Above pattern is sometimes called dependent pair:
- pattern matching on ONet yields the hidden existentially quantified type : hs
  and also a value whose type is based on it (Network i hs o)
- like hs “paired” with Network i hs o

Pattern match on the results to give both the type 'hs' and the data structure.

(could have implemented it as ONet :: Sing hs -> Network i hs o -> OpaqueNet i o)

key to making this work
- after pattern matching on ONet, must handle 'hs' in a completely polymorphic way
- cannot assume anything about 'hs'

e.g., not ok:

bad :: OpaqueNet i o -> Network i hs o
bad (ONet n) = n

because the type signature means the caller can decide what 'hs' can be.
but the function has a specific 'hs' : the network that ONet hides.
The caller must accommodate whatever is inside ONet.

------------------------------------------------------------------------------
The Universal and the Existential

central to using existential types
- who has the power to decide what the types will be instantiated as

Most polymorphic functions in Haskell are "universally qualified", e.g.,

    map :: (a -> b) -> [a] -> [b]

the caller of map decides what a and b are

map is defined such that it will work for any a and b that the caller chooses

But for

foo :: [Int] -> OpaqueNet i o

caller of 'foo' chooses i and o
'foo' implementation chooses 'hs'
To use what 'foo' returns, must being able to handle anything 'foo'returns.

summary:

universal   : caller chooses
              function’s implementation must accommodate any choice

existential : function’s implementation chooses
              caller must accommodate any choice

earlier showed to Network i hs o inside the OpaqueNet i o
must deal with it in a parametrically polymorphic way (i.e., handle any 'hs')

------------------------------------------------------------------------------
A familiar friend

OpaqueNet i o  is a "dependent pair"
- pairs 'hs' with 'Network i hs o'

another common term : DEPENDENT SUM

"sum types" are Either-like types that can be one thing or another.

Dependent pairs/existential types actually are similar to Either/sum types

working with existential types is no different than working with Either

Given:

foo :: String -> Either Int Bool

then must handle the result for both cases.
- the function gets to pick which

(can view OpaqueNet i o as being an infinite Either over '[], '[1], '[1,2], ...)

handle Either via pattern matching
OpaqueNet i o is the same
- do not know the type of Network i hs o until pattern matching (a "dependent pattern match")

------------------------------------------------------------------------------
REIFICATION

re-write randomNet with explicit singleton input style
-}
randomNet'
  :: forall m i hs o. (MonadRandom m, KnownNat i, KnownNat o)
  => Sing hs -> m (Network i hs o)
randomNet' = \case
  SNil            ->     O <$> randomWeights
  SNat `SCons` ss -> (:&~) <$> randomWeights <*> randomNet' ss

randomNet
  :: forall m i hs o. (MonadRandom m, KnownNat i, SingI hs, KnownNat o)
  => m (Network i hs o)
randomNet = randomNet' sing
{-
use sing :: SingI hs => Sing hs to call the Sing hs ->-style function
from the SingI hs => one

still need to get the list of integers to the type level
- to create a Network i hs o
- to put into ONet

Use 'SomeSing' (a lot like OpaqueNet)
- wrapping the Sing a inside an existential data constructor
- toSing
  - takes term-level value (here [Integer])
  - returns 'SomeSing' wrapping the type-level value (here [Nat])
- pattern match on SomeSing constructor brings 'a' into scope

SomeSing implementation:

data SomeSing :: * -> * where
    SomeSing :: Sing (a :: k) -> SomeSing k

so,

foo :: SomeSing Bool
foo  = SomeSing STrue

bar :: SomeSing Nat
bar  = SomeSing (SNat :: Sing 10)
-}

mSing :: IO ()
mSing  = do
  putStrLn "How many cats do you own?"
  c <- readLn :: IO Natural
  case toSing c of SomeSing sc -> print (fromSing sc + 2)

xIn :: Int -> SomeSing Nat
xIn  = toSing . intToNatural

xOut :: SomeSing Nat -> Int
xOut  = \case SomeSing n -> naturalToInt (fromSing n)

xio :: Int -> Int
xio  = xOut . xIn
{-
inside the case statement branch the type n :: Nat is in scope.
pattern matching on the SNat constructor gives a KnownNat n instance

toSing can convert value x of type a to a singleton representing type x with kind a

Can use this to write randomONet:
-}
randomONet
  :: (MonadRandom m, KnownNat i, KnownNat o)
  => [Natural]
  -> m (OpaqueNet i o)
randomONet hs = case toSing hs of
  SomeSing ss -> ONet <$> randomNet' ss
{-
converting term-level value to type level type known as reification

original goal is now within reach:
-}
main :: IO ()
main  = do
  putStrLn "What hidden layer structure do you want?"
  hs <- readLn
  n  <- randomONet hs
  case n of
    ONet (net :: Network 10 hs 3) ->
      print net
      -- more stuff with dynamically generated net
{-
------------------------------------------------------------------------------
The Boundary

existentially quantified types (e.g., SomeSing),
enables ability to work with types that depend on runtime results

toSing and SomeSing is a BOUNDARY between untyped world and typed world
reification : cleanly separates the two

Reification boundary similar to boundary between pure code and impure (IO, etc.) code.

working with typed/untyped worlds as the same thing
- write as as possible in typed world
- write only what absolutely must be in untyped world
  (e.g., values from runtime environment)

------------------------------------------------------------------------------
Continuation-Based Existentials

another way to work with existentials that can be more natural to use

when pattern matching on existential, must work with ANY values in constructor

oNetToFoo :: OpaqueNet i o -> Foo
oNetToFoo (ONet n) = f n

f take 'Network i hs o' and work for all 'hs'

could skip the constructor and represent an existential type
as something taking the continuation f and giving it what it needs
-}

type OpaqueNet' i o r = (forall hs. Network i hs o -> r) -> r

{-
"Tell me how you would make an r
- given a 'Network i hs o' that works for any 'hs'
- then I will make it for you”

Uses Rank-N types

Using above is similar to using constructor-style:
-}
runOpaqueNet'
  :: (KnownNat i, KnownNat o)
  => OpaqueNet' i o (R o)
  -> R i
  -> R o
runOpaqueNet' oN x = oN (`runNet` x)
--            :: ((forall hs. Network i hs o -> R o) -> R o)
--            -> R i
--            -> R o

numHiddens' :: OpaqueNet' i o Int -> Int
numHiddens' oN = oN go
 where
  go :: Network i hs o -> Int
  go  = \case
    O     _  -> 0
    _ :&~ n' -> 1 + go n'
--          :: ((forall hs. Network i hs o -> Int) -> Int)
--          -> Int
{-
Above "continuation transformation" formally known as skolemization

wrap
-}

oNet' :: Network i hs o -> OpaqueNet' i o r
oNet' n = \f -> f n

-- randomONet that returns a continuation-style existential:

withRandomONet'0
  :: (MonadRandom m, KnownNat i, KnownNat o)
  => [Natural]
  -> (forall hs. Network i hs o -> m r)
  -> m r
-- aka
-- => [Integer]
-- -> OpaqueNet' i o (m r)
withRandomONet'0 hs f = case toSing hs of
  SomeSing ss -> do
    net <- randomNet' ss
    f net
{-
-- version of `toSing` that returns a skolemized `SomeSing`
withSomeSing :: [Integer]
             -> (forall (hs :: [Nat]). Sing hs -> r)
             -> r
-}

withRandomONet'
  :: (MonadRandom m, KnownNat i, KnownNat o)
  => [Natural]
  -> (forall hs. Network i hs o -> m r)
  -> m r
--aka,
-- => [Integer]
-- -> OpaqueNet' i o (m r)
withRandomONet' hs f = withSomeSing hs $ \ss -> do
  net <- randomNet' ss
  f net

main' :: IO ()
main' = do
  putStrLn "What hidden layer structure do you want?"
  hs <- readLn
  withRandomONet' hs $ \(net :: Network 10 hs 3) ->
    print net
    -- blah blah stuff with our dynamically generated net
{-
Just like the case statement pattern match represented the lexical boundary
between the untyped and typed world when using constructor-style existentials,
the ... $ \net -> ... is also a BOUNDARY for continuation-style existentials.

------------------------------------------------------------------------------
A Tale of Two Styles

two equivalent ways representing/working with existential types
can always convert between them

which one to use or offer can make a difference in code clarity

PROS/CONS

continuation-style
- does not require a data type/constructor
  - can use newtypes, but they force users to learn those types/constructors
    for each existential type returned
  - if there are many existentials, then that means a lot of un/wrapping
- easier to use when functions return existentials
  - especially if intended to immediately use them
    - saves an extraneous pattern match
- when using several existentials at once, continuation-style is better
  because each nested existential doesn’t force another level of indentation:

foo = withSomeSing x $ \sx ->
      withSomeSing y $ \sy ->
      withSomeSing z $ \sz ->
        -- ...
vs.

foo = case toSing x of
        SomeSing sx ->
          case toSing y of
            SomeSing sy ->
              case toSing z of
                SomeSing sz ->
                  -- ...
In monadic context, use do notation and ScopedTypeVariables for a non-nested style

main = do
  ONet n1 <- randomONet [7,5,3] :: IO (OpaqueNet 10 1)
  ONet n2 <- randomONet [5,5,5] :: IO (OpaqueNet 10 1)
  ONet n3 <- randomONet [5,4,3] :: IO (OpaqueNet 10 1)
  hs <- readLn
  ONet (n4 :: Network 10 hs 1) <- randomONet hs
    -- ...

nicer than

main = withRandomONet' [7,5,3] $ \n1 ->
       withRandomONet' [5,5,5] $ \n2 ->
       withRandomONet' [5,4,3] $ \n3 -> do
         hs <- readLn
         withRandomONet' hs $ \(n4 :: Network 10 hs 1) -> do
           -- ...
trick less useful for functions like toSing where things are not returned in a monad
You could wrap it in Identity, but that’s kind of silly:

foo = runIdentity $ do
        SomeSing sx <- Identity $ toSing x
        SomeSing sy <- Identity $ toSing y
        SomeSing sz <- Identity $ toSing z
        return $ -- ...

Constructor-style
- necessary for writing typeclass instances
  - i.e., cannowrite a Show instance for (forall hs. Network i hs o -> r) -> r
    but can write one for OpaqueNet i o
  - i.e., Binary de/serialization instances

Haskell does not support using Rank-N types as arguments to type constructors,
so [OpaqueNet i o] is possible,
but not [OpaqueNet' i o r] or [(forall hs. Network i hs o -> r) -> r]

MVar (OpaqueNet i o) is OK
but not MVar ((forall hs. Network i hs o -> r) -> r)

(latter are known as impredicative types)

If the type constructor is a Monad, you can use ContT-style skolemization,
like (forall hs. Network i hs o -> [r]) -> [r] and (forall hs. Network i hs o -> IO r) -> IO r.
But doesn’t work for MVar and other useful type constructors.

When writing functions that take existentials as inputs,
the constructor-style is arguably more natural. But barely.

For example, we wrote a function to find the number of hidden layers in a network earlier:

numHiddens :: OpaqueNet i o -> Int

But the continuation-style version has a slightly messier type:

numHiddens' :: ((forall hs. Network i hs o -> Int) -> Int)
            -> Int

Even with with the type synonym, it’s still a little awkward:

numHiddens' :: OpaqueNet' i o Int -> Int

This is why you’ll encounter many more functions returning continuation-style existentials
in libraries than taking them, for the most part.

------------------------------------------------------------------------------
Serializing Networks

apply existential types and reification to serialization

------------------------------------------------------------------------------
Recap on the Binary Library

Serializing networks whose sizes are statically in
their types — is straightforward --- an advantages of having sizes in types

binary library : typeclass-based de/serializing

monads (Get, Put) for describing serialization schemes and also a
typeclass used to provide serialization instructions for different
types.

usually do not write instances
use GHC generics:
-}
data Weights i o = W
  { wBiases :: !(R o)
  , wNodes  :: !(L o i)
  } deriving (Show, Generic)

instance (KnownNat i, KnownNat o) => Binary (Weights i o)
{-
gives get and put:

get :: (KnownNat i, KnownNat o) => Get (Weights i o)
put :: (KnownNat i, KnownNat o) => Weights i o -> Put
However, for GADTs like Network, we have to things manually.

------------------------------------------------------------------------------
Serializing Network
-}
putNet :: (KnownNat i, KnownNat o) => Network i hs o -> Put
putNet  = \case
  O w     -> put w
  w :&~ n -> put w *> putNet n
{-
Usually a flag is needed  to tell what constructor is being serialized,
so the deserializer can know what constructor to deserialize.
But need needed for Network because it is known know how many :&~ layers to expect
via the type: deserializing Network 5 '[10,6,3] 2 means three (:&~)’s and one O.
-}
getNet :: forall i hs o. (KnownNat i, KnownNat o) => Sing hs -> Get (Network i hs o)
getNet  = \case
  SNil            ->     O <$> get
  SNat `SCons` ss -> (:&~) <$> get <*> getNet ss
{-
getNet similar to randomNet': pattern match on 'hs' to get constructors
and ollow what the singleton’s structure says.

a Network Binary instance : cannot have get take a Sing hs input
- change the arity/type of the function
- switch to SingI-style and have their Binary instances require a SingI hs constraint.
-}
instance (KnownNat i, SingI hs, KnownNat o) => Binary (Network i hs o) where
  put = putNet
  get = getNet sing
{-
------------------------------------------------------------------------------
Serializating OpaqueNet

to de/serialize, must know complete structure

the complete structure is not in the type
- need to encode it as a flag so deserializer will which constructors
-}
hiddenStruct :: Network i hs o -> [Natural]
hiddenStruct = \case
    O _                           -> []
    _ :&~ (n' :: Network h hs' o) -> natVal (Proxy @h) : hiddenStruct n'

{-
natVal :: KnownNat n => Proxy n -> Natural
returns the value-level Natural corresponding to type-level n :: Nat

natVal and hiddenStruct : turn type-level infor ('n', 'hs') into term-level values

inverse of reification functions (like toSing)

REFECTION: type level to the value level
-}
putONet :: (KnownNat i, KnownNat o) => OpaqueNet i o -> Put
putONet (ONet net) = do
  put (hiddenStruct net) -- put the structure (as a binary flag)
  putNet net             -- then put the network itself

-- deserialize: load the list of Integers and reify it
getONet :: (KnownNat i, KnownNat o) => Get (OpaqueNet i o)
getONet = do
  hs <- get              -- load flag
  withSomeSing hs        -- reify it
    (fmap ONet . getNet) -- create it

instance (KnownNat i, KnownNat o) => Binary (OpaqueNet i o) where
  put = putONet
  get = getONet
{-
used constructor-style existential (instead of continuation-style)
because typeclass instances for the latter are not possible

------------------------------------------------------------------------------
An Existence For All

learned
- reification : untyped world to typed world to create tyypes depend on runtime values
- reflection : de/serialization shows type safety in dynamic context

------------------------------------------------------------------------------
Exercises

1. Implement putONet' and getONet' using continuation-style existentials
https://github.com/mstksg/inCode/tree/master/code-samples/dependent-haskell/NetworkTyped2.hs#L188-L193
https://github.com/mstksg/inCode/tree/master/code-samples/dependent-haskell/NetworkTyped2.hs#L195-L203

2. Work with an existential wrapper over the entire network structure (inputs and outputs, too):

data SomeNet where
  SNet :: (KnownNat i, KnownNat o)
       => Network i hs o
       -> SomeNet

(Need KnownNat because of type erasure, to recover original input/output dimensions
back on pattern matching)

And write:

- A function to convert SomeNets to OpaqueNets.
  Return OpaqueNet with existentially quantified i and o in continuation-style.
  (Also try a data type to return it in constructor-style)
https://github.com/mstksg/inCode/tree/master/code-samples/dependent-haskell/NetworkTyped2.hs#L231-L234

- randomSNet, returning m SomeNet.
https://github.com/mstksg/inCode/tree/master/code-samples/dependent-haskell/NetworkTyped2.hs#L236-L245

While you’re at it, write it to return a random continuation-style SomeNet, too! (See the type of withRandomONet' for reference on how to write the type)

The binary instance for SomeNet.

Hint: Remember natVal :: KnownNat n => Proxy n -> Integer!

Hint: Remember that toSomeSing also works for Integers, to get Sings for Nats, too!

A bit of a stretch, because the set of all [Nat]s is non-enumerable and uncountable, but hopefully you get the picture!↩

Recall that I recommend (personally, and subjectively) a style where your external API functions and typeclass instances are implemented in SingI a => style, and your internal ones in Sing a -> style. This lets all of your internal functions fit together more nicely (Sing a -> style tends to be easier to write in, especially if you stay in it the entire time, because Sings are normal first-class values, unlike those global and magical typeclasses) while at the same time removing the burden of calling with explicit singletons from people using the functionality externally.↩

In older versions of singletons, before GHC 8 and TypeInType, we had to implement it using “kind proxies”. Don’t worry if you’re following along in 7.10; the basic usage of SomeSing is essentially identical either way.↩

Skolemization is probably one of the coolest words you’ll encounter when learning/using Haskell, and sometimes just knowing that you’re “skolemizing” something makes you feel cooler. Thank you Thoralf Skolem. If you ever see a “rigid, skolem” error in GHC, you can thank him for that too! He is also the inspiration behind my decision to name my first-born son Thoralf. (My second son’s name will be Curry)↩

It even lets you write Storable instances!↩
-}

max_stars_repo_name: haroldcarr/learn-haskell-coq-ml-etc

max_issues_repo_name: haroldcarr/learn-haskell-coq-ml-etc

max_forks_repo_name: haroldcarr/learn-haskell-coq-ml-etc


text:
module Mandroots where

import Data.Complex

data Qaudratic a = Qaudratic a a a deriving (Show, Eq)

solveQuadratic :: Qaudratic (Complex Double) -> (Complex Double, Complex Double)
solveQuadratic (Qaudratic a b c) = let d = sqrt (b*b - 4 * a * c) in ((-b + d)/(2 * a), (-b - d)/(2 * a))

evalQuadratic :: Num a => Qaudratic a -> a -> a
evalQuadratic (Qaudratic a b c) x =  x*x * a + x * b + c

rootsIteration :: Qaudratic (Complex Double) -> Qaudratic (Complex Double)
rootsIteration q@(Qaudratic a b c) = let (rt1, rt2) = solveQuadratic q in (Qaudratic rt1 rt2 c)

type QCD = Qaudratic (Complex Double)

timeToEscape :: (QCD -> QCD) -> QCD -> Int -> Int
timeToEscape f q max = count q 0
  where
    count p n
      | n ==  max = max
      | let a = evalQuadratic p (1 :+ 0) in realPart (a * conjugate a) > 25 = n
      | otherwise = count (f q) (n + 1)

ext: hs

max_stars_repo_name: AlwinHughes/Derivative

max_issues_repo_name: AlwinHughes/Derivative

max_forks_repo_name: AlwinHughes/Derivative


text:
{-# LANGUAGE CPP #-}
#if __GLASGOW_HASKELL__ > 707
{-# LANGUAGE AllowAmbiguousTypes #-}
#endif
{-# LANGUAGE OverlappingInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE QuasiQuotes #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE TypeFamilies #-}
module T1 where
import GHC.Exts (Constraint)
import DimMat
import Numeric.Units.Dimensional.Prelude
import Numeric.Units.Dimensional

import Text.PrettyPrint.ANSI.Leijen hiding (dot)

import GHC.TypeLits

import qualified Numeric.LinearAlgebra as H
import qualified Prelude as P

{- Example from http://ctms.engin.umich.edu/CTMS/index.php?example=InvertedPendulum&section=ControlStateSpace -}

-- not really stated on that page, but if you go back a couple of pages in their derivation
-- you can see that the type of u is a 1x1 matrix whose sole element is a force

massOfCart = (0.5 :: Double) *~ kilo gram
massOfPendulum = 0.2 *~ kilo gram
coefficientOfFrictionForCart = 0.1 *~ (newton / meter / second)
lengthToPendulumCenterOfMass = 0.3 *~ meter
massMomentOfInertiaOfPendulum = 0.006 *~ (kilo gram * meter * meter)
g = (9.8::Double) *~ (meter / second ^ pos2)

p = massMomentOfInertiaOfPendulum * (massOfCart + massOfPendulum)
  + (massOfCart*massOfPendulum*lengthToPendulumCenterOfMass*lengthToPendulumCenterOfMass)

a22 = negate (massMomentOfInertiaOfPendulum + massOfPendulum * lengthToPendulumCenterOfMass * lengthToPendulumCenterOfMass)
                  * coefficientOfFrictionForCart / p
a23 = (massOfPendulum * massOfPendulum * g * lengthToPendulumCenterOfMass * lengthToPendulumCenterOfMass) / p
a42 = negate (massOfPendulum * lengthToPendulumCenterOfMass * coefficientOfFrictionForCart) / p
a43 = massOfPendulum * g * lengthToPendulumCenterOfMass*(massOfCart + massOfPendulum)/p

{-

>>> inv aSmall `sub` aInvSmall
2><2     1                      m
s        -8.881784197001252e-16 -2.220446049250313e-16
m^-1 s^2 0.0                    0.0

>>> inv aSmall `multiply` aSmall
2><2   1                     m s^-1
1      1.0                   -6.732461810265988e-15
m^-1 s 5.014435047745458e-19 0.9999999999999999

-}
aSmall = [matD| a22, a23; a42, a43 |]

aInvSmall = scale (_1 / det aSmall)
    [matD| a43, negate a23; negate a42, a22 |]

-- id2 = Proxy :: Proxy (0 * a)

colMat :: D '(a, '[ra, '[]]) e -> D '(a, '[ra, '[]]) e
colMat x = x

{-
charEq lam1 lam2 x1 x2 =
        (det ( aSmall `sub` [matD| lam1, _0; _0, lam2 |]),
         colMat $ scale lam1 x1 `sub` dot aSmall x1,
         colMat $ scale lam2 x2 `sub` dot aSmall x2,
         hconcat x1 x2)
         --
         -- scale lam1 x1 `sub` dot aSmall x1, -- `asTypeOf` (undefined :: D '(s, ['[t],'[ ]]) Double),
         -- scale lam2 x2 `sub` dot aSmall x2) --  `asTypeOf` (undefined :: D '(s, ['[t],'[ ]]) Double),
         -- hconcat x1 x2 `asTypeOf` (undefined :: D '(t, ['[v],'[u]]) Double))
-- still has ambiguity. Do the SVD instead?
-}

{-
>>> cmap (Scale' ((2::Double) *~ second)) aSmall
2><2 1                    m s^-1
1    -0.36363636363636365 5.3454545454545475
m^-1 -0.9090909090909091  62.36363636363637

alternative to scale. Also more fancy examples (say swap metres and seconds
powers)...

-}
newtype Scale' f a = Scale' (Quantity f a)
instance (Num a, a ~ Double,
          x ~ Quantity d a,
          y ~ Quantity d' a,
          Mul f d d')
        => ApplyAB (Scale' f a) x y where
  applyAB (Scale' n) x = n * x

b21 = (massMomentOfInertiaOfPendulum + (massOfPendulum * lengthToPendulumCenterOfMass * lengthToPendulumCenterOfMass)) / p
b41 = massOfPendulum * lengthToPendulumCenterOfMass / p

-- example state value
x = [matD| 1.0 *~ meter;
           0.2 *~ (meter / second);
           _0 :: Dimensionless Double;
           0.1 *~ (second^neg1) |]

dx = scale (_1 / (1 *~ second)) x

-- example control input
u = [matD| (0 :: Double) *~ newton |]

--xDot = (a `multiply` x) `add` (b `multiply` u)

y = [matD| 1 *~ meter; _0 :: Dimensionless Double  |]

-- :t isLTI (1 *~ second) x u y, given that x u and y have monomorphic
-- types properly infers constraints on the a,b,c,d arguments!
isLTI time x u y a b c d =
    (scale (_1 /time) (colMat x) `add` dot a x `add` dot b (colMat u),
     colMat $ y `add` dot c x `add` dot d u)

testIsLTI =
  (\ a b c d -> case isLTI (1 *~ second) x u y a (colMat b) c (colMat d) of
   _ -> do
    print $ vsep
        [text "A = " </> indent 0 (pretty (zeroes `asTypeOf` a)),
         text "B = " </> indent 0 (pretty (zeroes `asTypeOf` b)),
         text "C = " </> indent 0 (pretty (zeroes `asTypeOf` c)),
         text "D = " </> indent 0 (pretty (zeroes `asTypeOf` d))]
    ) undefined undefined undefined undefined

{- | data type encoding units required by
http://en.wikibooks.org/wiki/Control_Systems/State-Space_Equations#State-Space_Equations

To refresh:

> dxs/dt = A xs + B us
> ys = C xs + D us

the units of d/dt are 1/iv
-}
class LiSystem (iv :: *) (xs :: (*, [[*]])) (ys :: (*, [[*]])) (us :: (*, [[*]]))
            (a :: (*, [[*]])) (b :: (*, [[*]])) (c :: (*, [[*]])) (d :: (*, [[*]]))
instance LiSystemCxt dxs iv xs ys us a b c d
    => LiSystem iv xs ys us a b c d

type LiSystemCxt dxs iv xs ys us a b c d =
   (Dot a xs dxs,
    Dot b us dxs,
    Dot c xs ys,
    Dot d us ys,
    DotS iv dxs xs
    )

-- | identity function but constrains types
isLiTuple ::
    (LiSystem iv xs ys us a b c d,
    t ~ (D a e, D b e, D c e, D d e)) =>
    t -> t
isLiTuple x = x

isExample ::
     (LiSystem DTime
        '(DLength , '[[DVelocity, DPlaneAngle, DAngularVelocity]])
        '(DLength, '[ '[DPlaneAngle]])
        '(DForce, '[ '[] ] )
        a b c d, t ~ (D a e, D b e, D c e, D d e)) =>
    t -> t
isExample x = x

pendulum = case (a',b',c',d') of
             sys -> isExample sys -- `const` evaluate sys
         where
           a' = [matD| _0, _1, _0, _0;
                       _0, a22, a23, _0;
                       _0, _0, _0, _1;
                       _0, a42, a43, _0 |]
           b' = [matD| _0;
                       b21;
                       _0;
                       b41 |]
           c' = [matD| _1, _0, _0, _0;
                       _0, _0, _1, _0 |]
           d' = zeroes

{-
-- poles (a,_,_,_) = eigenvalues a

-}
-- causes no problems for AV
-- evaluatePendulum = evaluate pendulum

{-
evaluate ::
    ( -- require all matrices to be at least 1x1
      -- really ought to be part of the MultiplyCxt
     a ~ [_1 ': _2, '[] ': _3],
     b ~ [_4 ': _5, '[] ': _6],
     c ~ [_7 ': _8, '[] ': _9],
     d ~ [_a ': _b, '[] ': _c],
     H.Field e,
     LiSystemCxt dxs iv xs ys us a b c d) =>
    (DimMat a e, DimMat b e, DimMat c e, DimMat d e)
    -> DimMat [xs, '[ '[] ]] e
    -> DimMat [us, '[ '[] ]] e
    -> (DimMat [dxs, '[ '[] ]] e, DimMat [ys, '[ '[] ]] e)
    -}
evaluate (a,b,c,d) x u = case (a `dot` x) `add` (b `dot` u) of
    xDot -> case (c `dot` x) `add` (d `dot` u) of
     y -> (xDot, y)

{-

ctrb = let (a,b,_,_) = pendulum
        in [blockD| b, a `multiply` b, a `multiply` a `multiply` b, a `multiply` a `multiply` a `multiply` b |]
        -}

max_stars_repo_name: aavogt/DimMat

max_issues_repo_name: aavogt/DimMat

max_forks_repo_name: aavogt/DimMat


text:
module Main where
import Layer
import Models
import Numeric.LinearAlgebra
import ImageProcess
import Optimization
import Trainer
import Data.Vector.Storable (Vector, length, toList)

main :: IO()
main =
  do
    let inputLayer = linearLayerInit 2 2
    let sigmoidLayer = sigmoidLayerInit 2 2
    let outputLayer = linearLayerInit 2 1 in
      let model = Model [inputLayer, sigmoidLayer, outputLayer] in
      let input = (1><2)[1.0, 2.0]::Matrix R in
      let output = forward input model in
      print output

testImageRegression :: FilePath -> IO()
testImageRegression filepath = do
  imageM <- getImageInfo filepath
  case imageM of
    Nothing -> putStrLn "Nohthing ! Guess what happened ? :) " >> return ()
    Just imagePixel8 ->
      let rgbs = getImageRGBs imagePixel8
          w = getImageWidth imagePixel8
          h =  getImageHeight imagePixel8
          rgbsList = Data.Vector.Storable.toList rgbs
          (rlist, glist, blist) = genRGBList rgbsList
          Just rMat = genMatrixFromList rlist w h
          Just gMat = genMatrixFromList glist w h
          Just bMat = genMatrixFromList blist w h in
          --- test
          print "prepare the model!" >>
          let rinputLayer = linearLayerInit w h
              rreluLayer1 = reluLayerInit w h
              rreluLayer2 = reluLayerInit w h
              rreluLayer3 = reluLayerInit w h
              routputLayer=  linearLayerInit w h
              rmodel = Model [rinputLayer, rreluLayer1, rreluLayer2, rreluLayer3, routputLayer]
              ginputLayer = linearLayerInit w h
              greluLayer1 = reluLayerInit w h
              greluLayer2 = reluLayerInit w h
              greluLayer3 = reluLayerInit w h
              goutputLayer=  linearLayerInit w h
              gmodel = Model [ginputLayer, greluLayer1, greluLayer2, greluLayer3, goutputLayer]
              binputLayer = linearLayerInit w h
              breluLayer1 = reluLayerInit w h
              breluLayer2 = reluLayerInit w h
              breluLayer3 = reluLayerInit w h
              boutputLayer=  linearLayerInit w h
              bmodel = Model [binputLayer, breluLayer1, breluLayer2, breluLayer3, boutputLayer] in
              print "begin to train the model" >>
              let trainConfig = TrainConfig {
                    trainConfigRowNum=w,
                    trainConfigColNum=h,
                    trainConfigLearningRate=0.95,
                    trainConfigMomentum=1.0,
                    trainConfigDecay=0.95,
                    trainConfigEpsilon=0.00001,
                    trainConfigBatchSize=100,
                    trainConfigOptimization="adadelta"
                    } in
              trainSingleData rmodel trainConfig rMat

max_stars_repo_name: neutronest/vortex

max_issues_repo_name: neutronest/vortex

max_forks_repo_name: neutronest/vortex


text:
{-# LANGUAGE DeriveDataTypeable #-}

module Statistics.Distribution.Combo
       (ComboDistribution(..), combo, combo', rs, create, comboStandard,
        getComboRandoms)
       where

import Control.Monad (replicateM)
import Control.Monad.Primitive
import Data.Typeable (Typeable)
import qualified Statistics.Distribution as D
import qualified Statistics.Distribution.Laplace as L
import qualified Statistics.Distribution.Normal as N
import System.Random.MWC
import qualified System.Random.MWC.Distributions as MWC

data ComboDistribution =
  CD {cdSplit :: Double -- probability of normal distribution
     ,cdMean :: Double -- (both normal mean and glue point)
     ,cdStd :: Double -- normal std
     ,cdLaplaceR :: Double -- exponetial deviation right-side
     ,cdLaplaceL :: Double -- exponential deviation left-side
     ,cdSplitLaplace :: Double   -- probability of right-side
     }
  deriving (Eq,Read,Show,Typeable)

instance D.Distribution ComboDistribution where
  cumulative = cumulative
  complCumulative = complCumulative

instance D.ContDistr ComboDistribution where
  density = density
  quantile = quantile

instance D.Mean ComboDistribution where
  mean (CD split m _ r l sl) =
    (split * m) +
    (1 - split) *
    (sl * r +
     (1 - sl) *
     l +
     m)
  {-# INLINE mean #-}

instance D.Variance ComboDistribution where
  variance (CD split _ s r l sl) =
    split * s +
    (1 - split) *
    (sl /
     ((1 / l) *
      (1 / l)) +
     (1 - sl) /
     ((1 / r) *
      (1 / r)))
  {-# INLINE variance #-}

instance D.MaybeMean ComboDistribution where
  maybeMean = Just . D.mean

instance D.MaybeVariance ComboDistribution where
  maybeStdDev = Just . D.stdDev
  maybeVariance = Just . D.variance

instance D.ContGen ComboDistribution where
  -- genContVar :: (PrimMonad m) => ComboDistribution -> Gen (PrimState m) -> m Double
  genContVar (CD split m s r l sl) g =
    do split' <- uniform g
       sl' <- uniform g
       n <- MWC.normal m s g
       e <- MWC.exponential 1 g
       return $
         if split' < split
            then n
            else if sl' < sl
                    then m + e / r
                    else m - e / l

density :: ComboDistribution -> Double -> Double
density (CD split m s r l sl) x =
  split *
  D.density (N.normalDistr m s) x +
  (1 - split) *
  D.density (L.laplace (1 / r)
                       (1 / l)
                       sl
                       m)
            x

cumulative :: ComboDistribution -> Double -> Double
cumulative (CD split m s r l sl) x =
  split *
  D.cumulative (N.normalDistr m s)
               x +
  (1 - split) *
  D.cumulative
    (L.laplace (1 / r)
               (1 / l)
               sl
               m)
    x

complCumulative :: ComboDistribution -> Double -> Double
complCumulative (CD split m s r l sl) x =
  split *
  D.complCumulative (N.normalDistr m s)
                    x +
  (1 - split) *
  D.complCumulative
    (L.laplace (1 / r)
               (1 / l)
               sl
               m)
    x

quantile :: ComboDistribution -> Double -> Double
quantile d p
  | p == 0 = -inf
  | p == 1 = inf
  | p == 0.5 = cdMean d
  | p > 0 && p < 1 = x
  | otherwise = error $
                 "Statistics.Distribution.Combo.quantile: p must be in [0,1] range. Got: " ++
                 show p
  where x =
          D.findRoot d
                     p
                     0
                     (-inf)
                     inf
        inf = 1 / 0

comboStandard :: ComboDistribution
comboStandard = CD 0.5 0.0 1.0 1.0 1.0 0.5

-- | Create the combo distribution.
combo :: Double             -- probability of normal
      -> Double             -- ^ mean
      -> Double             -- ^ stdev
      -> Double             -- ^ right-side lambda
      -> Double            -- ^ left-side lambda
      -> Double            -- ^ probability of right-side
      -> ComboDistribution
combo = CD

{-# INLINE combo #-}

-- | Create the combo distribution from a paramter list.
combo' :: [Double] -> ComboDistribution
combo' (split:m:s:r:l:sl:_) = CD split m s r l sl
combo' _ = comboStandard

{-# INLINE combo' #-}

rs :: PrimMonad m
   => Gen (PrimState m) -> [Double] -> Int -> m [Double]
rs g d n =
  replicateM n $
  D.genContVar (combo' d)
               g

getComboRandoms :: Int -> IO [Double]
getComboRandoms n =
  do g <- create
     let d = combo 0.5 0.0 1.2 1.2 1.2 0.5
     replicateM n $
       D.genContVar d g


text:
module Scheme.Data
  ( LispVal(..)
  , LispError(..)
  , showVal
  , ThrowsError
  , extractValue
  , trapError
  ) where

import Control.Monad.Error
import Data.Array
import Data.Complex (Complex(..))
import Text.ParserCombinators.Parsec hiding (spaces)

data LispVal
  = Atom String
  | List [LispVal]
  | DottedList [LispVal]
               LispVal
  | Number Integer
  | Ratio Rational
  | Float Double
  | Complex (Complex Double)
  | String String
  | Char Char
  | Bool Bool
  | Vector (Array Int LispVal)

instance Show LispVal where
  show = showVal

showVal :: LispVal -> String
showVal (String contents) = "\"" ++ contents ++ "\""
showVal (Atom name) = name
showVal (Number contents) = show contents
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (List contents) = "(" ++ unwordsList contents ++ ")"
showVal (DottedList head tail) =
  "(" ++ unwordsList head ++ " . " ++ showVal tail ++ ")"

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

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

instance Show LispError where
  show = showError

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

type ThrowsError = Either LispError

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

trapError action = catchError action (return . show)

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

max_stars_repo_name: underscorenico/Scheme

max_issues_repo_name: underscorenico/Scheme

max_forks_repo_name: underscorenico/Scheme


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeFamilyDependencies #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE GADTs #-}
{-# OPTIONS_GHC -Wwarn #-}
module Language.Test where

import Control.Arrow
import Data.Complex
import Data.Proxy
import Type.Reflection

import GHC.TypeLits

data INat = Z | S INat

data family Sing :: k -> *

data instance Sing (n :: INat) where
  SZ :: Sing 'Z
  SS :: Sing n -> Sing ('S n)

class IsSing a where sing :: Sing a
instance IsSing 'Z where sing = SZ
instance IsSing n => IsSing ('S n) where sing = SS sing
type SINat (n :: INat) = Sing n

newtype Unit a = MkUnit () deriving Show

unit :: Unit a
unit = MkUnit ()

type family Prod (n :: INat) (a :: *) = r | r -> n a where
  Prod 'Z a = Unit a
  Prod ('S n) a = (a, Prod n a)

fmapP' :: SINat n -> (a -> b) -> Prod n a -> Prod n b
fmapP' SZ     _ = const unit
fmapP' (SS m) f = f{- ALLOCATE HERE -} *** fmapP' m f

fmapPIx :: SINat n -> (Int -> a -> b) -> Int -> Prod n a -> Prod n b
fmapPIx SZ _ _ = const unit
fmapPIx (SS m) f k = f k *** fmapPIx m f (k+1)

fmapP :: IsSing n => (a -> b) -> Prod n a -> Prod n b
fmapP = fmapP' sing

type family Div2 (n :: INat) where
  Div2 'Z = 'Z
  Div2 ('S 'Z) = 'Z
  Div2 ('S ('S n)) = 'S (Div2 n)

div2 :: SINat n -> SINat (Div2 n)
div2 SZ = SZ
div2 (SS SZ) = SZ
div2 (SS (SS n)) = SS (div2 n)

assocL :: (a, (b, c)) -> ((a, b), c)
assocL (x, (y, z)) = ((x, y), z)

red :: SINat n -> ((a, a) -> a) -> Prod n a -> Prod (Div2 n) a
red SZ _ = id
red (SS SZ) _ = snd
red (SS (SS n)) f = assocL >>> f{- ALLOCATE HERE -} *** red n f

fold :: SINat n -> ((a, a) -> a) -> a -> Prod n a -> a
fold SZ _ z = const z
fold (SS SZ) _ _ = fst
fold n f z = red n f >>> fold (div2 n) f z

reduce :: (IsSing n, Monoid a) => Prod n a -> a
reduce = fold sing (uncurry mappend) mempty

data Fin (m :: INat) where
  FZ :: Fin ('S m)
  FS :: Fin m -> Fin ('S m)

data Le (n :: INat) (m :: INat) where
  LZ :: Le 'Z m
  LS :: Le n m -> Le ('S n) ('S m)

class (IsSing n, IsSing m) => IsLe (n :: INat) (m :: INat) where
  le :: Sing n -> Le n m
instance IsSing n => IsLe 'Z n where
  le SZ = LZ
instance IsLe n m => IsLe ('S n) ('S m) where
  le (SS n) = LS (le n)

get' :: Le n m -> Prod ('S m) a -> a
get' LZ     = fst
get' (LS n) = snd >>> get' n

get :: IsLe n m => SINat n -> Prod ('S m) a -> a
get n = get' (le n)

zw' :: forall a b c (n :: INat). SINat n -> ((a, b) -> c) -> (Prod n a, Prod n b) -> Prod n c
zw' SZ     _ = const unit
zw' (SS m) f = swap >>> f *** zw' m f

toNum :: forall (k :: INat) a. Num a => SINat k -> a
toNum SZ = 0
toNum (SS sm) = 1 + toNum sm

expn :: Float -> Float -> Complex Float -> Complex Float
expn n k c = cis (-2 * pi * k /  n) * c {- multiply exponential -}

add :: (Complex Float, Complex Float) -> Complex Float
add = uncurry (+)

sub :: (Complex Float, Complex Float) -> Complex Float
sub = uncurry (-)

type family Add (n :: INat) (m :: INat) where
  Add 'Z n = n
  Add ('S n) m = 'S (Add n m)

addn :: SINat n -> SINat m -> SINat (Add n m)
addn SZ n = n
addn (SS n) m = SS (addn n m)

type family Mul (n :: INat) (m :: INat) where
  Mul 'Z n = 'Z
  Mul ('S n) m = Add m (Mul n m)

mul :: SINat n -> SINat m -> SINat (Mul n m)
mul SZ _ = SZ
mul (SS n) m = addn m (mul n m)

type family Pow2 (n :: INat) where
  Pow2 'Z = 'S 'Z
  Pow2 ('S n) = Add (Pow2 n) (Pow2 n)

pow2 :: SINat n -> SINat (Pow2 n)
pow2 SZ = SS SZ
pow2 (SS n) = addn (pow2 n) (pow2 n)

-- Hack: proper way would be to traverse 'n' and 'm', but it is quite inefficient

szr :: SINat n -> Add n 'Z :~: n
szr SZ = Refl
szr (SS n) = case szr n of
               Refl -> Refl

ssr :: SINat n -> SINat m -> Add n ('S m) :~: 'S (Add n m)
ssr SZ _ = Refl
ssr (SS n) m = case ssr n m of
                 Refl -> Refl

addComm :: SINat n -> SINat m -> Add n m :~: Add m n
addComm SZ m = case szr m of Refl -> Refl -- unsafeCoerce Refl
addComm (SS n) m = case (addComm n m, ssr m n) of
                     (Refl, Refl) -> Refl

--deinterleave
split :: forall a (n :: INat). SINat n -> Prod (Add n n) a -> (Prod n a, Prod n a)
split SZ = id &&& id
split (SS n) =
  case addComm n (SS n) of
    Refl -> assocL >>> id *** split n >>> swap

assocR :: ((a, b), c) -> (a, (b, c))
assocR ((a, b), c) = (a, (b, c))

cat :: SINat n -> (Prod n a, Prod m a) -> Prod (Add n m) a
cat SZ = snd
cat (SS x) = assocR >>> id *** cat x

fromINat :: forall a n. Num a => SINat n -> a
fromINat SZ = 0
fromINat (SS x) = 1 + fromINat x

fromTo :: forall a n. Num a => SINat n -> Prod n a
fromTo = go 0

go :: forall a n. Num a => a -> SINat n -> Prod n a
go _ SZ = unit
go n (SS m) = (n, go (1 + n) m)

dft :: SINat n -> Prod (Pow2 n) (Complex Float) -> Prod (Pow2 n) (Complex Float)
dft SZ = id
dft (SS x)
   = split p2x >>>
     dft x {-EVENS-} *** dft x {-ODDS-} >>>
     id *** fmapPIx p2x (expn p2sx . fromIntegral) 0 {- Multiply by exponential -} >>>
     zw' p2x add {- Left side -} &&& zw' p2x sub {- right side -} >>>
     cat p2x
  where
    p2x = pow2 x
    p2sx :: Float
    p2sx = fromIntegral $ ((2 ^ (fromINat (SS x) :: Integer)) :: Integer)

--------------------------------------------------------------------------------
-- FFT BELOW -------------------------------------------------------------------

-- PRIM FUNCTIONS --------------------------------------------------------------

-- baseFFT is one of your PRIM funcs: you can copy from some online FFT implementation: e.g. Rosetta Code

baseFFT :: RealFloat a => [Complex a] -> [Complex a]
baseFFT [] = []
baseFFT [x] = [x]
baseFFT xs = zipWith (+) ys ts ++ zipWith (-) ys ts
    where n = length xs
          ys = baseFFT evens
          zs = baseFFT odds
          (evens, odds) = splitList xs
          ts = zipWith (\z k -> (exp' k n) * z) zs [0::Integer ..]

exp' :: (Floating a1, Integral a2, Integral a3) => a2 -> a3 -> Complex a1
exp' k n = cis $ -2 * pi * (fromIntegral k) / (fromIntegral n)

splitList :: [a] -> ([a], [a])
splitList [] = ([], [])
splitList [x] = ([x], [])
splitList (x:y:xs) = (x:xt, y:yt) where (xt, yt) = splitList xs

addPadding :: Num a => SINat n -> [a] -> [a]
addPadding sz l = l ++ replicate (padding $ length l) 0
  where
    padding k = 2 ^ (max ((fromINat sz) :: Integer) (ceiling (logBase 2 $ fromIntegral k :: Double))) - k

concatenate :: ([a], [a]) -> [a]
concatenate = uncurry (++)

{- Below should be your prim functions -}
mulExp :: (RealFloat a1, Integral a3) =>
                a3 -> Int -> [Complex a1] -> [Complex a1]
mulExp p2sx i l = zipWith (\k z -> exp' k (p2sx * fromIntegral len) * z) [i * len ..] l
  where
    len = length l
addc :: ([Complex Float], [Complex Float]) -> [Complex Float]
addc = uncurry $ zipWith (+)
subc :: ([Complex Float], [Complex Float]) -> [Complex Float]
subc = uncurry $ zipWith (-)

-- END PRIM FUNCTIONS ----------------------------------------------------------

-- ARROW DEFINITIONS -----------------------------------------------------------
---  Should not change, except types

type family Tree (n :: INat) (a :: *)  where
  Tree 'Z a = a
  Tree ('S n) a = (Tree n a, Tree n a)

splitL :: SINat m -> [Complex Float] -> Tree m [Complex Float]
splitL SZ = id
splitL (SS n) = splitList >>> splitL n *** splitL n

merge :: forall (n :: INat). SINat n -> Tree n [Complex Float] -> [Complex Float]
merge SZ  = id {- same: just concatenate all using a Prim -}
merge (SS n) = merge n *** merge n >>> concatenate

fmapTIx :: SINat n -> ((Int, [Complex Float]) -> [Complex Float]) -> Int -> Tree n [Complex Float] -> Tree n [Complex Float]
fmapTIx SZ f k = (\v -> f (k, v))
fmapTIx (SS x) f k = fmapTIx x f k *** fmapTIx x f (k + (2 ^ (fromINat x :: Integer)))

swap :: ((a,b), (c, d)) -> ((a, c), (b, d))
swap = ((fst >>> fst) &&& (snd >>> fst)) &&&
       ((fst >>> snd) &&& (snd >>> snd))

zwT :: SINat n -> (([Complex Float],[Complex Float]) -> [Complex Float])
    -> (Tree n [Complex Float], Tree n [Complex Float]) -> Tree n [Complex Float]
zwT SZ f = f
zwT (SS x) f = swap >>> zwT x f *** zwT x f

fft :: SINat n -> Tree n [Complex Float] -> Tree n [Complex Float]
fft SZ = baseFFT
fft (SS x)
  = fft x {-EVENS-} *** fft x {-ODDS-} >>>
    id *** fmapTIx x (uncurry $ mulExp p2sx) 0 {- Multiply by exponential -} >>>
    zwT x addc {- Left side -} &&& zwT x subc {- right side -}
  where
    p2sx :: Integer
    p2sx = 2 ^ (fromINat (SS x) :: Integer)

fastFourierR :: forall (n :: INat). SINat n -> [Complex Float] -> [Complex Float]
fastFourierR cores = addPadding cores >>> splitL cores >>> fft cores >>> merge cores

---------------------------------------------------------------------------------
-- UTILITY FUNCTIONS

fastFourier :: forall n. (KnownNat n, IsSing (FromNat n))
            => [Complex Float] -> [Complex Float]
fastFourier = fastFourierR (sing :: SINat (FromNat n))

fft8core :: [Complex Float] -> [Complex Float]
fft8core = fastFourier @3
fft16core :: [Complex Float] -> [Complex Float]
fft16core = fastFourier @4

type family ToNat (i :: INat) :: Nat where
  ToNat 'Z = 0
  ToNat ('S n) = 1 + ToNat n

type family FromNat (i :: Nat) :: INat where
  FromNat 0 = 'Z
  FromNat n = 'S (FromNat (n-1))

fmapT :: SINat n -> (a -> b) -> Tree n a -> Tree n b
fmapT SZ f = f
fmapT (SS x) f = fmapT x f *** fmapT x f

data SDict a where
  SDict :: Show a => SDict a

getDict :: Show a => SINat n -> Proxy a -> SDict (Tree n a)
getDict SZ _ = SDict
getDict (SS n) p = case getDict n p of
                        SDict -> SDict

test :: SINat n -> Tree n Int -> String
test n = case getDict n (Proxy :: Proxy Int) of
           SDict -> show

test2 :: Show (Tree n Int) => SINat n -> Tree n Int -> Bool
test2 _ _ =  False

test3 :: SINat n -> Tree n Int -> Bool
test3 n t = case getDict n (Proxy :: Proxy Int) of
              SDict -> test2 n t

max_stars_repo_name: session-arr/session-arr

max_issues_repo_name: session-arr/session-arr

max_forks_repo_name: dcastrop/SAlg


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE MonoLocalBinds #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Lib (benchmark) where

import Test.QuickCheck
import Criterion.Main
import Data.Proxy
import GHC.TypeLits
import Control.DeepSeq
import qualified Numeric.LinearAlgebra.Data as HM
import qualified Numeric.LinearAlgebra.HMatrix as HM
import qualified Data.Matrix as DM
import qualified LAoP.Matrix.Type as TM
import qualified Linear as LM
import qualified Linear.V as LM
import qualified Linear.Matrix as LM
import qualified Data.Vector as V
import LAoP.Utils hiding ((.))

randomMatrixHM :: Int -> Int -> Gen (HM.Matrix Double)
randomMatrixHM cols rows = do
  l <- vectorOf (cols * rows) arbitrary
  let lr = buildList l cols
      m  = HM.fromLists lr
  return m

randomMatrixDM :: Int -> Int -> Gen (DM.Matrix Double)
randomMatrixDM cols rows = do
  l <- vectorOf (cols * rows) arbitrary
  let lr = buildList l cols
      m  = DM.fromLists lr
  return m

randomMatrixTM :: forall a b . (TM.CountableDimensions a b, TM.FromListsN Double a b) => Gen (TM.Matrix Double a b)
randomMatrixTM = do
  let cols = fromInteger (natVal (Proxy :: Proxy (TM.Count a)))
      rows = fromInteger (natVal (Proxy :: Proxy (TM.Count b)))
  l <- vectorOf (cols * rows) arbitrary
  let lr = buildList l cols
      m  = TM.fromLists lr
  return m

randomMatrixLM :: forall a b . (KnownNat a, KnownNat b, LM.Dim a, LM.Dim b) => Gen (LM.V a (LM.V b Double))
randomMatrixLM = do
  let cols = fromInteger (natVal (Proxy :: Proxy a))
      rows = fromInteger (natVal (Proxy :: Proxy b))
  l <- vectorOf (cols * rows) (arbitrary :: Gen Double)
  let lr = buildList l cols
      v  = map (fromJust . LM.fromVector @b . V.fromList) lr
      (Just m)  = LM.fromVector @a (V.fromList v)
  return m

fromJust (Just a) = a

buildList [] _ = []
buildList l r  = take r l : buildList (drop r l) r

--
setupEnv1 = do
  -- HMatrix
  hm11 <- generate (resize 1 (randomMatrixHM 10 10))
  hm12 <- generate (resize 1 (randomMatrixHM 10 10))
  -- Matrix
  dm11 <- generate (resize 1 (randomMatrixDM 10 10))
  dm12 <- generate (resize 1 (randomMatrixDM 10 10))
  -- Linear
  lm11 <- generate (resize 1 (randomMatrixLM @10 @10))
  lm12 <- generate (resize 1 (randomMatrixLM @10 @10))
  -- Laop
  tm11 <- generate (resize 1 (randomMatrixTM @(Natural 0 10) @(Natural 0 10)))
  tm12 <- generate (resize 1 (randomMatrixTM @(Natural 0 10) @(Natural 0 10)))
  return (
   (hm11, hm12)
   ,(dm11, dm12)
   ,(lm11, lm12)
   ,(tm11, tm12)
         )

setupEnv2 = do
  -- HMatrix
  hm11 <- generate (resize 1 (randomMatrixHM 25 25))
  hm12 <- generate (resize 1 (randomMatrixHM 25 25))
  -- Matrix
  dm11 <- generate (resize 1 (randomMatrixDM 25 25))
  dm12 <- generate (resize 1 (randomMatrixDM 25 25))
  -- Linear
  lm11 <- generate (resize 1 (randomMatrixLM @25 @25))
  lm12 <- generate (resize 1 (randomMatrixLM @25 @25))
  -- Laop
  tm11 <- generate (resize 1 (randomMatrixTM @(Natural 0 25) @(Natural 0 25)))
  tm12 <- generate (resize 1 (randomMatrixTM @(Natural 0 25) @(Natural 0 25)))
  return (
   (hm11, hm12)
   ,(dm11, dm12)
   ,(lm11, lm12)
   ,(tm11, tm12)
         )

setupEnv3 = do
  -- HMatrix
  hm11 <- generate (resize 1 (randomMatrixHM 50 50))
  hm12 <- generate (resize 1 (randomMatrixHM 50 50))
  -- Matrix
  dm11 <- generate (resize 1 (randomMatrixDM 50 50))
  dm12 <- generate (resize 1 (randomMatrixDM 50 50))
  -- Linear
  lm11 <- generate (resize 1 (randomMatrixLM @50 @50))
  lm12 <- generate (resize 1 (randomMatrixLM @50 @50))
  -- Laop
  tm11 <- generate (resize 1 (randomMatrixTM @(Natural 0 50) @(Natural 0 50)))
  tm12 <- generate (resize 1 (randomMatrixTM @(Natural 0 50) @(Natural 0 50)))
  return (
   (hm11, hm12)
   ,(dm11, dm12)
   ,(lm11, lm12)
   ,(tm11, tm12)
         )

setupEnv4 = do
  -- HMatrix
  hm11 <- generate (resize 1 (randomMatrixHM 75 75))
  hm12 <- generate (resize 1 (randomMatrixHM 75 75))
  -- Matrix
  dm11 <- generate (resize 1 (randomMatrixDM 75 75))
  dm12 <- generate (resize 1 (randomMatrixDM 75 75))
  -- Linear
  lm11 <- generate (resize 1 (randomMatrixLM @75 @75))
  lm12 <- generate (resize 1 (randomMatrixLM @75 @75))
  -- Laop
  tm11 <- generate (resize 1 (randomMatrixTM @(Natural 0 75) @(Natural 0 75)))
  tm12 <- generate (resize 1 (randomMatrixTM @(Natural 0 75) @(Natural 0 75)))
  return (
   (hm11, hm12)
   ,(dm11, dm12)
   ,(lm11, lm12)
   ,(tm11, tm12)
         )

setupEnv5 = do
  -- HMatrix
  hm11 <- generate (resize 1 (randomMatrixHM 100 100))
  hm12 <- generate (resize 1 (randomMatrixHM 100 100))
  -- Matrix
  dm11 <- generate (resize 1 (randomMatrixDM 100 100))
  dm12 <- generate (resize 1 (randomMatrixDM 100 100))
  -- Linear
  lm11 <- generate (resize 1 (randomMatrixLM @100 @100))
  lm12 <- generate (resize 1 (randomMatrixLM @100 @100))
  -- Laop
  tm11 <- generate (resize 1 (randomMatrixTM @(Natural 0 100) @(Natural 0 100)))
  tm12 <- generate (resize 1 (randomMatrixTM @(Natural 0 100) @(Natural 0 100)))
  return (
   (hm11, hm12)
   ,(dm11, dm12)
   ,(lm11, lm12)
   ,(tm11, tm12)
         )

setupEnv6 = do
  -- HMatrix
  hm11 <- generate (resize 1 (randomMatrixHM 200 200))
  hm12 <- generate (resize 1 (randomMatrixHM 200 200))
  -- Matrix
  dm11 <- generate (resize 1 (randomMatrixDM 200 200))
  dm12 <- generate (resize 1 (randomMatrixDM 200 200))
  -- Linear
  lm11 <- generate (resize 1 (randomMatrixLM @200 @200))
  lm12 <- generate (resize 1 (randomMatrixLM @200 @200))
  -- Laop
  tm11 <- generate (resize 1 (randomMatrixTM @(Natural 0 200) @(Natural 0 200)))
  tm12 <- generate (resize 1 (randomMatrixTM @(Natural 0 200) @(Natural 0 200)))
  return (
   (hm11, hm12)
   ,(dm11, dm12)
   ,(lm11, lm12)
   ,(tm11, tm12)
         )

setupEnv7 = do
  -- HMatrix
  hm11 <- generate (resize 1 (randomMatrixHM 400 400))
  hm12 <- generate (resize 1 (randomMatrixHM 400 400))
  -- Matrix
  dm11 <- generate (resize 1 (randomMatrixDM 400 400))
  dm12 <- generate (resize 1 (randomMatrixDM 400 400))
  -- Linear
  lm11 <- generate (resize 1 (randomMatrixLM @400 @400))
  lm12 <- generate (resize 1 (randomMatrixLM @400 @400))
  -- Laop
  tm11 <- generate (resize 1 (randomMatrixTM @(Natural 0 400) @(Natural 0 400)))
  tm12 <- generate (resize 1 (randomMatrixTM @(Natural 0 400) @(Natural 0 400)))
  return (
   (hm11, hm12)
   ,(dm11, dm12)
   ,(lm11, lm12)
   ,(tm11, tm12)
         )

setupEnv8 = do
  -- HMatrix
  hm11 <- generate (resize 1 (randomMatrixHM 800 800))
  hm12 <- generate (resize 1 (randomMatrixHM 800 800))
  -- Matrix
  dm11 <- generate (resize 1 (randomMatrixDM 800 800))
  dm12 <- generate (resize 1 (randomMatrixDM 800 800))
  -- Linear
  lm11 <- generate (resize 1 (randomMatrixLM @800 @800))
  lm12 <- generate (resize 1 (randomMatrixLM @800 @800))
  -- Laop
  tm11 <- generate (resize 1 (randomMatrixTM @(Natural 0 800) @(Natural 0 800)))
  tm12 <- generate (resize 1 (randomMatrixTM @(Natural 0 800) @(Natural 0 800)))
  return (
   (hm11, hm12)
   ,(dm11, dm12)
   ,(lm11, lm12)
   ,(tm11, tm12)
         )

setupEnv9 = do
  -- HMatrix
  hm11 <- generate (resize 1 (randomMatrixHM 1600 1600))
  hm12 <- generate (resize 1 (randomMatrixHM 1600 1600))
  -- Matrix
  dm11 <- generate (resize 1 (randomMatrixDM 1600 1600))
  dm12 <- generate (resize 1 (randomMatrixDM 1600 1600))
  -- Linear
  lm11 <- generate (resize 1 (randomMatrixLM @1600 @1600))
  lm12 <- generate (resize 1 (randomMatrixLM @1600 @1600))
  -- Laop
  tm11 <- generate (resize 1 (randomMatrixTM @(Natural 0 1600) @(Natural 0 1600)))
  tm12 <- generate (resize 1 (randomMatrixTM @(Natural 0 1600) @(Natural 0 1600)))
  return (
   (hm11, hm12)
   ,(dm11, dm12)
   ,(lm11, lm12)
   ,(tm11, tm12)
         )

benchmark :: IO ()
benchmark = do
  print "Starting benchmarks..."
  defaultMain [
   env setupEnv1 $ \ ~((hm11, hm12), (dm11, dm12), (lm11, lm12), (tm11, tm12)) -> bgroup "Matrix Composition 1" [
   bgroup "10x10" [
     bench "hmatrix" $ nf (HM.mul hm11) hm12
   , bench "matrix" $ nf (DM.multStd2 dm11) dm12
   , bench "linear" $ nf (LM.!*! lm11) lm12
   , bench "laop" $ nf (TM.comp tm11) tm12
   ] ],
   env setupEnv2 $ \ ~((hm11, hm12), (dm11, dm12), (lm11, lm12), (tm11, tm12)) -> bgroup "Matrix Composition 2" [
   bgroup "25x25" [
     bench "hmatrix" $ nf (HM.mul hm11) hm12
   , bench "matrix" $ nf (DM.multStd2 dm11) dm12
   , bench "linear" $ nf (LM.!*! lm11) lm12
   , bench "laop" $ nf (TM.comp tm11) tm12
   ] ],
   env setupEnv3 $ \ ~((hm11, hm12), (dm11, dm12), (lm11, lm12), (tm11, tm12)) -> bgroup "Matrix Composition 3" [
   bgroup "50x50" [
     bench "hmatrix" $ nf (HM.mul hm11) hm12
   , bench "matrix" $ nf (DM.multStd2 dm11) dm12
   , bench "linear" $ nf (LM.!*! lm11) lm12
   , bench "laop" $ nf (TM.comp tm11) tm12
   ] ],
   env setupEnv4 $ \ ~((hm11, hm12), (dm11, dm12), (lm11, lm12), (tm11, tm12)) -> bgroup "Matrix Composition 4" [
   bgroup "75x75" [
     bench "hmatrix" $ nf (HM.mul hm11) hm12
   , bench "matrix" $ nf (DM.multStd2 dm11) dm12
   , bench "linear" $ nf (LM.!*! lm11) lm12
   , bench "laop" $ nf (TM.comp tm11) tm12
   ] ],
   env setupEnv5 $ \ ~((hm11, hm12), (dm11, dm12), (lm11, lm12), (tm11, tm12)) -> bgroup "Matrix Composition 5" [
   bgroup "100x100" [
     bench "hmatrix" $ nf (HM.mul hm11) hm12
   , bench "matrix" $ nf (DM.multStd2 dm11) dm12
   , bench "linear" $ nf (LM.!*! lm11) lm12
   , bench "laop" $ nf (TM.comp tm11) tm12
   ] ],
   env setupEnv6 $ \ ~((hm11, hm12), (dm11, dm12), (lm11, lm12), (tm11, tm12)) -> bgroup "Matrix Composition 6" [
   bgroup "200x200" [
     bench "hmatrix" $ nf (HM.mul hm11) hm12
   , bench "matrix" $ nf (DM.multStd2 dm11) dm12
   , bench "linear" $ nf (LM.!*! lm11) lm12
   , bench "laop" $ nf (TM.comp tm11) tm12
   ] ],
   env setupEnv7 $ \ ~((hm11, hm12), (dm11, dm12), (lm11, lm12), (tm11, tm12)) -> bgroup "Matrix Composition 7" [
   bgroup "400x400" [
     bench "hmatrix" $ nf (HM.mul hm11) hm12
   , bench "matrix" $ nf (DM.multStd2 dm11) dm12
   , bench "linear" $ nf (LM.!*! lm11) lm12
   , bench "laop" $ nf (TM.comp tm11) tm12
   ] ],
   env setupEnv8 $ \ ~((hm11, hm12), (dm11, dm12), (lm11, lm12), (tm11, tm12)) -> bgroup "Matrix Composition 8" [
   bgroup "800x800" [
     bench "hmatrix" $ nf (HM.mul hm11) hm12
   , bench "matrix" $ nf (DM.multStd2 dm11) dm12
   , bench "linear" $ nf (LM.!*! lm11) lm12
   , bench "laop" $ nf (TM.comp tm11) tm12
   ] ],
   env setupEnv9 $ \ ~((hm11, hm12), (dm11, dm12), (lm11, lm12), (tm11, tm12)) -> bgroup "Matrix Composition 9" [
   bgroup "1600x1600" [
     bench "hmatrix" $ nf (HM.mul hm11) hm12
   , bench "matrix" $ nf (DM.multStd2 dm11) dm12
   , bench "linear" $ nf (LM.!*! lm11) lm12
   , bench "laop" $ nf (TM.comp tm11) tm12
   ] ]
              ]


text:
module LexerSpec where

import Test.Hspec
import Data.Either
import Data.Complex
import Control.Exception (evaluate)

import qualified Lexer as L
import qualified Text.ParserCombinators.Parsec as PC
import qualified Text.ParserCombinators.Parsec.Error as PCError

instance Eq PC.ParseError where
         (==) x y = show x == show y

isLeft = null . rights . return

floatEqual x y = abs (x - y) < 1.0e-8
rightFloatEq y (Right x) = floatEqual x y
rightFloatEq _ _ = False

rightComplexEq (y1 :+ y2) (Right (x1 :+ x2)) = (floatEqual x1 y1) && (floatEqual x2 y2)
rightComplexEq _ _ = False

main = hspec $ do
  describe "rawParse" $ do
    it "parses identifiers" $ do
      L.parse L.identifier "a" `shouldBe` Right "a"
      L.parse L.identifier "_" `shouldBe` Right "_"
      L.parse L.identifier "_x9" `shouldBe` Right "_x9"
      L.parse L.identifier "ThisVariableIsExported" `shouldBe` Right "ThisVariableIsExported"
      L.parse L.identifier "αβ" `shouldBe` Right "αβ"
      L.parse L.identifier "foreach" `shouldBe` Right "foreach"
      L.parse L.identifier "for2" `shouldBe` Right "for2"
      L.parse L.identifier "breakage" `shouldBe` Right "breakage"

    it "parses keywords" $ do
      L.parse (L.keyword "for") "for" `shouldBe` Right ()
      L.parse (L.keyword "break") "break" `shouldBe` Right ()
      L.parse (L.keyword "continue") "continue" `shouldBe` Right ()
      L.parse (L.keyword "defer") "defer" `shouldBe` Right ()
      L.parse (L.keyword "for") "foreach" `shouldSatisfy` isLeft

    it "parses operators and delimiters" $ do
      L.parse (L.operator "+") "+" `shouldBe` Right ()
      L.parse (L.operator "+=") "+=" `shouldBe` Right ()
      L.parse (L.operator "&") "&" `shouldBe` Right ()
      L.parse (L.operator "&&") "&&" `shouldBe` Right ()
      L.parse (L.operator "&=") "&=" `shouldBe` Right ()
      L.parse (L.operator "&") "&&" `shouldSatisfy` isLeft

    it "parses integer literals" $ do
      L.parse L.integer "0" `shouldBe` Right 0
      L.parse L.integer "42" `shouldBe` Right 42
      L.parse L.integer "0600" `shouldBe` Right 0o600
      L.parse L.integer "0xBadFace3f" `shouldBe` Right 0xBadFace3f
      L.parse L.integer "894864549864587456180" `shouldBe` Right 894864549864587456180

    it "parses floating-point literals" $ do
      L.parse L.float "0." `shouldBe` Right 0.0
      L.parse L.float "72.40" `shouldBe` Right 72.4
      L.parse L.float "072.40" `shouldBe` Right 72.4
      L.parse L.float "2.71828" `shouldBe` Right 2.71828
      L.parse L.float "1.e+0" `shouldBe` Right 1.0
      L.parse L.float "6.67428e-11" `shouldBe` Right 6.67428e-11
      L.parse L.float "1E6" `shouldBe` Right 1.0e6
      L.parse L.float ".25" `shouldBe` Right 0.25
      L.parse L.float ".12345E+5" `shouldSatisfy` (rightFloatEq 0.12345e5)

    it "parses imaginary literals" $ do
      L.parse L.imaginary "0i" `shouldBe` Right (0 :+ 0)
      L.parse L.imaginary "011i" `shouldBe` Right (0 :+ 0o11)
      L.parse L.imaginary "0.i" `shouldBe` Right (0 :+ 0.0)
      L.parse L.imaginary "2.718128i" `shouldBe` Right (0 :+ 2.718128)
      L.parse L.imaginary "1.e+0i" `shouldBe` Right (0 :+ 1.0)
      L.parse L.imaginary "6.67428e-11i" `shouldBe` Right (0 :+ 6.67428e-11)
      L.parse L.imaginary "1E6i" `shouldBe` Right (0 :+ 1.0e6)
      L.parse L.imaginary ".25i" `shouldBe` Right (0 :+ 0.25)
      L.parse L.imaginary ".12345E+5i" `shouldSatisfy` rightComplexEq (0 :+ 0.12345e5)

    it "parses runes" $ do
      L.parse L.char "'a'" `shouldBe` Right 'a'
      L.parse L.char "'ä'" `shouldBe` Right 'ä'
      L.parse L.char "'本'" `shouldBe` Right '本'
      L.parse L.char "'\\t'" `shouldBe` Right '\t'
      L.parse L.char "'\\000'" `shouldBe` Right '\0'
      L.parse L.char "'\\007'" `shouldBe` Right '\7'
      L.parse L.char "'\\377'" `shouldBe` Right '\255'
      L.parse L.char "'\\x07'" `shouldBe` Right '\x7'
      L.parse L.char "'\\xff'" `shouldBe` Right '\xff'
      L.parse L.char "'\\u12e4'" `shouldBe` Right 'ዤ'
      L.parse L.char "'\\U0010FFFD'" `shouldBe` Right '\x10fffd'
      L.parse L.char "'\\\''" `shouldBe` Right '\''
      L.parse L.char "'aa'" `shouldSatisfy` isLeft
      L.parse L.char "'\\xa'" `shouldSatisfy` isLeft
      L.parse L.char "'\\0'" `shouldSatisfy` isLeft

    it "parses strings" $ do
      L.parse L.string "\"abc\"" `shouldBe` Right "abc"
      L.parse L.string "`abc`" `shouldBe` Right "abc"
      L.parse L.string "\"\\\\n\\n\\\\n\"" `shouldBe` Right "\\n\n\\n"
      L.parse L.string "`\\n\n\\n`" `shouldBe` Right "\\n\n\\n"
      L.parse L.string "\"\\n\"" `shouldBe` Right "\n"
      L.parse L.string "`\"`" `shouldBe` Right "\""
      L.parse L.string "\"\\\"\"" `shouldBe` Right "\""
      L.parse L.string "\"Hello, world!\\n\"" `shouldBe` Right "Hello, world!\n"
      L.parse L.string "\"日本語\"" `shouldBe` Right "日本語"
      L.parse L.string "\"\\u65e5本\\U00008a9e\"" `shouldBe` Right "\x65e5本\x00008a9e"
      L.parse L.string "\"\\xff\\u00FF\"" `shouldBe` Right "\xff\x00ff"
      L.parse L.string "\"\\uD800\"" `shouldBe` Right "\xd800"

max_stars_repo_name: waterlink/hgp

max_issues_repo_name: waterlink/hgp

max_forks_repo_name: waterlink/hgp


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

module School.Train.Test.ForwardPass
( forwardPassTest ) where

import Conduit ((.|), yield, liftIO, sinkList, await)
import Data.Either (isLeft)
import Numeric.LinearAlgebra (R)
import School.TestUtils (doCost, fromRight, randomAffineParams, randomMatrix, testState, weight1)
import School.Train.AppTrain (AppTrain)
import School.Train.ForwardPass
import School.Train.TrainState (TrainState(..), def)
import School.Types.PingPong (pingPongSingleton, reversePingPong, toPingPong)
import School.Types.Slinky (Slinky(..))
import School.Unit.Affine (affine)
import School.Unit.CostFunction (CostFunction)
import School.Unit.RecLin (recLin)
import School.Unit.Unit (Unit(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitBackward (BackwardStack)
import School.Unit.UnitParams (UnitParams(..))
import Test.Tasty (TestTree)
import Test.Tasty.QuickCheck hiding ((><))
import Test.Tasty.TH
import Test.QuickCheck.Monadic (assert, monadicIO)

weight :: CostFunction R (AppTrain R)
weight = weight1

prop_no_units :: Property
prop_no_units = monadicIO $ do
  let forward = forwardPass [] weight
  let source = yield ([], SNil)
  let pass = source .| forward .| sinkList
  result <- testState pass def
  assert $ isLeft result

prop_single_recLin :: (Positive Int) -> (Positive Int) -> Property
prop_single_recLin (Positive bSize) (Positive fSize) = monadicIO $ do
  let forward = forwardPass [recLin] weight
  input <- liftIO $ BatchActivation <$> (randomMatrix bSize fSize)
  let source = yield ([input], SNil)
  let pass = source .| forward .| await
  result <- testState pass def
  let out = apply recLin EmptyParams input
  let (cost, grad) = doCost weight out SNil
  let bParams = reversePingPong . paramList $ def
  let state = def { paramList = bParams }
  let stack = ([input], grad, cost)
  let check = Right $ (Just stack, state) :: Either String (Maybe (BackwardStack R), TrainState R)
  assert $ result == check

prop_aff_rl_aff_rl :: (Positive Int) -> (Positive Int) -> (Positive Int) -> (Positive Int) -> Property
prop_aff_rl_aff_rl (Positive b) (Positive f) (Positive h) (Positive o) = monadicIO $ do
  let units = [affine, recLin, affine, recLin]
  let forward = forwardPass units weight
  input <- liftIO $ BatchActivation <$> (randomMatrix b f)
  let source = yield ([input], SNil)
  let pass = source .| forward .| await
  params1 <- liftIO $ randomAffineParams f h
  params2 <- liftIO $ randomAffineParams h o
  let allParams = toPingPong [ params1
                             , EmptyParams
                             , params2
                             , EmptyParams
                             ]
  let paramList = fromRight (pingPongSingleton EmptyParams) allParams
  let initState = def { paramList }
  result <- testState pass initState
  let out1 = apply affine params1 input
  let out2 = apply recLin EmptyParams out1
  let out3 = apply affine params2 out2
  let out4 = apply recLin EmptyParams out3
  let (cost, grad) = doCost weight out4 SNil
  let bParams = reversePingPong paramList
  let state = def { paramList = bParams }
  let stack = ([out3, out2, out1, input], grad, cost)
  let check = Right $ (Just stack, state) :: Either String (Maybe (BackwardStack R), TrainState R)
  assert $ result == check

forwardPassTest:: TestTree
forwardPassTest = $(testGroupGenerator)

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
-- This file is auto-generated.  Do not edit directly.
-- |
--   Stability: Experimental
--
--   Generic interface to Blas using safe foreign calls.  Refer to the GHC documentation
--   for more information regarding appropriate use of safe and unsafe foreign calls.
--
--   The functions here are named in a similar fashion to the original Blas interface, with
--   the type-dependent letter(s) removed.  Some functions have been merged with
--   others to allow the interface to work on both real and complex numbers.  If you can't a
--   particular function, try looking for its corresponding complex equivalent (e.g.
--   @symv@ is a special case of 'hemv' applied to real numbers).
--
--   Note: although complex versions of @rot@ and @rotg@ exist in many implementations,
--   they are not part of the official Blas standard and therefore not included here.  If
--   you /really/ need them, submit a ticket so we can try to come up with a solution.
--
--   The documentation here is still incomplete.  Consult the
--   <http://netlib.org/blas/#_blas_routines official documentation> for more
--   information.
--
--   Notation:
--
--     * @⋅@ denotes dot product (without any conjugation).
--     * @*@ denotes complex conjugation.
--     * @⊤@ denotes transpose.
--     * @†@ denotes conjugate transpose (Hermitian conjugate).
--
--   Conventions:
--
--     * All scalars are denoted with lowercase Greek letters
--     * All vectors are denoted with lowercase Latin letters and are
--       assumed to be column vectors (unless transposed).
--     * All matrices are denoted with uppercase Latin letters.
{-# LANGUAGE FlexibleInstances, TypeFamilies #-}
module Blas.Generic.Safe (
    Numeric(..)
  , RealNumeric(..)
  , D.dsdot
  , S.sdsdot
  ) where
import Prelude (Floating, Double, Float, Int, IO)
import Data.Complex (Complex)
import Foreign (Ptr, Storable)
import Blas.Primitive.Types (Order, Transpose, Uplo, Diag, Side)
import qualified Blas.Specialized.Float.Safe as S
import qualified Blas.Specialized.Double.Safe as D
import qualified Blas.Specialized.ComplexFloat.Safe as C
import qualified Blas.Specialized.ComplexDouble.Safe as Z

-- | Blas operations that are applicable to real and complex numbers.
--
--   Instances are defined for the 4 types supported by Blas: the single- and
--   double-precision floating point types and their complex versions.
class (Floating a, Storable a) => Numeric a where

  -- | The corresponding real type of @a@.
  --
  --   In other words, @'RealType' ('Complex' a)@ is an alias for @a@.  For everything
  --   else, @'RealType' a@ is simply @a@.
  type RealType a :: *
  -- | Swap two vectors:
  --
  --   > (x, y) ← (y, x)
  swap :: Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Multiply a vector by a scalar.
  --
  --   > x ← α x
  scal :: Int
       -> a
       -> Ptr a
       -> Int
       -> IO ()

  -- | Copy a vector into another vector:
  --
  --   > y ← x
  copy :: Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Add a scalar-vector product to a vector.
  --
  --   > y ← α x + y
  axpy :: Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Calculate the bilinear dot product of two vectors:
  --
  --   > x ⋅ y ≡ ∑[i] x[i] y[i]
  dotu :: Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO a

  -- | Calculate the sesquilinear dot product of two vectors.
  --
  --   > x* ⋅ y ≡ ∑[i] x[i]* y[i]
  dotc :: Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO a

  -- | Calculate the Euclidean (L²) norm of a vector:
  --
  --   > ‖x‖₂ ≡ √(∑[i] x[i]²)
  nrm2 :: Int
       -> Ptr a
       -> Int
       -> IO (RealType a)

  -- | Calculate the Manhattan (L¹) norm, equal to the sum of the magnitudes of the elements:
  --
  --   > ‖x‖₁ = ∑[i] |x[i]|
  asum :: Int
       -> Ptr a
       -> Int
       -> IO (RealType a)

  -- | Calculate the index of the element with the maximum magnitude (absolute value).
  iamax :: Int
        -> Ptr a
        -> Int
        -> IO Int

  -- | Perform a general matrix-vector update.
  --
  --   > y ← α T(A) x + β y
  gemv :: Order
       -> Transpose
       -> Int
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform a general banded matrix-vector update.
  --
  --   > y ← α T(A) x + β y
  gbmv :: Order
       -> Transpose
       -> Int
       -> Int
       -> Int
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform a hermitian matrix-vector update.
  --
  --   > y ← α A x + β y
  hemv :: Order
       -> Uplo
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform a hermitian banded matrix-vector update.
  --
  --   > y ← α A x + β y
  hbmv :: Order
       -> Uplo
       -> Int
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform a hermitian packed matrix-vector update.
  --
  --   > y ← α A x + β y
  hpmv :: Order
       -> Uplo
       -> Int
       -> a
       -> Ptr a
       -> Ptr a
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> IO ()

  -- | Multiply a triangular matrix by a vector.
  --
  --   > x ← T(A) x
  trmv :: Order
       -> Uplo
       -> Transpose
       -> Diag
       -> Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Multiply a triangular banded matrix by a vector.
  --
  --   > x ← T(A) x
  tbmv :: Order
       -> Uplo
       -> Transpose
       -> Diag
       -> Int
       -> Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Multiply a triangular packed matrix by a vector.
  --
  --   > x ← T(A) x
  tpmv :: Order
       -> Uplo
       -> Transpose
       -> Diag
       -> Int
       -> Ptr a
       -> Ptr a
       -> Int
       -> IO ()

  -- | Multiply an inverse triangular matrix by a vector.
  --
  --   > x ← T(A⁻¹) x
  trsv :: Order
       -> Uplo
       -> Transpose
       -> Diag
       -> Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Multiply an inverse triangular banded matrix by a vector.
  --
  --   > x ← T(A⁻¹) x
  tbsv :: Order
       -> Uplo
       -> Transpose
       -> Diag
       -> Int
       -> Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Multiply an inverse triangular packed matrix by a vector.
  --
  --   > x ← T(A⁻¹) x
  tpsv :: Order
       -> Uplo
       -> Transpose
       -> Diag
       -> Int
       -> Ptr a
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform an unconjugated rank-1 update of a general matrix.
  --
  --   > A ← α x y⊤ + A
  geru :: Order
       -> Int
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform a conjugated rank-1 update of a general matrix.
  --
  --   > A ← α x y† + A
  gerc :: Order
       -> Int
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform a rank-1 update of a Hermitian matrix.
  --
  --   > A ← α x y† + A
  her :: Order
      -> Uplo
      -> Int
      -> RealType a
      -> Ptr a
      -> Int
      -> Ptr a
      -> Int
      -> IO ()

  -- | Perform a rank-1 update of a Hermitian packed matrix.
  --
  --   > A ← α x y† + A
  hpr :: Order
      -> Uplo
      -> Int
      -> RealType a
      -> Ptr a
      -> Int
      -> Ptr a
      -> IO ()

  -- | Perform a rank-2 update of a Hermitian matrix.
  --
  --   > A ← α x y† + y (α x)† + A
  her2 :: Order
       -> Uplo
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform a rank-2 update of a Hermitian packed matrix.
  --
  --   > A ← α x y† + y (α x)† + A
  hpr2 :: Order
       -> Uplo
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> IO ()

  -- | Perform a general matrix-matrix update.
  --
  --   > C ← α T(A) U(B) + β C
  gemm :: Order -- ^ Layout of all the matrices.
       -> Transpose -- ^ The operation @T@ to be applied to @A@.
       -> Transpose -- ^ The operation @U@ to be applied to @B@.
       -> Int -- ^ Number of rows of @T(A)@ and @C@.
       -> Int -- ^ Number of columns of @U(B)@ and @C@.
       -> Int -- ^ Number of columns of @T(A)@ and number of rows of @U(B)@.
       -> a -- ^ Scaling factor @α@ of the product.
       -> Ptr a -- ^ Pointer to a matrix @A@.
       -> Int -- ^ Stride of the major dimension of @A@.
       -> Ptr a -- ^ Pointer to a matrix @B@.
       -> Int -- ^ Stride of the major dimension of @B@.
       -> a -- ^ Scaling factor @β@ of the original @C@.
       -> Ptr a -- ^ Pointer to a mutable matrix @C@.
       -> Int -- ^ Stride of the major dimension of @C@.
       -> IO ()

  -- | Perform a symmetric matrix-matrix update.
  --
  --   > C ← α A B + β C    or    C ← α B A + β C
  --
  --   where @A@ is symmetric.  The matrix @A@ must be in an unpacked format, although the
  --   routine will only access half of it as specified by the @'Uplo'@ argument.
  symm :: Order -- ^ Layout of all the matrices.
       -> Side -- ^ Side that @A@ appears in the product.
       -> Uplo -- ^ The part of @A@ that is used.
       -> Int -- ^ Number of rows of @C@.
       -> Int -- ^ Number of columns of @C@.
       -> a -- ^ Scaling factor @α@ of the product.
       -> Ptr a -- ^ Pointer to a symmetric matrix @A@.
       -> Int -- ^ Stride of the major dimension of @A@.
       -> Ptr a -- ^ Pointer to a matrix @B@.
       -> Int -- ^ Stride of the major dimension of @B@.
       -> a -- ^ Scaling factor @α@ of the original @C@.
       -> Ptr a -- ^ Pointer to a mutable matrix @C@.
       -> Int -- ^ Stride of the major dimension of @C@.
       -> IO ()

  -- | Perform a Hermitian matrix-matrix update.
  --
  --   > C ← α A B + β C    or    C ← α B A + β C
  --
  --   where @A@ is Hermitian.  The matrix @A@ must be in an unpacked format, although the
  --   routine will only access half of it as specified by the @'Uplo'@ argument.
  hemm :: Order -- ^ Layout of all the matrices.
       -> Side -- ^ Side that @A@ appears in the product.
       -> Uplo -- ^ The part of @A@ that is used.
       -> Int -- ^ Number of rows of @C@.
       -> Int -- ^ Number of columns of @C@.
       -> a -- ^ Scaling factor @α@ of the product.
       -> Ptr a -- ^ Pointer to a Hermitian matrix @A@.
       -> Int -- ^ Stride of the major dimension of @A@.
       -> Ptr a -- ^ Pointer to a matrix @B@.
       -> Int -- ^ Stride of the major dimension of @B@.
       -> a -- ^ Scaling factor @α@ of the original @C@.
       -> Ptr a -- ^ Pointer to a mutable matrix @C@.
       -> Int -- ^ Stride of the major dimension of @C@.
       -> IO ()

  -- | Perform a symmetric rank-k update.
  --
  --   > C ← α A A⊤ + β C    or    C ← α A⊤ A + β C
  syrk :: Order
       -> Uplo
       -> Transpose
       -> Int
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform a Hermitian rank-k update.
  --
  --   > C ← α A A† + β C    or    C ← α A† A + β C
  herk :: Order
       -> Uplo
       -> Transpose
       -> Int
       -> Int
       -> RealType a
       -> Ptr a
       -> Int
       -> RealType a
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform a symmetric rank-2k update.
  --
  --   > C ← α A B⊤ + α* B A⊤ + β C    or    C ← α A⊤ B + α* B⊤ A + β C
  syr2k :: Order
        -> Uplo
        -> Transpose
        -> Int
        -> Int
        -> a
        -> Ptr a
        -> Int
        -> Ptr a
        -> Int
        -> a
        -> Ptr a
        -> Int
        -> IO ()

  -- | Perform a Hermitian rank-2k update.
  --
  --   > C ← α A B† + α* B A† + β C    or    C ← α A† B + α* B† A + β C
  her2k :: Order
        -> Uplo
        -> Transpose
        -> Int
        -> Int
        -> a
        -> Ptr a
        -> Int
        -> Ptr a
        -> Int
        -> RealType a
        -> Ptr a
        -> Int
        -> IO ()

  -- | Perform a triangular matrix-matrix multiplication.
  --
  --   > B ← α T(A) B    or    B ← α B T(A)
  --
  --   where @A@ is triangular.
  trmm :: Order
       -> Side
       -> Uplo
       -> Transpose
       -> Diag
       -> Int
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

  -- | Perform an inverse triangular matrix-matrix multiplication.
  --
  --   > B ← α T(A⁻¹) B    or    B ← α B T(A⁻¹)
  --
  --   where @A@ is triangular.
  trsm :: Order
       -> Side
       -> Uplo
       -> Transpose
       -> Diag
       -> Int
       -> Int
       -> a
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> IO ()

-- | Blas operations that are only applicable to real numbers.
class Numeric a => RealNumeric a where

  -- | Generate a Givens rotation. (Only available for real floating-point types.)
  rotg :: Ptr a
       -> Ptr a
       -> Ptr a
       -> Ptr a
       -> IO ()

  -- | Generate a modified Givens rotation. (Only available for real floating-point
  --   types.)
  rotmg :: Ptr a
        -> Ptr a
        -> Ptr a
        -> a
        -> Ptr a
        -> IO ()

  -- | Apply a Givens rotation. (Only available for real floating-point types.)
  rot :: Int
      -> Ptr a
      -> Int
      -> Ptr a
      -> Int
      -> a
      -> a
      -> IO ()

  -- | Apply a modified Givens rotation. (Only available for real floating-point types.)
  rotm :: Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> Int
       -> Ptr a
       -> IO ()

instance Numeric Float where
  type RealType Float = Float
  swap = S.swap
  scal = S.scal
  copy = S.copy
  axpy = S.axpy
  dotu = S.dotu
  dotc = S.dotc
  nrm2 = S.nrm2
  asum = S.asum
  iamax = S.iamax
  gemv = S.gemv
  gbmv = S.gbmv
  hemv = S.hemv
  hbmv = S.hbmv
  hpmv = S.hpmv
  trmv = S.trmv
  tbmv = S.tbmv
  tpmv = S.tpmv
  trsv = S.trsv
  tbsv = S.tbsv
  tpsv = S.tpsv
  geru = S.geru
  gerc = S.gerc
  her = S.her
  hpr = S.hpr
  her2 = S.her2
  hpr2 = S.hpr2
  gemm = S.gemm
  symm = S.symm
  hemm = S.hemm
  syrk = S.syrk
  herk = S.herk
  syr2k = S.syr2k
  her2k = S.her2k
  trmm = S.trmm
  trsm = S.trsm

instance Numeric Double where
  type RealType Double = Double
  swap = D.swap
  scal = D.scal
  copy = D.copy
  axpy = D.axpy
  dotu = D.dotu
  dotc = D.dotc
  nrm2 = D.nrm2
  asum = D.asum
  iamax = D.iamax
  gemv = D.gemv
  gbmv = D.gbmv
  hemv = D.hemv
  hbmv = D.hbmv
  hpmv = D.hpmv
  trmv = D.trmv
  tbmv = D.tbmv
  tpmv = D.tpmv
  trsv = D.trsv
  tbsv = D.tbsv
  tpsv = D.tpsv
  geru = D.geru
  gerc = D.gerc
  her = D.her
  hpr = D.hpr
  her2 = D.her2
  hpr2 = D.hpr2
  gemm = D.gemm
  symm = D.symm
  hemm = D.hemm
  syrk = D.syrk
  herk = D.herk
  syr2k = D.syr2k
  her2k = D.her2k
  trmm = D.trmm
  trsm = D.trsm

instance Numeric (Complex Float) where
  type RealType (Complex Float) = Float
  swap = C.swap
  scal = C.scal
  copy = C.copy
  axpy = C.axpy
  dotu = C.dotu
  dotc = C.dotc
  nrm2 = C.nrm2
  asum = C.asum
  iamax = C.iamax
  gemv = C.gemv
  gbmv = C.gbmv
  hemv = C.hemv
  hbmv = C.hbmv
  hpmv = C.hpmv
  trmv = C.trmv
  tbmv = C.tbmv
  tpmv = C.tpmv
  trsv = C.trsv
  tbsv = C.tbsv
  tpsv = C.tpsv
  geru = C.geru
  gerc = C.gerc
  her = C.her
  hpr = C.hpr
  her2 = C.her2
  hpr2 = C.hpr2
  gemm = C.gemm
  symm = C.symm
  hemm = C.hemm
  syrk = C.syrk
  herk = C.herk
  syr2k = C.syr2k
  her2k = C.her2k
  trmm = C.trmm
  trsm = C.trsm

instance Numeric (Complex Double) where
  type RealType (Complex Double) = Double
  swap = Z.swap
  scal = Z.scal
  copy = Z.copy
  axpy = Z.axpy
  dotu = Z.dotu
  dotc = Z.dotc
  nrm2 = Z.nrm2
  asum = Z.asum
  iamax = Z.iamax
  gemv = Z.gemv
  gbmv = Z.gbmv
  hemv = Z.hemv
  hbmv = Z.hbmv
  hpmv = Z.hpmv
  trmv = Z.trmv
  tbmv = Z.tbmv
  tpmv = Z.tpmv
  trsv = Z.trsv
  tbsv = Z.tbsv
  tpsv = Z.tpsv
  geru = Z.geru
  gerc = Z.gerc
  her = Z.her
  hpr = Z.hpr
  her2 = Z.her2
  hpr2 = Z.hpr2
  gemm = Z.gemm
  symm = Z.symm
  hemm = Z.hemm
  syrk = Z.syrk
  herk = Z.herk
  syr2k = Z.syr2k
  her2k = Z.her2k
  trmm = Z.trmm
  trsm = Z.trsm
instance RealNumeric Float where
  rotg = S.rotg
  rotmg = S.rotmg
  rot = S.rot
  rotm = S.rotm

instance RealNumeric Double where
  rotg = D.rotg
  rotmg = D.rotmg
  rot = D.rot
  rotm = D.rotm

max_stars_repo_name: Rufflewind/blas-hs

max_issues_repo_name: Rufflewind/blas-hs

max_forks_repo_name: Rufflewind/blas-hs


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE ScopedTypeVariables   #-}

module GrenadeExtras.Zip (
  Zip (..),
) where

import           Data.Serialize (Serialize, put, get)

--import           Data.Singletons
import           GHC.TypeLits

import           Grenade.Core (UpdateLayer, Layer, Shape(D1), S(S1D), Gradient,
                               runForwards, runBackwards, runUpdate, createRandom, Tape)
import           GrenadeExtras.GradNorm (GradNorm, normSquared)

import           Numeric.LinearAlgebra.Static ((#), split, R) -- row, (===), splitRows, unrow

data Zip :: Shape -> Shape -> * -> Shape -> Shape -> * -> * where
  Zip :: x -> y -> Zip ix ox x iy oy y

instance (Show x, Show y) => Show (Zip ix ox x iy oy y) where
  show (Zip x y) = "Zip\n" ++ show x ++ "\n" ++ show y

instance (UpdateLayer x, UpdateLayer y) => UpdateLayer (Zip ix ox x iy oy y) where
  type Gradient (Zip ix ox x iy oy y) = (Gradient x, Gradient y)
  runUpdate lr (Zip x y) (x', y') = Zip (runUpdate lr x x') (runUpdate lr y y')
  createRandom = Zip <$> createRandom <*> createRandom

instance ( Layer x ('D1 p) ('D1 m)
         , Layer y ('D1 q) ('D1 n)
         , KnownNat r
         , KnownNat p
         , KnownNat q
         , r ~ (p + q)
         , q ~ (r - p)
         , (p <=? r) ~ 'True
         , KnownNat o
         , KnownNat m
         , KnownNat n
         , o ~ (m + n)
         , n ~ (o - m)
         , (m <=? o) ~ 'True
         ) => Layer (Zip ('D1 p) ('D1 m) x ('D1 q) ('D1 n) y) ('D1 r) ('D1 o) where
   type Tape (Zip ('D1 p) ('D1 m) x ('D1 q) ('D1 n) y) ('D1 r) ('D1 o) = (Tape x ('D1 p) ('D1 m), Tape y ('D1 q) ('D1 n))

   runForwards (Zip x y) (S1D input) =
     let (inputX :: R p, inputY :: R q) = split input
         (xT, xOut :: S ('D1 m)) = runForwards x (S1D inputX)
         (yT, yOut :: S ('D1 n)) = runForwards y (S1D inputY)
     in case (xOut, yOut) of
         (S1D xOut', S1D yOut') ->
             ((xT, yT), S1D (xOut' # yOut'))

   runBackwards (Zip x y) (xTape, yTape) (S1D o) =
     let (ox :: R m , oy :: R n) = split o
         (x', xB :: S ('D1 p)) = runBackwards x xTape (S1D ox)
         (y', yB :: S ('D1 q)) = runBackwards y yTape (S1D oy)
      in case (xB, yB) of
           (S1D xB', S1D yB') ->
             ((x', y'), S1D (xB' # yB'))

instance (Serialize a, Serialize b) => Serialize (Zip sia sa a sib sb b) where
  put (Zip a b) = put a *> put b
  get = Zip <$> get <*> get

instance (GradNorm a, GradNorm b) => GradNorm (Zip sia sa a sib sb b) where
  normSquared (Zip a b) = normSquared a + normSquared b

instance (Num a, Num b) => Num (Zip sia sa a sib sb b) where
  (Zip w1 m1) + (Zip w2 m2) = Zip (w1+w2) (m1+m2)
  (Zip w1 m1) - (Zip w2 m2) = Zip (w1-w2) (m1-m2)
  (Zip w1 m1) * (Zip w2 m2) = Zip (w1*w2) (m1*m2)
  abs    (Zip w m) = Zip (abs w) (abs m)
  signum (Zip w m) = Zip (signum w) (signum m)
  negate (Zip w m) = Zip (negate w) (negate m)
  fromInteger i = Zip (fromInteger i) (fromInteger i)

max_stars_repo_name: helq/haskell-binary-classification

max_issues_repo_name: helq/haskell-binary-classification

max_forks_repo_name: helq/haskell-binary-classification