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