[FEAT] Implementing distributional function with matrix computation.

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

@@ -49,34 +49,71 @@ import Data.Array.Accelerate
 import Data.Array.Accelerate.Interpreter (run)
 import Gargantext.Core.Methods.Matrix.Accelerate.Utils
 import qualified Gargantext.Prelude as P
--- import Data.Array.Accelerate.LinearAlgebra (identity) TODO
------------------------------------------------------------------------
 
--- * Distributional Distance
+-- | `distributional m` returns the distributional distance between terms 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]
+--
 distributional :: Matrix Int -> Matrix Double
 distributional m' = run result
  where
     m = map fromIntegral $ use m'
     n = dim m'
 
-    d_m   = (.*) (matrixIdentity n) m
-
-    o_d_m = (#*#) (matrixOne n) d_m
-    d_m_o = transpose o_d_m
+    diag_m = diag m
 
-    mi    = (.*) ((./) m o_d_m) ((./) m d_m_o)
-    d_mi  = (.*) (matrixIdentity n) mi
+    d_1 = replicate (constant (Z :. n :. All)) diag_m
+    d_2 = replicate (constant (Z :. All :. n)) diag_m
 
-    w     = (.-) mi d_mi
+    mi    = (.*) ((./) m d_1) ((./) m d_2)
 
-    z     = (#*#) w (matrixOne n)
-    z'    = transpose z
+    -- w = (.-) mi d_mi
     
-    min_z_z' = zipWith min z z'
+    -- 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
 
-    result = (./) min_z_z' z
+    -- 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 (\(Z :. i :. j :. k) -> cond ((&&) ((/=) k i) ((/=) k j)) 1 0))
 
+    z_1 = sum ((.*) w' ii)
+    z_2 = sum ((.*) w_1 ii)
 
+    result = termDivNan z_1 z_2
 
 
 --

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

@@ -35,12 +35,9 @@ import Debug.Trace (trace)
 import Data.Array.Accelerate
 import Data.Array.Accelerate.Interpreter (run)
 import qualified Gargantext.Prelude as P
-import Data.Array.Accelerate.LinearAlgebra hiding (Matrix, transpose, Vector)
 
------------------------------------------------------------------------
--- | Main operators
--- Matrix Multiplication
-(#*#) :: ( Shape ix
+-- | Matrix cell by cell multiplication
+(.*) :: ( Shape ix
         , Slice ix
         , Elt a
         , P.Num (Exp a)
@@ -48,31 +45,32 @@ import Data.Array.Accelerate.LinearAlgebra hiding (Matrix, transpose, Vector)
      => Acc (Array ((ix :. Int) :. Int) a)
      -> Acc (Array ((ix :. Int) :. Int) a)
      -> Acc (Array ((ix :. Int) :. Int) a)
-(#*#) = multiplyMatrixMatrix
+(.*) = zipWith (*)
 
 
--- | Matrix cell by cell multiplication
-(.*) :: ( Shape ix
+(./) :: ( Shape ix
         , Slice ix
         , Elt 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 (/)
 
-(./) :: ( Shape ix
+-- | Term by term division where divisions by 0 produce 0 rather than NaN. 
+termDivNan :: ( 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 (/)
+termDivNan = zipWith (\i j -> cond ((==) j 0) 0 ((/) i j))
 
 (.-) :: ( Shape ix
         , Slice ix
@@ -399,11 +397,13 @@ theMatrixInt n = matrix n (dataMatrix n)
                                  , 4, 5, 3
                                  , 0, 3, 4
                                  ]
-                 | (P.==) x 4 =  [ 4, 4, 0, 0
-                                 , 4, 4, 0, 0
-                                 , 0, 0, 3, 3
-                                 , 0, 0, 3, 3
+                 | (P.==) x 4 =  [ 4, 1, 2, 1
+                                 , 1, 4, 0, 0
+                                 , 2, 0, 3, 3
+                                 , 1, 0, 3, 3
                                  ]
+
+
                  | P.otherwise = P.undefined
 
 {-