src/Gargantext/Core/Viz/Graph/Distances/Matrice.hs

   1 {-|
   2 Module      : Gargantext.Graph.Distances.Matrix
   3 Description :
   4 Copyright   : (c) CNRS, 2017-Present
   5 License     : AGPL + CECILL v3
   6 Maintainer  : team@gargantext.org
   7 Stability   : experimental
   8 Portability : POSIX
   9
  10 This module aims at implementig distances of terms context by context is
  11 the same referential of corpus.
  12
  13 Implementation use Accelerate library which enables GPU and CPU computation:
  14
  15   * Manuel M. T. Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell, and Vinod Grover.
  16     [Accelerating Haskell Array Codes with Multicore GPUs][CKLM+11].
  17     In _DAMP '11: Declarative Aspects of Multicore Programming_, ACM, 2011.
  18
  19   * Trevor L. McDonell, Manuel M. T. Chakravarty, Vinod Grover, and Ryan R. Newton.
  20     [Type-safe Runtime Code Generation: Accelerate to LLVM][MCGN15].
  21     In _Haskell '15: The 8th ACM SIGPLAN Symposium on Haskell_, ACM, 2015.
  22
  23 -}
  24
  25 {-# LANGUAGE TypeFamilies        #-}
  26 {-# LANGUAGE TypeOperators       #-}
  27 {-# LANGUAGE ScopedTypeVariables #-}
  28 {-# LANGUAGE ViewPatterns        #-}
  29
  30 module Gargantext.Core.Viz.Graph.Distances.Matrice
  31   where
  32
  33 import qualified Data.Foldable as P (foldl1)
  34 import Debug.Trace (trace)
  35 import Data.Array.Accelerate
  36 import Data.Array.Accelerate.Interpreter (run)
  37 import qualified Gargantext.Prelude as P
  38
  39
  40 -----------------------------------------------------------------------
  41 -- | Define a vector
  42 --
  43 -- >>> vector 3
  44 -- Vector (Z :. 3) [0,1,2]
  45 vector :: Elt c => Int -> [c] -> (Array (Z :. Int) c)
  46 vector n l = fromList (Z :. n) l
  47
  48 -- | Define a matrix
  49 --
  50 -- >>> matrix 3 ([1..] :: [Double])
  51 -- Matrix (Z :. 3 :. 3)
  52 --   [ 1.0, 2.0, 3.0,
  53 --     4.0, 5.0, 6.0,
  54 --     7.0, 8.0, 9.0]
  55 matrix :: Elt c => Int -> [c] -> Matrix c
  56 matrix n l = fromList (Z :. n :. n) l
  57
  58 -- | Two ways to get the rank (as documentation)
  59 --
  60 -- >>> rank (matrix 3 ([1..] :: [Int]))
  61 -- 2
  62 rank :: (Matrix a) -> Int
  63 rank m = arrayRank $ arrayShape m
  64
  65 -----------------------------------------------------------------------
  66 -- | Dimension of a square Matrix
  67 -- How to force use with SquareMatrix ?
  68 type Dim = Int
  69
  70 -- | Get Dimension of a square Matrix
  71 --
  72 -- >>> dim (matrix 3 ([1..] :: [Int]))
  73 -- 3
  74 dim :: Matrix a -> Dim
  75 dim m = n
  76   where
  77     Z :. _ :. n = arrayShape m
  78     -- indexTail (arrayShape m)
  79
  80 -----------------------------------------------------------------------
  81 -- TODO move to Utils
  82 runExp :: Elt e => Exp e -> e
  83 runExp e = indexArray (run (unit e)) Z
  84 -----------------------------------------------------------------------
  85
  86 -- | Sum of a Matrix by Column
  87 --
  88 -- >>> run $ matSumCol 3 (use $ matrix 3 [1..])
  89 -- Matrix (Z :. 3 :. 3)
  90 --   [ 12.0, 15.0, 18.0,
  91 --     12.0, 15.0, 18.0,
  92 --     12.0, 15.0, 18.0]
  93 matSumCol :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
  94 matSumCol r mat = replicate (constant (Z :. (r :: Int) :. All)) $ sum $ transpose mat
  95
  96 matSumCol' :: Matrix Double -> Matrix Double
  97 matSumCol' m = run $ matSumCol n m'
  98   where
  99     n  = dim m
 100     m' = use m
 101
 102
 103 -- | Proba computes de probability matrix: all cells divided by thee sum of its column
 104 -- if you need get the probability on the lines, just transpose it
 105 --
 106 -- >>> run $ matProba 3 (use $ matrix 3 [1..])
 107 -- Matrix (Z :. 3 :. 3)
 108 --   [ 8.333333333333333e-2, 0.13333333333333333, 0.16666666666666666,
 109 --       0.3333333333333333,  0.3333333333333333,  0.3333333333333333,
 110 --       0.5833333333333334,  0.5333333333333333,                 0.5]
 111 matProba :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
 112 matProba r mat = zipWith (/) mat (matSumCol r mat)
 113
 114 -- | Diagonal of the matrix
 115 --
 116 -- >>> run $ diag (use $ matrix 3 ([1..] :: [Int]))
 117 -- Vector (Z :. 3) [1,5,9]
 118 diag :: Elt e => Acc (Matrix e) -> Acc (Vector e)
 119 diag m = backpermute (indexTail (shape m))
 120                      (lift1 (\(Z :. x) -> (Z :. x :. (x :: Exp Int))))
 121                      m
 122
 123 -- | Divide by the Diagonal of the matrix
 124 --
 125 -- >>> run $ divByDiag 3 (use $ matrix 3 ([1..] :: [Double]))
 126 -- Matrix (Z :. 3 :. 3)
 127 --   [ 1.0, 0.4, 0.3333333333333333,
 128 --     4.0, 1.0, 0.6666666666666666,
 129 --     7.0, 1.6,                1.0]
 130 divByDiag :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
 131 divByDiag d mat = zipWith (/) mat (replicate (constant (Z :. (d :: Int) :. All)) $ diag mat)
 132
 133 -----------------------------------------------------------------------
 134 -- | Filters the matrix with the minimum of maximums
 135 --
 136 -- >>> run $ matMiniMax $ use $ matrix 3 [1..]
 137 -- Matrix (Z :. 3 :. 3)
 138 --   [ 0.0, 4.0, 7.0,
 139 --     0.0, 5.0, 8.0,
 140 --     0.0, 6.0, 9.0]
 141 matMiniMax :: (Elt a, Ord a, P.Num a) => Acc (Matrix a) -> Acc (Matrix a)
 142 matMiniMax m = filterWith' miniMax' (constant 0) m
 143   where
 144     miniMax' = the $ minimum $ maximum m
 145
 146
 147 -- | Filters the matrix with a constant
 148 --
 149 -- >>> run $ matFilter 5 $ use $ matrix 3 [1..]
 150 -- Matrix (Z :. 3 :. 3)
 151 --   [ 0.0, 0.0, 7.0,
 152 --     0.0, 0.0, 8.0,
 153 --     0.0, 6.0, 9.0]
 154 filter' :: Double -> Acc (Matrix Double) -> Acc (Matrix Double)
 155 filter' t m = filterWith t 0 m
 156
 157 filterWith :: Double -> Double -> Acc (Matrix Double) -> Acc (Matrix Double)
 158 filterWith t v m = map (\x -> ifThenElse (x > (constant t)) x (constant v)) (transpose m)
 159
 160 filterWith' :: (Elt a, Ord a) => Exp a -> Exp a -> Acc (Matrix a) -> Acc (Matrix a)
 161 filterWith' t v m = map (\x -> ifThenElse (x > t) x v) m
 162
 163
 164
 165
 166 -----------------------------------------------------------------------
 167 -- * Metrics of proximity
 168 -----------------------------------------------------------------------
 169 -- ** Conditional distance
 170
 171 -- *** Conditional distance (basic)
 172
 173 -- | Conditional distance (basic version)
 174 --
 175 -- 2 main metrics are actually implemented in order to compute the
 176 -- proximity of two terms: conditional and distributional
 177 --
 178 -- Conditional metric is an absolute metric which reflects
 179 -- interactions of 2 terms in the corpus.
 180 measureConditional :: Matrix Int -> Matrix Double
 181 --measureConditional m = run (matMiniMax $ matProba (dim m) $ map fromIntegral $ use m)
 182 measureConditional m = run $ matProba (dim m)
 183                            $ map fromIntegral
 184                            $ use m
 185
 186
 187 -- *** Conditional distance (advanced)
 188
 189 -- | Conditional distance (advanced version)
 190 --
 191 -- The conditional metric P(i|j) of 2 terms @i@ and @j@, also called
 192 -- "confidence" , is the maximum probability between @i@ and @j@ to see
 193 -- @i@ in the same context of @j@ knowing @j@.
 194 --
 195 -- If N(i) (resp. N(j)) is the number of occurrences of @i@ (resp. @j@)
 196 -- in the corpus and _[n_{ij}\] the number of its occurrences we get:
 197 --
 198 -- \[P_c=max(\frac{n_i}{n_{ij}},\frac{n_j}{n_{ij}} )\]
 199 conditional' :: Matrix Int -> (Matrix InclusionExclusion, Matrix SpecificityGenericity)
 200 conditional' m = ( run $ ie $ map fromIntegral $ use m
 201                  , run $ sg $ map fromIntegral $ use m
 202                  )
 203   where
 204     ie :: Acc (Matrix Double) -> Acc (Matrix Double)
 205     ie mat = map (\x -> x / (2*n-1)) $ zipWith (+) (xs mat) (ys mat)
 206     sg :: Acc (Matrix Double) -> Acc (Matrix Double)
 207     sg mat = map (\x -> x / (2*n-1)) $ zipWith (-) (xs mat) (ys mat)
 208
 209     n :: Exp Double
 210     n = P.fromIntegral r
 211
 212     r :: Dim
 213     r = dim m
 214
 215     xs :: Acc (Matrix Double) -> Acc (Matrix Double)
 216     xs mat = zipWith (-) (matSumCol r $ matProba r mat) (matProba r mat)
 217     ys :: Acc (Matrix Double) -> Acc (Matrix Double)
 218     ys mat = zipWith (-) (matSumCol r $ transpose $ matProba r mat) (matProba r mat)
 219
 220 -----------------------------------------------------------------------
 221 -- ** Distributional Distance
 222
 223 -- | Distributional Distance metric
 224 --
 225 -- Distributional metric is a relative metric which depends on the
 226 -- selected list, it represents structural equivalence of mutual information.
 227 --
 228 -- The distributional metric P(c) of @i@ and @j@ terms is: \[
 229 -- S_{MI} = \frac {\sum_{k \neq i,j ; MI_{ik} >0}^{} \min(MI_{ik},
 230 -- MI_{jk})}{\sum_{k \neq i,j ; MI_{ik}>0}^{}} \]
 231 --
 232 -- Mutual information
 233 -- \[S_{MI}({i},{j}) = \log(\frac{C{ij}}{E{ij}})\]
 234 --
 235 -- Number of cooccurrences of @i@ and @j@ in the same context of text
 236 --                        \[C{ij}\]
 237 --
 238 -- The expected value of the cooccurrences @i@ and @j@ (given a map list of size @n@)
 239 --            \[E_{ij}^{m} = \frac {S_{i} S_{j}} {N_{m}}\]
 240 --
 241 -- Total cooccurrences of term @i@ given a map list of size @m@
 242 --            \[S_{i} = \sum_{j, j \neq i}^{m} S_{ij}\]
 243 --
 244 -- Total cooccurrences of terms given a map list of size @m@
 245 --            \[N_{m} = \sum_{i,i \neq i}^{m} \sum_{j, j \neq j}^{m} S_{ij}\]
 246 --
 247 distributional :: Matrix Int -> Matrix Double
 248 distributional m = run -- $ matMiniMax
 249                        $ diagNull n
 250                        $ rIJ n
 251                        $ filterWith 0 100
 252                        $ filter' 0
 253                        $ s_mi
 254                        $ map fromIntegral
 255                           {- from Int to Double -}
 256                        $ use m
 257                           {- push matrix in Accelerate type -}
 258   where
 259
 260     _ri :: Acc (Matrix Double) -> Acc (Matrix Double)
 261     _ri mat = mat1 -- zipWith (/) mat1 mat2
 262       where
 263         mat1 = matSumCol n $ zipWith min (_myMin mat) (_myMin $ filterWith 0 100 $ diagNull n $ transpose mat)
 264         _mat2 = total mat
 265
 266     _myMin :: Acc (Matrix Double) -> Acc (Matrix Double)
 267     _myMin = replicate (constant (Z :. n :. All)) . minimum
 268
 269
 270     -- TODO fix NaN
 271     -- Quali TEST: OK
 272     s_mi :: Acc (Matrix Double) -> Acc (Matrix Double)
 273     s_mi m' = zipWith (\x y -> log (x / y)) (diagNull n m')
 274             $ zipWith (/) (crossProduct n m') (total m')
 275             -- crossProduct n m'
 276
 277
 278     total :: Acc (Matrix Double) -> Acc (Matrix Double)
 279     total = replicate (constant (Z :. n :. n)) . sum . sum
 280
 281     n :: Dim
 282     n = dim m
 283
 284 -- run $ (identityMatrix (DAA.constant (10::Int)) :: DAA.Acc (DAA.Matrix Int)) Matrix (Z :. 10 :. 10)
 285 identityMatrix :: Num a => Exp Int -> Acc (Matrix a)
 286 identityMatrix n =
 287         let zeros = fill (index2 n n) 0
 288             ones  = fill (index1 n)   1
 289         in
 290         permute const zeros (\(unindex1 -> i) -> index2 i i) ones
 291
 292
 293 eyeMatrix :: Num a => Dim -> Acc (Matrix a)
 294 eyeMatrix n' =
 295         let ones   = fill (index2 n n) 1
 296             zeros  = fill (index1 n)   0
 297             n = constant n'
 298         in
 299         permute const ones (\(unindex1 -> i) -> index2 i i) zeros
 300
 301 -- | TODO use Lenses
 302 data Direction = MatCol (Exp Int) | MatRow (Exp Int) | Diag
 303
 304 nullOf :: Num a => Dim -> Direction -> Acc (Matrix a)
 305 nullOf n' dir =
 306         let ones   = fill (index2 n n) 1
 307             zeros  = fill (index2 n n) 0
 308             n = constant n'
 309         in
 310         permute const ones ( lift1 ( \(Z :. (i :: Exp Int) :. (_j:: Exp Int))
 311                                                 -> case dir of
 312                                                      MatCol m -> (Z :. i :. m)
 313                                                      MatRow m -> (Z :. m :. i)
 314                                                      Diag     -> (Z :. i :. i)
 315                                    )
 316                            )
 317                            zeros
 318
 319 nullOfWithDiag :: Num a => Dim -> Direction -> Acc (Matrix a)
 320 nullOfWithDiag n dir = zipWith (*) (nullOf n dir) (nullOf n Diag)
 321
 322
 323 rIJ' :: Matrix Int -> Matrix Double
 324 rIJ' m = run $ sumRowMin (dim m) m'
 325   where
 326     m' = (map fromIntegral $ use m)
 327
 328 rIJ :: (Elt a, Ord a, P.Fractional (Exp a), P.Num a)
 329     => Dim -> Acc (Matrix a) -> Acc (Matrix a)
 330 rIJ n m = matMiniMax $ divide a b
 331   where
 332     a = sumRowMin n m
 333     b = sumColMin n m
 334
 335 divide :: (Elt a, Ord a, P.Fractional (Exp a), P.Num a)
 336     => Acc (Matrix a) -> Acc (Matrix a) -> Acc (Matrix a)
 337 divide = zipWith divide'
 338   where
 339     divide' a b = ifThenElse (b > (constant 0))
 340                              (a / b)
 341                              (constant 0)
 342
 343 -- | Nominator
 344 sumRowMin :: (Num a, Ord a) => Dim -> Acc (Matrix a) -> Acc (Matrix a)
 345 sumRowMin n m = {-trace (P.show $ run m') $-} m'
 346   where
 347     m' = reshape (shape m) vs
 348     vs = P.foldl1 (++)
 349        $ P.map (\z -> sumRowMin1 n (constant z) m) [0..n-1]
 350
 351 sumRowMin1 :: (Num a, Ord a) => Dim -> Exp Int -> Acc (Matrix a) -> Acc (Vector a)
 352 sumRowMin1 n x m = trace (P.show (run m,run $ transpose m)) $ m''
 353   where
 354     m'' = sum $ zipWith min (transpose m) m
 355     _m'  = zipWith (*) (zipWith (*) (nullOf n (MatCol x)) $ nullOfWithDiag n (MatRow x)) m
 356
 357
 358 -- | Denominator
 359 sumColMin :: (Num a, Ord a) => Dim -> Acc (Matrix a) -> Acc (Matrix a)
 360 sumColMin n m = reshape (shape m) vs
 361   where
 362     vs = P.foldl1 (++)
 363        $ P.map (\z -> sumColMin1 n (constant z) m) [0..n-1]
 364
 365
 366 sumColMin1 :: (Num a) => Dim -> Exp Int -> Acc (Matrix a) -> Acc (Matrix a)
 367 sumColMin1 n x m = zipWith (*) (nullOfWithDiag n (MatCol x)) m
 368
 369
 370
 371 {- | WIP fun with indexes
 372 selfMatrix :: Num a => Dim -> Acc (Matrix a)
 373 selfMatrix n' =
 374         let zeros = fill (index2 n n) 0
 375             ones  = fill (index2 n n) 1
 376             n = constant n'
 377         in
 378         permute const ones ( lift1 ( \(Z :. (i :: Exp Int) :. (_j:: Exp Int))
 379                                                 -> -- ifThenElse (i /= j)
 380                                                      --         (Z :. i :. j)
 381                                                               (Z :. i :. i)
 382                                     )) zeros
 383
 384 selfMatrix' :: (Elt a, P.Num (Exp a)) => Array DIM2 a -> Matrix a
 385 selfMatrix' m' = run $ selfMatrix n
 386   where
 387     n = dim m'
 388     m = use m'
 389 -}
 390 -------------------------------------------------
 391 diagNull :: Num a => Dim -> Acc (Matrix a) -> Acc (Matrix a)
 392 diagNull n m = zipWith (*) m eye
 393   where
 394     eye = eyeMatrix n
 395
 396 -------------------------------------------------
 397 crossProduct :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
 398 crossProduct n m = {-trace (P.show (run m',run m'')) $-} zipWith (*) m' m''
 399   where
 400     m'  = cross n m
 401     m'' = transpose $ cross n m
 402
 403
 404 crossT :: Matrix Double -> Matrix Double
 405 crossT  = run . transpose . use
 406
 407 crossProduct' :: Matrix Double -> Matrix Double
 408 crossProduct' m = run $ crossProduct n m'
 409   where
 410     n  = dim m
 411     m' = use m
 412
 413 runWith :: (Arrays c, Elt a1)
 414         => (Dim -> Acc (Matrix a1) -> a2 -> Acc c)
 415         -> Matrix a1
 416         -> a2
 417         -> c
 418 runWith f m = run . f (dim m) (use m)
 419
 420 -- | cross
 421 cross :: Dim -> Acc (Matrix Double) -> Acc (Matrix Double)
 422 cross n mat = diagNull n (matSumCol n $ diagNull n mat)
 423
 424 cross' :: Matrix Double -> Matrix Double
 425 cross' mat = run $ cross n mat'
 426   where
 427     mat' = use mat
 428     n = dim mat
 429
 430
 431 -----------------------------------------------------------------------
 432 -----------------------------------------------------------------------
 433 -- * Specificity and Genericity
 434
 435 {- | Metric Specificity and genericity: select terms
 436
 437 - let N termes and occurrences of i \[N{i}\]
 438
 439 - Cooccurrences of i and j \[N{ij}\]
 440 - Probability to get i given j : \[P(i|j)=N{ij}/N{j}\]
 441
 442 - Genericity of i  \[Gen(i)  = \frac{\sum_{j \neq i,j} P(i|j)}{N-1}\]
 443 - Specificity of j \[Spec(i) = \frac{\sum_{j \neq i,j} P(j|i)}{N-1}\]
 444
 445 - \[Inclusion (i) = Gen(i) + Spec(i)\)
 446 - \[GenericityScore = Gen(i)- Spec(i)\]
 447
 448 - References: Science mapping with asymmetrical paradigmatic proximity
 449 Jean-Philippe Cointet (CREA, TSV), David Chavalarias (CREA) (Submitted
 450 on 15 Mar 2008), Networks and Heterogeneous Media 3, 2 (2008) 267 - 276,
 451 arXiv:0803.2315 [cs.OH]
 452 -}
 453 type InclusionExclusion    = Double
 454 type SpecificityGenericity = Double
 455
 456 data SquareMatrix      = SymetricMatrix | NonSymetricMatrix
 457 type SymetricMatrix    = Matrix
 458 type NonSymetricMatrix = Matrix
 459
 460
 461 incExcSpeGen :: Matrix Int -> (Vector InclusionExclusion, Vector SpecificityGenericity)
 462 incExcSpeGen m = (run' inclusionExclusion m, run' specificityGenericity m)
 463   where
 464     run' fun mat = run $ fun $ map fromIntegral $ use mat
 465
 466     -- | Inclusion (i) = Gen(i)+Spec(i)
 467     inclusionExclusion :: Acc (Matrix Double) -> Acc (Vector Double)
 468     inclusionExclusion mat = zipWith (+) (pV mat) (pV mat)
 469
 470     -- | Genericity score = Gen(i)- Spec(i)
 471     specificityGenericity :: Acc (Matrix Double) -> Acc (Vector Double)
 472     specificityGenericity mat = zipWith (+) (pH mat) (pH mat)
 473
 474     -- | Gen(i)  : 1/(N-1)*Sum(j!=i, P(i|j)) : Genericity  of i
 475     pV :: Acc (Matrix Double) -> Acc (Vector Double)
 476     pV mat = map (\x -> (x-1)/(cardN-1)) $ sum $ p_ij mat
 477
 478     -- | Spec(i) : 1/(N-1)*Sum(j!=i, P(j|i)) : Specificity of j
 479     pH :: Acc (Matrix Double) -> Acc (Vector Double)
 480     pH mat = map (\x -> (x-1)/(cardN-1)) $ sum $ p_ji mat
 481
 482     cardN :: Exp Double
 483     cardN = constant (P.fromIntegral (dim m) :: Double)
 484
 485
 486 -- | P(i|j) = Nij /N(jj) Probability to get i given j
 487 --p_ij :: (Elt e, P.Fractional (Exp e)) => Acc (SymetricMatrix e) -> Acc (Matrix e)
 488 p_ij :: (Elt e, P.Fractional (Exp e)) => Acc (Matrix e) -> Acc (Matrix e)
 489 p_ij m = zipWith (/) m (n_jj m)
 490   where
 491     n_jj :: Elt e => Acc (SymetricMatrix e) -> Acc (Matrix e)
 492     n_jj myMat' = backpermute (shape m)
 493                          (lift1 ( \(Z :. (_ :: Exp Int) :. (j:: Exp Int))
 494                                    -> (Z :. j :. j)
 495                                 )
 496                          ) myMat'
 497
 498 -- | P(j|i) = Nij /N(ii) Probability to get i given j
 499 -- to test
 500 p_ji :: (Elt e, P.Fractional (Exp e)) => Acc (Array DIM2 e) -> Acc (Array DIM2 e)
 501 p_ji = transpose . p_ij
 502
 503
 504 -- | Step to ckeck the result in visual/qualitative tests
 505 incExcSpeGen_proba :: Matrix Int -> Matrix Double
 506 incExcSpeGen_proba m = run' pro m
 507   where
 508     run' fun mat = run $ fun $ map fromIntegral $ use mat
 509
 510     pro mat = p_ji mat
 511
 512 {-
 513 -- | Hypothesis to test maybe later (or not)
 514 -- TODO ask accelerate for instances to ease such writtings:
 515 p_ :: (Elt e, P.Fractional (Exp e)) => Acc (Array DIM2 e) -> Acc (Array DIM2 e)
 516 p_ m = zipWith (/) m (n_ m)
 517   where
 518     n_ :: Elt e => Acc (SymetricMatrix e) -> Acc (Matrix e)
 519     n_ m = backpermute (shape m)
 520                          (lift1 ( \(Z :. (i :: Exp Int) :. (j:: Exp Int))
 521                                    -> (ifThenElse (i < j) (lift (Z :. j :. j)) (lift (Z :. i :. i)) :: Exp DIM2)
 522                                 )
 523                          ) m
 524 -}
 525
 526 -- * For Tests (to be removed)
 527 -- | Test perfermance with this matrix
 528 -- TODO : add this in a benchmark folder
 529 distriTest :: Int -> Matrix Double
 530 distriTest n = distributional (theMatrix n)
 531
 532 theMatrix :: Int -> Matrix Int
 533 theMatrix n = matrix n (dataMatrix n)
 534   where
 535     dataMatrix :: Int -> [Int]
 536     dataMatrix x | (P.==) x 2 = [ 1, 1
 537                                 , 1, 2
 538                                 ]
 539
 540                  | (P.==) x 3 =  [ 1, 1, 2
 541                                  , 1, 2, 3
 542                                  , 2, 3, 4
 543                                  ]
 544                  | (P.==) x 4 =  [ 1, 1, 2, 3
 545                                  , 1, 2, 3, 4
 546                                  , 2, 3, 4, 5
 547                                  , 3, 4, 5, 6
 548                                  ]
 549                  | P.otherwise = P.undefined
 550
 551 {-
 552 theResult :: Int -> Matrix Double
 553 theResult n | (P.==) n 2 = let r = 1.6094379124341003 in [ 0, r, r, 0]
 554           | P.otherwise = [ 1, 1 ]
 555 -}
 556
 557
 558 colMatrix :: Elt e
 559           => Int -> [e] -> Acc (Array ((Z :. Int) :. Int) e)
 560 colMatrix n ns = replicate (constant (Z :. (n :: Int) :. All)) v
 561   where
 562     v = use $ vector (P.length ns) ns
 563
 564 -----------------------------------------------------------------------
 565