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