WIP: implement distributional in terms of massiv

gargantext.cabal

@@ -589,6 +589,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
@@ -749,6 +750,7 @@ common testDependencies
     , raw-strings-qq
     , resource-pool >= 0.4.0.0 && < 0.5
     , safe-exceptions >= 0.1.7.4 && < 0.2
+    , scientific < 0.4
     , servant-auth-client
     , servant-client >= 0.20 && < 0.21
     , servant-client-core >= 0.20 && < 0.21

src/Gargantext/Core/LinearAlgebra.hs

@@ -20,8 +20,10 @@ module Gargantext.Core.LinearAlgebra (
   , createIndices
 
   -- * Operations on matrixes
+  , (.*)
   , diag
   , termDivNan
+  , distributional
   ) where
 
 import Data.Bimap (Bimap)
@@ -32,7 +34,7 @@ import Data.Set qualified as S
 import Data.Set (Set)
 import Prelude
 import Data.Massiv.Array qualified as A
-import Data.Massiv.Array (D, Matrix, Vector)
+import Data.Massiv.Array (D, Matrix, Vector, Array, Ix3)
 
 newtype Index = Index { _Index :: Int }
   deriving newtype (Eq, Show, Ord, Num, Enum)
@@ -49,8 +51,8 @@ createIndices = set2indices . map2set
 -- | Computes the diagnonal matrix of the input one.
 diag :: Num e => Matrix D e -> Vector D e
 diag matrix =
-  let (A.Sz2 _rows cols) = A.size matrix
-      newSize = A.Sz1 cols
+  let (A.Sz2 rows _cols) = A.size matrix
+      newSize = A.Sz1 rows
   in A.backpermute' newSize (\j -> j A.:. j) matrix
 
 
@@ -91,37 +93,79 @@ diag matrix =
 --
 -- /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 D Int -> Matrix D Double
+distributional :: forall e. (Ord e, Fractional e, Num e)
+               => Matrix D Int
+               -> Matrix D e
 distributional m' = result
  where
-    m = map A.fromIntegral $ use m'
+    m :: Matrix D e
+    m = A.map fromIntegral m'
     n = dim m'
 
-    diag_m = A.diag m
+    diag_m :: Vector D e
+    diag_m = diag m
 
-    d_1 = replicate (constant (Z :. n :. All)) diag_m
-    d_2 = replicate (constant (Z :. All :. n)) diag_m
+    -- replicate (constant (Z :. n :. All)) diag_m
+    d_1 :: Matrix D e
+    d_1 = let (A.Sz1 r) = A.size diag_m
+          in A.backpermute' (A.Sz2 n r) (\(_ A.:. i) -> i) diag_m
 
-    mi    = (A..*) (termDivNan m d_1) (termDivNan m d_2)
+    -- replicate (constant (Z :. All :. n)) diag_m
+    d_2 :: Matrix D e
+    d_2 = let (A.Sz1 r) = A.size diag_m
+          in A.backpermute' (A.Sz2 r n) (\(i A.:. _) -> i) diag_m
+
+    mi :: Matrix D e
+    mi = (.*) (termDivNan m d_1) (termDivNan m d_2)
 
     -- The matrix permutations is taken care of below by directly replicating
     -- the matrix mi, making the matrix w unneccessary and saving one step.
-    w_1 = replicate (constant (Z :. All :. n :. All)) mi
-    w_2 = replicate (constant (Z :. n :. All :. All)) mi
-    w' = zipWith min w_1 w_2
+    -- replicate (constant (Z :. All :. n :. All)) mi
+    w_1 :: Array D Ix3 e
+    w_1 = let (A.Sz2 r c) = A.size mi
+          in A.backpermute' (A.Sz3 r n c) (\(x A.:> _y A.:. z) -> x A.:. z) mi
+
+    -- replicate (constant (Z :. n :. All :. All)) mi
+    w_2 :: Array D Ix3 e
+    w_2 = let (A.Sz2 r c) = A.size mi
+          in A.backpermute' (A.Sz3 n r c) (\(_x A.:> y A.:. z) -> y A.:. z) mi
+
+    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). 
-    ii = generate (constant (Z :. n :. n :. n))
-        (lift1 (\( i A.:. j A.:. k) -> cond ((&&) ((/=) k i) ((/=) k j)) 1 0))
+    -- 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 D Ix3 e
+    ii = A.makeArrayR A.D A.Par (A.Sz3 n n n) $ \(i A.:> j A.:. k) -> if k /= i && k /= j then 1 else 0
+
+    z_1 :: Matrix D e
+    z_1 = sumRows ((.*) w' ii)
 
-    z_1 = sum ((A..*) w' ii)
-    z_2 = sum ((A..*) w_1 ii)
+    z_2 :: Matrix D e
+    z_2 = sumRows ((.*) w_1 ii)
 
     result = termDivNan z_1 z_2
--}
 
 -- | Term by term division where divisions by 0 produce 0 rather than NaN.
 termDivNan :: (Eq a, Fractional a) => Matrix D a -> Matrix D a -> Matrix D a
 termDivNan = A.zipWith (\i j -> if j == 0 then 0 else i / j)
+
+sumRows :: Num e => Array D A.Ix3 e -> Array D A.Ix2 e
+sumRows 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.Index ix, Num a)
+     => Array D ix a
+     -> Array D ix a
+     -> Array D ix a
+(.*) = A.zipWith (*)
+
+-- | 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

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

@@ -143,7 +143,10 @@ import qualified Prelude
 distributional :: Matrix Int -> Matrix Double
 distributional = distributionalWith LLVM.run
 
-distributionalWith :: (forall a. Arrays a => Acc a -> a) -> Matrix Int -> Matrix Double
+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'

src/Gargantext/Orphans/Accelerate.hs

@@ -6,6 +6,7 @@ module Gargantext.Orphans.Accelerate where
 
 import Prelude
 import Test.QuickCheck
+import Data.Scientific ()
 import Data.Array.Accelerate (DIM2, Z (..), (:.) (..), Array, Elt, fromList, arrayShape)
 import Data.Array.Accelerate qualified as A
 import qualified Data.List.Split as Split
@@ -23,12 +24,9 @@ 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) -> [ resizeArray arr (max 1 (x - 1)) y  ]
-      (1,_) -> [ resizeArray arr x (max 1 ( y - 1)) ]
-      _     -> [ resizeArray arr (max 1 (x - 1)) (max 1 y)
-               , resizeArray arr (max 1 x) (max 1 (y - 1))
-               , resizeArray arr (max 1 (x - 1)) (max 1 (y - 1))
-               ]
+      (_,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 =

test/Test/Core/LinearAlgebra.hs

@@ -13,18 +13,21 @@ 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 Gargantext.Core.LinearAlgebra qualified as LA
 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 Gargantext.Prelude (headMay)
 import Prelude hiding ((^))
-import Test.Tasty
-import Test.Tasty.QuickCheck
-import Data.Massiv.Array qualified as Massiv
 import qualified Data.Array.Accelerate as A
 import qualified Data.List.Split.Internals as Split
-import Gargantext.Prelude (headMay)
+import Test.Tasty
+import Test.Tasty.QuickCheck
+import Data.Proxy
+import Data.Scientific
+
 
 --
 -- Utility types and functions
@@ -33,11 +36,11 @@ import Gargantext.Prelude (headMay)
 newtype SquareMatrix a = SquareMatrix { _SquareMatrix :: Matrix a }
   deriving newtype (Show, Eq)
 
-instance (Elt a, Show a, Arbitrary a) => Arbitrary (SquareMatrix a) where
+instance (Elt a, Show a, Prelude.Num a, Ord a, Arbitrary a) => Arbitrary (SquareMatrix a) where
   arbitrary = do
-    x <- choose (0,30)
+    x <- choose (1,30)
     let sh = Z :. x :. x
-    SquareMatrix . fromList sh <$> vectorOf (x*x) arbitrary
+    SquareMatrix . A.fromList sh <$> vectorOf (x*x) arbitrary
   shrink = map (SquareMatrix) . sliceArray . _SquareMatrix
 
 compareImplementations :: (Arbitrary a, Eq b, Show b)
@@ -63,8 +66,36 @@ mapCreateIndices (_m1, m2) = Bimap.fromList $ map (first LA.Index) $ M.toList m2
 
 type TermDivNanShape = Z :. Int :. Int
 
-twoByTwo :: Matrix Int
-twoByTwo = fromList (Z :. 2 :. 2) (Prelude.replicate 4 0)
+twoByTwo :: SquareMatrix Int
+twoByTwo = SquareMatrix $ fromList (Z :. 2 :. 2) (Prelude.replicate 4 0)
+
+testMatrix_01 :: SquareMatrix Int
+testMatrix_01 = SquareMatrix $ fromList (Z :. 14 :. 14) $
+  [  30,  36, -36, -16,   0,   7,  34,  -7,   5,  -4,   0,  21,   6, -35,
+      0, -31,  20, -15, -22,  -7, -22, -37, -29, -29,  23, -31, -29, -23,
+    -24, -29,  19,  -6,  16,   7,  15, -27, -27, -30,  -9, -33,  18, -23,
+      7, -36,  12,  26, -17,  -3,  -2, -15,  -4,  26,  24,   9,  -4,   4,
+     32,  28,  -2, -10,  34,  -3,  20,  -9, -22,  20, -26,  34,  18, -21,
+      7, -12,  12,  -2,  36,  10,  34, -37,  13,  -9, -28,  34,  33, -18,
+     -4, -32,  -1,  29,  29, -28,  24,  28,  35,  19,   8, -18,  25, -35,
+    -14,  -4, -24,  -1,   7,  34, -37, -28, -12, -32,  -5, -23,  27,  33,
+    -36, -28,  21, -29,  -2, -26,  -4, -31, -26, -21,  33, -11, -33,  20,
+     25,  14,   5,  -7,   5,  24,  37,   1,  -3,  23,  25, -16,  17,   5,
+    -35,  36,  -2,  -2,   1, -14,  34, -30, -10,  12,  25,  21,   0,  34,
+     17,  -1,  20, -19,  15,  20,  -5, -30, -35, -13,   5,  17, -10, -19,
+    -34, -11, -18,  26, -29, -28,   0,   3,  23,  -6,  36,   4,  16,  28,
+     13, -37, -16,   2,   7, -13,  21, -10, -33, -33, -26, -19,  -1,  29]
+
+testMatrix_02 :: SquareMatrix Int
+testMatrix_02 = SquareMatrix $ fromList (Z :. 7 :. 7) $
+  [  30,  36, -36, -16,   0,   7,  34,
+      0, -31,  20, -15, -22,  -7, -22,
+    -24, -29,  19,  -6,  16,   7,  15,
+      7, -36,  12,  26, -17,  -3,  -2,
+     32,  28,  -2, -10,  34,  -3,  20,
+      7, -12,  12,  -2,  36,  10,  34,
+     13, -37, -16,   2,   7, -13,  21]
+
 
 accelerate2MassivMatrix :: (Massiv.Unbox a, Elt a) => Matrix a -> Massiv.Matrix Massiv.D a
 accelerate2MassivMatrix m =
@@ -96,10 +127,10 @@ tests = testGroup "LinearAlgebra" [
     , testProperty "termDivNan"     compareTermDivNan
     , testProperty "diag"           compareDiag
     , testGroup "distributional" [
-        testProperty "2x2" (compareDistributional twoByTwo)
-      , testProperty "roundtrips" (compareImplementations' @(SquareMatrix Int) (Legacy.distributionalWith Naive.run . _SquareMatrix)
-                                                                               (Legacy.distributionalWith Naive.run . _SquareMatrix)
-                                                                               id)
+        testProperty "2x2" (compareDistributional (Proxy @Double) twoByTwo)
+      , testProperty "7x7" (compareDistributional (Proxy @Double) testMatrix_02)
+      , testProperty "14x14" (compareDistributional (Proxy @Double) testMatrix_01)
+      , testProperty "roundtrips" (compareDistributional (Proxy @Double))
       ]
     ]
 
@@ -121,8 +152,24 @@ compareDiag (SquareMatrix i1)
         accelerate =  Naive.run (Legacy.diag (use i1))
     in accelerate === massiv2AccelerateVector massiv
 
-compareDistributional :: Matrix Int
+compareDistributional :: forall e.
+                      ( Eq e
+                      , Show e
+                      , FromIntegral Int e
+                      , Prelude.RealFloat e
+                      , A.Ord e
+                      , Ord e
+                      , Prelude.Fractional (Exp e)
+                      , Prelude.Fractional e
+                      ) => Proxy e
+                      -> SquareMatrix Int
                       -> Property
-compareDistributional i1
-  = Legacy.distributionalWith Naive.run i1 === Legacy.distributionalWith Naive.run i1
-
+compareDistributional Proxy (SquareMatrix i1)
+  = let massiv     = Massiv.computeAs Massiv.B $ LA.distributional @e (accelerate2MassivMatrix i1)
+        accelerate = Legacy.distributionalWith 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