2 Module : Gargantext.Core.Methods.Matrix.Accelerate.Utils
4 Copyright : (c) CNRS, 2017-Present
5 License : AGPL + CECILL v3
6 Maintainer : team@gargantext.org
7 Stability : experimental
10 This module aims at implementig distances of terms context by context is
11 the same referential of corpus.
13 Implementation use Accelerate library which enables GPU and CPU computation:
15 * Manuel M. T. Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell, and Vinod Grover.
16 [Accelerating Haskell Array Codes with Multicore GPUs][CKLM+11].
17 In _DAMP '11: Declarative Aspects of Multicore Programming_, ACM, 2011.
19 * Trevor L. McDonell, Manuel M. T. Chakravarty, Vinod Grover, and Ryan R. Newton.
20 [Type-safe Runtime Code Generation: Accelerate to LLVM][MCGN15].
21 In _Haskell '15: The 8th ACM SIGPLAN Symposium on Haskell_, ACM, 2015.
25 {-# LANGUAGE TypeFamilies #-}
26 {-# LANGUAGE TypeOperators #-}
27 {-# LANGUAGE ScopedTypeVariables #-}
28 {-# LANGUAGE ViewPatterns #-}
30 module Gargantext.Core.Methods.Matrix.Accelerate.Utils
33 import qualified Data.Foldable as P (foldl1)
34 import Debug.Trace (trace)
35 import Data.Array.Accelerate
36 import Data.Array.Accelerate.Interpreter (run)
37 import qualified Gargantext.Prelude as P
38 import Data.Array.Accelerate.LinearAlgebra hiding (Matrix, transpose, Vector)
40 -----------------------------------------------------------------------
42 -- Matrix Multiplication
48 => Acc (Array ((ix :. Int) :. Int) a)
49 -> Acc (Array ((ix :. Int) :. Int) a)
50 -> Acc (Array ((ix :. Int) :. Int) a)
51 (#*#) = multiplyMatrixMatrix
54 -- | Matrix cell by cell multiplication
60 => Acc (Array ((ix :. Int) :. Int) a)
61 -> Acc (Array ((ix :. Int) :. Int) a)
62 -> Acc (Array ((ix :. Int) :. Int) a)
70 , P.Fractional (Exp a)
72 => Acc (Array ((ix :. Int) :. Int) a)
73 -> Acc (Array ((ix :. Int) :. Int) a)
74 -> Acc (Array ((ix :. Int) :. Int) a)
81 , P.Fractional (Exp a)
83 => Acc (Array ((ix :. Int) :. Int) a)
84 -> Acc (Array ((ix :. Int) :. Int) a)
85 -> Acc (Array ((ix :. Int) :. Int) a)
92 , P.Fractional (Exp a)
94 => Acc (Array ((ix :. Int) :. Int) a)
95 -> Acc (Array ((ix :. Int) :. Int) a)
96 -> Acc (Array ((ix :. Int) :. Int) a)
99 -----------------------------------------------------------------------
100 matrixOne :: Num a => Dim -> Acc (Matrix a)
103 ones = fill (index2 n n) 1
107 matrixIdentity :: Num a => Dim -> Acc (Matrix a)
109 let zeros = fill (index2 n n) 0
110 ones = fill (index1 n) 1
113 permute const zeros (\(unindex1 -> i) -> index2 i i) ones
116 matrixEye :: Num a => Dim -> Acc (Matrix a)
118 let ones = fill (index2 n n) 1
119 zeros = fill (index1 n) 0
122 permute const ones (\(unindex1 -> i) -> index2 i i) zeros
125 diagNull :: Num a => Dim -> Acc (Matrix a) -> Acc (Matrix a)
126 diagNull n m = zipWith (*) m (matrixEye n)
128 -----------------------------------------------------------------------
129 _runExp :: Elt e => Exp e -> e
130 _runExp e = indexArray (run (unit e)) Z
132 -----------------------------------------------------------------------
136 -- Vector (Z :. 3) [0,1,2]
137 vector :: Elt c => Int -> [c] -> (Array (Z :. Int) c)
138 vector n l = fromList (Z :. n) l
142 -- >>> matrix 3 ([1..] :: [Double])
143 -- Matrix (Z :. 3 :. 3)
147 matrix :: Elt c => Int -> [c] -> Matrix c
148 matrix n l = fromList (Z :. n :. n) l
150 -- | Two ways to get the rank (as documentation)
152 -- >>> rank (matrix 3 ([1..] :: [Int]))
154 rank :: (Matrix a) -> Int
155 rank m = arrayRank $ arrayShape m
157 -----------------------------------------------------------------------
158 -- | Dimension of a square Matrix
159 -- How to force use with SquareMatrix ?
162 -- | Get Dimension of a square Matrix
164 -- >>> dim (matrix 3 ([1..] :: [Int]))
166 dim :: Matrix a -> Dim
169 Z :. _ :. n = arrayShape m
170 -- indexTail (arrayShape m)
172 -----------------------------------------------------------------------
174 -- | Sum of a Matrix by Column
176 -- >>> run $ matSumCol 3 (use $ matrix 3 [1..])
177 -- Matrix (Z :. 3 :. 3)
178 -- [ 12.0, 15.0, 18.0,
181 matSumCol :: (Elt a, P.Num (Exp a)) => Dim -> Acc (Matrix a) -> Acc (Matrix a)
182 matSumCol r mat = replicate (constant (Z :. (r :: Int) :. All)) $ sum $ transpose mat
184 matSumCol' :: (Elt a, P.Num (Exp a)) => Matrix a -> Matrix a
185 matSumCol' m = run $ matSumCol n m'
191 -- | Proba computes de probability matrix: all cells divided by thee sum of its column
192 -- if you need get the probability on the lines, just transpose it
194 -- >>> run $ matProba 3 (use $ matrix 3 [1..])
195 -- Matrix (Z :. 3 :. 3)
196 -- [ 8.333333333333333e-2, 0.13333333333333333, 0.16666666666666666,
197 -- 0.3333333333333333, 0.3333333333333333, 0.3333333333333333,
198 -- 0.5833333333333334, 0.5333333333333333, 0.5]
199 matProba :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
200 matProba d mat = zipWith (/) mat (matSumCol d mat)
202 -- | Diagonal of the matrix
204 -- >>> run $ diag (use $ matrix 3 ([1..] :: [Int]))
205 -- Vector (Z :. 3) [1,5,9]
209 diag m = backpermute (indexTail (shape m))
210 (lift1 (\(Z :. x) -> (Z :. x :. (x :: Exp Int))))
213 -- | Divide by the Diagonal of the matrix
215 -- >>> run $ divByDiag 3 (use $ matrix 3 ([1..] :: [Double]))
216 -- Matrix (Z :. 3 :. 3)
217 -- [ 1.0, 0.4, 0.3333333333333333,
218 -- 4.0, 1.0, 0.6666666666666666,
220 divByDiag :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
221 divByDiag d mat = zipWith (/) mat (replicate (constant (Z :. (d :: Int) :. All)) $ diag mat)
223 -----------------------------------------------------------------------
224 -- | Filters the matrix with the minimum of maximums
226 -- >>> run $ matMiniMax $ use $ matrix 3 [1..]
227 -- Matrix (Z :. 3 :. 3)
231 matMiniMax :: (Elt a, Ord a, P.Num a)
234 matMiniMax m = filterWith' miniMax' (constant 0) m
236 miniMax' = the $ minimum $ maximum m
239 -- | Filters the matrix with a constant
241 -- >>> run $ matFilter 5 $ use $ matrix 3 [1..]
242 -- Matrix (Z :. 3 :. 3)
246 filter' :: Double -> Acc (Matrix Double) -> Acc (Matrix Double)
247 filter' t m = filterWith t 0 m
249 filterWith :: Double -> Double -> Acc (Matrix Double) -> Acc (Matrix Double)
250 filterWith t v m = map (\x -> ifThenElse (x > (constant t)) x (constant v)) (transpose m)
252 filterWith' :: (Elt a, Ord a) => Exp a -> Exp a -> Acc (Matrix a) -> Acc (Matrix a)
253 filterWith' t v m = map (\x -> ifThenElse (x > t) x v) m
256 ------------------------------------------------------------------------
257 ------------------------------------------------------------------------
262 data Direction = MatCol (Exp Int) | MatRow (Exp Int) | Diag
264 nullOf :: Num a => Dim -> Direction -> Acc (Matrix a)
266 let ones = fill (index2 n n) 1
267 zeros = fill (index2 n n) 0
270 permute const ones ( lift1 ( \(Z :. (i :: Exp Int) :. (_j:: Exp Int))
272 MatCol m -> (Z :. i :. m)
273 MatRow m -> (Z :. m :. i)
274 Diag -> (Z :. i :. i)
279 nullOfWithDiag :: Num a => Dim -> Direction -> Acc (Matrix a)
280 nullOfWithDiag n dir = zipWith (*) (nullOf n dir) (nullOf n Diag)
283 divide :: (Elt a, Ord a, P.Fractional (Exp a), P.Num a)
284 => Acc (Matrix a) -> Acc (Matrix a) -> Acc (Matrix a)
285 divide = zipWith divide'
287 divide' a b = ifThenElse (b > (constant 0))
292 sumRowMin :: (Num a, Ord a) => Dim -> Acc (Matrix a) -> Acc (Matrix a)
293 sumRowMin n m = {-trace (P.show $ run m') $-} m'
295 m' = reshape (shape m) vs
297 $ P.map (\z -> sumRowMin1 n (constant z) m) [0..n-1]
299 sumRowMin1 :: (Num a, Ord a) => Dim -> Exp Int -> Acc (Matrix a) -> Acc (Vector a)
300 sumRowMin1 n x m = trace (P.show (run m,run $ transpose m)) $ m''
302 m'' = sum $ zipWith min (transpose m) m
303 _m' = zipWith (*) (zipWith (*) (nullOf n (MatCol x)) $ nullOfWithDiag n (MatRow x)) m
306 sumColMin :: (Num a, Ord a) => Dim -> Acc (Matrix a) -> Acc (Matrix a)
307 sumColMin n m = reshape (shape m) vs
310 $ P.map (\z -> sumColMin1 n (constant z) m) [0..n-1]
313 sumColMin1 :: (Num a) => Dim -> Exp Int -> Acc (Matrix a) -> Acc (Matrix a)
314 sumColMin1 n x m = zipWith (*) (nullOfWithDiag n (MatCol x)) m
318 {- | WIP fun with indexes
319 selfMatrix :: Num a => Dim -> Acc (Matrix a)
321 let zeros = fill (index2 n n) 0
322 ones = fill (index2 n n) 1
325 permute const ones ( lift1 ( \(Z :. (i :: Exp Int) :. (_j:: Exp Int))
326 -> -- ifThenElse (i /= j)
331 selfMatrix' :: (Elt a, P.Num (Exp a)) => Array DIM2 a -> Matrix a
332 selfMatrix' m' = run $ selfMatrix n
337 -------------------------------------------------
338 -------------------------------------------------
339 crossProduct :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
340 crossProduct n m = {-trace (P.show (run m',run m'')) $-} zipWith (*) m' m''
343 m'' = transpose $ cross n m
346 crossT :: Matrix Double -> Matrix Double
347 crossT = run . transpose . use
349 crossProduct' :: Matrix Double -> Matrix Double
350 crossProduct' m = run $ crossProduct n m'
355 runWith :: (Arrays c, Elt a1)
356 => (Dim -> Acc (Matrix a1) -> a2 -> Acc c)
360 runWith f m = run . f (dim m) (use m)
363 cross :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
364 cross n mat = diagNull n (matSumCol n $ diagNull n mat)
366 cross' :: Matrix Double -> Matrix Double
367 cross' mat = run $ cross n mat'
374 -- | Hypothesis to test maybe later (or not)
375 -- TODO ask accelerate for instances to ease such writtings:
376 p_ :: (Elt e, P.Fractional (Exp e)) => Acc (Array DIM2 e) -> Acc (Array DIM2 e)
377 p_ m = zipWith (/) m (n_ m)
379 n_ :: Elt e => Acc (SymetricMatrix e) -> Acc (Matrix e)
380 n_ m = backpermute (shape m)
381 (lift1 ( \(Z :. (i :: Exp Int) :. (j:: Exp Int))
382 -> (ifThenElse (i < j) (lift (Z :. j :. j)) (lift (Z :. i :. i)) :: Exp DIM2)
387 theMatrixDouble :: Int -> Matrix Double
388 theMatrixDouble n = run $ map fromIntegral (use $ theMatrixInt n)
390 theMatrixInt :: Int -> Matrix Int
391 theMatrixInt n = matrix n (dataMatrix n)
393 dataMatrix :: Int -> [Int]
394 dataMatrix x | (P.==) x 2 = [ 1, 1
398 | (P.==) x 3 = [ 7, 4, 0
402 | (P.==) x 4 = [ 4, 4, 0, 0
407 | P.otherwise = P.undefined
410 theResult :: Int -> Matrix Double
411 theResult n | (P.==) n 2 = let r = 1.6094379124341003 in [ 0, r, r, 0]
412 | P.otherwise = [ 1, 1 ]
417 => Int -> [e] -> Acc (Array ((Z :. Int) :. Int) e)
418 colMatrix n ns = replicate (constant (Z :. (n :: Int) :. All)) v
420 v = use $ vector (P.length ns) ns
422 -----------------------------------------------------------------------