Attempt to speed up logDistributional2

We are now around 6/7 times slower than the LLVM code.

bench/Main.hs

@@ -17,7 +17,7 @@ import Gargantext.Prelude.Crypto.Auth (createPasswordHash)
 import Paths_gargantext
 import qualified Data.Array.Accelerate as A
 import qualified Data.Array.Accelerate as Accelerate
-import qualified Data.Array.Accelerate.Interpreter as LLVM
+import qualified Data.Array.Accelerate.LLVM.Native as LLVM
 import qualified Data.Array.Accelerate.Interpreter as Naive
 import qualified Data.List.Split as Split
 import qualified Data.Massiv.Array as Massiv
@@ -25,6 +25,7 @@ import qualified Data.Massiv.Array.Numeric as Massiv
 import qualified Gargantext.Core.LinearAlgebra as LA
 import qualified Gargantext.Core.Methods.Matrix.Accelerate.Utils as Accelerate
 import qualified Gargantext.Core.Methods.Similarities.Accelerate.Distributional as Accelerate
+import qualified Numeric.LinearAlgebra.Data as HM
 import Test.Tasty.Bench
 import Data.Array.Accelerate ((:.))
 
@@ -95,20 +96,37 @@ main = do
       , bgroup "Phylo" [
           bench "toPhylo (small)" $ nf toPhylo issue290PhyloSmall
       ]
+      , bgroup "logDistributional2" [
+          bench "Accelerate (Naive)" $ nf (Accelerate.logDistributional2With @Double Naive.run) accInput
+      ,   bench "Accelerate (LLVM)"  $ nf (Accelerate.logDistributional2With @Double LLVM.run) accInput
+      ,   bench "Massiv"  $ nf (LA.logDistributional2 @_ @Double) massivInput
+      ]
+      , bgroup "distributional" [
+          bench "Accelerate (Naive)" $ nf (Accelerate.distributionalWith @Double Naive.run) accInput
+      ,   bench "Accelerate (LLVM)"  $ nf (Accelerate.distributionalWith @Double LLVM.run) accInput
+      ,   bench "Massiv (reference implementation)" $ nf (LA.distributionalReferenceImplementation @_ @Double) massivInput
+      ,   bench "Massiv " $ nf (LA.distributional @_ @Double) massivInput
+      ]
       , bgroup "diag" [
           bench "Accelerate (Naive)" $ nf (Naive.run . Accelerate.diag . Accelerate.use) accInput
       ,   bench "Accelerate (LLVM)"  $ nf (LLVM.run  . Accelerate.diag . Accelerate.use) accInput
       ,   bench "Massiv " $ nf (LA.diag @_) massivInput
       ]
-      , bgroup "identityMatrix" [
+      , bgroup "matrixIdentity" [
           bench "Accelerate (Naive)" $ nf (Naive.run . Accelerate.matrixIdentity @Double) 1000
       ,   bench "Accelerate (LLVM)"  $ nf (LLVM.run  . Accelerate.matrixIdentity @Double) 1000
-      ,   bench "Massiv " $ nf (Massiv.compute @Massiv.U . Massiv.identityMatrix @Double . Massiv.Sz1) 1000
+      ,   bench "Massiv" $ nf (LA.matrixIdentity @Double) 1000
+      ,   bench "HMatrix" $ nf (HM.ident @Double) 1000
+      ]
+      , bgroup "matrixEye" [
+          bench "Accelerate (Naive)" $ nf (Naive.run . Accelerate.matrixEye @Double) 1000
+      ,   bench "Accelerate (LLVM)"  $ nf (LLVM.run  . Accelerate.matrixEye @Double) 1000
+      ,   bench "Massiv " $ nf (LA.matrixEye @Double) 1000
       ]
       , bgroup "matMaxMini" [
           bench "Accelerate (Naive)" $ nf (\v -> Naive.run . Accelerate.matMaxMini @Double . Accelerate.use) accDoubleInput
       ,   bench "Accelerate (LLVM)"  $ nf (\v -> LLVM.run  . Accelerate.matMaxMini @Double . Accelerate.use) accDoubleInput
-      ,   bench "Massiv " $ nf (Massiv.compute @Massiv.U . LA.matMaxMini) massivDoubleInput
+      ,   bench "Massiv " $ nf LA.matMaxMini massivDoubleInput
       ]
       , bgroup "(.*)" [
           bench "Accelerate (Naive)" $ nf (\v -> Naive.run $ (Accelerate.use v) Accelerate..* (Accelerate.use v)) accDoubleInput
@@ -132,16 +150,5 @@ main = do
             nf (\m -> LLVM.run $ Accelerate.termDivNan (Accelerate.use m) (Accelerate.use m)) accDoubleInput
       ,   bench "Massiv " $ nf (\m -> LA.termDivNan @Massiv.U m m) massivDoubleInput
       ]
-      , bgroup "distributional" [
-          bench "Accelerate (Naive)" $ nf (Accelerate.distributionalWith @Double Naive.run) accInput
-      ,   bench "Accelerate (LLVM)"  $ nf Accelerate.distributional accInput
-      ,   bench "Massiv (reference implementation)" $ nf (LA.distributionalReferenceImplementation @_ @Double) massivInput
-      ,   bench "Massiv " $ nf (LA.distributional @_ @Double) massivInput
-      ]
-      , bgroup "logDistributional2" [
-          bench "Accelerate (Naive)" $ nf (Accelerate.logDistributional2With @Double Naive.run) accInput
-      ,   bench "Accelerate (LLVM)"  $ nf Accelerate.logDistributional2 accInput
-      ,   bench "Massiv"  $ nf (LA.logDistributional2 @_ @Double) massivInput
-      ]
       ]
     ]

gargantext.cabal

@@ -907,7 +907,8 @@ benchmark garg-bench
   build-depends: base
                , bytestring ^>= 0.11.5.3
                , accelerate
-               , accelerate-llvm
+               , accelerate-llvm-native
+               , hmatrix
                , massiv
                , deepseq
                , gargantext

src/Gargantext/Core/LinearAlgebra/Distributional.hs

@@ -2,6 +2,7 @@
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeApplications #-}
 {-# LANGUAGE TypeOperators #-}
+{-# OPTIONS_GHC -Wall-missed-specialisations #-}
 {-|
 Module      : Gargantext.Core.LinearAlgebra.Distributional
 Description : The "distributional" algorithm, fast and slow implementations
@@ -227,6 +228,7 @@ logDistributional2 m = A.computeP
                      $ logDistributional' n m
   where
     n = dim m
+{-# SPECIALIZE logDistributional2 :: Matrix A.U Int -> Matrix A.U Double #-}
 
 logDistributional' :: forall r e.
                    ( A.Manifest r e
@@ -244,8 +246,8 @@ logDistributional' :: forall r e.
                    -> Matrix r e
 logDistributional' n m' = result
  where
-    m :: Matrix A.D e
-    m = A.map fromIntegral m'
+    m :: Matrix A.U e
+    m = A.compute $ A.map fromIntegral m'
 
     -- Scalar. Sum of all elements of m.
     to :: e
@@ -257,7 +259,7 @@ logDistributional' n m' = result
 
     -- Size n vector. s = [s_i]_i
     s :: Vector A.U e
-    s = A.computeP $ sumRowsD (m `subD` d_m)
+    s = A.compute $ sumRowsD (m `subD` d_m)
 
     -- Matrix nxn. Vector s replicated as rows.
     s_1 :: Matrix D e

src/Gargantext/Core/LinearAlgebra/Operations.hs

@@ -179,19 +179,19 @@ dim m = n
   where
     (A.Sz2 _ n) = A.size m
 
-matMaxMini :: (A.Source r a, Ord a, Num a, A.Shape r A.Ix2) => Matrix r a -> Matrix D a
-matMaxMini m = A.map (\x -> if x > miniMax then x else 0) m
+matMaxMini :: (A.Unbox a, A.Source r a, Ord a, Num a, A.Shape r A.Ix2) => Matrix r a -> Matrix A.U a
+matMaxMini m = A.compute $ A.map (\x -> if x > miniMax then x else 0) m
   where
     -- Convert the matrix to a list of rows, take the minimum of each row,
     -- and then the maximum of those minima.
     miniMax = maximum (map minimum (A.toLists m))
 
-sumM_go :: (A.Manifest r a, Num a, A.Load r A.Ix2 a) => Int -> Matrix r a -> Matrix r a
-sumM_go n mi = A.makeArray (A.getComp mi) (A.Sz2 n n) $ \(i A.:. j) ->
+sumM_go :: (A.Unbox a, A.Manifest r a, Num a, A.Load r A.Ix2 a) => Int -> Matrix r a -> Matrix A.U a
+sumM_go n mi = A.makeArrayR A.U A.Seq (A.Sz2 n n) $ \(i A.:. j) ->
   Prelude.sum [ if k /= i && k /= j then mi A.! (i A.:. k) else 0 | k <- [0 .. n - 1] ]
 
-sumMin_go :: (A.Manifest r a, Num a, Ord a, A.Load r A.Ix2 a) => Int -> Matrix r a -> Matrix r a
-sumMin_go n mi = A.makeArray (A.getComp mi) (A.Sz2 n n) $ \(i A.:. j) ->
+sumMin_go :: (A.Unbox a, A.Manifest r a, Num a, Ord a, A.Load r A.Ix2 a) => Int -> Matrix r a -> Matrix A.U a
+sumMin_go n mi = A.makeArrayR A.U A.Seq (A.Sz2 n n) $ \(i A.:. j) ->
   Prelude.sum
     [ if k /= i && k /= j
         then min (mi A.! (i A.:. k)) (mi A.! (j A.:. k))
@@ -199,11 +199,15 @@ sumMin_go n mi = A.makeArray (A.getComp mi) (A.Sz2 n n) $ \(i A.:. j) ->
     | k <- [0 .. n - 1]
     ]
 
-matrixEye :: Num e => Int -> Matrix D e
-matrixEye n = A.makeArrayR A.D A.Seq (A.Sz2 n n) $ \(i A.:. j) -> if i == j then 0 else 1
+matrixEye :: (A.Unbox e, Num e) => Int -> Matrix A.U e
+matrixEye n = A.makeArrayR A.U A.Seq (A.Sz2 n n) $ \(i A.:. j) -> if i == j then 0 else 1
+{-# INLINE matrixEye #-}
+{-# SPECIALIZE matrixEye :: Int -> Matrix A.U Double #-}
 
-matrixIdentity :: Num e => Int -> Matrix D e
-matrixIdentity n = A.makeArrayR A.D A.Seq (A.Sz2 n n) $ \(i A.:. j) -> if i == j then 1 else 0
+matrixIdentity :: (A.Unbox e, Num e) => Int -> Matrix A.U e
+matrixIdentity n = A.makeArrayR A.U A.Seq (A.Sz2 n n) $ \(i A.:. j) -> if i == j then 1 else 0
+{-# INLINE matrixIdentity #-}
+{-# SPECIALIZE matrixIdentity :: Int -> Matrix A.U Double #-}
 
-diagNull :: (A.Source r e, Num e) => Int -> Matrix r e -> Matrix D e
-diagNull n m = A.zipWith (*) m (matrixEye n)
+diagNull :: (A.Unbox e, A.Source r e, Num e) => Int -> Matrix r e -> Matrix A.U e
+diagNull n m = A.compute $ A.zipWith (*) m (matrixEye n)

test/Test/Core/LinearAlgebra.hs

@@ -235,19 +235,19 @@ compareLogDistributional2 Proxy (SquareMatrix i1)
 
 compareMatMaxMini :: SquareMatrix Int -> Property
 compareMatMaxMini (SquareMatrix i1)
-  = let massiv     = LA.matMaxMini @Massiv.U (LA.accelerate2MassivMatrix i1)
+  = let massiv     = LA.matMaxMini (LA.accelerate2MassivMatrix i1)
         accelerate =  Naive.run (A.matMaxMini (use i1))
     in accelerate === LA.massiv2AccelerateMatrix massiv
 
 compareSumMin_go :: SquareMatrix Int -> Property
 compareSumMin_go (SquareMatrix i1)
-  = let massiv     = LA.sumMin_go @Massiv.U (A.dim i1) (LA.accelerate2MassivMatrix i1)
+  = let massiv     = LA.sumMin_go (A.dim i1) (LA.accelerate2MassivMatrix i1)
         accelerate =  Naive.run (Legacy.sumMin_go (A.dim i1) (use i1))
     in accelerate === LA.massiv2AccelerateMatrix massiv
 
 compareSumM_go :: SquareMatrix Int -> Property
 compareSumM_go (SquareMatrix i1)
-  = let massiv     = LA.sumM_go @Massiv.U (A.dim i1) (LA.accelerate2MassivMatrix i1)
+  = let massiv     = LA.sumM_go (A.dim i1) (LA.accelerate2MassivMatrix i1)
         accelerate =  Naive.run (Legacy.sumM_go (A.dim i1) (use i1))
     in accelerate === LA.massiv2AccelerateMatrix massiv