git.iko.soy/grenade
on master help sign in
Add stupidly fast im2col in c
3830e55Huw Campbell
bench/bench.hs+4
compare with master
unifiedsplit
@@ -11,6 +11,10 @@ main = defaultMain [
11 , bench "im2col 28x28" $ whnf (im2colUnsafe 5 5 1 1) ((28><28) [1..])
11 , bench "im2col 28x28" $ whnf (im2colUnsafe 5 5 1 1) ((28><28) [1..])
12 , bench "im2col 100x100" $ whnf (im2colUnsafe 10 10 1 1) ((100><100) [1..])
12 , bench "im2col 100x100" $ whnf (im2colUnsafe 10 10 1 1) ((100><100) [1..])
13 ]
13 ]
14+ , bgroup "im2col_c" [ bench "im2col_c 3x4" $ whnf (im2col_c 2 2 1 1) ((3><4) [1..])
15+ , bench "im2col_c 28x28" $ whnf (im2col_c 5 5 1 1) ((28><28) [1..])
16+ , bench "im2col_c 100x100" $ whnf (im2col_c 10 10 1 1) ((100><100) [1..])
17+ ]
14 , bgroup "col2im" [ bench "col2im 3x4" $ whnf (col2imUnsafe 2 2 1 1 3 4) ((6><4) [1..])
18 , bgroup "col2im" [ bench "col2im 3x4" $ whnf (col2imUnsafe 2 2 1 1 3 4) ((6><4) [1..])
15 , bench "col2im 28x28" $ whnf (col2imUnsafe 5 5 1 1 28 28) ((576><25) [1..])
19 , bench "col2im 28x28" $ whnf (col2imUnsafe 5 5 1 1 28 28) ((576><25) [1..])
16 , bench "col2im 100x100" $ whnf (col2imUnsafe 10 10 1 1 100 100) ((8281><100) [1..])
20 , bench "col2im 100x100" $ whnf (col2imUnsafe 10 10 1 1 100 100) ((8281><100) [1..])
cbits/im2col.cnew file+30
@@ -0,0 +1,30 @@
1+#include "im2col.h"
2+
3+inline int is_a_ge_zero_and_a_lt_b(int a, int b) {
4+ return a >= 0 && a < b;
5+}
6+
7+void im2col_cpu(const double* data_im, int dataOffset, const int channels,
8+ const int height, const int width, const int kernel_h, const int kernel_w,
9+ const int stride_h, const int stride_w,
10+ double* data_col) {
11+
12+ data_im += dataOffset;
13+ const int output_h = (height - kernel_h) / stride_h + 1;
14+ const int output_w = (width - kernel_w) / stride_w + 1;
15+
16+ const int channel_size = height * width;
17+ for (int channel = channels; channel--; data_im += channel_size) {
18+ for (int fitting_height = 0; fitting_height <= (height - kernel_h); fitting_height += stride_h) {
19+ for (int fitting_width = 0; fitting_width <= (width - kernel_w); fitting_width += stride_w) {
20+ for (int kernel_row = 0; kernel_row < kernel_h; kernel_row++) {
21+ for (int kernel_col = 0; kernel_col < kernel_w; kernel_col++) {
22+ int input_row = fitting_height + kernel_row;
23+ int input_col = fitting_width + kernel_col;
24+ *(data_col++) = data_im[input_row * width + input_col];
25+ }
26+ }
27+ }
28+ }
29+ }
30+}
cbits/im2col.hnew file+7
@@ -0,0 +1,7 @@
1+#include <stdio.h>
2+#include <stdint.h>
3+
4+void im2col_cpu(const double* data_im, int dataOffset, const int channels,
5+ const int height, const int width, const int kernel_h, const int kernel_w,
6+ const int stride_h, const int stride_w,
7+ double* data_col);
grenade.cabal+6
@@ -10,6 +10,10 @@ cabal-version: >= 1.8
10build-type: Simple
10build-type: Simple
11description: grenade.
11description: grenade.
12
12
13+extra-source-files:
14+ cbits/im2col.h
15+ cbits/im2col.c
16+
13library
17library
14 build-depends:
18 build-depends:
15 base >= 4.8 && < 5
19 base >= 4.8 && < 5
@@ -53,6 +57,8 @@ library
53 Grenade.Layers.Internal.Convolution
57 Grenade.Layers.Internal.Convolution
54 Grenade.Layers.Internal.Pooling
58 Grenade.Layers.Internal.Pooling
55
59
60+ includes: cbits/im2col.h
61+ c-sources: cbits/im2col.c
56
62
57executable feedforward
63executable feedforward
58 ghc-options: -Wall -threaded -O2
64 ghc-options: -Wall -threaded -O2
src/Grenade/Layers/Convolution.hs+2−2
@@ -155,7 +155,7 @@ instance ( KnownNat kernelRows
155 sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
155 sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
156 ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
156 ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
157 oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
157 oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
158− c = im2colUnsafe kx ky sx sy ex
158+ c = im2col_c kx ky sx sy ex
159 mt = c LA.<> ek
159 mt = c LA.<> ek
160 r = col2vidUnsafe 1 1 1 1 ox oy mt
160 r = col2vidUnsafe 1 1 1 1 ox oy mt
161 rs = fmap (fromJust . create) r
161 rs = fmap (fromJust . create) r
@@ -172,7 +172,7 @@ instance ( KnownNat kernelRows
172 ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
172 ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
173 oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
173 oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)
174
174
175− c = im2colUnsafe kx ky sx sy ex
175+ c = im2col_c kx ky sx sy ex
176
176
177 eo = vecToList $ fmap extract dEdy
177 eo = vecToList $ fmap extract dEdy
178 ek = extract kernel
178 ek = extract kernel
src/Grenade/Layers/Internal/Convolution.hs+67−29
@@ -1,8 +1,11 @@
1+{-# LANGUAGE BangPatterns #-}
2+{-# LANGUAGE ForeignFunctionInterface #-}
1module Grenade.Layers.Internal.Convolution (
3module Grenade.Layers.Internal.Convolution (
2 col2vidUnsafe
4 col2vidUnsafe
3 , col2imUnsafe
5 , col2imUnsafe
4 , vid2colUnsafe
6 , vid2colUnsafe
5 , im2colUnsafe
7 , im2colUnsafe
8+ , im2col_c
6 , fittingStarts
9 , fittingStarts
7 , unsafeModifyMatrix
10 , unsafeModifyMatrix
8 ) where
11 ) where
@@ -13,11 +16,15 @@ import Data.STRef ( newSTRef, modifySTRef', writeSTRef, readSTRef )
13import Data.Foldable ( forM_ )
16import Data.Foldable ( forM_ )
14import Data.Traversable ( forM )
17import Data.Traversable ( forM )
15
18
19+import Foreign ( mallocForeignPtrArray0, withForeignPtr )
20+import Foreign.Ptr ( Ptr )
16import Foreign.Storable( Storable )
21import Foreign.Storable( Storable )
17
22
18import Numeric.LinearAlgebra hiding ( uniformSample, konst )
23import Numeric.LinearAlgebra hiding ( uniformSample, konst )
19import qualified Numeric.LinearAlgebra.Devel as U
24import qualified Numeric.LinearAlgebra.Devel as U
20
25
26+import System.IO.Unsafe ( unsafePerformIO )
27+
21-- This module provides provides im2col function and friends, ala caffe.
28-- This module provides provides im2col function and friends, ala caffe.
22--
29--
23-- /* From Caffe */
30-- /* From Caffe */
@@ -105,13 +112,13 @@ col2vidUnsafe kernelRows kernelColumns strideRows strideColumns destinationRows
105 offsetR <- newSTRef 0
112 offsetR <- newSTRef 0
106 offsetC <- newSTRef 0
113 offsetC <- newSTRef 0
107 forM_ [0 .. columnMatrixRows - 1] $ \ir -> do
114 forM_ [0 .. columnMatrixRows - 1] $ \ir -> do
108− inputColumn <- newSTRef 0
115+ inputColumn <- newSTRef offsetM
109− offsetR' <- readSTRef offsetR
116+ offsetR' <- readSTRef offsetR
110− offsetC' <- readSTRef offsetC
117+ offsetC' <- readSTRef offsetC
111 forM_ [offsetR' .. offsetR' + kernelRows -1] $ \kr ->
118 forM_ [offsetR' .. offsetR' + kernelRows -1] $ \kr ->
112 forM_ [offsetC' .. offsetC' + kernelColumns -1] $ \kc -> do
119 forM_ [offsetC' .. offsetC' + kernelColumns -1] $ \kc -> do
113 ic <- readSTRef inputColumn
120 ic <- readSTRef inputColumn
114− unsafeModifyMatrix dataIm kr kc (+ U.atM' columnMatrix ir (ic + offsetM))
121+ unsafeModifyMatrix dataIm kr kc (+ U.atM' columnMatrix ir ic)
115 modifySTRef' inputColumn (+1)
122 modifySTRef' inputColumn (+1)
116
123
117 if offsetC' + kernelColumns < destinationCols
124 if offsetC' + kernelColumns < destinationCols
@@ -120,6 +127,62 @@ col2vidUnsafe kernelRows kernelColumns strideRows strideColumns destinationRows
120
127
121 U.unsafeFreezeMatrix dataIm
128 U.unsafeFreezeMatrix dataIm
122
129
130+
131+im2col_c :: Int -> Int -> Int -> Int -> Matrix Double -> Matrix Double
132+im2col_c kernelRows kernelColumns strideRows strideColumns dataIm =
133+ let height = rows dataIm
134+ width = cols dataIm
135+ vec = flatten dataIm
136+ rowOut = (height - kernelRows) `div` strideRows + 1
137+ colOut = (width - kernelColumns) `div` strideColumns + 1
138+ kernelSize = kernelRows * kernelColumns
139+ numberOfPatches = rowOut * colOut
140+ in unsafePerformIO $ do
141+ outPtr <- mallocForeignPtrArray0 (numberOfPatches * kernelSize)
142+ let (inPtr, inOffset, _) = U.unsafeToForeignPtr vec
143+
144+ withForeignPtr inPtr $ \inPtr' ->
145+ withForeignPtr outPtr $ \outPtr' ->
146+ im2col_cpu inPtr' inOffset 1 height width kernelRows kernelColumns strideRows strideColumns outPtr'
147+
148+ let matVec = U.unsafeFromForeignPtr outPtr 0 (numberOfPatches * kernelSize)
149+ return $ U.matrixFromVector U.RowMajor numberOfPatches kernelSize matVec
150+
151+foreign import ccall safe
152+ im2col_cpu
153+ :: Ptr Double -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Ptr Double -> IO ()
154+
155+
156+unsafeModifyMatrix :: (Storable t) => U.STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()
157+unsafeModifyMatrix x r c f = U.unsafeReadMatrix x r c >>= U.unsafeWriteMatrix x r c . f
158+{-# INLINE unsafeModifyMatrix #-}
159+
160+
161+-- | Returns the starting sub matrix locations which fit inside the larger matrix for the
162+-- convolution. Takes into account the stride and kernel size.
163+fittingStarts :: Int -> Int -> Int -> Int -> Int -> Int -> [(Int,Int)]
164+fittingStarts nrows kernelrows steprows ncols kernelcols stepcolsh =
165+ let rs = fittingStart nrows kernelrows steprows
166+ cs = fittingStart ncols kernelcols stepcolsh
167+ in concatMap ( \r -> fmap (\c -> (r , c)) cs ) rs
168+
169+-- | Returns the starting sub vector which fit inside the larger vector for the
170+-- convolution. Takes into account the stride and kernel size.
171+fittingStart :: Int -> Int -> Int -> [Int]
172+fittingStart width kernel steps =
173+ let go left | left + kernel < width
174+ = left : go (left + steps)
175+ | left + kernel == width
176+ = [left]
177+ | otherwise
178+ = []
179+ in go 0
180+
181+
182+
183+
184+-- | Old functions (useful for sanity checking and benchmarking)
185+
123vid2colUnsafe :: Int -> Int -> Int -> Int -> Int -> Int -> [Matrix Double] -> Matrix Double
186vid2colUnsafe :: Int -> Int -> Int -> Int -> Int -> Int -> [Matrix Double] -> Matrix Double
124vid2colUnsafe kernelRows kernelColumns striderows stridecols vidrows vidcols dataVid = U.runSTMatrix $ do
187vid2colUnsafe kernelRows kernelColumns striderows stridecols vidrows vidcols dataVid = U.runSTMatrix $ do
125 let starts = fittingStarts vidrows kernelRows striderows vidcols kernelColumns stridecols
188 let starts = fittingStarts vidrows kernelRows striderows vidcols kernelColumns stridecols
@@ -168,28 +231,3 @@ im2colUnsafe kernelRows kernelColumns striderows stridecols dataIm = U.runSTMatr
168 modifySTRef' inputRowRef (+1)
231 modifySTRef' inputRowRef (+1)
169
232
170 return dataCol
233 return dataCol
171−
172−unsafeModifyMatrix :: (Storable t) => U.STMatrix s t -> Int -> Int -> (t -> t) -> ST s ()
173−unsafeModifyMatrix x r c f = U.unsafeReadMatrix x r c >>= U.unsafeWriteMatrix x r c . f
174−{-# INLINE unsafeModifyMatrix #-}
175−
176−
177−-- | Returns the starting sub matrix locations which fit inside the larger matrix for the
178−-- convolution. Takes into account the stride and kernel size.
179−fittingStarts :: Int -> Int -> Int -> Int -> Int -> Int -> [(Int,Int)]
180−fittingStarts nrows kernelrows steprows ncols kernelcols stepcolsh =
181− let rs = fittingStart nrows kernelrows steprows
182− cs = fittingStart ncols kernelcols stepcolsh
183− in concatMap ( \r -> fmap (\c -> (r , c)) cs ) rs
184−
185−-- | Returns the starting sub vector which fit inside the larger vector for the
186−-- convolution. Takes into account the stride and kernel size.
187−fittingStart :: Int -> Int -> Int -> [Int]
188−fittingStart width kernel steps =
189− let go left | left + kernel < width
190− = left : go (left + steps)
191− | left + kernel == width
192− = [left]
193− | otherwise
194− = []
195− in go 0
test/Test/Grenade/Layers/Convolution.hs+43
@@ -30,6 +30,21 @@ prop_im2col_no_stride = once $
30 out = im2colUnsafe 2 2 1 1 input
30 out = im2colUnsafe 2 2 1 1 input
31 in expected === out
31 in expected === out
32
32
33+prop_im2col_c = once $
34+ let input = (3><4)
35+ [ 1.0, 2.0, 3.0, 4.0
36+ , 5.0, 6.0, 7.0, 8.0
37+ , 9.0, 10.0, 11.0, 12.0 ]
38+ expected = (6><4)
39+ [ 1.0, 2.0, 5.0, 6.0
40+ , 2.0, 3.0, 6.0, 7.0
41+ , 3.0, 4.0, 7.0, 8.0
42+ , 5.0, 6.0, 9.0, 10.0
43+ , 6.0, 7.0, 10.0, 11.0
44+ , 7.0, 8.0, 11.0, 12.0 ]
45+ out = im2col_c 2 2 1 1 input
46+ in expected === out
47+
33prop_im2col_stride = once $
48prop_im2col_stride = once $
34 let input = (3><4)
49 let input = (3><4)
35 [ 1.0, 2.0, 3.0, 4.0
50 [ 1.0, 2.0, 3.0, 4.0
@@ -43,6 +58,19 @@ prop_im2col_stride = once $
43 out = im2colUnsafe 2 2 1 2 input
58 out = im2colUnsafe 2 2 1 2 input
44 in expected === out
59 in expected === out
45
60
61+prop_im2col_c_stride = once $
62+ let input = (3><4)
63+ [ 1.0, 2.0, 3.0, 4.0
64+ , 5.0, 6.0, 7.0, 8.0
65+ , 9.0, 10.0, 11.0, 12.0 ]
66+ expected = (4><4)
67+ [ 1.0, 2.0, 5.0, 6.0
68+ , 3.0, 4.0, 7.0, 8.0
69+ , 5.0, 6.0, 9.0, 10.0
70+ , 7.0, 8.0, 11.0, 12.0 ]
71+ out = im2col_c 2 2 1 2 input
72+ in expected === out
73+
46prop_im2col_other = once $
74prop_im2col_other = once $
47 let input = (3><4)
75 let input = (3><4)
48 [ 1.0, 2.0, 3.0, 4.0
76 [ 1.0, 2.0, 3.0, 4.0
@@ -54,6 +82,21 @@ prop_im2col_other = once $
54 out = im2colUnsafe 3 2 1 2 input
82 out = im2colUnsafe 3 2 1 2 input
55 in expected === out
83 in expected === out
56
84
85+prop_im2col_c_other = once $
86+ let input = (3><4)
87+ [ 1.0, 2.0, 3.0, 4.0
88+ , 5.0, 6.0, 7.0, 8.0
89+ , 9.0, 10.0, 11.0, 12.0 ]
90+ expected = (2><6)
91+ [ 1.0, 2.0, 5.0, 6.0 , 9.0, 10.0
92+ , 3.0, 4.0, 7.0, 8.0 , 11.0 ,12.0 ]
93+ out = im2col_c 3 2 1 2 input
94+ in expected === out
95+
96+prop_im2col_bigger = once $
97+ let input = (7><7) [ 1.0 .. ]
98+ in im2colUnsafe 5 5 2 2 input === im2col_c 5 5 2 2 input
99+
57-- If there's no overlap (stride is the same size as the kernel)
100-- If there's no overlap (stride is the same size as the kernel)
58-- then col2im . im2col should be symmetric.
101-- then col2im . im2col should be symmetric.
59prop_im2col_sym_on_same_stride = once $
102prop_im2col_sym_on_same_stride = once $