Persistence 1.0 → 1.1

raw patch · 7 files changed

+311/−45 lines, 7 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• ChangeLog.md +11/−1
• HasseDiagram.hs +146/−0
• Matrix.hs +31/−23
• Persistence.cabal +3/−3
• Persistence.hs +82/−3
• SimplicialComplex.hs +22/−11
• Util.hs +16/−4

@@ -1,5 +1,15 @@ # Revision history for Persistence -## 1.0  -- YYYY-mm-dd+## 1.0  -- 2018-05-11  * First version. Released on an unsuspecting world.++## 1.1  -- 2018-05-27++### Added+- Bottleneck distance, a way to compare the topology of different data sets.+- HasseDiagram module, will allow users to deduce topological features of information flow in networks.+- A function for encoding generic graphs as 2D vectors.++### Changed+- Improved documentation for all exposed modules.

@@ -0,0 +1,146 @@+{- |+Module     : Persistence.HasseDiagram+Copyright  : (c) Eben Cowley, 2018+License    : BSD 3 Clause+Maintainer : eben.cowley42@gmail.com+Stability  : experimental++This module implements algorithms for admissible Hasse diagrams. A Hasse diagram is admissible if it is stratified and oriented. A diagram is stratified if all the vertices can be arranged in rows such that all the sources of each vertex are in the next highest row and all the targets are in the next lowest row. A diagram is oriented if every vertex has a linear ordering on its targets.++A node in the diagram is represented as a tuple: the indices of the level 0 nodes in the diagram that are reachable from this node, the indices of targets in the next lowest level, and the indices of the sources in the next highest level. The entire diagram is simply an array where each entry is an array representing a particular level; index 0 represents level 0, etc.++Any directed graph can be encoded as an admissible Hasse diagram with 2 levels. The edges are level 1 and the vertices are level 0. The ordering on the targets of a node representing an edge is simply the terminal vertex first and the initial vertex second. This may be counterintuitive, but its helpful to interpret an arrow between two vertices as the "<" operator. This induces a linear ordering on the vertices of any acyclic complete subgraph - which is what the nodes in the Hasse diagram of the directed clique complex represent.++Any oriented simplicial complex can also be encoded as an admissible Hasse diagram. A node is a simplex, the targets are the faces of the simplex, and the sources are simplices of which the given simplex is a face.++The main feature of this module is an algorithm which takes the Hasse diagram of a directed graph and generates the Hasse diagram of the directed flag complex - the simplicial complex whose simplices are acyclic complete subgraphs of the given graph. Here acyclic refers to a directed graph without any sequence of arrows whose heads and tails match up an which has the same start and end vertex.++The idea is that, if your directed graph represents any kind of information flow, "sub-modules" in the network that simply take input, process it using its nodes, and then output it without spinning the information around at all. These "sub-modules" are the directed cliques/flags which I've been referring to as acyclic complete subgraphs up to this point. Constructing a simplicial complex out of them will allow you to both simplify the 1-dimensional topology of the network and possibly detect higher-dimensional topological features.++-}++module HasseDiagram+  ( Node+  , HasseDiagram+  , hsd2String+  , encodeDirectedGraph+  , directedFlagComplex+  , toSimplicialComplex+  ) where++import Util+import SimplicialComplex++import Data.List as L+import Data.Vector as V++{- |+  Type representing a node in a Hasse diagram.+  Hasse diagrams are being used to represent simplicial complexes so each node represents a simplex+  Contents of the tuple in order: Vector of references to vertices of the underlying directed graph,+  vector of references to the simplices faes in the next lowest level of the Hasse diagram,+  vector of references to "parent" simplices (simplices who have this simplex as a face) in the next highest level of the Hasse diagram.+-}+type Node = (Vector Int, Vector Int, Vector Int)++-- | Type representing an admissible Hasse diagram. Each entry in the vector represents a level in the Hasse diagram.+type HasseDiagram = Vector (Vector Node)++-- | Simple printing function for Hasse diagrams.+hsd2String :: HasseDiagram -> String+hsd2String = (L.intercalate "\n\n") . V.toList . (V.map (L.intercalate "\n" . V.toList . V.map show))++-- | Given the number of vertices in a directed graph, and pairs representing the direction of each edge (initial, terminal), construct a Hasse diagram representing the graph.+encodeDirectedGraph :: Int -> [(Int, Int)] -> HasseDiagram+encodeDirectedGraph numVerts cxns =+  let verts       = V.map (\n -> (n `cons` V.empty, V.empty, V.empty)) $ 0 `range` (numVerts - 1)++      encodeEdges _ vertices edges []          = V.empty `snoc` vertices `snoc` edges+      encodeEdges n vertices edges ((i, j):xs) =+        let v1 = vertices ! i; v2 = vertices ! j; edge = V.empty `snoc` j `snoc` i+        in encodeEdges (n + 1)+          (replaceElem i (one v1, two v1, (thr v1) `snoc` n) $+            replaceElem j (one v2, two v2, (thr v2) `snoc` n) vertices) +              (edges `snoc` (edge, edge, V.empty)) xs++  in encodeEdges 0 verts V.empty cxns++-- | Given a Hasse diagram representing a directed graph, construct the diagram representing the directed clique/flag complex of the graph.+directedFlagComplex :: HasseDiagram -> HasseDiagram+directedFlagComplex directedGraph =+  let edges    = V.last directedGraph+      fstSinks =+        V.map (\e ->+          V.map (\(e0, _) -> (two e0) ! 0) $+            findBothElems (\e1 e2 -> (two e1) ! 0 == (two e2) ! 0)+              (V.filter (\e0 -> (two e0) ! 1 == (two e) ! 1) edges) (V.filter (\e0 -> (two e0) ! 1 == (two e) ! 0) edges)) edges++      --take last level of nodes and their sinks. return modified last level, new level, and new sinks+      makeLevel :: Bool -> HasseDiagram -> Vector Node -> Vector (Vector Int) -> (Vector Node, Vector Node, Vector (Vector Int))+      makeLevel fstIter result oldNodes oldSinks =+        let maxindex = V.length oldNodes++            --given a node and a specific sink, construct a new node with new sinks that has the given index. Fst output is the modified input nodes, snd output is the new node, thrd output is the sinks of the new node+            makeNode :: Int -> Int -> Int -> Vector Node -> Vector Int -> (Vector Node, Node, Vector Int)+            makeNode newIndex oldIndex sinkIndex nodes sinks =+              let sink     = sinks ! sinkIndex+                  oldNode  = nodes ! oldIndex+                  verts    = sink `cons` (one oldNode) --the vertices of the new simplex are the vertices of the old simplex plus the sink+                  numFaces = V.length $ two oldNode++                  --find all the faces of the new node by looking at the faces of the old node+                  testTargets :: Int -> Node -> Vector Node -> Node -> Vector Int -> (Vector Node, Node, Vector Int)+                  testTargets i onode onodes newNode newSinks =+                    let faceVerts =+                          if fstIter then one $ (V.last $ V.init $ result) ! ((two onode) ! i)+                          else one $ (V.last $ result) ! ((two onode) ! i)+                    in+                      if i == numFaces then (onodes, newNode, newSinks)+                      else+                        case V.find (\(_, (v, _, _)) -> V.head v == sink && V.tail v == faceVerts) $ mapWithIndex (\j n -> (j, n)) onodes of+                          Just (j, n) ->+                            testTargets (i + 1) onode+                              (replaceElem j (one n, two n, (thr n) `smartSnoc` newIndex) onodes)+                                (one newNode, (two newNode) `snoc` j, thr newNode) (newSinks |^| (oldSinks ! j))+                          Nothing     -> error "Face not found, HasseDiagram.directedFlagComplex.makeDiagram.makeNode.testTargets"++              in testTargets 0 oldNode nodes (verts, oldIndex `cons` V.empty, V.empty) sinks++            loopSinks :: Int -> Int -> Vector Node -> (Vector Node, Vector Node, Vector (Vector Int), Int)+            loopSinks nodeIndex lastIndex nodes =+              let node     = oldNodes ! nodeIndex+                  sinks    = oldSinks ! nodeIndex+                  numSinks = V.length sinks++                  loop i (modifiedNodes, newNodes, newSinks) =+                    if i == numSinks then (modifiedNodes, newNodes, newSinks, i + lastIndex)+                    else+                      let (modNodes, newNode, ns) = makeNode (i + lastIndex) nodeIndex i modifiedNodes sinks+                      in loop (i + 1) (modNodes, newNodes `snoc` newNode, newSinks `snoc` ns)++              in loop 0 (nodes, V.empty, V.empty)++            loopNodes :: Int -> Int -> Vector Node -> Vector Node -> Vector (Vector Int) -> (Vector Node, Vector Node, Vector (Vector Int))+            loopNodes i lastIndex nodes newNodes newSinks =+              if i == maxindex then (nodes, newNodes, newSinks)+              else+                let (modifiedNodes, nnodes, nsinks, index) = loopSinks i lastIndex nodes+                in loopNodes (i + 1) lastIndex modifiedNodes (newNodes V.++ nnodes) (newSinks V.++ nsinks)++        in loopNodes 0 0 oldNodes V.empty V.empty++      loopLevels :: Int -> HasseDiagram -> Vector Node -> Vector (Vector Int) -> HasseDiagram+      loopLevels iter diagram nextNodes sinks =+        let (modifiedNodes, newNodes, newSinks) = makeLevel (iter < 2) diagram nextNodes sinks+            newDiagram                          = diagram `snoc` modifiedNodes+        in+          if V.null newNodes then newDiagram+          else loopLevels (iter + 1) newDiagram newNodes newSinks++  in loopLevels 0 directedGraph edges fstSinks++-- | Convert a Hasse diagram to a simplicial complex.+toSimplicialComplex :: HasseDiagram -> SimplicialComplex+toSimplicialComplex diagram =+  let sc = V.map (V.map not3) $ V.tail diagram+  in (V.length $ V.head diagram, (V.map (\(v, _) -> (v, V.empty)) $ sc ! 0) `cons` V.tail sc)

@@ -15,6 +15,7 @@  module Matrix   ( BMatrix+  , bMat2String   , getDiagonal   , getUnsignedDiagonal   , transposeMat@@ -26,22 +27,6 @@   , imgInKerBool   ) where -{--FOR DEVS-----------------------------------------------------------------Matrices are transformed by iterating through each row and selecting a pivot. Zero rows are skipped for finding column eschelon form but a row operation is performed (if possible) if there is a zero row for Smith normal form.--To get the smith normal form, the entire pivot row and column is eliminated before continuing. Also, the pivot is always a diagonal element.--To get column eschelon form, every element in the pivot row after the pivot is eliminated. To get the kernel, all column operations to get the matrix to this form are also performed on the identiy matrix. To get the image of one matrix inside the kernel of the one being put into column eschelon form, perform the inverse row operations on the matrix whose image is needed. See stanford paper or the blog post on simplicial homology.--To get the rank of a matrix, look at the number of non-zero columns in the column eschelon form. To get the kernel, look at the columns of the identity (after all of the same column operations have been performed on it) which correspond to zero columns of the column eschelon form.--Eliminating elements is a slighltly more complicated process since only integer operations are allowed. First, every element that must be eliminated is made divisible by the pivot by using the Bezout coefficients from the extended Euclidean algorithm. Once this is done, integer division and subtraction can be used to eliminate the elements.--Boolean matrices are much easier to work with, they are regular matrices with elements modulo 2. Bool is an instance of Num here and the instance is given in Util.----}- import Util  import Data.List as L@@ -56,6 +41,30 @@ -- | Matrix of integers modulo 2. Alternatively, matrix over the field wit h2 elements. type BMatrix = Vector (Vector Bool) +-- | Display an integer matrix.+iMat2String :: IMatrix -> String+iMat2String mat =+  let printVec vec =+        if V.null vec then ""+        else+          let x = V.head vec+          in (show x) L.++ ((if x < 0 then " " else "  ") L.++ (printVec $ V.tail vec))+      print m =+        if V.null m then ""+        else (printVec $ V.head m) L.++ ('\n':(print $ V.tail m))+  in print mat++-- | Display a boolean matrix (as 1's and 0's).+bMat2String :: BMatrix -> String+bMat2String mat =+  let printVec vec =+        if V.null vec then ""+        else (if V.head vec then "1 " else "0 ") L.++ (printVec $ V.tail vec)+      print m =+        if V.null m then ""+        else (printVec $ V.head m) L.++ ('\n':(print $ V.tail m))+  in print mat+ -- | Take the transpose a matrix (no fancy optimizations, yet). transposeMat :: Vector (Vector a) -> Vector (Vector a) transposeMat mat =@@ -148,15 +157,14 @@ --returns Nothing if the entire row is zero chooseGaussPivotBool :: (Int, Int) -> BMatrix -> Maybe (Bool, BMatrix, Maybe (Int, Int)) chooseGaussPivotBool (rowIndex, colIndex) mat =-  let row     = mat ! rowIndex --the following variable should be useful for quickly determining whether or not there are more entries to eleiminate-      indices = V.filter (\index -> index > colIndex) $ V.findIndices id row --but that method is not working for some reason+  let row     = mat ! rowIndex+      indices = V.filter (\index -> index > colIndex) $ V.findIndices id row   in     if not $ row ! colIndex then-      case indices of-        v | V.null v -> Nothing-        v            ->-          let j = V.head v-          in Just (V.length v > 0, V.map (switchElems colIndex j) mat, Just (colIndex, j))+      if V.null indices then Nothing+      else+        let j = V.head indices+        in Just (V.length indices > 0, V.map (switchElems colIndex j) mat, Just (colIndex, j))     else Just (V.length indices > 0, mat, Nothing)  --eliminates pivot row of a boolean matrix

@@ -2,9 +2,9 @@ -- documentation, see http://haskell.org/cabal/users-guide/  name:                Persistence-version:             1.0+version:             1.1 synopsis:            Quickly detect clusters and holes in data.-description:         Persistence is a topological data analysis library motivated by flexibility when it comes to the type of data being analyzed. If you have data that comes with a meaningful function into something of the Ord typeclass, you can use Persistence to detect clusters and holes in the data. You can also use the library to analyze undirected graphs, and interesting features for directed graphs will be added soon.+description:         Persistence is a topological data analysis library motivated by flexibility when it comes to the type of data being analyzed. If you have data that comes with a meaningful function into something of the Ord typeclass, you can use Persistence to detect clusters and holes in the data. You can also use the library to analyze the topology of directed\/undirected weighted\/unweighted graphs, and compare topologies of different data sets. license:             BSD3 license-file:        LICENSE author:              Eben Cowley@@ -16,7 +16,7 @@ cabal-version:       >=1.10  library-  exposed-modules:     SimplicialComplex, Persistence+  exposed-modules:     SimplicialComplex, HasseDiagram, Persistence   other-modules:       Util, Matrix   -- other-extensions:       build-depends:       base >=4.10 && <4.11, vector >=0.12 && <0.13, parallel >=3.2 && <3.3, containers >=0.5 && <0.6, maximal-cliques >=0.1 && <0.2

@@ -17,11 +17,14 @@  Persistent homology is the main event of topological data analysis. It allows one to identify clusters, holes, and voids that persist in the data throughout many scales. The output of the persistence algorithm is a barcode diagram, which currently have a somewhat crude representation in this library. A single barcode represents the filtration index where a feature appears and the index where it disappears (if it does). Thus, short barcodes are typically interpretted as noise and long barcodes are interpretted as actual features. +After you've got the persistent homology of a data set, you might want to compare it with that of a different data set. That's why this release includes two versions of "bottleneck distance," one works only if the number of features in each data set is the same and the other works regardless.+ -}  module Persistence   ( Filtration   , BarCode+  , Extended   , sim2String   , filtr2String   , getComplex@@ -32,6 +35,10 @@   , filterByWeightsLight   , makeVRFiltrationLight   , persistentHomology+  , bottelNeckDistance+  , bottelNeckDistances+  , safeBottelNeckDistance+  , safeBottelNeckDistances   ) where  import Util@@ -56,6 +63,27 @@ -- | (i, Just j) is a feature that appears at filtration index i and disappears at index j, (i, Nothing) begins at i and doesn't disappear. type BarCode = (Int, Maybe Int) +-- | Type for representing inifinite bottleneck distance.+data Extended a = Finite a | Infinity deriving Eq++instance (Ord a, Eq a) => Ord (Extended a) where+  Infinity > Infinity  = False+  Infinity > Finite _  = True+  Finite a > Finite b  = a > b+  Finite _ > Infinity  = False+  Infinity >= Infinity = True+  Infinity >= Finite _ = True+  Finite _ >= Infinity = False+  Finite a >= Finite b = a >= b+  Infinity < Infinity  = False+  Infinity < Finite a  = False+  Finite _ < Infinity  = True+  Finite a < Finite b  = a < b+  Infinity <= Infinity = True+  Infinity <= Finite _ = False+  Finite _ <= Infinity = True+  Finite a <= Finite b = a <= b+ -- | Shows all the information in a simplex. sim2String :: (Int, Vector Int, Vector Int) -> String sim2String (index, vertices, faces) =@@ -90,7 +118,7 @@ {- |   Given a list of scales, a simplicial complex, and a weighted graph (see SimplicialComplex) which encodes a metric on the vertices,   this function creates a filtration out of a simplicial complex by removing simplices that contain edges that are too long for each scale in the list.-  Really, this is a helper function to be called by makeVRFiltrationFast, but I decided to expose it in case you have a simplicial complex and weighted graph lying around.+  This is really a helper function to be called by makeVRFiltrationFast, but I decided to expose it in case you have a simplicial complex and weighted graph lying around.   The scales MUST be in decreasing order for this function. -} filterByWeightsFast :: Ord a => [a] -> (SimplicialComplex, Graph a) -> Filtration@@ -212,7 +240,7 @@           let maxindex = V.head d               begin    = one $ allSimplices ! (dim - 1) ! maxindex           in makeBarCodesAndMark dim (index + 1) marked (replaceElem maxindex (Just d) reduced) (V.tail simplices)-              (if begin == i then codes else (begin, Just i):codes, newMarked)+              ((begin, Just i):codes, newMarked)         where (i, v, f) = V.head simplices               d         = removePivotRows reduced $ removeUnmarked marked f @@ -222,7 +250,7 @@         | V.null d     =           makeEdgeCodes (index + 1) reduced (V.tail edges) (codes, marked `snoc` index)         | otherwise    =-          makeEdgeCodes (index + 1) (replaceElem (V.head d) (Just d) reduced) (V.tail edges) (if i == 0 then codes else (0, Just i):codes, marked)+          makeEdgeCodes (index + 1) (replaceElem (V.head d) (Just d) reduced) (V.tail edges) ((0, Just i):codes, marked)         where (i, v, f) = V.head edges               d         = removePivotRows reduced f @@ -251,3 +279,54 @@       verts = 0 `range` (numVerts - 1)    in makeInfiniteBarCodes $ makeFiniteBarCodes 1 (V.length allSimplices) [fstCodes] (verts `cons` (fstMarked `cons` V.empty)) (fstSlots `cons` V.empty)++{- |+  Return the maximum of minimum distances bewteen the bar codes.+  It's important to note that the function isn't "unsafe" in the sense that it will throw an exception,+  it will just give you a distance regardless of whether or not there is the same number of barcodes is in each list.+-}+bottelNeckDistance :: [BarCode] -> [BarCode] -> Extended Double+bottelNeckDistance diagram1 diagram2 =+  let v1 = V.fromList diagram1+      v2 = V.fromList diagram2++      metric (x1, Just y1) (x2, Nothing) = Infinity+      metric (x1, Nothing) (x2, Just y1) = Infinity+      metric (x1, Just y1) (x2, Just y2) =+        let dx = fromIntegral $ x2 - x1; dy = fromIntegral $ y2 - y1+        in Finite $ sqrt $ dx*dx + dy*dy++  in foldRelation (<) $ V.map (\p -> foldRelation (>) $ V.map (metric p) v2) v1++-- |  Get's all the bottle neck distances; a good way to determine the similarity of the topology of two filtrations.+bottelNeckDistances :: [[BarCode]] -> [[BarCode]] -> [Extended Double]+bottelNeckDistances diagrams1 diagrams2 =+  let d = (L.length diagrams1) - (L.length diagrams2)+  in+    if d >= 0 then (L.zipWith bottelNeckDistance diagrams1 diagrams2) L.++ (L.replicate d Infinity)+    else (L.zipWith bottelNeckDistance diagrams1 diagrams2) L.++ (L.replicate (-d) Infinity)+++-- | If the number of barcodes is the same, return the maximum of minimum distances bewteen the bar codes. Otherwise return nothing.+safeBottelNeckDistance :: [BarCode] -> [BarCode] -> Maybe (Extended Double)+safeBottelNeckDistance diagram1 diagram2 =+  if L.length diagram1 /= L.length diagram2 then Nothing+  else+    let v1 = V.fromList diagram1+        v2 = V.fromList diagram2++        metric (x1, Just y1) (x2, Nothing) = Infinity+        metric (x1, Nothing) (x2, Just y1) = Infinity+        metric (x1, Just y1) (x2, Just y2) =+          let dx = fromIntegral $ x2 - x1; dy = fromIntegral $ y2 - y1+          in Finite $ sqrt $ dx*dx + dy*dy++    in Just $ foldRelation (<) $ V.map (\p -> foldRelation (>) $ V.map (metric p) v2) v1++-- |  Safely get all the bottle neck distances; a good way to determine the similarity of the topology of two filtrations.+safeBottelNeckDistances :: [[BarCode]] -> [[BarCode]] -> [Maybe (Extended Double)]+safeBottelNeckDistances diagrams1 diagrams2 =+  let d = (L.length diagrams1) - (L.length diagrams2)+  in+    if d >= 0 then (L.zipWith safeBottelNeckDistance diagrams1 diagrams2) L.++ (L.replicate d Nothing)+    else (L.zipWith safeBottelNeckDistance diagrams1 diagrams2) L.++ (L.replicate (-d) Nothing)

@@ -9,11 +9,11 @@  An important thing to know about this module is the difference between "fast" and "light" functions. Fast functions encode the metric in a 2D vector, so that distances don't need to be computed over and over again and can instead be accessed in constant time. Unfortunately, this takes O(n^2) memory so I would strongly recomend against using it for larger data sets; this is why "light" functions exist. -A neighborhood graph is a graph where an edge exists between two vertices if and only if the vertices fall within a given distance of each other. Graphs here are more of a helper data type for constructing filtrations, which is why they have a somewhat odd representation.+A neighborhood graph is a graph where an edge exists between two vertices if and only if the vertices fall within a given distance of each other. Graphs here are more of a helper data type for constructing filtrations, which is why they have a somewhat odd representation. Not to worry though, I've supplied a way of encoding graphs of a more generic representation.  The clique complex of a graph is a simplicial complex whose simplices are the complete subgraphs of the given graph. The Vietoris-Rips complex is the clique complex of the neighborhood graph. -Betti numbers measure the number of n-dimensional holes in a given simplicial complex. The 0th Betti number is the number of connected components, or clusters; the 1st Betti number is the number of loops or tunnels; then 2nd Betti number measures voids or hollow volumes; and if you have high-dimensional data, you might have higher Betti numbers representing the number of high-dimensional holes.+Betti numbers measure the number of n-dimensional holes in a given simplicial complex. The 0th Betti number is the number of connected components, or clusters; the 1st Betti number is the number of loops or tunnels; the 2nd Betti number measures voids or hollow volumes; and if you have high-dimensional data, you might have higher Betti numbers representing the number of high-dimensional holes.  -} @@ -22,6 +22,7 @@   , Graph   , sc2String   , getDim+  , encodeWeightedGraph   , makeNbrhdGraph   , makeCliqueComplex   , makeVRComplexFast@@ -44,13 +45,21 @@ {- |   The first component of the pair is the number of vertices.   The second component is a vector whose entries are vectors of simplices of the same dimension.-  The dimension starts at 1 because there is no need to store individual vertices.-  A simplex is represented as a pair: the vector of its vertices in the original data set from which the complex was constructed,-  and the vector of the indices of the faces in the next lowest dimension. Edges are the exception - they do not store reference to their-  faces because it would be redundant with their vertices.+  Index 0 of the vecor corresponds to dimension 1 because there is no need to store individual vertices.+  A simplex is represented as a pair: the vector of its vertices (represent by their index in the original data set),+  and the vector of the indices of the faces in the next lowest dimension. Edges are the exception to the last part -+  they do not store reference to their faces because it would be redundant with their vertices. -} type SimplicialComplex = (Int, Vector (Vector (Vector Int, Vector Int))) +{- |+  This represents the (symmetric) adjacency matrix of some weighted undirected graph. The type `a` is whatever distance is in your data analysis regime.+  The reason graphs are represented like this is because their main function is too speed up the construction of simplicial complexes and filtrations.+  If the clique complex of this graph were to be constructed, only the adjacencies would matter. But if a filtration is constructed all the distances+  will be required over and over again - this allows them to be accessed in constant time.+-}+type Graph a = Vector (Vector (a, Bool))+ -- | Show all the information in a simplicial complex. sc2String :: SimplicialComplex -> String sc2String (v, allSimplices) =@@ -70,12 +79,14 @@ getDim = L.length . snd  {- |-  This represents the (symmetric) adjacency matrix of some undirected graph. The type `a` is whatever distance is in your data analysis regime.-  The reason graphs are represented like this is because their main function is too speed up the construction of simplicial complexes and filtrations.-  If the clique complex of this graph were to be constructed, only the adjacencies would matter. But if a filtration is constructed all the distances-  will be required over and over again - this allows them to be accessed in constant time.+  Takes the number of vertices and each edge paired with its weight to output the graph encoded as a 2D vector.+  If only you have an unweighted graph, you can still encode your graph by simply letting the type `a` be `()`.+  If you have a weighted graph but there isn't a distance between every vertex, you can use the `Extended` type (essentially the same as Maybe)+  from the Persistence module which is already an instance of Ord. -}-type Graph a = Vector (Vector (a, Bool))+encodeWeightedGraph :: Int -> (Int -> Int -> (a, Bool)) -> Graph a+encodeWeightedGraph numVerts edge =+  mapWithIndex (\i r -> mapWithIndex (\j _ -> edge i j) r) $ V.replicate numVerts (V.replicate numVerts ())  {- |   The first argument is a scale, the second is a metric, and the third is the data.

@@ -176,11 +176,16 @@ -- | Given a relation and two vectors, find all pairs of elements satisfying the relation. findBothElems :: (a -> b -> Bool) -> Vector a -> Vector b -> Vector (a, b) findBothElems rel vector1 vector2 =-  let calc i result =+  let len = V.length vector1+ +      calc i result =         let a = vector1 ! i-        in case V.find (\b -> rel a b) vector2 of-          Just b  -> calc (i + 1) $ result `snoc` (a, b)-          Nothing -> calc (i + 1) result+        in+          if i == len then result+          else case V.find (\b -> rel a b) vector2 of+            Just b  -> calc (i + 1) $ result `snoc` (a, b)+            Nothing -> calc (i + 1) result+   in calc 0 V.empty  {- |@@ -278,6 +283,13 @@             Just _  -> calc (x `cons` acc) xs             Nothing -> calc acc xs   in calc V.empty vector1++-- | snocs the element if an only if it isn't already in the vector.+smartSnoc :: Eq a => Vector a -> a -> Vector a+smartSnoc v e =+  case V.elemIndex e v of+    Just _  -> v+    Nothing -> v `snoc` e  -- | Returns whether or not there is an element that satisfies the predicate. existsVec :: (a -> Bool) -> Vector a -> Bool