Skip to content

haskell-gargantext
haskell-gargantext
Commit 1909fd3a authored by Alfredo Di Napoli's avatar Alfredo Di Napoli Committed by Alfredo Di Napoli
Browse files
Options

WIP: try to use away intermediate matrixes in distributional

parent 756aea9e
Hide whitespace changes
Inline Side-by-side
Showing with 438 additions and 337 deletions
bench/Main.hs
View file @ 1909fd3a
......@@ -130,6 +130,7 @@ main = do
, bgroup "distributional" [
bench "Accelerate (Naive)" $ nf (Accelerate.distributionalWith @Double Naive.run) accInput
, bench "Accelerate (LLVM)" $ nf Accelerate.distributional accInput
, bench "Massiv (reference implementation)" $ nf (LA.distributionalReferenceImplementation @_ @Double) massivInput
, bench "Massiv " $ nf (LA.distributional @_ @Double) massivInput
]
, bgroup "logDistributional2" [
......
gargantext.cabal
View file @ 1909fd3a
......@@ -189,6 +189,8 @@ library
Gargantext.Core.Config.Utils
Gargantext.Core.Config.Worker
Gargantext.Core.LinearAlgebra
Gargantext.Core.LinearAlgebra.Distributional
Gargantext.Core.LinearAlgebra.Operations
Gargantext.Core.Mail
Gargantext.Core.Mail.Types
Gargantext.Core.Methods.Matrix.Accelerate.Utils
......
src/Gargantext/Core/LinearAlgebra.hs
View file @ 1909fd3a
......@@ -21,41 +21,20 @@ module Gargantext.Core.LinearAlgebra (
-- * Functions
, createIndices
-- * Convertion functions
, accelerate2MassivMatrix
, accelerate2Massiv3DMatrix
, massiv2AccelerateMatrix
, massiv2AccelerateVector
-- * Operations on matrixes
, (.*)
, (.-)
, diag
, termDivNan
, distributional
--, logDistributional
, sumRows
-- * Internals for testing
, sumRowsReferenceImplementation
, matMaxMini
, sumM_go
, sumMin_go
-- * Handy re-exports
, module Gargantext.Core.LinearAlgebra.Operations
, module Gargantext.Core.LinearAlgebra.Distributional
) where
import Data.Array.Accelerate qualified as Acc
import Data.Bimap (Bimap)
import Data.Bimap qualified as Bimap
import Data.List.Split qualified as Split
import Data.Map.Strict (Map)
import Data.Map.Strict qualified as M
import Data.Massiv.Array (D, Matrix, Vector, Array, Ix3, U)
import Data.Massiv.Array qualified as A
import Data.Set qualified as S
import Data.Set (Set)
import Prelude
import Protolude.Safe (headMay)
import Data.Monoid
import Gargantext.Core.LinearAlgebra.Operations
import Gargantext.Core.LinearAlgebra.Distributional
newtype Index = Index { _Index :: Int }
deriving newtype (Eq, Show, Ord, Num, Enum)
......@@ -68,313 +47,3 @@ createIndices = set2indices . map2set
set2indices :: Ord t => Set t -> Bimap Index t
set2indices s = foldr (uncurry Bimap.insert) Bimap.empty (zip [0..] $ S.toList s)
-- | Converts an accelerate matrix into a Massiv matrix.
accelerate2MassivMatrix :: (A.Unbox a, Acc.Elt a) => Acc.Matrix a -> Matrix A.U a
accelerate2MassivMatrix m =
let (Acc.Z Acc.:. _r Acc.:. c) = Acc.arrayShape m
in A.fromLists' @A.U A.Par $ Split.chunksOf c (Acc.toList m)
-- | Converts a massiv matrix into an accelerate matrix.
massiv2AccelerateMatrix :: (Acc.Elt a, A.Source r a) => Matrix r a -> Acc.Matrix a
massiv2AccelerateMatrix m =
let m' = A.toLists2 m
r = Prelude.length m'
c = maybe 0 Prelude.length (headMay m')
in Acc.fromList (Acc.Z Acc.:. r Acc.:. c) (mconcat m')
-- | Converts a massiv vector into an accelerate one.
massiv2AccelerateVector :: (A.Source r a, Acc.Elt a) => A.Vector r a -> Acc.Vector a
massiv2AccelerateVector m =
let m' = A.toList m
r = Prelude.length m'
in Acc.fromList (Acc.Z Acc.:. r) m'
accelerate2Massiv3DMatrix :: (A.Unbox e, Acc.Elt e, A.Manifest r e)
=> Acc.Array (Acc.Z Acc.:. Int Acc.:. Int Acc.:. Int) e
-> A.Array r A.Ix3 e
accelerate2Massiv3DMatrix m =
let (Acc.Z Acc.:. _r Acc.:. _c Acc.:. _z) = Acc.arrayShape m
in A.fromLists' A.Par $ map (Split.chunksOf $ _z) $ Split.chunksOf (_c*_z) (Acc.toList m)
-- | Computes the diagnonal matrix of the input one.
diag :: (A.Unbox e, A.Manifest r e, A.Source r e, Num e) => Matrix r e -> Vector A.U e
diag matrix =
let (A.Sz2 rows _cols) = A.size matrix
newSize = A.Sz1 rows
in A.makeArrayR A.U A.Seq newSize $ (\(A.Ix1 i) -> matrix A.! (A.Ix2 i i))
-- | `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]
--
-- /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 :: forall r e.
( A.Manifest r e
, A.Unbox e
, A.Source r Int
, A.Size r
, Ord e
, Fractional e
, Num e
)
=> Matrix r Int
-> Matrix r e
distributional m' = result
where
mD :: Matrix D e
mD = A.map fromIntegral m'
m :: Matrix A.U e
m = A.compute mD
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_2 :: Matrix A.D e
d_2 = A.backpermute' (A.Sz2 diag_m_size n) (\(i A.:. _) -> i) diag_m
a :: Matrix D e
a = termDivNanD mD d_1
b :: Matrix D e
b = termDivNanD mD d_2
miDelayed :: Matrix D e
miDelayed = a `mulD` b
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
-- 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' :: 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).
-- generate (constant (Z :. n :. n :. n)) (lift1 (\( i A.:. j A.:. k) -> cond ((&&) ((/=) k i) ((/=) k j)) 1 0))
ii :: Array A.D Ix3 e
ii = A.makeArrayR A.D A.Seq (A.Sz3 n n n) $ \(i A.:> j A.:. k) -> if k /= i && k /= j then 1 else 0
z_1 :: Matrix A.D e
z_1 = sumRowsD (w' `mulD` ii)
z_2 :: Matrix A.D e
z_2 = sumRowsD (w_1 `mulD` ii)
result = A.computeP (termDivNanD z_1 z_2)
-- | Term by term division where divisions by 0 produce 0 rather than NaN.
termDivNan :: (A.Manifest r3 a, A.Source r1 a, A.Source r2 a, Eq a, Fractional a)
=> Matrix r1 a
-> Matrix r2 a
-> Matrix r3 a
termDivNan m1 = A.compute . termDivNanD m1
termDivNanD :: (A.Source r1 a, A.Source r2 a, Eq a, Fractional a)
=> Matrix r1 a
-> Matrix r2 a
-> Matrix D a
termDivNanD m1 m2 = A.zipWith (\i j -> if j == 0 then 0 else i / j) m1 m2
sumRows :: ( A.Load r A.Ix2 e
, 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
sumRows = A.compute . sumRowsD
sumRowsD :: ( A.Source r e
, Num e
) => Array r A.Ix3 e
-> Array D A.Ix2 e
sumRowsD matrix = A.map getSum $ A.foldlWithin' 1 (\(Sum s) n -> Sum $ s + n) mempty matrix
sumRowsReferenceImplementation :: ( A.Load r A.Ix2 e
, 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
sumRowsReferenceImplementation 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.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 . mulD m1
mulD :: (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
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/LinearAlgebra/Distributional.hs 0 → 100644
View file @ 1909fd3a
{-# LANGUAGE DerivingStrategies #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeOperators #-}
{-|
Module : Gargantext.Core.LinearAlgebra.Distributional
Description : The "distributional" algorithm, fast and slow implementations
Copyright : (c) CNRS, 2017-Present
License : AGPL + CECILL v3
Maintainer : team@gargantext.org
Stability : experimental
Portability : POSIX
-}
module Gargantext.Core.LinearAlgebra.Distributional (
distributional
-- * Internals for testing
, distributionalReferenceImplementation
) where
import Data.Massiv.Array (D, Matrix, Vector, Array, Ix3, U, Ix2 (..), IxN (..))
import Data.Massiv.Array qualified as A
import Gargantext.Core.LinearAlgebra.Operations
import Prelude
-- | `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]
--
-- /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 :: forall r e. ( A.Manifest r e
, A.Manifest r Int
, A.Unbox e
, A.Source r Int
, A.Size r
, Ord e
, Fractional e
, Num e
)
=> Matrix r Int
-> Matrix U e
distributional m' = A.computeP result
where
mD :: Matrix D e
mD = A.map fromIntegral m'
m :: Matrix A.U e
m = A.compute mD
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_2 :: Matrix A.D e
d_2 = A.backpermute' (A.Sz2 diag_m_size n) (\(i A.:. _) -> i) diag_m
a :: Matrix D e
a = termDivNanD mD d_1
b :: Matrix D e
b = termDivNanD mD d_2
miDelayed :: Matrix D e
miDelayed = a `mulD` b
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
-- 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' :: Array D Ix3 e
w' = A.zipWith min w_1 w_2
z_1 :: Matrix A.D e
z_1 = A.ifoldlWithin' 1 ( \(i :> j :. k) acc w'_val ->
let ii_val = if k /= i && k /= j then 1 else 0
z1_val = w'_val * ii_val
in acc + z1_val
) 0 w'
z_2 :: Matrix A.D e
z_2 = A.ifoldlWithin' 1 ( \(i :> j :. k) acc w1_val ->
let ii_val = if k /= i && k /= j then 1 else 0
z2_val = w1_val * ii_val
in acc + z2_val
) 0 w_1
result :: Matrix A.D e
result = termDivNanD z_1 z_2
-- | A reference implementation for \"distributional\" which is slower but
-- it's more declarative and can be used to assess the correctness of the
-- optimised version.
-- Same proviso about the shape of the matri applies for this function.
distributionalReferenceImplementation :: forall r e.
( A.Manifest r e
, A.Unbox e
, A.Source r Int
, A.Size r
, Ord e
, Fractional e
, Num e
)
=> Matrix r Int
-> Matrix r e
distributionalReferenceImplementation m' = result
where
mD :: Matrix D e
mD = A.map fromIntegral m'
m :: Matrix A.U e
m = A.compute mD
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_2 :: Matrix A.D e
d_2 = A.backpermute' (A.Sz2 diag_m_size n) (\(i A.:. _) -> i) diag_m
a :: Matrix D e
a = termDivNanD mD d_1
b :: Matrix D e
b = termDivNanD mD d_2
miDelayed :: Matrix D e
miDelayed = a `mulD` b
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
-- 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' :: 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).
-- generate (constant (Z :. n :. n :. n)) (lift1 (\( i A.:. j A.:. k) -> cond ((&&) ((/=) k i) ((/=) k j)) 1 0))
ii :: Array A.D Ix3 e
ii = A.makeArrayR A.D A.Seq (A.Sz3 n n n) $ \(i A.:> j A.:. k) -> if k /= i && k /= j then 1 else 0
z_1 :: Matrix A.D e
z_1 = sumRowsD (w' `mulD` ii)
z_2 :: Matrix A.D e
z_2 = sumRowsD (w_1 `mulD` ii)
result = A.computeP (termDivNanD z_1 z_2)
src/Gargantext/Core/LinearAlgebra/Operations.hs 0 → 100644
View file @ 1909fd3a
{-# LANGUAGE DerivingStrategies #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeOperators #-}
{-|
Module : Gargantext.Core.LinearAlgebra.Operations
Description : Operations on matrixes using massiv
Copyright : (c) CNRS, 2017-Present
License : AGPL + CECILL v3
Maintainer : team@gargantext.org
Stability : experimental
Portability : POSIX
-}
module Gargantext.Core.LinearAlgebra.Operations (
-- * Convertion functions
accelerate2MassivMatrix
, accelerate2Massiv3DMatrix
, massiv2AccelerateMatrix
, massiv2AccelerateVector
-- * Operations on matrixes
, (.*)
, (.-)
, diag
, diagD
, termDivNan
, sumRows
, dim
-- * Operations on delayed arrays
, mulD
, termDivNanD
, sumRowsD
, safeDiv
-- * Internals for testing
, sumRowsReferenceImplementation
, matMaxMini
, sumM_go
, sumMin_go
) where
import Data.Array.Accelerate qualified as Acc
import Data.List.Split qualified as Split
import Data.Massiv.Array (D, Matrix, Vector, Array)
import Data.Massiv.Array qualified as A
import Prelude
import Protolude.Safe (headMay)
import Data.Monoid
-- | Converts an accelerate matrix into a Massiv matrix.
accelerate2MassivMatrix :: (A.Unbox a, Acc.Elt a) => Acc.Matrix a -> Matrix A.U a
accelerate2MassivMatrix m =
let (Acc.Z Acc.:. _r Acc.:. c) = Acc.arrayShape m
in A.fromLists' @A.U A.Par $ Split.chunksOf c (Acc.toList m)
-- | Converts a massiv matrix into an accelerate matrix.
massiv2AccelerateMatrix :: (Acc.Elt a, A.Source r a) => Matrix r a -> Acc.Matrix a
massiv2AccelerateMatrix m =
let m' = A.toLists2 m
r = Prelude.length m'
c = maybe 0 Prelude.length (headMay m')
in Acc.fromList (Acc.Z Acc.:. r Acc.:. c) (mconcat m')
-- | Converts a massiv vector into an accelerate one.
massiv2AccelerateVector :: (A.Source r a, Acc.Elt a) => A.Vector r a -> Acc.Vector a
massiv2AccelerateVector m =
let m' = A.toList m
r = Prelude.length m'
in Acc.fromList (Acc.Z Acc.:. r) m'
accelerate2Massiv3DMatrix :: (A.Unbox e, Acc.Elt e, A.Manifest r e)
=> Acc.Array (Acc.Z Acc.:. Int Acc.:. Int Acc.:. Int) e
-> A.Array r A.Ix3 e
accelerate2Massiv3DMatrix m =
let (Acc.Z Acc.:. _r Acc.:. _c Acc.:. _z) = Acc.arrayShape m
in A.fromLists' A.Par $ map (Split.chunksOf $ _z) $ Split.chunksOf (_c*_z) (Acc.toList m)
-- | Computes the diagnonal matrix of the input one.
diag :: (A.Unbox e, A.Manifest r e, A.Source r e, Num e) => Matrix r e -> Vector A.U e
diag matrix =
let (A.Sz2 rows _cols) = A.size matrix
newSize = A.Sz1 rows
in A.makeArrayR A.U A.Seq newSize $ (\(A.Ix1 i) -> matrix A.! (A.Ix2 i i))
diagD :: (A.Source r e, A.Size r) => Matrix r e -> Vector A.D e
diagD matrix =
let (A.Sz2 rows _cols) = A.size matrix
newSize = A.Sz1 rows
in A.backpermute' newSize (\i -> i A.:. i) matrix
-- | Term by term division where divisions by 0 produce 0 rather than NaN.
termDivNan :: (A.Manifest r3 a, A.Source r1 a, A.Source r2 a, Eq a, Fractional a)
=> Matrix r1 a
-> Matrix r2 a
-> Matrix r3 a
termDivNan m1 = A.compute . termDivNanD m1
termDivNanD :: (A.Source r1 a, A.Source r2 a, Eq a, Fractional a)
=> Matrix r1 a
-> Matrix r2 a
-> Matrix D a
termDivNanD m1 m2 = A.zipWith safeDiv m1 m2
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
, 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
sumRows = A.compute . sumRowsD
sumRowsD :: ( A.Source r e
, Num e
) => Array r A.Ix3 e
-> Array D A.Ix2 e
sumRowsD matrix = A.map getSum $ A.foldlWithin' 1 (\(Sum s) n -> Sum $ s + n) mempty matrix
sumRowsReferenceImplementation :: ( A.Load r A.Ix2 e
, 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
sumRowsReferenceImplementation 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.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 . mulD m1
mulD :: (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
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))
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]
]
test/Test/Core/LinearAlgebra.hs
View file @ 1909fd3a
......@@ -112,7 +112,8 @@ tests = testGroup "LinearAlgebra" [
, testProperty "sumM_go" compareSumM_go
, testProperty "sumMin_go" compareSumMin_go
, testGroup "distributional" [
testProperty "2x2" (compareDistributional (Proxy @Double) twoByTwo)
testProperty "reference implementation roundtrips" compareDistributionalImplementations
, 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))
......@@ -145,6 +146,11 @@ compareSumRows i1
in counterexample "sumRows and reference implementation do not agree" (massiv === massiv') .&&.
accelerate === LA.massiv2AccelerateMatrix massiv
compareDistributionalImplementations :: SquareMatrix Int -> Property
compareDistributionalImplementations (SquareMatrix i1) =
let ma = LA.accelerate2MassivMatrix i1
in LA.distributional @Massiv.U @Double ma === LA.distributionalReferenceImplementation ma
compareDistributional :: forall e.
( Eq e
, Show e
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment