WIP: start the implementation of logDistributional

bench/Main.hs

@@ -21,6 +21,7 @@ import qualified Data.Array.Accelerate.Interpreter 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 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
@@ -99,6 +100,16 @@ main = do
       ,   bench "Accelerate (LLVM)"  $ nf (LLVM.run  . Accelerate.diag . Accelerate.use) accInput
       ,   bench "Massiv " $ nf (LA.diag @_) massivInput
       ]
+      , bgroup "identityMatrix" [
+          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
+      ]
+      , 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
+      ]
       , 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

src/Gargantext/Core/LinearAlgebra.hs

@@ -29,13 +29,18 @@ module Gargantext.Core.LinearAlgebra (
 
   -- * Operations on matrixes
   , (.*)
+  , (.-)
   , diag
   , termDivNan
   , distributional
+  --, logDistributional
   , sumRows
 
   -- * Internals for testing
   , sumRowsReferenceImplementation
+  , matMaxMini
+  , sumM_go
+  , sumMin_go
   ) where
 
 import Data.Array.Accelerate qualified as Acc
@@ -151,12 +156,12 @@ distributional :: forall r e.
 distributional m' = result
  where
     m :: Matrix A.U e
-    m = A.computeP $ A.map fromIntegral m'
+    m = A.compute $ A.map fromIntegral m'
     n = dim m'
 
     -- Computes the diagonal matrix of the input ..
     diag_m :: Vector A.U e
-    diag_m = A.computeP $ diag m
+    diag_m = diag m
 
     -- .. and its size.
     diag_m_size :: Int
@@ -252,8 +257,104 @@ mulD :: (A.Source r1 a, A.Source r2 a, A.Index ix, Num 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.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
+  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))
+
+{-
+logDistributional :: forall r e. (
+                     )
+                     => Matrix r Int
+                     -> Matrix r e
+logDistributional m =
+    interpreter
+  $ diagNull n
+  $ matMaxMini
+  $ logDistributional' n m
+  where
+    n = dim m
+
+logDistributional' :: forall r e. ( )
+                   => Int
+                   -> Matrix r Int
+                   -> Matrix r e
+logDistributional' n m' = result
+ where
+
+    m :: Matrix A.U Int
+    m = A.compute $ A.map fromIntegral m'
+
+    -- Scalar. Sum of all elements of m.
+    to = sum m
+
+    -- Diagonal matrix with the diagonal of m.
+    d_m = (.*) m (A.identityMatrix (A.Sz1 n))
+
+    -- Size n vector. s = [s_i]_i
+    s = sum ((.-) m d_m)
+
+    -- Matrix nxn. Vector s replicated as rows.
+    s_1 :: Matrix A.U Int
+    s_1 = A.makeArrayR A.U A.Seq (A.Sz2 n n) $ \(_ A.:. i) -> s A.! i
+
+    -- Matrix nxn. Vector s replicated as columns.
+    s_2 :: Matrix A.U Int
+    s_2 = A.makeArrayR A.U A.Seq (A.Sz2 n n) $ \(i A.:. _) -> s A.! i
+
+    -- 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 -> let x' = x * to in cond (x' < 0.5) 0 (log x'))) ((./) m ss))
+              (map (lift1 (\x -> let x' = x * to in cond (x' < 1) 0 (log x'))) ((./) m ss))
+
+    -- Matrix nxn.
+    sumMin :: Matrix A.U e
+    sumMin = sumMin_go n mi
+
+    -- Matrix nxn. All columns are the same.
+    sumM :: Matrix A.U e
+    sumM = sumM_go n mi
+
+    result = termDivNan sumMin sumM
+-}
+
+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) ->
+  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) ->
+  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]
+    ]
+

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

@@ -95,6 +95,8 @@ module Gargantext.Core.Methods.Similarities.Accelerate.Distributional
   -- internals for testing
   , distributionalWith
   , logDistributional2With
+  , sumMin_go
+  , sumM_go
   )
   where
 

test/Test/Core/LinearAlgebra.hs

@@ -8,23 +8,24 @@ module Test.Core.LinearAlgebra where
 
 import Data.Array.Accelerate hiding (Ord, Eq, map, (<=))
 import Data.Array.Accelerate.Interpreter qualified as Naive
+import Data.Array.Accelerate qualified as A
 import Data.Bifunctor (first)
 import Data.Bimap (Bimap)
 import Data.Bimap qualified as Bimap
 import Data.Map.Strict (Map)
 import Data.Map.Strict qualified as M
 import Data.Massiv.Array qualified as Massiv
+import Data.Proxy
+import Data.Scientific
 import Gargantext.Core.LinearAlgebra qualified as LA
+import Gargantext.Core.Methods.Matrix.Accelerate.Utils qualified as A
 import Gargantext.Core.Methods.Matrix.Accelerate.Utils qualified as Legacy
 import Gargantext.Core.Methods.Similarities.Accelerate.Distributional qualified as Legacy
 import Gargantext.Core.Viz.Graph.Index qualified as Legacy
 import Gargantext.Orphans.Accelerate (sliceArray)
 import Prelude hiding ((^))
-import qualified Data.Array.Accelerate as A
 import Test.Tasty
 import Test.Tasty.QuickCheck
-import Data.Proxy
-import Data.Scientific
 
 
 --
@@ -107,6 +108,9 @@ tests = testGroup "LinearAlgebra" [
     , testProperty "termDivNan"     compareTermDivNan
     , testProperty "diag"           compareDiag
     , testProperty "sumRows"        compareSumRows
+    , testProperty "matMaxMini"     compareMatMaxMini
+    , testProperty "sumM_go"        compareSumM_go
+    , testProperty "sumMin_go"      compareSumMin_go
     , testGroup "distributional" [
         testProperty "2x2" (compareDistributional (Proxy @Double) twoByTwo)
       , testProperty "7x7" (compareDistributional (Proxy @Double) testMatrix_02)
@@ -164,3 +168,21 @@ compareDistributional Proxy (SquareMatrix i1)
   where
     conv :: e -> Scientific
     conv = fromFloatDigits
+
+compareMatMaxMini :: SquareMatrix Int -> Property
+compareMatMaxMini (SquareMatrix i1)
+  = let massiv     = LA.matMaxMini @Massiv.U (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)
+        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)
+        accelerate =  Naive.run (Legacy.sumM_go (A.dim i1) (use i1))
+    in accelerate === LA.massiv2AccelerateMatrix massiv