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