text:
module Main where

import Svg
import Data.Complex
import Dft
import Data.Foldable
import System.Environment
import Data.Maybe
import Util
import qualified Data.Vector as V
import Numeric.FFT.Vector.Invertible

_N = 256 :: Int
_N_Minus1 = _N - 1 :: Int

square :: (Enum a, Floating a) => Int -> [a]
square n = [sum [sin ((fromIntegral (2 * (k * 2 + 1) * x)) * pi / (fromIntegral _N)) / (fromIntegral (k * 2 + 1)) | k <- [0..n-1]] | x <- [0.._N_Minus1]]

sawtooth = 0 : [1 - (2 / (fromIntegral _N)) * (fromIntegral x) | x <- [1.._N_Minus1]]

approxSawtooth :: (Enum a, Floating a) => Int -> [a]
approxSawtooth n = [foldr1 (+) [sin (fromIntegral (2 * k * x) * pi / (fromIntegral _N)) / (fromIntegral k) | k <- [1..n]] | x <- [0.._N_Minus1]]

signal = approxSawtooth 10

realToComplex :: (Floating a) => a -> Complex a
realToComplex a = a :+ 0

freqDomain :: [Complex Double]
freqDomain = let timeDomain = (V.fromList signal) in
  V.toList $ run dftR2C timeDomain

freqDomainAmp :: [Double]
freqDomainAmp = let (h : t) = freqDomain
                    (m, l) = splitAt (length t - 1) t in
  (magnitude h / 2) : (map magnitude m) ++ (map ((/2) . magnitude) l)

prefix :: Maybe String
prefix = Nothing

main :: IO ()
main = do
  args <- getArgs
  let f = (read . head) args
      graphStroke = Stroke (Just "black") (Just 0.03)
      graphScale = 8 / (fromIntegral _N)
      dasharray = fromMaybe "2" ((show . (*2)) <$> strokeWidthMaybe graphStroke) ++
                  fromMaybe " 1" (((' ' :) . show) <$> strokeWidthMaybe graphStroke)
      helperStroke = strokeWidthMap (/4) graphStroke
      signalGraphTransform = "translate(-4 3)"
  writeFile "background.svg" (
    prologue prefix "xlink" (Just (WidthHeight "600" "600")) (Just $ ViewBox "-5" "-5" "10" "10") ++
        "<" ++ maybePrefix prefix ++ "g transform=\"" ++ signalGraphTransform ++ "\">" ++
        foldr1 (++) (map (\i -> let x = fromIntegral (i * _N) / (fromIntegral f) in
                           "<" ++ maybePrefix prefix ++
                           "line" ++ strokeToString helperStroke ++
                           " stroke-dasharray=\"" ++
                           dasharray ++
                           "\" x1=\"0\" y1=\"1\" x2=\"0\" y2=\"-1\"" ++
                           " transform=\"translate(" ++ show (x * graphScale) ++ " 0)\" />") [0..f]) ++
        "<" ++ maybePrefix prefix ++ "line" ++ strokeToString helperStroke ++
        " stroke-dasharray=\"" ++ dasharray ++ "\" x1=\"" ++ show (-5 * graphScale) ++ "\" y1=\"0\"" ++
        " x2=\"" ++ show ((fromIntegral (_N + 5)) * graphScale) ++ "\" y2=\"0\" " ++
        " />" ++
        drawGraph prefix graphScale 1.0 graphStroke (zip [0..(fromIntegral _N_Minus1)] signal) ++
        "</" ++ maybePrefix prefix ++ "g>" ++
        "<" ++ maybePrefix prefix ++ "g transform=\"translate(-4 -4)\">" ++
        drawGraph prefix (graphScale * 2) freqDomainScale graphStroke (zip [0..(fromIntegral _N_Minus1)] freqDomainAmp) ++
        "<" ++ maybePrefix prefix ++ "circle cx=\"" ++ show ((fromIntegral f) * graphScale * 2) ++ "\" cy=\"" ++
        show ((freqDomainAmp !! f) * freqDomainScale)  ++ "\" r=\"" ++
        fromMaybe "1" (show <$> strokeWidthMaybe graphStroke) ++ "\"" ++
        (fromMaybe "" ((\s -> " fill=\"" ++ s ++ "\"") <$> (strokeColorMaybe graphStroke))) ++ "/>" ++
        "</" ++ maybePrefix prefix ++ "g>" ++
        epilogue prefix
    )
  for_ (zip3 (animateFreq (map realToComplex signal) f (_N `div` f)) [0..] signal)
    (\(content, num, s) ->
      writeFile ("graph" ++ show num ++ ".svg") (
        prologue prefix "xlink" (Just (WidthHeight "600" "600")) (Just $ ViewBox "-5" "-5" "10" "10") ++
        content ++
        "<" ++ maybePrefix prefix ++ "g transform=\"" ++ signalGraphTransform ++ "\">" ++
        "<" ++ maybePrefix prefix ++ "circle cx=\"" ++ show ((fromIntegral num) * graphScale) ++ "\" cy=\"" ++ show s ++ "\" r=\"" ++
        fromMaybe "1" (show <$> strokeWidthMaybe graphStroke) ++ "\"" ++
        (fromMaybe "" ((\s -> " fill=\"" ++ s ++ "\"") <$> (strokeColorMaybe graphStroke))) ++ "/>" ++
        "</" ++ maybePrefix prefix ++ "g>" ++
        epilogue prefix))

freqDomainScale :: RealFloat a => a
freqDomainScale = 1 / (2 * sqrt (fromIntegral (length signal)))

animateFreq :: (RealFloat a, Show a) => [Complex a] -> Int -> Int -> [String]
animateFreq signal f sumTraceN = let
  functionValues = ((map (\(s, wp, sum) -> (s, wp, sum * (freqDomainScale :+ 0)))) . tail . reverse) (correlateSignalAtFreq signal f)
  lastN = lastNWindow sumTraceN functionValues in
  map partToString (zip functionValues lastN) where
    partToString ((s, wavePart, sum), lastNTuples) =
      let opacities = [1.0 - (fromIntegral i) / (fromIntegral sumTraceN)  | i <- [0..sumTraceN - 1]]
          lastComplexes = zip lastNTuples opacities
          lastComplexesDraws =
            map
            (\((_, _, a :+ b), opacity) ->
              "<" ++ maybePrefix prefix ++ "circle fill=\"blue\" cx=\"" ++ show a ++
              "\" cy=\"" ++ show b ++ "\" r=\"0.01\" opacity=\"" ++ show opacity ++
              "\" />"
            )
            lastComplexes
      in
      drawComplex sum prefix (Stroke (Just "blue") (Just 0.03)) ++
      foldr1 (++) lastComplexesDraws ++
      drawComplex wavePart prefix (Stroke (Just "green") (Just 0.02)) ++
      drawComplex s prefix (Stroke (Just "red") (Just 0.01))

max_stars_repo_name: tilgalas/haskell-fourier

max_issues_repo_name: tilgalas/haskell-fourier

max_forks_repo_name: tilgalas/haskell-fourier


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
module Grenade.Layers.Dropout (
    Dropout (..)
  , randomDropout
  , SpecDropout (..)
  , specDropout
  , dropout
  , dropoutWithSeed
  ) where

import           Control.DeepSeq
import           Control.Monad.Primitive        (PrimBase, PrimState)
import           Control.Monad.Primitive        (PrimMonad)
import           Data.Maybe                     (fromJust)
import           Data.Proxy
import           Data.Reflection                (reifyNat)
import           Data.Serialize
import           Data.Singletons
import qualified Data.Vector.Generic            as VG
import qualified Data.Vector.Storable           as V
import qualified Data.Word                      as W
import           GHC.Generics                   hiding (R)
import           GHC.TypeLits
import           Numeric.LinearAlgebra.Static   hiding (Seed)
import           System.Random.MWC              hiding (create)

import           Grenade.Core
import           Grenade.Dynamic
import           Grenade.Dynamic.Internal.Build
import           Grenade.Types
import           Grenade.Utils.Vector

-- | Dropout layer help to reduce overfitting. Idea here is that the vector is a shape of 1s and 0s, which we multiply the input by. After backpropogation, we return a new matrix/vector, with
-- different bits dropped out. The provided argument is the proportion to drop in each training iteration (like 1% or 5% would be reasonable). The proportion is given as promille, i.e. 100 are 10%.
data Dropout (pct :: Nat) =
  Dropout
    { dropoutActive :: Bool     -- ^ Add possibility to deactivate dropout
    , dropoutSeed   :: !Int     -- ^ Seed
    }
  deriving (Generic)

instance NFData (Dropout pct) where rnf (Dropout a s) = rnf a `seq` rnf s
instance Show (Dropout pct) where show (Dropout _ _) = "Dropout"
instance Serialize (Dropout pct) where
  put (Dropout act seed) = put act >> put seed
  get = Dropout <$> get <*> get

instance UpdateLayer (Dropout pct) where
  type Gradient (Dropout pct) = ()
  runUpdate _ (Dropout act seed) _ = Dropout act (seed+1)
  runSettingsUpdate set (Dropout _ seed) = Dropout (setDropoutActive set) seed

instance RandomLayer (Dropout pct) where
  createRandomWith _ = randomDropout

randomDropout :: (PrimBase m) => Gen (PrimState m) -> m (Dropout pct)
randomDropout gen = Dropout True <$> uniform gen

instance (KnownNat pct, KnownNat i) => Layer (Dropout pct) ('D1 i) ('D1 i) where
  type Tape (Dropout pct) ('D1 i) ('D1 i) = V.Vector RealNum
  runForwards (Dropout act seed) (S1D x)
    | not act = (extract v, S1D $ dvmap (rate *) x) -- multily with rate to normalise throughput
    | otherwise = (extract v, S1D $ v * x)
    where
      rate = (/1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy pct)
      v = dvmap mask $ randomVector seed Uniform
      mask r
        | not act || r < rate = 1
        | otherwise = 0
  runForwards (Dropout act seed) (S1DV x)
    | not act = (v, S1DV $ mapVectorInPlace (rate *) x) -- multily with rate to normalise throughput
    | otherwise = (v, S1DV $ zipWithVectorInPlaceSnd (*) v x)
    where
      rate :: RealNum
      rate = (/1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy pct)
      v :: V.Vector RealNum
      v = mapVectorInPlace mask $ extract (randomVector seed Uniform :: R i)
      mask r
        | not act || r < rate = 1
        | otherwise = 0
  runBackwards _ v (S1D x)  = ((), S1D $ x * fromJust (create v))
  runBackwards _ v (S1DV x) = ((), S1DV $ x * v)

-------------------- DynamicNetwork instance --------------------

instance (KnownNat pct) => FromDynamicLayer (Dropout pct) where
  fromDynamicLayer inp _ (Dropout _ seed) = case tripleFromSomeShape inp of
    (rows, 1, 1) -> SpecNetLayer $ SpecDropout rows rate (Just seed)
    _            -> error "Dropout is only allows for vectors, i.e. 1D spaces."
    where rate = (/1000) $ fromIntegral $ max 0 $ min 1000 $ natVal (Proxy :: Proxy pct)

instance ToDynamicLayer SpecDropout where
  toDynamicLayer _ gen (SpecDropout rows rate mSeed) =
    reifyNat rows $ \(_ :: (KnownNat i) => Proxy i) ->
    reifyNat (round $ 1000 * rate) $ \(_ :: (KnownNat pct) => Proxy pct) ->
    case mSeed of
      Just seed -> return $ SpecLayer (Dropout True seed :: Dropout pct) (sing :: Sing ('D1 i)) (sing :: Sing ('D1 i))
      Nothing -> do
        layer <-  randomDropout gen
        return $ SpecLayer (layer :: Dropout pct) (sing :: Sing ('D1 i)) (sing :: Sing ('D1 i))

-- | Create a specification for a droput layer by providing the input size of the vector (1D allowed only!), a rate of nodes to keep (e.g. 0.95) and maybe a seed.
specDropout :: Integer -> RealNum -> Maybe Int -> SpecNet
specDropout i rate seed = SpecNetLayer $ SpecDropout i rate seed

-- | Create a dropout layer with the specified keep rate of nodes. The seed will be randomly initialized when the network is created. See also @dropoutWithSeed@.
dropout :: RealNum -> BuildM ()
dropout ratio = dropoutWithSeed ratio Nothing

-- | Create a dropout layer with the specified keep rate of nodes. The seed will be randomly initialized when the network is created. See also @dropoutWithSeed@.
dropoutWithSeed :: RealNum -> Maybe Int -> BuildM ()
dropoutWithSeed ratio mSeed = buildRequireLastLayerOut Is1D >>= \(i, _, _) -> buildAddSpec (specDropout i ratio mSeed)

-------------------- GNum instance --------------------

instance GNum (Dropout pct) where
  _ |* x = x
  _ |+ x = x

ext: hs

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
{-# LANGUAGE BangPatterns             #-}
{-# LANGUAGE ConstrainedClassMethods  #-}
{-# LANGUAGE DeriveGeneric            #-}
{-# LANGUAGE FlexibleContexts         #-}
{-# LANGUAGE FlexibleInstances        #-}
{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE TypeFamilies             #-}
{-# LANGUAGE TypeOperators            #-}
{-# LANGUAGE ViewPatterns             #-}

-- |
-- Module      :  Internal.Matrix
-- Copyright   :  (c) Alberto Ruiz 2007-15
-- License     :  BSD3
-- Maintainer  :  Alberto Ruiz
-- Stability   :  provisional
--
-- Internal matrix representation
--

module Internal.Matrix where

import           Control.DeepSeq       (NFData (..))
import           Data.Complex          (Complex)
import           Foreign.C.String      (CString, newCString)
import           Foreign.C.Types       (CInt (..))
import           Foreign.Marshal.Alloc (free)
import           Foreign.Marshal.Array (newArray)
import           Foreign.Ptr           (Ptr)
import           Foreign.Storable      (Storable)
import           Internal.Devel
import           Internal.Vector
import           Internal.Vectorized   hiding (( # ), ( #! ))
import           System.IO.Unsafe      (unsafePerformIO)
import           Text.Printf

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

data MatrixOrder = RowMajor | ColumnMajor deriving (Show,Eq)

-- | Matrix representation suitable for BLAS\/LAPACK computations.

data Matrix t = Matrix
    { irows :: {-# UNPACK #-} !Int
    , icols :: {-# UNPACK #-} !Int
    , xRow  :: {-# UNPACK #-} !Int
    , xCol  :: {-# UNPACK #-} !Int
    , xdat  :: {-# UNPACK #-} !(Vector t)
    }

rows :: Matrix t -> Int
rows = irows
{-# INLINE rows #-}

cols :: Matrix t -> Int
cols = icols
{-# INLINE cols #-}

size :: Matrix t -> (Int, Int)
size m = (irows m, icols m)
{-# INLINE size #-}

rowOrder :: Matrix t -> Bool
rowOrder m = xCol m == 1 || cols m == 1
{-# INLINE rowOrder #-}

colOrder :: Matrix t -> Bool
colOrder m = xRow m == 1 || rows m == 1
{-# INLINE colOrder #-}

is1d :: Matrix t -> Bool
is1d (size->(r,c)) = r==1 || c==1
{-# INLINE is1d #-}

-- data is not contiguous
isSlice :: Storable t => Matrix t -> Bool
isSlice m@(size->(r,c)) = r*c < dim (xdat m)
{-# INLINE isSlice #-}

orderOf :: Matrix t -> MatrixOrder
orderOf m = if rowOrder m then RowMajor else ColumnMajor

showInternal :: Storable t => Matrix t -> IO ()
showInternal m = printf "%dx%d %s %s %d:%d (%d)\n" r c slc ord xr xc dv
  where
    r  = rows m
    c  = cols m
    xr = xRow m
    xc = xCol m
    slc = if isSlice m then "slice" else "full"
    ord = if is1d m then "1d" else if rowOrder m then "rows" else "cols"
    dv = dim (xdat m)

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

-- | Matrix transpose.
trans :: Matrix t -> Matrix t
trans m@Matrix { irows = r, icols = c, xRow = xr, xCol = xc } =
             m { irows = c, icols = r, xRow = xc, xCol = xr }

cmat :: (Element t) => Matrix t -> Matrix t
cmat m
    | rowOrder m = m
    | otherwise  = extractAll RowMajor m

fmat :: (Element t) => Matrix t -> Matrix t
fmat m
    | colOrder m = m
    | otherwise  = extractAll ColumnMajor m

-- C-Haskell matrix adapters
{-# INLINE amatr #-}
amatr :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> Ptr a -> f) -> IO r
amatr x f g = unsafeWith (xdat x) (f . g r c)
  where
    r  = fi (rows x)
    c  = fi (cols x)

{-# INLINE amat #-}
amat :: Storable a => Matrix a -> (f -> IO r) -> (CInt -> CInt -> CInt -> CInt -> Ptr a -> f) -> IO r
amat x f g = unsafeWith (xdat x) (f . g r c sr sc)
  where
    r  = fi (rows x)
    c  = fi (cols x)
    sr = fi (xRow x)
    sc = fi (xCol x)

instance Storable t => TransArray (Matrix t)
  where
    type TransRaw (Matrix t) b = CInt -> CInt -> Ptr t -> b
    type Trans (Matrix t) b    = CInt -> CInt -> CInt -> CInt -> Ptr t -> b
    apply = amat
    {-# INLINE apply #-}
    applyRaw = amatr
    {-# INLINE applyRaw #-}

infixr 1 #
(#) :: TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
a # b = apply a b
{-# INLINE (#) #-}

(#!) :: (TransArray c, TransArray c1) => c1 -> c -> Trans c1 (Trans c (IO r)) -> IO r
a #! b = a # b # id
{-# INLINE (#!) #-}

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

copy :: Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy ord m = extractR ord m 0 (idxs[0,rows m-1]) 0 (idxs[0,cols m-1])

extractAll :: Element t => MatrixOrder -> Matrix t -> Matrix t
extractAll ord m = unsafePerformIO (copy ord m)

{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.

>>> flatten (ident 3)
[1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]
it :: (Num t, Element t) => Vector t

-}
flatten :: Element t => Matrix t -> Vector t
flatten m
    | isSlice m || not (rowOrder m) = xdat (extractAll RowMajor m)
    | otherwise                     = xdat m

-- | the inverse of 'Data.Packed.Matrix.fromLists'
toLists :: (Element t) => Matrix t -> [[t]]
toLists = map toList . toRows

-- | common value with \"adaptable\" 1
compatdim :: [Int] -> Maybe Int
compatdim [] = Nothing
compatdim [a] = Just a
compatdim (a:b:xs)
    | a==b = compatdim (b:xs)
    | a==1 = compatdim (b:xs)
    | b==1 = compatdim (a:xs)
    | otherwise = Nothing

-- | Create a matrix from a list of vectors.
-- All vectors must have the same dimension,
-- or dimension 1, which is are automatically expanded.
fromRows :: Element t => [Vector t] -> Matrix t
fromRows [] = emptyM 0 0
fromRows vs = case compatdim (map dim vs) of
    Nothing -> error $ "fromRows expects vectors with equal sizes (or singletons), given: " ++ show (map dim vs)
    Just 0  -> emptyM r 0
    Just c  -> matrixFromVector RowMajor r c . vjoin . map (adapt c) $ vs
  where
    r = length vs
    adapt c v
        | c == 0 = fromList[]
        | dim v == c = v
        | otherwise = constantD (v@>0) c

-- | extracts the rows of a matrix as a list of vectors
toRows :: Element t => Matrix t -> [Vector t]
toRows m
    | rowOrder m = map sub rowRange
    | otherwise  = map ext rowRange
  where
    rowRange = [0..rows m-1]
    sub k = subVector (k*xRow m) (cols m) (xdat m)
    ext k = xdat $ unsafePerformIO $ extractR RowMajor m 1 (idxs[k]) 0 (idxs[0,cols m-1])

-- | Creates a matrix from a list of vectors, as columns
fromColumns :: Element t => [Vector t] -> Matrix t
fromColumns m = trans . fromRows $ m

-- | Creates a list of vectors from the columns of a matrix
toColumns :: Element t => Matrix t -> [Vector t]
toColumns m = toRows . trans $ m

-- | Reads a matrix position.
(@@>) :: Storable t => Matrix t -> (Int,Int) -> t
infixl 9 @@>
m@Matrix {irows = r, icols = c} @@> (i,j)
    | i<0 || i>=r || j<0 || j>=c = error "matrix indexing out of range"
    | otherwise = atM' m i j
{-# INLINE (@@>) #-}

--  Unsafe matrix access without range checking
atM' :: Storable t => Matrix t -> Int -> Int -> t
atM' m i j = xdat m `at'` (i * (xRow m) + j * (xCol m))
{-# INLINE atM' #-}

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

matrixFromVector :: Storable t => MatrixOrder -> Int -> Int -> Vector t -> Matrix t
matrixFromVector _ 1 _ v@(dim->d) = Matrix { irows = 1, icols = d, xdat = v, xRow = d, xCol = 1 }
matrixFromVector _ _ 1 v@(dim->d) = Matrix { irows = d, icols = 1, xdat = v, xRow = 1, xCol = d }
matrixFromVector o r c v
    | r * c == dim v = m
    | otherwise = error $ "can't reshape vector dim = "++ show (dim v)++" to matrix " ++ shSize m
  where
    m | o == RowMajor = Matrix { irows = r, icols = c, xdat = v, xRow = c, xCol = 1 }
      | otherwise     = Matrix { irows = r, icols = c, xdat = v, xRow = 1, xCol = r }

-- allocates memory for a new matrix
createMatrix :: (Storable a) => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix ord r c = do
    p <- createVector (r*c)
    return (matrixFromVector ord r c p)

{- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns. (GNU-Octave groups by columns. To do it you can define @reshapeF r = tr' . reshape r@
where r is the desired number of rows.)

>>> reshape 4 (fromList [1..12])
(3><4)
 [ 1.0,  2.0,  3.0,  4.0
 , 5.0,  6.0,  7.0,  8.0
 , 9.0, 10.0, 11.0, 12.0 ]

-}
reshape :: Storable t => Int -> Vector t -> Matrix t
reshape 0 v = matrixFromVector RowMajor 0 0 v
reshape c v = matrixFromVector RowMajor (dim v `div` c) c v

-- | application of a vector function on the flattened matrix elements
liftMatrix :: (Element a, Element b) => (Vector a -> Vector b) -> Matrix a -> Matrix b
liftMatrix f m@Matrix { irows = r, icols = c, xdat = d}
    | isSlice m = matrixFromVector RowMajor r c (f (flatten m))
    | otherwise = matrixFromVector (orderOf m) r c (f d)

-- | application of a vector function on the flattened matrices elements
liftMatrix2 :: (Element t, Element a, Element b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
liftMatrix2 f m1@(size->(r,c)) m2
    | (r,c)/=size m2 = error "nonconformant matrices in liftMatrix2"
    | rowOrder m1 = matrixFromVector RowMajor    r c (f (flatten m1) (flatten m2))
    | otherwise   = matrixFromVector ColumnMajor r c (f (flatten (trans m1)) (flatten (trans m2)))

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

-- | Supported matrix elements.
class (Storable a) => Element a where
    constantD  :: a -> Int -> Vector a
    extractR :: MatrixOrder -> Matrix a -> CInt -> Vector CInt -> CInt -> Vector CInt -> IO (Matrix a)
    setRect  :: Int -> Int -> Matrix a -> Matrix a -> IO ()
    sortI    :: Ord a => Vector a -> Vector CInt
    sortV    :: Ord a => Vector a -> Vector a
    compareV :: Ord a => Vector a -> Vector a -> Vector CInt
    selectV  :: Vector CInt -> Vector a -> Vector a -> Vector a -> Vector a
    remapM   :: Matrix CInt -> Matrix CInt -> Matrix a -> Matrix a
    rowOp    :: Int -> a -> Int -> Int -> Int -> Int -> Matrix a -> IO ()
    gemm     :: Vector a -> Matrix a -> Matrix a -> Matrix a -> IO ()
    reorderV :: Vector CInt-> Vector CInt-> Vector a -> Vector a -- see reorderVector for documentation

instance Element Float where
    constantD  = constantAux cconstantF
    extractR   = extractAux c_extractF
    setRect    = setRectAux c_setRectF
    sortI      = sortIdxF
    sortV      = sortValF
    compareV   = compareF
    selectV    = selectF
    remapM     = remapF
    rowOp      = rowOpAux c_rowOpF
    gemm       = gemmg c_gemmF
    reorderV   = reorderAux c_reorderF

-- instance Element Double where
--     constantD  = constantAux cconstantR
--     extractR   = extractAux c_extractD
--     setRect    = setRectAux c_setRectD
--     sortI      = sortIdxD
--     sortV      = sortValD
--     compareV   = compareD
--     selectV    = selectD
--     remapM     = remapD
--     rowOp      = rowOpAux c_rowOpD
--     gemm       = gemmg c_gemmD
--     reorderV   = reorderAux c_reorderD

instance Element (Complex Float) where
    constantD  = constantAux cconstantQ
    extractR   = extractAux c_extractQ
    setRect    = setRectAux c_setRectQ
    sortI      = undefined
    sortV      = undefined
    compareV   = undefined
    selectV    = selectQ
    remapM     = remapQ
    rowOp      = rowOpAux c_rowOpQ
    gemm       = gemmg c_gemmQ
    reorderV   = reorderAux c_reorderQ

-- instance Element (Complex Double) where
--     constantD  = constantAux cconstantC
--     extractR   = extractAux c_extractC
--     setRect    = setRectAux c_setRectC
--     sortI      = undefined
--     sortV      = undefined
--     compareV   = undefined
--     selectV    = selectC
--     remapM     = remapC
--     rowOp      = rowOpAux c_rowOpC
--     gemm       = gemmg c_gemmC
--     reorderV   = reorderAux c_reorderC

instance Element (CInt) where
    constantD  = constantAux cconstantI
    extractR   = extractAux c_extractI
    setRect    = setRectAux c_setRectI
    sortI      = sortIdxI
    sortV      = sortValI
    compareV   = compareI
    selectV    = selectI
    remapM     = remapI
    rowOp      = rowOpAux c_rowOpI
    gemm       = gemmg c_gemmI
    reorderV   = reorderAux c_reorderI

instance Element Z where
    constantD  = constantAux cconstantL
    extractR   = extractAux c_extractL
    setRect    = setRectAux c_setRectL
    sortI      = sortIdxL
    sortV      = sortValL
    compareV   = compareL
    selectV    = selectL
    remapM     = remapL
    rowOp      = rowOpAux c_rowOpL
    gemm       = gemmg c_gemmL
    reorderV   = reorderAux c_reorderL

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

-- | reference to a rectangular slice of a matrix (no data copy)
subMatrix :: Element a
            => (Int,Int) -- ^ (r0,c0) starting position
            -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
            -> Matrix a -- ^ input matrix
            -> Matrix a -- ^ result
subMatrix (r0,c0) (rt,ct) m
    | rt <= 0 || ct <= 0 = matrixFromVector RowMajor (max 0 rt) (max 0 ct) (fromList [])
    | 0 <= r0 && 0 <= rt && r0+rt <= rows m &&
      0 <= c0 && 0 <= ct && c0+ct <= cols m = res
    | otherwise = error $ "wrong subMatrix "++show ((r0,c0),(rt,ct))++" of "++shSize m
  where
    p = r0 * xRow m + c0 * xCol m
    tot | rowOrder m = ct + (rt-1) * xRow m
        | otherwise  = rt + (ct-1) * xCol m
    res = m { irows = rt, icols = ct, xdat = subVector p tot (xdat m) }

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

maxZ :: (Num t1, Ord t1, Foldable t) => t t1 -> t1
maxZ xs = if minimum xs == 0 then 0 else maximum xs

conformMs :: Element t => [Matrix t] -> [Matrix t]
conformMs ms = map (conformMTo (r,c)) ms
  where
    r = maxZ (map rows ms)
    c = maxZ (map cols ms)

conformVs :: Element t => [Vector t] -> [Vector t]
conformVs vs = map (conformVTo n) vs
  where
    n = maxZ (map dim vs)

conformMTo :: Element t => (Int, Int) -> Matrix t -> Matrix t
conformMTo (r,c) m
    | size m == (r,c) = m
    | size m == (1,1) = matrixFromVector RowMajor r c (constantD (m@@>(0,0)) (r*c))
    | size m == (r,1) = repCols c m
    | size m == (1,c) = repRows r m
    | otherwise = error $ "matrix " ++ shSize m ++ " cannot be expanded to " ++ shDim (r,c)

conformVTo :: Element t => Int -> Vector t -> Vector t
conformVTo n v
    | dim v == n = v
    | dim v == 1 = constantD (v@>0) n
    | otherwise = error $ "vector of dim=" ++ show (dim v) ++ " cannot be expanded to dim=" ++ show n

repRows :: Element t => Int -> Matrix t -> Matrix t
repRows n x = fromRows (replicate n (flatten x))
repCols :: Element t => Int -> Matrix t -> Matrix t
repCols n x = fromColumns (replicate n (flatten x))

shSize :: Matrix t -> [Char]
shSize = shDim . size

shDim :: (Show a, Show a1) => (a1, a) -> [Char]
shDim (r,c) = "(" ++ show r ++"x"++ show c ++")"

emptyM :: Storable t => Int -> Int -> Matrix t
emptyM r c = matrixFromVector RowMajor r c (fromList[])

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

instance (Storable t, NFData t) => NFData (Matrix t)
  where
    rnf m | d > 0     = rnf (v @> 0)
          | otherwise = ()
      where
        d = dim v
        v = xdat m

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

extractAux :: (Eq t3, Eq t2, TransArray c, Storable a, Storable t1,
                Storable t, Num t3, Num t2, Integral t1, Integral t)
           => (t3 -> t2 -> CInt -> Ptr t1 -> CInt -> Ptr t
                  -> Trans c (CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt))
           -> MatrixOrder -> c -> t3 -> Vector t1 -> t2 -> Vector t -> IO (Matrix a)
extractAux f ord m moder vr modec vc = do
    let nr = if moder == 0 then fromIntegral $ vr@>1 - vr@>0 + 1 else dim vr
        nc = if modec == 0 then fromIntegral $ vc@>1 - vc@>0 + 1 else dim vc
    r <- createMatrix ord nr nc
    (vr # vc # m #! r) (f moder modec)  #|"extract"

    return r

type Extr x = CInt -> CInt -> CIdxs (CIdxs (OM x (OM x (IO CInt))))

-- foreign import ccall unsafe "extractD" c_extractD :: Extr Double
foreign import ccall unsafe "extractF" c_extractF :: Extr Float
-- foreign import ccall unsafe "extractC" c_extractC :: Extr (Complex Double)
foreign import ccall unsafe "extractQ" c_extractQ :: Extr (Complex Float)
foreign import ccall unsafe "extractI" c_extractI :: Extr CInt
foreign import ccall unsafe "extractL" c_extractL :: Extr Z

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

setRectAux :: (TransArray c1, TransArray c)
           => (CInt -> CInt -> Trans c1 (Trans c (IO CInt)))
           -> Int -> Int -> c1 -> c -> IO ()
setRectAux f i j m r = (m #! r) (f (fi i) (fi j)) #|"setRect"

type SetRect x = I -> I -> x ::> x::> Ok

-- foreign import ccall unsafe "setRectD" c_setRectD :: SetRect Double
foreign import ccall unsafe "setRectF" c_setRectF :: SetRect Float
-- foreign import ccall unsafe "setRectC" c_setRectC :: SetRect (Complex Double)
foreign import ccall unsafe "setRectQ" c_setRectQ :: SetRect (Complex Float)
foreign import ccall unsafe "setRectI" c_setRectI :: SetRect I
foreign import ccall unsafe "setRectL" c_setRectL :: SetRect Z

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

sortG :: (Storable t, Storable a)
      => (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt) -> Vector t -> Vector a
sortG f v = unsafePerformIO $ do
    r <- createVector (dim v)
    (v #! r) f #|"sortG"
    return r

-- sortIdxD :: Vector Double -> Vector CInt
-- sortIdxD = sortG c_sort_indexD
sortIdxF :: Vector Float -> Vector CInt
sortIdxF = sortG c_sort_indexF
sortIdxI :: Vector CInt -> Vector CInt
sortIdxI = sortG c_sort_indexI
sortIdxL :: Vector Z -> Vector I
sortIdxL = sortG c_sort_indexL

-- sortValD :: Vector Double -> Vector Double
-- sortValD = sortG c_sort_valD
sortValF :: Vector Float -> Vector Float
sortValF = sortG c_sort_valF
sortValI :: Vector CInt -> Vector CInt
sortValI = sortG c_sort_valI
sortValL :: Vector Z -> Vector Z
sortValL = sortG c_sort_valL

-- foreign import ccall unsafe "sort_indexD" c_sort_indexD :: CV Double (CV CInt (IO CInt))
foreign import ccall unsafe "sort_indexF" c_sort_indexF :: CV Float  (CV CInt (IO CInt))
foreign import ccall unsafe "sort_indexI" c_sort_indexI :: CV CInt   (CV CInt (IO CInt))
foreign import ccall unsafe "sort_indexL" c_sort_indexL :: Z :> I :> Ok

-- foreign import ccall unsafe "sort_valuesD" c_sort_valD :: CV Double (CV Double (IO CInt))
foreign import ccall unsafe "sort_valuesF" c_sort_valF :: CV Float  (CV Float (IO CInt))
foreign import ccall unsafe "sort_valuesI" c_sort_valI :: CV CInt   (CV CInt (IO CInt))
foreign import ccall unsafe "sort_valuesL" c_sort_valL :: Z :> Z :> Ok

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

compareG :: (TransArray c, Storable t, Storable a)
         => Trans c (CInt -> Ptr t -> CInt -> Ptr a -> IO CInt)
         -> c -> Vector t -> Vector a
compareG f u v = unsafePerformIO $ do
    r <- createVector (dim v)
    (u # v #! r) f #|"compareG"
    return r

-- compareD :: Vector Double -> Vector Double -> Vector CInt
-- compareD = compareG c_compareD
compareF :: Vector Float -> Vector Float -> Vector CInt
compareF = compareG c_compareF
compareI :: Vector CInt -> Vector CInt -> Vector CInt
compareI = compareG c_compareI
compareL :: Vector Z -> Vector Z -> Vector CInt
compareL = compareG c_compareL

-- foreign import ccall unsafe "compareD" c_compareD :: CV Double (CV Double (CV CInt (IO CInt)))
foreign import ccall unsafe "compareF" c_compareF :: CV Float (CV Float  (CV CInt (IO CInt)))
foreign import ccall unsafe "compareI" c_compareI :: CV CInt (CV CInt   (CV CInt (IO CInt)))
foreign import ccall unsafe "compareL" c_compareL :: Z :> Z :> I :> Ok

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

selectG :: (TransArray c, TransArray c1, TransArray c2, Storable t, Storable a)
        => Trans c2 (Trans c1 (CInt -> Ptr t -> Trans c (CInt -> Ptr a -> IO CInt)))
        -> c2 -> c1 -> Vector t -> c -> Vector a
selectG f c u v w = unsafePerformIO $ do
    r <- createVector (dim v)
    (c # u # v # w #! r) f #|"selectG"
    return r

-- selectD :: Vector CInt -> Vector Double -> Vector Double -> Vector Double -> Vector Double
-- selectD = selectG c_selectD
selectF :: Vector CInt -> Vector Float -> Vector Float -> Vector Float -> Vector Float
selectF = selectG c_selectF
selectI :: Vector CInt -> Vector CInt -> Vector CInt -> Vector CInt -> Vector CInt
selectI = selectG c_selectI
selectL :: Vector CInt -> Vector Z -> Vector Z -> Vector Z -> Vector Z
selectL = selectG c_selectL
-- selectC :: Vector CInt
--         -> Vector (Complex Double)
--         -> Vector (Complex Double)
--         -> Vector (Complex Double)
--         -> Vector (Complex Double)
-- selectC = selectG c_selectC
selectQ :: Vector CInt
        -> Vector (Complex Float)
        -> Vector (Complex Float)
        -> Vector (Complex Float)
        -> Vector (Complex Float)
selectQ = selectG c_selectQ

type Sel x = CV CInt (CV x (CV x (CV x (CV x (IO CInt)))))

-- foreign import ccall unsafe "chooseD" c_selectD :: Sel Double
foreign import ccall unsafe "chooseF" c_selectF :: Sel Float
foreign import ccall unsafe "chooseI" c_selectI :: Sel CInt
-- foreign import ccall unsafe "chooseC" c_selectC :: Sel (Complex Double)
foreign import ccall unsafe "chooseQ" c_selectQ :: Sel (Complex Float)
foreign import ccall unsafe "chooseL" c_selectL :: Sel Z

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

remapG :: (TransArray c, TransArray c1, Storable t, Storable a)
       => (CInt -> CInt -> CInt -> CInt -> Ptr t
                -> Trans c1 (Trans c (CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)))
       -> Matrix t -> c1 -> c -> Matrix a
remapG f i j m = unsafePerformIO $ do
    r <- createMatrix RowMajor (rows i) (cols i)
    (i # j # m #! r) f #|"remapG"
    return r

-- remapD :: Matrix CInt -> Matrix CInt -> Matrix Double -> Matrix Double
-- remapD = remapG c_remapD
remapF :: Matrix CInt -> Matrix CInt -> Matrix Float -> Matrix Float
remapF = remapG c_remapF
remapI :: Matrix CInt -> Matrix CInt -> Matrix CInt -> Matrix CInt
remapI = remapG c_remapI
remapL :: Matrix CInt -> Matrix CInt -> Matrix Z -> Matrix Z
remapL = remapG c_remapL
-- remapC :: Matrix CInt
--        -> Matrix CInt
--        -> Matrix (Complex Double)
--        -> Matrix (Complex Double)
-- remapC = remapG c_remapC
remapQ :: Matrix CInt -> Matrix CInt -> Matrix (Complex Float) -> Matrix (Complex Float)
remapQ = remapG c_remapQ

type Rem x = OM CInt (OM CInt (OM x (OM x (IO CInt))))

-- foreign import ccall unsafe "remapD" c_remapD :: Rem Double
foreign import ccall unsafe "remapF" c_remapF :: Rem Float
foreign import ccall unsafe "remapI" c_remapI :: Rem CInt
-- foreign import ccall unsafe "remapC" c_remapC :: Rem (Complex Double)
foreign import ccall unsafe "remapQ" c_remapQ :: Rem (Complex Float)
foreign import ccall unsafe "remapL" c_remapL :: Rem Z

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

rowOpAux :: (TransArray c, Storable a) =>
            (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Trans c (IO CInt))
         -> Int -> a -> Int -> Int -> Int -> Int -> c -> IO ()
rowOpAux f c x i1 i2 j1 j2 m = do
    px <- newArray [x]
    (m # id) (f (fi c) px (fi i1) (fi i2) (fi j1) (fi j2)) #|"rowOp"
    free px

type RowOp x = CInt -> Ptr x -> CInt -> CInt -> CInt -> CInt -> x ::> Ok

-- foreign import ccall unsafe "rowop_double"  c_rowOpD :: RowOp R
foreign import ccall unsafe "rowop_float"   c_rowOpF :: RowOp Float
-- foreign import ccall unsafe "rowop_TCD"     c_rowOpC :: RowOp C
foreign import ccall unsafe "rowop_TCS"     c_rowOpQ :: RowOp (Complex Float)
foreign import ccall unsafe "rowop_int32_t" c_rowOpI :: RowOp I
foreign import ccall unsafe "rowop_int64_t" c_rowOpL :: RowOp Z
foreign import ccall unsafe "rowop_mod_int32_t" c_rowOpMI :: I -> RowOp I
foreign import ccall unsafe "rowop_mod_int64_t" c_rowOpML :: Z -> RowOp Z

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

gemmg :: (TransArray c1, TransArray c, TransArray c2, TransArray c3)
      => Trans c3 (Trans c2 (Trans c1 (Trans c (IO CInt))))
      -> c3 -> c2 -> c1 -> c -> IO ()
gemmg f v m1 m2 m3 = (v # m1 # m2 #! m3) f #|"gemmg"

type Tgemm x = x :> x ::> x ::> x ::> Ok

-- foreign import ccall unsafe "gemm_double"  c_gemmD :: Tgemm R
foreign import ccall unsafe "gemm_float"   c_gemmF :: Tgemm Float
-- foreign import ccall unsafe "gemm_TCD"     c_gemmC :: Tgemm C
foreign import ccall unsafe "gemm_TCS"     c_gemmQ :: Tgemm (Complex Float)
foreign import ccall unsafe "gemm_int32_t" c_gemmI :: Tgemm I
foreign import ccall unsafe "gemm_int64_t" c_gemmL :: Tgemm Z
foreign import ccall unsafe "gemm_mod_int32_t" c_gemmMI :: I -> Tgemm I
foreign import ccall unsafe "gemm_mod_int64_t" c_gemmML :: Z -> Tgemm Z

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

reorderAux :: (TransArray c, Storable t, Storable a1, Storable t1, Storable a) =>
              (CInt -> Ptr a -> CInt -> Ptr t1
                    -> Trans c (CInt -> Ptr t -> CInt -> Ptr a1 -> IO CInt))
           -> Vector t1 -> c -> Vector t -> Vector a1
reorderAux f s d v = unsafePerformIO $ do
    k <- createVector (dim s)
    r <- createVector (dim v)
    (k # s # d # v #! r) f #| "reorderV"
    return r

type Reorder x = CV CInt (CV CInt (CV CInt (CV x (CV x (IO CInt)))))

-- foreign import ccall unsafe "reorderD" c_reorderD :: Reorder Double
foreign import ccall unsafe "reorderF" c_reorderF :: Reorder Float
foreign import ccall unsafe "reorderI" c_reorderI :: Reorder CInt
-- foreign import ccall unsafe "reorderC" c_reorderC :: Reorder (Complex Double)
foreign import ccall unsafe "reorderQ" c_reorderQ :: Reorder (Complex Float)
foreign import ccall unsafe "reorderL" c_reorderL :: Reorder Z

-- | Transpose an array with dimensions @dims@ by making a copy using @strides@. For example, for an array with 3 indices,
--   @(reorderVector strides dims v) ! ((i * dims ! 1 + j) * dims ! 2 + k) == v ! (i * strides ! 0 + j * strides ! 1 + k * strides ! 2)@
--   This function is intended to be used internally by tensor libraries.
reorderVector :: Element a
                    => Vector CInt -- ^ @strides@: array strides
                    -> Vector CInt -- ^ @dims@: array dimensions of new array @v@
                    -> Vector a    -- ^ @v@: flattened input array
                    -> Vector a    -- ^ @v'@: flattened output array
reorderVector = reorderV

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

foreign import ccall unsafe "saveMatrix" c_saveMatrix
    :: CString -> CString -> Double ::> Ok

{- | save a matrix as a 2D ASCII table
-}
saveMatrix
    :: FilePath
    -> String        -- ^ \"printf\" format (e.g. \"%.2f\", \"%g\", etc.)
    -> Matrix Double
    -> IO ()
saveMatrix name format m = do
    cname   <- newCString name
    cformat <- newCString format
    (m # id) (c_saveMatrix cname cformat) #|"saveMatrix"
    free cname
    free cformat
    return ()

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

max_stars_repo_name: schnecki/hmatrix-float

max_issues_repo_name: schnecki/hmatrix-float

max_forks_repo_name: schnecki/hmatrix-float