src/Gargantext/Core/Methods/Matrix/Accelerate/Utils.hs

   1 {-|
   2 Module      : Gargantext.Core.Methods.Matrix.Accelerate.Utils
   3 Description :
   4 Copyright   : (c) CNRS, 2017-Present
   5 License     : AGPL + CECILL v3
   6 Maintainer  : team@gargantext.org
   7 Stability   : experimental
   8 Portability : POSIX
   9
  10 This module aims at implementig distances of terms context by context is
  11 the same referential of corpus.
  12
  13 Implementation use Accelerate library which enables GPU and CPU computation:
  14
  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.
  18
  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.
  22
  23 -}
  24
  25 {-# LANGUAGE TypeFamilies        #-}
  26 {-# LANGUAGE TypeOperators       #-}
  27 {-# LANGUAGE ScopedTypeVariables #-}
  28 {-# LANGUAGE ViewPatterns        #-}
  29
  30 module Gargantext.Core.Methods.Matrix.Accelerate.Utils
  31   where
  32
  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
  39 import Debug.Trace (trace)
  40
  41 -- | Matrix cell by cell multiplication
  42 (.*) :: ( Shape ix
  43         , Slice ix
  44         , Elt a
  45         , P.Num (Exp a)
  46         )
  47      => Acc (Array ((ix :. Int) :. Int) a)
  48      -> Acc (Array ((ix :. Int) :. Int) a)
  49      -> Acc (Array ((ix :. Int) :. Int) a)
  50 (.*) = zipWith (*)
  51
  52
  53 (./) :: ( Shape ix
  54         , Slice ix
  55         , Elt a
  56         , P.Num (Exp a)
  57         , P.Fractional (Exp a)
  58         )
  59      => Acc (Array ((ix :. Int) :. Int) a)
  60      -> Acc (Array ((ix :. Int) :. Int) a)
  61      -> Acc (Array ((ix :. Int) :. Int) a)
  62 (./) = zipWith (/)
  63
  64 -- | Term by term division where divisions by 0 produce 0 rather than NaN.
  65 termDivNan :: ( Shape ix
  66         , Slice ix
  67         , Elt a
  68         , Eq a
  69         , P.Num (Exp a)
  70         , P.Fractional (Exp a)
  71         )
  72      => Acc (Array ((ix :. Int) :. Int) a)
  73      -> Acc (Array ((ix :. Int) :. Int) a)
  74      -> Acc (Array ((ix :. Int) :. Int) a)
  75 termDivNan = trace "termDivNan" $ zipWith (\i j -> cond ((==) j 0) 0 ((/) i j))
  76
  77 (.-) :: ( Shape ix
  78         , Slice ix
  79         , Elt a
  80         , P.Num (Exp a)
  81         , P.Fractional (Exp a)
  82         )
  83      => Acc (Array ((ix :. Int) :. Int) a)
  84      -> Acc (Array ((ix :. Int) :. Int) a)
  85      -> Acc (Array ((ix :. Int) :. Int) a)
  86 (.-) = zipWith (-)
  87
  88 (.+) :: ( Shape ix
  89         , Slice ix
  90         , Elt a
  91         , P.Num (Exp a)
  92         , P.Fractional (Exp a)
  93         )
  94      => Acc (Array ((ix :. Int) :. Int) a)
  95      -> Acc (Array ((ix :. Int) :. Int) a)
  96      -> Acc (Array ((ix :. Int) :. Int) a)
  97 (.+) = zipWith (+)
  98
  99 -----------------------------------------------------------------------
 100 matrixOne :: Num a => Dim -> Acc (Matrix a)
 101 matrixOne n' = ones
 102   where
 103     ones  = fill (index2 n n) 1
 104     n     = constant n'
 105
 106
 107 matrixIdentity :: Num a => Dim -> Acc (Matrix a)
 108 matrixIdentity n' =
 109         let zeros = fill (index2 n n) 0
 110             ones  = fill (index1 n)   1
 111             n = constant n'
 112         in
 113         permute const zeros (\(unindex1 -> i) -> Just_ $ index2 i i) ones
 114
 115
 116 matrixEye :: Num a => Dim -> Acc (Matrix a)
 117 matrixEye n' =
 118         let ones   = fill (index2 n n) 1
 119             zeros  = fill (index1 n)   0
 120             n = constant n'
 121         in
 122         permute const ones (\(unindex1 -> i) -> Just_ $ index2 i i) zeros
 123
 124
 125 diagNull :: Num a => Dim -> Acc (Matrix a) -> Acc (Matrix a)
 126 diagNull n m = trace ("diagNull") $ zipWith (*) m (matrixEye n)
 127
 128
 129 -- Returns an N-dimensional array with the values of x for the indices where
 130 -- the condition is true, 0 everywhere else.
 131 condOrDefault
 132   :: forall sh a. (Shape sh, Elt a)
 133   => (Exp sh -> Exp Bool) -> Exp a -> Acc (Array sh a) -> Acc (Array sh a)
 134 condOrDefault theCond def x = permute const zeros filterInd x
 135   where
 136     zeros = fill (shape x) (def)
 137     filterInd ix = (cond (theCond ix)) (Just_ ix) Nothing_
 138
 139 -----------------------------------------------------------------------
 140 _runExp :: Elt e => Exp e -> e
 141 _runExp e = indexArray (run (unit e)) Z
 142
 143 -----------------------------------------------------------------------
 144 -- | Define a vector
 145 --
 146 -- >>> vector 3
 147 -- Vector (Z :. 3) [0,1,2]
 148 vector :: Elt c => Int -> [c] -> (Array (Z :. Int) c)
 149 vector n l = fromList (Z :. n) l
 150
 151 -- | Define a matrix
 152 --
 153 -- >>> matrix 3 ([1..] :: [Double])
 154 -- Matrix (Z :. 3 :. 3)
 155 --   [ 1.0, 2.0, 3.0,
 156 --     4.0, 5.0, 6.0,
 157 --     7.0, 8.0, 9.0]
 158 matrix :: Elt c => Int -> [c] -> Matrix c
 159 matrix n l = fromList (Z :. n :. n) l
 160
 161 -- | Two ways to get the rank (as documentation)
 162 --
 163 -- >>> rank (matrix 3 ([1..] :: [Int]))
 164 -- 2
 165 rank :: (Matrix a) -> Int
 166 rank m = arrayRank m
 167
 168 -----------------------------------------------------------------------
 169 -- | Dimension of a square Matrix
 170 -- How to force use with SquareMatrix ?
 171 type Dim = Int
 172
 173 -- | Get Dimension of a square Matrix
 174 --
 175 -- >>> dim (matrix 3 ([1..] :: [Int]))
 176 -- 3
 177 dim :: Matrix a -> Dim
 178 dim m = n
 179   where
 180     Z :. _ :. n = arrayShape m
 181     -- indexTail (arrayShape m)
 182
 183 -----------------------------------------------------------------------
 184
 185 -- | Sum of a Matrix by Column
 186 --
 187 -- >>> run $ matSumCol 3 (use $ matrix 3 [1..])
 188 -- Matrix (Z :. 3 :. 3)
 189 --   [ 12.0, 15.0, 18.0,
 190 --     12.0, 15.0, 18.0,
 191 --     12.0, 15.0, 18.0]
 192 matSumCol :: (Elt a, P.Num (Exp a)) => Dim -> Acc (Matrix a) -> Acc (Matrix a)
 193 matSumCol r mat = replicate (constant (Z :. (r :: Int) :. All)) $ sum $ transpose mat
 194
 195 matSumLin :: (Elt a, P.Num (Exp a)) => Dim -> Acc (Matrix a) -> Acc (Matrix a)
 196 matSumLin r mat = replicate (constant (Z :. (r :: Int) :. All)) $ sum mat
 197
 198 matSumCol' :: (Elt a, P.Num (Exp a)) => Matrix a -> Matrix a
 199 matSumCol' m = run $ matSumCol n m'
 200   where
 201     n  = dim m
 202     m' = use m
 203
 204
 205 -- | Proba computes de probability matrix: all cells divided by thee sum of its column
 206 -- if you need get the probability on the lines, just transpose it
 207 --
 208 -- >>> run $ matProba 3 (use $ matrix 3 [1..])
 209 -- Matrix (Z :. 3 :. 3)
 210 --   [ 8.333333333333333e-2, 0.13333333333333333, 0.16666666666666666,
 211 --       0.3333333333333333,  0.3333333333333333,  0.3333333333333333,
 212 --       0.5833333333333334,  0.5333333333333333,                 0.5]
 213 matProba :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
 214 matProba d mat = zipWith (/) mat (matSumCol d mat)
 215
 216
 217
 218 -- | Diagonal of the matrix
 219 --
 220 -- >>> run $ diag (use $ matrix 3 ([1..] :: [Int]))
 221 -- Vector (Z :. 3) [1,5,9]
 222 diag :: Elt e
 223      => Acc (Matrix e)
 224      -> Acc (Vector e)
 225 diag m = backpermute (indexTail (shape m))
 226                      (lift1 (\(Z :. x) -> (Z :. x :. (x :: Exp Int))))
 227                      m
 228
 229 -- | Divide by the Diagonal of the matrix
 230 --
 231 -- >>> run $ divByDiag 3 (use $ matrix 3 ([1..] :: [Double]))
 232 -- Matrix (Z :. 3 :. 3)
 233 --   [ 1.0, 0.4, 0.3333333333333333,
 234 --     4.0, 1.0, 0.6666666666666666,
 235 --     7.0, 1.6,                1.0]
 236 divByDiag :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
 237 divByDiag d mat = zipWith (/) mat (replicate (constant (Z :. (d :: Int) :. All)) $ diag mat)
 238
 239 -----------------------------------------------------------------------
 240 -- | Filters the matrix with the minimum of maximums
 241 --
 242 -- >>> run $ matMiniMax $ use $ matrix 3 [1..]
 243 -- Matrix (Z :. 3 :. 3)
 244 --   [ 0.0, 4.0, 7.0,
 245 --     0.0, 5.0, 8.0,
 246 --     0.0, 6.0, 9.0]
 247 matMiniMax :: (Elt a, Ord a, P.Num a)
 248            => Acc (Matrix a)
 249            -> Acc (Matrix a)
 250 matMiniMax m = filterWith' (>=) miniMax' (constant 0) m
 251   where
 252     miniMax' = the $ minimum $ maximum m
 253
 254 matMaxMini :: (Elt a, Ord a, P.Num a)
 255            => Acc (Matrix a)
 256            -> Acc (Matrix a)
 257 matMaxMini m = filterWith' (>) miniMax' (constant 0) m
 258   where
 259     miniMax' = the $ maximum $ minimum m
 260
 261
 262
 263
 264
 265 -- | Filters the matrix with a constant
 266 --
 267 -- >>> run $ matFilter 5 $ use $ matrix 3 [1..]
 268 -- Matrix (Z :. 3 :. 3)
 269 --   [ 0.0, 0.0, 7.0,
 270 --     0.0, 0.0, 8.0,
 271 --     0.0, 6.0, 9.0]
 272 filter' :: Double -> Acc (Matrix Double) -> Acc (Matrix Double)
 273 filter' t m = filterWith t 0 m
 274
 275 filterWith :: Double -> Double -> Acc (Matrix Double) -> Acc (Matrix Double)
 276 filterWith t v m = map (\x -> ifThenElse (x > (constant t)) x (constant v)) (transpose m)
 277
 278 filterWith' :: (Elt a, Ord a) => (Exp a -> Exp a -> Exp Bool) -> Exp a -> Exp a -> Acc (Matrix a) -> Acc (Matrix a)
 279 filterWith' f t v m = map (\x -> ifThenElse (f x t) x v) m
 280
 281
 282 ------------------------------------------------------------------------
 283 ------------------------------------------------------------------------
 284 -- | TODO use Lenses
 285 data Direction = MatCol (Exp Int) | MatRow (Exp Int) | Diag
 286
 287 nullOf :: Num a => Dim -> Direction -> Acc (Matrix a)
 288 nullOf n' dir =
 289         let ones   = fill (index2 n n) 1
 290             zeros  = fill (index2 n n) 0
 291             n = constant n'
 292         in
 293         permute const ones ( Just_ . lift1 ( \(Z :. (i :: Exp Int) :. (_j:: Exp Int))
 294                                                 -> case dir of
 295                                                      MatCol m -> (Z :. i :. m)
 296                                                      MatRow m -> (Z :. m :. i)
 297                                                      Diag     -> (Z :. i :. i)
 298                                    )
 299                            )
 300                            zeros
 301
 302 nullOfWithDiag :: Num a => Dim -> Direction -> Acc (Matrix a)
 303 nullOfWithDiag n dir = zipWith (*) (nullOf n dir) (nullOf n Diag)
 304
 305
 306 divide :: (Elt a, Ord a, P.Fractional (Exp a), P.Num a)
 307     => Acc (Matrix a) -> Acc (Matrix a) -> Acc (Matrix a)
 308 divide = zipWith divide'
 309   where
 310     divide' a b = ifThenElse (b > (constant 0))
 311                              (a / b)
 312                              (constant 0)
 313
 314 -- | Nominator
 315 sumRowMin :: (Num a, Ord a) => Dim -> Acc (Matrix a) -> Acc (Matrix a)
 316 sumRowMin n m = {-trace (P.show $ run m') $-} m'
 317   where
 318     m' = reshape (shape m) vs
 319     vs = P.foldl1 (++)
 320        $ P.map (\z -> sumRowMin1 n (constant z) m) [0..n-1]
 321
 322 sumRowMin1 :: (Num a, Ord a) => Dim -> Exp Int -> Acc (Matrix a) -> Acc (Vector a)
 323 sumRowMin1 n x m = {-trace (P.show (run m,run $ transpose m)) $-} m''
 324   where
 325     m'' = sum $ zipWith min (transpose m) m
 326     _m'  = zipWith (*) (zipWith (*) (nullOf n (MatCol x)) $ nullOfWithDiag n (MatRow x)) m
 327
 328 -- | Denominator
 329 sumColMin :: (Num a, Ord a) => Dim -> Acc (Matrix a) -> Acc (Matrix a)
 330 sumColMin n m = reshape (shape m) vs
 331   where
 332     vs = P.foldl1 (++)
 333        $ P.map (\z -> sumColMin1 n (constant z) m) [0..n-1]
 334
 335
 336 sumColMin1 :: (Num a) => Dim -> Exp Int -> Acc (Matrix a) -> Acc (Matrix a)
 337 sumColMin1 n x m = zipWith (*) (nullOfWithDiag n (MatCol x)) m
 338
 339
 340
 341 {- | WIP fun with indexes
 342 selfMatrix :: Num a => Dim -> Acc (Matrix a)
 343 selfMatrix n' =
 344         let zeros = fill (index2 n n) 0
 345             ones  = fill (index2 n n) 1
 346             n = constant n'
 347         in
 348         permute const ones ( lift1 ( \(Z :. (i :: Exp Int) :. (_j:: Exp Int))
 349                                                 -> -- ifThenElse (i /= j)
 350                                                      --         (Z :. i :. j)
 351                                                               (Z :. i :. i)
 352                                     )) zeros
 353
 354 selfMatrix' :: (Elt a, P.Num (Exp a)) => Array DIM2 a -> Matrix a
 355 selfMatrix' m' = run $ selfMatrix n
 356   where
 357     n = dim m'
 358     m = use m'
 359 -}
 360 -------------------------------------------------
 361 -------------------------------------------------
 362 crossProduct :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
 363 crossProduct n m = {-trace (P.show (run m',run m'')) $-} zipWith (*) m' m''
 364   where
 365     m'  = cross n m
 366     m'' = transpose $ cross n m
 367
 368
 369 crossT :: Matrix Double -> Matrix Double
 370 crossT  = run . transpose . use
 371
 372 crossProduct' :: Matrix Double -> Matrix Double
 373 crossProduct' m = run $ crossProduct n m'
 374   where
 375     n  = dim m
 376     m' = use m
 377
 378 runWith :: (Arrays c, Elt a1)
 379         => (Dim -> Acc (Matrix a1) -> a2 -> Acc c)
 380         -> Matrix a1
 381         -> a2
 382         -> c
 383 runWith f m = run . f (dim m) (use m)
 384
 385 -- | cross
 386 cross :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
 387 cross n mat = diagNull n (matSumCol n $ diagNull n mat)
 388
 389 cross' :: Matrix Double -> Matrix Double
 390 cross' mat = run $ cross n mat'
 391   where
 392     mat' = use mat
 393     n    = dim mat
 394
 395
 396 {-
 397 -- | Hypothesis to test maybe later (or not)
 398 -- TODO ask accelerate for instances to ease such writtings:
 399 p_ :: (Elt e, P.Fractional (Exp e)) => Acc (Array DIM2 e) -> Acc (Array DIM2 e)
 400 p_ m = zipWith (/) m (n_ m)
 401   where
 402     n_ :: Elt e => Acc (SymetricMatrix e) -> Acc (Matrix e)
 403     n_ m = backpermute (shape m)
 404                          (lift1 ( \(Z :. (i :: Exp Int) :. (j:: Exp Int))
 405                                    -> (ifThenElse (i < j) (lift (Z :. j :. j)) (lift (Z :. i :. i)) :: Exp DIM2)
 406                                 )
 407                          ) m
 408 -}
 409
 410 theMatrixDouble :: Int -> Matrix Double
 411 theMatrixDouble n = run $ map fromIntegral (use $ theMatrixInt n)
 412
 413 theMatrixInt :: Int -> Matrix Int
 414 theMatrixInt n = matrix n (dataMatrix n)
 415   where
 416     dataMatrix :: Int -> [Int]
 417     dataMatrix x | (P.==) x 2 = [ 1, 1
 418                                 , 1, 2
 419                                 ]
 420
 421                  | (P.==) x 3 =  [ 7, 4, 0
 422                                  , 4, 5, 3
 423                                  , 0, 3, 4
 424                                  ]
 425                  | (P.==) x 4 =  [ 4, 1, 2, 1
 426                                  , 1, 1, 0, 0
 427                                  , 2, 0, 5, 3
 428                                  , 1, 0, 3, 4
 429                                  ]
 430
 431
 432                  | P.otherwise = P.undefined
 433
 434 {-
 435 theResult :: Int -> Matrix Double
 436 theResult n | (P.==) n 2 = let r = 1.6094379124341003 in [ 0, r, r, 0]
 437           | P.otherwise = [ 1, 1 ]
 438 -}
 439
 440
 441 colMatrix :: Elt e
 442           => Int -> [e] -> Acc (Array ((Z :. Int) :. Int) e)
 443 colMatrix n ns = replicate (constant (Z :. (n :: Int) :. All)) v
 444   where
 445     v = use $ vector (P.length ns) ns
 446
 447 -----------------------------------------------------------------------
 448