WIP: Implement distributional2 using massiv

There are some rounding errors in the tests, but otherwise
the implementation seems correct.

bench/Main.hs

@@ -120,6 +120,11 @@ main = do
       ,   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
@@ -136,6 +141,7 @@ main = do
       , 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
       ]
       ]
     ]

src/Gargantext/Core/LinearAlgebra/Distributional.hs

@@ -15,6 +15,7 @@ Portability : POSIX
 
 module Gargantext.Core.LinearAlgebra.Distributional (
     distributional
+  , logDistributional2
 
   -- * Internals for testing
   , distributionalReferenceImplementation
@@ -83,21 +84,14 @@ distributional m' = A.computeP result
     n :: Int
     n = dim m'
 
-    -- Computes the diagonal matrix of the input ..
     diag_m :: Vector A.U e
     diag_m = diag m
 
-    -- .. and its size.
-    diag_m_size :: Int
-    (A.Sz1 diag_m_size) = A.size 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 diag_m_size) (\(_ A.:. i) -> i) diag_m
+    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 diag_m_size n) (\(i A.:. _) -> i) diag_m
+    d_2 = A.backpermute' (A.Sz2 n n) (\(i A.:. _) -> i) diag_m
 
     a :: Matrix D e
     a = termDivNanD mD d_1
@@ -111,18 +105,11 @@ distributional m' = A.computeP result
     miMemo :: Matrix D e
     miMemo = A.delay (A.compute @U miDelayed)
 
-    mi_r, mi_c :: Int
-    (A.Sz2 mi_r mi_c) = A.size miMemo
-
-    -- 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 mi_r n mi_c) (\(x A.:> _y A.:. z) -> x A.:. z) miMemo
+    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 mi_r mi_c) (\(_x A.:> y A.:. z) -> y A.:. z) miMemo
+    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
@@ -174,17 +161,13 @@ distributionalReferenceImplementation m' = result
     diag_m :: Vector A.U e
     diag_m = diag m
 
-    -- .. and its size.
-    diag_m_size :: Int
-    (A.Sz1 diag_m_size) = A.size 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 diag_m_size) (\(_ A.:. i) -> i) diag_m
+    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 diag_m_size n) (\(i A.:. _) -> i) diag_m
+    d_2 = A.backpermute' (A.Sz2 n n) (\(i A.:. _) -> i) diag_m
 
     a :: Matrix D e
     a = termDivNanD mD d_1
@@ -198,18 +181,15 @@ distributionalReferenceImplementation m' = result
     miMemo :: Matrix D e
     miMemo = A.delay (A.compute @U miDelayed)
 
-    mi_r, mi_c :: Int
-    (A.Sz2 mi_r mi_c) = A.size miMemo
-
     -- 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 mi_r n mi_c) (\(x A.:> _y A.:. z) -> x A.:. z) miMemo
+    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 mi_r mi_c) (\(_x A.:> y A.:. z) -> y A.:. z) miMemo
+    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
@@ -227,3 +207,86 @@ distributionalReferenceImplementation m' = result
     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.D e
+    m = 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.computeP $ 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 / 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

@@ -24,12 +24,16 @@ module Gargantext.Core.LinearAlgebra.Operations (
   , (.*)
   , (.-)
   , diag
-  , diagD
   , termDivNan
   , sumRows
   , dim
+  , matrixEye
+  , matrixIdentity
+  , diagNull
 
   -- * Operations on delayed arrays
+  , diagD
+  , subD
   , mulD
   , termDivNanD
   , sumRowsD
@@ -110,20 +114,23 @@ safeDiv :: (Eq a, Fractional a) => a -> a -> a
 safeDiv i j = if j == 0 then 0 else i / j
 {-# INLINE safeDiv #-}
 
-sumRows :: ( A.Load r A.Ix2 e
+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 A.Ix3 e
-             -> Array r A.Ix2 e
+           ) => Array r ix e
+             -> Array r (A.Lower ix) e
 sumRows = A.compute . sumRowsD
 
-sumRowsD :: ( A.Source r e
+sumRowsD :: ( A.Index (A.Lower ix)
+            , A.Index ix
+            , A.Source r e
             , Num e
-            ) => Array r A.Ix3 e
-              -> Array D A.Ix2 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
@@ -192,3 +199,11 @@ 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
+
+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
+
+diagNull :: (A.Source r e, Num e) => Int -> Matrix r e -> Matrix D e
+diagNull n m = A.zipWith (*) m (matrixEye n)

test/Test/Core/LinearAlgebra.hs

@@ -3,6 +3,7 @@
 {-# LANGUAGE MultiWayIf #-}
 {-# LANGUAGE DerivingStrategies #-}
 {-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE ViewPatterns #-}
 
 module Test.Core.LinearAlgebra where
 
@@ -111,6 +112,8 @@ tests = testGroup "LinearAlgebra" [
     , testProperty "matMaxMini"     compareMatMaxMini
     , testProperty "sumM_go"        compareSumM_go
     , testProperty "sumMin_go"      compareSumMin_go
+    , testProperty "matrixEye"      compareMatrixEye
+    , testProperty "diagNull"       compareDiagNull
     , testGroup "distributional" [
         testProperty "reference implementation roundtrips" compareDistributionalImplementations
       , testProperty "2x2" (compareDistributional (Proxy @Double) twoByTwo)
@@ -118,6 +121,12 @@ tests = testGroup "LinearAlgebra" [
       , testProperty "14x14" (compareDistributional (Proxy @Double) testMatrix_01)
       , testProperty "roundtrips" (compareDistributional (Proxy @Double))
       ]
+    , testGroup "logDistributional2" [
+       testProperty "2x2" (compareLogDistributional2 (Proxy @Double) twoByTwo)
+      , testProperty "7x7" (compareLogDistributional2 (Proxy @Double) testMatrix_02)
+      , testProperty "14x14" (compareLogDistributional2 (Proxy @Double) testMatrix_01)
+      ,testProperty "roundtrips" (compareLogDistributional2 (Proxy @Double))
+    ]
     ]
 
 --
@@ -140,7 +149,7 @@ compareDiag (SquareMatrix i1)
 
 compareSumRows :: Array (Z :. Int :. Int :. Int) Int -> Property
 compareSumRows i1
-  = let massiv     = LA.sumRows @Massiv.U (LA.accelerate2Massiv3DMatrix i1)
+  = let massiv     = LA.sumRows @_ @Massiv.U (LA.accelerate2Massiv3DMatrix i1)
         massiv'    = LA.sumRowsReferenceImplementation @Massiv.U (LA.accelerate2Massiv3DMatrix i1)
         accelerate =  Naive.run (A.sum (use i1))
     in counterexample "sumRows and reference implementation do not agree" (massiv === massiv') .&&.
@@ -175,6 +184,32 @@ compareDistributional Proxy (SquareMatrix i1)
     conv :: e -> Scientific
     conv = fromFloatDigits
 
+compareLogDistributional2 :: forall e.
+                          ( Eq e
+                          , Show e
+                          , FromIntegral Int e
+                          , Prelude.RealFloat e
+                          , Massiv.Unbox e
+                          , A.Ord e
+                          , Ord e
+                          , Prelude.Fractional (Exp e)
+                          , Prelude.Fractional e
+                          , Prelude.Floating (Exp e)
+                          , Prelude.Floating e
+                          , Monoid e
+                          ) => Proxy e
+                          -> SquareMatrix Int
+                          -> Property
+compareLogDistributional2 Proxy (SquareMatrix i1)
+  = let massiv     = Massiv.computeAs Massiv.B $ LA.logDistributional2 @_ @e (LA.accelerate2MassivMatrix i1)
+        accelerate = Legacy.logDistributional2With Naive.run i1
+        expected   = map conv (A.toList accelerate)
+        actual     = map conv (mconcat (Massiv.toLists2 massiv))
+    in counterexample "size not equal" (Prelude.length expected === Prelude.length actual) .&&. expected === actual
+  where
+    conv :: e -> Scientific
+    conv = fromFloatDigits
+
 compareMatMaxMini :: SquareMatrix Int -> Property
 compareMatMaxMini (SquareMatrix i1)
   = let massiv     = LA.matMaxMini @Massiv.U (LA.accelerate2MassivMatrix i1)
@@ -192,3 +227,15 @@ 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
+
+compareMatrixEye :: Positive Int -> Property
+compareMatrixEye (getPositive -> n)
+  = let massiv     = Massiv.compute @Massiv.U $ LA.matrixEye @Int n
+        accelerate = Naive.run (Legacy.matrixEye n)
+    in accelerate === LA.massiv2AccelerateMatrix massiv
+
+compareDiagNull :: SquareMatrix Int -> Property
+compareDiagNull (SquareMatrix i1)
+  = let massiv     = Massiv.compute @Massiv.U $ LA.diagNull (A.dim i1) (LA.accelerate2MassivMatrix i1)
+        accelerate =  Naive.run (Legacy.diagNull (A.dim i1) (use i1))
+    in accelerate === LA.massiv2AccelerateMatrix massiv