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