Merge branch 'adinapoli/issue-291' into 'dev'

Replace performance-critical parts of the algorithm with `massiv` (and ditch `accelerate-llvm`)

See merge request !382

bench/Main.hs

 {-# LANGUAGE TypeApplications  #-}
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE NumericUnderscores #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeOperators #-}
 
 module Main where
 
@@ -9,9 +14,19 @@ import Gargantext.Core.Viz.Phylo.API.Tools (readPhylo)
 import Gargantext.Core.Viz.Phylo.PhyloMaker (toPhylo)
 import Gargantext.Core.Viz.Phylo.PhyloTools
 import Gargantext.Prelude.Crypto.Auth (createPasswordHash)
-
-import Test.Tasty.Bench
 import Paths_gargantext
+import qualified Data.Array.Accelerate as A
+import qualified Data.Array.Accelerate as Accelerate
+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
+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 ((:.))
 
 phyloConfig :: PhyloConfig
 phyloConfig = PhyloConfig {
@@ -37,10 +52,38 @@ phyloConfig = PhyloConfig {
   , exportFilter = [ByBranchSize {_branch_size = 3.0}]
   }
 
+matrixValues :: [Int]
+matrixValues = [ 1 .. 10_000 ]
+
+matrixDim :: Int
+matrixDim = 100
+
+testMatrix :: A.Matrix Int
+testMatrix = A.fromList (A.Z A.:. matrixDim A.:. matrixDim) $ matrixValues
+{-# INLINE testMatrix #-}
+
+testVector :: A.Array (A.Z :. Int :. Int :. Int) Int
+testVector = A.fromList (A.Z A.:. 20 A.:. 20 A.:. 20) $ matrixValues
+{-# INLINE testVector #-}
+
+testMassivMatrix :: Massiv.Matrix Massiv.U Int
+testMassivMatrix = Massiv.fromLists' Massiv.Par $ Split.chunksOf matrixDim $ matrixValues
+{-# INLINE testMassivMatrix #-}
+
+testMassivVector :: Massiv.Array Massiv.U Massiv.Ix3 Int
+testMassivVector = LA.accelerate2Massiv3DMatrix testVector
+{-# INLINE testMassivVector #-}
+
 main :: IO ()
 main = do
   _issue290Phylo     <- force . setConfig phyloConfig <$> (readPhylo =<< getDataFileName "bench-data/phylo/issue-290.json")
   issue290PhyloSmall <- force . setConfig phyloConfig <$> (readPhylo =<< getDataFileName "bench-data/phylo/issue-290-small.json")
+  let !accInput       = force testMatrix
+  let !accVector      = force testVector
+  let !massivVector   = force testMassivVector
+  let !(accDoubleInput :: Accelerate.Matrix Double) = force $ Naive.run $ Accelerate.map Accelerate.fromIntegral (Accelerate.use testMatrix)
+  let !massivInput    = force testMassivMatrix
+  let !(massivDoubleInput :: Massiv.Matrix Massiv.U Double) = force $ Massiv.computeP $ Massiv.map fromIntegral testMassivMatrix
   defaultMain
     [ bgroup "Benchmarks"
       [ bgroup "User creation" [
@@ -51,5 +94,59 @@ 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 "matrixIdentity" [
+          bench "Accelerate (Naive)" $ nf (Naive.run . Accelerate.matrixIdentity @Double) 1000
+      ,   bench "Accelerate (LLVM)"  $ nf (LLVM.run  . Accelerate.matrixIdentity @Double) 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 (Naive.run . Accelerate.matMaxMini @Double . Accelerate.use) accDoubleInput
+      ,   bench "Accelerate (LLVM)"  $ nf (LLVM.run  . Accelerate.matMaxMini @Double . Accelerate.use) accDoubleInput
+      ,   bench "Massiv " $ nf LA.matMaxMini massivDoubleInput
+      ]
+      , bgroup "(.*)" [
+          bench "Accelerate (Naive)" $ nf (\v -> Naive.run $ (Accelerate.use v) Accelerate..* (Accelerate.use v)) accDoubleInput
+      ,   bench "Accelerate (LLVM)"  $ nf (\v -> LLVM.run $ (Accelerate.use v) Accelerate..* (Accelerate.use v)) accDoubleInput
+      ,   bench "Massiv " $ nf (\v -> (v LA..* v) :: Massiv.Matrix Massiv.U Double) massivDoubleInput
+      ]
+      , bgroup "sumRows" [
+          bench "Accelerate (Naive)" $ nf (Naive.run . Accelerate.sum . Accelerate.use) accVector
+      ,   bench "Accelerate (LLVM)"  $ nf (LLVM.run . Accelerate.sum . Accelerate.use) accVector
+      ,   bench "Massiv "            $ nf LA.sumRows massivVector
+      ]
+      , bgroup "sumMin_go" [
+          bench "Accelerate (Naive)" $ nf (Naive.run . Accelerate.sumMin_go 100 . Accelerate.use) accDoubleInput
+      ,   bench "Accelerate (LLVM)"  $ nf (LLVM.run . Accelerate.sumMin_go 100 . Accelerate.use) accDoubleInput
+      ,   bench "Massiv "            $ nf (Massiv.compute @Massiv.U . LA.sumMin_go 100) massivDoubleInput
+      ]
+      , bgroup "termDivNan" [
+          bench "Accelerate (Naive)" $
+            nf (\m -> Naive.run $ Accelerate.termDivNan (Accelerate.use m) (Accelerate.use m)) accDoubleInput
+      ,   bench "Accelerate (LLVM)" $
+            nf (\m -> LLVM.run $ Accelerate.termDivNan (Accelerate.use m) (Accelerate.use m)) accDoubleInput
+      ,   bench "Massiv " $ nf (\m -> LA.termDivNan @Massiv.U m m) massivDoubleInput
+      ]
       ]
     ]

bin/update-project-dependencies

@@ -18,8 +18,8 @@ fi
 # with the `sha256sum` result calculated on the `cabal.project` and
 # `cabal.project.freeze`. This ensures the files stay deterministic so that CI
 # cache can kick in.
-expected_cabal_project_hash="ac293a4c66092996bc85fbf14ef34b7cce3ed5b0612ceb9e1a5f395059631e0b"
-expected_cabal_project_freeze_hash="32310c4d4e7b4679dcb90dcfcd0d6d1b175dbf885a77ffddca16d422998a521c"
+expected_cabal_project_hash="ae8a153144dbf82c8676cd494d5afced39122fca8eac4facf20d3aec1c7c5dcc"
+expected_cabal_project_freeze_hash="50f40d4fc0190b8e8645e2de3c3d40ddb3664339e7995f52f883d925685a9a49"
 
 
 cabal --store-dir=$STORE_DIR v2-build --dry-run

cabal.project

@@ -5,6 +5,7 @@ index-state: 2024-09-12T03:02:26Z
 
 with-compiler: ghc-9.4.8
 optimization: 2
+benchmarks: False
 
 packages:
     ./
@@ -14,13 +15,6 @@ source-repository-package
     location: https://github.com/AccelerateHS/accelerate.git
     tag: 334d05519436bb7f20f9926ec76418f5b8afa359
 
-source-repository-package
-    type: git
-    location: https://github.com/AccelerateHS/accelerate-llvm.git
-    tag: 2b5d69448557e89002c0179ea1aaf59bb757a6e3
-    subdir: accelerate-llvm-native/
-            accelerate-llvm/
-
 -- Patch for "Allow NOT to backtrack"
 source-repository-package
     type: git
@@ -32,13 +26,6 @@ source-repository-package
     location: https://gitlab.iscpif.fr/gargantext/opaleye-textsearch.git
     tag: 04b5c9044fef44393b66bffa258ca0b0f59c1087
 
-source-repository-package
-    type: git
-    location: https://github.com/adinapoli/llvm-hs.git
-    tag: 7533a9ccd3bfe77141745f6b61039a26aaf5c83b
-    subdir: llvm-hs
-            llvm-hs-pure
-
 source-repository-package
     type: git
     location: https://github.com/alpmestan/accelerate-arithmetic.git

cabal.project.freeze

@@ -26,8 +26,6 @@ constraints: any.Boolean ==0.2.4,
              any.accelerate ==1.3.0.0,
              accelerate +bounds-checks -debug -internal-checks -nofib -unsafe-checks,
              any.accelerate-arithmetic ==1.0.0.1,
-             any.accelerate-llvm ==1.3.0.0,
-             any.accelerate-llvm-native ==1.3.0.0,
              any.accelerate-utility ==1.0.0.1,
              any.adjunctions ==4.4.2,
              any.aeson ==2.1.2.1,
@@ -308,9 +306,6 @@ constraints: any.Boolean ==0.2.4,
              any.linear ==1.23,
              linear -herbie +template-haskell,
              any.list-t ==1.0.5.7,
-             any.llvm-hs ==12.0.0,
-             llvm-hs -debug -llvm-with-rtti +shared-llvm,
-             any.llvm-hs-pure ==12.0.0,
              any.lockfree-queue ==0.2.4,
              any.logict ==0.8.1.0,
              any.loop ==0.3.0,

gargantext.cabal

@@ -100,6 +100,10 @@ flag no-phylo-debug-logs
   default: False
   manual: True
 
+flag enable-benchmarks
+  default: False
+  manual: True
+
 library
   import:
     defaults
@@ -188,9 +192,14 @@ library
       Gargantext.Core.Config.Types
       Gargantext.Core.Config.Utils
       Gargantext.Core.Config.Worker
+      Gargantext.Core.LinearAlgebra
+      Gargantext.Core.LinearAlgebra.Distributional
+      Gargantext.Core.LinearAlgebra.Operations
       Gargantext.Core.Mail
       Gargantext.Core.Mail.Types
+      Gargantext.Core.Methods.Matrix.Accelerate.Utils
       Gargantext.Core.Methods.Similarities
+      Gargantext.Core.Methods.Similarities.Accelerate.Distributional
       Gargantext.Core.Methods.Similarities.Conditional
       Gargantext.Core.NLP
       Gargantext.Core.NodeStory
@@ -295,6 +304,9 @@ library
       Gargantext.Database.Schema.User
       Gargantext.Defaults
       Gargantext.MicroServices.ReverseProxy
+      Gargantext.Orphans
+      Gargantext.Orphans.Accelerate
+      Gargantext.Orphans.OpenAPI
       Gargantext.System.Logging
       Gargantext.Utils.Dict
       Gargantext.Utils.Jobs.Error
@@ -303,6 +315,7 @@ library
       Gargantext.Utils.SpacyNLP.Types
       Gargantext.Utils.Tuple
       Gargantext.Utils.Zip
+      Paths_gargantext
   other-modules:
       Gargantext.API.Admin.Auth
       Gargantext.API.Admin.FrontEnd
@@ -356,9 +369,7 @@ library
       Gargantext.Core.Flow.Ngrams
       Gargantext.Core.Flow.Types
       Gargantext.Core.Methods.Graph.MaxClique
-      Gargantext.Core.Methods.Matrix.Accelerate.Utils
       Gargantext.Core.Methods.Similarities.Accelerate.Conditional
-      Gargantext.Core.Methods.Similarities.Accelerate.Distributional
       Gargantext.Core.Methods.Similarities.Accelerate.SpeGen
       Gargantext.Core.Statistics
       Gargantext.Core.Text.Corpus
@@ -470,12 +481,9 @@ library
       Gargantext.Database.Schema.NodeNode
       Gargantext.Database.Schema.Prelude
       Gargantext.Database.Types
-      Gargantext.Orphans
-      Gargantext.Orphans.OpenAPI
       Gargantext.Utils.Aeson
       Gargantext.Utils.Servant
       Gargantext.Utils.UTCTime
-      Paths_gargantext
   ghc-options: -Wall -Wincomplete-uni-patterns -Wincomplete-record-updates -Wmissing-signatures -Wunused-binds -Wunused-imports -Wunused-packages -Werror -freduction-depth=300 -fprint-potential-instances
   hs-source-dirs:
       src
@@ -489,7 +497,6 @@ library
     , MonadRandom ^>= 0.6
     , QuickCheck ^>= 2.14.2
     , accelerate ^>= 1.3.0.0
-    , accelerate-llvm-native ^>= 1.3.0.0
     , aeson ^>= 2.1.2.1
     , ansi-terminal
     , array ^>= 0.5.4.0
@@ -554,6 +561,7 @@ library
     , json-stream ^>= 0.4.2.4
     , lens >= 5.2.2 && < 5.3
     , lens-aeson < 1.3
+    , massiv < 1.1
     , matrix ^>= 0.3.6.1
     , mime-mail >= 0.5.1
     , monad-control ^>= 1.0.3.1
@@ -586,6 +594,7 @@ library
     , replace-attoparsec ^>= 1.5.0.0
     , resource-pool >= 0.4.0.0 && < 0.5
     , safe-exceptions >= 0.1.7.4 && < 0.2
+    , scientific < 0.4
     , serialise ^>= 0.2.4.0
     , servant >= 0.20.1 && < 0.21
     , servant-auth ^>= 0.4.0.0
@@ -609,7 +618,7 @@ library
     , singletons ^>= 3.0.2
     , singletons-th >= 3.1 && < 3.2
     , smtp-mail >= 0.3.0.0
-    , split >= 0.2.0
+    , split >= 0.2.3.4
     , stemmer == 0.5.2
     , stm >= 2.5.1.0 && < 2.6
     , stm-containers >= 1.2.0.3 && < 1.3
@@ -694,14 +703,18 @@ executable gargantext
     , servant-routes < 0.2
     , servant-websockets >= 2.0.0 && < 2.1
     , shelly
-    , split ^>= 0.2.3.4
+    , split >= 0.2.3.4
     , text ^>= 2.0.2
     , toml-parser >= 2.0.1.0 && < 3
     , tree-diff
     , vector >= 0.12.3.0 && <= 0.13.1.0
   ghc-options: -Wall -Wincomplete-uni-patterns -Wincomplete-record-updates -Wmissing-signatures -Wunused-binds -Wunused-imports -Werror -freduction-depth=300 -fprint-potential-instances
 
-common testDependencies
+-- Dependencies needed by both test stanzas. Dependencies
+-- which are not shared must be placed in the `build-depends`
+-- of each separate stanza to not trigger any warning if `-Wunused-packages`
+-- is enabled.
+common commonTestDependencies
   build-depends:
       base >=4.7 && <5
     , QuickCheck ^>= 2.14.2
@@ -763,17 +776,20 @@ common testDependencies
 test-suite garg-test-tasty
   import:
       defaults
-    , testDependencies
+    , commonTestDependencies
   type: exitcode-stdio-1.0
   main-is: drivers/tasty/Main.hs
   build-depends:
       aeson-pretty ^>= 0.8.9
+    , accelerate >= 1.3.0.0
     , boolexpr ^>= 0.3
     , conduit ^>= 1.3.4.2
     , crawlerArxiv
     , cryptohash
     , directory ^>= 1.3.7.1
     , graphviz ^>= 2999.20.1.0
+    , massiv < 1.1
+    , scientific < 0.4
     , split
     , tasty >= 1.4.3 && < 1.6
     , tasty-golden
@@ -795,6 +811,7 @@ test-suite garg-test-tasty
       Test.API.Setup
       Test.API.Prelude
       Test.API.UpdateList
+      Test.Core.LinearAlgebra
       Test.Core.Notifications
       Test.Core.Orchestrator
       Test.Core.Similarity
@@ -845,7 +862,7 @@ test-suite garg-test-tasty
 test-suite garg-test-hspec
   import:
       defaults
-    , testDependencies
+    , commonTestDependencies
   type: exitcode-stdio-1.0
   main-is: drivers/hspec/Main.hs
   build-depends: process ^>= 1.6.18.0
@@ -889,11 +906,16 @@ benchmark garg-bench
   other-modules:
       Paths_gargantext
   build-depends: base
-               , bytestring ^>= 0.11.5.3
+               , accelerate
+               , accelerate-llvm-native
+               , hmatrix
+               , massiv
+               , deepseq
                , gargantext
                , gargantext-prelude
+               , split
                , tasty-bench
-  ghc-options:   "-with-rtsopts=-T -A32m"
+  ghc-options:   -threaded "-with-rtsopts=-N -T -A32m"
   if impl(ghc >= 8.6)
-    ghc-options: "-with-rtsopts=--nonmoving-gc"
+    ghc-options: -threaded "-with-rtsopts=-N --nonmoving-gc"
 

src/Gargantext/Core/LinearAlgebra.hs

+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
+{-# LANGUAGE TypeOperators #-}
+{-|
+Module      : Gargantext.Core.LinearAlgebra
+Description : Linear Algebra utility functions
+Copyright   : (c) CNRS, 2017-Present
+License     : AGPL + CECILL v3
+Maintainer  : team@gargantext.org
+Stability   : experimental
+Portability : POSIX
+
+Linear algebra utility functions to be used across all the Gargantext modules requiring it.
+-}
+
+module Gargantext.Core.LinearAlgebra (
+  -- * Handy re-exports
+    module Gargantext.Core.LinearAlgebra.Operations
+  , module Gargantext.Core.LinearAlgebra.Distributional
+  ) where
+
+import Gargantext.Core.LinearAlgebra.Operations
+import Gargantext.Core.LinearAlgebra.Distributional

src/Gargantext/Core/LinearAlgebra/Distributional.hs

+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
+{-|
+Module      : Gargantext.Core.LinearAlgebra.Distributional
+Description : The "distributional" algorithm, fast and slow implementations
+Copyright   : (c) CNRS, 2017-Present
+License     : AGPL + CECILL v3
+Maintainer  : team@gargantext.org
+Stability   : experimental
+Portability : POSIX
+
+-}
+
+module Gargantext.Core.LinearAlgebra.Distributional (
+    distributional
+  , logDistributional2
+
+  -- * Internals for testing
+  , distributionalReferenceImplementation
+  ) where
+
+import Data.Massiv.Array (D, Matrix, Vector, Array, Ix3, U, Ix2 (..), IxN (..))
+import Data.Massiv.Array qualified as A
+import Gargantext.Core.LinearAlgebra.Operations
+import Prelude
+
+-- | `distributional m` returns the distributional distance between each
+-- pair of terms as a matrix.  The argument m is the matrix $[n_{ij}]_{i,j}$
+-- where $n_{ij}$ is the coocccurrence between term $i$ and term $j$.
+--
+-- ## Basic example with Matrix of size 3: 
+--
+-- >>> theMatrixInt 3
+-- Matrix (Z :. 3 :. 3)
+--   [ 7, 4, 0,
+--     4, 5, 3,
+--     0, 3, 4]
+--
+-- >>> distributional $ theMatrixInt 3
+-- Matrix (Z :. 3 :. 3)
+--   [ 1.0, 0.0, 0.9843749999999999,
+--     0.0, 1.0,                0.0,
+--     1.0, 0.0,                1.0]
+--
+-- ## Basic example with Matrix of size 4: 
+--
+-- >>> theMatrixInt 4
+-- Matrix (Z :. 4 :. 4)
+--   [ 4, 1, 2, 1,
+--     1, 4, 0, 0,
+--     2, 0, 3, 3,
+--     1, 0, 3, 3]
+--
+-- >>> distributional $ theMatrixInt 4
+-- Matrix (Z :. 4 :. 4)
+--   [                  1.0,                   0.0, 0.5714285714285715, 0.8421052631578947,
+--                      0.0,                   1.0,                1.0,                1.0,
+--     8.333333333333333e-2,             4.6875e-2,                1.0,               0.25,
+--       0.3333333333333333, 5.7692307692307696e-2,                1.0,                1.0]
+--
+-- /IMPORTANT/: As this function computes the diagonal matrix in order to carry on the computation
+-- the input has to be a square matrix, or this function will fail at runtime.
+distributional :: forall r e. ( A.Manifest r e
+                              , A.Manifest r Int
+                              , A.Unbox e
+                              , A.Source r Int
+                              , A.Size r
+                              , Ord e
+                              , Fractional e
+                              , Num e
+                              )
+                              => Matrix r Int
+                              -> Matrix U e
+distributional m' = A.computeP result
+ where
+    mD :: Matrix D e
+    mD = A.map fromIntegral m'
+
+    m :: Matrix A.U e
+    m = A.compute mD
+
+    n :: Int
+    n = dim m'
+
+    diag_m :: Vector A.U e
+    diag_m = diag m
+
+    d_1 :: Matrix A.D e
+    d_1 = A.backpermute' (A.Sz2 n n) (\(_ A.:. i) -> i) diag_m
+
+    d_2 :: Matrix A.D e
+    d_2 = A.backpermute' (A.Sz2 n n) (\(i A.:. _) -> i) diag_m
+
+    a :: Matrix D e
+    a = termDivNanD mD d_1
+
+    b :: Matrix D e
+    b = termDivNanD mD d_2
+
+    miDelayed :: Matrix D e
+    miDelayed = a `mulD` b
+
+    miMemo :: Matrix D e
+    miMemo = A.delay (A.compute @U miDelayed)
+
+    w_1 :: Array D Ix3 e
+    w_1 = A.backpermute' (A.Sz3 n n n) (\(x A.:> _y A.:. z) -> x A.:. z) miMemo
+
+    w_2 :: Array D Ix3 e
+    w_2 = A.backpermute' (A.Sz3 n n n) (\(_x A.:> y A.:. z) -> y A.:. z) miMemo
+
+    w' :: Array D Ix3 e
+    w' = A.zipWith min w_1 w_2
+
+    z_1 :: Matrix A.D e
+    z_1 = A.ifoldlWithin' 1 ( \(i :> j :. k) acc w'_val ->
+      let ii_val = if k /= i && k /= j then 1 else 0
+          z1_val = w'_val * ii_val
+      in acc + z1_val
+      ) 0 w'
+
+    z_2 :: Matrix A.D e
+    z_2 = A.ifoldlWithin' 1 ( \(i :> j :. k) acc w1_val ->
+      let ii_val = if k /= i && k /= j then 1 else 0
+          z2_val = w1_val * ii_val
+      in acc + z2_val
+      ) 0 w_1
+
+    result :: Matrix A.D e
+    result = termDivNanD z_1 z_2
+
+-- | A reference implementation for \"distributional\" which is slower but
+-- it's more declarative and can be used to assess the correctness of the
+-- optimised version.
+-- Same proviso about the shape of the matri applies for this function.
+distributionalReferenceImplementation :: forall r e.
+                                      ( A.Manifest r e
+                                      , A.Unbox e
+                                      , A.Source r Int
+                                      , A.Size r
+                                      , Ord e
+                                      , Fractional e
+                                      , Num e
+                                      )
+                                      => Matrix r Int
+                                      -> Matrix r e
+distributionalReferenceImplementation m' = result
+ where
+    mD :: Matrix D e
+    mD = A.map fromIntegral m'
+
+    m :: Matrix A.U e
+    m = A.compute mD
+
+    n :: Int
+    n = dim m'
+
+    -- Computes the diagonal matrix of the input ..
+    diag_m :: Vector A.U e
+    diag_m = diag m
+
+    -- Then we create a matrix that contains the same elements of diag_m
+    -- for the rows and columns, to make it square again.
+    d_1 :: Matrix A.D e
+    d_1 = A.backpermute' (A.Sz2 n n) (\(_ A.:. i) -> i) diag_m
+
+    d_2 :: Matrix A.D e
+    d_2 = A.backpermute' (A.Sz2 n n) (\(i A.:. _) -> i) diag_m
+
+    a :: Matrix D e
+    a = termDivNanD mD d_1
+
+    b :: Matrix D e
+    b = termDivNanD mD d_2
+
+    miDelayed :: Matrix D e
+    miDelayed = a `mulD` b
+
+    miMemo :: Matrix D e
+    miMemo = A.delay (A.compute @U miDelayed)
+
+    -- The matrix permutations is taken care of below by directly replicating
+    -- the matrix mi, making the matrix w unneccessary and saving one step.
+    -- replicate (constant (Z :. All :. n :. All)) mi
+    w_1 :: Array D Ix3 e
+    w_1 = A.backpermute' (A.Sz3 n n n) (\(x A.:> _y A.:. z) -> x A.:. z) miMemo
+
+    -- replicate (constant (Z :. n :. All :. All)) mi
+    w_2 :: Array D Ix3 e
+    w_2 = A.backpermute' (A.Sz3 n n n) (\(_x A.:> y A.:. z) -> y A.:. z) miMemo
+
+    w' :: Array D Ix3 e
+    w' = A.zipWith min w_1 w_2
+
+    -- The matrix ii = [r_{i,j,k}]_{i,j,k} has r_(i,j,k) = 0 if k = i OR k = j
+    -- and r_(i,j,k) = 1 otherwise (i.e. k /= i AND k /= j).
+    -- generate (constant (Z :. n :. n :. n)) (lift1 (\( i A.:. j A.:. k) -> cond ((&&) ((/=) k i) ((/=) k j)) 1 0))
+    ii :: Array A.D Ix3 e
+    ii = A.makeArrayR A.D A.Seq (A.Sz3 n n n) $ \(i A.:> j A.:. k) -> if k /= i && k /= j then 1 else 0
+
+    z_1 :: Matrix A.D e
+    z_1 = sumRowsD (w' `mulD` ii)
+
+    z_2 :: Matrix A.D e
+    z_2 = sumRowsD (w_1 `mulD` ii)
+
+    result = A.computeP (termDivNanD z_1 z_2)
+
+
+logDistributional2 :: (A.Manifest r e
+                      , A.Unbox e
+                      , A.Source r Int
+                      , A.Shape r Ix2
+                      , Num e
+                      , Ord e
+                      , A.Source r e
+                      , Fractional e
+                      , Floating e
+                      )
+                   => Matrix r Int
+                   -> Matrix r e
+logDistributional2 m = A.computeP
+                     $ diagNull n
+                     $ matMaxMini
+                     $ logDistributional' n m
+  where
+    n = dim m
+
+logDistributional' :: forall r e.
+                   ( A.Manifest r e
+                   , A.Unbox e
+                   , A.Source r Int
+                   , A.Shape r Ix2
+                   , Num e
+                   , Ord e
+                   , A.Source r e
+                   , Fractional e
+                   , Floating e
+                   )
+                   => Int
+                   -> Matrix r Int
+                   -> Matrix r e
+logDistributional' n m' = result
+ where
+    m :: Matrix A.U e
+    m = A.compute $ A.map fromIntegral m'
+
+    -- Scalar. Sum of all elements of m.
+    to :: e
+    to = A.sum m
+
+    -- Diagonal matrix with the diagonal of m.
+    d_m :: Matrix A.D e
+    d_m = m `mulD` (matrixIdentity n)
+
+    -- Size n vector. s = [s_i]_i
+    s :: Vector A.U e
+    s = A.compute $ sumRowsD (m `subD` d_m)
+
+    -- Matrix nxn. Vector s replicated as rows.
+    s_1 :: Matrix D e
+    s_1 = A.backpermute' (A.Sz2 n n) (\(x :. _y) -> x) s
+
+    -- Matrix nxn. Vector s replicated as columns.
+    s_2 :: Matrix D e
+    s_2 = A.backpermute' (A.Sz2 n n) (\(_x :. y) -> y) s
+
+    -- Matrix nxn. ss = [s_i * s_j]_{i,j}. Outer product of s with itself.
+    ss :: Matrix A.D e
+    ss = s_1 `mulD` s_2
+
+    mi_divvy :: Matrix A.D e
+    mi_divvy = A.zipWith (\m_val ss_val ->
+      let x  = m_val `safeDiv` ss_val
+          x' = x * to
+      in if (x' < 1) then 0 else log x') m ss
+
+    -- Matrix nxn. mi = [m_{i,j}]_{i,j} where
+    -- m_{i,j} = 0 if n_{i,j} = 0 or i = j,
+    -- m_{i,j} = log(to * n_{i,j} / s_{i,j}) otherwise.
+    mi :: Matrix A.U e
+    mi = A.computeP $ mulD (matrixEye n) (mi_divvy)
+
+    sumMin :: Matrix A.U e
+    sumMin = sumMin_go n mi
+
+    sumM :: Matrix A.U e
+    sumM = sumM_go n mi
+
+    result :: Matrix r e
+    result = termDivNan sumMin sumM

src/Gargantext/Core/LinearAlgebra/Operations.hs

+{-# LANGUAGE DerivingStrategies #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeApplications #-}
+{-# LANGUAGE TypeOperators #-}
+{-|
+Module      : Gargantext.Core.LinearAlgebra.Operations
+Description : Operations on matrixes using massiv
+Copyright   : (c) CNRS, 2017-Present
+License     : AGPL + CECILL v3
+Maintainer  : team@gargantext.org
+Stability   : experimental
+Portability : POSIX
+
+-}
+
+module Gargantext.Core.LinearAlgebra.Operations (
+  -- * Convertion functions
+    accelerate2MassivMatrix
+  , accelerate2Massiv3DMatrix
+  , massiv2AccelerateMatrix
+  , massiv2AccelerateVector
+
+  -- * Operations on matrixes
+  , (.*)
+  , (.-)
+  , diag
+  , termDivNan
+  , sumRows
+  , dim
+  , matrixEye
+  , matrixIdentity
+  , diagNull
+
+  -- * Operations on delayed arrays
+  , diagD
+  , subD
+  , mulD
+  , termDivNanD
+  , sumRowsD
+  , safeDiv
+
+  -- * Internals for testing
+  , sumRowsReferenceImplementation
+  , matMaxMini
+  , sumM_go
+  , sumMin_go
+  ) where
+
+import Data.Array.Accelerate qualified as Acc
+import Data.List.Split qualified as Split
+import Data.Massiv.Array (D, Matrix, Vector, Array)
+import Data.Massiv.Array qualified as A
+import Prelude
+import Protolude.Safe (headMay)
+import Data.Monoid
+
+-- | Converts an accelerate matrix into a Massiv matrix.
+accelerate2MassivMatrix :: (A.Unbox a, Acc.Elt a) => Acc.Matrix a -> Matrix A.U a
+accelerate2MassivMatrix m =
+  let (Acc.Z Acc.:. _r Acc.:. c) = Acc.arrayShape m
+  in A.fromLists' @A.U A.Par $ Split.chunksOf c (Acc.toList m)
+
+-- | Converts a massiv matrix into an accelerate matrix.
+massiv2AccelerateMatrix :: (Acc.Elt a, A.Source r a) => Matrix r a -> Acc.Matrix a
+massiv2AccelerateMatrix m =
+  let m' = A.toLists2 m
+      r  = Prelude.length m'
+      c  = maybe 0 Prelude.length (headMay m')
+  in Acc.fromList (Acc.Z Acc.:. r Acc.:. c) (mconcat m')
+
+-- | Converts a massiv vector into an accelerate one.
+massiv2AccelerateVector :: (A.Source r a, Acc.Elt a) => A.Vector r a -> Acc.Vector a
+massiv2AccelerateVector m =
+  let m' = A.toList m
+      r  = Prelude.length m'
+  in Acc.fromList (Acc.Z Acc.:. r) m'
+
+accelerate2Massiv3DMatrix :: (A.Unbox e, Acc.Elt e, A.Manifest r e)
+                          => Acc.Array (Acc.Z Acc.:. Int Acc.:. Int Acc.:. Int) e
+                          -> A.Array r A.Ix3 e
+accelerate2Massiv3DMatrix m =
+  let (Acc.Z Acc.:. _r Acc.:. _c Acc.:. _z) = Acc.arrayShape m
+  in A.fromLists' A.Par $ map (Split.chunksOf $ _z) $ Split.chunksOf (_c*_z) (Acc.toList m)
+
+
+
+-- | Computes the diagnonal matrix of the input one.
+diag :: (A.Unbox e, A.Manifest r e, A.Source r e, Num e) => Matrix r e -> Vector A.U e
+diag matrix =
+  let (A.Sz2 rows _cols) = A.size matrix
+      newSize = A.Sz1 rows
+  in A.makeArrayR A.U A.Seq newSize $ (\(A.Ix1 i) -> matrix A.! (A.Ix2 i i))
+
+diagD :: (A.Source r e, A.Size r) => Matrix r e -> Vector A.D e
+diagD matrix =
+  let (A.Sz2 rows _cols) = A.size matrix
+      newSize = A.Sz1 rows
+  in A.backpermute' newSize (\i -> i A.:. i) matrix
+
+-- | Term by term division where divisions by 0 produce 0 rather than NaN.
+termDivNan :: (A.Manifest r3 a, A.Source r1 a, A.Source r2 a, Eq a, Fractional a)
+           => Matrix r1 a
+           -> Matrix r2 a
+           -> Matrix r3 a
+termDivNan m1 = A.compute . termDivNanD m1
+
+termDivNanD :: (A.Source r1 a, A.Source r2 a, Eq a, Fractional a)
+            => Matrix r1 a
+            -> Matrix r2 a
+            -> Matrix D a
+termDivNanD m1 m2 = A.zipWith safeDiv m1 m2
+
+safeDiv :: (Eq a, Fractional a) => a -> a -> a
+safeDiv i j = if j == 0 then 0 else i / j
+{-# INLINE safeDiv #-}
+
+sumRows :: ( A.Index (A.Lower ix)
+           , A.Index ix
+           , A.Source r e
+           , A.Manifest r e
+           , A.Strategy r
+           , A.Size r
+           , Num e
+           ) => Array r ix e
+             -> Array r (A.Lower ix) e
+sumRows = A.compute . sumRowsD
+
+sumRowsD :: ( A.Index (A.Lower ix)
+            , A.Index ix
+            , A.Source r e
+            , Num e
+            ) => Array r ix e
+              -> Array D (A.Lower ix) e
+sumRowsD matrix = A.map getSum $ A.foldlWithin' 1 (\(Sum s) n -> Sum $ s + n) mempty matrix
+
+sumRowsReferenceImplementation :: ( A.Load r A.Ix2 e
+                                  , A.Source r e
+                                  , A.Manifest r e
+                                  , A.Strategy r
+                                  , A.Size r
+                                  , Num e
+                                  ) => Array r A.Ix3 e
+                                    -> Array r A.Ix2 e
+sumRowsReferenceImplementation matrix =
+  let A.Sz3 rows cols z = A.size matrix
+  in A.makeArray (A.getComp matrix) (A.Sz2 rows cols) $ \(i A.:. j) ->
+       A.sum (A.backpermute' (A.Sz1 z) (\c -> i A.:> j A.:. c) matrix)
+
+-- | Matrix cell by cell multiplication
+(.*) :: (A.Manifest r3 a, A.Source r1 a, A.Source r2 a, A.Index ix, Num a)
+     => Array r1 ix a
+     -> Array r2 ix a
+     -> Array r3 ix a
+(.*) m1 = A.compute . mulD m1
+
+mulD :: (A.Source r1 a, A.Source r2 a, A.Index ix, Num a)
+     => Array r1 ix a
+     -> Array r2 ix a
+     -> Array D ix a
+mulD m1 m2 = A.zipWith (*) m1 m2
+
+-- | Matrix cell by cell substraction
+(.-) :: (A.Manifest r3 a, A.Source r1 a, A.Source r2 a, A.Index ix, Num a)
+     => Array r1 ix a
+     -> Array r2 ix a
+     -> Array r3 ix a
+(.-) m1 = A.compute . subD m1
+
+subD :: (A.Source r1 a, A.Source r2 a, A.Index ix, Num a)
+     => Array r1 ix a
+     -> Array r2 ix a
+     -> Array D ix a
+subD m1 m2 = A.zipWith (-) m1 m2
+
+
+-- | Get the dimensions of a /square/ matrix.
+dim :: A.Size r => Matrix r a -> Int
+dim m = n
+  where
+    (A.Sz2 _ n) = A.size 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.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.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))
+        else 0
+    | k <- [0 .. n - 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 :: (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.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)

src/Gargantext/Core/Methods/Matrix/Accelerate/Utils.hs

@@ -38,8 +38,6 @@ import Data.Array.Accelerate
 import Data.Array.Accelerate.Interpreter (run)
 import qualified Gargantext.Prelude as P
 
-import Debug.Trace (trace)
-
 -- | Matrix cell by cell multiplication
 (.*) :: ( Shape ix
         , Slice ix
@@ -55,26 +53,28 @@ import Debug.Trace (trace)
 (./) :: ( Shape ix
         , Slice ix
         , Elt a
+        , Eq a
         , P.Num (Exp a)
         , P.Fractional (Exp a)
         )
      => Acc (Array ((ix :. Int) :. Int) a)
      -> Acc (Array ((ix :. Int) :. Int) a)
      -> Acc (Array ((ix :. Int) :. Int) a)
-(./) = zipWith (/)
+(./) = zipWith safeDivCond
 
 -- | Term by term division where divisions by 0 produce 0 rather than NaN. 
-termDivNan :: ( Shape ix
-        , Slice ix
-        , Elt a
+termDivNan :: ( Elt a
               , Eq a
               , P.Num (Exp a)
               , P.Fractional (Exp a)
               )
-     => Acc (Array ((ix :. Int) :. Int) a)
-     -> Acc (Array ((ix :. Int) :. Int) a)
-     -> Acc (Array ((ix :. Int) :. Int) a)
-termDivNan = trace "termDivNan" $ zipWith (\i j -> cond ((==) j 0) 0 ((/) i j))
+           => Acc (Matrix a)
+           -> Acc (Matrix a)
+           -> Acc (Matrix a)
+termDivNan = zipWith safeDivCond
+
+safeDivCond :: (Eq a, P.Num (Exp a), P.Fractional (Exp a)) => Exp a -> Exp a -> Exp a
+safeDivCond i j = cond ((==) j 0) 0 ((/) i j)
 
 (.-) :: ( Shape ix
         , Slice ix

src/Gargantext/Core/Methods/Similarities.hs

@@ -20,10 +20,11 @@ import Data.Swagger
 import Data.Text (Text)
 import GHC.Generics (Generic)
 import Gargantext.Core.Methods.Similarities.Accelerate.Conditional (measureConditional')
-import Gargantext.Core.Methods.Similarities.Accelerate.Distributional (logDistributional2)
+import Gargantext.Core.LinearAlgebra.Operations (accelerate2MassivMatrix, massiv2AccelerateMatrix)
+import Gargantext.Core.LinearAlgebra.Distributional (logDistributional2)
 -- import Gargantext.Core.Text.Metrics.Count (coocOn)
 -- import Gargantext.Core.Viz.Graph.Index
-import Gargantext.Prelude (Ord, Eq, Int, Double, Show, map)
+import Gargantext.Prelude (Ord, Eq, Int, Double, Show, map, ($), (.))
 import Prelude (Enum, Bounded, minBound, maxBound)
 import Test.QuickCheck (elements)
 import Test.QuickCheck.Arbitrary
@@ -36,7 +37,7 @@ data Similarity = Conditional | Distributional
 
 measure :: Similarity -> Matrix Int -> Matrix Double
 measure Conditional    x = measureConditional' x
-measure Distributional x = logDistributional2 x
+measure Distributional x = massiv2AccelerateMatrix . logDistributional2 . accelerate2MassivMatrix $ x
 
 ------------------------------------------------------------------------
 withMetric :: GraphMetric -> Similarity

src/Gargantext/Core/Methods/Similarities/Accelerate/Distributional.hs

@@ -89,18 +89,22 @@ where $n_{ij}$ is the cooccurrence between term $i$ and term $j$
 {-# LANGUAGE ViewPatterns        #-}
 
 module Gargantext.Core.Methods.Similarities.Accelerate.Distributional
+  ( distributional
+  , logDistributional2
+
+  -- internals for testing
+  , distributionalWith
+  , logDistributional2With
+  , sumMin_go
+  , sumM_go
+  )
   where
 
 -- import qualified Data.Foldable as P (foldl1)
 -- import Debug.Trace (trace)
 import Data.Array.Accelerate as A
--- import Data.Array.Accelerate.Interpreter (run)
-import Data.Array.Accelerate.LLVM.Native (run) -- TODO: try runQ?
+import Data.Array.Accelerate.Interpreter qualified as Naive
 import Gargantext.Core.Methods.Matrix.Accelerate.Utils
-import qualified Gargantext.Prelude as P
-
-import Debug.Trace
-import Prelude (show, mappend{- , String, (<>), fromIntegral, flip -})
 
 import qualified Prelude
 
@@ -138,8 +142,16 @@ import qualified Prelude
 --     8.333333333333333e-2,             4.6875e-2,                1.0,               0.25,
 --       0.3333333333333333, 5.7692307692307696e-2,                1.0,                1.0]
 --
+-- /IMPORTANT/: As this function computes the diagonal matrix in order to carry on the computation
+-- the input has to be a square matrix, or this function will fail at runtime.
 distributional :: Matrix Int -> Matrix Double
-distributional m' = run $ result
+distributional = distributionalWith Naive.run
+
+distributionalWith :: (Elt e, FromIntegral Int e, Eq e, Prelude.Fractional (Exp e), Ord e)
+                   => (forall a. Arrays a => Acc a -> a)
+                   -> Matrix Int
+                   -> Matrix e
+distributionalWith interpret m' = interpret $ result
  where
     m = map A.fromIntegral $ use m'
     n = dim m'
@@ -149,7 +161,7 @@ distributional m' = run $ result
     d_1 = replicate (constant (Z :. n :. All)) diag_m
     d_2 = replicate (constant (Z :. All :. n)) diag_m
 
-    mi    = (.*) ((./) m d_1) ((./) m d_2)
+    mi    = (.*) (termDivNan m d_1) (termDivNan m d_2)
 
     -- w = (.-) mi d_mi
     
@@ -170,15 +182,36 @@ distributional m' = run $ result
     result = termDivNan z_1 z_2
 
 logDistributional2 :: Matrix Int -> Matrix Double
-logDistributional2 m = trace ("logDistributional2, dim=" `mappend` show n) . run
+logDistributional2 m = logDistributional2With Naive.run m
+
+logDistributional2With :: ( Elt e
+                          , Prelude.Num (Exp e)
+                          , Ord e
+                          , Prelude.Num e
+                          , FromIntegral Int e
+                          , Prelude.Fractional (Exp e)
+                          , Prelude.Floating (Exp e)
+                          )
+                       => (forall a. Arrays a => Acc a -> a)
+                       -> Matrix Int -> Matrix e
+logDistributional2With interpreter m = interpreter
                      $ diagNull n
                      $ matMaxMini
                      $ logDistributional' n m
   where
     n = dim m
 
-logDistributional' :: Int -> Matrix Int -> Acc (Matrix Double)
-logDistributional' n m' = trace ("logDistributional'") result
+logDistributional' :: ( Elt e
+                      , Prelude.Num (Exp e)
+                      , FromIntegral Int e
+                      , Eq e
+                      , Ord e
+                      , Prelude.Fractional (Exp e)
+                      , Prelude.Floating (Exp e)
+                      )
+                   => Int -> Matrix Int
+                   -> Acc (Matrix e)
+logDistributional' n m' = result
  where
     -- From Matrix Int to Matrix Double, i.e :
     -- m :: Matrix Int -> Matrix Double
@@ -236,10 +269,10 @@ logDistributional' n m' = trace ("logDistributional'") result
     -- k_diff_i_and_j = lift1 (\(Z :. i :. j :. k) -> ((&&) ((/=) k i) ((/=) k j)))
 
     -- Matrix nxn. 
-    sumMin = trace "sumMin" $ sumMin_go n mi -- sum (condOrDefault k_diff_i_and_j 0 w')
+    sumMin = sumMin_go n mi -- sum (condOrDefault k_diff_i_and_j 0 w')
 
     -- Matrix nxn. All columns are the same.
-    sumM = trace "sumM" $ sumM_go n mi -- trace "sumM" $ sum (condOrDefault k_diff_i_and_j 0 w_1)
+    sumM = sumM_go n mi -- trace "sumM" $ sum (condOrDefault k_diff_i_and_j 0 w_1)
 
     result = termDivNan sumMin sumM
 
@@ -264,103 +297,6 @@ logDistributional' n m' = trace ("logDistributional'") result
 --            \[N_{m} = \sum_{i,i \neq i}^{m} \sum_{j, j \neq j}^{m} S_{ij}\]
 --
 
-logDistributional :: Matrix Int -> Matrix Double
-logDistributional m' = run $ diagNull n $ result
- where
-    m = map fromIntegral $ use m'
-    n = dim m'
-
-    -- Scalar. Sum of all elements of m.
-    to = the $ sum (flatten m)
-
-    -- Diagonal matrix with the diagonal of m.
-    d_m = (.*) m (matrixIdentity n)
-
-    -- Size n vector. s = [s_i]_i
-    s = sum ((.-) m d_m)
-
-    -- Matrix nxn. Vector s replicated as rows. 
-    s_1 = replicate (constant (Z :. All :. n)) s
-    -- Matrix nxn. Vector s replicated as columns. 
-    s_2 = replicate (constant (Z :. n :. All)) s
-
-    -- Matrix nxn. ss = [s_i * s_j]_{i,j}. Outer product of s with itself. 
-    ss = (.*) s_1 s_2
-
-    -- Matrix nxn. mi = [m_{i,j}]_{i,j} where 
-    -- m_{i,j} = 0 if n_{i,j} = 0 or i = j, 
-    -- m_{i,j} = log(to * n_{i,j} / s_{i,j}) otherwise.
-    mi = (.*) (matrixEye n) 
-        (map (lift1 (\x -> cond (x == 0) 0 (log (x * to)))) ((./) m ss))
-
-    -- Tensor nxnxn. Matrix mi replicated along the 2nd axis.
-    w_1 = replicate (constant (Z :. All :. n :. All)) mi
-
-    -- Tensor nxnxn. Matrix mi replicated along the 1st axis.
-    w_2 = replicate (constant (Z :. n :. All :. All)) mi
-
-    -- Tensor nxnxn.
-    w' = zipWith min w_1 w_2
-
-    -- A predicate that is true when the input (i, j, k) satisfy 
-    -- k /= i AND k /= j
-    k_diff_i_and_j = lift1 (\(Z :. i :. j :. k) -> ((&&) ((/=) k i) ((/=) k j)))
-
-    -- Matrix nxn. 
-    sumMin = sum (condOrDefault k_diff_i_and_j 0 w')
-
-    -- Matrix nxn. All columns are the same.
-    sumM = sum (condOrDefault k_diff_i_and_j 0 w_1)
-
-    result = termDivNan sumMin sumM
-
-
-
-
-distributional'' :: Matrix Int -> Matrix Double
-distributional'' m = -- run {- $ matMaxMini -}
-                   run  $ diagNull n
-                       $ rIJ n
-                       $ filterWith 0 100
-                       $ filter' 0
-                       $ s_mi
-                       $ map A.fromIntegral
-                          {- from Int to Double -}
-                       $ use m
-                          {- push matrix in Accelerate type -}
-  where
-
-    _ri :: Acc (Matrix Double) -> Acc (Matrix Double)
-    _ri mat = mat1 -- zipWith (/) mat1 mat2
-      where
-        mat1 = matSumCol n $ zipWith min (_myMin mat) (_myMin $ filterWith 0 100 $ diagNull n $ transpose mat)
-        _mat2 = total mat
-
-    _myMin :: Acc (Matrix Double) -> Acc (Matrix Double)
-    _myMin = replicate (constant (Z :. n :. All)) . minimum
-
-
-    -- TODO fix NaN
-    -- Quali TEST: OK
-    s_mi :: Acc (Matrix Double) -> Acc (Matrix Double)
-    s_mi m' = zipWith (\x y -> log (x / y)) (diagNull n m')
-            $ zipWith (/) (crossProduct n m') (total m')
-            -- crossProduct n m'
-
-
-    total :: Acc (Matrix Double) -> Acc (Matrix Double)
-    total = replicate (constant (Z :. n :. n)) . sum . sum
-
-    n :: Dim
-    n = dim m
-
-rIJ :: (Elt a, Ord a, P.Fractional (Exp a), P.Num a)
-    => Dim -> Acc (Matrix a) -> Acc (Matrix a)
-rIJ n m = matMaxMini $ divide a b
-  where
-    a = sumRowMin n m
-    b = sumColMin n m
-
 -- * For Tests (to be removed)
 -- | Test perfermance with this matrix
 -- TODO : add this in a benchmark folder
@@ -376,25 +312,6 @@ distriTest n = logDistributional m == distributional m
 -- compact repr of "extend along an axis" op?
 -- general sparse repr ?
 
-type Extended sh = sh :. Int
-
-data Ext where
-  Along1 :: Int -> Ext
-  Along2 :: Int -> Ext
-
-along1 :: Int -> Ext
-along1 = Along1
-
-along2 :: Int -> Ext
-along2 = Along2
-
-type Delayed sh a = Exp sh -> Exp a
-
-data ExtArr sh a = ExtArr
-  { extSh  :: Extended sh
-  , extFun :: Delayed (Extended sh) a
-  }
-
 {-
 w_1_{i, j, k} = mi_{i, k}
 w_2_{i, j, k} = mi_{j, k}

src/Gargantext/Core/Viz/Graph/Index.hs

@@ -20,6 +20,7 @@ TODO:
 {-# LANGUAGE BangPatterns      #-}
 {-# LANGUAGE TypeOperators     #-}
 {-# LANGUAGE MonoLocalBinds    #-}
+{-# OPTIONS_GHC -Wno-unused-top-binds #-}
 
 module Gargantext.Core.Viz.Graph.Index
   where

src/Gargantext/Orphans.hs

@@ -5,6 +5,7 @@ module Gargantext.Orphans (
 
 import Data.Aeson qualified as JSON
 import Gargantext.Database.Admin.Types.Hyperdata (Hyperdata)
+import Gargantext.Orphans.Accelerate ()
 import Gargantext.Orphans.OpenAPI
 
 instance Hyperdata JSON.Value

src/Gargantext/Orphans/Accelerate.hs

+{-# OPTIONS_GHC -Wno-orphans #-}
+{-# LANGUAGE TypeApplications #-}
+{-# LANGUAGE MultiWayIf #-}
+
+module Gargantext.Orphans.Accelerate where
+
+import Prelude
+import Test.QuickCheck
+import Data.Scientific ()
+import Data.Array.Accelerate (DIM2, Z (..), (:.) (..), Array, Elt, fromList, arrayShape, DIM3)
+import Data.Array.Accelerate qualified as A
+import qualified Data.List.Split as Split
+
+instance (Show e, Elt e, Arbitrary e, Num e, Ord e) => Arbitrary (Array DIM3 e) where
+  arbitrary = do
+    x <- choose (1,10)
+    y <- choose (1,10)
+    z <- choose (1,10)
+    let sh = Z :. x :. y :. z
+    fromList sh <$> vectorOf (x * y * z) (getPositive <$> arbitrary)
+
+instance (Show e, Elt e, Arbitrary e) => Arbitrary (Array DIM2 e) where
+  arbitrary = do
+    x <- choose (1,128)
+    y <- choose (1,48)
+    let sh = Z :. x :. y
+    fromList sh <$> vectorOf (x * y) arbitrary
+  shrink arr = sliceArray arr
+
+-- Slice the array to the new shape, keeping the square dimensions.
+sliceArray :: (Elt e, Show e) => Array DIM2 e -> [Array DIM2 e]
+sliceArray arr =
+  case arrayShape arr of
+    (Z :. x :. y) -> case (x, y) of
+      (_,1) -> [ ]
+      (1,_) -> [ ]
+      _     -> [ resizeArray arr (max 1 (x - 1)) (max 1 (y - 1)) ]
+
+resizeArray :: (Show e, Elt e) => Array DIM2 e -> Int -> Int -> Array DIM2 e
+resizeArray arr rows cols =
+  let (Z :. _originRows :. originCols) = arrayShape arr
+      vals   = A.toList arr
+      chunks = map (take cols) $ Split.chunksOf originCols vals
+      m'     = mconcat $ take rows chunks
+  in A.fromList (Z :. rows :. cols) m'

stack.yaml

@@ -440,6 +440,7 @@ flags:
   formatting:
     "no-double-conversion": false
   gargantext:
+    "enable-benchmarks": false
     "no-phylo-debug-logs": true
     "test-crypto": false
   graphviz:
@@ -509,6 +510,8 @@ flags:
     "shared-llvm": true
   lzma:
     pkgconfig: true
+  massiv:
+    "unsafe-checks": false
   "math-functions":
     "system-erf": true
     "system-expm1": true

test/drivers/tasty/Main.hs

@@ -12,6 +12,7 @@ module Main where
 
 import Gargantext.Prelude
 
+import qualified Test.Core.LinearAlgebra         as LinearAlgebra
 import qualified Test.Core.Notifications         as Notifications
 import qualified Test.Core.Orchestrator          as Orchestrator
 import qualified Test.Core.Similarity            as Similarity
@@ -78,4 +79,5 @@ main = do
     , Notifications.qcTests
     , Orchestrator.qcTests
     , NgramsTerms.tests
+    , LinearAlgebra.tests
     ]