neural-network-hmatrix (empty) → 0.1.0.0
raw patch · 7 files changed
+739/−0 lines, 7 filesdep +basedep +deepseqdep +hmatrixsetup-changed
Dependencies added: base, deepseq, hmatrix, hmatrix-gsl, mtl, mwc-random, neural-network-base, parallel, vector
Files
- Data/NeuralNetwork/Backend/HMatrix.hs +91/−0
- Data/NeuralNetwork/Backend/HMatrix/Layers.hs +345/−0
- Data/NeuralNetwork/Backend/HMatrix/Utils.hs +177/−0
- LICENSE +29/−0
- Setup.hs +2/−0
- cbits/conv.c +70/−0
- neural-network-hmatrix.cabal +25/−0
+ Data/NeuralNetwork/Backend/HMatrix.hs view
@@ -0,0 +1,91 @@+{-# LANGUAGE MultiParamTypeClasses, FlexibleContexts, FlexibleInstances #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE TypeFamilies #-} +module Data.NeuralNetwork.Backend.HMatrix ( + module Data.NeuralNetwork.Backend.HMatrix.Layers, + ByHmatrix(..), + ErrCode(..) +) where + +import Data.NeuralNetwork +import Data.NeuralNetwork.Backend.HMatrix.Utils +import Data.NeuralNetwork.Backend.HMatrix.Layers +import Numeric.LinearAlgebra (Vector, Matrix) +import Control.Monad.Except +import Data.Functor.Identity + +data ErrCode = ErrMismatch +type Err = ExceptT ErrCode IO + +-- the backend type +data ByHmatrix = ByHmatrix + +-- with 1D input +instance (TranslateBody s, Component (RunLayer (SpecToTag s))) => + Backend ByHmatrix (SpecIn1D :++ s) where + type Env ByHmatrix = Err + type ConvertFromSpec (SpecIn1D :++ s) = RunLayer (SpecToTag s) + compile _ (a :++ l)= trans (size Nothing a) l + +-- with 2D input +instance (TranslateBody s, Component (RunLayer (SpecToTag s))) => + Backend ByHmatrix (SpecIn2D :++ s) where + type Env ByHmatrix = Err + type ConvertFromSpec (SpecIn2D :++ s) = RunLayer (SpecToTag s) + compile _ (a :++ l)= trans (size Nothing a) l + +instance RunInEnv Identity Err where + run = return . runIdentity + +-- It is necessary to propagate the size along the layers, +-- because fullconnect and convolution need to know +-- the previous size. +data Size = D1 Int | D2 Int Int Int + +class ComputeSize l where + size :: Maybe Size -> l -> Size +instance ComputeSize SpecIn1D where + size Nothing (In1D n) = D1 n +instance ComputeSize SpecIn2D where + size Nothing (In2D m n) = D2 1 m n +instance ComputeSize SpecReshape2DAs1D where + size (Just (D2 k m n)) _ = D1 (k*m*n) +instance ComputeSize SpecFullConnect where + size _ (FullConnect n) = D1 n +instance ComputeSize SpecConvolution where + size (Just (D2 _ m n)) (Convolution k f p) = D2 k (m+2*p-f+1) (n+2*p-f+1) +instance ComputeSize SpecMaxPooling where + size (Just (D2 k m n)) (MaxPooling s) = D2 k (m `div` s) (n `div` s) + +-- translate the body of specification +class TranslateBody s where + type SpecToTag s + trans :: Size -> s -> Err (RunLayer (SpecToTag s)) + +instance TranslateBody SpecFullConnect where + type SpecToTag SpecFullConnect = S F (T (SinglC :. Vector)) + trans (D1 s) (FullConnect n) = do u <- lift $ newFLayer s n + return $ Stack u (Activation (relu, relu')) + trans _ _ = throwError ErrMismatch + +instance TranslateBody SpecConvolution where + type SpecToTag SpecConvolution = S C (T (MultiC :. Matrix)) + trans (D2 k s t) (Convolution n f p) = do u <- lift $ newCLayer k n f p + return $ Stack u (Activation (relu, relu')) + trans _ _ = throwError ErrMismatch + +instance TranslateBody SpecReshape2DAs1D where + type SpecToTag SpecReshape2DAs1D = A + trans (D2 _ _ _) _ = return As1D + trans (D1 _) _ = throwError ErrMismatch + +instance TranslateBody SpecMaxPooling where + type SpecToTag SpecMaxPooling = M + trans (D2 _ _ _) (MaxPooling n) = return (MaxP n) + trans (D1 _) _ = throwError ErrMismatch + +instance (TranslateBody a, TranslateBody c, ComputeSize a) => TranslateBody (a :++ c) where + type SpecToTag (a :++ b) = S (SpecToTag a) (SpecToTag b) + trans s (a :++ c) = do u <- trans s a + v <- trans (size (Just s) a) c + return $ Stack u v
+ Data/NeuralNetwork/Backend/HMatrix/Layers.hs view
@@ -0,0 +1,345 @@+{-# LANGUAGE BangPatterns, TypeFamilies, TypeOperators, FlexibleInstances, FlexibleContexts, GADTs #-} +module Data.NeuralNetwork.Backend.HMatrix.Layers where + +import Numeric.LinearAlgebra hiding (R, C) +import Numeric.LinearAlgebra.Devel +import qualified Data.Vector as V +import qualified Data.Vector.Storable as SV +import qualified Data.Vector.Storable.Mutable as SVM +import System.Random.MWC +import System.Random.MWC.Distributions +import Control.Monad.ST +import Control.Monad (liftM2, forM_, when) +import GHC.Float +import Data.STRef +import Data.Functor.Identity +import Control.DeepSeq +import Data.NeuralNetwork +import Data.NeuralNetwork.Backend.HMatrix.Utils + +type R = Float + +-- We parameterise the activation layer T, where the parameter indicates how +-- elements are contained: +-- SinglC :. Vector, SinglC :. Matrix, MultiC :. Vector, MultiC :. Matrix +-- SinglC means the input has only one channel, while +-- MultiC means the input has more than one. +-- +-- type function composition +data (f :: * -> *) :. (g :: * -> *) :: * -> * +-- type function: Identity +data SinglC :: * -> * +data MultiC :: * -> * + +-- Tags for each form of layer +data F +data C +data A +data M +data T (c :: * -> *) +data S a b + +data RunLayer :: * -> * where + -- Densely connected layer + -- input: vector of size m + -- output: vector of size n + -- weights: matrix of size m x n + -- biases: vector of size n + Full :: !(Matrix R) -> !(Vector R) -> RunLayer F + -- convolutional layer + -- input: channels of 2D floats, of the same size (a x b), # of input channels: m + -- output: channels of 2D floats, of the same size (c x d), # of output channels: n + -- where c = a + 2*padding + 1 - s + -- d = b + 2*padding + 1 - t + -- feature: matrix of (s x t), # of features: m x n + -- padding: number of 0s padded at each side of channel + -- biases: bias for each output, # of biases: n + Conv :: !(V.Vector (V.Vector (Matrix R))) -> !(V.Vector R) -> Int -> RunLayer C + -- Reshape from channels of matrix to a single vector + -- input: m channels of 2D matrices + -- assuming that all matrices are of the same size a x b + -- output: 1D vector of the concatenation of all input channels + -- its size: m x a x b + As1D :: RunLayer A + -- max pooling layer + -- input: channels of 2D floats, of the same size (a x b), # of input channels: m + -- assuming that a and b are both multiple of stride + -- output: channels of 2D floats, of the same size (c x d), # of output channels: m + -- where c = a / stride + -- d = b / stride + MaxP :: Int -> RunLayer M + -- Activator + -- the input can be either a 1D vector, 2D matrix, or channels of either. + Activation :: (R->R, R->R) -> RunLayer (T c) + -- stacking two components a and b + -- the output of a should matches the input of b + Stack :: !(RunLayer a) -> !(RunLayer b) -> RunLayer (S a b) + +instance Component (RunLayer F) where + type Run (RunLayer F) = Identity + type Inp (RunLayer F) = Vector R + type Out (RunLayer F) = Vector R + -- trace is (input, weighted-sum) + newtype Trace (RunLayer F) = DTrace (Vector R, Vector R) + forwardT (Full w b) !inp = + let !bv = (inp <# w) `add` b + in return $ DTrace (inp,bv) + output (DTrace (_,!a)) = a + backward l (DTrace (!iv,!bv)) !odelta rate = + let Full w b = l + !d = scale (negate rate) odelta + !m = iv `outer` d + -- back-propagated error at input + !idelta = w #> odelta + -- update to weights + ---- for what reason, could this expression: w `add` (iv `outer` d) + ---- entails a huge space leak? especially, neither 'seq' nor + ---- 'deepseq' helps a bit. The only workaround is to expand the + ---- add function, and call SV.force on the result vector, which + ---- explcitly copy and drop reference to orignal computed result. + !w'= w `add` m + -- !w'= let (r,c) = size w + -- dat1 = flatten (tr' w) + -- dat2 = flatten (tr' m) + -- in matrixFromVector ColumnMajor r c $ SV.force $ dat1 `add` dat2 + !b'= b `add` d + -- !b'= SV.force $ b `add` d + in return $ (Full w' b', idelta) + +instance Component (RunLayer C) where + type Run (RunLayer C) = Identity + type Inp (RunLayer C) = V.Vector (Matrix R) + type Out (RunLayer C) = V.Vector (Matrix R) + -- trace is (input, convoluted output) + newtype Trace (RunLayer C) = CTrace (Inp (RunLayer C), V.Vector (Matrix R)) + forwardT (Conv fs bs p) !inp = + let !ov = parallel $ V.zipWith feature + (tr fs) -- feature matrix indexed majorly by each output + bs -- biases by each output + in return $ CTrace (inp,ov) + where + !osize = let (x,y) = size (V.head inp) + (u,v) = size (V.head $ V.head fs) + in (x+2*p-u+1, y+2*p-v+1) + -- transpose the features matrix + tr :: V.Vector (V.Vector a) -> V.Vector (V.Vector a) + tr uv = let n = V.length (V.head uv) + !vu = V.map (\i -> V.map (V.! i) uv) $ V.enumFromN 0 n + in vu + feature :: V.Vector (Matrix R) -> R -> Matrix R + feature f b = V.foldl1' add (V.zipWith (layerCorr2 p) f inp) `add` konst b osize + output (CTrace (_,a)) = a + backward l (CTrace (!iv,!av)) !odelta rate = + let Conv fs bs p = l + -- update to the feature matrix + m :: V.Vector (V.Vector (Matrix R)) + !m = parallel $ V.zipWith (\flts chn -> + -- chn: a single input channel + -- flts: all features used for chn + V.zipWith (\f d -> + let upd = scale (negate rate) (layerCorr2 p chn d) + in f `add` upd + ) flts odelta + ) fs iv + -- update to the biases + b :: V.Vector R + !b = V.zipWith (\b d -> b + (negate rate) * sumElements d) bs odelta + -- back-propagated error at input + idelta :: V.Vector (Matrix R) + !idelta = V.map (\f -> V.foldl1' add $ V.zipWith (layerConv2 p) f odelta) fs + in --trace ("CL:" ++ show odelta) + return $ (Conv m b p, idelta) + +instance Component (RunLayer A) where + type Run (RunLayer A) = Identity + type Inp (RunLayer A) = V.Vector (Matrix R) + type Out (RunLayer A) = Vector R + -- trace keeps information of (m, axb, b, output) + newtype Trace (RunLayer A) = ReshapeTrace (Int, Int, Int, Vector R) + forwardT _ !inp = + let !b = V.length inp + (!r,!c) = size (V.head inp) + !o = V.foldr' (\x y -> flatten x SV.++ y) SV.empty inp + in return $ ReshapeTrace (b, r*c, c, o) + output (ReshapeTrace (_,_,_,a)) = a + backward a (ReshapeTrace (b,n,c,_)) !odelta _ = + let !idelta = V.fromList $ map (reshape c) $ takesV (replicate b n) odelta + in return $ (a, idelta) + +instance Component (RunLayer M) where + type Run (RunLayer M) = Identity + type Inp (RunLayer M) = V.Vector (Matrix R) + type Out (RunLayer M) = V.Vector (Matrix R) + -- trace is (dimension of pools, index of max in each pool, pooled matrix) + -- for each channel. + newtype Trace (RunLayer M) = PTrace (V.Vector (IndexOf Matrix, Vector Int, Matrix R)) + -- forward is to divide the input matrix in stride x stride sub matrices, + -- and then find the max element in each sub matrices. + forwardT (MaxP stride) !inp = return $ PTrace $ parallel $ V.map mk inp + where + mk inp = let (!i,!v) = pool stride inp in (size v, i, v) + output (PTrace a) = V.map (\(_,_,!o) ->o) a + -- use the saved index-of-max in each pool to propagate the error. + backward l@(MaxP stride) (PTrace t) odelta _ = + let !idelta = V.zipWith gen t odelta in return $ (l, idelta) + where + gen (!si,!iv,_) od = unpool stride iv od + +instance (Component (RunLayer a), + Component (RunLayer b), + Run (RunLayer a) ~ Identity, + Run (RunLayer b) ~ Identity, + Out (RunLayer a) ~ Inp (RunLayer b) + ) => Component (RunLayer (S a b)) where + type Run (RunLayer (S a b)) = Identity + type Inp (RunLayer (S a b)) = Inp (RunLayer a) + type Out (RunLayer (S a b)) = Out (RunLayer b) + newtype Trace (RunLayer (S a b)) = TTrace (Trace (RunLayer b), Trace (RunLayer a)) + forwardT (Stack a b) !i = do + !tra <- forwardT a i + !trb <- forwardT b (output tra) + return $ TTrace (trb, tra) + output (TTrace !a) = output (fst a) + backward (Stack a b) (TTrace (!trb,!tra)) !odelta rate = do + (b', !odelta') <- backward b trb odelta rate + (a', !idelta ) <- backward a tra odelta' rate + return (Stack a' b', idelta) + +instance (Container c R) => Component (RunLayer (T (MultiC :. c))) where + type Run (RunLayer (T (MultiC :. c))) = Identity + type Inp (RunLayer (T (MultiC :. c))) = V.Vector (c R) + type Out (RunLayer (T (MultiC :. c))) = V.Vector (c R) + newtype Trace (RunLayer (T (MultiC :. c))) = TTraceM (V.Vector (Trace (RunLayer (T (SinglC :. c))))) + forwardT (Activation ac) !inp = + TTraceM <$> V.mapM (forwardT (Activation ac)) inp + output (TTraceM a) = V.map output a + backward a@(Activation ac) (TTraceM ts) !odelta r = do + idelta <- V.zipWithM (\t d -> snd <$> backward (Activation ac) t d r) ts odelta + return (a, idelta) + +instance (Container c R) => Component (RunLayer (T (SinglC :. c))) where + type Run (RunLayer (T (SinglC :. c))) = Identity + type Inp (RunLayer (T (SinglC :. c))) = c R + type Out (RunLayer (T (SinglC :. c))) = c R + newtype Trace (RunLayer (T (SinglC :. c))) = TTraceS (c R, c R) + forwardT (Activation (af,_)) !inp = return $ TTraceS (inp, cmap af inp) + output (TTraceS (_,!a)) = a + backward a@(Activation (_,ag)) (TTraceS (!iv,_)) !odelta _ = return $ (a, odelta `hadamard` cmap ag iv) + +newFLayer :: Int -- number of input values + -> Int -- number of neurons (output values) + -> IO (RunLayer F) -- new layer +newFLayer m n = + withSystemRandom . asGenIO $ \gen -> do + -- we build the weights in column major because in the back-propagation + -- algo, the computed update to weights is in column major. So it is + -- good for performance to keep the matrix always in column major. + w <- buildMatrix (normal 0 0.01 gen) ColumnMajor (m,n) + b <- return $ konst 1 n + return $ Full w b + +newCLayer :: Int -- number of input channels + -> Int -- number of output channels + -> Int -- size of each feature + -> Int -- size of padding + -> IO (RunLayer C) -- new layer +newCLayer inpsize outsize sfilter npadding = + withSystemRandom . asGenIO $ \gen -> do + fs <- V.replicateM inpsize $ V.replicateM outsize $ + buildMatrix (truncNormal 0 0.1 gen) RowMajor (sfilter, sfilter) + bs <- return $ V.replicate outsize 0.1 + return $ Conv fs bs npadding + where + truncNormal m s g = do + x <- standard g + if x >= 2.0 || x <= -2.0 + then truncNormal m s g + else return $! m + s * x + +buildMatrix g order (nr, nc) = do + vals <- SV.replicateM (nr*nc) (double2Float <$> g) + return $ matrixFromVector order nr nc vals + +layerCorr2 :: Int -> Matrix R -> Matrix R -> Matrix R +layerCorr2 p k m = c_corr2d_s k padded + where + padded = zeroPadded p m + (w,_) = size k + +layerConv2 :: Int -> Matrix R -> Matrix R -> Matrix R +layerConv2 p k m = c_conv2d_s k padded + where + padded = zeroPadded p m + (w,_) = size k + +-- max pool, picking out the maximum element +-- in each stride x stride sub-matrices. +-- assuming that the original matrix row and column size are +-- both multiple of stride +pool :: Int -> Matrix Float -> (Vector Int, Matrix Float) +pool 1 mat = let (r,c) = size mat in (SV.replicate (r*c) 0, mat) +-- pool 2 mat | orderOf mat == RowMajor = c_max_pool2_f mat +pool stride mat = runST $ do + ori <- unsafeThawMatrix mat + mxv <- newUndefinedMatrix RowMajor r' c' + mxi <- newUndefinedVector (r'*c') + forM_ [0..r'-1] $ \i -> do + forM_ [0..c'-1] $ \j -> do + (n,v) <- unsafeMaxIndEle ori (i*stride) (j*stride) stride stride + unsafeWriteVector mxi (i*c'+j) n + unsafeWriteMatrix mxv i j v + a <- unsafeFreezeVector mxi + b <- unsafeFreezeMatrix mxv + return (a,b) + where + (r,c) = size mat + r' = r `div` stride + c' = c `div` stride + unsafeMaxIndEle mm x y r c = do + mp <- newSTRef 0 + mv <- newSTRef (-10000.0) + forM_ [0..r-1] $ \ i -> do + forM_ [0..c-1] $ \ j -> do + v1 <- unsafeReadMatrix mm (x+i) (y+j) + v0 <- readSTRef mv + when (v1 > v0) $ do + writeSTRef mv v1 + writeSTRef mp (i*2+j) + p <- readSTRef mp + v <- readSTRef mv + return (p, v) + +-- the reverse of max pool. +-- assuming idx and mat are of the same size +unpool :: Int -> Vector Int -> Matrix Float -> Matrix Float +unpool stride idx mat = runSTMatrix $ do + mat' <- newMatrix' 0 r' c' + forM_ [0..r-1] $ \i -> do + forM_ [0..c-1] $ \j -> do + let pos = idx SV.! (i*c+j) + let (oi,oj) = pos `divMod` 2 + let val = mat `atIndex` (i,j) + unsafeWriteMatrix mat' (i*stride+oi) (j*stride+oj) val + return mat' + where + (r,c) = size mat + (r',c') = (r*stride, c*stride) + +-- a slightly faster way to pading the matrix +-- camparing to fromBlocks provided by hmatrix. +zeroPadded :: Int -> Matrix Float -> Matrix Float +zeroPadded p mat = runSTMatrix $ do + mat' <- newMatrix' 0 r' c' + setMatrix mat' p p mat + return mat' + where + (r,c) = size mat + (r',c') = (r+2*p,c+2*p) + +-- a slightly faster version of newMatrix, which based +-- directly on lower level Vector.Storable creation. +newMatrix' :: SVM.Storable t => t -> Int -> Int -> ST s (STMatrix s t) +newMatrix' v r c = do + vec <- SVM.replicate (r*c) v + vec <- SV.unsafeFreeze vec + unsafeThawMatrix $ reshape c vec
+ Data/NeuralNetwork/Backend/HMatrix/Utils.hs view
@@ -0,0 +1,177 @@+{-# LANGUAGE BangPatterns, FlexibleInstances, FlexibleContexts, ForeignFunctionInterface #-} +module Data.NeuralNetwork.Backend.HMatrix.Utils where +import Numeric.LinearAlgebra +import Numeric.LinearAlgebra.Devel +import Numeric.GSL.Fourier +import Data.Complex +import Control.Exception +import Control.DeepSeq +import Control.Parallel +import Control.Parallel.Strategies +import Control.Monad +import qualified Data.Vector as VecB +import qualified Data.Vector.Generic as VecGeneric +import qualified Data.Vector.Fusion.Bundle as VecFusion +import qualified Data.Vector.Fusion.Bundle.Monadic as VecFusionM +import qualified Data.Vector.Storable as VecS +import qualified Data.Vector.Storable.Mutable as VecM +import Foreign.Ptr ( Ptr ) +import Foreign.C.Types ( CInt(..) ) +import System.IO.Unsafe ( unsafePerformIO ) + +-- fft2d :: Matrix (Complex Double) -> Matrix (Complex Double) +-- fft2d m = let !x = fromRows $ map fft $ unsafeToRows m +-- !y = fromColumns $ map fft $ toColumns x +-- in y +-- ifft2d :: Matrix (Complex Double) -> Matrix (Complex Double) +-- ifft2d m = let !x = fromRows $ map ifft $ unsafeToRows m +-- !y = fromColumns $ map ifft $ toColumns x +-- in y + +-- fft2d :: Matrix (Complex Double) -> Matrix (Complex Double) +-- fft2d m = let rh:rr = map fft $ force $ toRows m +-- x = fromRows $ withStrategy (parList rdeepseq) rr `pseq` (rh : rr) +-- sh:ss = map fft $ force $ toColumns x +-- y = fromColumns $ withStrategy (parList rdeepseq) ss `pseq` (sh : ss) +-- in y +-- +-- ifft2d :: Matrix (Complex Double) -> Matrix (Complex Double) +-- ifft2d m = let rh:rr = map fft $ force $ toRows m +-- x = fromRows $ withStrategy (parList rdeepseq) rr `pseq` (rh : rr) +-- sh:ss = map fft $ force $ toColumns x +-- y = fromColumns $ withStrategy (parList rdeepseq) ss `pseq` (sh : ss) +-- in y + +-- conv2d_b :: (Numeric t, ConvFD t, Container Vector t, Container Matrix t) +-- => Matrix t -> Matrix t -> Matrix t +-- conv2d_b !k !m | w1 > w2 && h1 < h2 = error "convolution cannot be performed" +-- | w1 < w2 && h1 > h2 = error "convolution cannot be performed" +-- | w1 > w2 = conv2d_b m k +-- | otherwise = +-- -- convolution via FFT is actually cyclic conv. so we need to 0-pad the +-- -- matrix being convoluted by half the size of the kernel matrix. +-- -- the kernel matrix is also 0-padded to be the equal size. +-- -- finall, extra +-- let !hw = w1 `div` 2 +-- !hh = h1 `div` 2 +-- !z1 = konst 0 (w2-w1+hw,h2-h1+hh) +-- !z2 = konst 0 (hw, hh) +-- m1' = fft2d $ fromBlocks [[m1,0],[0,z1]] +-- m2' = fft2d $ fromBlocks [[m2,0],[0,z2]] +-- mr = ifft2d $ (force m1' `par` (force m2' `pseq` hadamard m1' m2')) +-- ms = subMatrix (w1-1,h1-1) (w2-w1+1,h2-h1+1) $ fst . fromComplex $ mr +-- in fromDouble $ ms +-- where +-- !m1 = complex $ toDouble k +-- !m2 = complex $ toDouble m +-- (w1,h1) = size m1 +-- (w2,h2) = size m2 +-- +-- corr2d_b k m | w > s && h < t = error "convolution cannot be performed" +-- | w < s && h > t = error "convolution cannot be performed" +-- | w > s = conv2d_b (rotate m) k +-- | otherwise = conv2d_b (rotate k) m +-- where +-- (w,h) = size k +-- (s,t) = size m +-- +-- corr2d_s :: (Numeric t, Container Vector t, Container Matrix t) +-- => Matrix t -> Matrix t -> Matrix t +-- corr2d_s k m | w > s && h < t = error "correlation cannot be performed" +-- | w < s && h > t = error "correlation cannot be performed" +-- | w > s = corr2d_s m k +-- | otherwise = final +-- where +-- (w,h) = size k +-- (s,t) = size m +-- (u,v) = (s-w+1, t-h+1) +-- subs = map (\s->subMatrix s (w,h) m) $ [(x,y) | x<-[0..u-1], y<-[0..v-1]] +-- -- use unsafe* methods to create the intermediate matrix fast. +-- t_rows = u*v +-- t_cols = w*h +-- !transformed = matrixFromVector RowMajor t_rows t_cols $ VecS.create $ do +-- mat <- VecM.new (t_rows*t_cols) +-- forM_ (zip [0..] subs) $ \(ri,rm) -> do +-- let bs = t_cols*ri +-- forM_ (zip [0..] $ unsafeToRows rm) $ \(ci, rv) -> do +-- let tv = VecM.unsafeSlice (bs+h*ci) h mat +-- {-# SCC "corr-data-copy" #-} VecS.unsafeCopy tv rv +-- return mat +-- !weights = flatten k +-- !final = reshape v $ transformed #> weights +-- +-- conv2d_s k m | w > s && h < t = error "convolution cannot be performed" +-- | w < s && h > t = error "convolution cannot be performed" +-- | w > s = corr2d_s (rotate m) k +-- | otherwise = corr2d_s (rotate k) m +-- where +-- (w,h) = size k +-- (s,t) = size m + +-- class ConvFD t where +-- fromDouble :: Matrix Double -> Matrix t +-- toDouble :: Matrix t -> Matrix Double +-- instance ConvFD Double where +-- fromDouble = id +-- toDouble = id +-- instance ConvFD Float where +-- fromDouble = single +-- toDouble = double +-- +-- rotate :: Element t => Matrix t -> Matrix t +-- rotate = fliprl . flipud + +foreign import ccall unsafe corr_sf_general :: + CInt -> + CInt -> CInt -> CInt -> Ptr Float -> + CInt -> CInt -> CInt -> Ptr Float -> + Ptr Float -> IO CInt + +data CConvType = CConv | CCorr + +c_corr2d_g :: CConvType -> Matrix Float -> Matrix Float -> Matrix Float +c_corr2d_g y k m | w > s = c_corr2d_s m k + | orderOf k == RowMajor && orderOf m == RowMajor = + let (r,c) = (s-w+1,t-h+1) + v = unsafePerformIO $ do + v <- VecM.unsafeNew (r * c) + VecM.unsafeWith v $ \rp -> + apply k id $ \kr kc ks kt kp -> + apply m id $ \mr mc ms mt mp -> + case y of + CCorr -> corr_sf_general 0 kr kc ks kp mr mc ms mp rp + CConv -> corr_sf_general 1 kr kc ks kp mr mc ms mp rp + VecS.unsafeFreeze v + in matrixFromVector RowMajor r c v + | otherwise = error "column major matrix not supported" + where + (w,h) = size k + (s,t) = size m + +c_corr2d_s = c_corr2d_g CCorr +c_conv2d_s = c_corr2d_g CConv + +-- parallel !vec = vec +parallel :: NFData a => VecB.Vector a -> VecB.Vector a +parallel vec = (VecB.tail vec `using` parvec) `pseq` vec + where + parvec = VecB.mapM (rparWith rdeepseq) + +-- foreign import ccall unsafe pool2_f :: CInt -> CInt -> CInt -> Ptr Float -> +-- Ptr Float -> Ptr CInt -> IO () +-- c_max_pool2_f :: Matrix Float -> (Vector Int, Matrix Float) +-- c_max_pool2_f mat +-- | orderOf mat == RowMajor = unsafePerformIO $ do +-- ind <- VecM.unsafeNew (r' * c') +-- mx <- VecM.unsafeNew (r' * c') +-- VecM.unsafeWith ind $ \pind -> +-- VecM.unsafeWith mx $ \pmax -> +-- apply mat id $ \mr mc ms mt mp -> +-- pool2_f mr mc ms mp pmax pind +-- ind <- VecS.unsafeFreeze ind +-- mx <- VecS.unsafeFreeze mx +-- return (VecS.map fromIntegral ind, matrixFromVector RowMajor r' c' mx) +-- where +-- (r,c) = size mat +-- r' = r `div` 2 +-- c' = c `div` 2
+ LICENSE view
@@ -0,0 +1,29 @@+BSD 3-Clause License + +Copyright (c) 2016, Jiasen Wu +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +* Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple +main = defaultMain
+ cbits/conv.c view
@@ -0,0 +1,70 @@+#include <stdlib.h> +#include <string.h> +#include "cblas.h" + +#define AT(p,s,x,y) ((p)+(s)*(x)+(y)) + + +// 2M byte working storage per thread. +// sufficient to handle correlation/convolution wthin +// the size (5x5, 128x128) or (7x7, 100x100) +#define WORKINGSTORAGE 524288 +static _Thread_local float workingSto[WORKINGSTORAGE]; + +/* +input: two matrices, assuming in row major, may not be in continuous. i.e. + one row has only column number of elements, but its length is of + the stride. + mat1 is the kernel + mat2 is the source + assuming that row1<=row2 and col1 <= col2 +output: continuous row-major matrix of size u x v + where u = row2-row1+1 + v = col2-co11+1 +*/ + +int corr_sf_general( + int reversemat1, int row1, int col1, int stride1, float *mat1, + int row2, int col2, int stride2, float *mat2, float *mat3) +{ + int u = row2-row1+1; + int v = col2-col1+1; + if ((u*v+1)*row1*col1 > WORKINGSTORAGE) return -1; + float *ws = workingSto; + for(int i=0;i<u;i++) { + for(int j=0;j<v;j++) { + for(int k=0;k<row1;k++) { + float *src=AT(mat2,stride2,i+k,j); + memcpy(ws, src, col1*sizeof(float)); + ws += col1; + } + } + } + float *vw = mat1; + if (stride1 > col1 || reversemat1) { + // we have expecting an extra row1*col1 elements at hand, pointed + // by ws at this point of time, and used for a continues + // copy of the kernel matrix. + vw = ws; + float *p1=vw, *p2=mat1; + for(int i=0;i<row1;i++) { + memcpy(p1,p2,col1*sizeof(float)); + p1 += col1; + p2 += stride1; + } + } + // reverse the kernel when do convolution. + if(reversemat1) { + int sz = row1*col1; + for(int i=0;i<sz/2;i++) { + float t = vw[i]; + vw[i] = vw[sz-1-i]; + vw[sz-1-i] = t; + } + } + // execute the matrix vector multiplication + ws = workingSto; + cblas_sgemv(CblasRowMajor, CblasNoTrans, u*v, row1*col1, + 1.0, ws, row1*col1, vw, 1, 0, mat3, 1); + return 0; +}
+ neural-network-hmatrix.cabal view
@@ -0,0 +1,25 @@+name: neural-network-hmatrix +version: 0.1.0.0 +license-file: LICENSE +license: BSD3 +author: Jiasen Wu +maintainer: jiasenwu@hotmail.com +Category: AI +Synopsis: Yet Another High Performance and Extendable Neural Network in Haskell +Description: Provides execution backend of neural network on top of hmatrix. +build-type: Simple +cabal-version: >=1.10 + +library + exposed-modules: Data.NeuralNetwork.Backend.HMatrix + other-modules: Data.NeuralNetwork.Backend.HMatrix.Utils + Data.NeuralNetwork.Backend.HMatrix.Layers + build-depends: base >= 4.7 && < 5, hmatrix, hmatrix-gsl, mwc-random, mtl, + vector, deepseq, parallel, neural-network-base + default-language: Haskell2010 + C-sources: cbits/conv.c + cc-options: -O3 -std=c11 -march=native + if os(windows) + extra-libraries: openblas + if os(linux) + extra-libraries: blas