text:
-- Copyright 2016 TensorFlow authors.
--
-- Licensed under the Apache License, Version 2.0 (the "License");
-- you may not use this file except in compliance with the License.
-- You may obtain a copy of the License at
--
--     http://www.apache.org/licenses/LICENSE-2.0
--
-- Unless required by applicable law or agreed to in writing, software
-- distributed under the License is distributed on an "AS IS" BASIS,
-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-- See the License for the specific language governing permissions and
-- limitations under the License.

{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE MonoLocalBinds #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
-- We use UndecidableInstances for type families with recursive definitions
-- like "\\".  Those instances will terminate since each equation unwraps one
-- cons cell of a type-level list.
{-# LANGUAGE UndecidableInstances #-}

module TensorFlow.Types
    ( TensorType(..)
    , TensorData(..)
    , TensorDataType(..)
    , Scalar(..)
    , Shape(..)
    , protoShape
    , Attribute(..)
    , DataType(..)
    , ResourceHandle
    , Variant
    -- * Lists
    , ListOf(..)
    , List
    , (/:/)
    , TensorTypeProxy(..)
    , TensorTypes(..)
    , TensorTypeList
    , fromTensorTypeList
    , fromTensorTypes
    -- * Type constraints
    , OneOf
    , type (/=)
    , OneOfs
    -- ** Implementation of constraints
    , TypeError
    , ExcludedCase
    , NoneOf
    , type (\\)
    , Delete
    , AllTensorTypes
    ) where

import Data.Functor.Identity (Identity(..))
import Data.Complex (Complex)
import Data.Default (def)
import Data.Int (Int8, Int16, Int32, Int64)
import Data.Maybe (fromMaybe)
import Data.Monoid ((<>))
import Data.ProtoLens.TextFormat (showMessageShort)
import Data.Proxy (Proxy(..))
import Data.String (IsString)
import Data.Word (Word8, Word16, Word32, Word64)
import Foreign.Storable (Storable)
import GHC.Exts (Constraint, IsList(..))
import Lens.Family2 (Lens', view, (&), (.~), (^..))
import Lens.Family2.Unchecked (iso)
import Text.Printf (printf)
import qualified Data.Attoparsec.ByteString as Atto
import Data.ByteString (ByteString)
import qualified Data.ByteString as B
import Data.ByteString.Builder (Builder)
import qualified Data.ByteString.Builder as Builder
import qualified Data.ByteString.Lazy as L
import qualified Data.Vector as V
import qualified Data.Vector.Storable as S
import Proto.Tensorflow.Core.Framework.AttrValue
    ( AttrValue(..)
    , AttrValue'ListValue(..)
    )
import Proto.Tensorflow.Core.Framework.AttrValue_Fields
    ( b
    , f
    , i
    , s
    , list
    , type'
    , shape
    , tensor
    )

import Proto.Tensorflow.Core.Framework.ResourceHandle
    (ResourceHandleProto)
import Proto.Tensorflow.Core.Framework.Tensor as Tensor
    (TensorProto(..))
import Proto.Tensorflow.Core.Framework.Tensor_Fields as Tensor
    ( boolVal
    , doubleVal
    , floatVal
    , intVal
    , int64Val
    , resourceHandleVal
    , stringVal
    , uint32Val
    , uint64Val
    )

import Proto.Tensorflow.Core.Framework.TensorShape
    (TensorShapeProto(..))
import Proto.Tensorflow.Core.Framework.TensorShape_Fields
    ( dim
    , size
    , unknownRank
    )
import Proto.Tensorflow.Core.Framework.Types (DataType(..))

import TensorFlow.Internal.VarInt (getVarInt, putVarInt)
import qualified TensorFlow.Internal.FFI as FFI

type ResourceHandle = ResourceHandleProto

-- | Dynamic type.
-- TensorFlow variants aren't supported yet. This type acts a placeholder to
-- simplify op generation.
data Variant

-- | The class of scalar types supported by tensorflow.
class TensorType a where
    tensorType :: a -> DataType
    tensorRefType :: a -> DataType
    tensorVal :: Lens' TensorProto [a]

instance TensorType Float where
    tensorType _ = DT_FLOAT
    tensorRefType _ = DT_FLOAT_REF
    tensorVal = floatVal

instance TensorType Double where
    tensorType _ = DT_DOUBLE
    tensorRefType _ = DT_DOUBLE_REF
    tensorVal = doubleVal

instance TensorType Int32 where
    tensorType _ = DT_INT32
    tensorRefType _ = DT_INT32_REF
    tensorVal = intVal

instance TensorType Int64 where
    tensorType _ = DT_INT64
    tensorRefType _ = DT_INT64_REF
    tensorVal = int64Val

integral :: Integral a => Lens' [Int32] [a]
integral = iso (fmap fromIntegral) (fmap fromIntegral)

instance TensorType Word8 where
    tensorType _ = DT_UINT8
    tensorRefType _ = DT_UINT8_REF
    tensorVal = intVal . integral

instance TensorType Word16 where
    tensorType _ = DT_UINT16
    tensorRefType _ = DT_UINT16_REF
    tensorVal = intVal . integral

instance TensorType Word32 where
    tensorType _ = DT_UINT32
    tensorRefType _ = DT_UINT32_REF
    tensorVal = uint32Val

instance TensorType Word64 where
    tensorType _ = DT_UINT64
    tensorRefType _ = DT_UINT64_REF
    tensorVal = uint64Val

instance TensorType Int16 where
    tensorType _ = DT_INT16
    tensorRefType _ = DT_INT16_REF
    tensorVal = intVal . integral

instance TensorType Int8 where
    tensorType _ = DT_INT8
    tensorRefType _ = DT_INT8_REF
    tensorVal = intVal . integral

instance TensorType ByteString where
    tensorType _ = DT_STRING
    tensorRefType _ = DT_STRING_REF
    tensorVal = stringVal

instance TensorType Bool where
    tensorType _ = DT_BOOL
    tensorRefType _ = DT_BOOL_REF
    tensorVal = boolVal

instance TensorType (Complex Float) where
    tensorType _ = DT_COMPLEX64
    tensorRefType _ = DT_COMPLEX64
    tensorVal = error "TODO (Complex Float)"

instance TensorType (Complex Double) where
    tensorType _ = DT_COMPLEX128
    tensorRefType _ = DT_COMPLEX128
    tensorVal = error "TODO (Complex Double)"

instance TensorType ResourceHandle where
    tensorType _ = DT_RESOURCE
    tensorRefType _ = DT_RESOURCE_REF
    tensorVal = resourceHandleVal

instance TensorType Variant where
    tensorType _ = DT_VARIANT
    tensorRefType _ = DT_VARIANT_REF
    tensorVal = error "TODO Variant"

-- | Tensor data with the correct memory layout for tensorflow.
newtype TensorData a = TensorData { unTensorData :: FFI.TensorData }

-- | Types that can be converted to and from 'TensorData'.
--
-- 'S.Vector' is the most efficient to encode/decode for most element types.
class TensorType a => TensorDataType s a where
    -- | Decode the bytes of a 'TensorData' into an 's'.
    decodeTensorData :: TensorData a -> s a
    -- | Encode an 's' into a 'TensorData'.
    --
    -- The values should be in row major order, e.g.,
    --
    --   element 0:   index (0, ..., 0)
    --   element 1:   index (0, ..., 1)
    --   ...
    encodeTensorData :: Shape -> s a -> TensorData a

-- All types, besides ByteString and Bool, are encoded as simple arrays and we
-- can use Vector.Storable to encode/decode by type casting pointers.

-- TODO(fmayle): Assert that the data type matches the return type.
simpleDecode :: Storable a => TensorData a -> S.Vector a
simpleDecode = S.unsafeCast . FFI.tensorDataBytes . unTensorData

simpleEncode :: forall a . (TensorType a, Storable a)
             => Shape -> S.Vector a -> TensorData a
simpleEncode (Shape xs) v =
    if product xs /= fromIntegral (S.length v)
        then error $ printf
            "simpleEncode: bad vector length for shape %v: expected=%d got=%d"
            (show xs) (product xs) (S.length v)
        else TensorData (FFI.TensorData xs dt (S.unsafeCast v))
  where
    dt = tensorType (undefined :: a)

instance TensorDataType S.Vector Float where
    decodeTensorData = simpleDecode
    encodeTensorData = simpleEncode

instance TensorDataType S.Vector Double where
    decodeTensorData = simpleDecode
    encodeTensorData = simpleEncode

instance TensorDataType S.Vector Int8 where
    decodeTensorData = simpleDecode
    encodeTensorData = simpleEncode

instance TensorDataType S.Vector Int16 where
    decodeTensorData = simpleDecode
    encodeTensorData = simpleEncode

instance TensorDataType S.Vector Int32 where
    decodeTensorData = simpleDecode
    encodeTensorData = simpleEncode

instance TensorDataType S.Vector Int64 where
    decodeTensorData = simpleDecode
    encodeTensorData = simpleEncode

instance TensorDataType S.Vector Word8 where
    decodeTensorData = simpleDecode
    encodeTensorData = simpleEncode

instance TensorDataType S.Vector Word16 where
    decodeTensorData = simpleDecode
    encodeTensorData = simpleEncode

-- TODO: Haskell and tensorflow use different byte sizes for bools, which makes
-- encoding more expensive. It may make sense to define a custom boolean type.
instance TensorDataType S.Vector Bool where
    decodeTensorData =
        S.convert . S.map (/= 0) . FFI.tensorDataBytes . unTensorData
    encodeTensorData (Shape xs) =
        TensorData . FFI.TensorData xs DT_BOOL . S.map fromBool . S.convert
      where
        fromBool x = if x then 1 else 0 :: Word8

instance {-# OVERLAPPABLE #-} (Storable a, TensorDataType S.Vector a, TensorType a)
    => TensorDataType V.Vector a where
    decodeTensorData = (S.convert :: S.Vector a -> V.Vector a) . decodeTensorData
    encodeTensorData x = encodeTensorData x . (S.convert :: V.Vector a -> S.Vector a)

instance {-# OVERLAPPING #-} TensorDataType V.Vector (Complex Float) where
    decodeTensorData = error "TODO (Complex Float)"
    encodeTensorData = error "TODO (Complex Float)"

instance {-# OVERLAPPING #-} TensorDataType V.Vector (Complex Double) where
    decodeTensorData = error "TODO (Complex Double)"
    encodeTensorData = error "TODO (Complex Double)"

instance {-# OVERLAPPING #-} TensorDataType V.Vector ByteString where
    -- Encoded data layout (described in third_party/tensorflow/c/c_api.h):
    --   table offsets for each element :: [Word64]
    --   at each element offset:
    --     string length :: VarInt64
    --     string data   :: [Word8]
    decodeTensorData tensorData =
        either (\err -> error $ "Malformed TF_STRING tensor; " ++ err) id $
            if expected /= count
                then Left $ "decodeTensorData for ByteString count mismatch " ++
                            show (expected, count)
                else V.mapM decodeString (S.convert offsets)
      where
        expected = S.length offsets
        count = fromIntegral $ product $ FFI.tensorDataDimensions
                    $ unTensorData tensorData
        bytes = FFI.tensorDataBytes $ unTensorData tensorData
        offsets = S.take count $ S.unsafeCast bytes :: S.Vector Word64
        dataBytes = B.pack $ S.toList $ S.drop (count * 8) bytes
        decodeString :: Word64 -> Either String ByteString
        decodeString offset =
            let stringDataStart = B.drop (fromIntegral offset) dataBytes
            in Atto.eitherResult $ Atto.parse stringParser stringDataStart
        stringParser :: Atto.Parser ByteString
        stringParser = getVarInt >>= Atto.take . fromIntegral
    encodeTensorData (Shape xs) vec =
        TensorData $ FFI.TensorData xs dt byteVector
      where
        dt = tensorType (undefined :: ByteString)
        -- Add a string to an offset table and data blob.
        addString :: (Builder, Builder, Word64)
                  -> ByteString
                  -> (Builder, Builder, Word64)
        addString (table, strings, offset) str =
            ( table <> Builder.word64LE offset
            , strings <> lengthBytes <> Builder.byteString str
            , offset + lengthBytesLen + strLen
            )
          where
            strLen = fromIntegral $ B.length str
            lengthBytes = putVarInt $ fromIntegral $ B.length str
            lengthBytesLen =
                fromIntegral $ L.length $ Builder.toLazyByteString lengthBytes
        -- Encode all strings.
        (table', strings', _) = V.foldl' addString (mempty, mempty, 0) vec
        -- Concat offset table with data.
        bytes = table' <> strings'
        -- Convert to Vector Word8.
        byteVector = S.fromList $ L.unpack $ Builder.toLazyByteString bytes

newtype Scalar a = Scalar {unScalar :: a}
    deriving (Show, Eq, Ord, Num, Fractional, Floating, Real, RealFloat,
              RealFrac, IsString)

instance (TensorDataType V.Vector a, TensorType a) => TensorDataType Scalar a where
    decodeTensorData = Scalar . headFromSingleton . decodeTensorData
    encodeTensorData x (Scalar y) = encodeTensorData x (V.fromList [y])

headFromSingleton :: V.Vector a -> a
headFromSingleton x
    | V.length x == 1 = V.head x
    | otherwise = error $
                  "Unable to extract singleton from tensor of length "
                  ++ show (V.length x)

-- | Shape (dimensions) of a tensor.
--
-- TensorFlow supports shapes of unknown rank, which are represented as
-- @Nothing :: Maybe Shape@ in Haskell.
newtype Shape = Shape [Int64] deriving Show

instance IsList Shape where
    type Item Shape = Int64
    fromList = Shape . fromList
    toList (Shape ss) = toList ss

protoShape :: Lens' TensorShapeProto Shape
protoShape = iso protoToShape shapeToProto
  where
    protoToShape p = fromMaybe (error msg) (view protoMaybeShape p)
      where msg = "Can't convert TensorShapeProto with unknown rank to Shape: "
                  ++ showMessageShort p
    shapeToProto s' = def & protoMaybeShape .~ Just s'

protoMaybeShape :: Lens' TensorShapeProto (Maybe Shape)
protoMaybeShape = iso protoToShape shapeToProto
  where
    protoToShape :: TensorShapeProto -> Maybe Shape
    protoToShape p =
        if view unknownRank p
            then Nothing
            else Just (Shape (p ^.. dim . traverse . size))
    shapeToProto :: Maybe Shape -> TensorShapeProto
    shapeToProto Nothing =
        def & unknownRank .~ True
    shapeToProto (Just (Shape ds)) =
        def & dim .~ fmap (\d -> def & size .~ d) ds

class Attribute a where
    attrLens :: Lens' AttrValue a

instance Attribute Float where
    attrLens = f

instance Attribute ByteString where
    attrLens = s

instance Attribute Int64 where
    attrLens = i

instance Attribute DataType where
    attrLens = type'

instance Attribute TensorProto where
    attrLens = tensor

instance Attribute Bool where
    attrLens = b

instance Attribute Shape where
    attrLens = shape . protoShape

instance Attribute (Maybe Shape) where
    attrLens = shape . protoMaybeShape

-- TODO(gnezdo): support generating list(Foo) from [Foo].
instance Attribute AttrValue'ListValue where
    attrLens = list

instance Attribute [DataType] where
    attrLens = list . type'

instance Attribute [Int64] where
    attrLens = list . i

-- | A heterogeneous list type.
data ListOf f as where
    Nil :: ListOf f '[]
    (:/) :: f a -> ListOf f as -> ListOf f (a ': as)

infixr 5 :/

type family All f as :: Constraint where
    All f '[] = ()
    All f (a ': as) = (f a, All f as)

type family Map f as where
    Map f '[] = '[]
    Map f (a ': as) = f a ': Map f as

instance All Eq (Map f as) => Eq (ListOf f as) where
    Nil == Nil = True
    (x :/ xs) == (y :/ ys) = x == y && xs == ys
    -- Newer versions of GHC use the GADT to tell that the previous cases are
    -- exhaustive.
#if __GLASGOW_HASKELL__ < 800
    _ == _ = False
#endif

instance All Show (Map f as) => Show (ListOf f as) where
    showsPrec _ Nil = showString "Nil"
    showsPrec d (x :/ xs) = showParen (d > 10)
                                $ showsPrec 6 x . showString " :/ "
                                    . showsPrec 6 xs

type List = ListOf Identity

-- | Equivalent of ':/' for lists.
(/:/) :: a -> List as -> List (a ': as)
(/:/) = (:/) . Identity

infixr 5 /:/

-- | A 'Constraint' specifying the possible choices of a 'TensorType'.
--
-- We implement a 'Constraint' like @OneOf '[Double, Float] a@ by turning the
-- natural representation as a conjunction, i.e.,
--
-- @
--    a == Double || a == Float
-- @
--
-- into a disjunction like
--
-- @
--     a \/= Int32 && a \/= Int64 && a \/= ByteString && ...
-- @
--
-- using an enumeration of all the possible 'TensorType's.
type OneOf ts a
    -- Assert `TensorTypes' ts` to make error messages a little better.
    = (TensorType a, TensorTypes' ts, NoneOf (AllTensorTypes \\ ts) a)

type OneOfs ts as = (TensorTypes as, TensorTypes' ts,
                        NoneOfs (AllTensorTypes \\ ts) as)

type family NoneOfs ts as :: Constraint where
    NoneOfs ts '[] = ()
    NoneOfs ts (a ': as) = (NoneOf ts a, NoneOfs ts as)

data TensorTypeProxy a where
    TensorTypeProxy :: TensorType a => TensorTypeProxy a

type TensorTypeList = ListOf TensorTypeProxy

fromTensorTypeList :: TensorTypeList ts -> [DataType]
fromTensorTypeList Nil = []
fromTensorTypeList ((TensorTypeProxy :: TensorTypeProxy t) :/ ts)
    = tensorType (undefined :: t) : fromTensorTypeList ts

fromTensorTypes :: forall as . TensorTypes as => Proxy as -> [DataType]
fromTensorTypes _ = fromTensorTypeList (tensorTypes :: TensorTypeList as)

class TensorTypes (ts :: [*]) where
    tensorTypes :: TensorTypeList ts

instance TensorTypes '[] where
    tensorTypes = Nil

-- | A constraint that the input is a list of 'TensorTypes'.
instance (TensorType t, TensorTypes ts) => TensorTypes (t ': ts) where
    tensorTypes = TensorTypeProxy :/ tensorTypes

-- | A simpler version of the 'TensorTypes' class, that doesn't run
-- afoul of @-Wsimplifiable-class-constraints@.
--
-- In more detail: the constraint @OneOf '[Double, Float] a@ leads
-- to the constraint @TensorTypes' '[Double, Float]@, as a safety-check
-- to give better error messages.  However, if @TensorTypes'@ were a class,
-- then GHC 8.2.1 would complain with the above warning unless @NoMonoBinds@
-- were enabled.  So instead, we use a separate type family for this purpose.
-- For more details: https://ghc.haskell.org/trac/ghc/ticket/11948
type family TensorTypes' (ts :: [*]) :: Constraint where
    -- Specialize this type family when `ts` is a long list, to avoid deeply
    -- nested tuples of constraints.  Works around a bug in ghc-8.0:
    -- https://ghc.haskell.org/trac/ghc/ticket/12175
    TensorTypes' (t1 ': t2 ': t3 ': t4 ': ts)
        = (TensorType t1, TensorType t2, TensorType t3, TensorType t4
              , TensorTypes' ts)
    TensorTypes' (t1 ': t2 ': t3 ': ts)
        = (TensorType t1, TensorType t2, TensorType t3, TensorTypes' ts)
    TensorTypes' (t1 ': t2 ': ts)
        = (TensorType t1, TensorType t2, TensorTypes' ts)
    TensorTypes' (t ': ts) = (TensorType t, TensorTypes' ts)
    TensorTypes' '[] = ()

-- | A constraint checking that two types are different.
type family a /= b :: Constraint where
    a /= a = TypeError a ~ ExcludedCase
    a /= b = ()

-- | Helper types to produce a reasonable type error message when the Constraint
-- "a /= a" fails.
-- TODO(judahjacobson): Use ghc-8's CustomTypeErrors for this.
data TypeError a
data ExcludedCase

-- | An enumeration of all valid 'TensorType's.
type AllTensorTypes =
    -- NOTE: This list should be kept in sync with
    -- TensorFlow.OpGen.dtTypeToHaskell.
    -- TODO: Add support for Complex Float/Double.
    '[ Float
     , Double
     , Int8
     , Int16
     , Int32
     , Int64
     , Word8
     , Word16
     , ByteString
     , Bool
     ]

-- | Removes a type from the given list of types.
type family Delete a as where
    Delete a '[] = '[]
    Delete a (a ': as) = Delete a as
    Delete a (b ': as) = b ': Delete a as

-- | Takes the difference of two lists of types.
type family as \\ bs where
    as \\ '[] = as
    as \\ (b ': bs) = Delete b as \\ bs

-- | A constraint that the type @a@ doesn't appear in the type list @ts@.
-- Assumes that @a@ and each of the elements of @ts@ are 'TensorType's.
type family NoneOf ts a :: Constraint where
    -- Specialize this type family when `ts` is a long list, to avoid deeply
    -- nested tuples of constraints.  Works around a bug in ghc-8.0:
    -- https://ghc.haskell.org/trac/ghc/ticket/12175
    NoneOf (t1 ': t2 ': t3 ': t4 ': ts) a
        = (a /= t1, a /= t2, a /= t3, a /= t4, NoneOf ts a)
    NoneOf (t1 ': t2 ': t3 ': ts) a = (a /= t1, a /= t2, a /= t3, NoneOf ts a)
    NoneOf (t1 ': t2 ': ts) a = (a /= t1, a /= t2, NoneOf ts a)
    NoneOf (t1 ': ts) a = (a /= t1, NoneOf ts a)
    NoneOf '[] a = ()

max_stars_repo_name: cmk/haskell

max_issues_repo_name: cmk/haskell

max_forks_repo_name: cmk/haskell


text:
module Main where

import Numeric.LinearAlgebra
import Hoqus.Dirac

main = do
  let dim = 4
  let v1 = ket dim 0
  let v2 = ket dim 4
  let prj = proj 4 0 0
  print v1
  print v2
  print $ v1 + v2
  print $ prj #> v1
  print $ prj #> v2

ext: hs

max_stars_repo_name: jmiszczak/hoqus

max_issues_repo_name: jmiszczak/hoqus

max_forks_repo_name: jmiszczak/hoqus


text:
{-|
Module      : ODMatrix
Description : Representación y procesamiento de matrices origen-destino.

Base de funciones para la operación sobre matrices origen-destino.

ODM Features
The ODM features are indicators to verify if two ODMs are elements of the same solution set.
-}
module ODMatrix (
    ODM
  , FeatVec
  --
  , boardingVector
  , alightingVector
  , onBoardVector
  --
  , inputVector
  --
  , getFeatures
  , compareFeatures
  --
  , toFlatList
  , fromFlatList
  ) where

  import Numeric.LinearAlgebra (Matrix, Vector,
                                rows, cols, flatten,
                                toColumns, toRows,
                                subMatrix, vector,
                                vjoin, subVector)
  import Numeric.LinearAlgebra.Data (toList, takeDiag, fromLists, toLists, assoc)

  import Control.Monad (join)

  import ODMatrix.Util.TriangularNumbers (isTrg, trglvl)

  -- | Representación de una matriz de viajes. Debe ser cuadrada.
  type ODM = Matrix Double

  -- | The feature represents the input information.
  type FeatVec = Vector Double

  -- * Cálculo de vectores de información derivados. Feature Vectors.

  -- | Dada una matriz de viajes, calcula el vector de ascensos.
  boardingVector :: ODM      -- ^ Matriz de viajes de entrada.
                 -> FeatVec  -- ^ Vector de ascensos
  boardingVector = vector . map (sum . toList) . toRows

  -- | Dada una matriz de viajes, calcula el vector de descensos.
  alightingVector :: ODM      -- ^ Matriz de viajes de entrada.
                  -> FeatVec  -- ^ Vector de ascensos
  alightingVector = vector . map (sum . toList) . toColumns

  -- | Dada una matriz de viajes, calcula el vector de pasajeros a bordo.
  onBoardVector :: ODM      -- ^ Matriz de viajes de entrada.
                -> FeatVec  -- ^ Vector de pasajeros a bordo.
  onBoardVector m = vector ((sum . toList) . flatten <$> blocks)
    where s = cols m
          blocks = map (\x -> subMatrix (0,x) (x+1,s-x) m) [0..(s-1)]

  -- | Get all the feature vectors. (boarding, alighting, onBoard)
  getFeatures :: ODM -> (FeatVec, FeatVec, FeatVec)
  getFeatures m = (boardingVector m, alightingVector m, onBoardVector m)

  -- | Given two ODMs, verify if they are members of the same solution set.
  compareFeatures :: ODM -> ODM -> Bool
  compareFeatures a b = getFeatures a == getFeatures b

  -- | Given a ODM, build the vector of independent terms
  inputVector :: ODM -> Vector Double
  inputVector odm = vjoin [a,b]
    where a = alightingVector odm
          b = subVector 1 (cols odm - 1) $ boardingVector odm

  -- * Convertions

  -- | Aplana una matriz de origen destino, por filas y solo el triangulo superior.
  toFlatList :: ODM
             -> [Double]
  toFlatList = join . zipWith drop [0..] . toLists

  -- | Convierte una lista en una matriz origen destino.
  fromFlatList :: [Double]
               -> ODM
  fromFlatList ls
    | (not . isTrg) (length ls) =
      error "fromFlatList: Invalid list. The length must be a triangular number."
    | otherwise =
      assoc (n,n) 0 . zip [(i,j)|i<-[0..n-1],j<-[i..n-1]] $ ls
        where n = trglvl $ length ls

max_stars_repo_name: renecura/odmatrix

max_issues_repo_name: renecura/odmatrix

max_forks_repo_name: renecura/odmatrix


text:
module Fuml.Base.PCA where

import Numeric.LinearAlgebra
import Data.List (foldl')

sumWith :: Num a => (Vector Double -> a) -> [Vector Double] -> a
sumWith f (v:vs) = foldl' (\acc v -> acc + f v) (f v) vs

meanVs :: [Vector Double] -> Vector Double
meanVs vecs = scale (recip $ realToFrac $ length vecs) $ sumWith id vecs

covVs :: Vector Double -> [Vector Double] -> Herm Double
covVs xmn xs
   = let k = length xs
         f xi = outerSame $ xi - xmn
     in sym $ scale (recip $ realToFrac k - 1 ) $ sumWith f xs

outerSame v = outer v v

--following https://github.com/albertoruiz/hmatrix/blob/master/examples/pca2.hs
type Stat = (Vector Double, [Double], Matrix Double)

statVs :: [Vector Double] -> Stat
statVs vs = (m, toList s, tr' v) where
    m = meanVs vs
    (s,v) = eigSH (covVs m vs)

pcaN :: Int -> Stat -> (Vector Double -> Vector Double , Vector Double -> Vector Double )
pcaN n (m,s,v) = (encode,decode)
  where
    encode x = vp #> (x - m)
    decode x = x <# vp + m
    vp = takeRows n v

max_stars_repo_name: ekalosak/open

max_issues_repo_name: ekalosak/open

max_forks_repo_name: ekalosak/open


text:
{-|
Module: Science.QuantumChemistry.NumericalTools.EigenValues
Description: Compute all the eigenvalues of the real symmetrical matrix
Copyright: @2016 Felipe Zapata
-}

module Science.QuantumChemistry.NumericalTools.EigenValues (
                   eigenSolve
                   ) where

-- =============================> Standard and third party libraries <===============================
import Control.Arrow ((***))
import Data.Array.Repa as R
import Data.Vector.Storable (convert)
import qualified Data.Vector.Unboxed as U
import qualified Numeric.LinearAlgebra.Data as ND
import Numeric.LinearAlgebra.HMatrix (eigSH')

-- =================> Internal Modules <======================
import Science.QuantumChemistry.GlobalTypes (VecUnbox)
import Science.QuantumChemistry.NumericalTools.JacobiMethod (jacobiP)
import Science.QuantumChemistry.NumericalTools.VectorTools (sortEigenData)

-- | Compute all the eigenvalues and eigenvectors of symmetrical matrix.
eigenSolve :: Array U DIM2 Double -> (VecUnbox, Array U DIM2 Double)
eigenSolve arr = (vs, fromUnboxed (extent arr) mtx )
  where (vs,mtx) = sortEigenData $ hv2VecUnbox *** hm2VecUnbox $ eigSH' $ repa2HM arr

vecUnbox2HV :: VecUnbox -> ND.Vector Double
vecUnbox2HV = convert

hv2VecUnbox :: ND.Vector Double -> VecUnbox
hv2VecUnbox = convert

repa2HM :: Array U DIM2 Double -> ND.Matrix Double
repa2HM arr = ND.reshape dim . convert . toUnboxed $ arr

 where (Z :. dim :. _) = extent arr

hm2Repa :: ND.Matrix Double -> Array U DIM2 Double
hm2Repa mtx = fromUnboxed (ix2 dim dim) . convert . ND.flatten $ mtx
 where dim = ND.rows mtx

hm2VecUnbox :: ND.Matrix Double -> VecUnbox
hm2VecUnbox = convert . ND.flatten

max_stars_repo_name: felipeZ/Haskell-abinitio

max_issues_repo_name: felipeZ/Haskell-abinitio

max_forks_repo_name: felipeZ/Haskell-abinitio


text:
module Main where

import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Cairo

import Statistics.Distribution.Normal
import Statistics.Distribution

-- pOak = class1
-- pBeech = class2

stdDevOak :: NormalDistribution
stdDevOak = normalDistr (- 2.25) 1

stdDevBeech :: NormalDistribution
stdDevBeech = normalDistr 1.596 1

type Priori = Double

pOak :: Priori
pOak = 0.25

pBeech :: Priori
pBeech = 0.75

oak :: [Double] -> [(Double,Double)]
oak xs = [ (x,(density stdDevOak x) * pOak) | x <- xs ]

beech :: [Double] -> [(Double, Double)]
beech xs = [ (x,(density stdDevBeech x) * pBeech) | x <- xs ]

identit :: [Double] -> [(Double, Double)]
identit xs = [(x, x) | x <- xs]

increments :: Double -> Double -> Double -> [Double]
increments from to inc
    | from > to = []
    | otherwise = (from):(increments (from + inc) to inc)

--main = print $ (increments (-5) 5 0.1)

main = toFile def "output.png" $ do
    layout_title .= "KMeans"
    setShapes [PointShapeCircle, PointShapePlus, PointShapeStar]
    setColors [opaque blue, opaque red]
    plot (points "class2" $ beech (increments (-5) 5 0.1))
    plot (line "class1" [oak (increments (-5) 5 0.1)])

max_stars_repo_name: Parrows/Parrows

max_issues_repo_name: Parrows/Parrows

max_forks_repo_name: Parrows/Parrows


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE BangPatterns #-}
module HLearn.Models.Distributions.Univariate.Binomial
    where

import HLearn.Algebra
import HLearn.Models.Distributions.Common
import HLearn.Models.Distributions.Univariate.Internal.Moments
import HLearn.Models.Distributions.Visualization.Gnuplot

import qualified Statistics.Distribution as S
import Statistics.Distribution.Binomial

-------------------------------------------------------------------------------
-- data types

newtype Binomial prob dp = Binomial {  bmoments :: (Moments3 dp) }
    deriving (Read,Show,Eq,Ord,Monoid,Group)

-------------------------------------------------------------------------------
-- Training

instance (Num dp) => HomTrainer (Binomial prob dp) where
    type Datapoint (Binomial prob dp) = dp
    train1dp dp = Binomial $ train1dp dp

-------------------------------------------------------------------------------
-- distribution

instance (Num dp) => Probabilistic (Binomial prob dp) where
    type Probability (Binomial prob dp) = prob

instance (Floating prob) => PDF (Binomial Double Int) where
    pdf (Binomial dist) dp = S.probability (binomial n p) dp
        where
            n = bin_n $ Binomial dist
            p = bin_p $ Binomial dist

bin_n :: Binomial Double Int -> Int
bin_n (Binomial dist) = round $ ((fromIntegral $ m1 dist :: Double) / (fromIntegral $ m0 dist)) / (bin_p $ Binomial dist)

bin_p :: Binomial Double Int -> Double
bin_p (Binomial dist) = ((fromIntegral $ m1 dist) / (fromIntegral $ m0 dist)) + 1 - (fromIntegral $ m2 dist)/(fromIntegral $ m1 dist)

instance
    ( PDF (Binomial prob dp)
--     , PlottableDataPoint dp
    , Show prob
    , Show dp
    , Ord dp
    , Ord prob
    , Num prob
    , Integral dp
    ) => PlottableDistribution (Binomial prob dp)
-- instance PlottableDistribution (Poisson Int Double)
        where

    plotType _ = Points

    samplePoints dist = [min..max]
        where
            min = 0
            max = maximum [20,floor $ 3*mu]
            mu = 5

-- instance (Fractional prob) => Mean (Binomial prob) where
--     mean (Binomial dist) =
--
-- instance (Fractional prob) => Variance (Binomiral prob) where
--     variance dist = mean dist

max_stars_repo_name: silky/HLearn

max_issues_repo_name: silky/HLearn

max_forks_repo_name: silky/HLearn


text:
{-# LANGUAGE MonoLocalBinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE Strict           #-}
{-# LANGUAGE StrictData       #-}
module FourierPinwheel.Array where

import qualified Data.Array.Accelerate                  as A
import qualified Data.Array.Accelerate.LLVM.PTX as A
import           Data.Array.Repa                        as R
import           Data.Complex
import           Data.List                              as L
import           Data.Vector.Generic                    as VG
import           Data.Vector.Storable                   as VS
import           DFT.Plan
import           Filter.Utils
import           Graphics.Rendering.Chart.Backend.Cairo
import           Graphics.Rendering.Chart.Easy
import           Image.IO
import           Text.Printf
import           Utils.Array
import           Utils.BLAS
import           Utils.Parallel
import FourierMethod.FourierSeries2D

data FPArray vector = FPArray
  { getFPArrayNumXFreq     :: Int
  , getFPArrayNumYFreq     :: Int
  , getFPArrayNumRFreq     :: Int
  , getFPArrayNumThetaFreq :: Int
  , getFPArrayNumRhoFreq   :: Int
  , getFPArrayNumPhiFreq   :: Int
  , getFPArray             :: [vector]
  }

{-# INLINE toMatrixAcc #-}
toMatrixAcc ::
     (VG.Vector vector e, A.Elt e) => FPArray (vector e) -> A.Acc (A.Array A.DIM2 e)
toMatrixAcc (FPArray numXFreq numYFreq numRFreq numThetaFreq _ _ arr) =
  A.use .
  A.fromList (A.Z A.:. (numRFreq * numThetaFreq) A.:. (numXFreq * numYFreq)) .
  VG.toList . VG.concat $
  arr

{-# INLINE parMapFPArray #-}
parMapFPArray ::
     (NFData vector) => (vector -> vector) -> FPArray vector -> FPArray vector
parMapFPArray f (FPArray numXFreq numYFreq numRFreq numThetaFreq numRhoFreq numPhiFreq vecs) =
  FPArray numXFreq numYFreq numRFreq numThetaFreq numRhoFreq numPhiFreq .
  parMap rdeepseq f $
  vecs

{-# INLINE parZipWithFPArray #-}
parZipWithFPArray ::
     (NFData vector)
  => (vector -> vector -> vector)
  -> FPArray vector
  -> FPArray vector
  -> FPArray vector
parZipWithFPArray f (FPArray numXFreq numYFreq numRFreq numThetaFreq numRhoFreq numPhiFreq vecs1) =
  FPArray numXFreq numYFreq numRFreq numThetaFreq numRhoFreq numPhiFreq .
  parZipWith rdeepseq f vecs1 . getFPArray

plotFPArray ::
     DFTPlan
  -> FilePath
  -> FPArray (VS.Vector (Complex Double))
  -> IO (R.Array R.U R.DIM4 (Complex Double))
plotFPArray plan filePath arr = do
  let planIDBackward =
        DFTPlanID
          IDFT1DG
          [ getFPArrayNumThetaFreq arr
          , getFPArrayNumXFreq arr
          , getFPArrayNumYFreq arr
          ]
          [1, 2]
  vecsR2 <- dftExecuteBatchP plan planIDBackward . getFPArray $ arr
  arrR2 <-
    R.computeUnboxedP .
    makeFilter2DInverse .
    R.fromUnboxed
      (R.Z R.:. getFPArrayNumRFreq arr R.:. getFPArrayNumThetaFreq arr R.:.
       getFPArrayNumXFreq arr R.:.
       getFPArrayNumYFreq arr) .
    VS.convert . VS.concat $
    vecsR2
  arr1 <-
    R.sumP . R.sumS . R.map (\x -> magnitude x Prelude.^ 2) . rotate4D2 $ arrR2
  plotImageRepa filePath .
    ImageRepa 8 .
    R.fromUnboxed
      (R.Z R.:. (1 :: Int) R.:. getFPArrayNumXFreq arr R.:.
       getFPArrayNumYFreq arr) .
    -- VG.map (\x -> log (x + 1)) .
    R.toUnboxed $ arr1
  return arrR2

plotFPArrayAcc ::
     (A.Elt (Complex Double))
  => [A.PTX]
  -> FilePath
  -> Int
  -> Double
  -> Double
  -> Int
  -> FPArray (VS.Vector (Complex Double))
  -> IO (R.Array R.U R.DIM4 (Complex Double))
plotFPArrayAcc ptxs filePath numPoints delta periodR2 numBatch arr@(FPArray numXFreq numYFreq numRFreq numThetaFreq _ _ vecs) = do
  let rows = numRFreq * numThetaFreq
      cols = numXFreq * numYFreq
      arrMat =
        A.transpose .
        A.use .
        A.fromList (A.Z A.:. rows A.:. cols) .
        R.toList .
        makeFilter2DInverse .
        fromUnboxed (Z :. rows :. numXFreq :. numXFreq) . VG.convert . VG.concat $
        vecs
  arrR2 <-
    computeFourierSeriesR2StreamAcc
      ptxs
      numXFreq
      numPoints
      rows
      periodR2
      delta
      numBatch
      arrMat
  (sumP . R.map (\x -> magnitude x ** 2) . rotate3D $ arrR2) >>=
    plotImageRepa filePath .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numPoints :. numPoints) . toUnboxed
  return .
    computeS .
    R.reshape (Z :. numRFreq :. numThetaFreq :. numPoints :. numPoints) $
    arrR2

plotFPArrayFreqency ::
     FilePath
  -> FPArray (VS.Vector (Complex Double))
  -> IO ()
plotFPArrayFreqency filePath arr = do
  arr1 <-
    R.sumP .
    R.sumS .
    R.map (\x -> magnitude x Prelude.^ 2) .
    rotate4D2 .
    makeFilter2DInverse .
    R.fromUnboxed
      (R.Z R.:. getFPArrayNumRFreq arr R.:. getFPArrayNumThetaFreq arr R.:.
       getFPArrayNumXFreq arr R.:.
       getFPArrayNumYFreq arr) .
    VS.convert . VS.concat . getFPArray $
    arr
  plotImageRepa filePath .
    ImageRepa 8 .
    R.fromUnboxed
      (R.Z R.:. (1 :: Int) R.:. getFPArrayNumXFreq arr R.:.
       getFPArrayNumYFreq arr) .
    -- VG.map (\x -> x^2) .
    R.toUnboxed $
    arr1

plotRThetaDist ::
     FilePath
  -> FilePath
  -> Int
  -> Int
  -> Int
  -> Double
  -> (Int, Int)
  -> R.Array R.U R.DIM4 (Complex Double)
  -> IO ()
plotRThetaDist filePathTheta filePathR numPoints numTheta numR periodEnv (x, y) arr = do
  let (Z :. numRFreq :. numThetaFreq :. _ :. _) = extent arr
      center = div numPoints 2
      centerR = div numR 2
      centerTheta = div numTheta 2
      centerRFreq = div numRFreq 2
      centerThetaFreq = div numThetaFreq 2
      vecRTheta =
        VG.convert . R.toUnboxed . R.computeUnboxedS . R.slice arr $
        (R.Z R.:. R.All R.:. R.All R.:. (x + center) R.:. (y + center))
      deltaTheta = 2 * pi / Prelude.fromIntegral numTheta
      deltaLogR = log periodEnv / Prelude.fromIntegral numR
      matrix =
        VG.convert .
        R.toUnboxed .
        R.computeUnboxedS .
        R.fromFunction (R.Z R.:. numR R.:. numTheta :. numRFreq :. numThetaFreq) $ \(Z :. r :. t :. rf :. tf) ->
          cis
            (fromIntegral ((r - centerR) * (rf - centerRFreq)) * 2 * pi /
             log periodEnv *
             deltaLogR +
             fromIntegral ((t - centerTheta) * (tf - centerThetaFreq)) *
             deltaTheta)
      m = numR * numTheta
      n = 1
      k = numRFreq * numThetaFreq
  vec <- VS.map magnitude <$> gemmBLAS m n k matrix vecRTheta
  let outputArr = fromUnboxed (Z :. numR :. numTheta) . VG.convert $ vec
      outputR = R.toList . sumS $ outputArr
      outputTheta =
        R.toList .
        sumS .
        R.backpermute (Z :. numTheta :. numR) (\(Z :. i :. j) -> (Z :. j :. i)) $
        outputArr
  toFile def filePathTheta $ do
    layout_title .= printf "Theta (%d , %d)" x y
    plot
      (line
         ""
         [ L.zip
             [ fromIntegral (i - centerTheta) * deltaTheta
             | i <- [0 .. numTheta - 1]
             ]
             outputTheta
         ])
  toFile def filePathR $ do
    layout_title .= printf "LogR (%d , %d)" x y
    plot
      (line
         ""
         [ L.zip
             [fromIntegral (i - centerR) * deltaLogR | i <- [0 .. numR - 1]]
             outputR
         ])

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE ScopedTypeVariables #-}
module Main
    ( main
    ) where

-------------------------------------------------------------------------------
import           Data.Foldable             as FT
import           Data.Ratio
import           Data.Sequence             as Seq
import           Test.Tasty
import           Test.Tasty.HUnit
import           Test.Tasty.QuickCheck
-------------------------------------------------------------------------------
import           Statistics.RollingAverage
-------------------------------------------------------------------------------

main :: IO ()
main = defaultMain tests

tests :: TestTree
tests = testGroup "rolling-average"
  [
    testProperty "window size of 0 always produces 0" $ \(samples :: [Int]) ->
      ravg (avgList 0 samples) === (0.0 :: Double)
  , testProperty "never grows samples larger than the limit" $ \(Positive lim) (samples :: [Int]) ->
      Seq.length (ravgSamples (avgList lim samples)) <= lim
  , testCase "can calculate rolling average as a Ratio" $ do
      ravg (avgList 3 ([0, 1, 2, 3, 4, 5] :: [Int])) @?= (((3 + 4 + 5) % 3) :: Rational)
  ]

-------------------------------------------------------------------------------
avgList :: Num a => Int -> [a] -> RollingAvg a
avgList lim samples = FT.foldl' ravgAdd (mkRavg lim) samples

max_stars_repo_name: Soostone/rolling-average

max_issues_repo_name: Soostone/rolling-average

max_forks_repo_name: Soostone/rolling-average


text:
--------------------------------------------------------------------------------
-- |
-- Module      : ArrayFire
-- Copyright   : David Johnson (c) 2019-2020
-- License     : BSD3
-- Maintainer  : David Johnson <djohnson.m@gmail.com>
-- Stability   : Experimental
-- Portability : GHC
--
-- <<https://user-images.githubusercontent.com/875324/59819703-0fbaf980-92f7-11e9-8f53-adebea590bfb.png>>
--
--------------------------------------------------------------------------------
module ArrayFire
  ( -- * Tutorial
    -- $tutorial

    -- ** Modules
    -- $modules

    -- ** Exceptions
    -- $exceptions

    -- ** Construction
    -- $construction

    -- ** Laws
    -- $laws

    -- ** Conversion
    -- $conversion

    -- ** Serialization
    -- $serialization

    -- ** Device
    -- $device
    module ArrayFire.Algorithm
  , module ArrayFire.Arith
  , module ArrayFire.Array
  , module ArrayFire.Backend
  , module ArrayFire.BLAS
  , module ArrayFire.Data
  , module ArrayFire.Device
  , module ArrayFire.Features
  , module ArrayFire.Graphics
  , module ArrayFire.Image
  , module ArrayFire.Index
  , module ArrayFire.LAPACK
  , module ArrayFire.Random
  , module ArrayFire.Signal
  , module ArrayFire.Sparse
  , module ArrayFire.Statistics
  , module ArrayFire.Types
  , module ArrayFire.Util
  , module ArrayFire.Vision
  , module Foreign.C.Types
  , module Data.Int
  , module Data.Word
  , module Data.Complex
  , module Foreign.Storable
  ) where

import ArrayFire.Algorithm
import ArrayFire.Arith
import ArrayFire.Array
import ArrayFire.Backend
import ArrayFire.BLAS
import ArrayFire.Data
import ArrayFire.Device
import ArrayFire.Features
import ArrayFire.Graphics
import ArrayFire.Image
import ArrayFire.Index
import ArrayFire.LAPACK
import ArrayFire.Random
import ArrayFire.Signal
import ArrayFire.Sparse
import ArrayFire.Statistics
import ArrayFire.Types
import ArrayFire.Util
import ArrayFire.Vision
import ArrayFire.Orphans ()
import Foreign.Storable
import Foreign.C.Types
import Data.Int
import Data.Complex
import Data.Word

-- $tutorial
--
-- [ArrayFire](http://arrayfire.org/docs/gettingstarted.htm) is a high performance parallel computing library that features modules for statistical and numerical methods.
-- Example usage is depicted below.
--
-- @
-- module Main where
--
-- import qualified ArrayFire as A
--
-- main :: IO ()
-- main = print $ A.matrix @Double (3,2) [[1,2,3],[4,5,6]]
-- @
--
-- Each 'Array' is constructed and displayed in column-major order.
--
-- @
-- ArrayFire Array
-- [3 2 1 1]
--    1.0000     4.0000
--    2.0000     5.0000
--    3.0000     6.0000
-- @

-- $modules
--
-- All child modules are re-exported top-level in the "ArrayFire" module.
-- We recommend importing "ArrayFire" qualified so as to avoid naming collisions.
--
-- >>> import qualified ArrayFire as A
--

-- $exceptions
--
-- @
-- {\-\# LANGUAGE TypeApplications \#\-}
-- module Main where
--
-- import qualified ArrayFire         as A
-- import           Control.Exception ( catch )
--
-- main :: IO ()
-- main = A.printArray action \`catch\` (\\(e :: A.AFException) -> print e)
--   where
--     action =
--       A.matrix \@Double (3,3) [[1..],[1..],[1..]]
--         \`A.mul\` A.matrix \@Double (2,2) [[1..],[1..]]
-- @
--
-- The above operation is invalid since the matrix multiply has improper dimensions. The caught exception produces the following error:
--
-- > AFException {afExceptionType = SizeError, afExceptionCode = 203, afExceptionMsg = "Invalid input size"}
--

-- $construction
-- An 'Array' can be constructed using the following smart constructors:
--
-- /Note/: All smart constructors (and ArrayFire internally) assume column-major order.
--
-- @
-- >>> scalar \@Double 2.0
-- ArrayFire Array
-- [1 1 1 1]
--    2.0000
-- @
--
-- @
-- >>> vector \@Double 10 [1..]
-- ArrayFire Array
-- [10 1 1 1]
--     1.0000
--     2.0000
--     3.0000
--     4.0000
--     5.0000
--     6.0000
--     7.0000
--     8.0000
--     9.0000
--    10.0000
-- @
--
-- @
-- >>> matrix \@Double (2,2) [[1,2],[3,4]]
-- ArrayFire Array
-- [2 2 1 1]
--     1.0000     3.0000
--     2.0000     4.0000
-- @
--
-- @
-- >>> cube \@Double (2,2,2) [[[2,2],[2,2]],[[2,2],[2,2]]]
-- ArrayFire Array
-- [2 2 2 1]
--    2.0000     2.0000
--    2.0000     2.0000
--
--    2.0000     2.0000
--    2.0000     2.0000
-- @
--
-- @
-- >>> tensor \@Double (2,2,2,2) [[[[2,2],[2,2]],[[2,2],[2,2]]], [[[2,2],[2,2]],[[2,2],[2,2]]]]
-- ArrayFire Array
-- [2 2 2 2]
--     2.0000     2.0000
--     2.0000     2.0000
--
--     2.0000     2.0000
--     2.0000     2.0000
--
--
--     2.0000     2.0000
--     2.0000     2.0000
--
--     2.0000     2.0000
--     2.0000     2.0000
-- @
--
-- Array construction can use Haskell's lazy lists, since 'take' is called on each dimension before sending to the C API.
--
-- >>> mkArray @Double [5,3] [1..]
-- ArrayFire Array
-- [5 3 1 1]
--    1.0000     6.0000    11.0000
--    2.0000     7.0000    12.0000
--    3.0000     8.0000    13.0000
--    4.0000     9.0000    14.0000
--    5.0000    10.0000    15.0000
--
-- Specifying up to 4 dimensions is allowed (anything higher is ignored).

-- $laws
-- Every 'Array' has an instance of 'Eq', 'Num', 'Fractional', 'Floating' and 'Show'
--
-- 'Num'
--
-- >>> 2.0 :: Array Double
-- ArrayFire Array
-- [1 1 1 1]
--    2.0000
--
-- >>> scalar @Int 1 + scalar @Int 1
-- ArrayFire Array
-- [1 1 1 1]
--         2
--
-- >>> scalar @Int 1 - scalar @Int 1
-- ArrayFire Array
-- [1 1 1 1]
--         0
--
-- >>> scalar @Double 10 / scalar @Double 10
-- ArrayFire Array
-- [1 1 1 1]
--   1.0000
--
-- >>> abs $ scalar @Double (-10)
-- ArrayFire Array
-- [1 1 1 1]
--   10.0000
--
-- >>> negate (scalar @Double 10)
-- ArrayFire Array
-- [1 1 1 1]
--   -10.0000
--
-- >>> fromInteger 1.0 :: Array Double
-- ArrayFire Array
-- [1 1 1 1]
--    1.0000
--
-- 'Eq'
--
-- >>> scalar @Double 1 [10] == scalar @Double 1 [10]
-- True
-- >>> scalar @Double 1 [10] /= scalar @Double 1 [10]
-- False
--
-- 'Floating'
--
-- >>> pi :: Array Double
-- ArrayFire Array
-- [1 1 1 1]
--    3.1416
--
-- >>> A.sqrt pi :: Array Double
-- ArrayFire Array
-- [1 1 1 1]
--    1.7725
--
-- 'Fractional'
--
-- >>> (pi :: Array Double) / pi
-- ArrayFire Array
-- [1 1 1 1]
--    1.000
--
-- >>> recip 0.5 :: Array Double
-- ArrayFire Array
-- [1 1 1 1]
--    2.000
--
-- 'Show'
--
-- >>> 0.0 :: Array Double
-- ArrayFire Array
-- [1 1 1 1]
--    0.000
--

-- $conversion
-- Any 'Array' can be exported into Haskell using `toVector'. This will create a Storable vector suitable for use in other C programs.
--
-- >>> vector :: Vector Double <- toVector <$> randu @Double [10,10]
--

-- $serialization
-- Each 'Array' can be serialized to disk and deserialized from disk efficiently.
--
-- @
-- import qualified ArrayFire as A
-- import           Control.Monad
--
-- main :: IO ()
-- main = do
--   let arr = A.'constant' [1,1,1,1] 10
--   idx <- A.'saveArray' "key" arr "file.array" False
--   foundIndex <- A.'readArrayKeyCheck' "file.array" "key"
--   when (idx == foundIndex) $ do
--     array <- A.'readArrayKey' "file.array" "key"
--     'print' array
--
-- -- ArrayFire Array
-- -- [ 1 1 1 1 ]
-- --         10
-- @
--

-- $device
-- The ArrayFire API is able to see which devices are present, and will by default use the GPU if available.
--
-- >>> afInfo
-- ArrayFire v3.6.4 (OpenCL, 64-bit Mac OSX, build 1b8030c5)
-- [0] APPLE: AMD Radeon Pro 555X Compute Engine, 4096 MB <-- brackets [] signify device being used.
-- -1- APPLE: Intel(R) UHD Graphics 630, 1536 MB
--

-- $visualization
-- The ArrayFire API is able to display visualizations using the Forge library
-- >>> window <- createWindow 800 600 "Histogram"
--

max_stars_repo_name: arrayfire/arrayfire-haskell

max_issues_repo_name: arrayfire/arrayfire-haskell

max_forks_repo_name: arrayfire/arrayfire-haskell


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE ScopedTypeVariables #-}
#if __GLASGOW_HASKELL__ >= 805
{-# LANGUAGE ExplicitNamespaces #-}
{-# LANGUAGE NoStarIsType #-}
#endif
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}

module Eigen.Matrix
  (
    -- * Types
    Matrix(..)
  , Vec(..)
  , MatrixXf
  , MatrixXd
  , MatrixXcf
  , MatrixXcd

    -- * Common API
  , Elem
  , C
  , natToInt
  , Row(..)
  , Col(..)

    -- * Encode/Decode a Matrix
  , encode
  , decode

    -- * Querying a Matrix
  , null
  , square
  , rows
  , cols
  , dims

    -- * Constructing a Matrix
  , empty
  , constant
  , zero
  , ones
  , identity
  , random
  , diagonal

  , (!)
  , coeff
  , generate
  , sum
  , prod
  , mean
  , trace
  , all
  , any
  , count
  , norm
  , squaredNorm
  , blueNorm
  , hypotNorm
  , determinant
  , add
  , sub
  , mul
  , map
  , imap
  , TriangularMode(..)
  , triangularView
  , filter
  , ifilter
  , length
  , foldl
  , foldl'
  , inverse
  , adjoint
  , transpose
  , conjugate
  , normalize
  , modify
  , block
  , unsafeFreeze
  , unsafeWith
  , fromList
  , toList
  ) where

import Control.Monad (when)
import Control.Monad.ST (ST)
import Prelude hiding
  (map, null, filter, length, foldl, any, all, sum)
import Control.Monad (forM_)
import Control.Monad.Primitive (PrimMonad(..))
import Data.Binary (Binary(..))
import qualified Data.Binary as Binary
import qualified Data.ByteString.Lazy as BSL
import Data.Complex (Complex)
import Data.Constraint.Nat
import Eigen.Internal
  ( Elem
  , Cast(..)
  , natToInt
  , Row(..)
  , Col(..)
  )
import qualified Eigen.Internal as Internal
import qualified Eigen.Matrix.Mutable as M
import qualified Data.List as List
import Data.Kind (Type)
import GHC.TypeLits (Nat, type (*), type (<=), KnownNat)
import Foreign.C.Types (CInt)
import Foreign.C.String (CString)
import Foreign.Marshal.Alloc (alloca)
import Foreign.Ptr (Ptr)
import Foreign.Storable (peek)

import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM

-- | Matrix to be used in pure computations.
--
--   * Uses column majour memory layout.
--
--   * Has a copy-free FFI using the <http://eigen.tuxfamily.org Eigen> library.
--
newtype Matrix :: Nat -> Nat -> Type -> Type where
  Matrix :: Vec (n * m) a -> Matrix n m a

-- | Used internally to track the size and corresponding C type of the matrix.
newtype Vec :: Nat -> Type -> Type where
  Vec :: VS.Vector (C a) -> Vec n a

instance forall n m a. (Elem a, Show a, KnownNat n, KnownNat m) => Show (Matrix n m a) where
  show m = List.concat
    [ "Matrix ", show (rows m), "x", show (cols m)
    , "\n", List.intercalate "\n" $ List.map (List.intercalate "\t" . List.map show) $ toList m, "\n"
    ]

instance forall n m a. (KnownNat n, KnownNat m, Elem a) => Binary (Matrix n m a) where
  put (Matrix (Vec vals)) = do
    put $ Internal.magicCode (undefined :: C a)
    put $ natToInt @n
    put $ natToInt @m
    put vals

  get = do
    get >>= (`when` fail "wrong matrix type") . (/= Internal.magicCode (undefined :: C a))
    Matrix . Vec <$> get

-- | Encode the sparse matrix as a lazy bytestring
encode :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> BSL.ByteString
encode = Binary.encode

-- | Decode the sparse matrix from a lazy bytestring
decode :: (Elem a, KnownNat n, KnownNat m) => BSL.ByteString -> Matrix n m a
decode = Binary.decode

-- | Alias for single precision matrix
type MatrixXf n m = Matrix n m Float
-- | Alias for double precision matrix
type MatrixXd n m = Matrix n m Double
-- | Alias for single precision matrix of complex numbers
type MatrixXcf n m = Matrix n m (Complex Float)
-- | Alias for double precision matrix of complex numbers
type MatrixXcd n m = Matrix n m (Complex Double)

-- | Construct an empty 0x0 matrix
empty :: Elem a => Matrix 0 0 a
{-# INLINE empty #-}
empty = Matrix (Vec (VS.empty))

-- | Is matrix empty?
null :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Bool
{-# INLINE null #-}
null m = cols m == 0 && rows m == 0

-- | Is matrix square?
--
square :: forall n m a. (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Bool
{-# INLINE square #-}
square _ = natToInt @n == natToInt @m

-- | Matrix where all coeffs are filled with the given value
constant :: forall n m a. (Elem a, KnownNat n, KnownNat m) => a -> Matrix n m a
{-# INLINE constant #-}
constant !val =
  let !cval = toC val
  in withDims $ \rs cs -> VS.replicate (rs * cs) cval

-- | Matrix where all coeffs are filled with 0
zero :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a
{-# INLINE zero #-}
zero = constant 0

-- | Matrix where all coeffs are filled with 1
ones :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a
{-# INLINE ones #-}
ones = constant 1

-- | The identity matrix (not necessarily square)
identity :: forall n m a. (Elem a, KnownNat n, KnownNat m) => Matrix n m a
identity =
  Internal.performIO $ do
     m :: M.IOMatrix n m a <- M.new
     Internal.call $ M.unsafeWith m Internal.identity
     unsafeFreeze m

-- | The random matrix of a given size
random :: forall n m a. (Elem a, KnownNat n, KnownNat m) => IO (Matrix n m a)
random = do
  m :: M.IOMatrix n m a <- M.new
  Internal.call $ M.unsafeWith m Internal.random
  unsafeFreeze m

withDims :: forall n m a. (Elem a, KnownNat n, KnownNat m) => (Int -> Int -> VS.Vector (C a)) -> Matrix n m a
{-# INLINE withDims #-}
withDims f =
  let !r = natToInt @n
      !c = natToInt @m
  in Matrix $ Vec $ f r c

-- | The number of rows in the matrix
rows :: forall n m a. KnownNat n => Matrix n m a -> Int
{-# INLINE rows #-}
rows _ = natToInt @n

-- | The number of colums in the matrix
cols :: forall n m a. KnownNat m => Matrix n m a -> Int
{-# INLINE cols #-}
cols _ = natToInt @m

-- | Return Matrix size as a pair of (rows, cols)
dims :: forall n m a. (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> (Int, Int)
{-# INLINE dims #-}
dims _ = (natToInt @n, natToInt @m)

-- | Return the value at the given position.
(!) :: forall n m a r c. (Elem a, KnownNat n, KnownNat r, KnownNat c, r <= n, c <= m) => Row r -> Col c -> Matrix n m a -> a
{-# INLINE (!) #-}
(!) = coeff

-- | Return the value at the given position.
coeff :: forall n m a r c. (Elem a, KnownNat n, KnownNat r, KnownNat c, r <= n, c <= m) => Row r -> Col c -> Matrix n m a -> a
{-# INLINE coeff #-}
coeff _ _ m@(Matrix (Vec vals)) =
  let !row  = natToInt @r
      !col  = natToInt @c
  in fromC $! VS.unsafeIndex vals $! col * rows m + row

unsafeCoeff :: (Elem a, KnownNat n) => Int -> Int -> Matrix n m a -> a
{-# INLINE unsafeCoeff #-}
unsafeCoeff row col m@(Matrix (Vec vals)) = fromC $! VS.unsafeIndex vals $! col * rows m + row

-- | Given a generation function `f :: Int -> Int -> a`, construct a Matrix of known size
--   using points in the matrix as inputs.
generate :: forall n m a. (Elem a, KnownNat n, KnownNat m) => (Int -> Int -> a) -> Matrix n m a
generate f = withDims $ \rs cs -> VS.create $ do
  vals :: VSM.MVector s (C a) <- VSM.new (rs * cs)
  forM_ [0 .. pred rs] $ \r ->
    forM_ [0 .. pred cs] $ \c ->
      VSM.write vals (c * rs + r) (toC $! f r c)
  pure vals

-- | The sum of all coefficients in the matrix
sum :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> a
sum = _prop Internal.sum

-- | The product of all coefficients in the matrix
prod :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> a
prod = _prop Internal.prod

-- | The arithmetic mean of all coefficients in the matrix
mean :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> a
mean = _prop Internal.mean

-- | The trace of a matrix is the sum of the diagonal coefficients.
--
--   'trace' m == 'sum' ('diagonal' m)
trace :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> a
trace = _prop Internal.trace

-- | Given a predicate p, determine if all values in the Matrix satisfy p.
all :: (Elem a, KnownNat n, KnownNat m) => (a -> Bool) -> Matrix n m a -> Bool
all f (Matrix (Vec vals)) = VS.all (f . fromC) vals

-- | Given a predicate p, determine if any values in the Matrix satisfy p.
any :: (Elem a, KnownNat n, KnownNat m) => (a -> Bool) -> Matrix n m a -> Bool
any f (Matrix (Vec vals)) = VS.any (f . fromC) vals

-- | Given a predicate p, determine how many values in the Matrix satisfy p.
count :: (Elem a, KnownNat n, KnownNat m) => (a -> Bool) -> Matrix n m a -> Int
count f (Matrix (Vec vals)) = VS.foldl' (\n x-> if f (fromC x) then (n + 1) else n) 0 vals

norm, squaredNorm, blueNorm, hypotNorm :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> a

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

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

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

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

-- | The determinant of the matrix
determinant :: forall n a. (Elem a, KnownNat n) => Matrix n n a -> a
determinant m = _prop Internal.determinant m

-- | Add two matrices.
add :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix n m a -> Matrix n m a
add m1 m2 = _binop Internal.add m1 m2

-- | Subtract two matrices.
sub :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix n m a -> Matrix n m a
sub m1 m2 = _binop Internal.sub m1 m2

-- | Multiply two matrices.
mul :: (Elem a, KnownNat p, KnownNat q, KnownNat r) => Matrix p q a -> Matrix q r a -> Matrix p r a
mul m1 m2 = _binop Internal.mul m1 m2

{- | Apply a given function to each element of the matrix.
Here is an example how to implement scalar matrix multiplication:

>>> let a = fromList [[1,2],[3,4]] :: MatrixXf 2 2
>>> a
Matrix 2x2
1.0 2.0
3.0 4.0
>>> map (*10) a
Matrix 2x2
10.0    20.0
30.0    40.0
-}
map :: Elem a => (a -> a) -> Matrix n m a -> Matrix n m a
map f (Matrix (Vec vals)) = Matrix $ Vec $ VS.map (toC . f . fromC) vals

{- | Apply a given function to each element of the matrix.
Here is an example how upper triangular matrix can be implemented:
>>> let a = fromList [[1,2,3],[4,5,6],[7,8,9]] :: MatrixXf
>>> a
Matrix 3x3
1.0 2.0 3.0
4.0 5.0 6.0
7.0 8.0 9.0
>>> imap (\row col val -> if row <= col then val else 0) a
Matrix 3x3
1.0 2.0 3.0
0.0 5.0 6.0
0.0 0.0 9.0
-}
imap :: (Elem a, KnownNat n, KnownNat m) => (Int -> Int -> a -> a) -> Matrix n m a -> Matrix n m a
imap f (Matrix (Vec vals)) =
  withDims $ \rs _ ->
    VS.imap (\n ->
      let (c,r) = divMod n rs
      in toC . f r c . fromC) vals

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

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

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

-- | Filter elements in the matrix with an indexed predicate. Filtered elements will be replaces by 0.
ifilter :: (Elem a, KnownNat n, KnownNat m) => (Int -> Int -> a -> Bool) -> Matrix n m a -> Matrix n m a
ifilter f = imap (\r c x -> if f r c x then x else 0)

-- | The length of the matrix.
length :: forall n m a r. (Elem a, KnownNat n, KnownNat m, r ~ (n * m), KnownNat r) => Matrix n m a -> Int
length _ = natToInt @r

-- | Left fold of a matrix, where accumulation is lazy.
foldl :: (Elem a, KnownNat n, KnownNat m) => (b -> a -> b) -> b -> Matrix n m a -> b
foldl f b (Matrix (Vec vals)) = VS.foldl (\a x -> f a (fromC x)) b vals

-- | Right fold of a matrix, where accumulation is strict.
foldl' :: Elem a => (b -> a -> b) -> b -> Matrix n m a -> b
foldl' f b (Matrix (Vec vals)) = VS.foldl' (\ !a x -> f a (fromC x)) b vals

-- | Return the diagonal of a matrix.
diagonal :: (Elem a, KnownNat n, KnownNat m, r ~ Min n m, KnownNat r) => Matrix n m a -> Matrix r 1 a
diagonal = _unop Internal.diagonal

{- | Inverse of the matrix
For small fixed sizes up to 4x4, this method uses cofactors. In the general case, this method uses PartialPivLU decomposition
-}
inverse :: forall n a. (Elem a, KnownNat n) => Matrix n n a -> Matrix n n a
inverse = _unop Internal.inverse

-- | Adjoint of the matrix
adjoint :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix m n a
adjoint = _unop Internal.adjoint

-- | Transpose of the matrix
transpose :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix m n a
transpose = _unop Internal.transpose

-- | Conjugate of the matrix
conjugate :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix n m a
conjugate = _unop Internal.conjugate

-- | Normalise the matrix by dividing it on its 'norm'
normalize :: forall n m a. (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix n m a
normalize (Matrix (Vec vals)) = Internal.performIO $ do
  vals' <- VS.thaw vals
  VSM.unsafeWith vals' $ \p ->
    let !rs = natToInt @n
        !cs = natToInt @m
    in Internal.call $ Internal.normalize p (toC rs) (toC cs)
  Matrix . Vec <$> VS.unsafeFreeze vals'

-- | Apply a destructive operation to a matrix. The operation will be performed in-place, if it is safe
--   to do so - otherwise, it will create a copy of the matrix.
modify :: (Elem a, KnownNat n, KnownNat m) => (forall s. M.MMatrix n m s a -> ST s ()) -> Matrix n m a -> Matrix n m a
modify f (Matrix (Vec vals)) = Matrix $ Vec $ VS.modify (f . M.fromVector ) vals

-- | Extract rectangular block from matrix defined by startRow startCol blockRows blockCols
block :: forall sr sc br bc n m a.
     (Elem a, KnownNat sr, KnownNat sc, KnownNat br, KnownNat bc, KnownNat n, KnownNat m)
  => (sr <= n, sc <= m, br <= n, bc <= m)
  => Row sr -- ^ starting row
  -> Col sc -- ^ starting col
  -> Row br -- ^ block of rows
  -> Col bc -- ^ block of cols
  -> Matrix n m a -- ^ extract from this
  -> Matrix br bc a -- ^ extraction
block _ _ _ _ m =
  let !startRow = natToInt @sr
      !startCol = natToInt @sc
  in generate $ \row col -> unsafeCoeff (startRow + row) (startCol + col) m

-- | Turn a mutable matrix into an immutable matrix without copying.
--   The mutable matrix should not be modified after this conversion.
unsafeFreeze :: (Elem a, KnownNat n, KnownNat m, PrimMonad p) => M.MMatrix n m (PrimState p) a -> p (Matrix n m a)
unsafeFreeze m = VS.unsafeFreeze (M.vals m) >>= pure . Matrix . Vec

-- | Pass a pointer to the matrix's data to the IO action. The data may not be modified through the pointer.
unsafeWith  :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> (Ptr (C a) -> CInt -> CInt -> IO b) -> IO b
unsafeWith m@(Matrix (Vec (vals))) f =
  VS.unsafeWith vals $ \p ->
    let !rs = toC $! rows m
        !cs = toC $! cols m
    in f p rs cs

_prop :: (Elem a, KnownNat n, KnownNat m) => (Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString) -> Matrix n m a -> a
{-# INLINE _prop #-}
_prop f m = fromC $ Internal.performIO $ alloca $ \p -> do
   Internal.call $ unsafeWith m (f p)
   peek p

_binop :: forall n m n1 m1 n2 m2 a. (Elem a, KnownNat n, KnownNat m, KnownNat n1, KnownNat m1, KnownNat n2, KnownNat m2)
  => (Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString)
  -> Matrix n m a
  -> Matrix n1 m1 a
  -> Matrix n2 m2 a
{-# INLINE _binop #-}
_binop g m1 m2 = Internal.performIO $ do
  m0 :: M.IOMatrix n2 m2 a <- M.new
  M.unsafeWith m0 $ \vals0 rows0 cols0 ->
      unsafeWith m1 $ \vals1 rows1 cols1 ->
          unsafeWith m2 $ \vals2 rows2 cols2 ->
              Internal.call $ g
                vals0 rows0 cols0
                vals1 rows1 cols1
                vals2 rows2 cols2
  unsafeFreeze m0

_unop :: forall n m n1 m1 a. (Elem a, KnownNat n, KnownNat m, KnownNat n1, KnownNat m1)
  => (Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString)
  -> Matrix n m a
  -> Matrix n1 m1 a
{-# INLINE _unop #-}
_unop g m1 = Internal.performIO $ do
  m0 :: M.IOMatrix n1 m1 a <- M.new
  M.unsafeWith m0 $ \vals0 rows0 cols0 ->
      unsafeWith m1 $ \vals1 rows1 cols1 ->
          Internal.call $ g
              vals0 rows0 cols0
              vals1 rows1 cols1
  unsafeFreeze m0

-- | Convert a matrix to a list.
toList :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> [[a]]
{-# INLINE toList #-}
toList m@(Matrix (Vec vals))
  | null m = []
  | otherwise = [[fromC $ vals `VS.unsafeIndex` (col * _rows + row) | col <- [0..pred _cols]] | row <- [0..pred _rows]]
  where
    !_rows = rows m
    !_cols = cols m

-- | Convert a list to a matrix. Returns 'Nothing' if the dimensions of the list do not match that
--   of the matrix.
fromList :: forall n m a. (Elem a, KnownNat n, KnownNat m) => [[a]] -> Maybe (Matrix n m a)
fromList list = do
  let myRows = natToInt @n
  let myCols = natToInt @m
  let _rows  = List.length list
  let _cols  = List.foldl' max 0 (List.map List.length list)
  if ((myRows /= _rows) || (myCols /= _cols))
    then Nothing
    else (Just . Matrix . Vec) $ VS.create $ do
      vm <- VSM.replicate (_rows * _cols) (toC (0 :: a))
      forM_ (zip [0..] list) $ \(row,vals) ->
        forM_ (zip [0..] vals) $ \(col, val) ->
          VSM.write vm (col * _rows + row) (toC val)
      pure vm

max_stars_repo_name: nilsalex/eigen

max_issues_repo_name: nilsalex/eigen

max_forks_repo_name: nilsalex/eigen


text:
module Mandelbrot where

import Data.List
import Data.Complex
import Data.Maybe

type Point = Complex Double

iterations = 1000

mandelbrot :: Point -> Point -> Point
mandelbrot c z = z ^^ 2 + c

series :: Point -> [Point]
series c = iterate (mandelbrot c) (0 :+ 0)

period :: Point -> Maybe Int
period c = helper 0 (take iterations (series c)) [] where
  helper count (z:zs) memory = if z `elem` memory
    then z `elemIndex` memory >>= (\i -> Just (i + 1))
    else helper (count + 1) zs (z:memory)
  helper count [] memory = Nothing

escaped :: Point -> Bool
escaped c = 4 < (realPart c * realPart c + imagPart c * imagPart c)

-- Does the series at 'c' eventually escape?
divergent :: Point -> Bool
divergent c = any escaped (take iterations (series c))

-- How many iterations does it take the series at 'c' to escape?
divergence :: Point -> Int
divergence c = helper 0 (0.0 :+ 0.0) where
  helper n acc = if n < iterations
    then
      let c' = mandelbrot c acc in
        if escaped c' then n else helper (n + 1) c'
    else
      n

max_stars_repo_name: michaelmp/mandelbrot

max_issues_repo_name: michaelmp/mandelbrot

max_forks_repo_name: michaelmp/mandelbrot


text:
-- H24.12.01 16:43 Eval Eval Eval!
module Main where
import System.Environment
import Text.ParserCombinators.Parsec hiding (spaces)
import Control.Monad
import Numeric
import Data.Ratio
import Data.Complex
import Data.Array
import Control.Monad.Error

main :: IO ()
main = do
  args <- getArgs
  result <- return $ liftM show $ readExpr (args !! 0) >>= eval
  putStrLn $ extractValue $ trapError result

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

spaces :: Parser ()
spaces = skipMany1 space

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

parseAtom :: Parser LispVal
parseAtom = do first <- letter <|> symbol1
               rest <- many (letter <|> digit <|> symbol)
               return . Atom $ first : rest

parseBool :: Parser LispVal
parseBool = liftM Bool ((char 't' >> return True)
                         <|> (char 'f' >> return False))

parseSigned :: (Real r) => ReadS r -> Parser [Char] -> Parser r
parseSigned reader lexer = liftM (fst . head . readSigned reader) (do
  sign <- char '-' <|> return '0'
  num <- lexer
  return (sign : num))

lexFloat :: Parser [Char]
lexFloat = do x <- many1 digit
              char '.'
              y <- many digit
              return $ x ++ "." ++ y

parseFloat :: Parser LispVal
parseFloat = (liftM Float $ parseSigned readFloat lexFloat)

parseDigits :: Parser LispVal
parseDigits = (liftM Number $ parseSigned readDec (many1 digit))

parseRatio :: Parser LispVal
parseRatio = do
  x <- parseSigned readDec (many1 digit)
  char '/'
  y <- parseSigned readDec (many1 digit)
  return $ Ratio (x % y)

toDouble :: LispVal -> Double
toDouble (Number n) = fromIntegral n
toDouble (Float f) = f

parseComplex :: Parser LispVal
parseComplex = do x <- (try parseFloat <|> parseDigits)
                  char '+'
                  y <- (try parseFloat <|> parseDigits)
                  char 'i'
                  return $ Complex (toDouble x :+ toDouble y)

parseNum :: Parser LispVal
parseNum = try parseComplex
           <|> try parseRatio
           <|> try parseFloat
           <|> parseDigits

parseNumber :: Parser LispVal
parseNumber = liftM Number
                (parseBin <|> parseOct <|> parseDig <|> parseHex)
  where
        parseBin = char 'b' >> parseSigned readBin (many1 $ oneOf "01")
        parseOct = char 'o' >> parseSigned readOct (many1 octDigit)
        parseDig = char 'd' >> parseSigned readDec (many1 digit)
        parseHex = char 'x' >> parseSigned readHex (many1 hexDigit)
        readBin  = readInt 2 (\x->x=='0'||x=='1') (\x->if x=='0' then 0 else 1)

parseSharpSyntax :: Parser LispVal
parseSharpSyntax = char '#' >> (parseBool
                     <|> parseNumber
                     <|> parseCharacter
                     <|> parseVector)

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

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

parseCharacter :: Parser LispVal
parseCharacter = liftM Character (char '\\' >>
                   option ' ' (try (string "space" >> return ' ')
                                <|> try (string "newline" >> return '\n')
                                <|> do { x <- anyChar;
                                         notFollowedBy alphaNum;
                                         return x}))

parseList :: Parser LispVal
parseList = do
  body <- sepEndBy parseExpr spaces
  if length body == 0 then return $ List []
  else parseTail >>= \tail-> return $
    case tail of
      Just (List l)         -> List (body ++ l)
      Just (DottedList h t) -> DottedList (body ++ h) t
      Just val              -> DottedList body val
      Nothing               -> List body

parseTail :: Parser (Maybe LispVal)
parseTail = liftM Just (char '.' >> skipMany space >> parseExpr)
            <|> return Nothing

makeSyntaxSugarParser :: Parser a -> String -> Parser LispVal
makeSyntaxSugarParser ch symbol =
  ch >> parseExpr >>= (\x -> return $ List [Atom symbol, x])

parseQuoted :: Parser LispVal
parseQuoted = makeSyntaxSugarParser (char '\'') "quote"

parseQusaiquote :: Parser LispVal
parseQusaiquote = makeSyntaxSugarParser (char '`') "qusaiquote"

parseUnquote :: Parser LispVal
parseUnquote = makeSyntaxSugarParser (char ',') "unquote"

parseUnquoteS :: Parser LispVal
parseUnquoteS = makeSyntaxSugarParser (string ",@") "unquote-splicing"

parens :: Parser a -> Parser a
parens = between (char '(' >> skipMany space) (char ')')

parseVector :: Parser LispVal
parseVector = do
  arrayValues <- parens (sepBy parseExpr spaces)
  return $ Vector (listArray (0,(length arrayValues - 1)) arrayValues)

parseExpr :: Parser LispVal
parseExpr = parseSharpSyntax
        <|> try parseNum
        <|> parseAtom
        <|> parseString
        <|> parseQuoted <|> parseQusaiquote
        <|> try parseUnquoteS <|> parseUnquote
        <|> parens parseList

readExpr :: String -> ThrowsError LispVal
readExpr input = case parse parseExpr "lisp" input of
                   Left err -> throwError $ Parser err
                   Right val -> return val

showVal :: LispVal -> String
showVal (String contents) = "\"" ++ contents ++ "\""
showVal (Atom name) = name
showVal (Number contents) = show contents
showVal (Float contents) = showFFloat Nothing contents ""
showVal (Ratio contents) = show (numerator contents) ++ "/"
                         ++ show (denominator contents)
showVal (Complex x) = showFFloat Nothing (realPart x) $ "+"
                    ++ showFFloat Nothing (imagPart x) "i"
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (Character '\n') = "#\\newline"
showVal (Character ' ') = "#\\space"
showVal (Character char) = "#\\" ++ [char]
showVal (List contents) = "(" ++ unwordsList contents ++ ")"
showVal (DottedList head tail) = "(" ++ unwordsList head
                               ++ " . " ++ showVal tail ++ ")"
showVal (Vector contents) = "#(" ++ unwordsList (elems contents) ++ ")"
unwordsList :: [LispVal] -> String
unwordsList = unwords . map showVal

instance Show LispVal where show = showVal

eval :: LispVal -> ThrowsError LispVal
eval val@(String _) = return val
eval val@(Number _) = return val
eval val@(Bool _) = return val
eval (List [Atom "quote", val]) = return val
eval (List (Atom func : args)) = mapM eval args >>= apply func
eval badForm = throwError $ BadSpecialForm "Unrecognized special form" badForm
--eval v= String$show v

apply :: String -> [LispVal] -> ThrowsError LispVal
apply func args = maybe (throwError $
                      NotFunction "Unrecognized primitive function args" func)
                    ($ args)
                    (lookup func primitives)

primitives :: [(String, [LispVal] -> ThrowsError LispVal)]
primitives = [ ("+", numericBinop (+)),
               ("-", numericBinop (-)),
               ("*", numericBinop (*)),
               ("/", numericBinop div),
               ("mod", numericBinop mod),
               ("quotient", numericBinop quot),
               ("remainder", numericBinop rem),
               ("boolean?", unaryOp booleanp),
               ("symbol?", unaryOp symbolp),
{-               ("char?", unaryOp charp),
               ("vector?", unaryOp vectorp),
--               ("procedure?", unaryOp procedurep),-}
               ("pair?", unaryOp pairp),
               ("number?", unaryOp numberp),{-
               ("complex?", unaryOp complexp),
               ("real?", unaryOp realp),
               ("rational?", unaryOp rationalp),
               ("integer?", unaryOp integerp),-}
               ("string?", unaryOp stringp),
--               ("port?", unaryOp portp),-}
               ("symbol->string", unaryOp symbol2string),
               ("string->symbol", unaryOp string2symbol)
             ]

numericBinop :: (Integer -> Integer -> Integer) -> [LispVal]
                -> ThrowsError LispVal
numericBinop op singleVal@[_] = throwError $ NumArgs 2 singleVal
numericBinop op params = mapM unpackNum params >>= return . Number . foldl1 op

unpackNum :: LispVal -> ThrowsError Integer
unpackNum (Number n) = return n
unpackNum nan = throwError $ TypeMismatch "number" nan

unaryOp :: (LispVal -> LispVal) -> [LispVal] -> ThrowsError LispVal
unaryOp f [v] = return $ f v

symbolp, numberp, stringp, booleanp, pairp :: LispVal -> LispVal
symbolp (Atom _)   = Bool True
symbolp _          = Bool False
numberp (Number _) = Bool True
numberp _          = Bool False
stringp (String _) = Bool True
stringp _          = Bool False
booleanp (Bool _)  = Bool True
booleanp   _       = Bool False
pairp (List _)     = Bool True
pairp (DottedList _ _) = Bool True
pairp _            = Bool False

symbol2string, string2symbol :: LispVal -> LispVal
symbol2string (Atom s)   = String s
symbol2string _          = String ""
string2symbol (String s) = Atom s
string2symbol _          = Atom ""
--booleanp, symbolp, charp, vectorp, {-procedurep,-}
--  pairp, numberp, complexp, realp, rationalp, integerp, stringp{-, portp-} :: LispVal -> Bool
--booleanp

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

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

instance Show LispError where show = showError
instance Error LispError where
  noMsg = Default "An error has occurred"
  strMsg = Default
type ThrowsError = Either LispError

--trapError ::
trapError action = catchError action (return . show)

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

-- ghc -package parsec -o 4errorcheck.exe --make 4errorcheck.hs
-- 4errorcheck "(+ 2 \"two\")"
-- 4errorcheck "(+ 2)"
-- 4errorcheck "(what? 2)"


text:
-- |
-- Module      : main
-- Description : Graphics module leveraging the Gloss package
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental
-- Portability : POSIX (not sure)
--
-- Jonatan H Sundqvist
-- June 2 2015
--

-- TODO | -
--        -

-- SPEC | -
--        -

module GlossGraphics where

---------------------------------------------------------------------------------------------------
-- We'll need these
---------------------------------------------------------------------------------------------------
-- import Graphics.Gloss.Data.Picture (line)
import qualified Graphics.Gloss as Gloss --(circleSolid, rectangleSolid)
import Graphics.Gloss.Interface.IO.Game hiding (Vector)
-- import Graphics.Gloss.Geometry.Angle (degToRad, radToDeg, normaliseAngle)
-- import Graphics.Gloss.Data.Vector

import Data.Complex

import Copernicus

---------------------------------------------------------------------------------------------------
-- Types
---------------------------------------------------------------------------------------------------
data World = World { bodies :: [Body Float],
                     grid   :: Bool
                   }

---------------------------------------------------------------------------------------------------
-- Interaction
---------------------------------------------------------------------------------------------------
simulate :: IO ()
simulate = playIO
	display 	-- Window mode
	white 		-- Background colour
	60			-- FPS (simulation steps per second, technically)
	world 		-- Initial world
	render		-- Converts world to Picture
	respond 	-- User interaction
	advance 	-- Advances the world to the next simulation step
	where
		display = InWindow "Simulator" (width, height) (25, 25)
		world   = World { bodies = map (\ (p', v', g') -> Body p' v' g') [(0.0:+0.0, v, g), (10.0:+0.0, (-20.0):+15.0, g), (35.0:+(-28.0), v, g), (20.5:+19.2, v, g)],
	                      grid   = True
	              }
		-- TODO: Refactor this ugly mess
		render w 	 		  = return . pictures $ zipWith drawBall colours (bodies w) ++ [drawGround (bodies w)] ++ if grid w then [renderGrid 45 45 width height] else []
		drawBall col (Body (x:+y) _ _) = color col . translate x y $ circleSolid 15 												-- TODO: Reorder arguments (?)
		drawGround _ 		  = translate 0 (30/2-fromIntegral height/2) . color green $ Gloss.rectangleSolid (fromIntegral width) 30
		respond e w           = return $ case e of
			EventKey (Char 'g') Down _ _ -> w { grid = not $ grid w }
			_                            -> w

		advance t w 		  = return $ w { bodies=map (animate t) $ bodies w }
		v 					  = 40.0:+20.0 --40.0 :+ 20.0
		g 					  = 0.0:+(-98.2) --0.0 :+ (-9.82)
		(width, height) 	  = (740, 540)
		colours = cycle [red, green, orange, makeColor 0.2 0.1 0.3 1.0]
		--parabola t v a p = let (px:+py) = p; (vx:+vy) = v; (ax:+ay) = a in (px + vx*t + 0.5*ax*t**2):+(py + vy*t + 0.5*ay*t**2)

--
--renderForces

--
-- Subgrids, colours, thickness, markings, hover for coordinates, snap-to-grid
renderGrid :: Float -> Float -> Int -> Int -> Picture
renderGrid dx dy w' h' = pictures $ rows ++ cols
	where
		cols = map (\cl -> line [(cl*dx-w/2, -h/2), (cl*dx-w/2, h/2)]) [1..cls]
		rows = map (\rw -> line [(-w/2, rw*dy-h/2), (w/2, rw*dy-h/2)]) [1..rws]
		rws  = h / dy -- Number of rows
		cls  = w / dx -- Number of columns
		w 	 = fromIntegral w'
		h 	 = fromIntegral h'

-- Transforms a vector from one coordinate space to another
-- by applying the given scaling and translation
-- Useful for converting between simulation and screen coordinates
-- TODO: Make pure (eg. use for pure coordinates so it doesn't depend on Gloss) (?)
transform :: Vector Float -> Vector Float -> Picture -> Picture
transform (sx:+sy) (dx:+dy) = scale sx sy . translate dx dy

max_stars_repo_name: SwiftsNamesake/Copernicus

max_issues_repo_name: SwiftsNamesake/Copernicus

max_forks_repo_name: SwiftsNamesake/Copernicus


text:
module Trainer where

import Numeric.LinearAlgebra
import Models
import Layer
import Optimization
import Data.Foldable
import Control.Monad.IO.Class
import Utils

data TrainConfig = TrainConfig {
  trainConfigRowNum :: Int,
  trainConfigColNum :: Int,
  trainConfigLearningRate :: Double,
  trainConfigMomentum :: Double,
  trainConfigDecay :: Double,
  trainConfigEpsilon :: Double,
  trainConfigBatchSize :: Double,
  trainConfigOptimization:: String
                               }

trainSingleData :: Model VLayer -> TrainConfig -> Matrix R -> IO ()
trainSingleData (Model []) _ _ = print "No models !" >> return ()
trainSingleData neuModel trainConfig inputData =
  let r = trainConfigRowNum trainConfig
      c = trainConfigColNum trainConfig
      paramsList = genParamsList neuModel
      gparamsList = genGParamsList neuModel
      accGradList = accMatrixInit paramsList
      accDeltaList = accMatrixInit paramsList
      optimize = Adadelta {
                adaDecay = trainConfigDecay trainConfig,
                adaEpsilon = trainConfigEpsilon trainConfig,
                adaBatchSize = (trainConfigBatchSize trainConfig),
                adaParams = paramsList,
                adaGparams = gparamsList,
                adaAccGrad = accGradList,
                adaAccDelta = accDeltaList
                            } in
  -- begin to train
  forM_ [1..10000] $ \idx -> do
    --print idx
    let output = forward inputData neuModel in
      --print "backpropagate the model" >>
      let neuModel = backward "mse" output inputData neuModel in
      --print "neumodel" >>
      let gparamsList = updateGParamsList neuModel gparamsList in
      --print "continue" >>
      if mod idx 100 == 0 && idx >= 100 then
        -- update the params in optimize
        -- update the params in model
        let optimize = paramUpdate optimize
            paramsList = adaParams optimize
            neuModel = updateModelParams neuModel paramsList in
        print (getAverageMat output) >>
        return()
      else return()

max_stars_repo_name: neutronest/vortex

max_issues_repo_name: neutronest/vortex

max_forks_repo_name: neutronest/vortex


text:
import Numeric.LinearAlgebra
import Graph
import Adjacency

main :: IO ()
main = do
  putStr . dispf 0 . edges'' . toAdjacency $ completeGraph [1..1000]

max_stars_repo_name: kimolas/causal-haskell

max_issues_repo_name: kimolas/causal-haskell

max_forks_repo_name: kimolas/causal-haskell


text:
module Math.IRT.MLE.Truncated
  ( DF (..)
  , MLEResult (..)
  , mleEst
  ) where

import Data.Default.Class

import Statistics.Distribution

import Math.IRT.Internal.Distribution
import Math.IRT.Internal.LogLikelihood
import Math.IRT.MLE.Internal.Generic

data DF = DF { steps         :: !Int
             , thetaEstimate :: !Double
             , lower_bound   :: !Double
             , upper_bound   :: !Double }

instance Default DF where
  def = DF 10 0.0 (-3.5) 3.5

mleEst :: (Distribution d, ContDistr d, DensityDeriv d, LogLikelihood d) => DF -> [Bool] -> [d] -> MLEResult
mleEst (DF n th lb ub) rs params =
    let res = generic_mleEst rs params n th
    in res { theta = min ub $ max lb $ theta res }

max_stars_repo_name: argiopetech/irt

max_issues_repo_name: argiopetech/irt

max_forks_repo_name: argiopetech/irt


text:
{-# LANGUAGE BangPatterns #-}
{-|
Provides a parser for text formatted MNIST data.
-}
module ML.MnistLoader (parseMnist) where

import           ML.Sample (Sample(..))

import           Data.Attoparsec.ByteString.Char8 (Parser, char, double, parseOnly, sepBy)
import qualified Data.ByteString.Char8 as BS
import qualified Data.Vector as V
import           Numeric.LinearAlgebra (R, assoc, vector)

makeSample :: ([R], [R]) -> Sample
makeSample (image, [expected]) = Sample (vector image) (assoc 10 0 [(floor expected, 1)])

-- | Break the list of samples into components.
--
-- These components are implicit in the file format.
toSampleSet :: V.Vector Sample -> (V.Vector Sample, V.Vector Sample, V.Vector Sample)
toSampleSet vec = (train, validation, test)
    where (train, rest) = V.splitAt 50000 vec
          (validation, test) = V.splitAt 10000 rest

-- | Parse an MNIST sample image.
sample :: Parser Sample
sample = do
    !x <- double `sepBy` (char ' ')
    let !s = makeSample $ splitAt 784 x
    pure s

-- | Parse the given bytestring into (training data, validation data, test data).
--
-- Data format:
--
--   * 70000 lines lines of 785 doubles.
--   * first 1     - 50000 lines are test data.
--   * lines 50001 - 60000 are validation data.
--   * lines 60001 - 70000 are test data.
--   * first 784 double of each line are the image.
--   * last double of each is the expected value.
parseMnist :: BS.ByteString -> (V.Vector Sample, V.Vector Sample, V.Vector Sample)
parseMnist = toSampleSet . V.unfoldr parseLines . BS.lines
    where parseLines :: [BS.ByteString] -> Maybe (Sample, [BS.ByteString])
          parseLines []     = Nothing
          parseLines (l:ls) = Just (fromRight $ parseOnly sample l, ls)

          fromRight (Right b) = b

max_stars_repo_name: m-renaud/ML

max_issues_repo_name: m-renaud/ML

max_forks_repo_name: m-renaud/ML