improve performance

app/Main.hs

@@ -16,8 +16,6 @@ import Prelude (String)
 import Protolude
 import qualified Data.IntMap                       as Dict
 import qualified Data.List as List
-import qualified Eigen.Matrix as Matrix
-import qualified Eigen.SparseMatrix as SMatrix
 import qualified Prelude as Prelude
 import qualified Data.Graph.Inductive.PatriciaTree as DGIP
 
@@ -31,9 +29,11 @@ setupEnv _    = getUnlabGraph Random
 main :: IO ()
 main = do
   fp <- getArgs
-
-  withArgs [] $
-    defaultMain
-      [ env (snd <$> setupEnv fp) $ \ ~g ->
-        bench "bench" (nf (\ x -> withG x (parts . clusteringOptim 3 Conf)) g)
-      ]
+  setupEnv fp >>= \(_, ~g) ->
+    evaluate $ withG g (clusteringOptim 3 Conf)
+  return ()
+  -- withArgs [] $
+  --   defaultMain
+  --     [ env (snd <$> setupEnv fp) $ \ ~g ->
+  --       bench "bench" (nf (\ x -> withG x (parts . clusteringOptim 3 Conf)) g)
+  --     ]

gargantext-graph.cabal

 cabal-version: 1.12
 
--- This file has been generated from package.yaml by hpack version 0.33.0.
+-- This file has been generated from package.yaml by hpack version 0.34.4.
 --
 -- see: https://github.com/sol/hpack
 --
@@ -9,8 +9,8 @@ cabal-version: 1.12
 name:           gargantext-graph
 version:        0.1.0.0
 description:    Please see the README on GitHub at <https://github.com/githubuser/gargantext-graph#readme>
-homepage:       https://github.com/https://gitlab.iscpif.fr/gargantext/gargantext-graph#readme
-bug-reports:    https://github.com/https://gitlab.iscpif.fr/gargantext/gargantext-graph/issues
+homepage:       https://github.com/gargantext/gargantext-graph#readme
+bug-reports:    https://github.com/gargantext/gargantext-graph/issues
 author:         Alexandre Delanoë
 maintainer:     alexandre+dev@delanoe.org
 copyright:      2021 CNRS / Alexandre Delanoë
@@ -23,12 +23,12 @@ extra-source-files:
 
 source-repository head
   type: git
-  location: https://github.com/https://gitlab.iscpif.fr/gargantext/gargantext-graph
+  location: https://github.com/gargantext/gargantext-graph
 
 library
   exposed-modules:
       Data.Array.Accelerate.Utils
-      Data.Eigen.Coeff
+      -- Data.Eigen.Coeff
       Graph.BAC.Clustering
       Graph.BAC.Proxemy
       Graph.BAC.ProxemyOptim
@@ -45,7 +45,10 @@ library
       Paths_gargantext_graph
   hs-source-dirs:
       src
-  default-extensions: DataKinds FlexibleInstances OverloadedStrings
+  default-extensions:
+      DataKinds
+      FlexibleInstances
+      OverloadedStrings
   build-depends:
       accelerate
     , accelerate-arithmetic
@@ -55,14 +58,16 @@ library
     , cassava
     , cereal
     , containers
-    , eigen
+    -- , eigen
     , fgl
     , haskell-igraph >=0.6.0
+    , hmatrix
     , lens
     , pretty-simple
     , protolude
     , reflection
-    , singletons
+    , singletons == 2.7
+    , sparse-linear
     , string-conversions
     , text
     , vector
@@ -86,7 +91,7 @@ executable gargantext-graph-exe
     , cereal
     , containers
     , criterion
-    , eigen
+    -- , eigen
     , fgl
     , gargantext-graph
     , haskell-igraph >=0.6.0
@@ -107,7 +112,10 @@ test-suite gargantext-graph-test
       Paths_gargantext_graph
   hs-source-dirs:
       test
-  default-extensions: DataKinds FlexibleInstances OverloadedStrings
+  default-extensions:
+      DataKinds
+      FlexibleInstances
+      OverloadedStrings
   ghc-options: -threaded -rtsopts -with-rtsopts=-N
   build-depends:
       QuickCheck
@@ -119,7 +127,7 @@ test-suite gargantext-graph-test
     , cassava
     , cereal
     , containers
-    , eigen
+    -- , eigen
     , fgl
     , gargantext-graph
     , haskell-igraph >=0.6.0

package.yaml

 name:                gargantext-graph
 version:             0.1.0.0
-github:              "https://gitlab.iscpif.fr/gargantext/gargantext-graph"
+github:              gargantext/gargantext-graph
 license:             AGPL
 author:              "Alexandre Delanoë"
 maintainer:          "alexandre+dev@delanoe.org"

src/Graph/BAC/ProxemyOptim.hs

+{-# LANGUAGE BangPatterns #-}
 {-| Module      : Gargantext.Core.Viz.Graph.ProxemyOptim
 Description : Proxemy
 Copyright   : (c) CNRS, 2017-Present
@@ -30,35 +31,33 @@ Gaume.
            , MultiParamTypeClasses
 #-}
 
+{-# LANGUAGE AllowAmbiguousTypes #-}
 module Graph.BAC.ProxemyOptim
   where
 
---import Debug.SimpleReflect
-import Data.Eigen.Coeff (coeffSM, coeffDM)
 import Data.IntMap (IntMap)
 import Data.Maybe (isJust)
 import Data.Proxy (Proxy(Proxy))
 import Data.Reflection
-import Eigen.Internal (CTriplet(..), Elem(..), toC, fromC, C(..), natToInt, Row(..), Col(..))
-import Eigen.Matrix (sum)
-import Eigen.SparseMatrix (SparseMatrix, SparseMatrixXd, (!), toMatrix)
 import GHC.TypeLits (KnownNat, Nat, SomeNat(SomeNat), type(+), natVal, sameNat, someNatVal)
 import Graph.FGL
 import Graph.Types
-import Prelude (String, readLn)
-import Protolude hiding (sum, natVal)
-import qualified Eigen.Matrix as DenseMatrix
+import Prelude (String, readLn, error)
+import Protolude hiding (sum, natVal, trace)
 import qualified Data.Graph.Inductive              as DGI
 import qualified Data.Graph.Inductive.PatriciaTree as DGIP
 import qualified Data.List                         as List
 import qualified Data.IntMap                       as Dict
 import qualified Data.Vector                       as V
 import qualified Data.Vector.Storable              as VS
-import qualified Eigen.Matrix                      as DMatrix
-import qualified Eigen.SparseMatrix                as SMatrix
 import qualified Prelude                           as Prelude
 import qualified Data.Set                          as Set
 
+import qualified Data.Matrix.Sparse.Static         as SMatrix
+import qualified Data.Vector.Sparse.Static         as SV
+import qualified Numeric.LinearAlgebra.Static      as DMatrix
+import Debug.Trace
+import qualified Data.Vector.Unboxed as UV
 
 ----------------------------------------------------------------
 
@@ -68,23 +67,25 @@ clusteringOptim :: forall n a b. (KnownNat n, Show a, Show b)
                 -> Similarity
                 -> FiniteGraph n a b
                 -> Clustering Node
-clusteringOptim l s fg@(FiniteGraph g) = -- traceShow ("mc", mc) $
+clusteringOptim l s fg@(FiniteGraph g) = trace "clusteringOptim" $ -- traceShow ("mc", mc) $
   make_clust_part sorted_edges matq
     where
-      adj    = adjacent fg False
-      tra    = transition adj
+      adj    = trace "adjacent False" $ adjacent fg False
+      tra    = trace "tra" $ transition adj
 
-      adj'   = adjacent fg True
-      tra'    = transition adj'
-      sorted_edges = sort_edges (natToInt @n)
-                   $ edges_confluence l adj' tra'
+      adj'   = trace "adjacent True" $ adjacent fg True
+      tra'    = trace "tra'" $ transition adj'
+      sorted_edges = adj `seq` tra `seq` adj' `seq` tra' `seq`
+                     (trace "sorted_edges"
+                   $ sort_edges (natToInt @n)
+                   $ edges_confluence l adj' tra')
 
       mc    = matconf False adj (proxemie l tra)
       md    = matmod fg
 
       matq = case s of
-        Conf -> SimConf mc -- $ matconf False adj (proxemie l tra)
-        Mod  -> SimMod  md -- $ matmod fg
+        Conf -> SimConf mc --  $ matconf False adj (proxemie l tra)
+        Mod  -> SimMod  md --  $ matmod fg
 
 -- | Where main Types are defined as
 data Similarity = Conf | Mod
@@ -114,8 +115,8 @@ type NeighborsFilter a b = DGI.Gr a b -> Node -> [Node]
 type RmEdge      = Bool
 ---------------------------------------------------------------
 -- Data Structure
-type VectorS          n = SparseMatrix 1 n Double
-type VectorD          n = DenseMatrix.Matrix 1 n Double
+type VectorS          n = SMatrix.Matrix 1 n Double
+type VectorD          n = DMatrix.L 1 n
 
 type AdjacencyMatrix  n = MatrixS n
 type TransitionMatrix n = MatrixS n
@@ -133,26 +134,24 @@ adjacent :: KnownNat n
          => FiniteGraph n a b
          -> IsReflexive
          -> AdjacencyMatrix n
-adjacent (FiniteGraph g) isReflexive =
-  SMatrix.fromVector $ VS.fromList $ triplets <> diag
+adjacent (FiniteGraph g) isReflexive = trace "adjacent"
+  SMatrix.fromList $ triplets <> diag
   where
-    triplets = [ CTriplet (toC i) (toC j) 1.0
+    triplets = [ (i, j, 1.0)
                | i <- nodes g
                , j <- neighbors g i
                , i /= j
                ]
 
     diag = case isReflexive of
-      True -> [ CTriplet (toC n) (toC n) 1.0
-             | n <- nodes g
-             ]
+      True  -> [ (n, n, 1.0) | n <- nodes g ]
       False -> []
 
 
-transition :: KnownNat         n
+transition :: (HasCallStack, KnownNat n)
            => AdjacencyMatrix  n
            -> TransitionMatrix n
-transition m = SMatrix.imap (\i j v -> v * (VS.!) s i) m
+transition m = trace "transition" $ SMatrix.imap (\i j v -> v * (SV.!) s i) m
   where
     s = sumWith Colonne (\s -> 1 / s) m
 
@@ -161,10 +160,28 @@ proxemie :: KnownNat         n
          => Length
          -> TransitionMatrix n
          -> ProxemyMatrix    n
-proxemie l tm = case l <= 1 of
-  True  -> toMatrix tm
-  False -> toMatrix $ foldl' (\m' _-> SMatrix.mul m' tm) tm [2 .. (l :: Int)]
-
+proxemie l tm = trace ("proxemie (l = " ++ show l ++ ")") $ case l <= 1 of
+  True  -> SMatrix.densify tm
+  False -> matInfo "proxemie tm" tm $
+    case matPow tm (l-1) of
+      tm' -> matInfo "proxemie tm'" tm' $ SMatrix.densify tm'
+
+
+-- https://en.wikipedia.org/wiki/Exponentiation_by_squaring
+matPow :: MatrixS n -> Int -> MatrixS n
+matPow mat n
+  | n <= 0 = error "matPow: positive exponent expected"
+  | otherwise = f mat n
+  where
+    f a i
+      | even i = f (a `mul` a) (i `quot` 2)
+      | i == 1 = a
+      | otherwise = g (a `mul` a) (i `quot` 2) a
+    g a i b
+      | even i = g (a `mul` a) (i `quot` 2) b
+      | i == 1 = a `mul` b
+      | otherwise = g (a `mul` a) (i `quot` 2) (a `mul` b)
+    mul = SMatrix.mul
 
 ---------------------------------------------------------------
 matconf :: forall n. KnownNat n
@@ -172,34 +189,35 @@ matconf :: forall n. KnownNat n
         -> AdjacencyMatrix  n
         -> ProxemyMatrix    n
         -> ConfluenceMatrix n
-matconf False a p = symmetry confmat
+matconf False a p = seq confmat $ trace "matconf" $ symmetry confmat
   where
     -- vcount  = natToInt @n
-    degs    = sumWith Colonne identity a
+    degs    = DMatrix.extract $ SV.toDense (sumWith Colonne identity a)
     sumdeg  = VS.sum degs
-    confmat = DMatrix.imap (\  x   y   v ->
-                            if x < y
-                              then let
-                                      prox_y_x_length = v
-                                      prox_y_x_infini = ((VS.!) degs x) / sumdeg
-                                 in 
-                                    (prox_y_x_length - prox_y_x_infini)
-                                  / (prox_y_x_length + prox_y_x_infini)
-                              else 0
-                        ) $ DMatrix.transpose p
+    tp      = DMatrix.tr p
+    confmat = seq tp . trace "confmat" $
+      DMatrix.imapL
+        (\(x, y) v ->
+            if x < y
+            then let prox_y_x_length = v
+                     prox_y_x_infini = ((VS.!) degs x) / sumdeg
+                 in
+                   (prox_y_x_length - prox_y_x_infini)
+                 / (prox_y_x_length + prox_y_x_infini)
+            else 0
+        ) tp
 matconf True _a _p = panic "MatConf True: TODO but not needed for now"
 
 
-
 matmod :: forall n a b. KnownNat n => FiniteGraph n a b -> ModularityMatrix n
-matmod fg = symmetry $ toMatrix modmat
+matmod fg = trace "matmod" $ symmetry $ SMatrix.densify modmat
   where
     n' = natToInt @n
     a  = adjacent fg False
-    sumRows = sumWith Ligne   identity a
-    sumCols = sumWith Colonne identity a
-    ecount = sum $ toMatrix a
-    modmat = SMatrix.imap (\x y v -> 
+    !sumRows = DMatrix.extract $ SV.toDense (sumWith Ligne   identity a)
+    !sumCols = DMatrix.extract $ SV.toDense (sumWith Colonne identity a)
+    !ecount = SMatrix.sum a
+    modmat = SMatrix.imap (\x y v ->
                                 if x < y
                                 then v - ((VS.!) sumRows x * (VS.!) sumCols y) / (2 * ecount)
                                 else 0
@@ -209,36 +227,60 @@ matmod fg = symmetry $ toMatrix modmat
 type UnsortedEdges = [(Node, Node, Double)]
 type SortedEdges   = [(Node, Node, Double)]
 
-edges_confluence :: forall           n a b
+matInfo
+  :: forall b (n :: Nat) (m :: Nat).
+     (KnownNat n, KnownNat m)
+  => String
+  -> SMatrix.Matrix n m Double
+  -> b -> b
+matInfo name mat = trace $
+  name ++ ": " ++ show nnz ++ "/" ++ show total ++
+  " (~" ++ pcstr ++ "%) nnz, " ++
+  szInfo ++
+  "\n\t= [ " ++ intercalate ", " snip ++ " , ... ]"
+  where total = n*m
+        nnz = SMatrix.nonZeros mat
+        pc = fromIntegral (100 * nnz) / fromIntegral total
+        pcstr = take 4 (show pc)
+        n = natToInt @n
+        m = natToInt @m
+        szInfo = show n ++ "x" ++ show m
+        snip = [ show (coeffSM 0 i mat) | i <- [1..10] ]
+
+edges_confluence :: forall           n
                   . KnownNat         n
                  => Length
                  -> AdjacencyMatrix  n
                  -> TransitionMatrix n
                  -> UnsortedEdges
-edges_confluence l am tm = -- traceShow ("degs", degs) $
+edges_confluence l am tm = trace ("edges_confluence (dim = " ++ show vcount ++ ")") $ -- traceShow ("degs", degs) $
+  matInfo "am" am $ matInfo "tm" tm $
   SMatrix.toList matconf'
     where
       vcount   = natToInt @n
-      degs     = sumWith Colonne identity am
-      sumdeg   = VS.sum degs
-      matconf' = SMatrix.imap
-                 (\x y _ ->
+      degs     = DMatrix.extract $ SV.toDense (sumWith Colonne identity am)
+      !sumdeg   = VS.sum degs
+      !tmam     = SMatrix.zipR 0 (,) am tm
+      matconf' = seq sumdeg $ trace "edges_confluence.matconf'" $ SMatrix.imap
+                 (\x y _v ->
                    if x < y
                       then
                         let
-                            deg_x = (VS.!) degs x - 1
-                            deg_y = (VS.!) degs y - 1
-
-                            tm' = SMatrix.imap (\i j v -> if (i == x && j == y) || ( i == y && j == x)
-                                                            then 0
-                                                            else
-                                                              if i == x && i /= y
-                                                                then (coeffSM i j am) / deg_x
-                                                                else
-                                                                   if i == y && i /= x
-                                                                     then (coeffSM i j am) / deg_y
-                                                                     else v
-                                               ) tm
+                            !deg_x = (VS.!) degs x - 1
+                            !deg_y = (VS.!) degs y - 1
+
+                            !tm' = {- trace ("tm': " ++ show (x, y)) $ -} SMatrix.imap
+                              (\i j (am_ij, tm_ij) ->
+                                 if (i == x && j == y) || ( i == y && j == x)
+                                 then 0
+                                 else
+                                   if i == x && i /= y
+                                   then am_ij / deg_x
+                                   else
+                                     if i == y && i /= x
+                                     then am_ij / deg_y
+                                     else tm_ij
+                              ) tmam
 
                             v  = fromList (Proxy @n) [(fromIntegral y,1)]
                             v' = doProx l v tm'
@@ -264,22 +306,18 @@ edges_confluence l am tm = -- traceShow ("degs", degs) $
           True  -> panic "doProx"
           False -> foldl' (\v'' _-> SMatrix.mul v'' tm'') v' [1 .. (l :: Int) ]
 
-
--- | TODO optimization
 sort_edges :: Int
            -> UnsortedEdges
            -> SortedEdges
-sort_edges n =  List.concat
-             . (map (List.sortOn (\(x,y,_) -> x * n + y)))
-             . (List.groupBy (\x y -> third x == third y))
-             . List.reverse
-             . (List.filter (\(x,y,_) -> x < y))
-             . (List.sortOn third)
+sort_edges n = trace "sort_edges"
+             . List.sortBy (\a b -> comparing third a b <> comparing xnpy a b)
+             . List.filter (\(x,y,_) -> x < y)
   where
     third
       :: forall a b c
        . (a,b,c) -> c
     third (_,_,c) = c
+    xnpy (x,y,_) = x*n+y
 
 
 ---------------------------------------------------------------
@@ -290,18 +328,18 @@ updateClustering c@(ClusteringIs parts idx currentScore _) f x y =
   let modX = fromMaybe 0 $ Dict.lookup x idx
       modY = fromMaybe 0 $ Dict.lookup y idx
    in case x == y || modX == modY of
-    True  -> c -- do case x' or y' are Nothing 
-    False -> let c' = updateWith c f (x,modX) (y,modY)
+    True  -> c -- do case x' or y' are Nothing
+    False -> let c' = updateWith c f x modX y modY
               in case score c' >= currentScore of
                    True  -> c'
                    False -> c
 
 updateWith :: Clustering Node
            -> (Int -> Int -> Double)
-           -> (Int,Int)
-           -> (Int,Int)
+           -> Int -> Int
+           -> Int -> Int
            -> Clustering Node
-updateWith c@(ClusteringIs parts idx _ _) f (x,modX) (y,modY) = -- traceShow ("score", x,y,score') $
+updateWith c@(ClusteringIs parts idx _ _) f x modX y modY = trace "updateWith" -- traceShow ("score", x,y,score') $
   ClusteringIs parts' idx' score' Part
     where
       parts' = Dict.filter (not . Set.null)
@@ -319,12 +357,21 @@ updateWith c@(ClusteringIs parts idx _ _) f (x,modX) (y,modY) = -- traceShow ("s
 
       px = Dict.elems parts'
 
-      score' = Prelude.sum $ List.concat
-                           $ map (\s -> [ if x'' < y'' then f x'' y'' else 0 -- to Validate with BG
-                                        | x'' <- Set.toList s
-                                        , y'' <- Set.toList s
-                                        ]
-                                  ) px
+      score' = getSum $
+        foldMap (\s -> fold [ Sum (f x'' y'')
+                            | x'' <- Set.toList s
+                            , y'' <- Set.toList s
+                            , x'' < y''
+                            ]
+                ) parts'
+
+        -- Prelude.sum $ List.concatMap
+        -- (\s -> [ f x'' y'' -- to Validate with BG
+        --        | x'' <- Set.toList s
+        --        , y'' <- Set.toList s
+        --        , x'' < y''
+        --        ]
+        -- ) px
 
 ---------------------------------------------------------------
 make_clust_part :: forall n
@@ -332,8 +379,8 @@ make_clust_part :: forall n
                 => SortedEdges
                 -> SimilarityMatrix n
                 -> Clustering Node
-make_clust_part se sm = foldl' (\c (e1,e2,_)
-                                 -> updateClustering c
+make_clust_part se sm = trace "make_clust_part" $ foldl' (\c (e1,e2,_)
+                                 -> {- trace ("foldl' step: " ++ show (e1, e2)) $ -} updateClustering c
                                     (\x y -> 2 * (coeffDM x y sm'))
                                     e1 e2
                                ) (ClusteringIs parts idx 0 Part) se
@@ -398,7 +445,7 @@ instance KnownNat n => Symmetry (MatrixS n) where
   symmetry = symmetryS
 
 symmetryD :: KnownNat n => MatrixD n -> MatrixD n
-symmetryD m = DMatrix.imap (\x y v -> if x < y then v else coeffDM y x m) m
+symmetryD m = DMatrix.imapL (\(x, y) v -> if x < y then v else coeffDM y x m) m
 
 symmetryS :: KnownNat n => MatrixS n -> MatrixS n
 symmetryS m = SMatrix.imap (\x y v -> if x < y then v else coeffSM y x m) m
@@ -421,21 +468,30 @@ vectorFromListS pn ns
        n = natVal pn
 
 vectorFromListD :: KnownNat n => Proxy n -> [(Integer, Double)] -> VectorD n
-vectorFromListD pn ns = toMatrix $ vectorFromListS pn ns
+vectorFromListD pn ns = SMatrix.densify $ vectorFromListS pn ns
 
 ------------------------------
 data Direction = Ligne | Colonne
 
-sumWith :: ( Elem a
-           , Elem t
-           , KnownNat n
-           ) => Direction -> (t -> a) -> SparseMatrix n n t -> VS.Vector a
-sumWith dir f m = VS.fromList
-       $ map (\v -> f v)
-       $ case dir of
-           Colonne -> somme $ SMatrix.getCols m
-           Ligne   -> somme $ SMatrix.getRows m
-          where
-            somme m' = map (sum . SMatrix.toMatrix) m'
+sumWith
+  :: (UV.Unbox a, UV.Unbox t, KnownNat n, Num t, HasCallStack, Show t)
+  => Direction -> (t -> a) -> SMatrix.Matrix n n t -> SV.V n a
+sumWith dir f m = SV.map f $
+  case dir of
+    Colonne -> SV.sum $ SMatrix.getCols m
+    Ligne   -> SV.sum $ SMatrix.getRows m
 
 --------------------------------------------
+
+natToInt :: forall (n :: Nat). KnownNat n => Int
+natToInt = fromIntegral $ natVal (Proxy :: Proxy n)
+
+coeffSM
+  :: (Num a, UV.Unbox a, Show a)
+  => Int -> Int -> SMatrix.Matrix n p a -> a
+coeffSM i j m = SMatrix.extractCol m j SV.! i
+
+coeffDM
+  :: (Num a, KnownNat n, KnownNat p)
+  => Int -> Int -> DMatrix.L n p -> Double
+coeffDM i j m = m `DMatrix.at` (i, j)

src/Graph/Tools/CSV.hs

@@ -22,7 +22,6 @@ import Data.Vector hiding (map, uniq)
 import Prelude (read)
 import Protolude
 import Graph.Types
-import qualified Eigen.SparseMatrix as SMatrix
 import qualified Data.Map                          as Map
 import qualified Data.ByteString.Lazy              as BL
 import qualified Data.Graph.Inductive              as DGI

src/Graph/Tools/Import.hs

+{-# LANGUAGE TypeApplications #-}
 {-| Module      : Graph.Tools
 Description :
 Copyright   : (c) CNRS, Alexandre Delanoë
@@ -25,7 +26,7 @@ import Graph.Types
 import Graph.Tools.Random
 import Graph.Tools.CSV
 import Graph.FGL
-import qualified Eigen.SparseMatrix as SMatrix
+import qualified Data.Matrix.Sparse.Static         as SMatrix
 import qualified Data.IntMap                       as Dict
 import qualified Data.Map                          as Map
 import qualified Data.ByteString.Lazy              as BL
@@ -44,7 +45,7 @@ data GraphData = LightGraph    { lightGraph :: Graph () () }
 
 
 getGraph :: GetGraph -> IO GraphData
-getGraph Random = randomAdjacency
+getGraph Random = randomAdjacency @100
                 >>= \m -> pure $ LightGraph
                                $ mkGraphUfromEdges
                                $ List.map (\(x,y,_) -> (x,y))

src/Graph/Tools/Random.hs

@@ -6,7 +6,6 @@ Maintainer  : alexandre+dev@delanoe.org
 Stability   : experimental
 Portability : POSIX
 
-
 -}
 
 {-# LANGUAGE ConstraintKinds #-}
@@ -22,18 +21,21 @@ import Data.Vector hiding (map, uniq)
 import Prelude (read)
 import Protolude
 import Graph.Types
-import qualified Eigen.SparseMatrix                as SMatrix
-import qualified Eigen.Matrix                      as DMatrix
+-- import qualified Eigen.SparseMatrix                as SMatrix
+-- import qualified Eigen.Matrix                      as DMatrix
 import qualified Data.Graph.Inductive.PatriciaTree as DGIP
 import qualified Data.Set                          as Set
 import qualified Data.Vector                       as Vector
 
+import qualified Data.Matrix.Sparse.Static    as SMatrix
+import qualified Numeric.LinearAlgebra.Static as DMatrix
+
 -- Random Matrix
 -- | Random Matrix && Graph
 -- TODO random matrix of any size for the tests
 
-randomMatrix :: IO (MatrixD 100)
-randomMatrix = DMatrix.random
+randomMatrix :: KnownNat n => IO (MatrixD n)
+randomMatrix = DMatrix.rand
 
 {-
 matrix2graph :: forall   n
@@ -44,10 +46,12 @@ matrix2graph m = withG (mkGraphUfromEdges $ map (\(x,y,_) -> (x,y)) $ SMatrix.to
                        identity
 -}
 
-randomAdjacency :: IO (MatrixS 100)
+randomAdjacency :: KnownNat n => IO (MatrixS n)
 randomAdjacency = do
   m1 <- randomMatrix
   m2 <- randomMatrix
-  pure $ SMatrix.fromMatrix
-       $ DMatrix.imap (\i j v -> if i < j && v > 0.9 then 1 else 0)
+  pure $ SMatrix.sparsify (\(i, j) v -> if i < j && v > 0.9 then Just 1 else Nothing)
        $ DMatrix.mul m1 m2
+  -- pure $ SMatrix.fromMatrix
+  --      $ DMatrix.imapL (\(i, j) v -> if i < j && v > 0.9 then 1 else 0)
+  --      $ DMatrix.mul m1 m2

src/Graph/Types.hs

@@ -15,8 +15,12 @@ module Graph.Types where
 
 import qualified Data.Graph.Inductive.PatriciaTree as DGIP
 import Data.IntMap (IntMap)
-import qualified Eigen.Matrix as DenseMatrix
-import Eigen.SparseMatrix (SparseMatrix)
+-- import qualified Eigen.Matrix as DenseMatrix
+-- import Eigen.SparseMatrix (SparseMatrix)
+
+import qualified Data.Matrix.Sparse.Static    as Sparse
+import qualified Data.Vector.Unboxed          as VU
+import qualified Numeric.LinearAlgebra.Static as Dense
 
 -- | Main Types use in this libray
 
@@ -26,8 +30,5 @@ type Dict = IntMap
 type Graph a b = DGIP.Gr a b
 
 -- | Type for Matrix computation optimizations (with Eigen)
-type MatrixD          n = DenseMatrix.Matrix n n Double
-type MatrixS          n = SparseMatrix n n Double
-
-
-
+type MatrixD          n = Dense.L n n
+type MatrixS          n = Sparse.Matrix n n Double

stack.yaml

@@ -44,8 +44,16 @@ extra-deps:
   commit: f5c0e0071ec7b6532f9a9cd3eb33d14f340fbcc9
 - git: https://gitlab.iscpif.fr/anoe/accelerate-utility.git
   commit: 83ada76e78ac10d9559af8ed6bd4064ec81308e4
+- git: https://github.com/alpmestan/sparse-linear.git
+  commit: 785f12c99615907b207f5499dd3b70a486a0249b
+  subdirs:
+    - sparse-linear
+- git: https://github.com/alpmestan/hmatrix.git
+  commit: 39155c60c97ba8cd5c8ab9c202b428fec42faa3c
+  subdirs:
+    - packages/base
 - accelerate-arithmetic-1.0.0.1@sha256:555639232aa5cad411e89247b27871d09352b987a754230a288c690b6de6d888,2096
-- eigen-3.3.7.0@sha256:7c24a86426b238e02ba0ac2179cc44997518b443a13828ed7a791fe49a6dffa5,82060
+# - eigen-3.3.7.0@sha256:7c24a86426b238e02ba0ac2179cc44997518b443a13828ed7a791fe49a6dffa5,82060
 
 
 # Override default flag values for local packages and extra-deps
@@ -64,7 +72,6 @@ extra-deps:
 # Override the architecture used by stack, especially useful on Windows
 # arch: i386
 # arch: x86_64
-#
 # Extra directories used by stack for building
 # extra-include-dirs: [/path/to/dir]
 # extra-lib-dirs: [/path/to/dir]