text:
{-# LANGUAGE ScopedTypeVariables, UndecidableInstances #-}
{-|
  Array access, including cell arrays and structures.

  Functions here are primarily thin wrappers to the underlying Matlab functions, and the same memory-management semantics hold.
  In particular, created arrays must be freed, 'copyMXArray' and 'freeMXArray' are deep operations, other set operations do not make copies.
-}

module Foreign.Matlab.Array (
    -- * Array manipulation
    anyMXArray,
    MNullArray, castMNull,
    mxArrayClass,
    mxArrayIsComplex,
    mxArraySize,
    mxArraySetSize,
    mxArrayLength,
    freeMXArray,
    copyMXArray,
    mIndexOffset,

    -- * Array element access
    MXArrayComponent (mxArrayGetOffset, mxArraySetOffset
      , mxArrayGetOffsetList, mxArraySetOffsetList
      , mxScalarGet, isMXScalar
      , createMXArray, createMXScalar
      , createColVector, createRowVector),
    castMXArray,
    -- | array element access
    mxArrayGet, mxArraySet,
    -- | array list access
    mxArrayGetList, mxArraySetList,
    mxArrayGetAll, mxArraySetAll,
    mxArrayGetOffsetSafe, mxArrayGetFirst, mxArrayGetLast,
    -- mxArrayGetSafe, --TODO--
    fromListIO, cellFromListsIO,
    isMNull,

    -- * Struct access
    -- |Structs in Matlab are always arrays, and so can be accessed using most array accessors.
    -- |However, the modification functions such as 'mxArraySet' are not implemented because they could disrupt field data in the entire struct array, and so some specialized functions are necessary.
    MStructArray,
    createStruct,
    mStructFields,
    mStructGet, mStructSet,
    mStructSetFields,
    mStructAddField, mStructRemoveField,
    mxCellGetAllOfType, mxCellGetArraysOfType,

    -- ** Object access
    -- |Some structs are also validated (blessed) user objects.
    mObjectGetClass, mObjectSetClass
  ) where

import           Control.Monad
import           Data.Foldable (toList)
import           Foreign
import           Foreign.C.String
import           Foreign.C.Types
import           Data.Complex
import qualified Data.Map.Strict as DM
import           Data.Maybe (catMaybes)
import           Foreign.Matlab.Util
import           Foreign.Matlab.Internal
import           Foreign.Matlab.Types

#include <matrix.h>

-- |(Un)cast an array to a generic type.
anyMXArray :: MXArray a -> MAnyArray
anyMXArray a = unsafeCastMXArray a

-- |A NULL (empty) array
type MNullArray = MXArray MNull
-- |Safely cast a generic array to a NULL array, or return Nothing if the array is not NULL
castMNull :: MAnyArray -> MIO (Maybe MNullArray)
castMNull a
  | isMNull a = pure $ Just (unsafeCastMXArray a)
  | otherwise = pure Nothing

foreign import ccall unsafe mxGetClassID :: MXArrayPtr -> IO MXClassID
-- |Return the representation of the type of the elements of an array
mxArrayClass :: MXArray a -> IO MXClass
mxArrayClass a
  | isMNull a = pure $ MXClassNull
  | otherwise = withMXArray a mxGetClassID >.= mx2hs

ndims :: MWSize -> Ptr MWSize -> IO MSize
ndims n s = map ii =.< peekArray (ii n) s
--nsubs = ndims

withNSubs :: With MSubs (MWSize, Ptr MWSize) (IO a)
withNSubs l f = withArrayLen (map ii l) (\l a -> f (ii l, a))
withNDims :: With MSize (MWSize, Ptr MWSize) (IO a)
withNDims = withNSubs . realMSize

foreign import ccall unsafe mxGetNumberOfDimensions :: MXArrayPtr -> IO MWSize
foreign import ccall unsafe mxGetDimensions :: MXArrayPtr -> IO (Ptr MWSize)
-- |Get the size (dimensions) of an array
mxArraySize :: MXArray a -> MIO MSize
mxArraySize a = withMXArray a $ \a -> do
  n <- mxGetNumberOfDimensions a
  s <- mxGetDimensions a
  ndims n s

foreign import ccall unsafe mxSetDimensions :: MXArrayPtr -> Ptr MWSize -> MWSize -> IO CInt
-- |Set dimension array and number of dimensions
mxArraySetSize :: MXArray a -> MSize -> IO ()
mxArraySetSize a s = do
  r <- withMXArray a (\a -> withNDims s (\(nd,d) -> mxSetDimensions a d nd))
  when (r /= 0) $ fail "mxArraySetSize"

foreign import ccall unsafe mxGetNumberOfElements :: MXArrayPtr -> IO CSize
-- |Like `numel` in MATLAB.
mxArrayLength :: MXArray a -> MIO Int
mxArrayLength a = ii =.< withMXArray a mxGetNumberOfElements

foreign import ccall unsafe mxCalcSingleSubscript :: MXArrayPtr -> MWSize -> Ptr MWIndex -> IO MWIndex
-- |Convert an array subscript into an offset
mIndexOffset :: MXArray a -> MIndex -> MIO Int
mIndexOffset _ (MSubs []) = pure 0
mIndexOffset _ (MSubs [i]) = pure i
mIndexOffset a (MSubs i) = ii =.< withMXArray a (withNSubs i . uncurry . mxCalcSingleSubscript)

foreign import ccall unsafe mxDuplicateArray :: MXArrayPtr -> IO MXArrayPtr
-- |Make a deep copy of an array
copyMXArray :: MXArray a -> MIO (MXArray a)
copyMXArray a = withMXArray a mxDuplicateArray >>= mkMXArray

foreign import ccall unsafe mxDestroyArray :: MXArrayPtr -> IO ()
-- |Destroy an array and all of its contents.
freeMXArray :: MXArray a -> MIO ()
freeMXArray a = do
  withMXArray a mxDestroyArray
  mxArraySetSize a [0, 0]

-- | Create and populate an MXArray in one go. Named without 'mx' due to possible
-- | conformity to a typeclass function.
fromListIO :: (Foldable t, MXArrayComponent a) => t a -> MIO (MXArray a)
fromListIO xs = do
  arr <- createMXArray [length xs]
  mxArraySetAll arr xsList
  pure arr
  where
    xsList = toList xs

-- | Like fromListIO but wraps elements in a cell. Most useful for converting a list of strings
-- | to a MATLAB cell array of strings. Named in conjunction with `fromListIO`, which is used
-- | as part of the implementation.
cellFromListsIO :: (Traversable s, Foldable t, MXArrayComponent a) => s (t a) -> MIO (MXArray MCell)
cellFromListsIO xss = do
  listOfStructArrays <- sequence $ fromListIO <$> xss
  arr <- createMXArray [length xss]
  mxArraySetAll arr (toList $ (MCell . anyMXArray) <$> listOfStructArrays)
  pure arr

-- |The class of standardly typeable array elements
class MXArrayComponent a where
  -- |Determine whether the given array is of the correct type
  isMXArray :: MXArray a -> MIO Bool
  -- |Create an array and initialize all its data elements to some default value, usually 0 or []
  createMXArray :: MSize -> MIO (MXArray a)

  -- |Determine if an array is singleton. Equivalent to
  --
  -- >  liftM (all (1 ==)) . mxArraySize
  isMXScalar :: MXArray a -> MIO Bool

  mxArrayGetOffset :: MXArray a -> Int -> MIO a
  mxArraySetOffset :: MXArray a -> Int -> a -> MIO ()

  mxArrayGetOffsetList :: MXArray a -> Int -> Int -> MIO [a]
  mxArraySetOffsetList :: MXArray a -> Int -> [a] -> MIO ()

  -- |Get the value of the first data element in an array or, more specifically, the value that the array will be interpreted as in scalar context
  mxScalarGet :: MXArray a -> MIO a

  -- |Create a singleton (scalar) array having the specified value
  createMXScalar :: a -> MIO (MXArray a)
  -- |Create a column vector from the given list.
  createColVector :: [a] -> MIO (MXArray a)
  -- |Create a row vector from the given list.
  createRowVector :: [a] -> MIO (MXArray a)

  isMXArray _ = pure False
  isMXScalar a = liftM2 (&&) (isMXArray a) (all (1 ==) =.< mxArraySize a)

  mxArrayGetOffsetList a o n = mapM (mxArrayGetOffset a) [o..o+n-1]
  mxArraySetOffsetList a o = zipWithM_ (mxArraySetOffset a . (o+)) [0..]

  mxScalarGet a = mxArrayGetOffset a 0

  createMXScalar x = do
    a <- createMXArray [1]
    mxArraySetOffset a 0 x
    pure a
  createRowVector l = do
    a <- createMXArray [1,length l]
    mxArraySetOffsetList a 0 l
    pure a
  createColVector l = do
    a <- createMXArray [length l]
    mxArraySetOffsetList a 0 l
    pure a

-- |Get the value of the specified array element.  Does not check bounds.
mxArrayGet :: MXArrayComponent a => MXArray a -> MIndex -> MIO a
mxArrayGet a i = mIndexOffset a i >>= mxArrayGetOffset a
-- |Set an element in an array to the specified value.  Does not check bounds.
mxArraySet :: MXArrayComponent a => MXArray a -> MIndex -> a -> MIO ()
mxArraySet a i v = do
  o <- mIndexOffset a i
  mxArraySetOffset a o v

-- |@'mxArrayGetList' a i n@ gets the sequential list of @n@ items from array @a@ starting at index @i@.  Does not check bounds.
mxArrayGetList :: MXArrayComponent a => MXArray a -> MIndex -> Int -> MIO [a]
mxArrayGetList a i n = do
  o <- mIndexOffset a i
  n <- if n == -1 then subtract o =.< mxArrayLength a else pure n
  mxArrayGetOffsetList a o n
-- |@'mxArraySetList' a i l@ sets the sequential items in array @a@ starting at index @i@ to @l@.  Does not check bounds.
mxArraySetList :: MXArrayComponent a => MXArray a -> MIndex -> [a] -> MIO ()
mxArraySetList a i l = do
  o <- mIndexOffset a i
  mxArraySetOffsetList a o l

-- |Get a flat list of all elements in the array.
mxArrayGetAll :: MXArrayComponent a => MXArray a -> IO [a]
mxArrayGetAll a = mxArrayGetList a mStart (-1)

-- |Set a flat list of all elements in the array.
mxArraySetAll :: MXArrayComponent a => MXArray a -> [a] -> IO ()
mxArraySetAll a = mxArraySetList a mStart

mxArrayGetFirst :: MXArrayComponent a => MXArray a -> MIO (Either String a)
mxArrayGetFirst arr = mxArrayGetOffsetSafe arr 0

mxArrayGetLast :: MXArrayComponent a => MXArray a -> MIO (Either String a)
mxArrayGetLast arr = do
  arrLen <- mxArrayLength arr
  mxArrayGetOffsetSafe arr (arrLen - 1)

-- |Like mxArrayGetOffset but safe.
mxArrayGetOffsetSafe :: forall a. MXArrayComponent a => MXArray a -> Int -> MIO (Either String a)
mxArrayGetOffsetSafe arr ix
  | isMNull arr = pure $ Left "Couldn't get element of null array"
  | otherwise = do
    arrLen <- mxArrayLength arr
    safeGetElem arrLen ix
  where
    safeGetElem :: Int -> Int -> MIO (Either String a)
    safeGetElem aLen aIx
      | aIx < aLen = Right <$> mxArrayGetOffset arr aIx
      | otherwise = pure $ Left $ "Couldn't get element at index "
        <> (show aIx) <> " of " <> (show aLen) <> "-length array"

-- |Safely cast a generic array to a type, or return Nothing if the array does not have the proper type
castMXArray :: forall a. MXArrayComponent a => MAnyArray -> MIO (Maybe (MXArray a))
castMXArray a
  | isMNull a = pure Nothing
  | otherwise = do
      y <- isMXArray b
      pure $ if y then Just b else Nothing
      where
        b :: MXArray a
        b = unsafeCastMXArray a

-- | Extract all arrays of a given type from a Cell Array.
mxCellGetArraysOfType :: MXArrayComponent a => MXArray MCell -> MIO ([MXArray a])
mxCellGetArraysOfType ca = do
  cellVals <- (fmap . fmap) mCell (mxArrayGetAll ca)
  mxaMays :: [Maybe (MXArray a)] <- sequence $ castMXArray <$> cellVals
  pure $ catMaybes mxaMays

-- | A convenience function to extract all arrays of a given type from a Cell Array;
-- | may have larger dimensions than the original Cell Array due to flattening.
mxCellGetAllOfType :: MXArrayComponent a => MXArray MCell -> MIO [a]
mxCellGetAllOfType ca = do
  as <- mxCellGetArraysOfType ca
  join <$> (sequence $ mxArrayGetAll <$> as)

foreign import ccall unsafe mxGetData :: MXArrayPtr -> IO (Ptr a)

class (MXArrayComponent a, MType mx a, Storable mx) => MXArrayData mx a where
  withArrayData :: MXArray a -> (Ptr mx -> IO b) -> IO b
  withArrayDataOff :: MXArray a -> Int -> (Ptr mx -> IO b) -> IO b
  arrayDataGet :: MXArray a -> Int -> IO a
  arrayDataSet :: MXArray a -> Int -> a -> IO ()

  arrayDataGetList :: MXArray a -> Int -> Int -> IO [a]
  arrayDataSetList :: MXArray a -> Int -> [a] -> IO ()

  withArrayData a f = withMXArray a (mxGetData >=> f)
  withArrayDataOff a o f = withArrayData a (\p -> f (advancePtr p o))
  arrayDataGet a o   = withArrayDataOff a o (mx2hs .=< peek)
  arrayDataSet a o v = withArrayDataOff a o (\p -> poke p (hs2mx v))
  arrayDataGetList a o n = withArrayDataOff a o (map mx2hs .=< peekArray n)
  arrayDataSetList a o l = withArrayDataOff a o (\p -> pokeArray p (map hs2mx l))
#let arrayDataComponent = "\
mxArrayGetOffset = arrayDataGet ;\
mxArraySetOffset = arrayDataSet ;\
mxArrayGetOffsetList = arrayDataGetList ;\
mxArraySetOffsetList = arrayDataSetList\
"
--"

foreign import ccall unsafe mxIsLogical :: MXArrayPtr -> IO CBool
foreign import ccall unsafe mxCreateLogicalArray :: MWSize -> Ptr MWSize -> IO MXArrayPtr
foreign import ccall unsafe mxGetLogicals :: MXArrayPtr -> IO (Ptr MXLogical)
foreign import ccall unsafe mxCreateLogicalScalar :: CBool -> IO MXArrayPtr
foreign import ccall unsafe mxIsLogicalScalar :: MXArrayPtr -> IO CBool
foreign import ccall unsafe mxIsLogicalScalarTrue :: MXArrayPtr -> IO CBool
instance MXArrayComponent MLogical where
  isMXArray a = boolC =.< withMXArray a mxIsLogical
  createMXArray s = withNDims s (uncurry mxCreateLogicalArray) >>= mkMXArray
  createMXScalar = mxCreateLogicalScalar . cBool >=> mkMXArray
  isMXScalar a = boolC =.< withMXArray a mxIsLogicalScalar
  mxScalarGet a = boolC =.< withMXArray a mxIsLogicalScalarTrue
  #arrayDataComponent
instance MXArrayData MXLogical MLogical where
  withArrayData a f = withMXArray a (mxGetLogicals >=> f)

foreign import ccall unsafe mxIsChar :: MXArrayPtr -> IO CBool
foreign import ccall unsafe mxCreateCharArray :: MWSize -> Ptr MWSize -> IO MXArrayPtr
foreign import ccall unsafe mxGetChars :: MXArrayPtr -> IO (Ptr MXChar)
foreign import ccall unsafe mxCreateStringFromNChars :: CString -> MWSize -> IO MXArrayPtr
instance MXArrayComponent MChar where
  isMXArray a = boolC =.< withMXArray a mxIsChar
  createMXArray s = withNDims s (uncurry mxCreateCharArray) >>= mkMXArray
  createRowVector s =
    mkMXArray =<< withCStringLen s (\(s,n) -> mxCreateStringFromNChars s (ii n))
  #arrayDataComponent
instance MXArrayData MXChar MChar where
  withArrayData a f = withMXArray a (mxGetChars >=> f)

foreign import ccall unsafe mxCreateNumericArray :: MWSize -> Ptr MWSize -> MXClassID -> (#type mxComplexity) -> IO MXArrayPtr
createNumericArray :: MXClass -> Bool -> MWSize -> Ptr MWSize -> IO MXArrayPtr
createNumericArray t c n s = mxCreateNumericArray n s (hs2mx t) (if c then (#const mxCOMPLEX) else (#const mxREAL))

#let numarray t = "\
foreign import ccall unsafe mxIs%s :: MXArrayPtr -> IO CBool\n\
instance MXArrayComponent M%s where\n\
  isMXArray a = boolC =.< withMXArray a mxIs%s\n\
  createMXArray s = withNDims s (uncurry $ createNumericArray (mxClassOf (undefined :: M%s)) False) >>= mkMXArray\n\
  \
mxArrayGetOffset = arrayDataGet ;\
mxArraySetOffset = arrayDataSet ;\
mxArrayGetOffsetList = arrayDataGetList ;\
mxArraySetOffsetList = arrayDataSetList\
  \n\
instance MXArrayData MX%s M%s\
", #t, #t, #t, #t, #t, #t

foreign import ccall unsafe mxIsDouble :: MXArrayPtr -> IO CBool
foreign import ccall unsafe mxCreateDoubleScalar :: MXDouble -> IO MXArrayPtr
foreign import ccall unsafe mxGetScalar :: MXArrayPtr -> IO MXDouble
instance MXArrayComponent MDouble where
  isMXArray a = boolC =.< withMXArray a mxIsDouble
  createMXScalar = mxCreateDoubleScalar . hs2mx >=> mkMXArray
  mxScalarGet a = withMXArray a mxGetScalar
  createMXArray s = withNDims s (uncurry $ createNumericArray (mxClassOf (undefined :: Double)) False) >>= mkMXArray
  #arrayDataComponent
instance MXArrayData MXDouble MDouble

#numarray Single
#numarray Int8
#numarray Int16
#numarray Int32
#numarray Int64
#numarray Uint8
#numarray Uint16
#numarray Uint32
#numarray Uint64

foreign import ccall unsafe mxIsCell :: MXArrayPtr -> IO CBool
foreign import ccall unsafe mxCreateCellArray :: MWSize -> Ptr MWSize -> IO MXArrayPtr
foreign import ccall unsafe mxGetCell :: MXArrayPtr -> MWIndex -> IO MXArrayPtr
foreign import ccall unsafe mxSetCell :: MXArrayPtr -> MWIndex -> MXArrayPtr -> IO ()
instance MXArrayComponent MCell where
  isMXArray a = boolC =.< withMXArray a mxIsCell
  createMXArray s = withNDims s (uncurry mxCreateCellArray) >>= mkMXArray
  -- Get a the specified cell element (not a copy).
  mxArrayGetOffset a o = withMXArray a (\a -> mxGetCell a (ii o) >>= mkMXArray >.= MCell)
  -- Set an element in a cell array to the specified value. The cell takes ownership of the array: and no copy is made. Any existing value should be freed first.
  mxArraySetOffset a o (MCell v) = withMXArray a (\a -> withMXArray v (mxSetCell a (ii o)))

-- |A (array of) structs
type MStructArray = MXArray MStruct

foreign import ccall unsafe mxIsStruct :: MXArrayPtr -> IO CBool
foreign import ccall unsafe mxIsObject :: MXArrayPtr -> IO CBool

foreign import ccall unsafe mxCreateStructArray :: MWSize -> Ptr MWSize -> CInt -> Ptr CString -> IO MXArrayPtr
-- |Create an N-Dimensional structure array having the specified fields; initialize all values to 'MNullArray'
createStruct :: MSize -> [String] -> MIO MStructArray
createStruct s f =
  withNDims s (\(nd,d) ->
    mapWithArrayLen withCString f (\(f,nf) ->
      mxCreateStructArray nd d (ii nf) f))
  >>= mkMXArray

foreign import ccall unsafe mxGetNumberOfFields :: MXArrayPtr -> IO CInt
foreign import ccall unsafe mxGetFieldNameByNumber :: MXArrayPtr -> CInt -> IO CString
foreign import ccall unsafe mxGetFieldNumber :: MXArrayPtr -> CString -> IO CInt
-- |Get the names of the fields
mStructFields :: MStructArray -> MIO [String]
mStructFields a = withMXArray a $ \a -> do
  n <- mxGetNumberOfFields a
  forM [0..pred n] (mxGetFieldNameByNumber a >=> peekCString)

foreign import ccall unsafe mxGetField :: MXArrayPtr -> MWIndex -> CString -> IO MXArrayPtr
foreign import ccall unsafe mxSetField :: MXArrayPtr -> MWIndex -> CString -> MXArrayPtr -> IO ()
foreign import ccall unsafe mxGetFieldByNumber :: MXArrayPtr -> MWIndex -> CInt -> IO MXArrayPtr
foreign import ccall unsafe mxSetFieldByNumber :: MXArrayPtr -> MWIndex -> CInt -> MXArrayPtr -> IO ()

-- |Return the contents of the named field for the given element.
-- |Returns 'MNullArray' on no such field or if the field itself is NULL
mStructGet :: MStructArray -> MIndex -> String -> MIO MAnyArray
-- |Sets the contents of the named field for the given element. The input is stored in the array -- no copy is made.
mStructSet :: MStructArray -> MIndex -> String -> MXArray a -> MIO ()
mStructGet a i f = do
  o <- mIndexOffset a i
  withMXArray a (\a -> withCString f (mxGetField a (ii o) >=> mkMXArray))
mStructSet a i f v = do
  o <- mIndexOffset a i
  withMXArray a (\a -> withCString f (withMXArray v . mxSetField a (ii o)))

foreign import ccall unsafe mxAddField :: MXArrayPtr -> CString -> IO CInt
foreign import ccall unsafe mxRemoveField :: MXArrayPtr -> CInt -> IO ()
-- |Add a field to a structure array.
mStructAddField :: MStructArray -> String -> MIO ()
-- |Remove a field from a structure array. Does nothing if no such field exists.
mStructRemoveField :: MStructArray -> String -> MIO ()
mStructAddField a f = do
  i <- withMXArray a (withCString f . mxAddField)
  when (i < 0) $ fail "mxAddField"
mStructRemoveField a f = withMXArray a $ \a -> do
  i <- withCString f (mxGetFieldNumber a)
  if i < 0
    then fail "mxRemoveField"
    else mxRemoveField a i

structGetOffsetFields :: MStructArray -> [String] -> Int -> IO MStruct
structGetOffsetFields a f o =
  MStruct =.< withMXArray a (\a -> DM.fromList <$>
    (zipWithM (\f -> ((,) f) .=< (mxGetFieldByNumber a (ii o) >=> mkMXArray)) f [0..]))

-- |Set the fields of a struct index to the given value list.  The list corresponds to the field list and must match in size.
mStructSetFields :: MStructArray -> MIndex -> [MXArray a] -> MIO ()
mStructSetFields a i v = do
  o <- mIndexOffset a i
  withMXArray a (\a -> zipWithM_ (\v -> withMXArray v . mxSetFieldByNumber a (ii o)) v [0..])

instance MXArrayComponent MStruct where
  isMXArray a = liftM2 (||) (boolC =.< withMXArray a mxIsStruct) (boolC =.< withMXArray a mxIsObject)
  createMXArray s = createStruct s []
  mxArrayGetOffset a o = do
    f <- mStructFields a
    structGetOffsetFields a f o
  mxArraySetOffset = error "mxArraySet undefined for MStruct: use mStructSet"
  mxArrayGetOffsetList a o n = do
    f <- mStructFields a
    mapM (structGetOffsetFields a f) [o..o+n-1]
  createMXScalar (MStruct fv) = do
    a <- createStruct [1] f
    withMXArray a $ \a -> zipWithM_ (\i v -> withMXArray v (mxSetFieldByNumber a 0 i)) [0..] v
    pure a
    where
      (f,v) = unzip $ DM.toList fv

foreign import ccall unsafe mxGetClassName :: MXArrayPtr -> IO CString
-- |Determine if a struct array is a user defined object, and return its class name, if any.
mObjectGetClass :: MStructArray -> IO (Maybe String)
mObjectGetClass a = do
  b <- boolC =.< withMXArray a mxIsObject
  if b
    then Just =.< withMXArray a (mxGetClassName >=> peekCString)
    else pure Nothing

foreign import ccall unsafe mxSetClassName :: MXArrayPtr -> CString -> IO CInt
-- |Set classname of an unvalidated object array.  It is illegal to call this function on a previously validated object array.
mObjectSetClass :: MStructArray -> String -> IO ()
mObjectSetClass a c = do
  r <- withMXArray a (withCString c . mxSetClassName)
  when (r /= 0) $ fail "mObjectSetClass"

castReal :: MXArray (Complex a) -> MXArray a
castReal = unsafeCastMXArray

foreign import ccall unsafe mxGetImagData :: MXArrayPtr -> IO (Ptr a)
withRealDataOff :: MXArrayData mx a => MXArray (Complex a) -> Int -> (Ptr mx -> IO b) -> IO b
withRealDataOff = withArrayDataOff . castReal
withImagDataOff :: MXArrayData mx a => MXArray (Complex a) -> Int -> (Ptr mx -> IO b) -> IO b
withImagDataOff a o f = withMXArray a (mxGetImagData >=> \p -> f (advancePtr p o))

foreign import ccall unsafe mxIsComplex :: MXArrayPtr -> IO CBool
mxArrayIsComplex :: MXArray a -> IO Bool
mxArrayIsComplex a = boolC =.< withMXArray a mxIsComplex

-- |Complex array access.
instance (RealFloat a, MNumeric a, MXArrayData mx a) => MXArrayComponent (MComplex a) where
  isMXArray a = liftM2 (&&) (isMXArray (castReal a)) (mxArrayIsComplex a)
  createMXArray s = withNDims s (uncurry $ createNumericArray (mxClassOf (undefined :: a)) True) >>= mkMXArray

  mxArrayGetOffset a o = do
    r <- withRealDataOff a o (mx2hs .=< peek)
    c <- withImagDataOff a o (mx2hs .=< peek)
    pure $ r :+ c
  mxArraySetOffset a o (r :+ c) = do
    withRealDataOff a o (\p -> poke p (hs2mx r))
    withImagDataOff a o (\p -> poke p (hs2mx c))
  mxArrayGetOffsetList a o n = do
    r <- withRealDataOff a o (map mx2hs .=< peekArray n)
    c <- withImagDataOff a o (map mx2hs .=< peekArray n)
    pure $ zipWith (:+) r c
  mxArraySetOffsetList a o v = do
    withRealDataOff a o (\p -> pokeArray p (map hs2mx r))
    withImagDataOff a o (\p -> pokeArray p (map hs2mx c))
    where (r,c) = unzip $ map (\(r:+c) -> (r,c)) v

ext: hsc

max_stars_repo_name: CornellCAC/haskell-matlab

max_issues_repo_name: CornellCAC/haskell-matlab

max_forks_repo_name: CornellCAC/haskell-matlab


text:
module Carlson
  (carlsonRF, carlsonRF',
  carlsonRD, carlsonRD',
  carlsonRJ, carlsonRJ',
  carlsonRC, carlsonRC',
  carlsonRG, carlsonRG')
  where
import           Data.Complex
import           Internal

rf_ :: Cplx -> Cplx -> Cplx -> Double -> ((Double,Double,Double), Cplx)
rf_ x y z err =
  let a = (x+y+z)/3 in
  let delta = map (\u -> magnitude(1-u/a)) [x,y,z] in
  if maximum delta < err
    then ((delta !! 0, delta !! 1, delta !! 2), a)
    else
      let (sqrtx, sqrty, sqrtz) = (sqrt x, sqrt y, sqrt z) in
      let lambda = sqrtx*sqrty + sqrty*sqrtz + sqrtz*sqrtx in
      rf_ ((x+lambda)/4) ((y+lambda)/4) ((z+lambda)/4) err

carlsonRF' :: Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRF' err x y z =
  if zeros > 1
    then error "At most one of x, y, z can be 0"
    else
      let ((dx,dy,dz), a) = rf_ x y z err in
      let (e2,e3) = (dx*dy + dy*dz + dz*dx, dx*dy*dz) in
      toCplx(1 - e2/10 + e3/14 + e2*e2/24 - 3*e2*e3/44 - 5*e2*e2*e2/208 +
          3*e3*e3/104 + e2*e2*e3/16) / sqrt a
  where
    zeros = sum (map (\u -> fromEnum (u == 0)) [x,y,z])

carlsonRF :: Cplx -> Cplx -> Cplx -> Cplx
carlsonRF = carlsonRF' 1e-15

rd_ :: Cplx -> Cplx -> Cplx -> Cplx -> Cplx -> Double ->
       ((Double,Double,Double), Cplx, Cplx, Cplx)
rd_ x y z s fac err =
  let a = (x+y+z+z+z)/5 in
  let delta = map (\u -> magnitude(1-u/a)) [x,y,z] in
  if maximum delta < err
    then ((delta !! 0, delta !! 1, delta !! 2), a, s, fac)
    else
      let (sqrtx, sqrty, sqrtz) = (sqrt x, sqrt y, sqrt z) in
      let lambda = sqrtx*sqrty + sqrty*sqrtz + sqrtz*sqrtx in
      let s' = s + fac / (sqrt z * (z + lambda)) in
      rd_ ((x+lambda)/4) ((y+lambda)/4) ((z+lambda)/4) s' (fac/4) err

carlsonRD' :: Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRD' err x y z =
  if zeros > 1
    then error "At most one of x, y, z can be 0"
    else
      let ((dx,dy,dz), a, s, fac) = rd_ x y z 0 1 err in
      let
        (e2,e3,e4,e5) = (dx*dy + dy*dz + 3*dz*dz + 2*dz*dx + dx*dz + 2*dy*dz,
          dz*dz*dz + dx*dz*dz + 3*dx*dy*dz + 2*dy*dz*dz + dy*dz*dz + 2*dx*dz*dz,
          dy*dz*dz*dz + dx*dz*dz*dz + dx*dy*dz*dz + 2*dx*dy*dz*dz,
          dx*dy*dz*dz*dz) in
      3*s + fac * toCplx(1 - 3*e2/14 + e3/6 + 9*e2*e2/88 - 3*e4/22 - 9*e2*e3/52 +
        3*e5/26 - e2*e2*e2/16 + 3*e3*e3/40 + 3*e2*e4/20 + 45*e2*e2*e3/272 -
        9*(e3*e4 + e2*e5)/68) / a / sqrt a
  where
    zeros = sum (map (\u -> fromEnum (u == 0)) [x,y,z])

carlsonRD :: Cplx -> Cplx -> Cplx -> Cplx
carlsonRD = carlsonRD' 1e-15

rj_ :: Cplx -> Cplx -> Cplx -> Cplx -> Cplx -> Double -> Cplx -> Int ->
       Double -> [Cplx] -> [Cplx] -> Double -> (Cplx, Int, [Cplx], [Cplx])
rj_ x y z p a maxmagns delta f fac d e err =
  let q = (4/err)**(1/6) * maxmagns / fromIntegral f in
  if magnitude a > q
    then (a, f, d, e)
    else
      let dnew = (sqrt p + sqrt x)*(sqrt p + sqrt y)*(sqrt p + sqrt z)
          d' = (fromIntegral f * dnew) : d
          e' = (toCplx fac * delta / dnew / dnew) : e
          lambda = sqrt x * sqrt y + sqrt y * sqrt z + sqrt z * sqrt x
          x' = (x + lambda) / 4
          y' = (y + lambda) / 4
          z' = (z + lambda) / 4
          p' = (p + lambda) / 4
          a' = (a + lambda) / 4
      in
      rj_ x' y' z' p' a' maxmagns delta (4*f) (fac/64) d' e' err

carlsonRJ' :: Double -> Cplx -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRJ' err x y z p =
  if zeros > 1
    then error "At most one of x, y, z, p can be 0"
    else
      let a0 = (x + y + z + p + p) / 5
          maxmagns = maximum $ map (\u -> magnitude(a0-u)) [x, y, z, p]
          delta = (p-x)*(p-y)*(p-z)
      in
      let (a, f, d, e) = rj_ x y z p a0 maxmagns delta 1 1 [] [] err
          f' = fromIntegral f
      in
      let x' = (a0 - x) / f' / a
          y' = (a0 - y) / f' / a
          z' = (a0 - z) / f' / a
          p' = -(x'+y'+z') / 2
          e2 = x'*y' + x'*z' + y'*z' - 3*p'*p'
          e3 = x'*y'*z' + 2*e2*p' + 4*p'*p'*p'
          e4 = p'*(2*x'*y'*z' + e2*p' + 3*p'*p'*p')
          e5 = x'*y'*z'*p'*p'
          h = zipWith (\u v -> atanx_over_x(sqrt u) / v) e d
      in
      (1 - 3*e2/14 + e3/6 + 9*e2*e2/88 - 3*e4/22 - 9*e2*e3/52 + 3*e5/26) /
        f' / a / sqrt a + 6 * sum h
  where
    zeros = sum (map (\u -> fromEnum (u == 0)) [x,y,z,p])
    atanx_over_x w = if w == 0 then 1 else atanC w / w

carlsonRJ :: Cplx -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRJ = carlsonRJ' 1e-15

rc_ :: Cplx -> Cplx -> Cplx -> Double -> Int -> Double -> (Cplx, Int)
rc_ x y a magn f err =
  let q = (1/3/err)**(1/8) * magn / fromIntegral f in
  if magnitude a > q
    then (a, f)
    else
      let lambda = 2 * sqrt x * sqrt y + y
          a' = (a + lambda) / 4
          x' = (x + lambda) / 4
          y' = (y + lambda) / 4
      in
      rc_ x' y' a' magn (4*f) err

carlsonRC' :: Double -> Cplx -> Cplx -> Cplx
carlsonRC' err x y =
  if y == 0
    then error "y cannot be 0"
    else
      let a0 = (x + y + y) / 3
          magn = magnitude(a0-x)
      in
      let (a, f) = rc_ x y a0 magn 1 err
          f' = fromIntegral f
      in
      let s = (y - a0) / f' / a in
      (1 + 3*s*s/10 + s*s*s/7 + 3*s*s*s*s/8 + 9*s*s*s*s*s/22 +
        159*s*s*s*s*s*s/208 + 9*s*s*s*s*s*s*s/8) / sqrt a

carlsonRC :: Cplx -> Cplx -> Cplx
carlsonRC = carlsonRC' 1e-15

carlsonRG' :: Double -> Cplx -> Cplx -> Cplx -> Cplx
carlsonRG' err x y z =
  if zeros > 1
    then sqrt(x+y+z) / 2
    else
      if z == 0
        then carlsonRG' err z x y
        else
          (z * carlsonRF' err x y z -
            (x-z) * (y-z) * carlsonRD' err x y z / 3 +
            sqrt x * sqrt y / sqrt z) / 2
  where
    zeros = sum (map (\u -> fromEnum (u == 0)) [x,y,z])

carlsonRG :: Cplx -> Cplx -> Cplx -> Cplx
carlsonRG = carlsonRG' 1e-15

ext: hs

max_stars_repo_name: stla/elliptic

max_issues_repo_name: stla/elliptic

max_forks_repo_name: stla/elliptic


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

module BPTT where

import Control.Monad
import Control.Monad.Random
import Data.Char
import Data.List
import Data.Maybe
import Data.Singletons
import Data.Singletons.Prelude
import Data.Singletons.TypeLits
import GHC.Show
import Numeric (showIntAtBase)
import qualified Numeric.LinearAlgebra as LA
import Numeric.LinearAlgebra.Static

-- dependently typed weights
data Weights i o = W {
    wBiases :: !(R o)
  , wNodes  :: !(L o i)
  }

-- an simple 3-layer RNN, including one recurrent layer,
-- is not like the Network type discussed in Justin Le's
-- blog series (because not all layers are just stacked
-- one on top of the other.)
-- Instead:
data RNN i h o = MkRNN {
    getU    :: !(Weights i h)
  , getW    :: !(Weights h h)
  , getV    :: !(Weights h o)
  }

-- If one were to leverage a.g. ad the Hackage library for
-- autodifferentiable functions, then functions being used for
-- activation could be provided in a symbolic form amenable to
-- being automatically differentiated (for use in the backprop
-- process.) Instead, each activation function used & its derivative
-- are provided separately.

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

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

softmax :: KnownNat v => R v -> R v
softmax v =
  let vexp = exp v
  in  vexp / (konst $ norm_1 vexp)

-- The formula for the (x, y)-entry of softmax',
-- which is a matrix, is
--
-- softmax'[x, y] = -softmax[x]    * softmax[y],    x /= y
--                  (1-softmax[x]) * softmax[y],    x == y
--
-- which is the same as adding the vector
-- 'softmax v' diagonally down the negative of
-- said vector's outer product with itself.
-- I would prefer to build this using 'build', but
-- I don't know how right now.
softmax' :: KnownNat v => R v -> L v v
softmax' v =
  let s = softmax v
  in  (diag s) + konst (-1) * (s `outer` s)

-- runLayer is for performing the linear mapping
-- due to a particular layer's set of parameters,
-- on an appropriately sized vector
runLayer :: (KnownNat i, KnownNat o)
         => Weights i o
         -> R i -> R o
runLayer (W wB wN) v_in =
  wB + wN #> v_in

-- Verbose 'run RNN', for exposing all intermediate computations.
-- (R i, R h,  R h, R o, R o)
--  x_t  h_t-1 h_t  z_t  y_t

-- runRNN can be expressed in terms of this, if I locate / find
-- generalized tuple accessors
runRNN_verbose  :: (KnownNat i, KnownNat h, KnownNat o)
                => RNN i h o
                -> R i -> R h
                -> (R i, R h, R h, R o, R o)
runRNN_verbose (MkRNN rnn_U rnn_W rnn_V)
  x_t h_t_1 =
    let h_t = logistic $
                -- recurrent parameters, weighting the prior hidden state
                (runLayer rnn_W h_t_1) +
                -- input parameters, weighting the current input
                (runLayer rnn_U x_t)
        -- output parameters, weighting the current hidden state
        z_t = runLayer rnn_V h_t
        -- the final activation provides the prediction
        y_t = logistic z_t
    in  (x_t, h_t_1, h_t, z_t, y_t)

-- Forward phase of training an RNN.
-- No adjustments to the parameters (weights) are made.
-- Instead, foldl is used to accumulate a list of tuples,
-- each tuple bearing all input, intermediate and output vectors
-- for one of the steps (t := 1, 2, ..., n) in the forward direction.
forward_phase :: (KnownNat i, KnownNat h, KnownNat o)
              => RNN i h o    -- RNN
              -> [(R i, R o)] -- list of training inputs paired with target outputs
              -> R h          -- initial hidden state

                              -- All tuples of resulting vectors, from every step carried out:
                              -- For t ranging from 1 to n, n the number of training samples:
              -> [( R i,      -- x_t
                    R h,      -- h_t-1
                    R h,      -- h_t
                    R o,      -- z_t
                    R o,      -- y_t
                    R o)]     -- tgt_t
forward_phase
  rnn
  pairs
  h_0 =
    snd $ foldl unroll (h_0, []) pairs
    where
      -- unroll :: (R h, [(R i, R h, R h, R o, R o, R o)]) -> (R i, R o)
      --        -> (R h, [(R i, R h, R h, R o, R o, R o)])

      unroll (h_t_1, results) pair@(x_t, tgt_t) =
        --                               ^ tgt_t plays no role in the forward direction,
        --                                 but we pattern-match on it so that we can re-package it
        --                                 up into this step's tuple for the result list
        let result@(_, _, h_t, z_t, y_t) = runRNN_verbose rnn x_t h_t_1
        in  (h_t, results ++ [(x_t, h_t_1, h_t, z_t, y_t, tgt_t)])

-- Backward phase of training an RNN.
-- Backpropagation is required.
-- There's a lot that goes into backprop!
-- (Basically, everything that is used or
-- produced at any point in the forward step.)
backward_phase  :: (KnownNat i, KnownNat h, KnownNat o)
                => RNN i h o   -- an RNN to train

                                -- All tuples of resulting vectors, from every step carried out:
                                -- For t ranging from 1 to n, n the number of training samples:
                -> [( R i,      -- x_t
                      R h,      -- h_t-1
                      R h,      -- h_t
                      R o,      -- z_t
                      R o,      -- y_t
                      R o)]     -- tgt_t

                -- At every step #t of BPTT, there are *two* sources of "delta-h_t" the bump
                -- to h_t required for the sake of trying to bump down the error at that same step:
                -- a) a bump due to what W is at the current time, "delta-h given W", and
                -- b) a bump *propagated backwards* from the next step #(t+1).
                -- But in the base case, this dh_t_dh_t' is chosen as 0 (because it's undefined.)

                -- Note the *symmetry* with forward_phase:
                --    just as we 'primed' forward_phase  with an initial hidden state,    h_0        : 0,
                --    we also    'prime'  backward_phase with an initial back-prop'ed dh, dh_n_dh_n' : 0.
                -> R h        -- dh_n_dh_n'

                -> Double     -- a training rate
                -> RNN i h o  -- the trained RNN
backward_phase rnn intermediates dh_n_dh_n' rate =
  snd $ foldr bptt (dh_n_dh_n', rnn) intermediates
    where
      -- -- A BPTT step (at time t) consumes the following data:
      -- --       x_t, h_t, h_t_1, z_t,  y_t, tgt_t
      -- bptt :: (R i, R h, R h,   R o,  R o, R o)
      -- -- It also "accumulates" by swapping out,
      -- -- based upon time-t data  <dh_t   | dh_t+1>, rnn
      --          -> (R h, RNN i h o)
      -- -- the time-(t-1) data     <dh_t_1 | dh_t>,   rnn';
      -- -- in this way, the rnn' "accumulates" training bumps to the parameters
      --          -> (R h, RNN i h o)

      bptt    (x_t, h_t, h_t_1, z_t, y_t, tgt_t)
              (dh_t_dh_t', MkRNN
                -- We'll pattern match into each Weights term u, w, v for convenience
                -- later in deriving the desired infinitesimal bumps.
                u@(W uB uN)
                w@(W wB wN)
                v@(W vB vN)) =
                -- We are hard-coding the use of squared-error for the error function:
                -- e(y_t) := (y_t - tgt_t) ^ 2
                -- whose derivative with respect to the pre-activation z_t is:

                {-
                  de/dz = de/dy           * dy/dz
                        = (e' at y_t)     * (y' at z_t)
                        = 2 (y_t - tgt_t) * (logistic' z_t)
                  These 2 factors appear (in reverse order) as 'dedz' below;
                  the factor of 2 is discarded & we defer all control of the
                  desired training rate to the scalar input param 'rate'.
                -}
            let dedz          = logistic' z_t * (y_t - tgt_t)
                              -- = softmax' z_t #> (y_t - tgt_t)

                {-
                  Informally:
                  dedz is a vector that reports,
                    for a unit "nudge" to each component z_t(j) ->
                      size / direction of resulting "nudge" for error component e_t(j).
                  The choice of dz_t below corresponds to the following training strategy:
                    a) bump each z-component in the direction that will *decrease* its
                       matching e-component (thus, factor of '-1')
                    b) but re-scaled based on the learning rate desired (thus, 'scale rate')
                -}
                dz_t          = -1 * konst rate * dedz

                {-
                  Informally:
                  Having set out to make these infinitesimal changes to the z_t components,
                  we need to figure out
                    what changes to V,    *given* what h_t is,     and
                    what changes to h_t,  *given* what V is
                  will help us accomplish that.

                  Observe that
                    z  = vN h              + vB
                  therefore (product rule on 1st summand)
                    dz = (dvN h) + (vN dh) + dvB

                  we do *not* try to formally solve for dvN, dh, dvB in here, but rather use these
                  as instructions for how to recombine the vectors dz_t / h_t and matrix V
                  so as to arrive at dvN_h ("bump to v's matrix part, given h") and
                  dh_vN ("bump to h, given v's matrix part").

                  As for v's bias-vector part, vB, we simply nudge it directly by the dz (z-nudge)
                  we worked out above. Why? Look at the infinitesimals equation above - a nudge
                  dvB performed on the bias vector will show up directly in dz.
                -}
                dvB           = dz_t
                {-
                  Instead of trying to solve for dvN, dimensional considerations lead us to this formula
                  for dvN in terms of dz_t, h_t:
                    vN   goes    from h-vectors  to o-vectors; i.e.    o * h -matrix
                    h_t            an h-vector
                    dz_t                         an o-vector
                    so the way to combine these to get a "bumps matrix" compatible with vN is
                    dz_t `outer` h_t   <- "column" o-vector times "row" h-vector
                -}
                dvN_h         = dz_t `outer` h_t
                {-
                  Similar to the above, but for building a "bumps h-vector" out of o * h matrix
                  and o-vector.
                -}
                dh_vN         = tr vN #> dz_t
                {-
                  See the note in the type signature of 'backward_phase' about the 2 sources of
                  training bumps to h_t at every step t.
                  Read dh_t_dh_t' as (bump to h at current time t,
                                      propped backwards from bump to h at future time t+1)
                -}
                dh_t          = dh_vN + dh_t_dh_t'

                -- Having now characterized in what way we are going to bump the
                -- entries in the current hidden state h_t, so as to achieve
                -- a desired reduction in error (i.e. dh_t), we now use it to
                -- determine how we should modify the intermediate preactivation,
                -- calling it q, q_t := (W * h_t_1 + U * x_t) to achieve the same goal:

                -- de / dq = (de / dh) * (dh / dq)  -- chain rule
                --         = dh_t      * (dh / dq)  -- back propagation
                --         = dh_t      * (g' at q)  -- h is the activation (g) of q

                dedq          = dh_t * logistic' q
                                  where q = runLayer w h_t_1 + runLayer u x_t

                -- and work out a desired perturbation to q, based upon scaling dedq
                dq_t          = (-1) * konst rate * dedq

                {-
                  q  = W * h_t_1 + U * x_t
                  dq = (dW * h_t_1) + (W * dh_t_1)
                        + (dU * x_t) + [the term in dx_t -> 0, since x constant]
                -}

                {-
                  Recurrent contribution to time-t hidden state:
                  q_t  = wN h_t_1                   + wB        + ...
                  dq_t = (dwN h_t_1) + (wN dh_t_1)  + dwB       + ...
                -}
                dwB           = dq_t
                -- dq_t and h_t_1 both happen to be h-vectors having the same
                -- dimensionality, so here the driving intuition is that such
                -- "matrix-bumps" take the form [destination] `outer` [origin]
                dwN_h         = dq_t `outer` h_t_1

                -- This is what will be given back from the current fold instance
                -- to the next one
                dh_t_1_dh_t   = tr wN #> dq_t
                {-
                  Input contribution to time-t hidden state:
                  q_t  = ... + uN x_t                 + uB
                  dq_t = ... + (duN x_t) + (uN dx_t)  + duB
                       = ... + (duN x_t) + ZERO       + duB
                                        -- ^ x_t, the training inputs, are constant!
                -}
                duB           = dq_t
                duN_x         = dq_t `outer` x_t

                -- Don't forget to apply all the RNN-training nudges!
                vB'           = vB + dvB
                vN'           = vN + dvN_h
                wB'           = wB + dwB
                wN'           = wN + dwN_h
                uB'           = uB + duB
                uN'           = uN + duN_x

            in  (dh_t_1_dh_t, MkRNN (W uB' uN')
                                    (W wB' wN')
                                    (W vB' vN'))

{-
 Making some simple training data.

 Here's a simple sentential grammar for producing "sentences"
 of ints, that later can go in an overall "document":

  <sentence>        ::= <sentence-start> <sentence-middle> <sentence-finish>
  <sentence-start>  ::= a multiple of 20, 20 * s, for some s : 1 <= s <= 10
  <sentence-middle> ::= a sequence of l integers, for some l : 0 <= l <= 10
                        that are derived from the sentence-start as
                          s - m, (s - m) + n, (s - m + n) - m, (s - m + n - m) + n, ...
                        for some choice of parameters      m, n : 1 <= m, n <= 5
  <sentence-finish> ::= a multiple of 17, 17 * f, for some f : 1 <= f <= 10
-}

-- Generate a "biarithmetic" sequence for sentence-middle, given
-- a base number s, a down size m, an up size n
biarithmetic :: Int -> Int -> Int -> [Int]
biarithmetic    s      m      n   =  concat $
 [ [ s - (j * m) + (k * n) | k <- [j-1 .. j] ]
                           | j <- [1..] ]

-- For use as a multiple of 20
genStartParam :: MonadRandom m => m Int
genStartParam = getRandomR (3, 10)

-- For use as a mid-sequence length
genMiddleLenParam :: MonadRandom m => m Int
genMiddleLenParam = getRandomR (0, 10)

-- For use as the down size and up size of a biarithmetic sequence
genMiddleDownUpParams :: MonadRandom m => m (Int, Int)
genMiddleDownUpParams = do
  m <- getRandomR (1, 5)
  n <- getRandomR (1, 5)
  return (m, n)

-- For use as a multiple of 17
genFinishParam :: MonadRandom m => m Int
genFinishParam = getRandomR (1, 10)

-- For use as a number of "sentences".
genSentenceCount :: MonadRandom m => m Int
genSentenceCount = getRandomR (500, 1000)

type Sentence = [Int]

genSentence :: MonadRandom m => m Sentence
genSentence = do
  -- 1. Get sentence-start
  s       <- genStartParam
  let sent_0 = s * 20
  -- 2. Get sentence-middle length
  l       <- genMiddleLenParam
  -- 3. Get sentence up & down size params
  (m, n)  <- genMiddleDownUpParams
  -- 4. Get sentence-middle
  let sent_mid = take l $ biarithmetic sent_0 m n
  -- 5. Get sentence-finish
  f       <- genFinishParam
  let sent_f = f * 17
  -- 6. Return random sentence
  return (sent_0 : sent_mid ++ [sent_f])

type Document = [Sentence]

genDocument :: MonadRandom m => m Document
genDocument = do
  -- 1. Get number of sentences
  sent_count <- genSentenceCount
  -- 2. For as many sentences as 'sent_count', generate one
  --    random sentence each & then combine the results inside
  --    MonadRandom context
  let rand_sent_stream = repeat genSentence
  sentences <- sequenceA $ take sent_count rand_sent_stream
  return sentences

{-
 Now that we have a way of generating "documents" made up out of
  numerical "sentences", we need to convert this into a form that
  would be characteristically used to train a neural net.
  Convert each integer to a (possibly sparse) vector. How?
  Every integer present in a (d :: Document) should always lie
  within the range (0 <= j <= 256), so we could use a binary representation
  taking up 8 bits (tossing in an extra 2 to get 10, for good measure)
-}

-- The full source module citation is probably forgivable if I only do it this once.
-- Best practice is to 'import qualified' with an alias intToVector :: Int -> R 10
intToVector m =
  let binary_string = reverse $ showIntAtBase 2 intToDigit m "" -- :: [Char]
      binary_padded = binary_string ++
                        (take (10 - length binary_string) $ repeat '0')
      doubles_list  = (fromIntegral . digitToInt) <$> binary_padded
      doubles_vect  = Numeric.LinearAlgebra.Static.fromList doubles_list
  in  doubles_vect

-- The other way round, for casting RNN predictions back into Int's.
-- This is unsafe in case the entries of any y_t end up rounding to
-- numbers other than 0, 1.
-- (Does the hmatrix API support a more direct way of casting from
-- R 10 to the list of rounded numbers?)
vectorToInt :: R 10 -> Int
vectorToInt v =
  let bins = (fmap toInteger) <$> LA.toList
                                $ LA.toZ
                                $ LA.roundVector
                                $ extract v
  in  fromIntegral $
        foldr (\p@(bit,  power)  total ->
                   bit * power + total)
        0
        (zip bins [ 2^j | j <- [0..] ])

zero_vec :: KnownNat n => R n
zero_vec = 0

-- via Justin Le's blog series on neural nets
randomWeights :: (MonadRandom m, KnownNat i, KnownNat o)
              => m (Weights i o)
randomWeights = do
    s1 :: Int <- getRandom
    s2 :: Int <- getRandom
    let wB = randomVector  s1 Uniform * 2 - 1
        wN = uniformSample s2 (-1) 1
    return $ W wB wN

randomRNN :: (MonadRandom m, KnownNat i, KnownNat h, KnownNat o)
          => m (RNN i h o)
randomRNN =
  (pure MkRNN)  <*> randomWeights
                <*> randomWeights
                <*> randomWeights

trainRNN :: (KnownNat i, KnownNat h, KnownNat o)
         => RNN i h o     -- an RNN
         -> [(R i, R o)]  -- the training data: pairs of inputs & training outputs
         -> Double
         -> RNN i h o     -- an RNN with trained parameters
trainRNN rnn train_data rate =
  let intermediates = forward_phase rnn train_data zero_vec
      trained       = backward_phase rnn intermediates zero_vec rate
  in  trained

-- Before re-formatting into a vector form suitable for use
-- in training, we have Document ~ [Sentence] ~ [[Int]];
-- now, turn each Int into a binary-repr vector, R 10
genSampleData :: MonadRandom m => m [[ R 10 ]]
genSampleData = do
  doc <- genDocument
  return $ (intToVector <$>) <$> doc

-- Every sentence (sublist) in the document (list of lists) will
-- furnish one sequence of training data for the RNN to "unfold"
-- over and perform BPTT on. So, BPTT is performed once per sentence
-- in the document overall.
pairEachWithNext :: [x] -> [(x, x)]
pairEachWithNext xs =
  zip (init xs) (tail xs)

-- In the special case of an 'RNN i h o' where o = i,
-- we can pipe back a prediction y_t due to x_t and an
-- initial state h, to furnish the next x_(t+1) for further
-- predictions. We'll use this to observe the way that a
-- trained RNN behaves.
rnnStep :: (KnownNat i, KnownNat h)
        => RNN i h i
        -> (R h, R i)
        -> (R h, R i)
rnnStep rnn (h_t_1, x_t) =
  let full_output = runRNN_verbose rnn x_t h_t_1
      third five@(_, _, x3, _, _)  = x3
      fifth five@(_, _, _,  _, x5) = x5
  in  (third full_output, fifth full_output)

getDocumentIO :: IO Document
getDocumentIO = evalRandIO genDocument

takePredictions  :: (KnownNat i, KnownNat h)
                 => RNN i h i
                 -> (R h, R i)
                 -> Int
                 -> [R i]
takePredictions rnn (h_0, x_1) n =
  let hxs_stream = take (n+1) $ iterate (rnnStep rnn) (h_0, x_1)
  in  snd <$> hxs_stream

takePredsFromN :: KnownNat h
               => RNN 10 h 10
               -> Int      -- word of a sentence
               -> Int      -- number of predictions
               -> [R 10]
takePredsFromN rnn x_1 n =
  takePredictions rnn (zero_vec, intToVector x_1) n

trainRNNOnSentence  :: KnownNat h
                    => RNN 10 h 10
                    -> Sentence
                    -> Double
                    -> RNN 10 h 10
trainRNNOnSentence rnn sent rate =
  let training_data = pairEachWithNext (intToVector <$> sent)
  in  trainRNN rnn training_data rate

-- trainRNNOnDocument has to be structured like a fold
-- so that the trained RNN due to each sentence can be
-- passed to the following sentence (at the same time,
-- I think that hidden state "restarts" with the beginning
-- of a new sentence.)
trainRNNOnDocument  :: KnownNat h
                    => RNN 10 h 10
                    -> Document
                    -> Double
                    -> RNN 10 h 10
trainRNNOnDocument rnn sents rate =
  foldl (\rnn sent -> trainRNNOnSentence rnn sent rate)
        rnn
        sents

observeUntrainedPredictions :: Int    -- initial word
                            -> Int    -- number of predictions
                            -> IO ()  -- action printing info to console
observeUntrainedPredictions n_0 num = do
  -- get an RNN of random parameters
  putStrLn "Getting a random RNN ..."
  (rnn :: RNN 10 300 10) <- evalRandIO randomRNN
  -- get a series of predictions from the untrained RNN
  putStrLn $ "Producing: " ++ show num ++
             " predictions, starting from: " ++ show n_0
  let vs = takePredsFromN rnn n_0 num
  let ns = vectorToInt <$> vs
  putStrLn "The RNN has predicted: "
  print ns

observeTrainedPredictions :: Int    -- initial word
                          -> Int    -- number of predictions
                          -> Double -- learning rate
                          -> IO ()  -- action printing info to console
observeTrainedPredictions n_0 num rate = do
  -- get an RNN of random parameters
  putStrLn "Getting a random RNN ..."
  (rnn :: RNN 10 300 10) <- evalRandIO randomRNN
  -- get a random document
  putStrLn "Getting a random document to train on ..."
  document <- evalRandIO genDocument
  -- train RNN on the document
  putStrLn "Training the RNN on the document ..."
  let rnn' = trainRNNOnDocument rnn document rate
  -- get a series of predictions from the trained RNN
  putStrLn $ "Producing: " ++ show num ++
             " predictions, starting from: " ++ show n_0 ++
             " and using a trained RNN"
  let vs = takePredsFromN rnn' n_0 num
  let ns = vectorToInt <$> vs
  putStrLn "The trained RNN has predicted: "
  print ns

observeUntrainedVsTrained :: Int    -- initial word
                          -> Int    -- number of predictions
                          -> Double -- learning rate
                          -> IO ()  -- action printing info to console
observeUntrainedVsTrained n_0 num rate = do
  -- get an RNN of random parameters
  putStrLn "Getting a random RNN ..."
  (rnn :: RNN 10 300 10) <- evalRandIO randomRNN
  -- get a series of predictions from the untrained RNN
  putStrLn $ "Producing: " ++ show num ++
             " predictions, starting from: " ++ show n_0
  let vs = takePredsFromN rnn n_0 num
  let ns = vectorToInt <$> vs
  putStrLn "The RNN has predicted: "
  print ns
  -- get a random document
  putStrLn "Getting a random document to train on ..."
  document <- evalRandIO genDocument
  -- train RNN on the document
  putStrLn "Training the RNN on the document ..."
  let rnn' = trainRNNOnDocument rnn document rate
  -- get a series of predictions from the trained RNN
  putStrLn $ "Producing: " ++ show num ++
             " predictions, starting from: " ++ show n_0 ++
             " and using a trained RNN"
  let vs' = takePredsFromN rnn' n_0 num
  let ns' = vectorToInt <$> vs'
  putStrLn "The trained RNN has predicted: "
  print ns'

max_stars_repo_name: enlambdment/ffnn-vs-rnn

max_issues_repo_name: enlambdment/ffnn-vs-rnn

max_forks_repo_name: enlambdment/ffnn-vs-rnn


text:
{-# LANGUAGE BangPatterns    #-}
{-# LANGUAGE NamedFieldPuns  #-}
{-# LANGUAGE RecordWildCards #-}
module Network
    ( NeuralNet (..),
      train,
      initWB,
      loadData,
      loadTestData,
    ) where

import           Codec.Compression.GZip (decompress)
import           Control.Monad
import qualified Data.Array.IO          as A
import qualified Data.ByteString        as BS
import qualified Data.ByteString.Lazy   as BL
import           Data.Int               (Int64)
import           Data.List
import           GHC.Base               (build)
import           Numeric.LinearAlgebra  hiding (build, normalize)
import           Prelude                hiding ((<>))
import           System.Directory       (getCurrentDirectory)
import qualified System.Random          as R
import           Text.Printf
import           Util

type Image        = Vector Double
type Label        = Vector Double
type Bias         = Vector Double
type Weight       = Matrix Double
type Activation   = Matrix Double
type Nabla        = Matrix Double
type ImagesMat    = Matrix Double
type LabelsMat    = Matrix Double

data NeuralNet = NN
    { weights   :: ![Weight]   -- ^ List of all weights
    , biases    :: ![Bias]    -- ^ List of all biases
    , eta       :: !Double    -- ^ Learning rate
    , epochs    :: !Int       -- ^ No. of epochs
    , layers    :: !Int       -- ^ No. of layers
    , layerSize :: !Int       -- ^ Size of hidden layers
    , batchSize :: !Int       -- ^ Size of mini-batch
    , trainData :: ![(Image, Label)]
    , testData  :: ![(ImagesMat, Int)]
    } deriving (Show, Eq)

test :: NeuralNet -> IO ()
test NN{..} = do
    let (imgs, labels) = unzip testData
    let guesses = [guess img (zip weights $ toMatrix <$> biases) | img <- imgs]
    let correct = sum $ (\(a,b) -> if a == b then 1 else 0) <$> zip guesses labels
    printf "Epoch %d: %.2f%%\n" (30-epochs :: Int)  (correct/100 :: Float)
  where
    guess :: ImagesMat -> [(Weight, Matrix Double)] -> Int
    guess i wb = maxIndex . flatten . head $ feedforward i wb

train :: NeuralNet -> IO NeuralNet
train net@NN{..} = do
    shuffledData <- shuffle trainData
    let miniBatches = getMiniBatch <$> chunkList batchSize shuffledData
    let newNet = foldl trainBatch net miniBatches
    test newNet
    case epochs of
        0 -> return newNet
        _ -> train newNet { epochs=epochs-1 }

trainBatch :: NeuralNet -> (ImagesMat, LabelsMat) -> NeuralNet
trainBatch net@NN{weights=ws, biases=bs, ..} (mX, mY) = newNet
    where newNet = let activations = feedforward mX (zip ws $ bTm bs batchSize)
                      in updateMiniBatch net $
                          backprop mX mY net activations

-- | Make biases suitable for full-matrix backprop.
bTm :: [Bias] -> Int -> [Matrix Double]
bTm biases batchSize = fmap (fromColumns . replicate batchSize) biases

-- | Turn mini-batch matrices into one matrix.
getMiniBatch :: [(Image, Label)] -> (ImagesMat, LabelsMat)
getMiniBatch mBatchData = (fromColumns imgs, fromColumns labels)
  where
    (imgs, labels) = unzip mBatchData

feedforward :: ImagesMat -> [(Weight, Matrix Double)] -> [Activation]
feedforward = scanr (\(w, b) a -> sigmoid $ w <> a + b)

--  | Stochastic gradient descent
sgd = id

--  | Updates mini-batch
--  NOTE: use 'seq' or 'deepseq' for full evaluation,
--  instead of just reducing to WHNF
updateMiniBatch :: NeuralNet -> ([Vector Double], [Matrix Double])  -> NeuralNet
updateMiniBatch net@NN{..} (nablaB, nablaW) = net {weights=wnew, biases=bnew}
    where
        wnew :: [Weight]
        !wnew = [ w - scale (eta/fromIntegral batchSize) nw |
            (w, nw)  <- zip weights nablaW ]
        bnew :: [Bias]
        !bnew = [ b - scale (eta/fromIntegral batchSize) nb |
            (b, nb)  <- zip biases nablaB ]

--  | Returns nablaW and nablaB
backprop :: ImagesMat -> LabelsMat -> NeuralNet -> [Activation] -> ([Vector Double], [Matrix Double])
backprop x y NN{..} (aL:as) = (nablaB, nablaW)
    where
        nablaW :: [Matrix Double]
        nablaW = zipWithSafe (<>) delta (tr <$> as) -- THIS IS WHERE THE MEMORY LEAK IS IF I DO JUST zipWith
        nablaB :: [Vector Double]
        nablaB = sumRows <$> delta
        delta :: [Matrix Double]
        delta = init $ deltas [crossEntropy y aL] y as weights

deltas :: [Matrix Double] -> LabelsMat -> [Activation] -> [Weight] -> [Matrix Double]
deltas deltal _ [] [] = deltal
deltas deltal y (a:as) (wlp1:ws) = deltas (deltal++[(tr wlp1 <> last deltal) * sigmoid' a]) y as ws

quadraticLoss :: LabelsMat -> Activation -> Matrix Double
quadraticLoss y aL = (aL - y) * sigmoid' aL

crossEntropy :: LabelsMat -> Activation -> Matrix Double
crossEntropy y aL = (aL - y)

{-- Arithmetic Functions --}
cost :: Matrix Double -> Matrix Double -> Matrix Double
cost a y = ((a-y)**2) / 2

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

sigmoid' :: Floating a => a -> a
sigmoid' x = x * (1 - x)

normalize :: (Integral a, Floating b) => a -> b
normalize x = (fromIntegral x) / 255

{-- Data Processing Functions --}
initWB :: NeuralNet -> IO NeuralNet
initWB net@NN{batchSize = bs, layerSize = ls, layers = n} = do
    headW <- randn ls 784
    midW  <- replicateM (n-2) (randn ls ls)
    lastW <- randn 10 ls

    initB <- replicateM (n-1) (randn ls 1)
    lastB <- randn 10 1

    let weights = reverse $ headW : midW ++ [lastW]
    let biases = reverse . fmap flatten $ initB ++ [lastB]

    return net { weights = weights
               , biases  = biases }

loadData :: IO [(Vector Double, Vector Double)]
loadData = do
    trainImgs <- getData "train-images-idx3-ubyte.gz"
    trainLabels <- getData "train-labels-idx1-ubyte.gz"
    let !labels = BL.toStrict trainLabels
    let !imgs = BL.toStrict trainImgs
    let !l = [vectorizeLabel $ getLabel n labels | n <- [0..49999]] :: [Label]
    let !i = [getImage n imgs | n <- [0..49999]] :: [Image]
    return $ zip i l

loadTestData :: IO [(Matrix Double, Int)]
loadTestData = do
    testImgs <- getData "t10k-images-idx3-ubyte.gz"
    testLabels <- getData "t10k-labels-idx1-ubyte.gz"
    let !labels = BL.toStrict testLabels
    let !imgs = BL.toStrict testImgs
    let !l = [getLabel n labels | n <- [0..9999]] :: [Int]
    let !i = asColumn <$> [getImage n imgs | n <- [0..9999]] :: [Matrix Double]
    return $ zip i l

getData :: FilePath -> IO BL.ByteString
getData path = do
    currentDir <- getCurrentDirectory
    fileData <- decompress <$> BL.readFile (currentDir ++ "/data/mnist_dataset/" ++ path)
    return fileData

getImage :: Int -> BS.ByteString -> Image
getImage n imgs = fromList [normalize $ BS.index imgs (16 + n*784 + s) | s <- [0..783]]

-- | `getImage` but lazy
getImage' :: Int64 -> BL.ByteString -> Image
getImage' n imgs = fromList [normalize $ BL.index imgs (16 + n*784 + s) | s <- [0..783]]

getLabel :: Num a => Int -> BS.ByteString -> a
getLabel n labels = fromIntegral $ BS.index labels (n+8)

getGuess :: [Activation] -> Int
getGuess = maxIndex . flatten . head

vectorizeLabel :: Int -> Vector Double
vectorizeLabel l = fromList $ x ++ 1 : y
    where (x,y) = splitAt l $ replicate 9 0

{-- HELPER FUNCTIONS --}
toMatrix :: Vector Double -> Matrix Double
toMatrix = reshape 1

-- | Means of rows in matrix as a column vector
meanRows :: Matrix Double -> Vector Double
meanRows = fst . meanCov . tr

sumRows :: (Num a, Element a) => Matrix a -> Vector a
sumRows m = fromList $ sum <$> toLists m

-- | zipWith, but throws error if lengths don't match
zipWithSafe :: (a -> b -> c) -> [a] -> [b] -> [c]
zipWithSafe f as bs
    | length as == length bs = zipWith f as bs
    | otherwise = error $ "zipWith: lengths don't match (" ++ (show $ length as) ++ " " ++ (show $ length bs) ++ ")"

-- | Randomly shuffle a list
--   /O(N)/
shuffle :: [a] -> IO [a]
shuffle xs = do
        ar <- newArray n xs
        forM [1..n] $ \i -> do
            j <- R.randomRIO (i,n)
            vi <- A.readArray ar i
            vj <- A.readArray ar j
            A.writeArray ar j vi
            return vj
  where
    n = length xs
    newArray :: Int -> [a] -> IO (A.IOArray Int a)
    newArray n xs =  A.newListArray (1,n) xs

max_stars_repo_name: ttesmer/haskell-mnist

max_issues_repo_name: ttesmer/haskell-mnist

max_forks_repo_name: ttesmer/haskell-mnist


text:
module Eris.Pantheon.Minkowski where

import Eris.Meta.DataTypes
import Numeric.LinearAlgebra hiding (Vector)
import qualified Data.Vector.Storable as DV

-- import Numeric.LinearAlgebra.Data
{-
   This module contains different derivation of Minkowski distance.
   These functions define distance between two high dimensional vectors X and Y.
   One distance definition may have many different names .

   1. sad (sub absolute difference, taxicab metric, manhattan distance, l1-norm of (X-Y) )
-}
{- L1 norm family -}
-- Sum of Absolute Difference: norm(l1) of (v1 - v2)
-- sad , taxicab, manhattan distance are different names of the same implementation.
sad :: Vector -> Vector -> Double
sad v1 v2 = norm_1 (v1 - v2)

taxicab :: Vector -> Vector -> Double
taxicab = sad

manhattanDistance :: Vector -> Vector -> Double
manhattanDistance = sad

-- | Mean Absolute Difference
mad :: Vector -> Vector -> Double
mad v1 v2 = norm1 / n
  where
    norm1 = sad v1 v2
    n = fromIntegral . DV.length $ v1

{- L2 norm based -}
-- | Euclidean Distance
--
euclideanDistance :: Vector -> Vector -> Double
euclideanDistance v1 v2 = norm_2 (v1 - v2)

-- | (L2 norm based) Sum Suqred Difference
-- Squred euclidean distance
ssd :: Vector -> Vector -> Double
ssd v1 v2 = v <.> v
  where
    v = v1 - v2

-- | (L2 norm based) Mean Sqaured Error (MSE)
--  ssd averaged over the vector
mse :: Vector -> Vector -> Double
mse v1 v2 = squaredL2 / n
  where
    squaredL2 = ssd v1 v2
    n = fromIntegral . DV.length $ v1

-- | (L2 norm based) Root mean squared Error (rmse)
-- euclidean distance averged over the vector
rmse :: Vector -> Vector -> Double
rmse v1 v2 = euclideanDistance v1 v2 / sqrtN
  where
    sqrtN = sqrt . fromIntegral $ DV.length v1

-- | Minkowski Distance
-- This is a degenerate case of Minkowski distance
-- in this case p is an Integral number for sure.
minkowskiDistance :: Int -> Vector -> Vector -> Double
minkowskiDistance p v1 v2 =
  let vDiff = v1 - v2
      absV = DV.map abs vDiff
  in lpnorm p absV
  where
    lpnorm :: Int -> Vector -> Double
    lpnorm 0 vec =
      let vl = -DV.length vec
      in norm_1 $ DV.map (f0 vl) vec
    lpnorm p vec = nroot p $ norm_1 $ DV.map (\n -> n ^ p) vec
    nroot
      :: (Integral a, Floating b)
      => a -> b -> b
    nroot 0 _ = 1
    nroot n f = f ** (1 / fromIntegral n)
    f0 :: Int -> Double -> Double
    f0 l v = (2 ^ l) * (v / (1 + v))

{- Minkowski based (p = 0, p -> infinity) -}
-- | Chebyshev Distance
-- The limiting case of Minkowski distance when p reaches infinity
cbv :: Vector -> Vector -> Double
cbv v1 v2 =
  let absV = DV.map abs $ v1 - v2
  in DV.maximum absV

-- | Canberra Distance
--
cbr :: Vector -> Vector -> Double
cbr v1 v2 =
  let sum1 = DV.sum . (DV.map abs) $ v1
      sum2 = DV.sum . (DV.map abs) $ v2
  in (sad v1 v2) / (sum1 + sum2)

max_stars_repo_name: emmettng/eris

max_issues_repo_name: emmettng/eris

max_forks_repo_name: emmettng/eris


text:
{-# language GADTs #-}
{-# language ConstraintKinds #-}
{-# language FlexibleContexts #-}

{-# language ScopedTypeVariables #-}
{-# LANGUAGE QuasiQuotes #-}
{-# language TypeApplications #-}

module FFT where

import Feldspar
import Feldspar.Software hiding (Arr)
import Feldspar.Software.Verify
import Feldspar.Software.Compile

import Feldspar.Array.Vector
import Feldspar.Array.Buffered

import Data.Bits (Bits)
import Data.Complex (Complex)

import Data.Selection
import Data.Default.Class

import Control.Monad

-- language-c-quote
import Language.C.Quote.GCC
import qualified Language.C.Syntax as C

-- imperative-edsl
import qualified Language.Embedded.Backend.C  as Imp

import Prelude hiding ((==), (/=), (>=), length)

--------------------------------------------------------------------------------
-- * FFT
--------------------------------------------------------------------------------

type SRef    a = Reference Software a
type SArr    a = Array  Software a
type SIArr   a = IArray Software a
type SStore  a = Store  Software a
type SSyntax a = Syntax SExp a

-- | Collection of constraints for immutable arrays.
type Immutable arr a = (
    Manifestable Software arr a
  , Finite   SExp arr
  , Indexed  SExp arr, ArrElem arr ~ a
  , Slicable SExp arr
  )

--------------------------------------------------------------------------------
-- ** Helper functions.

-- | Creates a value of type @a@ where the lowest @n@ bits are ones.
oneBits :: (Bits a, Num (SExp a), SType' a)
  => SExp Int32 -> SExp a
oneBits n = complement (ones .<<. n)

-- | Check if the @i@'th bit of @a@ is a one.
testBit :: (Integral a, Bits a, Num (SExp a), SType' a) =>
  SExp a -> SExp Index -> SExp Bool
testBit a i = i2b (a .&. (1 .<<. i2n i))

-- | Flip the @i@'th bit of @a@.
flipBit :: (Bits a, Num (SExp a), SType' a) =>
  SExp a -> SExp Index -> SExp a
flipBit a i = a `xor` (1 .<<. (i2n i))
--  SExp Index -> SExp a -> SExp a
--flipBit i a = a `xor` (1 .<<. (i2n i))

-- |
zeroBit :: (Bits a, Num (SExp a), SType' a) =>
  SExp Index -> SExp a -> SExp a
zeroBit i a = a + (a .&. (ones .<<. (i2n i)))

-- | Zeroes all but the lowest @i@ bits in @a@.
leastBits :: (Bits a, Num (SExp a), SType' a) =>
  SExp Int32 -> SExp a -> SExp a
leastBits i a = a .&. oneBits i

-- | Two to the power of @n@, i.e. @2^n@.
pow2 :: (Num (SExp a), Bits a, SType' a)
  => SExp Index -> SExp a
pow2  n = 1 .<<. i2n n
twoTo n = pow2 n

--------------------------------------------------------------------------------
-- ** Riffle network.

-- | Riffle indices.
rotBit :: SExp Index -> SExp Index -> SExp Index
rotBit k i = lefts .|. rights
  where
    k'     = i2n k
    ir     = i .>>. 1
    rights = ir .&. oneBits k'
    lefts  = (((ir .>>. k') .<<. 1) .|. (i .&. 1)) .<<. k'

-- | Riffle network for a pully array.
riffle :: (Immutable arr a, SSyntax a)
  => SExp Index -> arr -> SPull a
riffle k arr = Pull (length arr) $ \i -> arr ! rotBit k i

-- |
revBit :: (Immutable arr a, SSyntax a)
  => SStore a -> SExp Length -> arr -> Software (SManifest a)
revBit st n vec = loopStore st 1 1 (n-1) (step) vec
  where
    step i arr = return $ riffle i arr
    --step :: (Immutable arr a, SSyntax a) => SExp Index -> arr -> Software (SPush a)
--    step i arr = return $ unroll 2 $ riffle i arr
--    step k arr = return $ pairwise @Software (\i -> (i, rotBit k i)) $ riffle k arr

--------------------------------------------------------------------------------
-- ** Twiddle factors.

-- |
tw :: (Floating a, SType' a, SType' (Complex a))
  => SExp Index -> SExp Index -> SExp (Complex a)
tw n k = polar 1 (-2 * pi * i2n k / i2n n)

-- |
twids
  :: ( Immutable ts  (SExp (Complex a))
     , Immutable vec (SExp (Complex a))
     , SType' a
     , SType' (Complex a)
     , RealFloat a
     )
  => ts
  -> SExp Index
  -> SExp Index
  -> SExp Length
  -> vec
  -> SPull (SExp (Complex a))
twids ts n k l vec = Pull l $ \i ->
  let
    j = (leastBits (i2n k) i) .<<. (n'-1-k')
  in
    (testBit i k) ? ((ts!j) * (vec!i)) $ (vec!i)
  where
    n' = i2n n
    k' = i2n k

--------------------------------------------------------------------------------
-- ** Butterfly.

-- | Butterfly network.
bfly
  :: ( Immutable vec (SExp (Complex a))
     , RealFloat a
     , SType' a
     , SType' (Complex a)
     )
  => SExp Index -> vec -> SPull (SExp (Complex a))
bfly k as = Pull (length as) $ \i ->
    let a = as ! i
        b = as ! flipBit i k
    in  (testBit i k) ? (b-a) $ (a+b)

--------------------------------------------------------------------------------
-- ** FFT Core.

-- | Core of the FFT
fftCore
  :: ( Immutable ts  (SExp (Complex a))
     , Immutable vec (SExp (Complex a))
     , RealFloat a
     , SType' a
     , SType' (Complex a)
     )
  => SStore (SExp (Complex a))
  -> ts
  -> SExp Length
  -> vec
  -> Software (SManifest (SExp (Complex a)))
fftCore st ts n vec =
  let
    step k = return . twids ts n k (length vec) . bfly k
--    step k = return . pairwise @Software (\i -> (i, flipBit i k)) . twids ts n k (length vec) . bfly k
  in
    do arr <- loopStore st ((i2n n :: SExp Int32)-1) (-1) 0 (step . i2n) vec
--       revBit st n arr
       return arr

-- | Radix-2 Decimation-In-Frequency Fast Fourier Transformation of the given
-- complex vector. The given vector must be power-of-two sized, (for example 2,
-- 4, 8, 16, 32, etc.) The output is non-normalized.
fft
  :: ( Immutable vec (SExp (Complex a))
     , RealFloat a
     , SType' a
     , SType' (Complex a)
     )
  => SExp Length
  -> SStore (SExp (Complex a))
  -> vec
  -> Software (SManifest (SExp (Complex a)))
fft n st vec =
  do n  <- shareM (ilog2 (length vec))
     ts <- manifestFresh $ Pull (pow2 (n-1)) (tw (pow2 n))
     fftCore st ts n vec

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

fftCore2 ::
  SStore (SExp (Complex Double)) ->
  SIArr (SExp (Complex Double)) ->
  SExp Length ->
  SPull (SExp (Complex Double)) ->
  Software (SIArr (SExp (Complex Double)))
fftCore2 store iarr exp pull = do
  let iarrm :: SManifest (SExp (Complex Double))
      iarrm = M iarr
  (M ibrr) <- fftCore store iarrm exp pull
  return ibrr

fft2 ::
  SExp Length ->
  SStore (SExp (Complex Double)) ->
  SPull (SExp (Complex Double)) ->
  Software (SIArr (SExp (Complex Double)))
fft2 n store pull = do
  (M ts) <- manifestFresh $ Pull (twoTo (n-1)) (tw (twoTo n))
  fftCore2 store ts n pull

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

example :: Software ()
example = do
  size :: SExp Length <- fget stdin
  n :: SExp Length <- fget stdin
  assert (1 `shiftL` n == size) "2^n == size"
  assert (size >= 2) "not too small"
  st   :: SStore (SExp (Complex Double)) <- newInPlaceStore size
  arr  :: SArr   (SExp (Complex Double)) <- newArr size
  iarr :: SIArr  (SExp (Complex Double)) <- freezeArr arr
  fft n st iarr
  assert (value False) "oh no"
  return ()

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

bad :: Software ()
bad = do
  size :: SExp Length <- fget stdin
  fprintf stdout "%d" (share size id)
  return ()

--------------------------------------------------------------------------------
-- FFT Bench. Copy of https://github.com/Feldspar/raw-feldspar/blob/master/examples/FFT_bench.hs
--------------------------------------------------------------------------------

printTime_def = [cedecl|
void printTime(typename clock_t start, typename clock_t end)
{
  printf("CPU time (sec): %f\n", (double)(end-start) / CLOCKS_PER_SEC);
}
|]

sizeOf_double_complex :: SExp Length
sizeOf_double_complex = 16

-- | Measure the time for 100 runs of 'fftCore' (excluding initialization) for
-- arrays of the given size
benchmark :: SExp Length -> Software ()
benchmark n = do
  addInclude "<stdio.h>"
  addInclude "<string.h>"
  addInclude "<time.h>"

  addDefinition printTime_def

  start <- newObject "clock_t" False
  end   <- newObject "clock_t" False

  st :: SStore (SExp (Complex Double)) <- newStore n
  inp <- unsafeFreezeStore n st
  callProc "memset"
      [ iarrArg (manifest inp)
      , valArg (0 :: SExp Index)
      , valArg (n*sizeOf_double_complex)
      ]

  n  <- shareM (ilog2 (length inp))
  ts <- manifestFresh $ Pull (twoTo (n-1)) (tw (twoTo n))
  -- Change `manifestFresh` to `return` to avoid pre-computing twiddle factors

  callProcAssign start "clock" []

  for 0 1 99 $ \(_ :: SExp Index) ->
    void $ fftCore st ts n inp

  callProcAssign end "clock" []
  callProc "printTime" [objArg start, objArg end]

runBenchmark n = runCompiled'
    (def :: CompilerOpts) --{compilerAssertions = select []}
    -- Note: important to turn off assertions when running the benchmarks
    --       (in old Raw-Feldspar that is, Co-Feldspar gets rid of them)
    def {Imp.externalFlagsPre = ["-O3"], Imp.externalFlagsPost = ["-lm"]}
    (benchmark n)

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

dummy :: Software ()
dummy = do
  size :: SExp Length <- fget stdin
  n :: SExp Length <- fget stdin
  assert (1 `shiftL` n == size) "2^n == size"
  assert (size >= 2) "not too small"
  st   :: SStore (SExp (Complex Double)) <- newInPlaceStore size
  arr  :: SArr   (SExp (Complex Double)) <- newArr size
  iarr :: SIArr  (SExp (Complex Double)) <- freezeArr arr
  fft n st iarr
  fft n st iarr
  assert (value False) "oh no"
  return ()

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

max_stars_repo_name: markus-git/co-feldspar

max_issues_repo_name: markus-git/co-feldspar

max_forks_repo_name: markus-git/co-feldspar


text:
module Main where

import Lib
import Types
import Numeric.Vector
import Wavefront.Wavefront

cam :: Camera
cam = Camera { c_getPos = vec3 3 2 (-5)
             , c_getForward = vec3 0 0 (1)
             , c_getUp = vec3 0 1 0
             , c_getRight = vec3 1 0 0
             , c_getFOV = 90
             }

checkerShader :: BRDF
checkerShader = let spacing = VecConst 1 1 1
                    quarterSpacing = VecMath VDiv spacing (VecConst 4 4 4)
                    cs = VecMath
                           VSub
                           (VecConst 1 1 1)
                           (VecMath
                             VMul
                             (VecConst 2 2 2)
                             (VecMath
                               VDiv
                               (VecMath
                                 VMod
                                 (VecMath
                                   VAdd
                                   quarterSpacing
                                   (VecMath VAbs (VecMath VAdd UV quarterSpacing) (VecConst 0 0 0))
                                 )
                                 spacing
                               )
                               spacing
                             )
                           )
                    lt = ValMath LessThan (ValMath Mul (SeparateX cs) (SeparateY cs)) (ValConst 0)
                    value = ValMath Add (ValConst 0.7) (ValMath Mul (ValConst 0.3) lt)
                    color = CombineRGB value value value
                in Diffuse color

shader :: BRDF
shader = Diffuse $ ColConst 0.8 0.8 0.8

sun :: Light
sun = Directional { l_getAngle = 2.5
                  , l_getStrength = 1
                  , l_getDirection = normalized $ vec3 0.5 (-1) (-0.4)
                  , l_getColor = vec3 1 1 1
                  }

scene :: Scene
scene = Scene { s_getCamera = cam
              , s_getObjects = []
              , s_getLights = [sun]
              , s_getSkyColor = vec3 0.6 0.7 0.9
              , s_getWidth = 640
              , s_getHeight = 480
              , s_getSeed = 1
              , s_getSamples = 128
              , s_getBounces = 2
              }

main :: IO ()
main = do
  maybeCube <- readObj "cube.obj" <$> readFile "/home/craig/haskell/raytracer/app/cube.obj"
  case maybeCube of
    Left e -> print e
    Right cube -> do
      -- print cube
      let scene' = scene { s_getObjects = [Object cube shader] }
      writeImg "testimg.png" $ parRender scene'


text:
{-# LANGUAGE TypeFamilies #-}
module Data.Running
  ( Coord (Coord)
  , PointRecord (PointRecord)
  , pointDistance
  , pointCoord
  , pointLat
  , pointLng
  , pointTime
  , pointSpeed
) where
import Data.Time.Clock (UTCTime)
import Data.Time.Format ()
import Statistics.Sample (mean, stdDev)
import Data.AdditiveGroup
import Data.VectorSpace

newtype Coord = Coord (Double, Double) deriving (Eq, Ord, Show)

data PointRecord = PointRecord
  { pointTime :: UTCTime
  , pointCoord :: Coord
  , pointDistance :: Double
  , pointSpeed :: Double
  } deriving (Eq, Ord, Show)

pointLat (PointRecord _ (Coord (a, _)) _ _) = a
pointLng (PointRecord _ (Coord (_, b)) _ _) = b

instance AdditiveGroup Coord where
  zeroV = Coord (0, 0)
  Coord (a, b) ^+^ Coord (c, d) = Coord (a + c, b + d)
  negateV (Coord (a, b)) = Coord (-a, -b)

instance InnerSpace Coord where
  Coord (a, b) <.> Coord (c, d) = a * c + b * d

instance VectorSpace Coord where
  type Scalar Coord = Double
  s *^ Coord (a, b) = Coord (s * a, s * b)

--instance VectorSpace PointRecord where
--  type Scalar =


text:
{-# LANGUAGE FlexibleContexts #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Statistics.ICA
-- Copyright   :  (c) A. V. H. McPhail 2010, 2014
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Independent Components Analysis
--
--  implements the FastICA algorithm found in:
--
-- * Aapo Hyvärinen and Erkki Oja,
--   Independent Component Analysis: Algorithms and Applications,
--   /Neural Networks/, 13(4-5):411-430, 2000
--
--   <http://www.google.com/url?sa=t&source=web&cd=2&ved=0CBgQFjAB&url=http%3A%2F%2Fciteseerx.ist.psu.edu%2Fviewdoc%2Fdownload%3Fdoi%3D10.1.1.79.7003%26rep%3Drep1%26type%3Dpdf&ei=RQozTJb6L4_fcbCV6cMD&usg=AFQjCNGClLIB9MAvbrEj45SyUx9cYubLyA&sig2=hg5Wnfy3dLPkoIc1hqSfjg>
--
-----------------------------------------------------------------------------

module Numeric.Statistics.ICA (
                               sigmoid, sigmoid',
                               demean, whiten,
                               ica, icaDefaults,
                               --
                               NormType(..)
                          ) where

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

import qualified Data.Array.IArray as I

import Prelude hiding ((<>))

import Numeric.LinearAlgebra

import qualified Data.Vector.Generic as GV

import Numeric.GSL.Statistics

import Numeric.Statistics

import System.Random

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

data NormType = NormZero | NormOne | NormTwo | NormInf

pnorm :: Normed (Vector a) => NormType -> Vector a -> R
pnorm NormZero = norm_0
pnorm NormOne  = norm_1
pnorm NormTwo  = norm_2
pnorm NormInf  = norm_Inf

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

-- | sigmoid transfer function
sigmoid :: Double -> Double
sigmoid u = u * exp((-u**2)/2)

-- | derivative of sigmoid transfer function
sigmoid' :: Double -> Double
sigmoid' u = -u**2 * exp((-u**2)/2)

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

-- preprocessing:
--   demean
--   whiten
--        eigenvalue decomposition of covariance matrix E{xx^T} = EDE^T
--                   E orthogonal matrix of eigenvectors
--                   D diagonal matrix of eigenvalues, D = diag(d_1,...,d_n)
--                   x_white = ED^{-1/2}E^Tx
--                   D^{-1/2} = diag{d_1^{-1/2},...}
--

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

-- | remove the mean from data
demean :: I.Array Int (Vector Double)                  -- ^ the data
       -> (I.Array Int (Vector Double),Vector Double)  -- ^ (demeaned data,mean)
demean d = let u = I.elems $ fmap mean d
               d' = I.listArray (I.bounds d) (zipWith (-) (I.elems d) (map scalar u))
               u' = fromList u
               in (d',u')

-- | whiten data
whiten :: I.Array Int (Vector Double)                 -- ^ the data
       -> Double                                      -- ^ eigenvalue threshold
       -> (I.Array Int (Vector Double),Matrix Double) -- ^ (whitened data,transform)
whiten d q = let cv = covarianceMatrix d
                 (val',vec') = eigSH $ trustSym cv  -- the covariance matrix is real symmetric
                 val = toList val'
                 vec = toColumns vec'
                 v' = zip val vec
                 v = filter ((> q) . fst) v'        -- keep only eigens > than parameter
                 (dd',e') = unzip v
                 dd = diag $ (** (-0.5)) $ fromList dd'  -- square root of eigenvalues diagonalised
                 e = fromColumns e'
                 x = fromRows $ I.elems d
                 t = e <> dd <> tr' e          -- the actual mathematics
                 x' = t <> x                     -- the actual mathematics
                 d' = I.listArray (I.bounds d) (toRows x')
             in (d',t)

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

-- assuming that a weight vector is a row

-- algorithm:
-- 1  initial random weight vectors w_i
-- 2  w_i^+ = E{xg(w^Tx)} - E{g'(w^Tx)}w     (newton phase)
-- 3  W = W = (WW^T)^{-1/2)W                 (decorrelation) W = ( ..., w_i, ...)^T
--                                           WW^T = FDF^T (eigenvalue decomposition)
-- 4  w_i = w^+/norm(w^+)                    (normalisation) (almost any norm but not Frobenius)
-- 5  if not converged (dot w w^+ ~ 1 implies convergence) go to step 2
--
-- in matrix form, 2 becomes:
--    W^+ = W + (diag a_i)[(diag b_i) + E{g(y)y^T}]W
--
--      where
--            y = Wx
--            b_i = -E{y_ig(y_i)}
--            a_i = -1/(b_i-E{g'(y_i)})
--
--    g(u) = tanh(au) 0<=a<=2, often a = 1
--    g(u) = u exp(-u^2/2)

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

unconcat 0 _ _  = []
unconcat r c xs = [take c xs] ++ unconcat (r-1) c (drop c xs)

random_vector :: Int -> (Int,Int) -> Matrix Double
random_vector s (r,c) = fromLists $ unconcat r c $ randomRs (-1,1) (mkStdGen s)

-- g g' w x -> w'
update :: (Double -> Double) -> (Double -> Double) -> Matrix Double -> Matrix Double -> Matrix Double
update g g' w x = let y = w <> x
                      ys = toRows y
                      bis = map (\y' -> - mean (y' * (GV.map g y'))) ys
                      ais = zipWith (\b y' -> -1 / (b - mean (GV.map g y'))) bis ys
                      r = rows y
                      ix = ((1,1),(r,r))
                      cov = fromArray2D $ I.listArray ix $ map (\(m,n) -> covariance (GV.map g' (ys!!(m-1))) (ys!!(n-1))) $ I.range ix
                  in w + (diag $ fromList ais) <> ((diag $ fromList bis) + cov) <> w

decorrelate :: Matrix Double -> Matrix Double
decorrelate m = let (d',v') = eig m
                    d = fst $ fromComplex d'
                    v = fst $ fromComplex v'
                in v <> (diag (d ** (-0.5))) <> tr' v <> m
{-decorrelate n t w = let w' = w / (scalar $ sqrt $ pnorm n (w <> tr' w))
                    in decorrelate' t w w'
    where decorrelate' t' m m'
              | converged t' m m' = m'
              | otherwise         = decorrelate' t' m' ((scale 1.5 m') - (scale 0.5 (m' <> tr' m' <> m')))
-}

normalise :: NormType -> Matrix Double -> Matrix Double
normalise t m = fromRows $ map (\v -> v / (scalar $ pnorm t v)) (toRows m)

converged :: Double -> Matrix Double -> Matrix Double -> Bool
converged t m m' = let d' = map ((-) 1) $ zipWith dot (toRows m) (toRows m')
                   in maximum d' <= t

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

ica' :: (Double -> Double)          -- ^ transfer function (tanh,u exp(u^2/2), etc...)
     -> (Double -> Double)          -- ^ derivative of transfer function
     -> NormType                    -- ^ type of normalisation: Infinity, PNorm1, PNorm2
     -> Double                      -- ^ convergence tolerance for feature vectors
     -> Matrix Double               -- ^ weight matrix
     -> [Matrix Double]             -- ^ input data in chunks
     -> Matrix Double               -- ^ ica transform (weight matrix)
ica' _ _  _ _ _ []     = error "no sample data"
ica' g g' n t w (x:xs) = let w' = normalise n $ decorrelate $ update g g' w x
                             in if converged t w w'
                                then w'
                                else ica' g g' n t w' (xs ++ [x])

-- | perform an ICA transform
ica :: Int                         -- ^ random seed
    -> (Double -> Double)          -- ^ transfer function (tanh,u exp(u^2/2), etc...)
    -> (Double -> Double)          -- ^ derivative of transfer function
    -> NormType                    -- ^ type of normalisation: Infinity, PNorm1, PNorm2
    -> Double                      -- ^ convergence tolerance for feature vectors
--    -> Int                         -- ^ output dimensions
    -> Int                         -- ^ sampling size (must be smaller than length of data)
    -> I.Array Int (Vector Double) -- ^ data
    -> (I.Array Int (Vector Double),Matrix Double) -- ^ transformed data, ica transform
ica r g g' n t s a = let j = I.rangeSize $ I.bounds a
                         w = random_vector r (j,j)
                         x' = fromRows $ I.elems a
                         -- next line is BAD if distribution not stationary
                         x = concat $ toBlocksEvery j s x'
                         w' = ica' g g' n t w x
                         y = w' <> x'
                     in (I.listArray (1,1) $ toRows y,w')

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

-- | ICA with default values: no dimension reduction, euclidean norms, 16 sample groups, sigmoid
icaDefaults :: Int                         -- ^ random seed
            -> I.Array Int (Vector Double) -- ^ data
            -> (I.Array Int (Vector Double),Matrix Double) -- ^ transformed data, ica transform
icaDefaults r a = let c = I.rangeSize $ I.bounds a
                      s = (GV.length $ (a I.! 1)) `div` 16
                  in ica r sigmoid sigmoid' NormInf 0.0000001 s a

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

max_stars_repo_name: amcphail/hstatistics

max_issues_repo_name: amcphail/hstatistics

max_forks_repo_name: amcphail/hstatistics


text:
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeOperators       #-}
{-# LANGUAGE TypeApplications    #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
module Test.Grenade.Layers.LRN where

import           Data.Singletons              ()
import           GHC.TypeLits
import           Hedgehog
import           Data.Proxy
import           Data.Kind (Type)
import           Test.Hedgehog.Hmatrix
import           Test.Grenade.Layers.Internal.Reference
import qualified System.Random.MWC as MWC
import           Control.DeepSeq
import           Control.Exception (evaluate)
import           Test.Hedgehog.Compat
import           Data.Serialize
import           Data.Either

import           Grenade.Core
import           Grenade.Layers.LRN
import           Grenade.Types
import           Numeric.LinearAlgebra.Data   as NLD hiding ((===))
import           Numeric.LinearAlgebra.Static as H hiding ((===))

expectedDepth1
  = [
      [0.99992500, 1.99940020, 2.99797659],
      [3.99520671, 4.99064546, 5.98385086],
      [6.97438480, 7.96181378, 8.94570965],
      [1.09990018, 2.09930569, 3.09776755],
      [4.09483851, 5.09007376, 6.08303166],
      [7.07327453, 8.06036937, 9.04388861],
      [1.19987041, 2.19920173, 3.19754459],
      [4.19445196, 5.18947928, 6.18218531],
      [7.17213277, 8.15888920, 9.14202759]
    ]

expectedDepth2
  = [
      [0.99992500, 1.99940020, 2.99797659],
      [3.99520671, 4.99064546, 5.98385086],
      [6.97438480, 7.96181378, 8.94570965],

      [1.09981771, 2.09867639, 3.09568020],
      [4.08993980, 5.08057535, 6.06671955],
      [7.04752047, 8.02214440, 8.98977845],

      [1.19976155, 2.19847498, 3.19524425],
      [4.18918085, 5.17940607, 6.16505399],
      [7.14527441, 8.11923556, 9.08612680]
    ]

expectedDepth3
  = [
      [0.99983428, 1.99873942, 2.99581955],
      [3.99018517, 4.98095623, 5.96726513],
      [6.94825962, 7.92310558, 8.89098967],

      [1.09969897, 2.09791554, 3.09330926],
      [4.08455479, 5.07034831, 6.04941408],
      [7.02051070, 7.98243693, 8.93403711],

      [1.19976155, 2.19847498, 3.19524425],
      [4.18918085, 5.17940607, 6.16505399],
      [7.14527441, 8.11923556, 9.08612680]
    ]

prop_lrn_forwards_fixed :: Property
prop_lrn_forwards_fixed = property $ do
  let ch1 = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
      ch2 = map (map (0.1 +)) ch1
      ch3 = map (map (0.1 +)) ch2
      dat = concat [ch1, ch2, ch3]
      mat = (NLD.fromLists dat) :: NLD.Matrix RealNum

      lrn   = LRN :: (LRN "0.0001" "0.75" "1" 1)
      lrn'  = LRN :: (LRN "0.0001" "0.75" "1" 2)
      lrn'' = LRN :: (LRN "0.0001" "0.75" "1" 3)

      Just inp = (H.create mat) :: Maybe (H.L 9 3)
      (out,   expMat)   = getOutput lrn   expectedDepth1 inp
      (out',  expMat')  = getOutput lrn'  expectedDepth2 inp
      (out'', expMat'') = getOutput lrn'' expectedDepth3 inp

  out   `isSimilarMatrixTo` expMat
  out'  `isSimilarMatrixTo` expMat'
  out'' `isSimilarMatrixTo` expMat''
  where
    -- Helper function to
    getOutput :: (KnownSymbol a, KnownSymbol b, KnownSymbol k, KnownNat n) => LRN a b k n -> [[RealNum]] -> H.L 9 3 -> (Matrix RealNum, Matrix RealNum)
    getOutput l e inp = (H.extract o, NLD.fromLists e)
      where
      (_, S3D o :: (S ('D3 3 3 3))) = runForwards l ((S3D inp) :: (S ('D3 3 3 3)))

data OpaqueLRN :: Type where
  OpaqueLRN :: (KnownNat i) => LRN "0.0001" "0.75" "1" i -> OpaqueLRN

genOpaqueLRN :: Gen OpaqueLRN
genOpaqueLRN = do
  s :: Integer <- choose 1 7
  let Just s' = someNatVal s
  case s' of
      SomeNat (Proxy :: Proxy i') ->
          return . OpaqueLRN $ (LRN :: LRN "0.0001" "0.75" "1" i')

prop_lrn_forwards :: Property
prop_lrn_forwards = property $ do
  OpaqueLRN (lrn :: LRN "0.0001" "0.75" "1" n) <- blindForAll genOpaqueLRN
  let n' = natVal (Proxy :: Proxy n)
  source <- forAll $ genLists3D 5 11 7 -- 5 channels, 11 rows, 7 columns
  let inp = S3D (H.fromList $ (concat . concat) source) :: S ('D3 11 7 5)
      (_, o :: (S ('D3 11 7 5))) = runForwards lrn inp
      out = naiveLRNForwards 0.0001 0.75 1 (fromIntegral n') source

  assert $ allClose o (S3D (H.fromList $ (concat . concat) out))

prop_lrn_backwards :: Property
prop_lrn_backwards = property $ do
  OpaqueLRN (lrn :: LRN "0.0001" "0.75" "1" n) <- blindForAll genOpaqueLRN
  let n' = natVal (Proxy :: Proxy n)
  source <- forAll $ genLists3D 5 11 7 -- 5 channels, 11 rows, 7 columns
  let inp = S3D (H.fromList $ (concat . concat) source) :: S ('D3 11 7 5)
      (tape, o :: (S ('D3 11 7 5))) = runForwards lrn inp
      ((),    d :: (S ('D3 11 7 5))) = runBackwards lrn tape o
      out = naiveLRNForwards 0.0001 0.75 1 (fromIntegral n') source
      back = naiveLRNBackwards 0.0001 0.75 1 (fromIntegral n') source out

  assert $ allClose d (S3D (H.fromList $ (concat . concat) back))

prop_lrn_show :: Property
prop_lrn_show = withTests 1 $ property $ do
  let lrn = LRN :: (LRN "0.0001" "0.75" "1" 1)
  (show lrn) `seq` success

prop_lrn_rnf :: Property
prop_lrn_rnf = withTests 1 $ property $ do
  let lrn = LRN :: (LRN "0.0001" "0.75" "1" 1)
  (r :: ()) <- evalIO $ evaluate $ rnf lrn
  r `seq` success

prop_lrn_run_update :: Property
prop_lrn_run_update = withTests 1 $ property $ do
  source <- forAll $ genLists 12 17
  gen <- evalIO MWC.create
  (lrn :: LRN "0.0001" "0.75" "1" 1) <- evalIO $ createRandomWith UniformInit gen
  let mat = (NLD.fromLists source) :: NLD.Matrix RealNum
      Just inp = (H.create mat) :: Maybe (H.L 12 17)
      (tape, o :: (S ('D3 3 17 4))) = runForwards lrn ((S3D inp) :: (S ('D3 3 17 4)))
      (grad, _ :: (S ('D3 3 17 4))) = runBackwards lrn tape o
  (runUpdate defSGD lrn grad) `seq` (runUpdate defAdam lrn grad) `seq` success
  (reduceGradient @(LRN "0.0001" "0.75" "1" 1) [grad]) `seq` success

prop_lrn_serializable :: Property
prop_lrn_serializable = withTests 1 $ property $ do
  gen <- evalIO MWC.create
  (lrn :: LRN "0.0001" "0.75" "1" 1) <- evalIO $ createRandomWith UniformInit gen
  let bs  = encode lrn
      dec = decode bs :: Either String (LRN "0.0001" "0.75" "1" 1)
  assert $ isRight dec
  (put lrn) `seq` success

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

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE FlexibleContexts #-}

module SVM
  ( multiplyByItself
  , testmatrix
  , testvector
  , replicateVector
  , replicateVectorSquared
  , makeLabelAgreementMatrix
  , makeQuadraticMatrix
  , testmatrixsize
  , alphasQuadProg
  , trainPLAandSVM
  , countSV
  ) where

import           LinearRegression                (createMatrixX,
                                                  createRandomPoints,
                                                  createVectorY,
                                                  makeTargetFunction, y)
import           Numeric.LinearAlgebra
import           Numeric.LinearAlgebra.Data
import           Numeric.LinearAlgebra.HMatrix   as HMatrix
import           Numeric.Minimization.QuadProgPP
import           PLA

----------------------------------------------------
-- Solutions to Homework 7 on SVM
----------------------------------------------------
-- | a test matrix - not part of homework
testmatrix :: Matrix R
testmatrix = matrix 2 [1, 1, 1, -1, -1, -1, 3, 4, 5, 6]

testmatrixsize = size testmatrix

-- | a test vector - not part of homework
testvector :: Vector R
testvector = vector [-1, 1, 1, -1, 1]

-- | create a matrix where every entry x_n_m is the dot product of the nth vector with the mth vector
multiplyByItself :: Matrix R -> Matrix R
multiplyByItself matrixA = matrixA HMatrix.<> (tr matrixA)

-- | make a matrix where every entry corresponds to the product of the labels of the corresponding data - it will be 1 if they agree and -1 otherwise
makeLabelAgreementMatrix :: Vector R -> Matrix R
makeLabelAgreementMatrix labelvector =
  diag labelvector HMatrix.<> replicateVectorSquared labelvector

-- | make a matrix form a vector by repeating the vector n times and arranging into n rows, each corresponding to the vector
replicateVector :: Int -> Vector R -> Matrix R
replicateVector n v = fromRows $ replicate n v

-- | make a squared matrix from a vector, such that each rows corresponds to the vector
replicateVectorSquared :: Vector R -> Matrix R
replicateVectorSquared v = replicateVector (size v) v

-- | makes the matrix to pass to quadratix programming
makeQuadraticMatrix :: Matrix R -> Vector R -> Matrix R
makeQuadraticMatrix dataMatrix labelVector =
  multiplyByItself dataMatrix * makeLabelAgreementMatrix labelVector

makeQuadraticMatrix' :: Matrix R -> Vector R -> Matrix R
makeQuadraticMatrix' dataMatrix labelVector = quadMatrix + smallValtoDiag
  where
    quadMatrix = makeQuadraticMatrix dataMatrix labelVector
    smallValtoDiag = diag $ konst 0.00000000001 (size labelVector)

-- | solves quadratic Programming on the input data matrix and label Vector
alphasQuadProg ::
     Matrix R -> Vector R -> Either QuadProgPPError (Vector R, Double)
alphasQuadProg dataMatrix labelVector =
  solveQuadProg
    (matrixA, vectorB)
    (Just (labelVecMatrix, zeros))
    (Just (idMatrix, zeros))
  where
    matrixA = makeQuadraticMatrix' dataMatrix labelVector
    vectorB = konst (-1) (size labelVector)
    labelVecMatrix = diag labelVector
    idMatrix = ident (size labelVector)
    zeros = konst 0 (size labelVector)

-- | Creates the data matrix from a list of data points
createDataMatrix :: [(R, R)] -> Matrix R
createDataMatrix listOfPoints = matrix 2 listOfNumbers
  where
    listOfNumbers = concat [[a, b] | (a, b) <- listOfPoints]

setSmallTo0 x
  | x < (10 ** (-4)) = 0
  | otherwise = x

countSV vector = length . filter (/= 0) $ toList vector

-- | trains a perceptron and SVM
trainPLAandSVM :: Int -> IO ()
trainPLAandSVM n
    --making target function, training and testing points
    -----------------------------------------------------
 = do
  target <- makeTargetFunction
  trainpoints <- createRandomPoints n
  testpoints <- createRandomPoints 1000
  let trainX = createDataMatrix trainpoints
  let trainY = createVectorY (y target) trainpoints
  let testX = createDataMatrix testpoints
  let testY = createVectorY (y target) testpoints
    --Perceptron
    ---------------------------------
  let initialWeights = (0, 0, 0)
  (finalWeights, epochs) <-
    trainPLAWithInitial initialWeights trainpoints target
  let misclassifiedTestPoints =
        length $ filter (isMisclassified target finalWeights) testpoints
  let prob = fromIntegral misclassifiedTestPoints / 1000
    --SVM
    --------------------------
  let quadmatrix = makeQuadraticMatrix trainX trainY
  let labelVecMatrix = asColumn trainY
  let alphas = alphasQuadProg trainX trainY
  let alphasModified = fmap (cmap setSmallTo0 . fst) alphas
  let numSVs = fmap countSV alphasModified
    --Printing the output
    ---------------------------
  putStrLn "Perceptron:"
  putStrLn "-------------------------"
  putStr "Weights: "
  print finalWeights
  putStr "Epochs: "
  print epochs
  putStr "Probability f /= g: "
  print prob
  putStrLn "Support Vector Machine:"
  putStrLn "-------------------------"
  putStr "Alphas:"
  print alphas
  putStr "Alphas very small set ot 0:"
  print alphasModified
  putStr "number of SVs:"
  print numSVs