text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Data.Eigen.Matrix (
    -- * Matrix type
    -- | Matrix aliases follows Eigen naming convention
    Matrix(..),
    MatrixXf,
    MatrixXd,
    MatrixXcf,
    MatrixXcd,
    I.Elem,
    I.CComplex,
    valid,
    -- * Matrix conversions
    fromList,
    toList,
    fromFlatList,
    toFlatList,
    generate,
    -- * Standard matrices and special cases
    empty,
    null,
    square,
    zero,
    ones,
    identity,
    constant,
    random,
    -- * Accessing matrix data
    cols,
    rows,
    dims,
    (!),
    coeff,
    unsafeCoeff,
    col,
    row,
    block,
    topRows,
    bottomRows,
    leftCols,
    rightCols,
    -- * Matrix properties
    sum,
    prod,
    mean,
    minCoeff,
    maxCoeff,
    trace,
    norm,
    squaredNorm,
    blueNorm,
    hypotNorm,
    determinant,
    -- * Generic reductions
    fold,
    fold',
    ifold,
    ifold',
    fold1,
    fold1',
    -- * Boolean reductions
    all,
    any,
    count,
    -- * Basic matrix algebra
    add,
    sub,
    mul,
    -- * Mapping over elements
    map,
    imap,
    filter,
    ifilter,
    -- * Matrix transformations
    diagonal,
    transpose,
    inverse,
    adjoint,
    conjugate,
    normalize,
    modify,
    convert,
    TriangularMode(..),
    triangularView,
    lowerTriangle,
    upperTriangle,
    -- * Matrix serialization
    encode,
    decode,
    -- * Mutable matrices
    thaw,
    freeze,
    unsafeThaw,
    unsafeFreeze,
    -- * Raw pointers
    unsafeWith,
) where

import qualified Prelude as P
import qualified Data.List as L
import Prelude hiding (null, sum, all, any, map, filter)
import Data.Tuple
import Data.Complex hiding (conjugate)
import Data.Binary hiding (encode, decode)
import qualified Data.Binary as B
import Foreign.Ptr
import Foreign.C.Types
import Foreign.C.String
import Foreign.Storable
import Foreign.Marshal.Alloc
import Text.Printf
import Control.Monad
import Control.Monad.ST
import Control.Monad.Primitive
#if __GLASGOW_HASKELL__ >= 710
#else
import Control.Applicative hiding (empty)
#endif
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM
import qualified Data.Eigen.Internal as I
import qualified Data.Eigen.Matrix.Mutable as M
import qualified Data.ByteString.Lazy as BSL

-- | Matrix to be used in pure computations, uses column major memory layout, features copy-free FFI with C++ <http://eigen.tuxfamily.org Eigen> library.

data Matrix a b where
    Matrix :: I.Elem a b => !Int -> !Int -> !(VS.Vector b) -> Matrix a b

-- | Alias for single precision matrix
type MatrixXf = Matrix Float CFloat
-- | Alias for double precision matrix
type MatrixXd = Matrix Double CDouble
-- | Alias for single previsiom matrix of complex numbers
type MatrixXcf = Matrix (Complex Float) (I.CComplex CFloat)
-- | Alias for double prevision matrix of complex numbers
type MatrixXcd = Matrix (Complex Double) (I.CComplex CDouble)

-- | Pretty prints the matrix
instance (I.Elem a b, Show a) => Show (Matrix a b) where
    show m@(Matrix rows cols _) = concat [
        "Matrix ", show rows, "x", show cols,
        "\n", L.intercalate "\n" $ P.map (L.intercalate "\t" . P.map show) $ toList m, "\n"]

-- | Basic matrix math exposed through Num instance: @(*)@, @(+)@, @(-)@, `fromInteger`, `signum`, `abs`, `negate`
instance I.Elem a b => Num (Matrix a b) where
    (*) = mul
    (+) = add
    (-) = sub
    fromInteger = constant 1 1 . fromInteger
    signum = map signum
    abs = map abs
    negate = map negate

-- | Matrix binary serialization
instance I.Elem a b => Binary (Matrix a b) where
    put (Matrix rows cols vals) = do
        put $ I.magicCode (undefined :: b)
        put rows
        put cols
        put vals

    get = do
        get >>= (`when` fail "wrong matrix type") . (/= I.magicCode (undefined :: b))
        Matrix <$> get <*> get <*> get

-- | Encode the matrix as a lazy byte string
encode :: I.Elem a b => Matrix a b -> BSL.ByteString
encode = B.encode

-- | Decode matrix from the lazy byte string
decode :: I.Elem a b => BSL.ByteString -> Matrix a b
decode = B.decode

-- | Empty 0x0 matrix
{-# INLINE empty #-}
empty :: I.Elem a b => Matrix a b
empty = Matrix 0 0 VS.empty

-- | Is matrix empty?
{-# INLINE null #-}
null :: I.Elem a b => Matrix a b -> Bool
null (Matrix rows cols _) = rows == 0 && cols == 0

-- | Is matrix square?
{-# INLINE square #-}
square :: I.Elem a b => Matrix a b -> Bool
square (Matrix rows cols _) = rows == cols

-- | Matrix where all coeffs are filled with given value
{-# INLINE constant #-}
constant :: I.Elem a b => Int -> Int -> a -> Matrix a b
constant rows cols val = Matrix rows cols $ VS.replicate (rows * cols) (I.cast val)

-- | Matrix where all coeff are 0
{-# INLINE zero #-}
zero :: I.Elem a b => Int -> Int -> Matrix a b
zero rows cols = constant rows cols 0

-- | Matrix where all coeff are 1
{-# INLINE ones #-}
ones :: I.Elem a b => Int -> Int -> Matrix a b
ones rows cols = constant rows cols 1

-- | The identity matrix (not necessarily square).
identity :: I.Elem a b => Int -> Int -> Matrix a b
identity rows cols = I.performIO $ do
    m <- M.new rows cols
    I.call $ M.unsafeWith m I.identity
    unsafeFreeze m

-- | The random matrix of a given size
random :: I.Elem a b => Int -> Int -> IO (Matrix a b)
random rows cols = do
    m <- M.new rows cols
    I.call $ M.unsafeWith m I.random
    unsafeFreeze m

-- | Number of rows for the matrix
{-# INLINE rows #-}
rows :: I.Elem a b => Matrix a b -> Int
rows (Matrix rows _ _) = rows

-- | Number of columns for the matrix
{-# INLINE cols #-}
cols :: I.Elem a b => Matrix a b -> Int
cols (Matrix _ cols _) = cols

-- | Mtrix size as (rows, cols) pair
{-# INLINE dims #-}
dims :: I.Elem a b => Matrix a b -> (Int, Int)
dims (Matrix rows cols _) = (rows, cols)

-- | Matrix coefficient at specific row and col
{-# INLINE (!) #-}
(!) :: I.Elem a b => Matrix a b -> (Int, Int) -> a
(!) m (row,col) = coeff row col m

-- | Matrix coefficient at specific row and col
{-# INLINE coeff #-}
coeff :: I.Elem a b => Int -> Int -> Matrix a b -> a
coeff row col m@(Matrix rows cols _)
    | not (valid m) = error "matrix is not valid"
    | row < 0 || row >= rows = error $ printf "Matrix.coeff: row %d is out of bounds [0..%d)" row rows
    | col < 0 || col >= cols = error $ printf "Matrix.coeff: col %d is out of bounds [0..%d)" col cols
    | otherwise = unsafeCoeff row col m

-- | Unsafe version of coeff function. No bounds check performed so SEGFAULT possible
{-# INLINE unsafeCoeff #-}
unsafeCoeff :: I.Elem a b => Int -> Int -> Matrix a b -> a
unsafeCoeff row col (Matrix rows _ vals) = I.cast $ VS.unsafeIndex vals $ col * rows + row

-- | List of coefficients for the given col
{-# INLINE col #-}
col :: I.Elem a b => Int -> Matrix a b -> [a]
col c m@(Matrix rows _ _) = [coeff r c m | r <- [0..pred rows]]

-- | List of coefficients for the given row
{-# INLINE row #-}
row :: I.Elem a b => Int -> Matrix a b -> [a]
row r m@(Matrix _ cols _) = [coeff r c m | c <- [0..pred cols]]

-- | Extract rectangular block from matrix defined by startRow startCol blockRows blockCols
block :: I.Elem a b => Int -> Int -> Int -> Int -> Matrix a b -> Matrix a b
block startRow startCol blockRows blockCols m =
    generate blockRows blockCols $ \row col ->
        coeff (startRow + row) (startCol + col) m

-- | Verify matrix dimensions and memory layout
{-# INLINE valid #-}
valid :: I.Elem a b => Matrix a b -> Bool
valid (Matrix rows cols vals) = rows >= 0 && cols >= 0 && VS.length vals == rows * cols

-- | The maximum coefficient of the matrix
{-# INLINE maxCoeff #-}
maxCoeff :: (I.Elem a b, Ord a) => Matrix a b -> a
maxCoeff = fold1' max

-- | The minimum coefficient of the matrix
{-# INLINE minCoeff #-}
minCoeff :: (I.Elem a b, Ord a) => Matrix a b -> a
minCoeff = fold1' min

-- | Top @N@ rows of matrix
{-# INLINE topRows #-}
topRows :: I.Elem a b => Int -> Matrix a b -> Matrix a b
topRows n m@(Matrix _ cols _) = block 0 0 n cols m

-- | Bottom @N@ rows of matrix
{-# INLINE bottomRows #-}
bottomRows :: I.Elem a b => Int -> Matrix a b -> Matrix a b
bottomRows n m@(Matrix rows cols _) = block (rows - n) 0 n cols m

-- | Left @N@ columns of matrix
{-# INLINE leftCols #-}
leftCols :: I.Elem a b => Int -> Matrix a b -> Matrix a b
leftCols n m@(Matrix rows _ _) = block 0 0 rows n m

-- | Right @N@ columns of matrix
{-# INLINE rightCols #-}
rightCols :: I.Elem a b => Int -> Matrix a b -> Matrix a b
rightCols n m@(Matrix rows cols _) = block 0 (cols - n) rows n m

-- | Construct matrix from a list of rows, column count is detected as maximum row length. Missing values are filled with 0
fromList :: I.Elem a b => [[a]] -> Matrix a b
fromList list = Matrix rows cols vals where
    rows = length list
    cols = L.foldl' max 0 $ P.map length list
    vals = VS.create $ do
        vm <- VSM.replicate (rows * cols) (I.cast (0 `asTypeOf` (head (head list))))
        forM_ (zip [0..] list) $ \(row, vals) ->
            forM_ (zip [0..] vals) $ \(col, val) ->
                VSM.write vm (col * rows + row) (I.cast val)
        return vm

-- | Convert matrix to a list of rows
toList :: I.Elem a b => Matrix a b -> [[a]]
toList m@(Matrix rows cols vals)
    | not (valid m) = error "matrix is not valid"
    | otherwise = [[I.cast $ vals `VS.unsafeIndex` (col * rows + row) | col <- [0..pred cols]] | row <- [0..pred rows]]

-- | Build matrix of given dimensions and values from given list split on rows. Invalid list length results in error.
fromFlatList :: I.Elem a b => Int -> Int -> [a] -> Matrix a b
fromFlatList rows cols list
    | not (rows * cols == (length list)) = error $ concat ["cannot construct ", show rows, "x", show cols, " matrix from ", show $ length list, " values"]
    | otherwise = Matrix rows cols vals where
        vals = VS.create $ do
            vm <- VSM.replicate (rows * cols) (I.cast (0 `asTypeOf` (head list)))
            forM_ (zip [(col * rows + row) | row <- [0..pred rows], col <- [0..pred cols]] list) $ \(idx, val) ->
                VSM.write vm idx (I.cast val)
            return vm

-- | Convert matrix to a list by concatenating rows
toFlatList :: I.Elem a b => Matrix a b -> [a]
toFlatList m@(Matrix rows cols vals)
    | not (valid m) = error "matrix is not valid"
    | otherwise = [I.cast $ vals `VS.unsafeIndex` (col * rows + row) | row <- [0..pred rows], col <- [0..pred cols]]

-- | [generate rows cols (λ row col -> val)]
--
-- Create matrix using generator function @λ row col -> val@
--
generate :: I.Elem a b => Int -> Int -> (Int -> Int -> a) -> Matrix a b
generate rows cols f = Matrix rows cols $ VS.create $ do
    vals <- VSM.new (rows * cols)
    forM_ [0..pred rows] $ \row ->
        forM_ [0..pred cols] $ \col ->
            VSM.write vals (col * rows + row) (I.cast $ f row col)
    return vals

-- | The sum of all coefficients of the matrix
sum :: I.Elem a b => Matrix a b -> a
sum = _prop I.sum

-- | The product of all coefficients of the matrix
prod :: I.Elem a b => Matrix a b -> a
prod = _prop I.prod

-- | The mean of all coefficients of the matrix
mean :: I.Elem a b => Matrix a b -> a
mean = _prop I.mean

-- | The trace of a matrix is the sum of the diagonal coefficients and can also be computed as sum (diagonal m)
trace :: I.Elem a b => Matrix a b -> a
trace = _prop I.trace

-- | Applied to a predicate and a matrix, all determines if all elements of the matrix satisfies the predicate
all :: I.Elem a b => (a -> Bool) -> Matrix a b -> Bool
all f = VS.all (f . I.cast) . _vals

-- | Applied to a predicate and a matrix, any determines if any element of the matrix satisfies the predicate
any :: I.Elem a b => (a -> Bool) -> Matrix a b -> Bool
any f = VS.any (f . I.cast) . _vals

-- | Returns the number of coefficients in a given matrix that evaluate to true
count :: I.Elem a b => (a -> Bool) -> Matrix a b -> Int
count f = VS.foldl' (\n x -> if f (I.cast x) then succ n else n) 0 . _vals

{-| For vectors, the l2 norm, and for matrices the Frobenius norm.
    In both cases, it consists in the square root of the sum of the square of all the matrix entries.
    For vectors, this is also equals to the square root of the dot product of this with itself.
-}
norm :: I.Elem a b => Matrix a b -> a
norm = _prop I.norm

-- | For vectors, the squared l2 norm, and for matrices the Frobenius norm. In both cases, it consists in the sum of the square of all the matrix entries. For vectors, this is also equals to the dot product of this with itself.
squaredNorm :: I.Elem a b => Matrix a b -> a
squaredNorm = _prop I.squaredNorm

-- | The l2 norm of the matrix using the Blue's algorithm. A Portable Fortran Program to Find the Euclidean Norm of a Vector, ACM TOMS, Vol 4, Issue 1, 1978.
blueNorm :: I.Elem a b => Matrix a b -> a
blueNorm = _prop I.blueNorm

-- | The l2 norm of the matrix avoiding undeflow and overflow. This version use a concatenation of hypot calls, and it is very slow.
hypotNorm :: I.Elem a b => Matrix a b -> a
hypotNorm = _prop I.hypotNorm

-- | The determinant of the matrix
determinant :: I.Elem a b => Matrix a b -> a
determinant m
    | square m = _prop I.determinant m
    | otherwise = error "Matrix.determinant: non-square matrix"

-- | Adding two matrices by adding the corresponding entries together. You can use @(+)@ function as well.
add :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b
add m1 m2
    | dims m1 == dims m2 = _binop const I.add m1 m2
    | otherwise = error "Matrix.add: matrices should have the same size"

-- | Subtracting two matrices by subtracting the corresponding entries together. You can use @(-)@ function as well.
sub :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b
sub m1 m2
    | dims m1 == dims m2 = _binop const I.sub m1 m2
    | otherwise = error "Matrix.add: matrices should have the same size"

-- | Matrix multiplication. You can use @(*)@ function as well.
mul :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b
mul m1 m2
    | cols m1 == rows m2 = _binop (\(rows, _) (_, cols) -> (rows, cols)) I.mul m1 m2
    | otherwise = error "Matrix.mul: number of columns for lhs matrix should be the same as number of rows for rhs matrix"

{- | Apply a given function to each element of the matrix.

Here is an example how to implement scalar matrix multiplication:

>>> let a = fromList [[1,2],[3,4]] :: MatrixXf

>>> a
Matrix 2x2
1.0 2.0
3.0 4.0

>>> map (*10) a
Matrix 2x2
10.0    20.0
30.0    40.0

-}
map :: I.Elem a b => (a -> a) -> Matrix a b -> Matrix a b
map f (Matrix rows cols vals) = Matrix rows cols (VS.map (I.cast . f . I.cast) vals)

{- | Apply a given function to each element of the matrix.

Here is an example how upper triangular matrix can be implemented:

>>> let a = fromList [[1,2,3],[4,5,6],[7,8,9]] :: MatrixXf

>>> a
Matrix 3x3
1.0 2.0 3.0
4.0 5.0 6.0
7.0 8.0 9.0

>>> imap (\row col val -> if row <= col then val else 0) a
Matrix 3x3
1.0 2.0 3.0
0.0 5.0 6.0
0.0 0.0 9.0

-}

imap :: I.Elem a b => (Int -> Int -> a -> a) -> Matrix a b -> Matrix a b
imap f (Matrix rows cols vals) = Matrix rows cols (VS.imap (\n -> let (c, r) = divMod n rows in I.cast . f r c . I.cast) vals)

data TriangularMode
    -- | View matrix as a lower triangular matrix.
    = Lower
    -- | View matrix as an upper triangular matrix.
    | Upper
    -- | View matrix as a lower triangular matrix with zeros on the diagonal.
    | StrictlyLower
    -- | View matrix as an upper triangular matrix with zeros on the diagonal.
    | StrictlyUpper
    -- | View matrix as a lower triangular matrix with ones on the diagonal.
    | UnitLower
    -- | View matrix as an upper triangular matrix with ones on the diagonal.
    | UnitUpper deriving (Eq, Enum, Show, Read)

-- | Triangular view extracted from the current matrix
triangularView :: I.Elem a b => TriangularMode -> Matrix a b -> Matrix a b
triangularView Lower         = imap $ \row col val -> case compare row col of { LT -> 0; _ -> val }
triangularView Upper         = imap $ \row col val -> case compare row col of { GT -> 0; _ -> val }
triangularView StrictlyLower = imap $ \row col val -> case compare row col of { GT -> val; _ -> 0 }
triangularView StrictlyUpper = imap $ \row col val -> case compare row col of { LT -> val; _ -> 0 }
triangularView UnitLower     = imap $ \row col val -> case compare row col of { GT -> val; LT -> 0; EQ -> 1 }
triangularView UnitUpper     = imap $ \row col val -> case compare row col of { LT -> val; GT -> 0; EQ -> 1 }

-- | Lower trinagle of the matrix. Shortcut for @triangularView Lower@
lowerTriangle :: I.Elem a b => Matrix a b -> Matrix a b
lowerTriangle = triangularView Lower

-- | Upper trinagle of the matrix. Shortcut for @triangularView Upper@
upperTriangle :: I.Elem a b => Matrix a b -> Matrix a b
upperTriangle = triangularView Upper

-- | Filter elements in the matrix. Filtered elements will be replaced by 0
filter :: I.Elem a b => (a -> Bool) -> Matrix a b -> Matrix a b
filter f = map (\x -> if f x then x else 0)

-- | Filter elements in the matrix. Filtered elements will be replaced by 0
ifilter :: I.Elem a b => (Int -> Int -> a -> Bool) -> Matrix a b -> Matrix a b
ifilter f = imap (\r c x -> if f r c x then x else 0)

-- | Reduce matrix using user provided function applied to each element.
fold :: I.Elem a b => (c -> a -> c) -> c -> Matrix a b -> c
fold f a (Matrix _ _ vals) = VS.foldl (\a x -> f a (I.cast x)) a vals

-- | Reduce matrix using user provided function applied to each element. This is strict version of 'fold'
fold' :: I.Elem a b => (c -> a -> c) -> c -> Matrix a b -> c
fold' f a (Matrix _ _ vals) = VS.foldl' (\a x -> f a (I.cast x)) a vals

-- | Reduce matrix using user provided function applied to each element and it's index
ifold :: I.Elem a b => (Int -> Int -> c -> a -> c) -> c -> Matrix a b -> c
ifold f a (Matrix rows _ vals) = VS.ifoldl (\a n x -> let (c,r) = divMod n rows in f r c a (I.cast x)) a vals

-- | Reduce matrix using user provided function applied to each element and it's index. This is strict version of 'ifold'
ifold' :: I.Elem a b => (Int -> Int -> c -> a -> c) -> c -> Matrix a b -> c
ifold' f a (Matrix rows _ vals) = VS.ifoldl' (\a n x -> let (c,r) = divMod n rows in f r c a (I.cast x)) a vals

-- | Reduce matrix using user provided function applied to each element.
fold1 :: I.Elem a b => (a -> a -> a) -> Matrix a b -> a
fold1 f = foldl1 f . P.map I.cast . VS.toList . _vals

-- | Reduce matrix using user provided function applied to each element. This is strict version of 'fold'
fold1' :: I.Elem a b => (a -> a -> a) -> Matrix a b -> a
fold1' f = L.foldl1' f . P.map I.cast . VS.toList . _vals

-- | Diagonal of the matrix
diagonal :: I.Elem a b => Matrix a b -> Matrix a b
diagonal = _unop (\(rows, cols) -> (min rows cols, 1)) I.diagonal

{- | Inverse of the matrix

For small fixed sizes up to 4x4, this method uses cofactors. In the general case, this method uses PartialPivLU decomposition
-}
inverse :: I.Elem a b => Matrix a b -> Matrix a b
inverse m
    | square m = _unop id I.inverse m
    | otherwise = error "Matrix.inverse: non-square matrix"

-- | Adjoint of the matrix
adjoint :: I.Elem a b => Matrix a b -> Matrix a b
adjoint = _unop swap I.adjoint

-- | Transpose of the matrix
transpose :: I.Elem a b => Matrix a b -> Matrix a b
transpose = _unop swap I.transpose

-- | Conjugate of the matrix
conjugate :: I.Elem a b => Matrix a b -> Matrix a b
conjugate = _unop id I.conjugate

-- | Nomalize the matrix by deviding it on its 'norm'
normalize :: I.Elem a b => Matrix a b -> Matrix a b
normalize (Matrix rows cols vals) = I.performIO $ do
    vals <- VS.thaw vals
    VSM.unsafeWith vals $ \p ->
        I.call $ I.normalize p (I.cast rows) (I.cast cols)
    Matrix rows cols <$> VS.unsafeFreeze vals

-- | Apply a destructive operation to a matrix. The operation will be performed in place if it is safe to do so and will modify a copy of the matrix otherwise.
modify :: I.Elem a b => (forall s. M.MMatrix a b s -> ST s ()) -> Matrix a b -> Matrix a b
modify f (Matrix rows cols vals) = Matrix rows cols (VS.modify (f . M.MMatrix rows cols) vals)

-- | Convert matrix to different type using user provided element converter
convert :: (I.Elem a b, I.Elem c d) => (a -> c) -> Matrix a b -> Matrix c d
convert f (Matrix rows cols vals) = Matrix rows cols $ VS.map (I.cast . f . I.cast) vals

-- | Yield an immutable copy of the mutable matrix
freeze :: I.Elem a b => PrimMonad m => M.MMatrix a b (PrimState m) -> m (Matrix a b)
freeze (M.MMatrix mrows mcols mvals) = VS.freeze mvals >>= return . Matrix mrows mcols

-- | Yield a mutable copy of the immutable matrix
thaw :: I.Elem a b => PrimMonad m => Matrix a b -> m (M.MMatrix a b (PrimState m))
thaw (Matrix rows cols vals) = VS.thaw vals >>= return . M.MMatrix rows cols

-- | Unsafe convert a mutable matrix to an immutable one without copying. The mutable matrix may not be used after this operation.
unsafeFreeze :: I.Elem a b => PrimMonad m => M.MMatrix a b (PrimState m) -> m (Matrix a b)
unsafeFreeze (M.MMatrix mrows mcols mvals) = VS.unsafeFreeze mvals >>= return . Matrix mrows mcols

-- | Unsafely convert an immutable matrix to a mutable one without copying. The immutable matrix may not be used after this operation.
unsafeThaw :: I.Elem a b => PrimMonad m => Matrix a b -> m (M.MMatrix a b (PrimState m))
unsafeThaw (Matrix rows cols vals) = VS.unsafeThaw vals >>= return . M.MMatrix rows cols

-- | Pass a pointer to the matrix's data to the IO action. The data may not be modified through the pointer.
unsafeWith :: I.Elem a b => Matrix a b -> (Ptr b -> CInt -> CInt -> IO c) -> IO c
unsafeWith m@(Matrix rows cols vals) f
    | not (valid m) = fail "Matrix.unsafeWith: matrix layout is invalid"
    | otherwise = VS.unsafeWith vals $ \p -> f p (I.cast rows) (I.cast cols)

{-# INLINE _prop #-}
_prop :: I.Elem a b => (Ptr b -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> a
_prop f m = I.cast $ I.performIO $ alloca $ \p -> do
    I.call $ unsafeWith m (f p)
    peek p

{-# INLINE _binop #-}
_binop :: I.Elem a b => ((Int, Int) -> (Int, Int) -> (Int, Int)) -> (Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> Matrix a b -> Matrix a b
_binop f g m1 m2 = I.performIO $ do
    m0 <- uncurry M.new $ f (dims m1) (dims m2)
    M.unsafeWith m0 $ \vals0 rows0 cols0 ->
        unsafeWith m1 $ \vals1 rows1 cols1 ->
            unsafeWith m2 $ \vals2 rows2 cols2 ->
                I.call $ g
                    vals0 rows0 cols0
                    vals1 rows1 cols1
                    vals2 rows2 cols2
    unsafeFreeze m0

{-# INLINE _unop #-}
_unop :: I.Elem a b => ((Int,Int) -> (Int,Int)) -> (Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> Matrix a b
_unop f g m1 = I.performIO $ do
    m0 <- uncurry M.new $ f (dims m1)
    M.unsafeWith m0 $ \vals0 rows0 cols0 ->
        unsafeWith m1 $ \vals1 rows1 cols1 ->
            I.call $ g
                vals0 rows0 cols0
                vals1 rows1 cols1
    unsafeFreeze m0

{-# INLINE _vals #-}
_vals :: I.Elem a b => Matrix a b -> VS.Vector b
_vals (Matrix _ _ vals) = vals

max_stars_repo_name: osidorkin/haskell-eigen

max_issues_repo_name: osidorkin/haskell-eigen

max_forks_repo_name: osidorkin/haskell-eigen


text:
{-# LANGUAGE AllowAmbiguousTypes, DataKinds, RankNTypes, TypeFamilies,
             TypeOperators #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}
module TestMnistRNN (testTrees, shortTestForCITrees) where

import Prelude

import           Control.Monad (foldM)
import qualified Data.Array.DynamicS as OT
import           Data.Array.Internal (valueOf)
import           Data.List (foldl', unfoldr)
import           Data.Proxy (Proxy (Proxy))
import qualified Data.Vector.Generic as V
import           GHC.TypeLits (KnownNat)
import           Numeric.LinearAlgebra (Matrix, Vector)
import qualified Numeric.LinearAlgebra as HM
import           System.IO (hPutStrLn, stderr)
import           System.Random
import           Test.Tasty
import           Test.Tasty.HUnit hiding (assert)
import           Text.Printf

-- until stylish-haskell accepts NoStarIsType
import qualified GHC.TypeLits

import HordeAd
import HordeAd.Core.OutdatedOptimizer
import HordeAd.Tool.MnistRnnShaped
import HordeAd.Tool.MnistTools

testTrees :: [TestTree]
testTrees = [ sinRNNTests
            , mnistRNNTestsShort
            , mnistRNNTestsLong
            ]

shortTestForCITrees :: [TestTree]
shortTestForCITrees = [ sinRNNTests
                      , mnistRNNTestsShort
                      ]

-- * A recurrent net and an autoregressive model for sine, following
-- https://blog.jle.im/entry/purely-functional-typed-models-2.html
-- and obtaining matching results

-- A version written using matrices

hiddenLayerSinRNN :: DualMonad d r m
                  => r
                  -> DualNumber d (Vector r)
                  -> DualNumberVariables d r
                  -> m (DualNumber d (Vector r), DualNumber d (Vector r))
hiddenLayerSinRNN x s variables = do
  let wX = var2 variables 0
      wS = var2 variables 1
      b = var1 variables 0
  y <- returnLet $ wX #>!! V.singleton x + wS #>! s + b
  yLogistic <- logisticAct y
  return (y, yLogistic)

outputLayerSinRNN :: DualMonad d r m
                  => DualNumber d (Vector r)
                  -> DualNumberVariables d r
                  -> m (DualNumber d r)
outputLayerSinRNN vec variables = do
  let w = var1 variables 1
      b = var0 variables 0
  returnLet $ w <.>! vec + b

fcfcrnn :: DualMonad d r m
        => r
        -> DualNumber d (Vector r)
        -> DualNumberVariables d r
        -> m (DualNumber d r, DualNumber d (Vector r))
fcfcrnn x s variables = do
  (hiddenLayer, sHiddenLayer) <- hiddenLayerSinRNN x s variables
  outputLayer <- outputLayerSinRNN hiddenLayer variables
  return (outputLayer, sHiddenLayer)

unrollLast' :: forall d r m. DualMonad d r m
            => (r
                -> DualNumber d (Vector r)
                -> DualNumberVariables d r
                -> m (DualNumber d r, DualNumber d (Vector r)))
            -> (Vector r
                -> DualNumber d (Vector r)
                -> DualNumberVariables d r
                -> m (DualNumber d r, DualNumber d (Vector r)))
unrollLast' f xs s0 variables =
  let g :: (DualNumber d r, DualNumber d (Vector r)) -> r
        -> m (DualNumber d r, DualNumber d (Vector r))
      g (_, s) x = f x s variables
  in V.foldM' g (undefined, s0) xs

zeroState :: DualMonad d r m
          => Int
          -> (a
              -> DualNumber d (Vector r)
              -> DualNumberVariables d r
              -> m (DualNumber d r2, DualNumber d (Vector r)))
          -> (a
              -> DualNumberVariables d r
              -> m (DualNumber d r2))
zeroState k f xs variables =
  fst <$> f xs (constant $ HM.konst 0 k) variables

nnSinRNN :: DualMonad d r m
         => Vector r
         -> DualNumberVariables d r
         -> m (DualNumber d r)
nnSinRNN = zeroState 30 (unrollLast' fcfcrnn)

nnSinRNNLoss :: DualMonad d r m
             => (Vector r, r)
             -> DualNumberVariables d r
             -> m (DualNumber d r)
nnSinRNNLoss (xs, target) variables = do
  result <- nnSinRNN xs variables
  lossSquared target result

series :: [Double]
series = [sin (2 * pi * t / 25) | t <- [0 ..]]

samples :: [(Vector Double, Double)]
samples  = [(V.fromList $ init c, last c) | c <- chunksOf 19 series]

sgdShow :: HasDelta r
        => (a -> DualNumberVariables 'DModeGradient r -> DualMonadGradient r (DualNumber 'DModeGradient r))
        -> [a]
        -> Domains r
        -> r
sgdShow f trainData parameters =
  let result = fst $ sgd 0.1 f trainData parameters
  in snd $ dReverse 1 (f $ head trainData) result

sgdTestCase :: String
            -> (a
                -> DualNumberVariables 'DModeGradient Double
                -> DualMonadGradient Double (DualNumber 'DModeGradient Double))
            -> (Int, [Int], [(Int, Int)], [OT.ShapeL])
            -> IO [a]
            -> Double
            -> TestTree
sgdTestCase prefix f nParameters trainDataIO expected =
  let ((nParams0, nParams1, nParams2, _), totalParams, range, parameters0) =
        initializerFixed 44 0.05 nParameters
      name = prefix ++ " "
             ++ unwords [ show nParams0, show nParams1, show nParams2
                        , show totalParams, show range ]
  in testCase name $ do
       trainData <- trainDataIO
       sgdShow f trainData parameters0
         @?= expected

sgdTestCaseAlt :: String
            -> (a
                -> DualNumberVariables 'DModeGradient Double
                -> DualMonadGradient Double (DualNumber 'DModeGradient Double))
            -> (Int, [Int], [(Int, Int)], [OT.ShapeL])
            -> IO [a]
            -> [Double]
            -> TestTree
sgdTestCaseAlt prefix f nParameters trainDataIO expected =
  let ((nParams0, nParams1, nParams2, _), totalParams, range, parameters0) =
        initializerFixed 44 0.05 nParameters
      name = prefix ++ " "
             ++ unwords [ show nParams0, show nParams1, show nParams2
                        , show totalParams, show range ]
  in testCase name $ do
       trainData <- trainDataIO
       let res = sgdShow f trainData parameters0
       assertBool ("wrong result: " ++ show res ++ " is expected to be a member of " ++ show expected) $ res `elem` expected

prime :: IsScalar 'DModeGradient r
      => (r
          -> DualNumber 'DModeGradient (Vector r)
          -> DualNumberVariables 'DModeGradient r
          -> DualMonadValue r (DualNumber 'DModeGradient r, DualNumber 'DModeGradient (Vector r)))
      -> Domains r
      -> Vector r
      -> [r]
      -> Vector r
prime f parameters =
  foldl' (\s x -> primalValue (fmap snd . f x (constant s)) parameters)

feedback :: IsScalar 'DModeGradient r
         => (r
             -> DualNumber 'DModeGradient (Vector r)
             -> DualNumberVariables 'DModeGradient r
             -> DualMonadValue r (DualNumber 'DModeGradient r, DualNumber 'DModeGradient (Vector r)))
         -> Domains r
         -> Vector r
         -> r
         -> [r]
feedback f parameters s0 x0 =
  let go (x, s) =
        let (D y _, sd') = primalValueGeneral (f x s) parameters
        in Just (x, (y, sd'))
  in unfoldr go (x0, constant s0)

feedbackTestCase :: String
                 -> (Double
                     -> DualNumber 'DModeGradient (Vector Double)
                     -> DualNumberVariables 'DModeGradient Double
                     -> DualMonadValue Double
                                        ( DualNumber 'DModeGradient Double
                                        , DualNumber 'DModeGradient (Vector Double) ))
                 -> (a
                     -> DualNumberVariables 'DModeGradient Double
                     -> DualMonadGradient Double (DualNumber 'DModeGradient Double))
                 -> (Int, [Int], [(Int, Int)], [OT.ShapeL])
                 -> [a]
                 -> [Double]
                 -> TestTree
feedbackTestCase prefix fp f nParameters trainData expected =
  let ((nParams0, nParams1, nParams2, _), totalParams, range, parameters0) =
        initializerFixed 44 0.05 nParameters
      name = prefix ++ " "
             ++ unwords [ show nParams0, show nParams1, show nParams2
                        , show totalParams, show range ]
      trained = fst $ sgd 0.1 f trainData parameters0
      primed = prime fp trained (HM.konst 0 30) (take 19 series)
      output = feedback fp trained primed (series !! 19)
  in testCase name $
       take 30 output @?= expected

-- A version written using vectors

hiddenLayerSinRNNV :: DualMonad d r m
                   => r
                   -> DualNumber d (Vector r)
                   -> DualNumberVariables d r
                   -> m (DualNumber d (Vector r), DualNumber d (Vector r))
hiddenLayerSinRNNV x s variables = do
  let wX = var1 variables 0
      b = var1 variables 31
  y <- returnLet
       $ scale (HM.konst x 30) wX + sumTrainableInputsL s 1 variables 30 + b
  yLogistic <- logisticAct y
  return (y, yLogistic)

outputLayerSinRNNV :: DualMonad d r m
                   => DualNumber d (Vector r)
                   -> DualNumberVariables d r
                   -> m (DualNumber d r)
outputLayerSinRNNV vec variables = do
  let w = var1 variables 32
      b = var0 variables 0
  returnLet $ w <.>! vec + b

fcfcrnnV :: DualMonad d r m
         => r
         -> DualNumber d (Vector r)
         -> DualNumberVariables d r
         -> m (DualNumber d r, DualNumber d (Vector r))
fcfcrnnV x s variables = do
  (hiddenLayer, sHiddenLayer) <- hiddenLayerSinRNNV x s variables
  outputLayer <- outputLayerSinRNNV hiddenLayer variables
  return (outputLayer, sHiddenLayer)

nnSinRNNLossV :: DualMonad d r m
              => (Vector r, r)
              -> DualNumberVariables d r
              -> m (DualNumber d r)
nnSinRNNLossV (xs, target) variables = do
  result <- zeroState 30 (unrollLast' fcfcrnnV) xs variables
  lossSquared target result

-- Autoregressive model with degree 2

ar2Sin :: DualMonad d r m
       => r
       -> DualNumber d (Vector r)
       -> DualNumberVariables d r
       -> m (DualNumber d r, DualNumber d (Vector r))
ar2Sin yLast s variables = do
  let c = var0 variables 0
      phi1 = var0 variables 1
      phi2 = var0 variables 2
      yLastLast = index0 s 0  -- dummy vector for compatibility
  y <- returnLet $ c + scale yLast phi1 + phi2 * yLastLast
  return (y, constant $ V.singleton yLast)

ar2SinLoss :: DualMonad d r m
           => (Vector r, r)
           -> DualNumberVariables d r
           -> m (DualNumber d r)
ar2SinLoss (xs, target) variables = do
  result <- zeroState 30 (unrollLast' ar2Sin) xs variables
  lossSquared target result

sinRNNTests :: TestTree
sinRNNTests = testGroup "Sine RNN tests"
  [ sgdTestCase "train" nnSinRNNLoss (1, [30, 30], [(30, 1), (30, 30)], [])
                (return $ take 30000 samples) 5.060827754123346e-5
  , feedbackTestCase "feedback" fcfcrnn nnSinRNNLoss
                     (1, [30, 30], [(30, 1), (30, 30)], [])
                     (take 10000 samples)
                     [-0.9980267284282716,-0.9655322144631203,-0.8919588317267176,-0.7773331580548076,-0.6212249872512189,-0.4246885094957385,-0.19280278430361192,6.316924614971235e-2,0.3255160857644734,0.5731149496491759,0.7872840563791541,0.957217059407527,1.0815006200684472,1.1654656874016613,1.2170717188563214,1.2437913143303263,1.251142657837598,1.2423738174804864,1.2186583377053681,1.1794148708577938,1.1226117988569018,1.0450711676413071,0.9428743310020188,0.8120257428038534,0.6495453130357101,0.45507653540664667,0.23281831228915612,-6.935736916677385e-3,-0.24789484923780786,-0.4705527193222155]
  , sgdTestCase "trainVV" nnSinRNNLossV (1, replicate 33 30, [], [])
                (return $ take 30000 samples) 4.6511403967229306e-5
      -- different random initial paramaters produce a worse result;
      -- matrix implementation faster, because the matrices still fit in cache
  , feedbackTestCase "feedbackVV" fcfcrnnV nnSinRNNLossV
                     (1, replicate 33 30, [], [])
                     (take 10000 samples)
                     [-0.9980267284282716,-0.9660899403337656,-0.8930568599923028,-0.7791304201898077,-0.6245654477568863,-0.4314435277698684,-0.2058673183484546,4.0423225394292085e-2,0.29029630688547203,0.5241984159992963,0.7250013011527577,0.8820730400055012,0.9922277361823716,1.057620382863504,1.08252746840241,1.070784986731554,1.0245016946328942,0.9438848015250431,0.827868146535437,0.6753691437632174,0.48708347071773117,0.26756701680655437,2.6913747557207532e-2,-0.21912614372802072,-0.45154893423928943,-0.6525638736434227,-0.8098403108946983,-0.9180866488182939,-0.9775459850131992,-0.9910399864230198]
  , sgdTestCase "trainAR" ar2SinLoss (3, [], [], [])
                (return $ take 30000 samples) 6.327978161031336e-23
  , feedbackTestCase "feedbackAR" ar2Sin ar2SinLoss
                     (3, [], [], [])
                     (take 10000 samples)
                     [-0.9980267284282716,-0.9510565162972417,-0.8443279255081759,-0.6845471059406962,-0.48175367412103653,-0.24868988719256901,-3.673766846290505e-11,0.24868988711894977,0.4817536740469978,0.6845471058659982,0.8443279254326351,0.9510565162207472,0.9980267283507953,0.9822872506502898,0.9048270523889208,0.7705132427021685,0.5877852522243431,0.3681245526237731,0.12533323351198067,-0.1253332336071494,-0.36812455271766376,-0.5877852523157643,-0.7705132427900961,-0.9048270524725681,-0.9822872507291605,-0.9980267284247174,-0.9510565162898851,-0.844327925497479,-0.6845471059273313,-0.48175367410584324]
  ]

-- * A 1 recurrent layer net with 128 neurons for MNIST, based on
-- https://medium.com/machine-learning-algorithms/mnist-using-recurrent-neural-network-2d070a5915a2
-- *Not* LSTM.
-- Doesn't train without Adam, regardless of whether mini-batch sgd
-- is used and whether a second recurrent layer. It does train with Adam,
-- but only after very carefully tweaking initialization. This is
-- extremely sensitive to initial parameters, more than to anything
-- else. Probably, gradient is vanishing if parameters are initialized
-- with a probability distribution that doesn't have the right variance. See
-- https://stats.stackexchange.com/questions/301285/what-is-vanishing-gradient.

hiddenLayerMnistRNNL :: DualMonad d r m
                     => Vector r
                     -> DualNumber d (Vector r)
                     -> DualNumberVariables d r
                     -> m (DualNumber d (Vector r), DualNumber d (Vector r))
hiddenLayerMnistRNNL x s variables = do
  let wX = var2 variables 0  -- 128x28
      wS = var2 variables 1  -- 128x128
      b = var1 variables 0  -- 128
      y = wX #>!! x + wS #>! s + b
  yTanh <- tanhAct y
  return (yTanh, yTanh)  -- tanh in both, as per https://github.com/keras-team/keras/blob/v2.8.0/keras/layers/legacy_rnn/rnn_cell_impl.py#L468

middleLayerMnistRNNL :: DualMonad d r m
                     => DualNumber d (Vector r)
                     -> DualNumber d (Vector r)
                     -> DualNumberVariables d r
                     -> m (DualNumber d (Vector r), DualNumber d (Vector r))
middleLayerMnistRNNL vec s variables = do
  let wX = var2 variables 3  -- 128x128
      wS = var2 variables 4  -- 128x128
      b = var1 variables 2  -- 128
      y = wX #>! vec + wS #>! s + b
  yTanh <- tanhAct y
  return (yTanh, yTanh)

outputLayerMnistRNNL :: DualMonad d r m
                     => DualNumber d (Vector r)
                     -> DualNumberVariables d r
                     -> m (DualNumber d (Vector r))
outputLayerMnistRNNL vec variables = do
  let w = var2 variables 2  -- 10x128
      b = var1 variables 1  -- 10
  returnLet $ w #>! vec + b  -- I assume there is no activations, as per https://www.tensorflow.org/api_docs/python/tf/compat/v1/layers/dense

fcfcrnnMnistL :: DualMonad d r m
              => Vector r
              -> DualNumber d (Vector r)
              -> DualNumberVariables d r
              -> m (DualNumber d (Vector r), DualNumber d (Vector r))
fcfcrnnMnistL = hiddenLayerMnistRNNL

fcfcrnnMnistL2 :: DualMonad d r m
               => Vector r
               -> DualNumber d (Vector r)
               -> DualNumberVariables d r
               -> m (DualNumber d (Vector r), DualNumber d (Vector r))
fcfcrnnMnistL2 x s@(D u _) variables = do
  let len = V.length u `div` 2
      s1 = slice1 0 len s
      s2 = slice1 len len s
  (vec1, s1') <- hiddenLayerMnistRNNL x s1 variables
  (vec2, s2') <- middleLayerMnistRNNL vec1 s2 variables
  s3 <- returnLet $ append1 s1' s2'
  return (vec2, s3)

unrollLastG :: forall d a b c m r. DualMonad d r m
            => (a -> b -> DualNumberVariables d r -> m (c, b))
            -> ([a] -> b -> DualNumberVariables d r -> m (c, b))
unrollLastG f xs s0 variables =
  let g :: (c, b) -> a -> m (c, b)
      g (_, s) x = f x s variables
  in foldM g (undefined, s0) xs

nnMnistRNNL :: forall d r m. DualMonad d r m
            => Int
            -> [Vector r]
            -> DualNumberVariables d r
            -> m (DualNumber d (Vector r))
nnMnistRNNL width x variables = do
  rnnLayer <- zeroState width (unrollLastG fcfcrnnMnistL) x variables
  outputLayerMnistRNNL rnnLayer variables

nnMnistRNNL2 :: DualMonad d r m
             => Int
             -> [Vector r]
             -> DualNumberVariables d r
             -> m (DualNumber d (Vector r))
nnMnistRNNL2 width x variables = do
  rnnLayer <- zeroState (2 * width) (unrollLastG fcfcrnnMnistL2) x variables
  outputLayerMnistRNNL rnnLayer variables

nnMnistRNNLossL :: forall d r m. DualMonad d r m
                => Int
                -> ([Vector r], Vector r)
                -> DualNumberVariables d r
                -> m (DualNumber d r)
nnMnistRNNLossL width (xs, target) variables = do
  result <- nnMnistRNNL width xs variables
  lossSoftMaxCrossEntropyV target result

nnMnistRNNLossL2 :: DualMonad d r m
                 => Int
                 -> ([Vector r], Vector r)
                 -> DualNumberVariables d r
                 -> m (DualNumber d r)
nnMnistRNNLossL2 width (xs, target) variables = do
  result <- nnMnistRNNL2 width xs variables
  lossSoftMaxCrossEntropyV target result

testMnistRNNL :: forall r. IsScalar 'DModeGradient r
              => Int -> [([Vector r], Vector r)] -> Domains r -> r
testMnistRNNL width inputs parameters =
  let matchesLabels :: ([Vector r], Vector r) -> Bool
      matchesLabels (glyph, label) =
        let nn = nnMnistRNNL width glyph
            value = primalValue nn parameters
        in V.maxIndex value == V.maxIndex label
  in fromIntegral (length (filter matchesLabels inputs))
     / fromIntegral (length inputs)

testMnistRNNL2 :: forall r. IsScalar 'DModeGradient r
               => Int -> [([Vector r], Vector r)] -> Domains r -> r
testMnistRNNL2 width inputs parameters =
  let matchesLabels :: ([Vector r], Vector r) -> Bool
      matchesLabels (glyph, label) =
        let nn = nnMnistRNNL2 width glyph
            value = primalValue nn parameters
        in V.maxIndex value == V.maxIndex label
  in fromIntegral (length (filter matchesLabels inputs))
     / fromIntegral (length inputs)

-- A version written using vectors

hiddenLayerMnistRNNV :: DualMonad d r m
                     => Int
                     -> Vector r
                     -> DualNumber d (Vector r)
                     -> DualNumberVariables d r
                     -> m (DualNumber d (Vector r), DualNumber d (Vector r))
hiddenLayerMnistRNNV width x s variables = do
  let b = var1 variables (width + width)  -- 128
      y = sumConstantDataL x 0 variables width
          + sumTrainableInputsL s width variables width
          + b
  yTanh <- tanhAct y
  return (yTanh, yTanh)

outputLayerMnistRNNV :: DualMonad d r m
                     => Int
                     -> DualNumber d (Vector r)
                     -> DualNumberVariables d r
                     -> m (DualNumber d (Vector r))
outputLayerMnistRNNV width vec variables = do
  let b = var1 variables (width + width + 1 + 10)  -- 10
  returnLet $ sumTrainableInputsL vec (width + width + 1) variables 10 + b

fcfcrnnMnistV :: DualMonad d r m
              => Int
              -> Vector r
              -> DualNumber d (Vector r)
              -> DualNumberVariables d r
              -> m (DualNumber d (Vector r), DualNumber d (Vector r))
fcfcrnnMnistV = hiddenLayerMnistRNNV

nnMnistRNNV :: DualMonad d r m
            => Int
            -> [Vector r]
            -> DualNumberVariables d r
            -> m (DualNumber d (Vector r))
nnMnistRNNV width x variables = do
  rnnLayer <- zeroState width (unrollLastG $ fcfcrnnMnistV width) x variables
  outputLayerMnistRNNV width rnnLayer variables

nnMnistRNNLossV :: DualMonad d r m
                => Int
                -> ([Vector r], Vector r)
                -> DualNumberVariables d r
                -> m (DualNumber d r)
nnMnistRNNLossV width (xs, target) variables = do
  result <- nnMnistRNNV width xs variables
  lossSoftMaxCrossEntropyV target result

testMnistRNNV :: forall r. IsScalar 'DModeGradient r
              => Int -> [([Vector r], Vector r)] -> Domains r -> r
testMnistRNNV width inputs parameters =
  let matchesLabels :: ([Vector r], Vector r) -> Bool
      matchesLabels (glyph, label) =
        let nn = nnMnistRNNV width glyph
            value = primalValue nn parameters
        in V.maxIndex value == V.maxIndex label
  in fromIntegral (length (filter matchesLabels inputs))
     / fromIntegral (length inputs)

lenMnistRNNL :: Int -> Int -> (Int, [Int], [(Int, Int)], [OT.ShapeL])
lenMnistRNNL width nLayers =
  ( 0
  , [width, 10] ++ replicate (nLayers - 1) width
  , [(width, 28), (width, width), (10, width)]
    ++ concat (replicate (nLayers - 1) [(width, width), (width, width)])
  , []
  )

lenMnistRNNV :: Int -> Int -> (Int, [Int], [(Int, Int)], [OT.ShapeL])
lenMnistRNNV width nLayers =
  ( 0
  , replicate width 28 ++ replicate width width ++ [width]
    ++ replicate 10 width ++ [10]
    ++ concat (replicate (nLayers - 1)
                (replicate width width ++ replicate width width ++ [width]))
  , []
  , []
  )

mnistTestCaseRNN
  :: String
  -> Int
  -> Int
  -> (Int
      -> ([Vector Double], Vector Double)
      -> DualNumberVariables 'DModeGradient Double
      -> DualMonadGradient Double (DualNumber 'DModeGradient Double))
  -> (Int -> [([Vector Double], Vector Double)] -> Domains Double -> Double)
  -> (Int -> Int -> (Int, [Int], [(Int, Int)], [OT.ShapeL]))
  -> Int
  -> Int
  -> Double
  -> TestTree
mnistTestCaseRNN prefix epochs maxBatches f ftest flen width nLayers
                 expected =
  let ((nParams0, nParams1, nParams2, _), totalParams, range, parameters0) =
        initializerFixed 44 0.2 (flen width nLayers)
      name = prefix ++ ": "
             ++ unwords [ show epochs, show maxBatches
                        , show width, show nLayers
                        , show nParams0, show nParams1, show nParams2
                        , show totalParams, show range ]
  in testCase name $ do
       hPutStrLn stderr $ printf "\n%s: Epochs to run/max batches per epoch: %d/%d"
              prefix epochs maxBatches
       let rws (input, target) =
             ( map (\k -> V.slice (k * 28) 28 input) [0 .. 27]
             , target )
       trainData <- map rws <$> loadMnistData trainGlyphsPath trainLabelsPath
       testData <- map rws <$> loadMnistData testGlyphsPath testLabelsPath
       -- There is some visual feedback, because some of these take long.
       let runBatch :: (Domains Double, StateAdam Double)
                    -> (Int, [([Vector Double], Vector Double)])
                    -> IO (Domains Double, StateAdam Double)
           runBatch (parameters@(!_, !_, !_, !_), stateAdam) (k, chunk) = do
             let res@(parameters2, _) =
                   sgdAdamBatch 150 (f width) chunk parameters stateAdam
                 !trainScore = ftest width chunk parameters2
                 !testScore = ftest width testData parameters2
                 !lenChunk = length chunk
             hPutStrLn stderr $ printf "\n%s: (Batch %d with %d points)" prefix k lenChunk
             hPutStrLn stderr $ printf "%s: Training error:   %.2f%%" prefix ((1 - trainScore) * 100)
             hPutStrLn stderr $ printf "%s: Validation error: %.2f%%" prefix ((1 - testScore ) * 100)
             return res
           runEpoch :: Int
                    -> (Domains Double, StateAdam Double)
                    -> IO (Domains Double)
           runEpoch n (params2, _) | n > epochs = return params2
           runEpoch n paramsStateAdam = do
             hPutStrLn stderr $ printf "\n%s: [Epoch %d]" prefix n
             let trainDataShuffled = shuffle (mkStdGen $ n + 5) trainData
                 chunks = take maxBatches
                          $ zip [1 ..] $ chunksOf 5000 trainDataShuffled
             !res <- foldM runBatch paramsStateAdam chunks
             runEpoch (succ n) res
       res <- runEpoch 1 (parameters0, initialStateAdam parameters0)
       let testErrorFinal = 1 - ftest width testData res
       testErrorFinal @?= expected

-- * A version written using matrices to express mini-batches of data
-- and so using matrix multiplication to run the neural net

hiddenLayerMnistRNNB :: DualMonad d r m
                     => Matrix r  -- the mini-batch of data 28x150
                     -> DualNumber d (Matrix r)  -- state for mini-batch 128x150
                     -> DualNumberVariables d r
                     -> m (DualNumber d (Matrix r), DualNumber d (Matrix r))
hiddenLayerMnistRNNB x s variables = do
  let wX = var2 variables 0  -- 128x28
      wS = var2 variables 1  -- 128x128
      b = var1 variables 0  -- 128
      batchSize = HM.cols x
      y = wX <>!! x + wS <>! s + asColumn2 b batchSize
  yTanh <- returnLet $ tanh y
  return (yTanh, yTanh)

middleLayerMnistRNNB :: DualMonad d r m
                     => DualNumber d (Matrix r)  -- 128x150
                     -> DualNumber d (Matrix r)  -- 128x150
                     -> DualNumberVariables d r
                     -> m (DualNumber d (Matrix r), DualNumber d (Matrix r))
middleLayerMnistRNNB batchOfVec@(D u _) s variables = do
  let wX = var2 variables 3  -- 128x128
      wS = var2 variables 4  -- 128x128
      b = var1 variables 2  -- 128
      batchSize = HM.cols u
      y = wX <>! batchOfVec + wS <>! s + asColumn2 b batchSize
  yTanh <- returnLet $ tanh y
  return (yTanh, yTanh)

outputLayerMnistRNNB :: DualMonad d r m
                     => DualNumber d (Matrix r)  -- 128x150
                     -> DualNumberVariables d r
                     -> m (DualNumber d (Matrix r))
outputLayerMnistRNNB batchOfVec@(D u _) variables = do
  let w = var2 variables 2  -- 10x128
      b = var1 variables 1  -- 10
      batchSize = HM.cols u
  returnLet $ w <>! batchOfVec + asColumn2 b batchSize

fcfcrnnMnistB :: DualMonad d r m
              => Matrix r
              -> DualNumber d (Matrix r)
              -> DualNumberVariables d r
              -> m (DualNumber d (Matrix r), DualNumber d (Matrix r))
fcfcrnnMnistB = hiddenLayerMnistRNNB

fcfcrnnMnistB2 :: DualMonad d r m
               => Matrix r  -- 28x150
               -> DualNumber d (Matrix r)  -- 256x150
               -> DualNumberVariables d r
               -> m (DualNumber d (Matrix r), DualNumber d (Matrix r))
fcfcrnnMnistB2 x s@(D u _) variables = do
  let len = HM.rows u `div` 2
      s1 = rowSlice2 0 len s
      s2 = rowSlice2 len len s
  (vec1, s1') <- hiddenLayerMnistRNNB x s1 variables
  (vec2, s2') <- middleLayerMnistRNNB vec1 s2 variables
  return (vec2, rowAppend2 s1' s2')

zeroStateB :: DualMonad d r m
           => (Int, Int)
           -> (a
               -> DualNumber d (Matrix r)
               -> DualNumberVariables d r
               -> m (DualNumber d r2, DualNumber d (Matrix r)))
           -> (a
               -> DualNumberVariables d r
               -> m (DualNumber d r2))
zeroStateB ij f xs variables =
  fst <$> f xs (constant $ HM.konst 0 ij) variables

nnMnistRNNB :: DualMonad d r m
            => Int
            -> [Matrix r]
            -> DualNumberVariables d r
            -> m (DualNumber d (Matrix r))
nnMnistRNNB width xs variables = do
  let batchSize = HM.cols $ head xs
  rnnLayer <- zeroStateB (width, batchSize) (unrollLastG fcfcrnnMnistB)
                         xs variables
  outputLayerMnistRNNB rnnLayer variables

nnMnistRNNB2 :: DualMonad d r m
             => Int
             -> [Matrix r]
             -> DualNumberVariables d r
             -> m (DualNumber d (Matrix r))
nnMnistRNNB2 width xs variables = do
  let batchSize = HM.cols $ head xs
  rnnLayer <- zeroStateB (2 * width, batchSize) (unrollLastG fcfcrnnMnistB2)
                         xs variables
  outputLayerMnistRNNB rnnLayer variables

nnMnistRNNLossB :: DualMonad d r m
                => Int
                -> ([Matrix r], Matrix r)
                -> DualNumberVariables d r
                -> m (DualNumber d r)
nnMnistRNNLossB width (xs, target) variables = do
  result <- nnMnistRNNB width xs variables
  vec@(D u _) <- lossSoftMaxCrossEntropyL target result
  returnLet $ scale (recip $ fromIntegral $ V.length u) $ sumElements0 vec

nnMnistRNNLossB2 :: DualMonad d r m
                 => Int
                 -> ([Matrix r], Matrix r)
                 -> DualNumberVariables d r
                 -> m (DualNumber d r)
nnMnistRNNLossB2 width (xs, target) variables = do
  result <- nnMnistRNNB2 width xs variables
  vec@(D u _) <- lossSoftMaxCrossEntropyL target result
  returnLet $ scale (recip $ fromIntegral $ V.length u) $ sumElements0 vec

mnistTestCaseRNNB
  :: String
  -> Int
  -> Int
  -> (Int
      -> ([Matrix Double], Matrix Double)
      -> DualNumberVariables 'DModeGradient Double
      -> DualMonadGradient Double (DualNumber 'DModeGradient Double))
  -> (Int -> [([Vector Double], Vector Double)] -> Domains Double -> Double)
  -> (Int -> Int -> (Int, [Int], [(Int, Int)], [OT.ShapeL]))
  -> Int
  -> Int
  -> Double
  -> TestTree
mnistTestCaseRNNB prefix epochs maxBatches f ftest flen width nLayers
                  expected =
  let ((nParams0, nParams1, nParams2, _), totalParams, range, parameters0) =
        initializerFixed 44 0.2 (flen width nLayers)
      name = prefix ++ ": "
             ++ unwords [ show epochs, show maxBatches
                        , show width, show nLayers
                        , show nParams0, show nParams1, show nParams2
                        , show totalParams, show range ]
  in testCase name $ do
       hPutStrLn stderr $ printf "\n%s: Epochs to run/max batches per epoch: %d/%d"
              prefix epochs maxBatches
       let rws (input, target) =
             ( map (\k -> V.slice (k * 28) 28 input) [0 .. 27]
             , target )
       trainData <- map rws <$> loadMnistData trainGlyphsPath trainLabelsPath
       testData <- map rws <$> loadMnistData testGlyphsPath testLabelsPath
       let packChunk :: [([Vector Double], Vector Double)]
                     -> ([Matrix Double], Matrix Double)
           packChunk chunk =
             let (inputs, targets) = unzip chunk
                 behead !acc ([] : _) = reverse acc
                 behead !acc l = behead (HM.fromColumns (map head l) : acc)
                                        (map tail l)
             in (behead [] inputs, HM.fromColumns targets)
           -- There is some visual feedback, because some of these take long.
           runBatch :: (Domains Double, StateAdam Double)
                    -> (Int, [([Vector Double], Vector Double)])
                    -> IO (Domains Double, StateAdam Double)
           runBatch (parameters@(!_, !_, !_, !_), stateAdam) (k, chunk) = do
             let res@(parameters2, _) =
                   sgdAdam (f width) (map packChunk $ chunksOf 150 chunk)
                           parameters stateAdam
                 !trainScore = ftest width chunk parameters2
                 !testScore = ftest width testData parameters2
                 !lenChunk = length chunk
             hPutStrLn stderr $ printf "\n%s: (Batch %d with %d points)" prefix k lenChunk
             hPutStrLn stderr $ printf "%s: Training error:   %.2f%%" prefix ((1 - trainScore) * 100)
             hPutStrLn stderr $ printf "%s: Validation error: %.2f%%" prefix ((1 - testScore ) * 100)
             return res
           runEpoch :: Int
                    -> (Domains Double, StateAdam Double)
                    -> IO (Domains Double)
           runEpoch n (params2, _) | n > epochs = return params2
           runEpoch n paramsStateAdam = do
             hPutStrLn stderr $ printf "\n%s: [Epoch %d]" prefix n
             let trainDataShuffled = shuffle (mkStdGen $ n + 5) trainData
                 chunks = take maxBatches
                          $ zip [1 ..] $ chunksOf 5000 trainDataShuffled
             !res <- foldM runBatch paramsStateAdam chunks
             runEpoch (succ n) res
       res <- runEpoch 1 (parameters0, initialStateAdam parameters0)
       let testErrorFinal = 1 - ftest width testData res
       testErrorFinal @?= expected

-- * A version written using shaped tensors

mnistTestCaseRNNS
  :: forall out_width batch_size d r m.
     ( KnownNat out_width, KnownNat batch_size
     , r ~ Double, d ~ 'DModeGradient, m ~ DualMonadGradient Double )
  => String
  -> Int
  -> Int
  -> (forall out_width' batch_size'.
      (DualMonad d r m, KnownNat out_width', KnownNat batch_size')
      => Proxy out_width'
      -> MnistDataBatchS batch_size' r
      -> DualNumberVariables d r
      -> m (DualNumber d r))
  -> (forall out_width' batch_size'.
      (IsScalar d r, KnownNat out_width', KnownNat batch_size')
      => Proxy r -> Proxy out_width'
      -> MnistDataBatchS batch_size' r
      -> Domains r
      -> r)
  -> (forall out_width'. KnownNat out_width'
      => Proxy out_width' -> (Int, [Int], [(Int, Int)], [OT.ShapeL]))
  -> Double
  -> TestTree
mnistTestCaseRNNS prefix epochs maxBatches trainWithLoss ftest flen expected =
  let proxy_out_width = Proxy @out_width
      batch_size = valueOf @batch_size
      ((_, _, _, nParamsX), totalParams, range, parametersInit) =
        initializerFixed 44 0.2 (flen proxy_out_width)
      name = prefix ++ ": "
             ++ unwords [ show epochs, show maxBatches
                        , show (valueOf @out_width :: Int), show batch_size
                        , show nParamsX, show totalParams, show range ]
  in testCase name $ do
    hPutStrLn stderr $ printf "\n%s: Epochs to run/max batches per epoch: %d/%d"
           prefix epochs maxBatches
    trainData <- map shapeBatch
                 <$> loadMnistData trainGlyphsPath trainLabelsPath
    testData <- map shapeBatch
                <$> loadMnistData testGlyphsPath testLabelsPath
    let testDataS = packBatch @LengthTestData testData
        -- There is some visual feedback, because some of these take long.
        runBatch :: (Domains r, StateAdam r)
                 -> (Int, [MnistDataS r])
                 -> IO (Domains r, StateAdam r)
        runBatch (parameters@(!_, !_, !_, !_), stateAdam) (k, chunk) = do
          let f = trainWithLoss proxy_out_width
              chunkS = map (packBatch @batch_size)
                       $ filter (\ch -> length ch >= batch_size)
                       $ chunksOf batch_size chunk
              res@(parameters2, _) = sgdAdam f chunkS parameters stateAdam
              !trainScore =
                ftest (Proxy @r) proxy_out_width
                      (packBatch @(10 GHC.TypeLits.* batch_size) chunk)
                      parameters2
              !testScore = ftest (Proxy @r) proxy_out_width
                                testDataS parameters2
              !lenChunk = length chunk
          hPutStrLn stderr $ printf "\n%s: (Batch %d with %d points)" prefix k lenChunk
          hPutStrLn stderr $ printf "%s: Training error:   %.2f%%" prefix ((1 - trainScore) * 100)
          hPutStrLn stderr $ printf "%s: Validation error: %.2f%%" prefix ((1 - testScore ) * 100)
          return res
        runEpoch :: Int -> (Domains r, StateAdam r) -> IO (Domains r)
        runEpoch n (params2, _) | n > epochs = return params2
        runEpoch n paramsStateAdam = do
          hPutStrLn stderr $ printf "\n%s: [Epoch %d]" prefix n
          let trainDataShuffled = shuffle (mkStdGen $ n + 5) trainData
              chunks = take maxBatches
                       $ zip [1 ..]
                       $ chunksOf (10 * batch_size) trainDataShuffled
          !res <- foldM runBatch paramsStateAdam chunks
          runEpoch (succ n) res
    res <- runEpoch 1 (parametersInit, initialStateAdam parametersInit)
    let testErrorFinal = 1 - ftest (Proxy @r) proxy_out_width testDataS res
    testErrorFinal @?= expected

mnistRNNTestsLong :: TestTree
mnistRNNTestsLong = testGroup "MNIST RNN long tests"
  [ mnistTestCaseRNN "99LL 1 epoch, all batches" 1 99
                     nnMnistRNNLossL testMnistRNNL lenMnistRNNL 128 1
                     8.209999999999995e-2
  , mnistTestCaseRNNB "99BB 1 epoch, all batches" 1 99
                      nnMnistRNNLossB testMnistRNNL lenMnistRNNL 128 1
                      8.209999999999995e-2
  , mnistTestCaseRNN "99LL2 1 epoch, all batches" 1 99
                     nnMnistRNNLossL2 testMnistRNNL2 lenMnistRNNL 128 2
                     6.259999999999999e-2
  , mnistTestCaseRNNB "99BB2 1 epoch, all batches" 1 99
                      nnMnistRNNLossB2 testMnistRNNL2 lenMnistRNNL 128 2
                      6.259999999999999e-2
  , mnistTestCaseRNN "99VV 1 epoch, all batches" 1 99
                     nnMnistRNNLossV testMnistRNNV lenMnistRNNV 128 1
                     6.740000000000002e-2
  , mnistTestCaseRNNS @128 @150 "1S 1 epoch, 1 batch" 1 1
                      rnnMnistLossFusedS rnnMnistTestS rnnMnistLenS
                      0.4375
  ]

mnistRNNTestsShort :: TestTree
mnistRNNTestsShort = testGroup "MNIST RNN short tests"
  [ let glyph = V.unfoldrExactN sizeMnistGlyph (uniformR (0, 1))
        label = V.unfoldrExactN sizeMnistLabel (uniformR (0, 1))
        rws v = map (\k -> V.slice (k * 28) 28 v) [0 .. 27]
        trainData = map ((\g -> (rws (glyph g), label g)) . mkStdGen) [1 .. 140]
    in sgdTestCaseAlt "randomLL 140"
                      (nnMnistRNNLossL 128)
                      (lenMnistRNNL 128 1)
                      (return trainData)
                      [39.26529871965807, 39.26529500592892]
  , let rws (input, target) =
          (map (\k -> V.slice (k * 28) 28 input) [0 .. 27], target)
    in sgdTestCase "firstLL 100 trainset samples only"
                   (nnMnistRNNLossL 128)
                   (lenMnistRNNL 128 1)
                   (map rws . take 100
                    <$> loadMnistData trainGlyphsPath trainLabelsPath)
                   2.7790856895965272
  , mnistTestCaseRNN "1LL 1 epoch, 1 batch" 1 1
                     nnMnistRNNLossL testMnistRNNL lenMnistRNNL 128 1
                     0.2845
  , mnistTestCaseRNNB "1BB 1 epoch, 1 batch" 1 1
                      nnMnistRNNLossB testMnistRNNL lenMnistRNNL 128 1
                      0.2845
  , let glyph = V.unfoldrExactN sizeMnistGlyph (uniformR (0, 1))
        label = V.unfoldrExactN sizeMnistLabel (uniformR (0, 1))
        rws v = map (\k -> V.slice (k * 28) 28 v) [0 .. 27]
        trainData = map ((\g -> (rws (glyph g), label g)) . mkStdGen) [1 .. 140]
    in sgdTestCaseAlt "randomLL2 140"
                      (nnMnistRNNLossL2 128)
                      (lenMnistRNNL 128 2)
                      (return trainData)
                      [30.061871495723956, 30.06187089990927]
  , let rws (input, target) =
          (map (\k -> V.slice (k * 28) 28 input) [0 .. 27], target)
    in sgdTestCase "firstLL2 99 trainset samples only"
                   (nnMnistRNNLossL2 128)
                   (lenMnistRNNL 128 2)
                   (map rws . take 99
                    <$> loadMnistData trainGlyphsPath trainLabelsPath)
                   2.772595855528805
  , mnistTestCaseRNN "1LL2 1 epoch, 1 batch" 1 1
                     nnMnistRNNLossL2 testMnistRNNL2 lenMnistRNNL 128 2
                     0.2945
  , mnistTestCaseRNNB "1BB2 1 epoch, 1 batch" 1 1
                      nnMnistRNNLossB2 testMnistRNNL2 lenMnistRNNL 128 2
                      0.2945
  , let glyph = V.unfoldrExactN sizeMnistGlyph (uniformR (0, 1))
        label = V.unfoldrExactN sizeMnistLabel (uniformR (0, 1))
        rws v = map (\k -> V.slice (k * 28) 28 v) [0 .. 27]
        trainData = map ((\g -> (rws (glyph g), label g)) . mkStdGen) [1 .. 100]
    in sgdTestCase "randomVV 100"
                   (nnMnistRNNLossV 128)
                   (lenMnistRNNV 128 1)
                   (return trainData)
                   48.93543453250378
  , let rws (input, target) =
          (map (\k -> V.slice (k * 28) 28 input) [0 .. 27], target)
    in sgdTestCase "firstVV 100 trainset samples only"
                   (nnMnistRNNLossV 128)
                   (lenMnistRNNV 128 1)
                   (map rws . take 100
                    <$> loadMnistData trainGlyphsPath trainLabelsPath)
                   2.749410768938081
  , mnistTestCaseRNN "1VV 1 epoch, 1 batch" 1 1
                     nnMnistRNNLossV testMnistRNNV lenMnistRNNV 128 1
                     0.3024
  , mnistTestCaseRNNS @120 @15 "1S 1 epoch, 1 batch" 1 1
                      rnnMnistLossFusedS rnnMnistTestS rnnMnistLenS
                      0.8418
  ]

ext: hs

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE TupleSections #-}
module Main where

import Map
import Sector
-- import RepaMap
import Map.PixelMap
-- import SectorTree
-- import SectorMap
import Map.Dimension

import Linear.V as L
import Data.Complex
import qualified Data.Vector as V
import Codec.Picture
-- import Data.Kind

-- mandelbrot :: Int -> Complex Double -> Bool
{-# INLINE mandelbrot #-}
mandelbrot :: (Ord t, Num t, RealFloat a) => t -> Complex a -> Bool
mandelbrot n x = mandelbrot' x x n where
  mandelbrot' z _ _ | (sqr . magnitude $ z) > 4 = False
  mandelbrot' _ _ i | i <= 0 = True
  mandelbrot' z c i = mandelbrot' (z*z + c) c (i - 1)
  sqr a = a * a

{-# INLINE toComplexCoords #-}
toComplexCoords :: V 2 Double -> Complex Double
toComplexCoords (V v) = v V.! 1 :+ v V.! 0

{-# INLINE mandelMap #-}
mandelMap :: Int -> PixelMap Double
mandelMap n = fmap toGrayScale
            . changeCoordinates toComplexCoords
            $ getPoint >>= return . mandelbrot n

{-# INLINE mandelMap' #-}
mandelMap' :: Int -> Sector 2 Double -> Int -> PixelMap Double
mandelMap' n = optimizeMap (mandelMap n)

viewSector :: Zone -> Double -> Sector 2 Double
viewSector (Zone x y r) ratio = V . V.fromList $ [(x - r, x + r), (y - r * ratio, y + r * ratio)]

data Zone = Zone { zoneX :: Double
                 , zoneY :: Double
                 , zoneR :: Double }

interestingZones :: [Zone]
interestingZones =
  [ Zone 0 0 2
  , Zone 0.1 (-0.7) (5e-4)
  , Zone (0.109)  (-0.7443) (0.005)
  -- , Zone (0.1102) (-0.7463) (0.005)
  , Zone (0.1127) (-0.7453) (6.5e-4)
  , Zone (1.0405) (-0.16)   (0.026)
  ]

main :: IO ()
main = do
  let resX          = 1920 * resMultiplier -- x and y dimensions
      resY          = 1080 * resMultiplier
      resMultiplier = 1 -- increment this to get 4k, 8k, etc.
      ratio         = (fromIntegral resX) / (fromIntegral resY)
      zone          = Zone 0 0 2
      depth         = 400
      sec           = viewSector zone ratio
      img           = bakePixelMap
                        sec
                        (resolution resX resY)
                        (mandelMap depth)
  putStrLn "Execution started."
  savePngImage ("/tmp/mandelbrot.png") img
  putStrLn "Execution finished"

max_stars_repo_name: zephyrys/myworld-rewrite

max_issues_repo_name: zephyrys/myworld-rewrite

max_forks_repo_name: zephyrys/myworld-rewrite


text:
{-# LANGUAGE OverloadedStrings #-}

module Releaser.Demand.Ops
  ( demandUniformIn3To15FixedDds
  , demandConst9FixedDds
  , demandConst10FixedDds
  , demandConst11FixedDds
  , demandUnif95_175
  , demandUnif78_158
  , demandExp135
  , demandExp118
  , demandExp105
  ) where

import           Statistics.Distribution.Exponential
import           Statistics.Distribution.Uniform

import           SimSim

import           Releaser.Demand.Type
import           Releaser.SettingsPeriod
import           Releaser.Util

dds :: Integer
dds = 7

dueDateSlack :: Time
dueDateSlack = fromInteger dds * periodLength

demandUniformIn3To15FixedDds :: ConfigDemand
demandUniformIn3To15FixedDds = ConfigDemand ("U(3,15) with DDS=" <> tshow dueDateSlack) dds (\sim -> generateOrdersUniform sim 3 15 dueDateSlack)

demandConst9FixedDds :: ConfigDemand
demandConst9FixedDds = ConfigDemand ("Const(9) with DDS=" <> tshow dueDateSlack) dds (\sim -> generateOrdersUniform sim 9 9 dueDateSlack)

demandConst10FixedDds :: ConfigDemand
demandConst10FixedDds = ConfigDemand ("Const(10) with DDS=" <> tshow dueDateSlack) dds (\sim -> generateOrdersUniform sim 10 10 dueDateSlack)

demandConst11FixedDds :: ConfigDemand
demandConst11FixedDds = ConfigDemand ("Const(11) with DDS=" <> tshow dueDateSlack) dds (\sim -> generateOrdersUniform sim 11 11 dueDateSlack)

demandUnif95_175 :: ConfigDemand
demandUnif95_175 =
  ConfigDemand
    ("Uniform interarrival-time w/ Unif(95,175) with DDS=" <> tshow dueDateSlack)
    dds
    (\sim ->
       let pts = productTypes sim
        in generateOrdersFixedDueDateSlack sim (uniformDistr 0.098958333 0.182291667) (uniformDistr 0.5001 (fromIntegral (length pts) + 0.4999)) dueDateSlack)

demandUnif78_158 :: ConfigDemand
demandUnif78_158 =
  ConfigDemand
    ("Uniform interarrival-time w/ Unif(78,158) with DDS=" <> tshow dueDateSlack)
    dds
    (\sim ->
       let pts = productTypes sim
        in generateOrdersFixedDueDateSlack sim (uniformDistr (78 / 960) (158 / 960)) (uniformDistr 0.5001 (fromIntegral (length pts) + 0.4999)) dueDateSlack)

demandExp118 :: ConfigDemand
demandExp118 =
  ConfigDemand
    ("Exponential interarrival-time w/ Exp(118) with DDS=" <> tshow dueDateSlack)
    dds
    (\sim ->
       let pts = productTypes sim
        in generateOrdersFixedDueDateSlack sim (exponential (960 / 118)) (uniformDistr 0.5001 (fromIntegral (length pts) + 0.4999)) dueDateSlack)

demandExp135 :: ConfigDemand
demandExp135 =
  ConfigDemand
    ("Exponential interarrival-time w/ Exp(135) with DDS=" <> tshow dueDateSlack)
    dds
    (\sim ->
       let pts = productTypes sim
        in generateOrdersFixedDueDateSlack sim (exponential (960 / 135)) (uniformDistr 0.5001 (fromIntegral (length pts) + 0.4999)) dueDateSlack)

demandExp105 :: ConfigDemand
demandExp105 =
  ConfigDemand
    ("Exponential interarrival-time w/ Exp(105) with DDS=" <> tshow dueDateSlack)
    dds
    (\sim ->
       let pts = productTypes sim
           -- test :: IO ()
           -- test = do
           --   nrs <- foldM (\xs _ -> do
           --                  x <- generateOrdersFixedDueDateSlack sim (exponential (960 / 105)) (uniformDistr 0.5001 (fromIntegral (length pts) + 0.4999)) dueDateSlack
           --                  return $ length x : xs
           --              ) [] [0..10000]
           --   putStrLn $ "res: " ++ show (fromIntegral (sum nrs) / fromIntegral (length nrs))
        in generateOrdersFixedDueDateSlack sim (exponential (960 / 105)) (uniformDistr 0.5001 (fromIntegral (length pts) + 0.4999)) dueDateSlack)

-- interArrivalTimeDistribution :: UniformDistribution
-- interArrivalTimeDistribution = uniformDistr (61.935483871/960) (274.285714286/960)

-- productTypeDistribution :: UniformDistribution
-- productTypeDistribution = uniformDistr 1 2

-- demandExponential

  -- generateOrdersFixedDueDateSlack sim interArrivalTimeDistribution productTypeDistribution dueDateSlack

max_stars_repo_name: schnecki/borl-releaser

max_issues_repo_name: schnecki/borl-releaser

max_forks_repo_name: schnecki/borl-releaser


text:
-- Gaussian arithmetic
-- http://www.codewars.com/kata/55182e34a5808b90f4001335/

module Data.Complex.Gaussian.Instances (Gaussian (..)) where

import Data.Complex.Gaussian (Gaussian (..))

instance Num Gaussian where
    (+) (Gaussian a b) (Gaussian c d) = Gaussian (a + c) (b + d)
    (*) (Gaussian a b) (Gaussian c d) = Gaussian (a*c-b*d) (a*d+b*c)
    negate (Gaussian a b) = Gaussian (negate a) (negate b)
    signum (Gaussian 0 0) = Gaussian 0 0
    signum (Gaussian a b) | a  > 0 && b >= 0 = Gaussian   1   0
                          | a <= 0 && b >  0 = Gaussian   0   1
                          | a  < 0 && b <= 0 = Gaussian (-1)  0
                          | a >= 0 && b <  0 = Gaussian   0 (-1)
    abs (Gaussian 0 0) = Gaussian 0 0
    abs (Gaussian a b) | a  > 0 && b >= 0 = Gaussian   a    b
                       | a <= 0 && b >  0 = Gaussian   b  (-a)
                       | a  < 0 && b <= 0 = Gaussian (-a) (-b)
                       | a >= 0 && b <  0 = Gaussian (-b)   a
    fromInteger n = Gaussian n 0

instance Integral Gaussian where
    toInteger (Gaussian a b) = a*a + b*b
    quotRem (Gaussian a b) (Gaussian c d) = (q, (Gaussian a b) - q * (Gaussian c d))
        where q = Gaussian ((a*c+b*d) `quot` (c*c+d*d)) ((b*c-a*d) `quot` (c*c+d*d))

max_stars_repo_name: gafiatulin/codewars

max_issues_repo_name: gafiatulin/codewars

max_forks_repo_name: gafiatulin/codewars


text:
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PackageImports        #-}

module FPNLA.Matrix.Instances.HMatrix () where

import           Control.DeepSeq               (NFData (rnf))
import           Control.Parallel.Strategies   (rseq, withStrategy)
import "hmatrix" Numeric.LinearAlgebra hiding  (Matrix, Vector)
--import           Numeric.LinearAlgebra.HMatrix
--import "hmatrix" Numeric.LinearAlgebra         (buildVector)
import           Numeric.LinearAlgebra.HMatrix hiding (Matrix, Vector)
import qualified Numeric.LinearAlgebra.HMatrix  as HM (Matrix, Vector)
import           FPNLA.Matrix                  (Matrix (generate_m, fromList_m, transpose_m, dim_m, elem_m, subMatrix_m, fromBlocks_m, toBlocks_m, map_m), MatrixVector (fromCols_vm, toCols_vm, row_vm, col_vm), Vector (generate_v, fromList_v, concat_v, elem_v, length_v, foldr_v, map_v, zipWith_v))
{-
instance (NFData e) => NFData (Vector e) where
    -- Asumo que Vector es estricto
    -- http://haskell.1045720.n5.nabble.com/NFData-instance-for-Numeric-LinearAlgebra-Matrix-td4265725.html
    rnf v = (withStrategy rseq v) `seq` ()

instance (NFData e) => NFData (Matrix e) where
    -- Matrix es estricto
    -- http://haskell.1045720.n5.nabble.com/NFData-instance-for-Numeric-LinearAlgebra-Matrix-td4265725.html
    rnf m = withStrategy rseq m `seq` ()
-}

buildVector :: Int -> (Int -> e) -> HM.Vector e
buildVector size f = undefined

buildMatrix :: Int -> Int -> (Int -> Int -> e) -> HM.Matrix e
buildMatrix rows cols f = undefined

instance (Element e) => Vector HM.Vector e where
    generate_v = buildVector
    fromList_v = fromList
    concat_v = undefined -- joinsas
    elem_v pos v = undefined -- v @> pos
    length_v = undefined -- dim
    foldr_v = undefined -- foldVector
    map_v = undefined -- mapVector
    zipWith_v = undefined -- zipVectorWith

instance (Element e) => Matrix HM.Matrix e where
    generate_m = buildMatrix
    fromList_m = (><)
    transpose_m = undefined -- trans
    dim_m m = (rows m, cols m)
    elem_m i j m = undefined -- m @@> (i, j)
    map_m = undefined -- mapMatrix
    --zipWith_m
    subMatrix_m posI posJ cantRows cantCols = subMatrix (posI, posJ) (cantRows, cantCols)
    fromBlocks_m = fromBlocks
    toBlocks_m = toBlocksEvery

instance (Element e) => MatrixVector HM.Matrix HM.Vector e where
    row_vm pos m = toRows m !! pos
    col_vm pos m = toColumns m !! pos
    fromCols_vm = fromColumns
    toCols_vm = toColumns

max_stars_repo_name: mauroblanco/fpnla-examples

max_issues_repo_name: mauroblanco/fpnla-examples

max_forks_repo_name: mauroblanco/fpnla-examples


text:
{-|
Module: MachineLearning.SoftmaxClassifier
Description: Softmax Classifier.
Copyright: (c) Alexander Ignatyev, 2017-2018.
License: BSD-3
Stability: experimental
Portability: POSIX

Softmax Classifier (Multiclass Logistic Regression).
-}

module MachineLearning.SoftmaxClassifier
(
  module MachineLearning.Model
  , module MachineLearning.Classification.MultiClass
  , SoftmaxClassifier(..)
)

where

import Prelude hiding ((<>))
import qualified Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra((<>), (<.>), (|||))
import qualified Data.Vector.Storable as V

import qualified MachineLearning as ML
import MachineLearning.Types (R, Vector, Matrix)
import MachineLearning.Utils (reduceByRows, sumByRows)
import MachineLearning.Model
import MachineLearning.Classification.MultiClass

-- | Softmax Classifier, takes number of classes.
data SoftmaxClassifier = Softmax Int

instance Classifier SoftmaxClassifier where
  cscore (Softmax _) x theta = scores - reduceByRows V.maximum scores
    where scores = x <> (LA.tr theta)

  chypothesis m x theta = V.fromList predictions
    where scores = cscore m x theta
          scores' = LA.toRows scores
          predictions = map (fromIntegral . LA.maxIndex) scores'

  ccost m lambda x y theta =
    let nSamples = fromIntegral $ LA.rows x
        scores = cscore m x theta
        sum_probs = sumByRows $ exp scores
        loss = LA.sumElements $ (log sum_probs) - remap scores y
        reg = ccostReg lambda theta
    in (loss + reg) / nSamples

  cgradient m lambda x y theta =
    let nSamples = fromIntegral $ LA.rows x
        ys = processOutput m y
        scores = cscore m x theta
        sum_probs = sumByRows $ exp scores
        probs = (exp scores) / sum_probs
        probs' = probs - ys
        dw =  (LA.tr probs') <> x
        reg = cgradientReg lambda theta
    in (dw + reg)/ nSamples

  cnumClasses (Softmax nLabels) = nLabels

remap :: Matrix -> Vector -> Matrix
remap m v = LA.remap cols rows m
  where cols = LA.asColumn $ V.fromList [0..(fromIntegral $ LA.rows m)-1]
        rows = LA.toInt $ LA.asColumn v

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}

module Main where

import Codec.Compression.GZip (decompress)
import Numeric.LinearAlgebra as V
import qualified Data.ByteString.Lazy as BS
import qualified Data.ByteString as B
import Data.Bifunctor
import Data.Functor
import Data.Word8
import qualified Data.Vector.Storable as V
import Control.Monad
import Control.Arrow
import Control.Applicative
import Pipes as P
import qualified Pipes.Prelude as PP
import System.Random
import Parser
import Data.Network
import Data.Network.FullyConnected as FC
import Algorithm.Training
import Algorithm.GradientDescent
import Data.Network.Activation
import Data.Network.Utils (makeBatchP, repeatM, squeeze, unzipP, integersP, toGaussian, averageP, accuracyBatch)

-- import System.Random
-- import Data.Binary.Get
-- import Control.Monad.Except as E
-- import qualified Data.Vector as V

main = do
  putStrLn "Extracting MNIST dataset"
  s <- decompress <$> BS.readFile "train-images-idx3-ubyte.gz" -- images data lazy bytestring
  l <- decompress <$> BS.readFile "train-labels-idx1-ubyte.gz" -- labels data lazy bytestring

  s' <- decompress <$> BS.readFile "t10k-images-idx3-ubyte.gz" -- images data lazy bytestring
  l' <- decompress <$> BS.readFile "t10k-labels-idx1-ubyte.gz" -- labels data lazy bytestring

  let mnisttraining = parseMNISTSet s l
  case mnisttraining of
    Right ((nrows, ncols), nsamples, samples) -> do
      putStrLn "MNIST training dataset succesfully parsed!"
      putStrLn $ show nsamples ++ " training samples of " ++ show nrows ++ "x" ++ show ncols ++ " labeled images."
      -- let s = take 100 samples
      -- forM_ s $ \(im, l) -> do
      --   putStrLn $ renderImage (nrows, ncols) im
      --   putStrLn $ "Target: " ++ show l
      --   putStrLn $ "Label: " ++ show (V.maxIndex l)
      let mnisttest = parseMNISTSet s' l'
      case mnisttest of
        Right ((nrows', ncols'), nsamples', samples') -> do
          putStrLn "MNIST test dataset succesfully parsed!"
          putStrLn $ show nsamples' ++ " test samples of " ++ show nrows' ++ "x" ++ show ncols' ++ " labeled images."

          let batchsize = 20
              -- Sample producer
              trainingSamplesP = each samples :: Producer MNISTSample IO ()
              (imageP, labelP) = unzipP trainingSamplesP :: (Producer Image IO (), Producer Label IO ())
              -- Batch sample producer(pairs of matrix)
              trainingBatches = PP.zip (imageP >-> (makeBatchP batchsize)) (labelP  >-> (makeBatchP batchsize))
              testSamplesP    = each samples' :: Producer MNISTSample IO ()
              testBatches     = uncurry PP.zip $ bimap (>-> makeBatchP batchsize) (>-> makeBatchP batchsize) (unzipP testSamplesP)
          g <- getStdGen
          -- let xavier = fmap (\(FC w b af) -> FC (scale (sqrt $ recip $ fromIntegral $ (uncurry (+)) (size w)) w) b af)
          let !(model :: FeedForwardFC Batch) = (makeRandomNetwork g id id 784 [(1000, logistic)] (10, softmaxMat))
              -- !model' = fmap (\(FC w b af) -> let (u,_,_) = svd w in (FC u b af)) model
              lf = mse
              lrf = const 0.05 -- \n -> exp (-(1.0 + fromIntegral n))
              tolerance = 0.0001
              test m = do
                (i,t) <- await
                let r = runNetwork m i
                    loss = evaluate lf (t, fst r)
                    a    = accuracyBatch (fst r, t)
                yield ((toDouble loss), a)
              reporter (i,t) ts = do
                putStrLn $ show ts
              interval = 10
              stopper = (uncurry (liftM2 (<|>)) . ((closeEnough tolerance) &&& (longEnough 100)))
              loop m = do -- the actual training steps
                s@(sid, (i,t)) <- await
                -- lift $ putStrLn $ "Input " ++ show (fst s)
                r <- lift $ runEffect $ (stepP m lf 0 lrf s) >-> reportP interval (reporter (i,t)) >->
                  (descendP (stopper))
                -- lift $ putStrLn ""
                yield r
                loop (tsmodel r)

              trainingChain m = PP.zip (integersP 0) trainingBatches >-> (PP.take 100) >-> (loop m)
              -- Train over all training samples, then test with testing dataset, for each iteration
              -- One iteration = one epoch
              loop' m n = do
                lift $ putStrLn $ "Epoch: " ++ show n
                (Just r) <- lift $ PP.last (trainingChain m)
                yield r
                let (lossP, accuracyP) = unzipP $ testBatches >-> test (tsmodel r)
                (loss, accuracy) <- lift $ liftM2 (,) (averageP lossP) (averageP accuracyP)
                lift $ putStrLn $ "Test loss: " ++ show (fst loss)
                lift $ putStrLn $ "Test accuracy: " ++ show (fst accuracy)
                loop' (tsmodel r) (succ n)

          runEffect $ for ((loop' model 0) >-> (PP.take 100)) $ \ts -> do
            lift $ putStrLn ""
            lift $ print ts

      putStrLn "Done!"
    Left er -> putStrLn er

max_stars_repo_name: DrPyser/NeuralNetwork

max_issues_repo_name: DrPyser/NeuralNetwork

max_forks_repo_name: DrPyser/NeuralNetwork


text:
-------------------------------------------------------------------------------
-- |
-- Module    :  Spaces.Action
-- Copyright :  (c) Sentenai 2017
-- License   :  BSD3
-- Maintainer:  sam@sentenai.com
-- Stability :  experimental
-- Portability: non-portable
--
-- typeclass for a discrete action space, as well as helper functions
-------------------------------------------------------------------------------
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Reinforce.Spaces.Action
  ( DiscreteActionSpace(..)
  -- , oneHot
  , oneHot'
  , allActions
  , randomChoice
  ) where

import Control.Monad.IO.Class
-- import Numeric.LinearAlgebra.Static (R)
-- import qualified Numeric.LinearAlgebra.Static as LA

import Control.MonadMWCRandom
import GHC.TypeLits
import Data.Proxy
import Data.Vector (Vector)
import qualified Data.Vector as V

-- | Mostly tags around an Enum, but includes information about the size of
-- an action space and is used in helper functions.
class (Bounded a, Enum a) => DiscreteActionSpace a where
  type Size a :: Nat

  toAction :: Int -> a
  toAction = toEnum

  fromAction :: a -> Int
  fromAction = fromEnum

-- | one-hot encode a bounded enumerable. Doesn't care if minBound is < or > 0
-- oneHot :: forall a . (KnownNat (Size a), DiscreteActionSpace a) => a -> R (Size a)
-- oneHot e = LA.vector . V.toList
--   $ V.unsafeUpd (replicateZeros (Proxy :: Proxy a)) [(fromEnum e, 1)]

-- | one-hot encode a bounded enumerable
oneHot' :: forall a . (DiscreteActionSpace a) => a -> Vector Double
oneHot' e = V.unsafeUpd (replicateZeros (Proxy :: Proxy a)) [(fromEnum e, 1)]

-- | helper function to initialize a one-hot vector
replicateZeros :: forall a . (Enum a, Bounded a) => Proxy a -> Vector Double
replicateZeros _ = V.fromList $ replicate (fromEnum (maxBound :: a) + 1) 0

-- | helper function to get all actions in a discrete action space
allActions :: DiscreteActionSpace a => [a]
allActions = [minBound..maxBound]

-- | make a uniform-random selection of an Action in a discrete action space
randomChoice
  :: forall m a . (MonadIO m , MonadMWCRandom m, DiscreteActionSpace a)
  => m a
randomChoice = toEnum . fst <$> sampleFrom uniformDist
  where
    uniformDist :: [Double]
    uniformDist = fmap (\a -> convert a / total) allActions
      where
        convert :: a -> Double
        convert = fromIntegral . fromEnum

        total :: Double
        total = sum (fmap convert allActions)

max_stars_repo_name: juliendehos/reinforce

max_issues_repo_name: sentenai/reinforce

max_forks_repo_name: sentenai/reinforce


text:
module DeterministicSpec where

import           Test.Hspec

import           Data.Complex
import           Data.List                      ( sortBy )
import           Data.Ord                       ( compare )
import           Data.Ratio
import           Pineappl                       ( Prob(P)
                                                , LogProb(LP)
                                                )
import           Pineappl.Deterministic         ( WrappedBDDist(BDDist)
                                                , bddist
                                                , condition
                                                , factor
                                                , runBDDist
                                                , sample
                                                , uniform
                                                )
import           Pineappl.Quantum               ( fromAmplitude )

spec :: Spec
spec = do
  describe "DIPPL tutorial milestones" $ do
    let fairCoin = uniform [0, 1]
        coinFlip = sample fairCoin
    it "should get correct result for 3 coin flips"
      $ BDDist ((\a b c -> a + b + c) <$> coinFlip <*> coinFlip <*> coinFlip)
      `shouldBe` bddist
                   [ (0, P (1 % 8))
                   , (1, P (3 % 8))
                   , (2, P (3 % 8))
                   , (3, P (1 % 8))
                   ]

    it "should get correct result for 3 coin flips (monadic)"
      $          BDDist
                   (do
                     a <- coinFlip
                     b <- coinFlip
                     c <- coinFlip
                     return (a + b + c)
                   )
      `shouldBe` bddist
                   [ (0, P (1 % 8))
                   , (1, P (3 % 8))
                   , (2, P (3 % 8))
                   , (3, P (1 % 8))
                   ]

    it "should get correct result for 'funny binomial'"
      $          BDDist
                   (do
                     a <- coinFlip
                     b <- coinFlip
                     c <- coinFlip
                     factor $ P (if a + b > 0 then 1 else 1 % 4)
                     return (a + b + c)
                   )
      `shouldBe` bddist
                   [ (0, P (1 % 32))
                   , (1, P (9 % 32))
                   , (2, P (3 % 8))
                   , (3, P (1 % 8))
                   ]

  describe "Monty Hall" $ do
    it "should get correct answer with the switch"
      $          BDDist
                   (do
                     initialChoice <- sample $ uniform [1 .. 3]
                     revealed      <- case initialChoice of
                       1 -> sample $ uniform [2, 3]
                       2 -> return 3
                       3 -> return 2
                       _ -> error "oops"
                     case (initialChoice, revealed) of
                       (1, 2) -> return 3
                       (1, 3) -> return 2
                       (2, 3) -> return 1
                       (3, 2) -> return 1
                       _      -> error "oops"
                   )
      `shouldBe` bddist [(1, P (1 % 2)), (2, P (1 % 4)), (3, P (1 % 4))]

  describe "Bayes theorem" $ do
    let bernoulli p = bddist [(True, p), (False, 1 - p)]

    it "should get correct answer for flu testing problem"
      $          BDDist
                   (do
                     hasFlu     <- sample . bernoulli $ P (1 % 10)
                     testResult <- if hasFlu
                       then sample . bernoulli $ P (7 % 10)
                       else sample . bernoulli $ P (1 % 10)
                     condition testResult
                     return hasFlu
                   )
      `shouldBe` bernoulli (P (7 % 16))

    it "should get correct answer for traffic problem"
      $ let president = bernoulli $ P (1 % 100)
            accident  = bernoulli $ P (1 % 10)
            traffic p a = bernoulli . P $ case (p, a) of
              (False, False) -> 1 % 10
              (False, True ) -> 1 % 2
              (True , False) -> 3 % 5
              (True , True ) -> 9 % 10
            joint = do
              p <- sample president
              a <- sample accident
              t <- sample $ traffic p a
              return (p, a, t)
        in  do
              BDDist (joint >>= \(p, a, t) -> condition t >> return a)
                `shouldBe` bernoulli (P (8 % 23))
              BDDist (joint >>= \(p, a, t) -> condition (t && p) >> return a)
                `shouldBe` bernoulli (P (1 % 7))

  describe "Quantum mechanics" $ do
    it "should get correct answer for double-slit experiment"
      $ let
          ampl y = LP $ exp (i * ds)
           where
            i  = 0 :+ 1
            ds = sqrt $ 1 + fromIntegral y ^ 2
          outcome = BDDist
            (do
              slit <- sample $ bddist [ (y, ampl y) | y <- [-1, 1] ]
              sample $ bddist [ (y, ampl (y - slit)) | y <- [-10 .. 10] ]
            )
          intensity = runBDDist $ fromAmplitude outcome
          sorted =
            fmap snd
              . sortBy (flip compare)
              . fmap (\(a, b) -> (b, a))
              $ intensity
        in
          {- Looks like:
              -10 | #
              -9 | ####
              -8 | ######
              -7 | ######
              -6 | ###
              -5 | #
              -4 |
              -3 | ###
              -2 | ######
              -1 | #########
               0 | ###########
               1 | #########
               2 | ######
               3 | ###
               4 |
               5 | #
               6 | ###
               7 | ######
               8 | ######
               9 | ####
              10 | # -}
          take 5 sorted `shouldBe` [0, 1, -1, 8, -8]

max_stars_repo_name: mcwitt/pineappl

max_issues_repo_name: mcwitt/pineappl

max_forks_repo_name: mcwitt/pineappl


text:
{-# LANGUAGE MagicHash, UnboxedTuples #-}

module VecUtil (unpack2, unpack3, tangentBitangent, vec3Of, vec2Of) where

import Numeric.Vector

unpack2 :: Vec2d -> (Double, Double)
unpack2 v = case unpackV2# v of
              (# x, y #) -> (x, y)

unpack3 :: Vec3d -> (Double, Double, Double)
unpack3 v = case unpackV3# v of
              (# x, y, z #) -> (x, y, z)

tangentBitangent :: Vec3d -> (Vec3d, Vec3d)
tangentBitangent normal =
  -- Cross product with an arbitrary vector to produce the tangent. The if is for the rare case that
  -- the normal is the arbitrary vector. Is there a better approach?
  let tangent = normalized $ normal × (if normal == (vec3 1 0 0) then vec3 0 1 0 else vec3 1 0 0)
      bitangent = normal × tangent -- Don't need to normalize, magnitude is mathematically 1
  in (tangent, bitangent)

vec3Of :: Double -> Vec3d
vec3Of = realToFrac

vec2Of :: Double -> Vec2d
vec2Of = realToFrac


text:
{-# LANGUAGE RecordWildCards #-}

module Arpack.Foreign.Class where

import Data.Vector.Storable.Mutable (IOVector)
import Numeric.LinearAlgebra (Matrix, Vector)

import Arpack.Options

class Arpack t where
  arpack :: Options t -> Int -> (IOVector t -> IOVector t -> IO ())
         -> IO (Vector t, Matrix t)

max_stars_repo_name: ttuegel/arpack

max_issues_repo_name: ttuegel/arpack

max_forks_repo_name: ttuegel/arpack


text:
-- | Module for operating on vectors and matrices. It includes some functions
-- for using bra-ket notation.
module DiracLists where

import Data.Complex
import Data.List (transpose)

-- This can be used to hide the standard operators.
-- import Prelude hiding ((+))
-- import qualified Prelude

-- | Function 'ket' provides basic functionality for producing vectors from the
-- canonical basis. Vectors and matrices are represented as lists and are stored
-- in column-order (as in Mathematica). Operations on vectors and matrices have
-- to be defined using new operators.
ket :: Int -> Int -> [Complex Double]
ket d i = (take i zeros) ++ [1] ++ (take (d-1-i) zeros)
  where zeros = replicate (d-1) 0

-- | The 'norm' function calculates the Euclidean norm of a vector.
norm :: [Complex Double] -> Double
norm a = sqrt $ realPart $ sum [ x^2 | x<-a ]

-- | The 'proj' function builds d-dimensional operator |i><j|.
proj :: Int -> Int -> Int -> [Complex Double]
proj d i j = outerWith (*) (ket d i) (ket d j)

-- | The 'outer' defines the form of outer product resulting in paris of elements.
outer :: [a] -> [b] -> [(a,b)]
outer a b = [ (x,y) | x<-a, y<-b ]

-- | General form of the outer product resulting in paris of elements.
outerWith :: (a -> b -> c) -> [a] -> [b] -> [c]
outerWith f a b = [ f x y | x<-a, y<-b ]

-- | Kronecker product is defined in terms of the outer product.
kron ::  Num a => [a] -> [a] -> [a]
kron a b = outerWith (*) a b

-- | Function 'overlap' for calculating the inner product of two vectors. In
-- Dirac notation this is represented as <v|w>.
overlap :: Num a => [a] -> [a] -> a
overlap vec1 vec2 = sum $ zipWith (*) vec1 vec2

-- | The scalar product is defined to be consistent with the other operators. It
-- is identical with the 'overlap' function and it suppose to resemble the
-- braket.
(<>) :: [Complex Double] -> [Complex Double] -> Complex Double
(<>) = overlap

-- | Definition of the vector addition for two complex vectors. The standard '+'
-- operator does not work, and it is not a good idead to hide it.
(+>) :: [Complex Double] -> [Complex Double] -> [Complex Double]
(+>) = zipWith (+)

-- | The multiplication of a vector by a scalar.
(.*) :: Complex Double -> [Complex Double] -> [Complex Double]
(.*) c v = [ c*x | x<-v ]

-- | Matrix-vector multiplication.
(#>) :: [[Complex Double]] -> [Complex Double] -> [Complex Double]
(#>) m v = [ overlap w v | w<-tm ]
  where tm = transpose m

max_stars_repo_name: jmiszczak/hoqus

max_issues_repo_name: jmiszczak/hoqus

max_forks_repo_name: jmiszczak/hoqus


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

module Network.Trainer
( Trainer(..)
, CostFunction
, CostFunction'
, TrainingData
, Selection
, StopCondition

, quadraticCost
, quadraticCost'
, minibatch
, online

, trainNTimes
, trainUntilErrorLessThan
, trainUntil
) where

import Network.Network
import Network.Neuron
import Network.Layer

import System.Random
import System.Random.Shuffle (shuffle')
import Data.List.Split (chunksOf)
import Numeric.LinearAlgebra

-- | Trainer is a typeclass for all trainer types - a trainer will take in
--   an instance of itself, a network, a list of training data, and return a
--   new network trained on the data.
class (Network n) => Trainer a n where
  fit :: Selection -> a -> n -> [TrainingData] -> n
  evaluate :: a -> n -> TrainingData -> Double

-- | A CostFunction is used for evaluating a network's performance on a given
--   input
type CostFunction = Vector Double -> Vector Double -> Double

-- | A CostFunction' (derivative) is used in backPropagation
type CostFunction' = Vector Double -> Vector Double -> Vector Double

-- | A tuple of (input, expected output)
type TrainingData = (Vector Double, Vector Double)

-- | A selection function for performing gradient descent
type Selection = [TrainingData] -> [[TrainingData]]

-- | A predicate (given a network, trainer, a list of training
--   data, and the number of [fit]s performed) that
--   tells the trainer to stop training
type StopCondition t n = n -> t -> [TrainingData] -> Int -> Bool

-- | The quadratic cost function (1/2) * sum (y - a) ^ 2
quadraticCost :: Vector Double -> Vector Double -> Double
quadraticCost y a = sumElements $ 0.5 * (a - y) ** 2

-- | The derivative of the quadratic cost function sum (y - a)
quadraticCost' :: Vector Double -> Vector Double -> Vector Double
quadraticCost' y a = a - y

-- | The minibatch function becomes a Selection when partially applied
--   with the minibatch size
minibatch :: Int -> [TrainingData] -> [[TrainingData]]
minibatch size = chunksOf size

-- | If we want to train the network online
online :: [TrainingData] -> [[TrainingData]]
online = minibatch 1

-- | This function returns true if the error of the network is less than
--   a given error value, given a network, a trainer, a list of
--   training data, and a counter (should start with 0)
--   Note: Is there a way to have a counter with a recursive function
--         without providing 0?
networkErrorLessThan :: (Trainer t n) => Double -> n -> t -> [TrainingData] -> Int -> Bool
networkErrorLessThan err network trainer dat _ = meanError < err
  where meanError = (sum errors) / fromIntegral (length errors)
        errors = map (evaluate trainer network) dat

-- | Given a network, a trainer, a list of training data,
--   and N, this function trains the network with the list of
--   training data N times
trainNTimes :: (Trainer t n, RandomGen g) => g -> n -> t -> Selection -> [TrainingData] -> Int -> n
trainNTimes g network trainer s dat n =
  trainUntil g network trainer s dat completion 0
  where completion _ _ _ n' = (n == n')

-- | Given a network, a trainer, a list of training data,
--   and an error value, this function trains the network with the list of
--   training data until the error of the network (calculated
--   by averaging the errors of each training data) is less than
--   the given error value
trainUntilErrorLessThan :: (Trainer t n, RandomGen g) => g -> n -> t -> Selection -> [TrainingData] -> Double -> n
trainUntilErrorLessThan g network trainer s dat err =
  trainUntil g network trainer s dat (networkErrorLessThan err) 0

-- | This function trains a network until a given TrainCompletionPredicate
--   is satisfied.
trainUntil :: (Trainer t n, RandomGen g) => g -> n -> t -> Selection -> [TrainingData] -> StopCondition t n -> Int -> n
trainUntil g network trainer s dat completion n =
  if completion network trainer dat n
    then network
    else trainUntil g' network' trainer s (shuffle' dat (length dat) g'') completion (n+1)
    where network' = fit s trainer network dat
          (g', g'') = split g

max_stars_repo_name: AkatsukiSirius/LambdaNet

max_issues_repo_name: world-admin/LambdaNet

max_forks_repo_name: world-admin/LambdaNet