Persistence 1.1.4.2 → 2.0

raw patch · 17 files changed

+3754/−1433 lines, 17 filesdep +Persistencedep ~basedep ~containersdep ~maximal-cliquessetup-changedPVP ok

version bump matches the API change (PVP)

Dependencies added: Persistence

Dependency ranges changed: base, containers, maximal-cliques, vector

API changes (from Hackage documentation)

Files

• ChangeLog.md +0/−41
• HasseDiagram.hs +0/−146
• LICENSE +1/−1
• Matrix.hs +0/−279
• Persistence.cabal +61/−21
• Persistence.hs +0/−334
• README.md +60/−0
• Setup.hs +0/−2
• SimplicialComplex.hs +0/−282
• Util.hs +0/−327
• changelog.md +67/−0
• src/Persistence/Filtration.hs +895/−0
• src/Persistence/HasseDiagram.hs +203/−0
• src/Persistence/Matrix.hs +932/−0
• src/Persistence/SimplicialComplex.hs +608/−0
• src/Persistence/Util.hs +329/−0
• testing/Testing.hs +598/−0

@@ -1,41 +0,0 @@-# Revision history for Persistence--## 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.     --## 1.1.3 -- 2018-07-30--### Changed--- Fixed a major bug with persistent homology; high dimensional holes were being detected in low dimensional data sets.-- Persistent homology now filters out bar codes of the form (i, Just i), as they say nothing about the topology of the underlying complexes.--## 1.1.4 -- 2018-09-15--### Changed--- Fixed spelling error.-- Persistence now exports the constructors for `Extended a`.--## 1.1.4.1 -- 2018-09-15--### Changed--- Fixed all spelling errors, should actually build now.--## 1.1.4.2 -- 2018-09-15--### Changed--- Fixed non-exhaustive pattern match in `BottleNeckDistance` functions.

@@ -1,146 +0,0 @@-{- |-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)

@@ -1,4 +1,4 @@-Copyright (c) 2018, Eben Cowley+Copyright (c) 2019, Eben Kadile  All rights reserved. 

@@ -1,279 +0,0 @@-{- |-Module     : Persistence.Matrix-Copyright  : (c) Eben Cowley, 2018-License    : BSD 3 Clause-Maintainer : eben.cowley42@gmail.com-Stability  : experimental--This module contains a variety of matrix utility functions, used in the computation of Betti numbers and simplicial homology groups.--Most importantly, it includes functions for computing the rank, normal form, and kernel of matrices. For the computation of homology groups and Betti numbers, one must perform column operations on one matrix to get it into column echelon form and find its kernel while also performing the inverse row operations on the next matrix to be operated on.--Bool is an instance of Num here (instance given in Util) so that functions can be somewhat generalized to act on both integers and integers modulo 2.---}--module Matrix-  ( BMatrix-  , bMat2String-  , getDiagonal-  , getUnsignedDiagonal-  , transposeMat-  , transposePar-  , multiply-  , multiplyPar-  , rankBool-  , kernelBool-  , imgInKerBool-  ) where--import Util--import Data.List as L-import Data.Vector as V-import Control.Parallel.Strategies----BASIC STUFF----------------------------------------------------------------- | Matrix of integers.-type IMatrix = Vector (Vector Int)---- | 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 =-  V.map (\i -> V.map (\row -> row ! i) mat) $ 0 `range` ((V.length $ V.head mat) - 1)---- | Take the transpose of a matrix using parallel evaluation of rows.-transposePar :: Vector (Vector a) -> Vector (Vector a)-transposePar mat =-  parMapVec (\i -> V.map (\row -> row ! i) mat) $ 0 `range` ((V.length $ V.head mat) - 1)---- | Multiply two matrices-multiply :: Num a => Vector (Vector a) -> Vector (Vector a) -> Vector (Vector a)-multiply mat1 mat2 =-  let t = transposeMat mat2-  in V.map (\row -> V.map (dotProduct row) t) mat1---- | Multiply matrices, evaluate rows in parallel if processors are available-multiplyPar :: Num a => Vector (Vector a) -> Vector (Vector a) -> Vector (Vector a)-multiplyPar mat1 mat2 = runEval $ do-  let t = transposeMat mat2-  rseq t-  return $ parMapVec (\row -> V.map (dotProduct row) t) mat1---- | Get the diagonal elements.-getDiagonal :: Vector (Vector a) -> [a]-getDiagonal matrix =-  if V.null matrix then []-  else L.map (\i -> matrix ! i ! i) [0..(min (V.length matrix) (V.length $ V.head matrix)) - 1]---- | Get the absolute value of each of the diagonal elements in a list.-getUnsignedDiagonal :: Num a => Vector (Vector a) -> [a]-getUnsignedDiagonal matrix =-  if V.null matrix then []-  else L.map (\i -> abs $ matrix ! i ! i) [0..(min (V.length matrix) (V.length $ V.head matrix)) - 1]----assumes index1 < index2-colOperation :: Int -> Int -> (Int, Int, Int, Int) -> IMatrix -> IMatrix-colOperation index1 index2 (c11, c12, c21, c22) matrix =-  let calc row =-        let elem1  = row ! index1-            elem2  = row ! index2-            first  = V.take index1 row-            second = V.drop (index1 + 1) (V.take index2 row)-            third  = V.drop (index2 + 1) row-        in first V.++ (cons (c11*elem1 + c12*elem2) second) V.++ (cons (c22*elem2 + c21*elem1) third)-  in V.map calc matrix--colOperationPar :: Int -> Int -> (Int, Int, Int, Int) -> IMatrix -> IMatrix-colOperationPar index1 index2 (c11, c12, c21, c22) matrix =-  let calc row =-        let elem1  = row ! index1-            elem2  = row ! index2-            first  = V.take index1 row-            second = V.drop (index1 + 1) (V.take index2 row)-            third  = V.drop (index2 + 1) row-        in first V.++ (cons (c11*elem1 + c12*elem2) second) V.++ (cons (c22*elem2 + c21*elem1) third)-  in parMapVec calc matrix----assumes index1 < index2-rowOperation :: Int -> Int -> (Int, Int, Int, Int) -> IMatrix -> IMatrix-rowOperation index1 index2 (c11, c12, c21, c22) matrix =-  let row1   = matrix ! index1-      row2 = matrix ! index2-      first  = V.take index1 matrix-      second = V.drop (index1 + 1) $ V.take index2 matrix-      third  = V.drop (index2 + 1) matrix-  in first V.++ (cons ((c11 `mul` row1) `add` (c12 `mul` row2)) second)-    V.++ (cons ((c22 `mul` row2) `add` (c21 `mul` row1)) third)--rowOperationPar :: Int -> Int -> (Int, Int, Int, Int) -> IMatrix -> IMatrix-rowOperationPar index1 index2 (c11, c12, c21, c22) matrix =-  let row1   = matrix ! index1-      row2 = matrix ! index2-      first  = V.take index1 matrix-      second = V.drop (index1 + 1) (V.take index2 matrix)-      third  = V.drop (index2 + 1) matrix-  in runEval $ do-     a <- rpar $ (c11 `mul` row1) `add` (c12 `mul` row2)-     b <- rpar $ (c21 `mul` row1) `add` (c22 `mul` row2)-     rseq (a,b)-     return $ first V.++ (a `cons` second) V.++ (b `cons` third)----BOOLEAN MATRICES------------------------------------------------------------RANK------------------------------------------------------------------------given the index of the pivot row and the matrix---determines whether there is a non-zero element in the row, does necessary rearranging---and returns the column operation that was performed if there was one---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-      indices = V.filter (\index -> index > colIndex) $ V.findIndices id row-  in-    if not $ row ! colIndex then-      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-elimRowBool :: (Int, Int) -> Int -> BMatrix -> BMatrix-elimRowBool (rowIndex, colIndex) numCols elems =-  let row = elems ! rowIndex-      elim i mat-        | i == numCols  = mat-        | not $ row ! i = elim (i + 1) mat-        | otherwise     = elim (i + 1) $ V.map (\row -> replaceElem i ((row ! i) + (row ! colIndex)) row) mat-  in elim (colIndex + 1) elems---- | Find the rank of a mod 2 matrix (number of linearly independent columns).-rankBool :: BMatrix -> Int-rankBool matrix =-  let rows  = V.length matrix-      cols  = V.length $ V.head matrix-      cols1 = cols - 1--      doColOps (rowIndex, colIndex) mat =-        if rowIndex == rows || colIndex == cols then mat else-          case chooseGaussPivotBool (rowIndex, colIndex) mat of-            Just (True, mx, _)  -> doColOps (rowIndex + 1, colIndex + 1) $ elimRowBool (rowIndex, colIndex) cols mx-            Just (False, mx, _) -> doColOps (rowIndex + 1, colIndex + 1) mat-            Nothing             -> doColOps (rowIndex + 1, colIndex) mat--      countNonZeroCols mat =-        V.sum $ V.map (\i ->-           if existsVec (\j -> mat ! j ! i /= 0) (0 `range` (rows - 1)) then 1 else 0) $ 0 `range` cols1-  in countNonZeroCols $ doColOps (0, 0) matrix----KERNEL----------------------------------------------------------------------eliminates all the entries in the pivot row that come after the pivot, after the matrix has been improved---returns the new matrix (fst) paired with the identity with whatever column operations were performed (snd)-elimRowBoolWithId :: (Int, Int) -> Int -> BMatrix -> BMatrix -> (BMatrix, BMatrix)-elimRowBoolWithId (rowIndex, colIndex) numCols elems identity =-  let row = elems ! rowIndex-      elim i mat ide-        | i == numCols  = (mat, ide)-        | not $ row ! i = elim (i + 1) mat ide-        | otherwise     =-          let transform = V.map (\row -> replaceElem i ((row ! i) + (row ! colIndex)) row)-          in elim (i + 1) (transform mat) (transform ide)-  in elim (colIndex + 1) elems identity---- | Finds the basis of the kernel of a matrix, arranges the basis vectors into the rows of a matrix.-kernelBool :: BMatrix -> BMatrix-kernelBool matrix =-  let rows     = V.length matrix-      cols     = V.length $ V.head matrix-      cols1    = cols - 1-      identity = V.map (\i -> (V.replicate i False) V.++ (cons True (V.replicate (cols1 - i) False))) $ 0 `range` cols1--      doColOps (rowIndex, colIndex) (ker, ide) =-        if rowIndex == rows || colIndex == cols then (ker, ide)-        else-          case chooseGaussPivotBool (rowIndex, colIndex) ker of-            Just (True, _, Nothing)      ->-              doColOps (rowIndex + 1, colIndex + 1) $-                elimRowBoolWithId (rowIndex, colIndex) cols ker ide-            Just (True, mx, Just (i, j)) ->-              doColOps (rowIndex + 1, colIndex + 1) $-                elimRowBoolWithId (rowIndex, colIndex) cols mx $ V.map (switchElems i j) ide-            Just (False, _, Just (i, j)) -> doColOps (rowIndex + 1, colIndex + 1) (ker, V.map (switchElems i j) ide)-            Just (False, _, _)           -> doColOps (rowIndex + 1, colIndex + 1) (ker, ide)-            Nothing                      -> doColOps (rowIndex + 1, colIndex) (ker, ide)--      result = doColOps (0, 0) (matrix, identity)-      ker    = fst result-      img    = snd result-  in V.map (\i -> img ! i) $ V.filter (\i -> forallVec (\row -> not $ row ! i) ker) $ 0 `range` cols1----IMAGE IN BASIS OF KERNEL--------------------------------------------------elimRowBoolWithInv :: (Int, Int) -> Int -> BMatrix -> BMatrix -> (BMatrix, BMatrix)-elimRowBoolWithInv (rowIndex, colIndex) numCols toColEch toImage =-  let row = toColEch ! rowIndex-      elim i ech img-          | i == numCols  = (ech, img)-          | not $ row ! i = elim (i + 1) ech img-          | otherwise     =-            let transform1 = V.map (\r -> replaceElem i ((r ! i) + (r ! colIndex)) r)-                transform2 = \mat -> replaceElem colIndex ((mat ! i) `add` (mat ! colIndex)) mat-            in elim (i + 1) (transform1 ech) (transform2 img)-  in elim (colIndex + 1) toColEch toImage---- | Calculates the image of the second matrix represented in the basis of the kernel of the first matrix.-imgInKerBool :: BMatrix -> BMatrix -> BMatrix-imgInKerBool toColEch toImage =-  let rows  = V.length toColEch-      cols  = V.length $ V.head toColEch-      cols1 = cols - 1--      doColOps (rowIndex, colIndex) (ech, img) =-        if rowIndex == rows || colIndex == cols then (ech, img)-        else-          case chooseGaussPivotBool (rowIndex, colIndex) ech of-            Just (True, _, Nothing)       ->-              doColOps (rowIndex + 1, colIndex + 1) $-                elimRowBoolWithInv (rowIndex, colIndex) cols ech img-            Just (True, mx, Just (i, j))  ->-              doColOps (rowIndex + 1, colIndex + 1) $-                elimRowBoolWithInv (rowIndex, colIndex) cols mx $ switchElems i j img-            Just (False, mx, Just (i, j)) -> doColOps (rowIndex + 1, colIndex + 1) (mx, switchElems i j img)-            Just (False, _, _)            -> doColOps (rowIndex + 1, colIndex + 1) (ech, img)-            Nothing                       -> doColOps (rowIndex + 1, colIndex) (ech, img)--      result = doColOps (0, 0) (toColEch, toImage)-      ker    = fst result-      img    = snd result-  in V.map (\i -> img ! i) $ V.filter (\i -> forallVec (\row -> not $ row ! i) ker) $ 0 `range` cols1

@@ -1,24 +1,64 @@--- Initial Persistence.cabal generated by cabal init.  For further --- documentation, see http://haskell.org/cabal/users-guide/+cabal-version: 1.12 -name:                Persistence-version:             1.1.4.2-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 the topology of directed\/undirected weighted\/unweighted graphs, and compare topologies of different data sets.-license:             BSD3-license-file:        LICENSE-author:              Eben Cowley-maintainer:          eben.cowley42@gmail.com--- copyright:           -category:            Data-build-type:          Simple-extra-source-files:  ChangeLog.md-cabal-version:       >=1.10+-- This file has been generated from package.yaml by hpack version 0.31.2.+--+-- see: https://github.com/sol/hpack+--+-- hash: 8b7fcc7836b916f046653da257033f5fbdeb90a49069ff85be71bca358f359c0 +name:           Persistence+version:        2.0+category:       Data, Math+synopsis:       A versatile library for topological data anlysis.+description:    A topological data anlysis library motivated by flexibility when it comes to the type of data being analyzed. If your data comes with a meaningful binary function into into an ordered set, you can use Persistence to analyze your data. The library also provides functions for analyzing directed\/undirected, weighted\/unweighted graphs. See the README for resources on learning about topological data anlysis.+homepage:       https://github.com/Ebanflo42/Persistence+bug-reports:    https://github.com/Ebanflo42/Persistence/issues+author:         Eben Kadile+maintainer:     eben.cowley42@gmail.com+copyright:      2019 Eben Kadile+license:        BSD3+license-file:   LICENSE+build-type:     Simple+extra-source-files:+    README.md+    changelog.md++source-repository head+  type: git+  location: https://github.com/Ebanflo42/Persistence+ library-  exposed-modules:     SimplicialComplex, HasseDiagram, Persistence-  other-modules:       Util, Matrix-  -- other-extensions:    -  build-depends:       base >=4.9 && <4.12, vector >=0.12 && <0.13, parallel >=3.2 && <3.3, containers >=0.5 && <0.6, maximal-cliques >=0.1 && <0.2-  -- hs-source-dirs:      -  default-language:    Haskell2010+  exposed-modules:+      Persistence.Filtration+      Persistence.HasseDiagram+      Persistence.SimplicialComplex+  other-modules:+      Paths_Persistence+      Persistence.Matrix+      Persistence.Util+  hs-source-dirs:+      src+  build-depends:+      base >=4.0 && <4.13+    , containers >0.5+    , maximal-cliques >=0.1+    , parallel ==3.2.*+    , vector >=0.12+  default-language: Haskell2010++test-suite persistence-test+  type: exitcode-stdio-1.0+  main-is: Testing.hs+  other-modules:+      Paths_Persistence+  hs-source-dirs:+      testing+  ghc-options: -threaded -rtsopts -with-rtsopts=-N+  build-depends:+      Persistence+    , base >=4.0 && <4.11.2+    , containers ==0.5.*+    , maximal-cliques >=0.1+    , parallel ==3.2.*+    , vector >=0.12+  default-language: Haskell2010

@@ -1,334 +0,0 @@-{- |-Module     : Persistence.Persistence-Copyright  : (c) Eben Cowley, 2018-License    : BSD 3 Clause-Maintainer : eben.cowley42@gmail.com-Stability  : experimental--This module contains functions for constructing filtrations and computing persistent homology, as well as a few utility functions for working with filtrations.--A filtration is a finite sequence of simplicial complexes where each complex is a subset of the next. This means that a filtration can be thought of as a single simplicial complex where each of the simplices is labeled with a "filtration index" that represents the index in the sequence where that simplex enters the filtration.--One way to create a filtration, given a simplicial complex, a metric for the vertices, and a list of distances, is to loop through the distances from greatest to least: create a simplicial complex each iteration which excludes simplices that contain pairs of vertices which are further than the current distance apart. This method will produce a filtration of Vietoris-Rips complexes - each filtration index will correspond to a VR complex whose scale is the corresponding distance.--An essential thing to note about the way this library is set up is the distinction between "fast" and "light" functions. Light functions call the metric every time distance between two points is required, which is a lot. Fast functions store the distances between points and access them in constant time, BUT this means they use O(n^2) memory with respect to the number of data points, so it's a really bad idea to use this optimization on substantially large data.--IMPORTANT NOTE: This library currently only handles filtrations where all the vertices have filtration index 0. Since I haven't thought of any way filtrations that don't respect this property could come up in applications, I have chosen to exclude them in favor of less memory usage and a slight speedup.--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 (Finite, Infinity)-  , sim2String-  , filtr2String-  , getComplex-  , getNumSimplices-  , getDimension-  , filterByWeightsFast-  , makeVRFiltrationFast-  , filterByWeightsLight-  , makeVRFiltrationLight-  , persistentHomology-  , bottleNeckDistance-  , bottleNeckDistances-  , safeBottleNeckDistance-  , safeBottleNeckDistances-  ) where--import Util-import Matrix-import SimplicialComplex--import Data.List as L-import Data.Vector as V-import Control.Parallel.Strategies-import Data.Algorithm.MaximalCliques----DATA TYPES----------------------------------------------------------------{- |-  The first component is simply the number of vertices (all vertices are assumed to have filtration index 0).-  The second component is a vector with an entry for each dimension of simplices, starting at dimension 1 for edges.-  Each simplex is represented as a triple: its filtration index, the indices of its vertices in the original data, and the indices of its faces in the next lowest dimension.-  Edges do not have reference to their faces, as it would be redundant with their vertices. All simplices are sorted according to filtration index upon construction of the filtration.--}-type Filtration = (Int, Vector (Vector (Int, Vector Int, Vector Int)))---- | (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) =-  "Filtration index: " L.++ (show index) L.++-    "; Vertex indices: " L.++ (show vertices) L.++-      "; Boundary indices: " L.++ (show faces) L.++ "\n"---- | Shows all the information in a filtration.-filtr2String :: Filtration -> String-filtr2String = (intercalate "\n") . toList . (V.map (L.concat . toList . (V.map sim2String))) . snd---- | Gets the simplicial complex specified by the filtration index. This is O(n) with respect to the number of simplices.-getComplex :: Int -> Filtration -> SimplicialComplex-getComplex index (n, simplices) = (n, V.map (V.map not1 . V.filter (\(i, _, _) -> i == index)) simplices)--{- |-  The first argument is a list of dimensions, the second argument is a list of filtration indices.-  The function returns the number of simplices in the filtration whose dimension and index exist in the respective lists.--}-getNumSimplices :: [Int] -> [Int] -> Filtration -> Int-getNumSimplices dimensions indices (_, simplices) =-  L.length $ L.concat $-    L.map (\d -> V.toList $ V.filter (\(i, _, _) ->-      L.elemIndex i indices /= Nothing) $ simplices ! (d - 1)) dimensions---- | Return the dimension of the highest dimensional simplex in the filtration (constant time).-getDimension :: Filtration -> Int-getDimension = V.length . snd----FILTRATION CONSTRUCTION---------------------------------------------------{- |-  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.-  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-filterByWeightsFast scales ((numVerts, simplices), graph) =-  let edgeInSimplex edge simplex = (existsVec (\x -> V.head edge == x) simplex) && (existsVec (\x -> V.last edge == x) simplex)-      edgeTooLong scale edge     = scale <= (fst $ graph ! (edge ! 0) ! (edge ! 1))-      maxIndex                   = (L.length scales) - 1--      calcIndices 0 [] sc         = sc-      calcIndices i (scl:scls) sc =-        let longEdges = V.filter (edgeTooLong scl) $ V.map (\(i, v, f) -> v) $ V.head sc --find edges excluded by this scale-        in calcIndices (i - 1) scls $ V.map (V.map (\(j, v, f) ->-          if j == 0 then --if the simplex has not yet been assigned a fitration index-            if existsVec (\edge -> edgeInSimplex edge v) longEdges then (i, v, f) --if a long edge is in the simplex, assign it the current index-            else (0, v, f) --otherwise wait until next iteration-          else (j, v, f))) sc --otherwise leave it alone--      sortFiltration simplices =-        let sortedSimplices =-              V.map (quicksort (\((i, _, _), _) ((j, _, _), _) -> i > j)) $-                V.map (mapWithIndex (\i s -> (s, i))) simplices-            newFaces dim (i, v, f) =-              let findNew j =-                    case V.findIndex (\x -> snd x == j) $ sortedSimplices ! (dim - 1) of-                      Just k  -> k-                      Nothing -> error "Persistence.sortFiltration.newFaces.findNew"-              in (i, v, (V.map findNew f))-        in-          if V.null simplices then simplices-          else mapWithIndex (\i ss -> V.map ((newFaces i) . fst) ss) sortedSimplices--  in (numVerts, sortFiltration $ --sort the simplices by filtration index-      calcIndices maxIndex (L.tail scales) $-        V.map (V.map (\(v, f) -> (0, v, f))) $ simplices)--{- |-  Given a list of scales, a metric, and a data set, this function constructs a filtration of the Vietoris-Rips complexes associated with the scales.-  The scales MUST be in decreasing order. Note that this a fast function, meaning it uses O(n^2) memory to quickly access distances where n is the number of data points.--}-makeVRFiltrationFast :: (Ord a, Eq b) => [a] -> (b -> b -> a) -> [b] -> Filtration-makeVRFiltrationFast scales metric dataSet = filterByWeightsFast scales $ makeVRComplexFast (L.head scales) metric dataSet---- | The same as filterbyWeightsFast except it uses far less memory at the cost of speed. Note that the scales must be in decreasing order.-filterByWeightsLight :: Ord a => [a] -> (b -> b -> a) -> [b] -> SimplicialComplex -> Filtration-filterByWeightsLight scales metric dataSet (numVerts, simplices) =-  let edgeInSimplex edge simplex = (existsVec (\x -> V.head edge == x) simplex) && (existsVec (\x -> V.last edge == x) simplex)-      vector                     = V.fromList dataSet-      edgeTooLong scale edge     = scale <= (metric (vector ! (edge ! 0)) (vector ! (edge ! 1)))-      maxIndex                   = (L.length scales) - 1--      calcIndices 0 [] sc         = sc-      calcIndices i (scl:scls) sc =-        let longEdges = V.filter (edgeTooLong scl) $ V.map (\(i, v, f) -> v) $ V.head sc --find edges excluded by this scale-        in calcIndices (i - 1) scls $ V.map (V.map (\(j, v, f) ->-          if j == 0 then --if the simplex has not yet been assigned a fitration index-            if existsVec (\edge -> edgeInSimplex edge v) longEdges then (i, v, f) --if a long edge is in the simplex, assign it the current index-            else (0, v, f) --otherwise wait until next iteration-          else (j, v, f))) sc --otherwise leave it alone--      sortFiltration simplices =-        let sortedSimplices =-              V.map (quicksort (\((i, _, _), _) ((j, _, _), _) -> i > j)) $-                V.map (mapWithIndex (\i s -> (s, i))) simplices-            newFaces dim (i, v, f) =-              let findNew j =-                    case V.findIndex (\x -> snd x == j) $ sortedSimplices ! (dim - 1) of-                      Just k  -> k-                      Nothing -> error "Persistence.filterByWeightsLight.sortFiltration.newFaces.findNew"-              in (i, v, (V.map findNew f))-        in-          if V.null simplices then simplices-          else mapWithIndex (\i ss -> V.map ((newFaces i) . fst) ss) sortedSimplices--  in (numVerts, sortFiltration $ --sort the simplices by filtration index-      calcIndices maxIndex (L.tail scales) $-        V.map (V.map (\(v, f) -> (0, v, f))) $ simplices)---- | Given a list of scales in decreasing order, a metric, and a data set, this constructs the filtration of Vietoris-Rips complexes corresponding to the scales.-makeVRFiltrationLight :: (Ord a, Eq b) => [a] -> (b -> b -> a) -> [b] -> Filtration-makeVRFiltrationLight scales metric dataSet = filterByWeightsLight scales metric dataSet $ makeVRComplexLight (L.head scales) metric dataSet----PERSISTENT HOMOLOGY-------------------------------------------------------{- |-  Each index in the list is a list of barcodes whose dimension corresponds to the index.-  That is, the first list will represent clusters, the second list will represent tunnels or punctures,-  the third will represent hollow volumes, and the nth index list will represent n-dimensional holes in the data.--}-persistentHomology :: Filtration -> [[BarCode]]-persistentHomology (numVerts, allSimplices) =-  let-      --union minus intersection-      uin :: Ord a => Vector a -> Vector a -> Vector a-      u `uin` v =-        let findAndInsert elem vec =-              let (vec1, vec2) = biFilter (\x -> x > elem) vec-              in-                if V.null vec2 then vec `snoc` elem-                else if elem == V.head vec2 then vec1 V.++ V.tail vec2-                else vec1 V.++ (elem `cons` vec2)-        in-          if V.null u then v-          else (V.tail u) `uin` (findAndInsert (V.head u) v)-          -      removeUnmarked marked = V.filter (\x -> existsVec (\y -> y == x) marked)--      removePivotRows reduced chain =-        if V.null chain then chain-        else-          case reduced ! (V.head chain) of-            Nothing -> chain-            Just t  -> removePivotRows reduced (chain `uin` t) --eliminate the element corresponding to the pivot in a different chain--      makeBarCodesAndMark :: Int -> Int -> Vector Int -> Vector (Maybe (Vector Int)) -> Vector (Int, Vector Int, Vector Int) -> ([BarCode], Vector Int) -> ([BarCode], Vector Int, Vector Int)-      makeBarCodesAndMark dim index marked reduced simplices (codes, newMarked)-        | V.null simplices = (codes, newMarked, V.findIndices (\x -> x == Nothing) reduced)-        | V.null d         = makeBarCodesAndMark dim (index + 1) marked reduced (V.tail simplices) (codes, newMarked `snoc` index)-        | otherwise        =-          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)-              ((begin, Just i):codes, newMarked)-        where (i, v, f) = V.head simplices-              d         = removePivotRows reduced $ removeUnmarked marked f--      makeEdgeCodes :: Int -> Vector (Maybe (Vector Int)) -> Vector (Int, Vector Int, Vector Int) -> ([BarCode], Vector Int) -> ([BarCode], Vector Int, Vector Int)-      makeEdgeCodes index reduced edges (codes, marked)-        | V.null edges = (codes, marked, V.findIndices (\x -> x == Nothing) reduced)-        | 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) ((0, Just i):codes, marked)-        where (i, v, f) = V.head edges-              d         = removePivotRows reduced f--      makeFiniteBarCodes :: Int -> Int -> [[BarCode]] -> Vector (Vector Int) -> Vector (Vector Int) -> ([[BarCode]], Vector (Vector Int), Vector (Vector Int))-      makeFiniteBarCodes dim maxdim barcodes marked slots =-        if dim == maxdim then (barcodes, marked, slots)-        else-          let (newCodes, newMarked, unusedSlots) = makeBarCodesAndMark dim 0 (V.last marked) (V.replicate (V.length $ allSimplices ! (dim - 1)) Nothing) (allSimplices ! dim) ([], V.empty)-          in makeFiniteBarCodes (dim + 1) maxdim (barcodes L.++ [newCodes]) (marked `snoc` newMarked) (slots `snoc` unusedSlots)--      makeInfiniteBarCodes :: ([[BarCode]], Vector (Vector Int), Vector (Vector Int)) -> [[BarCode]]-      makeInfiniteBarCodes (barcodes, marked, unusedSlots) =-        let makeCodes i codes =-              let slots = unusedSlots ! i; marks = marked ! i-              in codes L.++ (V.toList $ V.map (\j -> (one $ allSimplices ! (i - 1) ! j, Nothing)) $ slots |^| marks)-            loop _ []     = []-            loop 0 (x:xs) = (x L.++ (V.toList $ V.map (\j -> (0, Nothing)) $ (unusedSlots ! 0) |^| (marked ! 0))):(loop 1 xs)-            loop i (x:xs) = (makeCodes i x):(loop (i + 1) xs)-        in loop 0 barcodes--      edges    = V.map (\(i, v, f) -> (i, v, (V.reverse v))) $ V.head allSimplices-      numEdges = V.length edges--      (fstCodes, fstMarked, fstSlots) = makeEdgeCodes 0 (V.replicate numVerts Nothing) edges ([], V.empty)--      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.--}-bottleNeckDistance :: [BarCode] -> [BarCode] -> Extended Double-bottleNeckDistance diagram1 diagram2 =-  let v1 = V.fromList diagram1-      v2 = V.fromList diagram2--      metric (x1, Nothing) (x2, Nothing) = Finite $ fromIntegral $ abs $ x1 - x2-      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.-bottleNeckDistances :: [[BarCode]] -> [[BarCode]] -> [Extended Double]-bottleNeckDistances diagrams1 diagrams2 =-  let d = (L.length diagrams1) - (L.length diagrams2)-  in-    if d >= 0 then (L.zipWith bottleNeckDistance diagrams1 diagrams2) L.++ (L.replicate d Infinity)-    else (L.zipWith bottleNeckDistance 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.-safeBottleNeckDistance :: [BarCode] -> [BarCode] -> Maybe (Extended Double)-safeBottleNeckDistance diagram1 diagram2 =-  if L.length diagram1 /= L.length diagram2 then Nothing-  else-    let v1 = V.fromList diagram1-        v2 = V.fromList diagram2--        metric (x1, Nothing) (x2, Nothing) = Finite $ fromIntegral $ abs $ x1 - x2-        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.-safeBottleNeckDistances :: [[BarCode]] -> [[BarCode]] -> [Maybe (Extended Double)]-safeBottleNeckDistances diagrams1 diagrams2 =-  let d = (L.length diagrams1) - (L.length diagrams2)-  in-    if d >= 0 then (L.zipWith safeBottleNeckDistance diagrams1 diagrams2) L.++ (L.replicate d Nothing)-    else (L.zipWith safeBottleNeckDistance diagrams1 diagrams2) L.++ (L.replicate (-d) Nothing)

@@ -0,0 +1,60 @@+# Persistence+A topological data analysis library for Haskell.++This library is motivated by flexibility when it comes to the type of data being analyzed. If your data comes with a meaningful binary function into into an ordered set, you can use Persistence to analyze your data. The library also provides functions for analyzing directed\/undirected, weighted\/unweighted graphs.++Visit https://hackage.haskell.org/package/Persistence to see the documentation for the stable version. There is also documentation in each module.++GitHub: https://github.com/Ebanflo42/Persistence++If you have the Haskell `stack` tool installed, compile and run tests with `stack test` (you may have to expose the modules `Util` and `Matrix` in the file `Persistence.cabal` to get it to work).++# Learning about Topological Data Analysis++Computing simplicial homology:++https://jeremykun.com/2013/04/10/computing-homology/++Constructing the Vietoris-Rips complex:++https://pdfs.semanticscholar.org/e503/c24dcc7a8110a001ae653913ccd064c1044b.pdf++Constructing the Cech complex:++https://www.academia.edu/15228439/Efficient_construction_of_the_%C4%8Cech_complex++Computing persistent homology:++http://geometry.stanford.edu/papers/zc-cph-05/zc-cph-05.pdf++The algorithm for finding the directed clique complex is inspired by the pseudocode in the supplementary materials of this paper:++https://www.frontiersin.org/articles/10.3389/fncom.2017.00048/full++Computing and working with persistence landscapes:++https://academic.csuohio.edu/bubenik_p/papers/persistenceLandscapes.pdf++# Major TODOs:++`Testing.hs`:++1) More tests for Persistence landscape functions.++2) Make some filtrations whose vertices don't all have index 0 and test persistent homology on them.++`SimplicialComplex.hs`:++1) Fix simplicial homology over the integers.++2) Implement construction of the Cech complex (n points form an (n-1)-simplex if balls of a certain radius centered at each of the points intersect).++3) Implement construction of the alpha-complex (sub-complex of the Delaunay triangulation where the vertices of every simplex are within a certain distance).++General:++1) Update documentation for `Filtration.hs` once more changes have been finalized.++2) A more consistent, well-motivated, and concise philosophy for parallelism needs to be implemented.++See each of the files for an overview of its inner-workings.

@@ -1,2 +0,0 @@-import Distribution.Simple-main = defaultMain

@@ -1,282 +0,0 @@-{- |-Module     : Persistence.SimplicialComplex-Copyright  : (c) Eben Cowley, 2018-License    : BSD 3 Clause-Maintainer : eben.cowley42@gmail.com-Stability  : experimental--This module provides functions for constructing neighborhood graphs, clique complexes, Vietoris-Rips complexes, as well as the computation of Betti numbers.--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. 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; 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.---}--module SimplicialComplex-  ( SimplicialComplex-  , Graph-  , sc2String-  , getDim-  , encodeWeightedGraph-  , makeNbrhdGraph-  , makeCliqueComplex-  , makeVRComplexFast-  , makeVRComplexLight-  , bettiNumbers-  , bettiNumbersPar-  ) where--import Util-import Matrix--import Data.List as L-import Data.Vector as V-import Data.IntSet as S-import Control.Parallel.Strategies-import Data.Algorithm.MaximalCliques----BASIC STUFF---------------------------------------------------------------{- |-  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.-  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) =-  if V.null allSimplices then (show v) L.++ " vertices"-  else-  let edges     = V.head allSimplices-      simplices = V.tail allSimplices-      showSimplex s =-        '\n':(intercalate "\n" $ V.toList $ V.map show s)-      showAll sc    =-        if V.null sc then '\n':(show v) L.++ " vertices"-        else showSimplex (V.head sc) L.++ ('\n':(showAll (V.tail sc)))-  in (intercalate "\n" $ V.toList $ V.map (show . fst) edges) L.++ ('\n':(showAll simplices))---- | Get the dimension of the highest dimensional simplex (constant time).-getDim :: SimplicialComplex -> Int-getDim = L.length . snd--{- |-  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.--}-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.-  This function records the distance between every element of the data and whether or not it is smaller than the given scale.--}-makeNbrhdGraph :: (Ord a, Eq b) => a -> (b -> b -> a) -> [b] -> Graph a-makeNbrhdGraph scale metric dataSet =-  let vector      = V.fromList dataSet-  in V.map (\x -> V.map (\y -> let d = metric x y in (d, d <= scale)) vector) vector----CONSTRUCTION--------------------------------------------------------------{- |-  Makes a simplicial complex where the simplices are the complete subgraphs (cliques) of the given graph.-  Mainly a helper function for makeVRComplexFast, but it might be useful if you happen to have a graph you want to analyze.-  This utilizes any available processors in parallel because the construction is quite expensive.--}-makeCliqueComplex :: Graph a -> SimplicialComplex-makeCliqueComplex graph =-  let numVerts = V.length graph-      --make a list with an entry for every dimension of simplices-      organizeCliques 1 _       = []-      organizeCliques i cliques =-        let helper = biFilter (\simplex -> i == V.length simplex) cliques --find the simplices with the given number of vertices-        in (fst helper):(organizeCliques (i - 1) $ snd helper) --append them to the next iteration of the function--      makePair simplices = --pair the organized maximal cliques with the dimension of the largest clique-        case simplices of-          (x:_) ->-            let dim = V.length x-            in (dim, organizeCliques dim $ V.fromList simplices)-          []    -> (1, []) --if there are no maximal cliques this acts as a flag so that the algorithm doesn't try to generate all the other simplices--      maxCliques :: (Int, Vector (Vector (Vector Int)))-      maxCliques = --find all maximal cliques and sort them from largest to smallest (excludes maximal cliques which are single points)-        (\(x, y) -> (x, V.fromList y)) $ makePair $ sortVecs $ L.map V.fromList $-          L.filter (\c -> L.length c > 1) $ getMaximalCliques (\i j -> snd $ graph ! i ! j) [0..numVerts - 1]--      --generates faces of simplices and records the boundary indices-      combos :: Int -> Int -> Vector (Vector (Vector Int)) -> Vector (Vector (Vector Int, Vector Int)) -> Vector (Vector (Vector Int, Vector Int))-      combos i max sc result =-        if i == max then-          (V.map (\s -> (s, V.empty)) $ V.last sc) `cons` result --don't record boundary indices for edges-        else-          let i1        = i + 1 --sc is in reverse order, so sc !! i1 is the array of simplices one dimension lower-              current   = sc ! i-              next      = sc ! i1-              len       = V.length next-              allCombos = V.map getCombos current --get all the faces of every simplex-              uCombos   = bigU allCombos --union of faces-              --the index of the faces of each simplex can be found by adding the number of (n-1)-simplices to the index of each face in the union of faces-              indices   = parMapVec (V.map (\face -> len + (elemIndexUnsafe face uCombos))) allCombos-          in combos i1 max (replaceElem i1 (next V.++ uCombos) sc) $ (V.zip current indices) `cons` result--      fstmc2 = fst maxCliques - 2-  in-    if fstmc2 == (-1) then (numVerts, V.empty) --if there are no maximal cliques, the complex is just a bunch of points-    else-      let sc = combos 0 fstmc2 (snd maxCliques) V.empty-      in (numVerts, sc)--{- |-  Constructs the Vietoris-Rips complex given a scale, metric, and data set.-  Also uses O(n^2) memory (where n is the number of data points) for a graph storing all the distances between data points.-  This utilizes any available processors in parallel because the construction is quite expensive.--}-makeVRComplexFast :: (Ord a, Eq b) => a -> (b -> b -> a) -> [b] -> (SimplicialComplex, Graph a)-makeVRComplexFast scale metric dataSet =-  let graph = makeNbrhdGraph scale metric dataSet-      sc    = makeCliqueComplex graph-  in (sc, graph)--{- |-  Constructs the Vietoris-Rips complex given a scale, metric, and data set.-  This utilizes any available processors in parallel because the construction is quite expensive.--}-makeVRComplexLight :: (Ord a, Eq b) => a -> (b -> b -> a) -> [b] -> SimplicialComplex-makeVRComplexLight scale metric dataSet =-  let numVerts = L.length dataSet-      vector   = V.fromList dataSet--      --make a list with an entry for every dimension of simplices-      organizeCliques 1 _       = []-      organizeCliques i cliques =-        let helper = biFilter (\simplex -> i == V.length simplex) cliques --find the simplices with the given number of vertices-        in (fst helper):(organizeCliques (i - 1) $ snd helper) --append them to the next iteration of the function--      makePair simplices = --pair the organized maximal cliques with the dimension of the largest clique-        case simplices of-          (x:_) ->-            let dim = V.length x-            in (dim, organizeCliques dim $ V.fromList simplices)-          []    -> (1, []) --if there are no maximal cliques this acts as a flag so that the algorithm doesn't try to generate all the other simplices--      maxCliques :: (Int, Vector (Vector (Vector Int)))-      maxCliques = --find all maximal cliques and sort them from largest to smallest (excludes maximal cliques which are single points)-        (\(x, y) -> (x, V.fromList y)) $ makePair $ sortVecs $ L.map V.fromList $-          L.filter (\c -> L.length c > 1) $ getMaximalCliques (\i j -> metric (vector ! i) (vector ! j) <= scale) [0..numVerts - 1]--      --generates faces of simplices and records the boundary indices-      combos :: Int -> Int -> Vector (Vector (Vector Int)) -> Vector (Vector (Vector Int, Vector Int)) -> Vector (Vector (Vector Int, Vector Int))-      combos i max sc result =-        if i == max then-          (V.map (\s -> (s, V.empty)) $ V.last sc) `cons` result --don't record boundary indices for edges-        else-          let i1        = i + 1 --sc is in reverse order, so sc !! i1 is the array of simplices one dimension lower-              current   = sc ! i-              next      = sc ! i1-              len       = V.length next-              allCombos = V.map getCombos current --get all the faces of every simplex-              uCombos   = bigU allCombos --union of faces-              --the index of the faces of each simplex can be found by adding the number of (n-1)-simplices to the index of each face in the union of faces-              indices   = parMapVec (V.map (\face -> len + (elemIndexUnsafe face uCombos))) allCombos-          in combos i1 max (replaceElem i1 (next V.++ uCombos) sc) $ (V.zip current indices) `cons` result--      fstmc2 = fst maxCliques - 2-  in-    if fstmc2 == (-1) then (numVerts, V.empty) --if there are no maximal cliques, the complex is just a bunch of points-    else-      let sc = combos 0 fstmc2 (snd maxCliques) V.empty-      in (numVerts, sc)----BOOLEAN HOMOLOGY------------------------------------------------------------gets the first boundary operator (because edges don't need to point to their subsimplices)-makeEdgeBoundariesBool :: SimplicialComplex -> BMatrix-makeEdgeBoundariesBool sc =-  transposeMat $ V.map (\edge ->-    V.map (\vert -> vert == V.head edge || vert == V.last edge) $ 0 `range` (fst sc - 1)) $-      V.map fst $ V.head $ snd sc----gets the boundary coefficients for a simplex of dimension 2 or greater---first argument is dimension of the simplex---second argument is the simplicial complex---third argument is the simplex paired with the indices of its faces-makeSimplexBoundaryBool :: Int -> SimplicialComplex -> (Vector Int, Vector Int) -> Vector Bool-makeSimplexBoundaryBool dim simplices (simplex, indices) =-  mapWithIndex (\i s -> V.elem i indices) (V.map fst $ (snd simplices) ! (dim - 2))----makes boundary operator for all simplices of dimension 2 or greater---first argument is the dimension of the boundary operator, second is the simplicial complex-makeBoundaryOperatorBool :: Int -> SimplicialComplex -> BMatrix-makeBoundaryOperatorBool dim sc = transposeMat $ V.map (makeSimplexBoundaryBool dim sc) $ (snd sc) ! (dim - 1)----makes all the boundary operators-makeBoundaryOperatorsBool :: SimplicialComplex -> Vector BMatrix-makeBoundaryOperatorsBool sc =-  let dim = getDim sc-      calc i-        | i > dim   = V.empty-        | i == 1    = (makeEdgeBoundariesBool sc) `cons` (calc 2)-        | otherwise = (makeBoundaryOperatorBool i sc) `cons` (calc (i + 1))-  in calc 1--{- |-  The nth index of the output corresponds to the nth Betti number.-  The zeroth Betti number is the number of connected components, the 1st is the number of tunnels/punctures,-  the 2nd is the number of hollow volumes, and so on for higher-dimensional holes.--}-bettiNumbers :: SimplicialComplex -> [Int]-bettiNumbers sc =-  let dim      = (getDim sc) + 1-      boundOps = makeBoundaryOperatorsBool sc-      ranks    = --dimension of image paired with dimension of kernel-        (0, V.length $ V.head boundOps) `cons`-          (V.map (\op -> let rank = rankBool op in (rank, (V.length $ V.head op) - rank)) boundOps)-      calc 1   = [(snd $ ranks ! 0) - (fst $ ranks ! 1)]-      calc i   =-        let i1 = i - 1-        in-          if i == dim then (snd $ V.last ranks):(calc i1)-          else ((snd $ ranks ! i1) - (fst $ ranks ! i)):(calc i1)-  in-    if L.null $ snd sc then [fst sc]-    else L.reverse $ calc dim---- | Calculates all of the Betti numbers in parallel.-bettiNumbersPar :: SimplicialComplex -> [Int]-bettiNumbersPar sc =-  let dim      = (getDim sc) + 1-      boundOps = makeBoundaryOperatorsBool sc-      ranks    = --dimension of image paired with dimension of kernel-        (0, V.length $ V.head boundOps) `cons`-          (parMapVec (\op -> let rank = rankBool op in (rank, (V.length $ V.head op) - rank)) boundOps)-      calc 1   = [(snd $ ranks ! 0) - (fst $ ranks ! 1)]-      calc i   =-        let i1 = i - 1-        in-          if i == dim then evalPar (snd $ V.last ranks) (calc i1) --see Util for evalPar-          else evalPar ((snd $ ranks ! i1) - (fst $ ranks ! i)) (calc i1)-  in-    if L.null $ snd sc then [fst sc]-    else L.reverse $ calc dim

@@ -1,327 +0,0 @@-{- |-Module     : Persistence.Util-Copyright  : (c) Eben Cowley, 2018-License    : BSD 3 Clause-Maintainer : eben.cowley42@gmail.com-Stability  : experimental--This module contains miscellaneous utility functions used throughout the Persistence library.---}--module Util where--import Data.List as L-import Data.Vector as V-import Control.Parallel.Strategies---- | Simple instance of Num where True is 1 and False is 0, all operations work like arithmetic modulo 2.-instance Num Bool where-  p + q  = p `xor` q-  p * q  = p && q-  p - q  = p `xor` (not q)-  negate = not-  abs    = id-  fromInteger 0 = False-  fromInteger _ = True-  signum bool   = if bool then 1 else 0---- | Exclusive or.-xor :: Bool -> Bool -> Bool-xor False False = False-xor True False  = True-xor False True  = True-xor True True   = False---- | First element of a triple.-one (a, _, _) = a--- | Second element of a triple.-two (_, b, _) = b--- | Third element of a triple.-thr (_, _, c) = c--- | Last two elements of a triple.-not1 (_, b, c) = (b, c)--- | First and last elements of a triple.-not2 (a, _, c) = (a, c)--- | First two elements of a triple.-not3 (a, b, _) = (a, b)---- | Concatenate a vector of vectors.-flatten :: Vector (Vector a) -> Vector a-flatten = V.foldl1 (V.++)---- | Multiply a vector by a scalar.-mul :: Num a => a -> Vector a -> Vector a-mul s = V.map (*s)--{- |-  Add two vectors together component-wise.-  WARNING: If one vector is longer than the other,-  the longer vector will simply be cut off.--}-add :: Num a => Vector a -> Vector a -> Vector a-add = V.zipWith (+)--{- |-  Subtract the second vector from the first vector component-wise.-  WARNING: If one vector is longer than the other,-  the longer vector will simply be cut off.--}-subtr :: Num a => Vector a -> Vector a -> Vector a-subtr = V.zipWith (\x y -> x - y)---- | Dot product.-dotProduct :: Num a => Vector a -> Vector a -> a-dotProduct vec1 vec2-  | a && b = fromIntegral 0-  | a      = error "First vector passed to dotProduct too short."-  | b      = error "Second vector passed to dotProduct too short."-  | otherwise   = (V.head vec1)*(V.head vec2) + (dotProduct (V.tail vec1) (V.tail vec2))-    where a = V.null vec1; b = V.null vec2---- | Extended Euclidean algorithm. Finds the gcd of the two inputs plus the coefficients that multiply each input and sum to give the gcd.-extEucAlg :: Integral a => a -> a -> (a, a, a)-extEucAlg a b =-  let eeaHelper r s t =-        case snd r of-          0 -> (fst r, fst s, fst t)-          _ ->-            let r1    = fst r-                r2    = snd r-                s2    = snd s-                t2    = snd t-                q     = r1 `div` r2-                nextr = r1 - q*r2-                nexts = fst s - q*s2-                nextt = fst t - q*t2-            in eeaHelper (r2, nextr) (s2, nexts) (t2, nextt)-  in (\(x, y, z) -> if x < 0 then (-x, -y, -z) else (x, y, z)) $ eeaHelper (a, b) (0, 1) (1, 0)---- | Returns whether or not the first number divides the second number.-divides :: Int -> Int -> Bool-0 `divides` b = False-a `divides` b-  | b < 0     = False-  | b == 0    = True-  | otherwise = a `divides` (b - (abs a))---- | Switches the elements of the vector at the given indices.-switchElems ::Int -> Int -> Vector a -> Vector a-switchElems i j vector-  | j == i    = vector-  | j < i     =-    let first  = V.take j vector-        second = V.drop (j + 1) (V.take i vector)-        third  = V.drop (i + 1) vector-    in first V.++ (cons (vector ! i) second) V.++ (cons (vector ! j) third)-  | otherwise =-    let first  = V.take i vector-        second = V.drop (i + 1) (V.take j vector)-        third  = V.drop (j + 1) vector-    in first V.++ (cons (vector ! j) second) V.++ (cons (vector ! i) third)---- | Return all vectors missing exactly one element from the original vector.-getCombos :: Vector a -> Vector (Vector a)-getCombos vector =-  let len    = V.length vector-      calc i =-        if i == len then empty-        else-          let i1 = i + 1-          in ((V.take i vector) V.++ (V.drop i1 vector)) `cons` (calc i1)-  in calc 0---- | Returns whether or not every element satisfies the predicate.-forallVec :: (a -> Bool) -> Vector a -> Bool-forallVec p vector =-  if V.null vector then True-  else (p $ V.head vector) && (forallVec p $ V.tail vector)---- | Map a function that takes into account the index of each element.-mapWithIndex :: (Int -> a -> b) -> Vector a -> Vector b-mapWithIndex f vector =-  let helper i vec =-        if V.null vec then empty-        else cons (f i $ V.head vec) $ helper (i + 1) (V.tail vec)-  in helper 0 vector---- | Map a function that takes into account the index of each element in parallel.-parMapWithIndex :: (Int -> a -> b) -> Vector a -> Vector b-parMapWithIndex f vector =-  let helper i vec = runEval $-        if V.null vec then return empty-        else-          let current = f i $ V.head vec; rest = helper (i + 1) $ V.tail vec-          in rpar current >> rseq rest >> (return $ current `cons` rest)-  in helper 0 vector---- | Return the element satisfying the predicate and its index if it exists.-elemAndIndex :: (a -> Bool) -> Vector a -> Maybe (a, Int)-elemAndIndex p vector =-  let helper i vec-        | V.null vec     = Nothing-        | p $ V.head vec = Just (V.head vec, i)-        | otherwise      = helper (i + 1) $ V.tail vec-  in helper 0 vector---- | Return the elements satisfying the predicate and their indices.-elemAndIndices :: (a -> Bool) -> Vector a -> [(a, Int)]-elemAndIndices p vector =-  let helper i vec-        | V.null vec     = []-        | p $ V.head vec = (V.head vec, i) : (helper (i + 1) $ V.tail vec)-        | otherwise      = helper (i + 1) $ V.tail vec-  in helper 0 vector---- | 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 len = V.length vector1- -      calc i result =-        let a = vector1 ! i-        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--{- |-  Return the vector of elements that satisfy the predicate in the first component and-  the vector of elements that don't satisfy the predicate in the second component.--}-biFilter :: (a -> Bool) -> Vector a -> (Vector a, Vector a)-biFilter p vector =-  let calc true false v-        | V.null v  = (true, false)-        | p x       = calc (true `snoc` x) false $ V.tail v-        | otherwise = calc true (false `snoc` x) $ V.tail v-        where x = V.head v-  in calc V.empty V.empty vector---- | Orders a list of vectors from greatest to least length.-sortVecs :: [Vector a] -> [Vector a]-sortVecs []     = []-sortVecs (v:vs) =-  let len  = V.length v-      less = sortVecs $ L.filter (\u -> V.length u < len) vs-      more = sortVecs $ L.filter (\u -> V.length u >= len) vs-  in more L.++ [v] L.++ less---- | Parallel map a function over a vector.-parMapVec :: (a -> b) -> Vector a -> Vector b-parMapVec f v = runEval $ evalTraversable rpar $ V.map f v---- | Filter a vector with a predicate that takes into account the index of the element.-filterWithIndex :: (Int -> a -> Bool) -> Vector a -> Vector a-filterWithIndex p vector =-  let maxIndex = V.length vector - 1-      calc i-        | i == maxIndex    = V.empty-        | p i (vector ! i) = (vector ! i) `cons` calc (i + 1)-        | otherwise        = calc (i + 1)-  in calc 0---- | Generate a range of integers in vector form.-range :: Int -> Int -> Vector Int-range x y-  | x == y = x `cons` empty-  | x < y  = x `cons` (range (x + 1) y)-  | x > y  = (range x (y + 1)) `snoc` y---- | Replace the element at the given index with the given element.-replaceElem :: Int -> a -> Vector a -> Vector a-replaceElem i e v = (V.take i v) V.++ (e `cons` (V.drop (i + 1) v))---- | Quicksort treating the given predicate as the < operator. Works like this because its more convenient to make a lambda instead of a complete instance of Ord.-quicksort :: (a -> a -> Bool) -> Vector a -> Vector a-quicksort rel vector = --rel is the > operator-  if V.null vector then empty-  else-    let x       = V.head vector-        xs      = V.tail vector-        lesser  = V.filter (rel x) xs-        greater = V.filter (not . (rel x)) xs-    in (quicksort rel lesser) V.++ (x `cons` (quicksort rel greater))---- | Takes the union of all of the vectors.-bigU :: Eq a => Vector (Vector a) -> Vector a-bigU =-  let exists x v-        | V.null v      = False-        | V.head v == x = True-        | otherwise     = exists x (V.tail v)-      union v1 v2 =-        if V.null v1 then v2-        else-          let x = V.head v1-          in-            if exists x v2 then union (V.tail v1) v2-            else union (V.tail v1) (x `cons` v2)-  in V.foldl1 union---- | The element being searched for, the vector being searched, and the lower and upper bounds on the indices.-binarySearch :: Ord a => a -> Vector a -> Int -> Int -> Maybe Int-binarySearch value xs low high-  | high < low        = Nothing-  | xs ! mid > value  = binarySearch value xs low (mid - 1)-  | xs ! mid < value  = binarySearch value xs (mid + 1) high-  | otherwise         = Just mid-  where mid = low + ((high - low) `div` 2)---- | Intersection of SORTED vectors.-(|^|) :: Ord a => Vector a -> Vector a -> Vector a-vector1 |^| vector2 =-  let len          = V.length vector2 - 1-      calc acc v =-        if V.null v then acc-        else-          let x = V.head v; xs = V.tail v-          in case binarySearch x vector2 0 len of-            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-existsVec p v-  | V.null v     = False-  | p $ V.head v = True-  | otherwise    = existsVec p $ V.tail v---- | If the relation were the "greater than" operator, this would find the minimum element of the vector.-foldRelation :: (a -> a -> Bool) -> Vector a -> a-foldRelation rel vec =-  let calc w v-        | V.null v  = w-        | rel w x   = calc x xs-        | otherwise = calc w xs-        where x = V.head v; xs = V.tail v-  in calc (V.head vec) (V.tail vec)---- | Unsafe index finding.-elemIndexUnsafe :: Eq a => a -> Vector a -> Int-elemIndexUnsafe elem vector =-  let find i v-        | V.null v         = error "Element isn't here, Persistence.Util.elemIndexUnsafe"-        | V.head v == elem = i-        | otherwise        = find (i + 1) $ V.tail v-  in find 0 vector--{- |-  Spark the first argument for parallel evaluation and force evaluation of the second argument,-  then return the first argument concatenated to the second. This is useful especially if the second-  argument is a recursive call that calls evalPar again, so that every elemment of the list will be -  sparked for parallelism.--}-evalPar :: a -> [a] -> [a]-evalPar c r = runEval $ rpar c >> rseq r >> return (c:r)

@@ -0,0 +1,67 @@+# Revision history for Persistence++## 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.++## 1.1.3 -- 2018-07-30++### Changed++- Fixed a major bug with persistent homology; high dimensional holes were being detected in low dimensional data sets.+- Persistent homology now filters out bar codes of the form (i, Just i), as they say nothing about the topology of the underlying complexes.++## 1.1.4 -- 2018-09-15++### Changed++- Fixed spelling error.+- Persistence now exports the constructors for `Extended a`.++## 1.1.4.1 -- 2018-09-15++### Changed++- Fixed all spelling errors, should actually build now.++## 1.1.4.2 -- 2018-09-15++### Changed++- Fixed non-exhaustive pattern match in `BottleNeckDistance` functions.++## 2.0 -- To be determined++### Changed++- The module `Persistence` has been renamed to `Filtration,` and all modules now exist within `Persistence`.++- Bar codes now take a type as a parameter so that they can not only represent the filtration indices at which features appear but also the scales.++- Persistent homology functions now take data sets as either vectors or lists.++- Fixed minor issue with bottleneck distance and removed the unsafe functions.++- Improved documentation for all exposed modules++- `Extended` type now exports its constructors.++### Added++- Simplex type synonyms for making other type synonyms and signatures more readable.++- Data structures for filtrations both with and without all vertices having filtration index equal to zero, and persistent homology functions for processing each structure.++- Persistent homology functions that return Bar codes in terms of scales.++- Functions for constructing and manipulating persistence landscapes.

@@ -0,0 +1,895 @@+{- |+Module     : Persistence.Filtration+Copyright  : (c) Eben Kadile, 2018+License    : BSD 3 Clause+Maintainer : eben.cowley42@gmail.com+Stability  : experimental++This module contains functions for constructing filtrations and computing persistent homology, persistence landscapes, and computing bottleneck distance between barcode diagrams.++A filtration is a finite sequence of simplicial complexes where each complex is a subset of the next. This means that a filtration can be thought of as a single simplicial complex where each of the simplices is labeled with a "filtration index" that represents the index in the sequence where that simplex enters the filtration.++One way to create a filtration, given a simplicial complex, a metric for the vertices, and a list of distances, is to loop through the distances from greatest to least: create a simplicial complex each iteration which excludes simplices that contain pairs of vertices which are further than the current distance apart. This method will produce a filtration of Vietoris-Rips complexes - each filtration index will correspond to a Rips complex whose scale is the corresponding distance. This filtration represents the topology of the data at each of the scales with which it was constructed.++NOTE: It's important that, even though the smallest filtration index represents the smallest scale at which the data is being anaylzed, all functions in this library receive your list of scales sorted in *decreasing* order.++An essential thing to note in this library is the distinction between "fast" and "light" functions. Light functions call the metric every time distance between two points is required, which is a lot. Fast functions store the distances between points and access them in constant time, BUT this means they use O(n^2) memory with respect to the number of data points, so it's a really bad idea to use this optimization on substantially large data if you don't have a lot of RAM.++Persistent homology is the main event of topological data analysis. It allows one to identify clusters, tunnels, cavities, and higher dimensional holes that persist in the data throughout many scales. The output of the persistence algorithm is a barcode diagram. A single barcode represents the filtration index where a feature appears and the index where it disappears (if it does). Alternatively, a barcode can represent the scale at which a feature and the scale at which it ends. Thus, short barcodes are typically interpretted as sampling irregularities and long barcodes are interpretted as actual features of whatever the underlying data set represents. In this context, what a feature *is* depends on which dimension the barcode diagram is; 0-dimensional features are connected components, 1-dimensional features are loops or tunnels, 2-dimensional features are hollow volumes, and higher dimensional features correspond to heigher-dimensional cavities.++After you've got the barcodes of a data set, you might want to compare it with that of a different data set. This is the purpose of bottleneck distance, which corresponds to the Hausdorff distance between barcode diagrams.++Another way to compare barcode diagrams is by using persistence landscapes. The peristence landscape of a barcode diagram is a finite sequence of piecewise-linear, real-valued functions. This means they can be used to take averages and compute distances between barcode diagrams. See "A Persistence Landscapes Toolbox For Topological Statistics" by Bubenik and Dlotko for more information.++WARNING: The persistence landscape functions have not been fully tested. Use them with caution. If you get any errors or unexpected output, please don't hesitate to email me.++-}++module Persistence.Filtration (+  -- * Types+    FilterSimplex+  , SimpleFiltration+  , Filtration+  , Extended (MinusInfty, Finite, Infinity)+  , BarCode+  , Landscape+  -- * Utilities+  , sim2String+  , filtr2String+  , getComplex+  , getDimension+  , simple2Filtr+  -- * Construction+  , filterByWeightsFast+  , makeRipsFiltrationFast+  , makeRipsFiltrationFastPar+  , filterByWeightsLight+  , makeRipsFiltrationLight+  , makeRipsFiltrationLightPar+  -- * Persistent homology+  , indexBarCodes+  , indexBarCodesSimple+  , scaleBarCodes+  , scaleBarCodesSimple+  -- * Comparing barcode diagrams+  , indexMetric+  , bottleNeckDistance+  , bottleNeckDistances+  , calcLandscape+  , evalLandscape+  , evalLandscapeAll+  , linearComboLandscapes+  , avgLandscapes+  , diffLandscapes+  , normLp+  , metricLp+  ) where++import Persistence.Util+import Persistence.SimplicialComplex++import Data.Maybe+import Data.List as L+import Data.Vector as V+import Data.Algorithm.MaximalCliques++import Control.Parallel.Strategies++-- * Types++{- |+  This type synonym exists to make other synonyms more concise.+  Each simplex in a filtration is represented as a triple: its filtration index,+  the indices of its vertices in the original data, and the indices of its faces in the next lowest dimension.+  Edges do not have reference to their faces, as it would be redundant with their vertices.+  All simplices are sorted according to filtration index upon construction of the filtration.+  In each dimension, all simplices are sorted in increasing order of filtration index,+  and every simplices face indices are sorted in decreasing order;+  both of these facts are critical to the computation of persistent homology.+-}+type FilterSimplex = (Int, Vector Int, Vector Int)++{- |+  A type representing a filtration whose vertices all have filtration index 0.+  Slightly faster and slightly less memory usage. The first component is simply the number of vertices.+  The second component is a vector with an entry for each dimension of simplices, starting at dimension 1 for edges.+-}+type SimpleFiltration = (Int, Vector (Vector FilterSimplex))++{- |+  Representation of a filtration which, unlike `SimpleFiltration`, can cope with vertices that have a non-zero filtration index. Vertices of the filtration are represented like all other simplices except that they don't their own have vertices or faces.++  Note that, since this library currently only deals with static pointcloud data, all of the filtration construction functions produce vertices whose filtration index is 0. Thus, if you want to use this type you will have to construct the instances yourself.+-}+type Filtration = Vector (Vector FilterSimplex)++{- |+  A type extending the number line to positive and negative infinity.+  Used for representing infinite barcodes, bottleneck distance, and persistence landscapes.+-}+data Extended a = Finite a+                | Infinity+                | MinusInfty+                deriving (Eq, Show)++{- |+  The ordering is inherited from the type a,+  Infinity is greater than everything else and MinusInfty is less than everything else.+-}+instance (Ord a, Eq a) => Ord (Extended a) where++  _ > Infinity        = False+  Infinity > _        = True+  Finite a > Finite b = a > b+  MinusInfty > _      = False+  _ > MinusInfty      = True++  Infinity >= _            = True+  _ >= Infinity            = False+  Finite a >= Finite b     = a >= b+  MinusInfty >= MinusInfty = True+  MinusInfty >= _          = False++  Infinity < _            = False+  _ < Infinity            = True+  Finite a < Finite b     = a < b+  MinusInfty < MinusInfty = False+  MinusInfty < _          = True++  _ <= Infinity        = True+  Infinity <= _        = False+  Finite a <= Finite b = a <= b+  MinusInfty <= _      = True+  _ <= MinusInfty      = False++-- | Arithmetic is defined in the canonical way based on the arithmetic of `a`.+instance Num a => Num (Extended a) where++  _ + Infinity        = Infinity+  _ + MinusInfty      = MinusInfty+  Infinity   + _      = Infinity+  MinusInfty + _      = MinusInfty+  Finite x + Finite y = Finite (x + y)++  _ - Infinity        = MinusInfty+  _ - MinusInfty      = Infinity+  Infinity   - _      = Infinity+  MinusInfty - _      = MinusInfty+  Finite x - Finite y = Finite (x - y)++  _ * Infinity        = Infinity+  _ * MinusInfty      = MinusInfty+  Infinity   * _      = Infinity+  MinusInfty * _      = MinusInfty+  Finite x * Finite y = Finite (x * y)++  abs Infinity    = Infinity+  abs MinusInfty  = Infinity+  abs (Finite x)  = Finite $ abs x++  fromInteger = Finite . fromInteger++  signum Infinity   = Finite (fromInteger 1)+  signum MinusInfty = Finite (fromInteger (-1))+  signum (Finite x) = Finite (signum x)++fi :: (Integral a, Floating b) => Extended a -> Extended b+fi (Finite i) = Finite $ fromIntegral i+fi Infinity   = Infinity+fi MinusInfty = MinusInfty++unExtend :: Floating a => Extended a -> a+unExtend Infinity   = 1.0/0.0+unExtend (Finite x) = x+unExtend MinusInfty = -1.0/0.0++-- | (x, Finite y) is a topological feature that appears at the index or scale x and disappears at the index or scale y. (x, Infinity) begins at x and doesn't disappear.+type BarCode a = (a, Extended a)++{- |+  A Persistence landscape is a certain type of piecewise linear function based on a barcode diagram.+  It can be represented as a list of critical points paired with critical values.+  Useful for taking averages and differences between barcode diagrams.+-}+type Landscape = Vector (Vector (Extended Double, Extended Double))++-- * Utilities++-- | Shows all the information in a simplex.+sim2String :: FilterSimplex -> String+sim2String (index, vertices, faces) =+  "Filtration index: " L.++ (show index) L.+++    "; Vertex indices: " L.++ (show vertices) L.+++      "; Boundary indices: " L.++ (show faces) L.++ "\n"++-- | Shows all the information in a filtration.+filtr2String :: Either SimpleFiltration Filtration -> String+filtr2String (Left f)  =+  "Simple filtration:\n" L.++ ((intercalate "\n") $ toList+    $ V.map (L.concat . toList . (V.map sim2String)) $ snd f)+filtr2String (Right f) =+  (intercalate "\n") $ toList $ V.map (L.concat . toList . (V.map sim2String)) f++{- |+  Gets the simplicial complex specified by the filtration index.+  This is O(n) with respect to the number of simplices.+-}+getComplex :: Int -> Either SimpleFiltration Filtration -> SimplicialComplex+getComplex index (Left (n, simplices)) =+  (n, dropRightWhile V.null+    $ V.map (V.map not1 . V.filter (\(i, _, _) -> i <= index)) simplices)+getComplex index (Right simplices)     =+  (V.length $ V.filter (\v ->+    one v <= index) (V.head simplices),+      dropRightWhile V.null $ V.map (V.map not1+        . V.filter (\(i, _, _) -> i <= index)) (V.tail simplices))++-- | Return the dimension of the highest dimensional simplex in the filtration (constant time).+getDimension :: Either SimpleFiltration Filtration -> Int+getDimension (Left sf) = V.length $ snd sf+getDimension (Right f) = V.length f - 1++-- | Convert a simple filtration into an ordinary filtration.+simple2Filtr :: SimpleFiltration -> Filtration+simple2Filtr (n, x) =+  let x' = (V.map (\(i, v, _) -> (i, v, V.reverse v)) $ V.head x) `cons` (V.tail x)+  in (mapWithIndex (\i (a,b,c) ->+       (a,i `cons` V.empty,c)) $ V.replicate n (0, V.empty, V.empty)) `cons` x'++-- * Construction++{- |+  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.+  This is really a helper function to be called by makeRipsFiltrationFast,+  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.+-}+filterByWeightsFast :: Ord a+                    => Either (Vector a) [a] -- ^Scales in decreasing order+                    -> (SimplicialComplex, Graph a) -- ^Simplicial complex and a graph encoding the distance between every data point as well as whether or not they are within the largest scale of each other.+                    -> SimpleFiltration+filterByWeightsFast scales' ((numVerts, simplices), graph) =+  let scales = case scales' of Left v -> V.toList v; Right l -> l+      edgeInSimplex edge simplex =+        (V.any (\x -> V.head edge == x) simplex) && (V.any (\x -> V.last edge == x) simplex)+      edgeTooLong scale edge     = scale <= (fst $ graph ! (edge ! 0) ! (edge ! 1))+      maxIndex                   = (L.length scales) - 1++      calcIndices 0 [] sc         = sc+      calcIndices i (scl:scls) sc =+        --find edges excluded by this scale+        let longEdges = V.filter (edgeTooLong scl) $ V.map (\(i, v, f) -> v) $ V.head sc+        in calcIndices (i - 1) scls $ V.map (V.map (\(j, v, f) ->+          --if the simplex has not yet been assigned a fitration index+          if j == 0 then+            --if a long edge is in the simplex, assign it the current index+            if V.any (\edge -> edgeInSimplex edge v) longEdges then (i, v, f)+            --otherwise wait until next iteration+            else (0, v, f)+          --otherwise leave it alone+          else (j, v, f))) sc++      sortFiltration simplices =+        let sortedSimplices =+              --sorted in reverse order+              V.map (quickSort (\((i, _, _), _) ((j, _, _), _) -> i > j)) $+                V.map (mapWithIndex (\i s -> (s, i))) simplices+            newFaces dim (i, v, f) =+              let findNew j =+                    case V.findIndex (\x -> snd x == j) $ sortedSimplices ! (dim - 1) of+                      Just k  -> k+                      Nothing -> error "Persistence.Filtration.sortFiltration.newFaces.findNew. This is a bug. Please email the Persistence maintainers."+              in (i, v, (V.map findNew f))+        in+          if V.null simplices then simplices+          else mapWithIndex (\i ss -> V.map ((newFaces i) . fst) ss) sortedSimplices++      sortBoundaries = V.map (V.map (\(i, v, f) -> (i, v, quickSort (\a b -> a <= b) f)))++  --sort the simplices by filtration index,+  --then sort boundaries so that the boundary chains can be acquired easily+  in (numVerts, sortBoundaries $ sortFiltration $+      calcIndices maxIndex (L.tail scales) $+        V.map (V.map (\(v, f) -> (0, v, f))) $ simplices)++{- |+  This function constructs a filtration of the Vietoris-Rips complexes associated with the scales.+  Note that this a fast function, meaning it uses O(n^2) memory to quickly access distances where n is the number of data points.+-}+makeRipsFiltrationFast :: (Ord a, Eq b)+                       => Either (Vector a) [a] -- ^Scales in decreasing order+                       -> (b -> b -> a) -- ^Metric+                       -> Either (Vector b) [b] -- ^Data set+                       -> SimpleFiltration+makeRipsFiltrationFast scales metric =+  let scale = case scales of Left v -> V.head v; Right l -> L.head l+  in (filterByWeightsFast scales) . (makeRipsComplexFast scale metric)++{- |+  Same as above except it uses parallelism when computing the Vietoris-Rips complex of the largest scale.+-}+makeRipsFiltrationFastPar :: (Ord a, Eq b)+                          => Either (Vector a) [a] -- ^Scales in decreasing order+                          -> (b -> b -> a) -- ^Metric+                          -> Either (Vector b) [b] -- ^Data set+                          -> SimpleFiltration+makeRipsFiltrationFastPar scales metric =+  let scale = case scales of Left v -> V.head v; Right l -> L.head l+  in (filterByWeightsFast scales) . (makeRipsComplexFastPar scale metric)++{- |+  The same as filterbyWeightsFast except it uses far less memory at the cost of speed.+  Note that the scales must be in decreasing order.+-}+filterByWeightsLight :: Ord a+                     => Either (Vector a) [a] -- ^Scales in decreasing order+                     -> (b -> b -> a) -- ^Metric+                     -> Either (Vector b) [b] -- ^Data set+                     -> SimplicialComplex -- ^Vietoris-Rips complex of the data at the largest scale.+                     -> SimpleFiltration+filterByWeightsLight scales' metric dataSet (numVerts, simplices) =+  let scales = case scales' of Left v -> V.toList v; Right l -> l++      edgeInSimplex edge simplex =+        (V.any (\x -> V.head edge == x) simplex) && (V.any (\x -> V.last edge == x) simplex)+      vector                     = case dataSet of Left v -> v; Right l -> V.fromList l+      edgeTooLong scale edge     = scale <= (metric (vector ! (edge ! 0)) (vector ! (edge ! 1)))+      maxIndex                   = (L.length scales) - 1++      calcIndices 0 [] sc         = sc+      calcIndices i (scl:scls) sc =+        --find edges excluded by this scale+        let longEdges = V.filter (edgeTooLong scl) $ V.map (\(i, v, f) -> v) $ V.head sc+        in calcIndices (i - 1) scls $ V.map (V.map (\(j, v, f) ->+          --if the simplex has not yet been assigned a fitration index+          if j == 0 then+            --if a long edge is in the simplex, assign it the current index+            if V.any (\edge -> edgeInSimplex edge v) longEdges then (i, v, f)+            --otherwise wait until next iteration+            else (0, v, f)+          --otherwise leave it alone+          else (j, v, f))) sc++      sortFiltration simplices =+        let sortedSimplices =+              --sorted in increasing order+              V.map (quickSort (\((i, _, _), _) ((j, _, _), _) -> i > j)) $+                V.map (mapWithIndex (\i s -> (s, i))) simplices+            newFaces dim (i, v, f) =+              let findNew j =+                    case V.findIndex (\x -> snd x == j) $ sortedSimplices ! (dim - 1) of+                      Just k  -> k+                      Nothing -> error "Persistence.Filtration.filterByWeightsLight.sortFiltration.newFaces.findNew. This is a bug. Please email the Persistence maintainers."+              in (i, v, (V.map findNew f))+        in+          if V.null simplices then simplices+          else mapWithIndex (\i ss -> V.map ((newFaces i) . fst) ss) sortedSimplices++  in (numVerts, sortFiltration $ --sort the simplices by filtration index+      calcIndices maxIndex (L.tail scales) $+        V.map (V.map (\(v, f) -> (0, v, f))) $ simplices)++{- |+  Constructs the filtration of Vietoris-Rips complexes corresponding to each of the scales.+-}+makeRipsFiltrationLight :: (Ord a, Eq b)+                        => Either (Vector a) [a] -- ^List of scales in decreasing order+                        -> (b -> b -> a) -- ^Metric+                        -> Either (Vector b) [b] -- ^Data set+                        -> SimpleFiltration+makeRipsFiltrationLight scales metric dataSet =+  let scale = case scales of Left v -> V.head v; Right l -> L.head l+  in filterByWeightsLight scales metric dataSet $ makeRipsComplexLight scale metric dataSet++{- |+  Same as above except it uses parallelism when computing the Vietoris-Rips complex of the largest scale.+-}+makeRipsFiltrationLightPar :: (Ord a, Eq b)+                        => Either (Vector a) [a] -- ^List of scales in decreasing order+                        -> (b -> b -> a) -- ^Metric+                        -> Either (Vector b) [b] -- ^Data set+                        -> SimpleFiltration+makeRipsFiltrationLightPar scales metric dataSet =+  let scale = case scales of Left v -> V.head v; Right l -> L.head l+  in filterByWeightsLight scales metric dataSet $ makeRipsComplexLightPar scale metric dataSet++-- * Persistent Homology++type Chain = Vector Int --indices of the simplices in the sum++{- |+  The nth entry in the list will describe the n-dimensional topology of the filtration.+  That is, the first list will represent clusters, the second list will represent tunnels or punctures, the third will represent hollow volumes,+  and the nth index list will represent n-dimensional holes in the data.+  Features are encoded by the filtration indices where they appear and disappear.+-}+indexBarCodes :: Filtration -> Vector (Vector (BarCode Int))+indexBarCodes filtration =+  let maxdim = getDimension (Right filtration)++      --given a vector of indices of simplices which are marked+      --and a vector of boundary chains paired with the indices of their simplices+      --remove the unmarked simplices from the chain+      removeUnmarked :: Vector Int -> Vector (Int, Chain) -> Vector (Int, Chain)+      removeUnmarked marked = V.map (\(i, c) -> (i, V.filter (\j -> V.elem j marked) c))++      --eliminate monomials in the boundary chain until it is no longer+      --or there is a monomial which can't be eliminated+      removePivotRows :: Vector (Maybe Chain) -> Chain -> Chain+      removePivotRows slots chain =+        if V.null chain then V.empty+        else+          case slots ! (V.head chain) of+            Nothing -> chain+            Just c  -> removePivotRows slots (chain `uin` c)++      --given the indices of the marked simplices from the last iteration,+      --slots from the last iteration,and boundary chains+      --mark the appropriate simplices, fill in the appropriate slots, and identify bar codes+      --boundary chains are paired with the index of their coresponding simplex+      makeFiniteBarCodes :: Int+                         -> Vector Int+                         -> Vector (Maybe Chain)+                         -> Vector (Int, Chain)+                         -> Vector (BarCode Int)+                         -> (Vector Int, Vector (Maybe Chain), Vector (BarCode Int))+      makeFiniteBarCodes dim newMarked slots boundaries barcodes =+        if V.null boundaries then (newMarked, slots, barcodes)+        else+          let boundary = V.head boundaries+              reduced  = removePivotRows slots $ snd boundary+          in+            --mark the simplex if its boundary chain is reduced to null+            if V.null reduced then+              makeFiniteBarCodes dim+                (newMarked `snoc` (fst boundary)) slots (V.tail boundaries) barcodes+            else+              let pivot = V.head reduced+              --put the pivot chain in the pivot's slot, add the new barcode to the list+              in makeFiniteBarCodes dim newMarked (replaceElem pivot+                (Just reduced) slots) (V.tail boundaries) ((one $ filtration ! (dim - 1) ! pivot,+                  Finite $ one $ filtration ! dim ! (fst boundary)) `cons` barcodes)++      --get the finite bar codes for each dimension+      loopFiniteBarCodes :: Int+                         -> Vector (Vector Int)+                         -> Vector (Vector (Maybe Chain))+                         -> Vector (Vector (BarCode Int))+                         -> ( Vector (Vector Int), Vector (Vector (Maybe Chain))+                            , Vector (Vector (BarCode Int)))+      loopFiniteBarCodes dim marked slots barcodes =+        --the slots vector made when looping over the vertices will be null+        if dim > maxdim then (marked, V.tail slots, (V.tail barcodes) V.++ (V.empty `cons` V.empty))+        else+          let numSlots = if dim == 0 then 0 else V.length $ filtration ! (dim - 1) --see above+              boundaries =+                removeUnmarked (V.last marked)+                  $ mapWithIndex (\i (_, _, f) -> (i, f)) $ filtration ! dim+              (newMarked, newSlots, newCodes) =+                makeFiniteBarCodes dim V.empty (V.replicate numSlots Nothing) boundaries V.empty+          in loopFiniteBarCodes (dim + 1) (marked `snoc` newMarked)+            (slots `snoc` newSlots) (barcodes V.++ (newCodes `cons` V.empty))++      --if a simplex isn't marked and has an empty slot,+      --an infinite bar code begins at it's filtration index+      makeInfiniteBarCodes :: Int -> Vector Int -> Vector (Maybe Chain) -> Vector (BarCode Int)+      makeInfiniteBarCodes dim marked slots =+        V.map (\i -> (one $ filtration ! dim ! i, Infinity))+          $ V.filter (\i -> slots ! i == Nothing) marked++      --add the infinite bar codes to the list of bar codes in each dimension+      loopInfiniteBarCodes :: Int+                           -> ( Vector (Vector Int), Vector (Vector (Maybe Chain))+                              , Vector (Vector (BarCode Int)))+                           -> Vector (Vector (BarCode Int))+      loopInfiniteBarCodes dim (marked, slots, barcodes) =+        if dim > maxdim then barcodes+        else+          loopInfiniteBarCodes (dim + 1) (marked, slots, replaceElem dim ((barcodes ! dim)+            V.++ (makeInfiniteBarCodes dim (marked ! dim) (slots ! dim))) barcodes)++  in V.map (V.filter (\(a, b) -> b /= Finite a))+    $ loopInfiniteBarCodes 0 $ loopFiniteBarCodes 0 V.empty V.empty V.empty++-- | Same as above except this function acts on filtrations whose vertices all have filtration index zero (for a very slight speedup).+indexBarCodesSimple :: SimpleFiltration -> Vector (Vector (BarCode Int))+indexBarCodesSimple (numVerts, allSimplices) =+  let removeUnmarked marked = V.filter (\x -> V.elem x marked)++      removePivotRows reduced chain =+        if V.null chain then chain+        else+          case reduced ! (V.head chain) of+            Nothing -> chain+            Just t  -> removePivotRows reduced (chain `uin` t)++      makeBarCodesAndMark :: Int+                          -> Int+                          -> Vector Int+                          -> Vector (Maybe (Vector Int))+                          -> Vector (Int, Vector Int, Vector Int)+                          -> (Vector (BarCode Int), Vector Int)+                          -> (Vector (BarCode Int), Vector Int, Vector Int)+      makeBarCodesAndMark dim index marked reduced simplices (codes, newMarked)+        | V.null simplices = (codes, newMarked, V.findIndices (\x -> x == Nothing) reduced)+        | V.null d         =+          makeBarCodesAndMark dim (index + 1) marked reduced+            (V.tail simplices) (codes, newMarked `snoc` index)+        | otherwise        =+          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)+              ((begin, Finite i) `cons` codes, newMarked)+        where (i, v, f) = V.head simplices+              d         = removePivotRows reduced $ removeUnmarked marked f++      makeEdgeCodes :: Int+                    -> Vector (Maybe (Vector Int))+                    -> Vector (Int, Vector Int, Vector Int)+                    -> (Vector (BarCode Int), Vector Int)+                    -> (Vector (BarCode Int), Vector Int, Vector Int)+      makeEdgeCodes index reduced edges (codes, marked)+        | V.null edges = (codes, marked, V.findIndices (\x -> x == Nothing) reduced)+        | 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) ((0, Finite i) `cons` codes, marked)+        where (i, v, f) = V.head edges+              d         = removePivotRows reduced f++      makeFiniteBarCodes :: Int+                         -> Int+                         -> Vector (Vector (BarCode Int))+                         -> Vector (Vector Int)+                         -> Vector (Vector Int)+                         -> ( Vector (Vector (BarCode Int))+                            , Vector (Vector Int)+                            , Vector (Vector Int)+                            )+      makeFiniteBarCodes dim maxdim barcodes marked slots =+        if dim == maxdim then (barcodes, marked, slots)+        else+          let (newCodes, newMarked, unusedSlots) =+                makeBarCodesAndMark dim 0 (V.last marked) (V.replicate (V.length+                  $ allSimplices ! (dim - 1)) Nothing) (allSimplices ! dim) (V.empty, V.empty)+          in makeFiniteBarCodes (dim + 1) maxdim+            (barcodes V.++ (newCodes `cons` V.empty))+              (marked `snoc` newMarked) (slots `snoc` unusedSlots)++      makeInfiniteBarCodes :: ( Vector (Vector (BarCode Int))+                              , Vector (Vector Int)+                              , Vector (Vector Int)+                              )+                           -> Vector (Vector (BarCode Int))+      makeInfiniteBarCodes (barcodes, marked, unusedSlots) =+        let+            makeCodes :: Int -> Vector (BarCode Int) -> Vector (BarCode Int)+            makeCodes i codes =+              let slots = unusedSlots ! i; marks = marked ! i+              in codes V.++ (V.map (\j -> (one+                   $ allSimplices ! (i - 1) ! j, Infinity)) $ slots |^| marks)+            loop :: Int -> Vector (Vector (BarCode Int)) -> Vector (Vector (BarCode Int))+            loop i v+              | V.null v  = V.empty+              | i == 0    =+                ((V.head v) V.++ (V.map (\j -> (0, Infinity))+                  $ (unusedSlots ! 0) |^| (marked ! 0))) `cons` (loop 1 $ V.tail v)+              | otherwise = (makeCodes i $ V.head v) `cons` (loop (i + 1) $ V.tail v)+        in loop 0 barcodes++      edges    = V.map (\(i, v, f) -> (i, v, (V.reverse v))) $ V.head allSimplices+      numEdges = V.length edges++      (fstCodes, fstMarked, fstSlots) =+        makeEdgeCodes 0 (V.replicate numVerts Nothing) edges (V.empty, V.empty)++      verts = 0 `range` (numVerts - 1)++  in+    V.map (V.filter (\(a, b) -> b /= Finite a))+      $ makeInfiniteBarCodes $ makeFiniteBarCodes 1+        (V.length allSimplices) (fstCodes `cons` V.empty)+          (verts `cons` (fstMarked `cons` V.empty)) (fstSlots `cons` V.empty)++{- |+  The nth entry in the list will describe the n-dimensional topology of the filtration.+  However, features are encoded by the scales where they appear and disappear. For consistency,+  scales must be in decreasing order.+-}+scaleBarCodes :: Either (Vector a) [a] -> Filtration -> Vector (Vector (BarCode a))+scaleBarCodes scales filtration =+  let s = V.reverse $ (\a -> case a of Left v -> v; Right l -> V.fromList l) scales++      translateBarCode (i, Infinity) = (s ! i, Infinity)+      translateBarCode (i, Finite j) = (s ! i, Finite $ s ! j)++  in V.map (V.map translateBarCode) $ indexBarCodes filtration++{- |+  Same as above except acts only on filtrations whose vertices all have filtration index 0.+  Note that scales must be in decreasing order.+-}+scaleBarCodesSimple :: Either (Vector a) [a] -> SimpleFiltration -> Vector (Vector (BarCode a))+scaleBarCodesSimple scales filtration =+  let s = V.reverse $ (\a -> case a of Left v -> v; Right l -> V.fromList l) scales++      translateBarCode (i, Infinity) = (s ! i, Infinity)+      translateBarCode (i, Finite j) = (s ! i, Finite $ s ! j)++  in V.map (V.map translateBarCode) $ indexBarCodesSimple filtration++-- * Comparing barcode diagrams++{- |+  The standard (Euclidean) metric between index barcodes.+  The distance between infinite and finite barcodes is infinite,+  and the distance between two infinite barcodes is the absolute value of the+  difference of their fst component.+-}+indexMetric :: BarCode Int -> BarCode Int -> Extended Double+indexMetric (_, Finite _) (_, Infinity) = Infinity+indexMetric (_, Infinity) (_, Finite _) = Infinity+indexMetric (i, Infinity) (j, Infinity) =+  Finite $ fromIntegral $ abs $ i - j+indexMetric (i, Finite j) (k, Finite l) =+  let x = i - k; y = j - l+  in Finite $ sqrt $ fromIntegral $ x*x + y*y++{- |+  Given a metric, return the Hausdorff distance+  (referred to as bottleneck distance in TDA) between the two sets.+  Returns nothing if either list of barcodes is empty.+-}+bottleNeckDistance :: Ord b+                   => (BarCode a -> BarCode a -> Extended b)+                   -> Vector (BarCode a)+                   -> Vector (BarCode a)+                   -> Maybe (Extended b)+bottleNeckDistance metric diagram1 diagram2+  | V.null diagram1 = Nothing+  | V.null diagram2 = Nothing+  | otherwise       =+    let first  = V.maximum $ V.map (\p -> V.minimum $ V.map (metric p) diagram2) diagram1+        second = V.maximum $ V.map (\p -> V.minimum $ V.map (metric p) diagram1) diagram2+    in Just $ max first second++{- |+  Get's all the bottleneck distances;+  a good way to determine the similarity of the topology of two filtrations.+-}+bottleNeckDistances :: Ord b => (BarCode a -> BarCode a -> Extended b)+                    -> Vector (Vector (BarCode a))+                    -> Vector (Vector (BarCode a))+                    -> Vector (Maybe (Extended b))+bottleNeckDistances metric diagrams1 diagrams2 =+  let d = (L.length diagrams1) - (L.length diagrams2)+  in+    if d >= 0+    then (V.zipWith (bottleNeckDistance metric) diagrams1 diagrams2) V.++ (V.replicate d Nothing)+    else (V.zipWith (bottleNeckDistance metric) diagrams1 diagrams2) V.++ (V.replicate (-d) Nothing)++-- | Compute the persistence landscape of the barcodes for a single dimension.+calcLandscape :: Vector (BarCode Int) -> Landscape+calcLandscape brcds =+  let half = Finite 0.5++      (i,j) `leq` (k,l) = i > k || j <= l++      innerLoop :: (Extended Double, Extended Double)+                -> Vector (Extended Double, Extended Double)+                -> Landscape+                -> Landscape+      innerLoop (b, d) barcodes result =+        case V.findIndex (\(b', d') -> d' > d) barcodes of+          Nothing ->+            outerLoop barcodes (((V.fromList [(0, b), (Infinity, 0)])+              V.++ (V.head result)) `cons` (V.tail result))+          Just i  -> let (b', d') = barcodes ! i in+            if b' >= d then+              if b == d then+                let new = [(Finite 0.0, b')]+                in+                  if d' == Infinity then+                    outerLoop (rmIndex i barcodes) (((V.fromList ((Infinity, Infinity):new))+                      V.++ (V.head result)) `cons` (V.tail result))+                  else+                    innerLoop (b', d') (rmIndex i barcodes)+                      ((V.fromList ((half*(b' + d'), half*(d' - b')):new)+                        V.++ (V.head result)) `cons` (V.tail result))+              else+                let new = [(Finite 0.0, d), (Finite 0.0, b')]+                in+                  if d' == Infinity then+                    outerLoop (rmIndex i barcodes) (((V.fromList ((Infinity, Infinity):new))+                      V.++ (V.head result)) `cons` (V.tail result))+                  else+                    innerLoop (b', d') (rmIndex i barcodes)+                      (((V.fromList ((half*(b' + d'), half*(d' - b')):new))+                        V.++ (V.head result)) `cons` (V.tail result))+            else+              let newbr = (half*(b' + d), half*(d - b'))+              in+                if d' == Infinity then+                  outerLoop (orderedInsert leq newbr barcodes)+                    (((V.fromList [(Infinity, Infinity), newbr])+                      V.++ (V.head result)) `cons` (V.tail result))+                else+                  innerLoop (b', d') (orderedInsert leq newbr barcodes)+                    (((V.fromList [(half*(b' + d'), half*(d' - b')), newbr])+                      V.++ (V.head result)) `cons` (V.tail result))++      outerLoop :: Vector (Extended Double, Extended Double) -> Landscape -> Landscape+      outerLoop barcodes result =+        if not $ V.null barcodes then+          let (b, d) = V.head barcodes+          in+            if (b, d) == (MinusInfty, Infinity)+            then+              outerLoop (V.tail barcodes)+                ((V.fromList [(MinusInfty, Infinity),+                  (Infinity, Infinity)]) `cons` result)+            else if d == Infinity+            then+              outerLoop (V.tail barcodes) ((V.fromList+                [(MinusInfty, Finite 0.0),(b, Finite 0.0),(Infinity,Infinity)]) `cons` result)+            else+              let newCritPoints =+                    if b == Infinity+                    then [(MinusInfty, Infinity)]+                    else [(MinusInfty, Finite 0.0), (half*(b + d), half*(d - b))]+              in innerLoop (b, d) (V.tail barcodes) ((V.fromList newCritPoints) `cons` result)+        else result++  in V.map (quickSort (\(x1, _) (x2, _) -> x1 > x2))+       $ outerLoop (quickSort leq $ V.map (\(i, j) -> (fi $ Finite i, fi j)) brcds) V.empty++-- | Evaluate the nth function in the landscape for the given point.+evalLandscape :: Landscape -> Int -> Extended Double -> Extended Double+evalLandscape landscape i arg =+  let fcn = landscape ! i++      findPointNeighbors :: Ord a => Int -> a -> Vector a -> (Int, Int)+      findPointNeighbors helper x vector =+        let len = V.length vector+            i   = len `div` 2+            y   = vector ! i+        in+          if x == y+          then (helper + i, helper + i)+          else if x > y+          then+            case vector !? (i + 1) of+              Nothing -> (helper + i, helper + i)+              Just z  ->+                if x < z+                then (helper + i, helper + i + 1)+                else findPointNeighbors (helper + i) x $ V.drop i vector+          else+            case vector !? (i - 1) of+              Nothing -> (helper + i, helper + i)+              Just z  ->+                if x > z+                then (helper + i - 1, helper + i)+                else findPointNeighbors helper x $ V.take i vector++      (i1, i2) = findPointNeighbors 0 arg $ V.map fst fcn+      (x1, x2) = (fst $ fcn ! i1, fst $ fcn ! i2)+      (y1, y2) = (snd $ fromMaybe (error "Persistence.Filtration.evalLandscape. This is a bug. Please email the Persistence mainstainers.") $ V.find (\a -> x1 == fst a) fcn, snd $ fromMaybe (error "Persistence.Filtration.evalLandscape. This is a bug. Please email the Persistence mainstainers.") $ V.find (\a -> x2 == fst a) fcn)++  in+    if x1 == x2+    then y1+    else+      case (x1, x2) of+        (MinusInfty, Infinity)   -> arg+        (MinusInfty, Finite _)   -> y1+        (Finite a, Finite b)     ->+          case arg of+            Finite c   ->+              let t = Finite $ (c - a)/(b - a)+              in t*y2 + ((Finite 1.0) - t)*y1+            _          -> error "Persistence.Filtration.evalLandscape.findPointNeighbors. This is a bug. Please email the Persistence maintainers."+        (Finite a, Infinity)     ->+          case arg of+            Infinity   -> y2+            Finite c   ->+              case y2 of+                Infinity   -> Finite $ c - a+                Finite 0.0 -> Finite 0.0+                _          -> error $ "Persistence.Filtration.evalLandscape: y2 = " L.++ (show y2) L.++ ". This is a bug. Please email the Persistence maintainers."+            _          -> error "Persistence.Filtration.evalLandscape.findPointNeighbors: bad argument. This is a bug. Please email the Persistence maintainers."+        anything                 -> error $ "Persistence.Filtration.evalLandscape.findPointNeighbors: " L.++ (show anything) L.++ ". This is a bug. Please email the Persistence maintainers."++-- | Evaluate all the real-valued functions in the landscape.+evalLandscapeAll :: Landscape -> Extended Double -> Vector (Extended Double)+evalLandscapeAll landscape arg =+  if V.null landscape then V.empty+  else (evalLandscape landscape 0 arg) `cons` (evalLandscapeAll (V.tail landscape) arg)++{- |+  Compute a linear combination of the landscapes.+  If the coefficient list is too short, the rest of the coefficients are assumed to be zero.+  If it is too long, the extra coefficients are discarded.+-}+linearComboLandscapes :: [Double] -> [Landscape] -> Landscape+linearComboLandscapes coeffs landscapes =+  let maxlen      = L.maximum $ L.map V.length landscapes+      emptylayer  = V.fromList [(MinusInfty, Finite 0.0), (Infinity, Finite 0.0)]+      landscapes' = L.map (\l -> l V.++ (V.replicate (maxlen - V.length l) emptylayer)) landscapes++      myconcat v1 v2+        | V.null v1 = v2+        | V.null v2 = v1+        | otherwise = ((V.head v1) V.++ (V.head v2)) `cons` (myconcat (V.tail v1) (V.tail v2))++      xs        = L.map (V.map (V.map fst)) landscapes'+      concatted = L.foldl myconcat V.empty xs+      unionXs   = V.map ((quickSort (>)) . V.fromList . L.nub . V.toList) concatted+      yVals     = L.map (\landscape ->+                    mapWithIndex (\i v -> V.map (evalLandscape landscape i) v) unionXs) landscapes'+      yVals'    = L.zipWith (\coeff yvals ->+                    V.map (V.map ((Finite coeff)*)) yvals) coeffs yVals+      finalY    = L.foldl1 (\acc new -> V.zipWith (V.zipWith (+)) acc new) yVals'+  in V.zipWith V.zip unionXs finalY++-- | Average the persistence landscapes.+avgLandscapes :: [Landscape] -> Landscape+avgLandscapes landscapes =+  let numscapes = L.length landscapes+      coeffs    = L.replicate numscapes (1.0/(fromIntegral numscapes))+  in linearComboLandscapes coeffs landscapes++-- | Subtract the second landscape from the first.+diffLandscapes :: Landscape -> Landscape -> Landscape+diffLandscapes scape1 scape2 = linearComboLandscapes [1, -1] [scape1, scape2]++{- |+  If p>=1 then it will compute the L^p norm on the given interval.+  Uses trapezoidal approximation.+  You should ensure that the stepsize partitions the interval evenly.+-}+normLp :: Extended Double -- ^p, the power of the norm+       -> (Double, Double) -- ^Interval to compute the integral over+       -> Double -- ^Step size+       -> Landscape -- ^Persistence landscape whose norm is to be computed+       -> Maybe Double+normLp p interval step landscape =+  let len = V.length landscape+      a   = fst interval+      b   = snd interval++      fcn x =+        let vals = V.map (\n ->+                     abs $ unExtend $ evalLandscape landscape n (Finite x)) $ 0 `range` (len - 1)+        in+          case p of+            Infinity    -> V.maximum vals+            Finite 1.0  -> V.sum vals+            Finite 2.0  -> sqrt $ V.sum $ V.map (\a -> a*a) vals+            Finite p'   -> (**(1.0/p')) $ V.sum $ V.map (**p') vals++      computeSum :: Double -> Double -> Double+      computeSum currentX result =+        let nextX = currentX + step+        in+          if nextX > b then result + (fcn nextX)+          else computeSum nextX (result + 2.0*(fcn nextX))++  in+    if p < (Finite 1.0) then Nothing+    else Just $ 0.5*step*(computeSum a $ fcn a)++{- |+  Given the same information as above, computes the L^p distance between the two landscapes.+  One way to compare the topologies of two filtrations.+-}+metricLp :: Extended Double -- ^p, power of the metric+         -> (Double, Double) -- ^Interval on which the integral will be computed+         -> Double -- ^Step size+         -> Landscape -- ^First landscape+         -> Landscape -- ^Second landscape+         -> Maybe Double+metricLp p interval step scape1 scape2 = normLp p interval step $ diffLandscapes scape1 scape2

@@ -0,0 +1,203 @@+{- |+Module     : Persistence.HasseDiagram+Copyright  : (c) Eben Kadile, 2018+License    : BSD 3 Clause+Maintainer : eben.cowley42@gmail.com+Stability  : experimental++This module allows one to do computations involving directed graphs. Namely, it allows you to convert a directed graph (presented in a generic way) into a simplicial complex and provides functions for constructing the "directed clique complex," see below.++This module uses 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 of arrays representing each 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 and 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 are groups of nodes that simply take input, process it, 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.++The algorithm for constructing the directed clique complex comes from this paper by Markram et al: https://www.frontiersin.org/articles/10.3389/fncom.2017.00048/full.++-}++module Persistence.HasseDiagram+  ( Node+  , HasseDiagram+  , hsd2String+  , dGraph2sc+  , encodeDirectedGraph+  , directedFlagComplex+  , hDiagram2sc+  ) where++import Persistence.Util+import Persistence.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,+  construct a 1-dimensional simplicial complex in the canonical way.+  Betti numbers of this simplicial complex can be used to count cycles and connected components.+-}+dGraph2sc :: Int -> [(Int, Int)] -> SimplicialComplex+dGraph2sc v edges =+  (v, V.fromList [V.fromList $ L.map (\(i, j) -> (i `cons` (j `cons` V.empty), V.empty)) edges])++{- |+  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.+Algorithm adapted from the one shown in the supplementary materials of this paper: https://www.frontiersin.org/articles/10.3389/fncom.2017.00048/full+-}+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+                  --the vertices of the new simplex are+                  --the vertices of the old simplex plus the sink+                  verts    = sink `cons` (one oldNode)+                  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 "HasseDiagram.directedFlagComplex.makeDiagram.makeNode.testTargets. This is a bug. Please email the Persistence maintainers."++              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.+hDiagram2sc :: HasseDiagram -> SimplicialComplex+hDiagram2sc 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)

@@ -0,0 +1,932 @@+{- |+Module     : Persistence.Matrix+Copyright  : (c) Eben Kadile, 2018+License    : BSD 3 Clause+Maintainer : eben.cowley42@gmail.com+Stability  : experimental++This module contains a variety of matrix utility functions, used in the computation of Betti numbers and simplicial homology groups.++Most importantly, it includes functions for computing the rank, normal form, and kernel of matrices. For the computation of homology groups and Betti numbers, one must perform column operations on one matrix to get it into column echelon form and find its kernel while also performing the inverse row operations on the next matrix to be operated on.++Bool is an instance of Num here (instance given in Util) so that functions can be somewhat generalized to act on both integers and integers modulo 2.++-}++module Persistence.Matrix (+  -- * Types+    IMatrix+  , BMatrix+  -- * Utilities+  , getDiagonal+  , getUnsignedDiagonal+  , transposeMat+  , transposePar+  , multiply+  , multiplyPar+  -- Int matrices+  , rankInt+  , rankIntPar+  , normalFormInt+  , normalFormIntPar+  , kernelInt+  , kernelIntPar+  , imgInKerInt+  , imgInKerIntPar+  -- Bool matrices+  , rankBool+  , kernelBool+  , 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 Persistence.Util++import Data.List as L+import Data.Vector as V++import Control.Parallel.Strategies++-- * Types++-- | Matrix of integers.+type IMatrix = Vector (Vector Int)++-- | Matrix of integers modulo 2. Alternatively, matrix over the field with 2 elements.+type BMatrix = Vector (Vector Bool)++-- * Utilities++isMatrix :: Vector (Vector a) -> Bool+isMatrix mat =+  let rowLen = V.length $ V.head mat+  in V.all (\r -> V.length r == rowLen) mat++-- | Take the transpose a matrix (no fancy optimizations, yet).+transposeMat :: Vector (Vector a) -> Maybe (Vector (Vector a))+transposeMat mat =+  if isMatrix mat+  then Just $ V.map (\i -> V.map (\row -> row ! i) mat) $ 0 `range` ((V.length $ V.head mat) - 1)+  else Nothing++-- | Take the transpose of a matrix using parallel evaluation of rows.+transposePar :: Vector (Vector a) -> Vector (Vector a)+transposePar mat =+  parMapVec (\i -> V.map (\row -> row ! i) mat) $ 0 `range` ((V.length $ V.head mat) - 1)++-- | Multiply two matrices+multiply :: Num a => Vector (Vector a) -> Vector (Vector a) -> Vector (Vector a)+multiply mat1 mat2 =+  let t =+        case transposeMat mat2 of+          Just m  -> m+          Nothing -> error "error in multiply"+  in V.map (\row -> V.map (dotProduct row) t) mat1++-- | Multiply matrices, evaluate rows in parallel if processors are available+multiplyPar :: Num a => Vector (Vector a) -> Vector (Vector a) -> Vector (Vector a)+multiplyPar mat1 mat2 = runEval $ do+  let t =+        case transposeMat mat2 of+          Just m  -> m+          Nothing -> error "error in multiplyPar"+  rseq t+  return $ parMapVec (\row -> V.map (dotProduct row) t) mat1++-- | Get the diagonal elements.+getDiagonal :: Vector (Vector a) -> [a]+getDiagonal matrix =+  if V.null matrix then []+  else L.map (\i -> matrix ! i ! i) [0..(min (V.length matrix) (V.length $ V.head matrix)) - 1]++-- | Get the absolute value of each of the diagonal elements in a list.+getUnsignedDiagonal :: Num a => Vector (Vector a) -> [a]+getUnsignedDiagonal matrix =+  if V.null matrix then []+  else L.map (\i -> abs $ matrix ! i ! i) [0..(min (V.length matrix) (V.length $ V.head matrix)) - 1]++--assumes index1 < index2+colOperation :: Int -> Int -> (Int, Int, Int, Int) -> IMatrix -> IMatrix+colOperation index1 index2 (c11, c12, c21, c22) matrix =+  let calc row =+        let elem1  = row ! index1+            elem2  = row ! index2+            first  = V.take index1 row+            second = V.drop (index1 + 1) (V.take index2 row)+            third  = V.drop (index2 + 1) row+        in first V.++ (cons (c11*elem1 + c12*elem2) second) V.++ (cons (c22*elem2 + c21*elem1) third)+  in V.map calc matrix++colOperationPar :: Int -> Int -> (Int, Int, Int, Int) -> IMatrix -> IMatrix+colOperationPar index1 index2 (c11, c12, c21, c22) matrix =+  let calc row =+        let elem1  = row ! index1+            elem2  = row ! index2+            first  = V.take index1 row+            second = V.drop (index1 + 1) (V.take index2 row)+            third  = V.drop (index2 + 1) row+        in first V.++ (cons (c11*elem1 + c12*elem2) second)+             V.++ (cons (c22*elem2 + c21*elem1) third)+  in parMapVec calc matrix++--assumes index1 < index2+rowOperation :: Int -> Int -> (Int, Int, Int, Int) -> IMatrix -> IMatrix+rowOperation index1 index2 (c11, c12, c21, c22) matrix =+  let row1   = matrix ! index1+      row2 = matrix ! index2+      first  = V.take index1 matrix+      second = V.drop (index1 + 1) $ V.take index2 matrix+      third  = V.drop (index2 + 1) matrix+  in first V.++ (cons ((c11 `mul` row1) `add` (c12 `mul` row2)) second)+    V.++ (cons ((c22 `mul` row2) `add` (c21 `mul` row1)) third)++rowOperationPar :: Int -> Int -> (Int, Int, Int, Int) -> IMatrix -> IMatrix+rowOperationPar index1 index2 (c11, c12, c21, c22) matrix =+  let row1   = matrix ! index1+      row2 = matrix ! index2+      first  = V.take index1 matrix+      second = V.drop (index1 + 1) (V.take index2 matrix)+      third  = V.drop (index2 + 1) matrix+  in runEval $ do+     a <- rpar $ (c11 `mul` row1) `add` (c12 `mul` row2)+     b <- rpar $ (c21 `mul` row1) `add` (c22 `mul` row2)+     rseq (a,b)+     return $ first V.++ (a `cons` second) V.++ (b `cons` third)++-- * Int matrices++--RANK--------------------------------------------------------------------++--finds the pivot in a given row for Gaussian elimination given the index of the pivot row and the matrix+--returns whether or not the row needs to be eliminated with the rearranged matrix and+--the column switch performed (if there was one), returns Nothing if the row is all zeroes+chooseGaussPivotInt :: (Int, Int) -> IMatrix -> Maybe (Bool, IMatrix, Maybe (Int, Int))+chooseGaussPivotInt (rowIndex, colIndex) mat =+  let row     = mat ! rowIndex+      indices = V.filter (\index -> index > colIndex) $ V.findIndices (\x -> x /= 0) row+  in+    if row ! colIndex == 0 then+      if V.null indices then Nothing+      else+        let j = V.head indices+        in+          if row ! j == 0 then error "Persistence.Matrix.chooseGaussPivotInt. This is a bug. Please email the persistence maintainers."+          else Just (V.length indices > 1, V.map (switchElems colIndex j) mat, Just (colIndex, j))+    else Just (V.length indices > 0, mat, Nothing)++--does gaussian elimination on the pivot row of an integer matrix+improveRowInt :: (Int, Int) -> Int -> IMatrix -> IMatrix+improveRowInt (rowIndex, colIndex) numCols matrix =+  let improve i mat =+        if i == numCols then mat+        else+          let row   = mat ! rowIndex+              pivot = row ! colIndex+              x     = row ! i+              next  = i + 1+          in+            --boundary operators have lots of zeroes+            --better to catch that instead of doing unnecessary %+            if pivot == 0 then+              if V.all (\a -> a == 0) row then mat+              else error "Persistence.Matrix.improveRowInt. This is a bug. Please email the Persistence maintainers."+            else+              if x == 0 || (x `mod` pivot == 0) then+                improve next mat+              else+                let gcdTriple = extEucAlg pivot x+                    gcd       = one gcdTriple+                in improve next $ colOperation colIndex i (thr gcdTriple, two gcdTriple, x `div` gcd, -(pivot `div` gcd)) mat++  in improve (colIndex + 1) matrix++--given pivot index and pivot paired with matrix whose pivot row has been improved+--eliminates the entries in the pivot row+--the kinds of matrices that the functions work on will have lots of zeroes+--better to catch that with a condition than perform an unnecessary division+elimRowInt :: (Int, Int) -> IMatrix -> IMatrix+elimRowInt (rowIndex, colIndex) elems =+  let pCol  = V.map (\row -> row ! colIndex) elems+      pivot = pCol ! rowIndex+      c1    = colIndex + 1++      makeCoeffs i v =+        if V.null v then empty+        else let x = V.head v; xs = V.tail v in+          if x == 0 then makeCoeffs (i + 1) xs+          else (i, x `div` pivot) `cons` (makeCoeffs (i + 1) xs)++      calc :: IMatrix -> Vector (Int, Int) -> IMatrix+      calc mat ops =+        if V.null ops then mat+        else+          let (i, coeff) = V.head ops+          in calc (mapWithIndex (\j row ->+               replaceElem i ((row ! i) - coeff*(pCol ! j)) row) mat) (V.tail ops)++  in+    if pivot == 0 then error "Persistence.Matrix.elimRowInt. This is a bug. Please email the Persistence maintainers."+    else calc elems $ makeCoeffs c1 $ V.drop c1 $ elems ! rowIndex++-- | Finds the rank of integer matrix (number of linearly independent columns).+rankInt :: IMatrix -> Int+rankInt matrix =+  let rows     = V.length matrix+      cols     = V.length $ V.head matrix+      cols1    = cols - 1++      doColOps (rowIndex, colIndex) mat =+        if rowIndex == rows || colIndex == cols then mat else+          case chooseGaussPivotInt (rowIndex, colIndex) mat of+            Just (True, mx, _)  ->+              doColOps (rowIndex + 1, colIndex + 1)+                $ elimRowInt (rowIndex, colIndex) $ improveRowInt (rowIndex, colIndex) cols mx+            Just (False, mx, _) -> doColOps (rowIndex + 1, colIndex + 1) mx+            Nothing             -> doColOps (rowIndex + 1, colIndex) mat++      countNonZeroCols mat =+        V.sum $ V.map (\i -> if V.any (\j -> mat ! j ! i /= 0)+          $ 0 `range` (rows - 1) then 1 else 0) $ 0 `range` cols1++  in countNonZeroCols $ doColOps (0, 0) matrix++--does gaussian elimination on the pivot row of an integer matrix in parallel+improveRowIntPar :: (Int, Int) -> Int -> IMatrix -> IMatrix+improveRowIntPar (rowIndex, colIndex) numCols matrix =+  let improve i mat =+        if i == numCols then mat else+          let row   = mat ! rowIndex+              pivot = row ! colIndex+              x     = row ! i+              next  = i + 1+          in+            --boundary operators have lots of zeroes+            --better to catch that instead of doing unnecessary %+            if x == 0 || (x `mod` pivot == 0) then+              improve next mat+            else+              let gcdTriple = extEucAlg pivot x+                  gcd       = one gcdTriple+              in improve next+                   $ colOperationPar colIndex i+                     (thr gcdTriple, two gcdTriple, x `div` gcd, -(pivot `div` gcd)) mat+  in improve (colIndex + 1) matrix++--eliminates a row in parallel+elimRowIntPar :: (Int, Int) -> IMatrix -> IMatrix+elimRowIntPar (rowIndex, colIndex) elems =+  let pCol  = V.map (\row -> row ! colIndex) elems+      pivot = pCol ! rowIndex+      c1    = colIndex + 1++      makeCoeffs i v =+        if V.null v then empty+        else let x = V.head v; xs = V.tail v in+          if x == 0 then makeCoeffs (i + 1) xs+          else (i, x `div` pivot) `cons` (makeCoeffs (i + 1) xs)++      calc :: IMatrix -> Vector (Int, Int) -> IMatrix+      calc mat ops =+        if V.null ops then mat+        else+          let (i, coeff) = V.head ops+          in+            calc (parMapWithIndex+              (\j row -> replaceElem i ((row ! i) - coeff*(pCol ! j)) row) mat) (V.tail ops)++  in+    if pivot == 0+    then error "Persistence.Matrix.elimRowIntPar. This is a bug. Please email the Persistence maintainers."+    else calc elems $ makeCoeffs c1 $ V.drop c1 $ elems ! rowIndex++-- | Calculates the rank of a matrix by operating on multiple rows in parallel.+rankIntPar :: IMatrix -> Int+rankIntPar matrix =+  let rows     = V.length matrix+      cols     = V.length $ V.head matrix+      cols1    = cols - 1++      doColOps (rowIndex, colIndex) mat =+        if rowIndex == rows || colIndex == cols then mat else+          case chooseGaussPivotInt (rowIndex, colIndex) mat of+            Just (True, mx, _)  ->+              doColOps (rowIndex + 1, colIndex + 1) $ elimRowIntPar (rowIndex, colIndex) $+                improveRowIntPar (rowIndex, colIndex) cols mx+            Just (False, mx, _) -> doColOps (rowIndex + 1, colIndex + 1) mx+            Nothing             -> doColOps (rowIndex + 1, colIndex) mat++      countNonZeroCols mat =+        V.sum $ parMapVec (\i -> if V.any (\j -> mat ! j ! i /= 0) $+          0 `range` (rows - 1) then 1 else 0) $ 0 `range` cols1++  in countNonZeroCols $ doColOps (0, 0) matrix++--NORMAL FORM-------------------------------------------------------------++--rearranges matrix so that the pivot entry is in the correct position+--returns true if more elimination is necessary+--returns Nothing if there is nothing but zeroes after the current pivot position+chooseRowPivotInt :: (Int, Int) -> Int -> Int -> IMatrix -> Maybe (Bool, IMatrix)+chooseRowPivotInt (rowIndex, colIndex) numRows numCols mat =+  let row      = mat ! rowIndex+      rIndices = V.toList $ V.findIndices (\x -> x /= 0) row+  in+    if 0 == row ! colIndex then+      case rIndices of+        (i:is)  -> Just ((L.length is) > 0, V.map (switchElems i colIndex) mat)+        []      ->+          case V.toList $ V.findIndices (\x -> x /= 0) $ V.map (\r -> r ! colIndex) mat of+            (i:_)  -> Just (True, switchElems i rowIndex mat)+            []     -> Nothing+    else Just ((L.length rIndices) > 1, mat)++--given pivot index and pivot paired with matrix, improves pivot column with row operations+improveColInt :: Int -> Int -> IMatrix -> IMatrix+improveColInt pIndex maxIndex matrix =+  let improve i mat =+        if i == maxIndex then mat else+          let pivot = matrix ! pIndex ! pIndex+              x     = matrix ! i ! pIndex+              next  = i + 1+          in+            --boundary operators have lots of zeroes+            --better to catch that instead of doing unnecessary %+            if x == 0 || (x `mod` pivot == 0) then+              improve next mat+            else+              let gcdTriple = extEucAlg pivot x+                  gcd       = one gcdTriple+              in improve next+                   $ rowOperation pIndex i+                     (thr gcdTriple, two gcdTriple, x `div` gcd, -(pivot `div` gcd)) mat+  in improve (pIndex + 1) matrix++--eliminates the pivot column of a matrix to obtain normal form+elimColInt :: (Int, Int) -> IMatrix -> IMatrix+elimColInt (rowIndex, colIndex) elems =+  let pRow  = elems ! rowIndex+      pivot = pRow ! colIndex+      ri1   = rowIndex + 1+      makeCoeffs i v =+        if V.null v then empty+        else let x = V.head v; xs = V.tail v in+          if x == 0 then makeCoeffs (i + 1) xs+          else (i, x `div` pivot) `cons` (makeCoeffs (i + 1) xs)+      calc :: IMatrix -> Vector (Int, Int) -> IMatrix+      calc mat ops =+        if V.null ops then mat+        else let (i, coeff) = V.head ops in+          calc (replaceElem i ((mat ! i) `subtr` (coeff `mul` pRow)) mat) (V.tail ops)+  in+    if pivot == 0+    then error "Persistence.Matrix.elimColInt. This is a bug. Please email the Persistence maintainters."+    else calc elems $ makeCoeffs ri1 $ V.drop ri1 $ V.map (\row -> row ! colIndex) elems++finish :: Int -> IMatrix -> IMatrix+finish diagLen matrix =+  let calc i mat =+        let i1    = i + 1+            row   = mat ! i+            entry = row ! i+            nextR = mat ! i1+            nextE = nextR ! i1+        in+          if entry == 0 || i1 == diagLen then mat+          else if entry `divides` nextE then calc i1 mat+          else+            let mat'      = replaceElem i (replaceElem i1 nextE row) mat+                gcdTriple = extEucAlg entry nextE; gcd = one gcdTriple+                improve   =+                  colOperation i i1+                    (thr gcdTriple, two gcdTriple, nextE `div` gcd, entry `div` gcd)+                cleanup   = \m -> elimColInt (i, i) $ elimRowInt (i, i) m+            in calc i1 $ cleanup $ improve mat'++      absDiag i mat =+        if i == V.length mat+        then mat+        else absDiag (i + 1)+               $ replaceElem i (replaceElem i (abs $ mat!i!i) $ mat!i) mat++      filtered = V.partition (\row -> V.any (\x -> x /= 0) row) matrix+  in absDiag 0 $ calc 0 $ (fst filtered) V.++ (snd filtered)++-- | Get the Smith normal form of an integer matrix.+normalFormInt :: IMatrix -> IMatrix+normalFormInt matrix =+  let rows = V.length matrix+      cols = V.length $ V.head matrix+      diag = min rows cols++      calc (rowIndex, colIndex) mat =+        if rowIndex == rows || colIndex == cols then mat+        else+          case chooseRowPivotInt (rowIndex, colIndex) rows cols mat of+            Just (True, mx)  ->+              calc (rowIndex + 1, colIndex + 1) $+                elimColInt (rowIndex, colIndex) $ improveColInt rowIndex rows $+                  elimRowInt (rowIndex, colIndex) $ improveRowInt (rowIndex, colIndex) cols mx+            Just (False, mx) ->+              calc (rowIndex + 1, colIndex + 1) $+                elimColInt (rowIndex, colIndex) $ improveColInt rowIndex rows mx+            Nothing          -> calc (rowIndex + 1, colIndex) mat++  in+    if V.null matrix+    then empty+    else finish diag $ calc (0, 0) matrix++--improves the pivot column of a matrix in parallel+improveColIntPar :: Int -> Int -> IMatrix -> IMatrix+improveColIntPar pIndex maxIndex matrix =+  let improve i mat =+        if i == maxIndex then mat else+          let col   = V.map (\row -> row ! pIndex) mat+              pivot = col ! pIndex+              x     = col ! i+              next  = i + 1+          in+            --boundary operators have lots of zeroes+            --better to catch that instead of doing unnecessary %+            if x == 0 || (x `mod` pivot == 0) then+              improve next mat+            else+              let gcdTriple = extEucAlg pivot x+                  gcd       = one gcdTriple+              in improve next+                   $ rowOperationPar pIndex i+                     (thr gcdTriple, two gcdTriple, x `div` gcd, -(pivot `div` gcd)) mat+  in improve (pIndex + 1) matrix++--NEEDS TO BE PARALLELIZED+--eliminates pivot column in parallel+elimColIntPar :: (Int, Int) -> IMatrix -> IMatrix+elimColIntPar (rowIndex, colIndex) elems =+  let pRow  = elems ! rowIndex+      pivot = pRow ! colIndex+      ri1   = rowIndex + 1++      makeCoeffs i v =+        if V.null v then empty+        else let x = V.head v; xs = V.tail v in+          if x == 0 then makeCoeffs (i + 1) xs+          else (i, x `div` pivot) `cons` (makeCoeffs (i + 1) xs)++      calc :: IMatrix -> Vector (Int, Int) -> IMatrix+      calc mat ops =+        if V.null ops then mat+        else let (i, coeff) = V.head ops in+          calc (replaceElem i ((mat ! i) `subtr` (coeff `mul` pRow)) mat) (V.tail ops)++  in calc elems $ makeCoeffs ri1 $ V.drop ri1 $ V.map (\row -> row ! colIndex) elems++-- | Gets the Smith normal form of a matrix, uses lots of parallelism if processors are available.+normalFormIntPar :: IMatrix -> IMatrix+normalFormIntPar matrix =+  let rows = V.length matrix+      cols = V.length $ V.head matrix+      diag = min rows cols++      calc (rowIndex, colIndex) mat =+        if rowIndex == rows || colIndex == cols then mat+        else+          case chooseRowPivotInt (rowIndex, colIndex) rows cols mat of+            Just (True, mx)  ->+              calc (rowIndex + 1, colIndex + 1) $+                elimColIntPar (rowIndex, colIndex) $ improveColIntPar rowIndex rows $+                  elimRowIntPar (rowIndex, colIndex) $ improveRowIntPar (rowIndex, colIndex) cols mx+            Just (False, mx) ->+              calc (rowIndex + 1, colIndex + 1) $+                elimColIntPar (rowIndex, colIndex) $ improveColIntPar rowIndex rows mx+            Nothing          -> calc (rowIndex + 1, colIndex) mat++  in+    if V.null matrix+    then empty+    else finish diag $ calc (0, 0) matrix++--KERNEL------------------------------------------------------------------++--improves the pivot row of an integer matrix+--performs the same column operations on the identity+improveRowIntWithId :: (Int, Int) -> Int -> IMatrix -> IMatrix -> (IMatrix, Int, IMatrix)+improveRowIntWithId (rowIndex, colIndex) numCols elems identity =+  let improve i mat ide =+        if i == numCols then (mat, mat ! rowIndex ! colIndex, ide) else+          let row   = mat ! rowIndex+              pivot = row ! colIndex+              x     = row ! i+              next  = i + 1+          in+            if x == 0 || (x `mod` pivot == 0) then+              improve next mat ide+            else+            let gcdTriple = extEucAlg pivot x+                gcd       = one gcdTriple+                transform =+                  colOperation colIndex i+                    (thr gcdTriple, two gcdTriple, x `div` gcd, -(pivot `div` gcd))+            in improve next (transform mat) (transform ide)++  in improve (colIndex + 1) elems identity++--eliminates all the entries in the pivot row that come after the pivot+--after the matrix has been improved+--returns the new matrix (fst) and+--the identity with whatever column operations were performed (snd)+elimRowIntWithId :: (Int, Int) -> Int -> (IMatrix, Int, IMatrix) -> (IMatrix, IMatrix)+elimRowIntWithId (rowIndex, colIndex) numCols (elems, pivot, identity) =+  let row            = elems ! rowIndex+      elim i mat ide =+        if i == numCols then (mat, ide)+        else+          let coeff     = (row ! i) `div` pivot+              transform =+                V.map (\r -> (V.take i r)+                  V.++ (cons ((r ! i) - coeff*(r ! colIndex)) (V.drop (i + 1) r)))+          in elim (i + 1) (transform mat) (transform ide)+  in elim (colIndex + 1) elems identity++-- | Finds a basis for the kernel of a matrix, arranges the basis vectors into the rows of a matrix.+kernelInt :: IMatrix -> IMatrix+kernelInt matrix =+  let rows     = V.length matrix+      cols     = V.length $ V.head matrix+      cols1    = cols - 1+      identity = V.map (\i -> (V.replicate i 0)+                   V.++ (cons 1 (V.replicate (cols1 - i) 0))) $ 0 `range` cols1++      doColOps (rowIndex, colIndex) (elems, ide) =+        if rowIndex == rows || colIndex == cols then (elems, ide) else+          case chooseGaussPivotInt (rowIndex, colIndex) elems of+            Just (True, mx, Just (i, j))  ->+              doColOps (rowIndex + 1, colIndex + 1) $ elimRowIntWithId (rowIndex, colIndex) cols $+                improveRowIntWithId (rowIndex, colIndex) cols mx $ V.map (switchElems i j) ide+            Just (True, _, Nothing)       ->+              doColOps (rowIndex + 1, colIndex + 1) $ elimRowIntWithId (rowIndex, colIndex) cols $+                improveRowIntWithId (rowIndex, colIndex) cols elems ide+            Just (False, mx, Just (i, j)) ->+              doColOps (rowIndex + 1, colIndex + 1) (mx, V.map (switchElems i j) ide)+            Just (False, _, _)            -> doColOps (rowIndex + 1, colIndex + 1) (elems, ide)+            Nothing                       -> doColOps (rowIndex + 1, colIndex) (elems, ide)++      result   = doColOps (0, 0) (matrix, identity)+      elems    = fst result+      ide      = snd result+  in V.map (\i -> V.map (\row -> row ! i) ide)+       $ V.filter (\i -> V.all (\row -> row ! i == 0) elems) $ 0 `range` cols1++--improves row in parallel and does the same thing to the identity matrix in parallel+improveRowIntWithIdPar :: (Int, Int) -> Int -> IMatrix -> IMatrix -> (IMatrix, Int, IMatrix)+improveRowIntWithIdPar (rowIndex, colIndex) numCols elems identity =+  let improve i mat ide =+        if i == numCols then (mat, mat ! rowIndex ! colIndex, ide) else+          let row   = mat ! rowIndex+              pivot = row ! colIndex+              x     = row ! i+              next  = i + 1+          in+            --boundary operators have lots of zeroes+            --better to catch that instead of doing unnecessary %+            if x == 0 || (x `mod` pivot == 0) then+              improve next mat ide+            else+            let gcdTriple = extEucAlg pivot x+                gcd       = one gcdTriple+                transform = colOperationPar colIndex i+                              (thr gcdTriple, two gcdTriple, x `div` gcd, -(pivot `div` gcd))+            in improve next (transform mat) (transform ide)+  in improve (colIndex + 1) elems identity++--eliminates entries in the pivot row in parallel and does the same to the identity+elimRowIntWithIdPar :: (Int, Int) -> Int -> (IMatrix, Int, IMatrix) -> (IMatrix, IMatrix)+elimRowIntWithIdPar (rowIndex, colIndex) numCols (elems, pivot, identity) =+  let row            = elems ! rowIndex+      elim i mat ide =+        if i == numCols then (mat, ide)+        else+          let coeff     = (row ! i) `div` pivot+              transform = parMapVec (\r -> (V.take i r)+                            V.++ (cons ((r ! i) - coeff*(r ! colIndex)) (V.drop (i + 1) r)))+          in elim (i + 1) (transform mat) (transform ide)+  in elim (colIndex + 1) elems identity++-- | Computes basis vectors for the kernel of an integer matrix and arranges them into the rows of a matrix using lots of parallelism if processors are available.+kernelIntPar :: IMatrix -> IMatrix+kernelIntPar matrix =+  let rows     = V.length matrix+      cols     = V.length $ V.head matrix+      cols1    = cols - 1+      identity = V.map (\i -> (V.replicate i 0)+                   V.++ (cons 1 (V.replicate (cols1 - i) 0))) $ 0 `range` cols1++      doColOps (rowIndex, colIndex) (elems, ide) =+        if rowIndex == rows || colIndex == cols then (elems, ide) else+          case chooseGaussPivotInt (rowIndex, colIndex) elems of+            Just (True, mx, Just (i, j))  ->+              doColOps (rowIndex + 1, colIndex + 1)+                $ elimRowIntWithIdPar (rowIndex, colIndex) cols+                  $ improveRowIntWithIdPar (rowIndex, colIndex) cols mx+                    $ V.map (switchElems i j) ide+            Just (True, _, Nothing)       ->+              doColOps (rowIndex + 1, colIndex + 1)+                $ elimRowIntWithIdPar (rowIndex, colIndex) cols+                  $ improveRowIntWithIdPar (rowIndex, colIndex) cols elems ide+            Just (False, mx, Just (i, j)) ->+              doColOps (rowIndex + 1, colIndex + 1) (mx, V.map (switchElems i j) ide)+            Just (False, _, _)            -> doColOps (rowIndex + 1, colIndex + 1) (elems, ide)+            Nothing                       -> doColOps (rowIndex + 1, colIndex) (elems, ide)++      result   = doColOps (0, 0) (matrix, identity)+      elems    = fst result+      ide      = snd result+  in V.map (\i -> V.map (\row -> row ! i) ide)+       $ V.filter (\i -> V.all (\row -> row ! i == 0) elems) $ 0 `range` cols1++--FIND IMAGE IN BASIS OF KERNEL-------------------------------------------++--improves the row of the first matrix+--and performs the inverse column operations on the second matrix+improveRowIntWithInv :: (Int, Int) -> Int -> IMatrix -> IMatrix -> (IMatrix, Int, IMatrix)+improveRowIntWithInv (rowIndex, colIndex) numCols kernel image =+  let improve i ker img =+        if i == numCols then (ker, ker ! rowIndex ! colIndex, img) else+          let row   = ker ! rowIndex+              pivot = row ! colIndex+              x     = row ! i+              next  = i + 1+          in+            if x == 0 || (x `mod` pivot == 0) then+              improve next ker img+            else+              let gcdTriple  = extEucAlg pivot x+                  gcd        = one gcdTriple+                  q1         = pivot `div` gcd+                  q2         = x `div` gcd+                  transform1 =+                    colOperationPar colIndex i (thr gcdTriple, two gcdTriple, q2, -q1)+                  transform2 =+                    rowOperationPar colIndex i (-q1, -(two gcdTriple), -q2, thr gcdTriple)+              in improve next (transform1 ker) (transform2 img)+  in improve (colIndex + 1) kernel image++--eliminates the row if the first matrix+--performs inverse column operations on the second matrix+elimRowIntWithInv :: (Int, Int) -> Int -> (IMatrix, Int, IMatrix) -> (IMatrix, IMatrix)+elimRowIntWithInv (rowIndex, colIndex) numCols (kernel, pivot, image) =+  let elim i ker img+        | i == numCols            = (ker, img)+        | row ! i == 0 = elim (i + 1) ker img+        | otherwise               =+          let coeff      = (row ! i) `div` pivot+              transform1 = V.map (\r -> replaceElem i ((r ! i) - coeff*(r ! colIndex)) r)+              transform2 = \mat -> replaceElem colIndex+                             ((coeff `mul` (mat ! i)) `add` (mat ! colIndex)) mat+          in elim (i + 1) (transform1 ker) (transform2 img)+        where row = ker ! rowIndex+  in elim (colIndex + 1) kernel image++-- | Calculates the image of the second matrix represented in the basis of the kernel of the first matrix.+imgInKerInt :: IMatrix -> IMatrix -> IMatrix+imgInKerInt toColEsch toImage =+  let rows     = V.length toColEsch+      cols     = V.length $ V.head toColEsch++      doColOps (rowIndex, colIndex) (ker, img) =+        if rowIndex == rows || colIndex == cols then (ker, img)+        else case chooseGaussPivotInt (rowIndex, colIndex) ker of+          Just (True, _, Nothing)       ->+            doColOps (rowIndex + 1, colIndex + 1) $ elimRowIntWithInv (rowIndex, colIndex) cols $+              improveRowIntWithInv (rowIndex, colIndex) cols ker img+          Just (True, mx, Just (i, j))  ->+            doColOps (rowIndex + 1, colIndex + 1) $ elimRowIntWithInv (rowIndex, colIndex) cols $+              improveRowIntWithInv (rowIndex, colIndex) cols mx $ switchElems i j img+          Just (False, mx, Just (i, j)) ->+            doColOps (rowIndex + 1, colIndex + 1) (mx, switchElems i j img)+          Just (False, _, _)            -> doColOps (rowIndex + 1, colIndex + 1) (ker, img)+          Nothing                       -> doColOps (rowIndex + 1, colIndex) (ker, img)++      result  = doColOps (0, 0) (toColEsch, toImage)+      colEsch = fst result+      image   = snd result+  in V.map (\i -> image ! i)+       $ V.filter (\i -> V.all (\row -> row ! i == 0) colEsch) $ 0 `range` (cols - 1)++--improves row and does inverse operations in parallel+improveRowIntWithInvPar :: (Int, Int) -> Int -> IMatrix -> IMatrix -> (IMatrix, Int, IMatrix)+improveRowIntWithInvPar (rowIndex, colIndex) numCols kernel image =+  let improve i ker img =+        if i == numCols then (ker, ker ! rowIndex ! colIndex, img) else+          let row   = ker ! rowIndex+              pivot = row ! colIndex+              x     = row ! i+              next  = i + 1+          in+            if x == 0 || (x `mod` pivot == 0) then+              improve next ker img+            else+              let gcdTriple  = extEucAlg pivot x+                  gcd        = one gcdTriple+                  q1         = pivot `div` gcd+                  q2         = x `div` gcd+                  transform1 =+                    colOperationPar colIndex i (thr gcdTriple, two gcdTriple, q2, -q1)+                  transform2 =+                    rowOperationPar colIndex i (-q1, -(two gcdTriple), -q2, thr gcdTriple)+              in improve next (transform1 ker) (transform2 img)+  in improve (colIndex + 1) kernel image++--eliminates row in parallel+--INVERSE OPERATIONS NEED TO BE PARALLELIZED+elimRowIntWithInvPar :: (Int, Int) -> Int -> (IMatrix, Int, IMatrix) -> (IMatrix, IMatrix)+elimRowIntWithInvPar (rowIndex, colIndex) numCols (kernel, pivot, image) =+  let elim i ker img+        | i == numCols            = (ker, img)+        | row ! i == 0 = elim (i + 1) ker img+        | otherwise               =+          let coeff      = (row ! i) `div` pivot+              transform1 = parMapVec (\r -> replaceElem i ((r ! i) - coeff*(r ! colIndex)) r)+              transform2 = \mat -> replaceElem colIndex+                             ((coeff `mul` (mat ! i)) `add` (mat ! colIndex)) mat+          in elim (i + 1) (transform1 ker) (transform2 img)+        where row = ker ! rowIndex+  in elim (colIndex + 1) kernel image++-- | Calculates the image of the second matrix represented in the basis of the kernel of the first matrix. Uses lots of parallelism if processors are available.+imgInKerIntPar :: IMatrix -> IMatrix -> IMatrix+imgInKerIntPar toColEsch toImage =+  let rows     = V.length toColEsch+      cols     = V.length $ V.head toColEsch++      doColOps (rowIndex, colIndex) (ker, img) =+        if rowIndex == rows || colIndex == cols then (ker, img)+        else case chooseGaussPivotInt (rowIndex, colIndex) ker of+          Just (True, _, Nothing)       ->+            doColOps (rowIndex + 1, colIndex + 1) $ elimRowIntWithInvPar (rowIndex, colIndex) cols $+              improveRowIntWithInvPar (rowIndex, colIndex) cols ker img+          Just (True, mx, Just (i, j))  ->+            doColOps (rowIndex + 1, colIndex + 1) $ elimRowIntWithInvPar (rowIndex, colIndex) cols $+              improveRowIntWithInvPar (rowIndex, colIndex) cols mx $ switchElems i j img+          Just (False, mx, Just (i, j)) ->+            doColOps (rowIndex + 1, colIndex + 1) (mx, switchElems i j img)+          Just (False, _, _)            -> doColOps (rowIndex + 1, colIndex + 1) (ker, img)+          Nothing                       -> doColOps (rowIndex + 1, colIndex) (ker, img)++      result = doColOps (0, 0) (toColEsch, toImage)+      ker    = fst result+      img    = snd result+  in V.map (\i -> img ! i)+       $ V.filter (\i -> V.all (\row -> row ! i == 0) ker) $ 0 `range` (cols - 1)++-- * Bool matrices++--RANK--------------------------------------------------------------------++--given the index of the pivot row and the matrix+--determines whether there is a non-zero element in the row, does necessary rearranging+--and returns the column operation that was performed if there was one+--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+      indices = V.filter (\index -> index > colIndex) $ V.findIndices id row+  in+    if not $ row ! colIndex then+      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+elimRowBool :: (Int, Int) -> Int -> BMatrix -> BMatrix+elimRowBool (rowIndex, colIndex) numCols elems =+  let row = elems ! rowIndex+      elim i mat+        | i == numCols  = mat+        | not $ row ! i = elim (i + 1) mat+        | otherwise     = elim (i + 1)+                            $ V.map (\row -> replaceElem i ((row ! i) + (row ! colIndex)) row) mat+  in elim (colIndex + 1) elems++-- | Find the rank of a mod 2 matrix (number of linearly independent columns).+rankBool :: BMatrix -> Int+rankBool matrix =+  let rows  = V.length matrix+      cols  = V.length $ V.head matrix+      cols1 = cols - 1++      doColOps (rowIndex, colIndex) mat =+        if rowIndex == rows || colIndex == cols then mat else+          case chooseGaussPivotBool (rowIndex, colIndex) mat of+            Just (True, mx, _)  ->+              doColOps (rowIndex + 1, colIndex + 1) $ elimRowBool (rowIndex, colIndex) cols mx+            Just (False, mx, _) -> doColOps (rowIndex + 1, colIndex + 1) mat+            Nothing             -> doColOps (rowIndex + 1, colIndex) mat++      countNonZeroCols mat =+        V.sum $ V.map (\i ->+           if V.any (\j -> mat ! j ! i /= 0) (0 `range` (rows - 1)) then 1 else 0) $ 0 `range` cols1+  in countNonZeroCols $ doColOps (0, 0) matrix++--KERNEL------------------------------------------------------------------++--eliminates all the entries in the pivot row that come after the pivot+--after the matrix has been improved+--returns the new matrix (fst) paired+--and the identity with whatever column operations were performed (snd)+elimRowBoolWithId :: (Int, Int) -> Int -> BMatrix -> BMatrix -> (BMatrix, BMatrix)+elimRowBoolWithId (rowIndex, colIndex) numCols elems identity =+  let row = elems ! rowIndex+      elim i mat ide+        | i == numCols  = (mat, ide)+        | not $ row ! i = elim (i + 1) mat ide+        | otherwise     =+          let transform = V.map (\row -> replaceElem i ((row ! i) + (row ! colIndex)) row)+          in elim (i + 1) (transform mat) (transform ide)+  in elim (colIndex + 1) elems identity++-- | Finds the basis of the kernel of a matrix, arranges the basis vectors into the rows of a matrix.+kernelBool :: BMatrix -> BMatrix+kernelBool matrix =+  let rows     = V.length matrix+      cols     = V.length $ V.head matrix+      cols1    = cols - 1+      identity = V.map (\i -> (V.replicate i False)+                   V.++ (cons True (V.replicate (cols1 - i) False))) $ 0 `range` cols1++      doColOps (rowIndex, colIndex) (ker, ide) =+        if rowIndex == rows || colIndex == cols then (ker, ide)+        else+          case chooseGaussPivotBool (rowIndex, colIndex) ker of+            Just (True, _, Nothing)      ->+              doColOps (rowIndex + 1, colIndex + 1) $+                elimRowBoolWithId (rowIndex, colIndex) cols ker ide+            Just (True, mx, Just (i, j)) ->+              doColOps (rowIndex + 1, colIndex + 1) $+                elimRowBoolWithId (rowIndex, colIndex) cols mx $ V.map (switchElems i j) ide+            Just (False, _, Just (i, j)) ->+              doColOps (rowIndex + 1, colIndex + 1) (ker, V.map (switchElems i j) ide)+            Just (False, _, _)           -> doColOps (rowIndex + 1, colIndex + 1) (ker, ide)+            Nothing                      -> doColOps (rowIndex + 1, colIndex) (ker, ide)++      result = doColOps (0, 0) (matrix, identity)+      ker    = fst result+      img    = snd result+  in V.map (\i -> img ! i) $ V.filter (\i -> V.all (\row -> not $ row ! i) ker) $ 0 `range` cols1++--IMAGE IN BASIS OF KERNEL------------------------------------------------++elimRowBoolWithInv :: (Int, Int) -> Int -> BMatrix -> BMatrix -> (BMatrix, BMatrix)+elimRowBoolWithInv (rowIndex, colIndex) numCols toColEch toImage =+  let row = toColEch ! rowIndex+      elim i ech img+          | i == numCols  = (ech, img)+          | not $ row ! i = elim (i + 1) ech img+          | otherwise     =+            let transform1 = V.map (\r -> replaceElem i ((r ! i) + (r ! colIndex)) r)+                transform2 = \mat -> replaceElem colIndex ((mat ! i) `add` (mat ! colIndex)) mat+            in elim (i + 1) (transform1 ech) (transform2 img)+  in elim (colIndex + 1) toColEch toImage++-- | Calculates the image of the second matrix represented in the basis of the kernel of the first matrix.+imgInKerBool :: BMatrix -> BMatrix -> BMatrix+imgInKerBool toColEch toImage =+  let rows  = V.length toColEch+      cols  = V.length $ V.head toColEch+      cols1 = cols - 1++      doColOps (rowIndex, colIndex) (ech, img) =+        if rowIndex == rows || colIndex == cols then (ech, img)+        else+          case chooseGaussPivotBool (rowIndex, colIndex) ech of+            Just (True, _, Nothing)       ->+              doColOps (rowIndex + 1, colIndex + 1) $+                elimRowBoolWithInv (rowIndex, colIndex) cols ech img+            Just (True, mx, Just (i, j))  ->+              doColOps (rowIndex + 1, colIndex + 1) $+                elimRowBoolWithInv (rowIndex, colIndex) cols mx $ switchElems i j img+            Just (False, mx, Just (i, j)) ->+              doColOps (rowIndex + 1, colIndex + 1) (mx, switchElems i j img)+            Just (False, _, _)            -> doColOps (rowIndex + 1, colIndex + 1) (ech, img)+            Nothing                       -> doColOps (rowIndex + 1, colIndex) (ech, img)++      result = doColOps (0, 0) (toColEch, toImage)+      ker    = fst result+      img    = snd result+  in V.map (\i -> img ! i) $ V.filter (\i -> V.all (\row -> not $ row ! i) ker) $ 0 `range` cols1

@@ -0,0 +1,608 @@+{- |+Module     : Persistence.SimplicialComplex+Copyright  : (c) Eben Kadile, 2018+License    : BSD 3 Clause+Maintainer : eben.cowley42@gmail.com+Stability  : experimental++This module provides functions for constructing neighborhood graphs, clique complexes, and Vietoris-Rips complexes, as well as the computation of Betti numbers.++A simplicial complex is like a generalization of a graph, where you can have any number of n-dimensional simplices (line-segments, triangles, tetrahedrons), glued along their boundaries so that the intersection of any two simplices is a simplicial complex. These objects are fundamental to topological data analysis.++An important thing to note about this module, and the library in general, 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 unless you have a lot of RAM. "Light" functions do not use this speed optimization.++The neighborhood graph of point cloud data is a graph where an edge exists between two data points 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; 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.++-}++module Persistence.SimplicialComplex (+  -- * Types+    Simplex+  , SimplicialComplex+  , Graph+  -- * Utilities+  , sc2String+  , getDim+  , encodeWeightedGraph+  , wGraph2sc+  -- * Construction+  , makeNbrhdGraph+  , makeNbrhdGraphPar+  , makeCliqueComplex+  , makeCliqueComplexPar+  , makeRipsComplexFast+  , makeRipsComplexFastPar+  , makeRipsComplexLight+  , makeRipsComplexLightPar+  -- * Homology+  , bettiNumbers+  , bettiNumbersPar+  --, simplicialHomology+  --, simplicialHomologyPar+  ) where++import Persistence.Util+import Persistence.Matrix++import Data.List as L+import Data.Vector as V+import Data.IntSet as S+import Control.Parallel.Strategies+import Data.Algorithm.MaximalCliques++-- * Types++{- |+  This is type synonym exists to make other synonyms more concise.+  A simplex is represented as a pair: the vector of its vertices (their index in the original data set),+  and the vector of the indices of the faces in the next lowest dimension of the simplicial complex.+  1-simplices are the exception:+  they do not store reference to their faces because it would be redundant with their vertices.+-}+type Simplex = (Vector Int, Vector Int)++{- |+  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.+  Index 0 of the vecor corresponds to dimension 1 because there is no need to store individual vertices.+-}+type SimplicialComplex = (Int, Vector (Vector Simplex))++{- |+  This represents the (symmetric) adjacency matrix of some weighted, undirected graph. The type a  is whatever distance is in your data analysis procedure.+  The reason graphs are represented like this is because their main purpose 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))++-- * Simplicial complex utilities++-- | Show all the information in a simplicial complex.+sc2String :: SimplicialComplex -> String+sc2String (v, allSimplices) =+  if V.null allSimplices then (show v) L.++ " vertices"+  else+  let edges     = V.head allSimplices+      simplices = V.tail allSimplices+      showSimplex s =+        '\n':(intercalate "\n" $ V.toList $ V.map show s)+      showAll sc    =+        if V.null sc then '\n':(show v) L.++ " vertices"+        else showSimplex (V.head sc) L.++ ('\n':(showAll (V.tail sc)))+  in (intercalate "\n" $ V.toList $ V.map (show . fst) edges) L.++ ('\n':(showAll simplices))++-- | Get the dimension of the highest dimensional simplex (constant time).+getDim :: SimplicialComplex -> Int+getDim = L.length . snd++{- |+  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 Filtration module which is already an instance of Ord.+-}+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 ())++-- | Given a weighted graph, construct a 1-dimensional simplicial complex in the canonical way. Betti numbers of this simplicial complex can be used to count cycles and connected components.+wGraph2sc :: Graph a -> SimplicialComplex+wGraph2sc graph =+  let numVerts = V.length graph+      getEdges index result =+        if index == numVerts then result+        else+          V.map (\(i, b) ->+            (V.fromList [index, i], V.empty)) $ V.filter snd+              $ mapWithIndex (\i (_, b) -> (i, b)) $ V.take (numVerts - index) $ graph ! index+  in (numVerts, (getEdges 0 V.empty) `cons` V.empty)++{- |+  The first argument is a scale, the second is a metric, and the third is the data.+  This function records the distance between every element of the data and whether or not it is smaller than the given scale.+-}+makeNbrhdGraph :: (Ord a, Eq b) => a -> (b -> b -> a) -> Either (Vector b) [b] -> Graph a+makeNbrhdGraph scale metric dataSet =+  let vector = case dataSet of Left v -> v; Right l -> V.fromList l+  in V.map (\x -> V.map (\y -> let d = metric x y in (d, d <= scale)) vector) vector++-- | Parallel version of the above.+makeNbrhdGraphPar :: (Ord a, Eq b) => a -> (b -> b -> a) -> Either (Vector b) [b] -> Graph a+makeNbrhdGraphPar scale metric dataSet =+  let vector = case dataSet of Left v -> v; Right l -> V.fromList l+  in parMapVec (\x -> V.map (\y -> let d = metric x y in (d, d <= scale)) vector) vector++-- * Simplicial complex construction++{- |+  Makes a simplicial complex where the simplices are the complete subgraphs (cliques) of the given graph.+  Mainly a helper function for makeRipsComplexFast,+  but it might be useful if you happen to have a graph you want to analyze.+  I highly recomend using the parallel version, as this process is very expensive.+-}+makeCliqueComplex :: Graph a -> SimplicialComplex+makeCliqueComplex graph =+  let numVerts = V.length graph++      --make a list with an entry for every dimension of simplices+      organizeCliques 1 _       = []+      organizeCliques i cliques =+        --find the simplices with the given number of vertices+        let helper = V.partition (\simplex -> i == V.length simplex) cliques+        --append them to the next iteration of the function+        in (fst helper):(organizeCliques (i - 1) $ snd helper)++      --pair the organized maximal cliques with the dimension of the largest clique+      makePair simplices =+        case simplices of+          (x:_) ->+            let dim = V.length x+            in (dim, organizeCliques dim $ V.fromList simplices)+          --if there are no maximal cliques this acts as a flag+          --so that the algorithm doesn't try to generate all the other simplices+          []    -> (1, [])++      --find all maximal cliques and sort them from largest to smallest+      --excludes maximal cliques which are single points+      maxCliques :: (Int, Vector (Vector (Vector Int)))+      maxCliques =+        (\(x, y) -> (x, V.fromList y)) $ makePair+          $ sortVecs $ L.map V.fromList $ L.filter (\c -> L.length c > 1)+            $ getMaximalCliques (\i j -> snd $ graph ! i ! j) [0..numVerts - 1]++      --generates faces of simplices and records the boundary indices+      combos :: Int+             -> Int+             -> Vector (Vector (Vector Int))+             -> Vector (Vector (Vector Int, Vector Int))+             -> Vector (Vector (Vector Int, Vector Int))+      combos i max sc result =+        if i == max then+          --don't record boundary indices for edges+          (V.map (\s -> (s, V.empty)) $ V.last sc) `cons` result+        else+          let i1        = i + 1+              --sc is in reverse order, so sc !! i1 is the array of simplices one dimension lower+              current   = sc ! i+              next      = sc ! i1+              len       = V.length next+              allCombos = V.map getCombos current --get all the faces of every simplex+              uCombos   = bigU allCombos --union of faces+              --the index of the faces of each simplex can be found by+              --adding the number of (n-1)-simplices to the index of each face in the union of faces+              indices   = V.map (V.map (\face -> len + (elemIndexUnsafe face uCombos))) allCombos+          in combos i1 max (replaceElem i1 (next V.++ uCombos) sc)+               $ (V.zip current indices) `cons` result++      fstmc2 = fst maxCliques - 2+  in+    --if there are no maximal cliques, the complex is just a bunch of points+    if fstmc2 == (-1) then (numVerts, V.empty)+    else+      let sc = combos 0 fstmc2 (snd maxCliques) V.empty+      in (numVerts, sc)++-- | Parallel version of the above.+makeCliqueComplexPar :: Graph a -> SimplicialComplex+makeCliqueComplexPar graph =+  let numVerts = V.length graph++      --make a list with an entry for every dimension of simplices+      organizeCliques 1 _       = []+      organizeCliques i cliques =+        --find the simplices with the given number of vertices+        let helper = V.partition (\simplex -> i == V.length simplex) cliques+        --append them to the next iteration of the function+        in (fst helper):(organizeCliques (i - 1) $ snd helper)++      --pair the organized maximal cliques with the dimension of the largest clique+      makePair simplices =+        case simplices of+          (x:_) ->+            let dim = V.length x+            in (dim, organizeCliques dim $ V.fromList simplices)+          --if there are no maximal cliques this acts as a flag+          --so that the algorithm doesn't try to generate all the other simplices+          []    -> (1, [])++      --find all maximal cliques and sort them from largest to smallest+      --excludes maximal cliques which are single points+      maxCliques :: (Int, Vector (Vector (Vector Int)))+      maxCliques =+        (\(x, y) -> (x, V.fromList y)) $ makePair+          $ sortVecs $ L.map V.fromList $ L.filter (\c -> L.length c > 1)+            $ getMaximalCliques (\i j -> snd $ graph ! i ! j) [0..numVerts - 1]++      --generates faces of simplices and records the boundary indices+      combos :: Int+             -> Int+             -> Vector (Vector (Vector Int))+             -> Vector (Vector (Vector Int, Vector Int))+             -> Vector (Vector (Vector Int, Vector Int))+      combos i max sc result =+        if i == max then+          --don't record boundary indices for edges+          (V.map (\s -> (s, V.empty)) $ V.last sc) `cons` result+        else+          let i1        = i + 1+              --sc is in reverse order, so sc !! i1 is the array of simplices one dimension lower+              current   = sc ! i+              next      = sc ! i1+              len       = V.length next+              allCombos = V.map getCombos current --get all the faces of every simplex+              uCombos   = bigU allCombos --union of faces+              --the index of the faces of each simplex can be found by adding+              --the number of (n-1)-simplices to the index of each face in the union of faces+              indices   =+                parMapVec (V.map (\face -> len + (elemIndexUnsafe face uCombos))) allCombos+          in combos i1 max (replaceElem i1 (next V.++ uCombos) sc)+               $ (V.zip current indices) `cons` result++      fstmc2 = fst maxCliques - 2+  in+    --if there are no maximal cliques, the complex is just a bunch of points+    if fstmc2 == (-1) then (numVerts, V.empty)+    else+      let sc = combos 0 fstmc2 (snd maxCliques) V.empty+      in (numVerts, sc)++{- |+  Constructs the Vietoris-Rips complex given a scale, metric, and data set.+  Also uses O(n^2) memory (where n is the number of data points)+  for a graph storing all the distances between data points.+-}+makeRipsComplexFast :: (Ord a, Eq b)+                    => a+                    -> (b -> b -> a)+                    -> Either (Vector b) [b]+                    -> (SimplicialComplex, Graph a)+makeRipsComplexFast scale metric dataSet =+  let graph = makeNbrhdGraph scale metric dataSet+      sc    = makeCliqueComplex graph+  in (sc, graph)++-- | Parallel version of the above.+makeRipsComplexFastPar :: (Ord a, Eq b)+                       => a+                       -> (b -> b -> a)+                       -> Either (Vector b) [b]+                       -> (SimplicialComplex, Graph a)+makeRipsComplexFastPar scale metric dataSet =+  let graph = makeNbrhdGraphPar scale metric dataSet+      sc    = makeCliqueComplexPar graph+  in (sc, graph)++-- | Constructs the Vietoris-Rips complex given a scale, metric, and data set.+makeRipsComplexLight :: (Ord a, Eq b)+                     => a+                     -> (b -> b -> a)+                     -> Either (Vector b) [b]+                     -> SimplicialComplex+makeRipsComplexLight scale metric dataSet =+  let vector   = case dataSet of Left v -> v; Right l -> V.fromList l+      numVerts = L.length vector++      --make a list with an entry for every dimension of simplices+      organizeCliques 1 _       = []+      organizeCliques i cliques =+        --find the simplices with the given number of vertices+        let helper = V.partition (\simplex -> i == V.length simplex) cliques+        --append them to the next iteration of the function+        in (fst helper):(organizeCliques (i - 1) $ snd helper)++      --pair the organized maximal cliques with the dimension of the largest clique+      makePair simplices =+        case simplices of+          (x:_) ->+            let dim = V.length x+            in (dim, organizeCliques dim $ V.fromList simplices)+          --if there are no maximal cliques this acts as a flag+          --so that the algorithm doesn't try to generate all the other simplices+          []    -> (1, [])++      --find all maximal cliques and sort them from largest to smallest+      --excludes maximal cliques which are single points+      maxCliques :: (Int, Vector (Vector (Vector Int)))+      maxCliques =+        (\(x, y) -> (x, V.fromList y)) $ makePair $ sortVecs+          $ L.map V.fromList $ L.filter (\c -> L.length c > 1)+            $ getMaximalCliques (\i j -> metric (vector ! i) (vector ! j) <= scale) [0..numVerts-1]++      --generates faces of simplices and records the boundary indices+      combos :: Int+             -> Int+             -> Vector (Vector (Vector Int))+             -> Vector (Vector (Vector Int, Vector Int))+             -> Vector (Vector (Vector Int, Vector Int))+      combos i max sc result =+        if i == max then+          --don't record boundary indices for edges+          (V.map (\s -> (s, V.empty)) $ V.last sc) `cons` result+        else+          let i1        = i + 1+              --sc is in reverse order, so sc !! i1 is the array of simplices one dimension lower+              current   = sc ! i+              next      = sc ! i1+              len       = V.length next+              allCombos = V.map getCombos current --get all the faces of every simplex+              uCombos   = bigU allCombos --union of faces+              --the index of the faces of each simplex can be found by+              --adding the number of (n-1)-simplices to the index of each face in the union of faces+              indices   = V.map (V.map (\face ->+                            len + (elemIndexUnsafe face uCombos))) allCombos+          in combos i1 max (replaceElem i1 (next V.++ uCombos) sc)+               $ (V.zip current indices) `cons` result++      fstmc2 = fst maxCliques - 2+  in+    --if there are no maximal cliques, the complex is just a bunch of points+    if fstmc2 == (-1) then (numVerts, V.empty)+    else+      let sc = combos 0 fstmc2 (snd maxCliques) V.empty+      in (numVerts, sc)++-- | Parallel version of the above.+makeRipsComplexLightPar :: (Ord a, Eq b)+                        => a+                        -> (b -> b -> a)+                        -> Either (Vector b) [b]+                        -> SimplicialComplex+makeRipsComplexLightPar scale metric dataSet =+  let numVerts = L.length dataSet+      vector   = case dataSet of Left v -> v; Right l -> V.fromList l++      --make a list with an entry for every dimension of simplices+      organizeCliques 1 _       = []+      organizeCliques i cliques =+        --find the simplices with the given number of vertices+        let helper = V.partition (\simplex -> i == V.length simplex) cliques+        --append them to the next iteration of the function+        in (fst helper):(organizeCliques (i - 1) $ snd helper)++      --pair the organized maximal cliques with the dimension of the largest clique+      makePair simplices =+        case simplices of+          (x:_) ->+            let dim = V.length x+            in (dim, organizeCliques dim $ V.fromList simplices)+          --if there are no maximal cliques this acts as a flag+          --so that the algorithm doesn't try to generate all the other simplices+          []    -> (1, [])++      --find all maximal cliques and sort them from largest to smallest+      --excludes maximal cliques which are single points+      maxCliques :: (Int, Vector (Vector (Vector Int)))+      maxCliques =+        (\(x, y) -> (x, V.fromList y)) $ makePair+          $ sortVecs $ L.map V.fromList $ L.filter (\c -> L.length c > 1)+            $ getMaximalCliques (\i j -> metric (vector ! i) (vector ! j) <= scale) [0..numVerts-1]++      --generates faces of simplices and records the boundary indices+      combos :: Int+             -> Int+             -> Vector (Vector (Vector Int))+             -> Vector (Vector (Vector Int, Vector Int))+             -> Vector (Vector (Vector Int, Vector Int))+      combos i max sc result =+        if i == max then+          --don't record boundary indices for edges+          (V.map (\s -> (s, V.empty)) $ V.last sc) `cons` result+        else+          let i1        = i + 1+              --sc is in reverse order, so sc !! i1 is the array of simplices one dimension lower+              current   = sc ! i+              next      = sc ! i1+              len       = V.length next+              allCombos = V.map getCombos current --get all the faces of every simplex+              uCombos   = bigU allCombos --union of faces+              --the index of the faces of each simplex can be found by+              --adding the number of (n-1)-simplices to the index of each face in the union of faces+              indices   = parMapVec (V.map (\face ->+                            len + (elemIndexUnsafe face uCombos))) allCombos+          in combos i1 max (replaceElem i1 (next V.++ uCombos) sc)+               $ (V.zip current indices) `cons` result++      fstmc2 = fst maxCliques - 2+  in+    --if there are no maximal cliques, the complex is just a bunch of points+    if fstmc2 == (-1) then (numVerts, V.empty)+    else+      let sc = combos 0 fstmc2 (snd maxCliques) V.empty+      in (numVerts, sc)++-- * Homology++--gets the first boundary operator (because edges don't need to point to their subsimplices)+makeEdgeBoundariesBool :: SimplicialComplex -> BMatrix+makeEdgeBoundariesBool sc =+  let mat =+        V.map (\edge ->+          V.map (\vert -> vert == V.head edge || vert == V.last edge) $ 0 `range` (fst sc - 1))+            $ V.map fst $ V.head $ snd sc+  in+    case transposeMat mat of+      Just m  -> m+      Nothing -> error "error in makeEdgeBoundariesBool"++--gets the boundary coefficients for a simplex of dimension 2 or greater+--first argument is dimension of the simplex+--second argument is the simplicial complex+--third argument is the simplex paired with the indices of its faces+makeSimplexBoundaryBool :: Int -> SimplicialComplex -> (Vector Int, Vector Int) -> Vector Bool+makeSimplexBoundaryBool dim simplices (simplex, indices) =+  mapWithIndex (\i s -> V.elem i indices) (V.map fst $ (snd simplices) ! (dim - 2))++--makes boundary operator for all simplices of dimension 2 or greater+--first argument is the dimension of the boundary operator, second is the simplicial complex+makeBoundaryOperatorBool :: Int -> SimplicialComplex -> BMatrix+makeBoundaryOperatorBool dim sc =+  let mat = V.map (makeSimplexBoundaryBool dim sc) $ (snd sc) ! (dim - 1)+  in+    case transposeMat mat of+      Just m  -> m+      Nothing -> error "error in makeBoundaryOperatorBool"++--makes all the boundary operators+makeBoundaryOperatorsBool :: SimplicialComplex -> Vector BMatrix+makeBoundaryOperatorsBool sc =+  let dim = getDim sc+      calc i+        | i > dim   = V.empty+        | i == 1    = (makeEdgeBoundariesBool sc) `cons` (calc 2)+        | otherwise = (makeBoundaryOperatorBool i sc) `cons` (calc (i + 1))+  in calc 1++{- |+  The nth index of the output corresponds to the nth Betti number.++  The zeroth Betti number is the number of connected components, the 1st is the number of tunnels/  punctures, the 2nd is the number of hollow volumes, and so on for higher-dimensional holes.+-}+bettiNumbers :: SimplicialComplex -> [Int]+bettiNumbers sc =+  let dim      = (getDim sc) + 1+      boundOps = makeBoundaryOperatorsBool sc+      --dimension of image paired with dimension of kernel+      ranks    = (0, V.length $ V.head boundOps)+                   `cons` (V.map (\op -> let rank = rankBool op+                     in (rank, (V.length $ V.head op) - rank)) boundOps)+      calc 1   = [(snd $ ranks ! 0) - (fst $ ranks ! 1)]+      calc i   =+        let i1 = i - 1+        in+          if i == dim then (snd $ V.last ranks):(calc i1)+          else ((snd $ ranks ! i1) - (fst $ ranks ! i)):(calc i1)+  in+    if L.null $ snd sc then [fst sc]+    else L.reverse $ calc dim++-- | Calculates all of the Betti numbers in parallel.+bettiNumbersPar :: SimplicialComplex -> [Int]+bettiNumbersPar sc =+  let dim      = (getDim sc) + 1+      boundOps = makeBoundaryOperatorsBool sc+      --dimension of image paired with dimension of kernel+      ranks    = (0, V.length $ V.head boundOps)+                   `cons` (parMapVec (\op -> let rank = rankBool op+                     in (rank, (V.length $ V.head op) - rank)) boundOps)+      calc 1   = [(snd $ ranks ! 0) - (fst $ ranks ! 1)]+      calc i   =+        let i1 = i - 1+        in+          if i == dim then evalPar (snd $ V.last ranks) (calc i1) --see Util for evalPar+          else evalPar ((snd $ ranks ! i1) - (fst $ ranks ! i)) (calc i1)+  in+    if L.null $ snd sc then [fst sc]+    else L.reverse $ calc dim++--gets the first boundary operator+--because edges don't need to point to their faces+makeEdgeBoundariesInt :: SimplicialComplex -> IMatrix+makeEdgeBoundariesInt sc =+  let mat =+        V.map (\e -> let edge = fst e in+            replaceElem (V.head edge) (-1) $+              replaceElem (V.last edge) 1 $+                V.replicate (fst sc) 0) $+                  V.head $ snd sc+  in+    case transposeMat mat of+      Just m  -> m+      Nothing -> error "error in makeEdgeBoundariesInt"++--gets the boundary coefficients for a simplex of dimension 2 or greater+--first argument is the simplices of one dimension lower+--second argument is the simplex paired with the indices of its faces+makeSimplexBoundaryInt :: Vector (Vector Int, Vector Int) -> (Vector Int, Vector Int) -> Vector Int+makeSimplexBoundaryInt simplices (_, indices) =+  let calc1 ixs result =+        if V.null ixs then result+        else calc2 (V.tail ixs) $ replaceElem (V.head ixs) (-1) result+      calc2 ixs result =+        if V.null ixs then result+        else calc1 (V.tail ixs) $ replaceElem (V.head ixs) 1 result+  in calc1 indices $ V.replicate (V.length simplices) 0++--makes boundary operator for all simplices of dimension 2 or greater+--first argument is the dimension of the boundary operator, second is the simplicial complex+makeBoundaryOperatorInt :: Int -> SimplicialComplex -> IMatrix+makeBoundaryOperatorInt dim sc =+  let mat = V.map (makeSimplexBoundaryInt ((snd sc) ! (dim - 2))) $ (snd sc) ! (dim - 1)+  in+    case transposeMat mat of+      Just m  -> m+      Nothing -> error $ show mat++--makes all the boundary operators+makeBoundaryOperatorsInt :: SimplicialComplex -> Vector IMatrix+makeBoundaryOperatorsInt sc =+  let dim = getDim sc+      calc 1 = (makeEdgeBoundariesInt sc) `cons` (calc 2)+      calc i =+        if i > dim then V.empty+        else (makeBoundaryOperatorInt i sc) `cons` (calc $ i + 1)+  in calc 1++{- |+  Calculates all homology groups of the complex. An int list represents a single homology group.++  Specifically, an int list represents a direct product of cyclic groups - a 0 entry corresponds to the infinite cylic group (also known as the integers), and an entry of n corresponds to a cyclic group of order n (also known as the integers modulo n). The nth index of the output represents the nth homology group.++  If you don't have a rigorous understanding of what homology is, I recomend using BettiNumbers on large data-sets instead, as the output of those functions is more easily interprettable.+-}+simplicialHomology :: SimplicialComplex -> [[Int]]+simplicialHomology sc =+  let dim      = getDim sc+      boundOps = makeBoundaryOperatorsInt sc+      calc 0   = [getUnsignedDiagonal $ normalFormInt (boundOps ! 0)]+      calc i   =+        if i == dim then+          let op = V.last boundOps+          in (L.replicate ((V.length $ V.head op) - (rankInt op)) 0):(calc $ i - 1)+        else+          let i1 = i - 1+          in (getUnsignedDiagonal $ normalFormInt $ imgInKerInt (boundOps ! i1) (boundOps ! i)):(calc i1)+  in+    if L.null $ snd sc then [L.replicate (fst sc) 0]+    else L.reverse $ L.map (L.filter (/=1)) $ calc dim++-- | Same as simplicialHomology except it computes each of the groups in parallel and uses parallel matrix computations.+simplicialHomologyPar :: SimplicialComplex -> [[Int]]+simplicialHomologyPar sc =+  let dim      = getDim sc+      boundOps = makeBoundaryOperatorsInt sc+      calc 0   = [getUnsignedDiagonal $ normalFormInt (boundOps ! 0)]+      calc i   =+        if i == dim then+          let op = V.last boundOps+          in evalPar (L.replicate ((V.length $ V.head op) - (rankInt op)) 0) $ calc $ i - 1+        else+          let i1 = i - 1+          in evalPar (getUnsignedDiagonal $ normalFormIntPar $ --see Util for evalPar+            imgInKerIntPar (boundOps ! i1) (boundOps ! i)) $ calc i1+  in+    if L.null $ snd sc then [L.replicate (fst sc) 0]+    else L.reverse $ L.map (L.filter (/=1)) $ calc dim

@@ -0,0 +1,329 @@+{- |+Module     : Persistence.Util+Copyright  : (c) Eben Kadile, 2018+License    : BSD 3 Clause+Maintainer : eben.cowley42@gmail.com+Stability  : experimental++This module contains miscellaneous utility functions used throughout the Persistence library.++-}++module Persistence.Util where++import Data.List as L+import Data.Vector as V+import Control.Parallel.Strategies++-- | Simple instance of Num where True is 1 and False is 0, all operations work like arithmetic modulo 2.+instance Num Bool where+  p + q  = p `xor` q+  p * q  = p && q+  p - q  = p `xor` (not q)+  negate = not+  abs    = id+  fromInteger 0 = False+  fromInteger _ = True+  signum bool   = if bool then 1 else 0++-- | Exclusive or.+xor :: Bool -> Bool -> Bool+xor False False = False+xor True False  = True+xor False True  = True+xor True True   = False++-- | First element of a triple.+one (a, _, _) = a++-- | Second element of a triple.+two (_, b, _) = b++-- | Third element of a triple.+thr (_, _, c) = c++-- | Last two elements of a triple.+not1 (_, b, c) = (b, c)++-- | First and last elements of a triple.+not2 (a, _, c) = (a, c)++-- | First two elements of a triple.+not3 (a, b, _) = (a, b)++-- | Concatenate a vector of vectors.+flatten :: Vector (Vector a) -> Vector a+flatten = V.foldl1 (V.++)++-- | Multiply a vector by a scalar.+mul :: Num a => a -> Vector a -> Vector a+mul s = V.map (*s)++{- |+  Add two vectors together component-wise.+  WARNING: If one vector is longer than the other,+  the longer vector will simply be cut off.+-}+add :: Num a => Vector a -> Vector a -> Vector a+add = V.zipWith (+)++{- |+  Subtract the second vector from the first vector component-wise.+  WARNING: If one vector is longer than the other,+  the longer vector will simply be cut off.+-}+subtr :: Num a => Vector a -> Vector a -> Vector a+subtr = V.zipWith (\x y -> x - y)++-- | Dot product. WARNING: Vectors must be the same length+dotProduct :: Num a => Vector a -> Vector a -> a+dotProduct vec1 vec2+  | a && b = fromIntegral 0+  | a      = error "Persistence.Util.dotProduct (first argument too short). This is a bug. Please email the Persistence maintainers."+  | b      = error "Persistence.Util.dotProduct (second argument too short). This is a bug. Please email the Persistence maintainers."+  | otherwise   = (V.head vec1)*(V.head vec2) + (dotProduct (V.tail vec1) (V.tail vec2))+    where a = V.null vec1; b = V.null vec2++-- | Extended Euclidean algorithm. Finds the gcd of the two inputs plus the coefficients that multiply each input and sum to give the gcd.+extEucAlg :: Integral a => a -> a -> (a, a, a)+extEucAlg a b =+  let eeaHelper r s t =+        case snd r of+          0 -> (fst r, fst s, fst t)+          _ ->+            let r1    = fst r+                r2    = snd r+                s2    = snd s+                t2    = snd t+                q     = r1 `div` r2+                nextr = r1 - q*r2+                nexts = fst s - q*s2+                nextt = fst t - q*t2+            in eeaHelper (r2, nextr) (s2, nexts) (t2, nextt)+  in (\(x, y, z) -> if x < 0 then (-x, -y, -z) else (x, y, z)) $ eeaHelper (a, b) (0, 1) (1, 0)++-- | Returns whether or not the first number divides the second number. Returns false if the first number is 0.+divides :: Int -> Int -> Bool+0 `divides` b = False+a `divides` b+  | b < 0     = False+  | b == 0    = True+  | otherwise = a `divides` (b - (abs a))++-- | Switches the elements of the vector at the given indices.+switchElems ::Int -> Int -> Vector a -> Vector a+switchElems i j vector+  | j == i    = vector+  | j < i     =+    let first  = V.take j vector+        second = V.drop (j + 1) (V.take i vector)+        third  = V.drop (i + 1) vector+    in first V.++ (cons (vector ! i) second) V.++ (cons (vector ! j) third)+  | otherwise =+    let first  = V.take i vector+        second = V.drop (i + 1) (V.take j vector)+        third  = V.drop (j + 1) vector+    in first V.++ (cons (vector ! j) second) V.++ (cons (vector ! i) third)++-- | Remove the element at the given index.+rmIndex :: Int -> Vector a -> Vector a+rmIndex i v = (V.take i v) V.++ (V.drop (i + 1) v)++-- | Generate a range of integers in vector form.+range :: Int -> Int -> Vector Int+range x y+  | x == y = x `cons` empty+  | x < y  = x `cons` (range (x + 1) y)+  | x > y  = (range x (y + 1)) `snoc` y++-- | Return all vectors missing exactly one element from the original vector.+getCombos :: Vector a -> Vector (Vector a)+getCombos vector = V.map (\i -> rmIndex i vector) $ 0 `range` (V.length vector - 1)++-- | Drop all elements from the right until the predicate is unsatisfied.+dropRightWhile :: (a -> Bool) -> Vector a -> Vector a+dropRightWhile p v = if p (V.last v) then dropRightWhile p (V.init v) else v++-- | Filter a vector with a predicate that takes into account the index of the element.+filterWithIndex :: (Int -> a -> Bool) -> Vector a -> Vector a+filterWithIndex p vector =+  let maxIndex = V.length vector - 1+      calc i+        | i == maxIndex    = V.empty+        | p i (vector ! i) = (vector ! i) `cons` calc (i + 1)+        | otherwise        = calc (i + 1)+  in calc 0++-- | Map a function that takes into account the index of each element.+mapWithIndex :: (Int -> a -> b) -> Vector a -> Vector b+mapWithIndex f vector =+  let helper i vec =+        if V.null vec then empty+        else cons (f i $ V.head vec) $ helper (i + 1) (V.tail vec)+  in helper 0 vector++-- | Parallel map a function over a vector.+parMapVec :: (a -> b) -> Vector a -> Vector b+parMapVec f = runEval . (evalTraversable rpar) . (V.map f)++-- | Map a function that takes into account the index of each element in parallel.+parMapWithIndex :: (Int -> a -> b) -> Vector a -> Vector b+parMapWithIndex f = runEval . (evalTraversable rpar) . (mapWithIndex f)++-- | Return the element satisfying the predicate and its index if it exists.+elemAndIndex :: (a -> Bool) -> Vector a -> Maybe (a, Int)+elemAndIndex p vector =+  let helper i vec+        | V.null vec     = Nothing+        | p $ V.head vec = Just (V.head vec, i)+        | otherwise      = helper (i + 1) $ V.tail vec+  in helper 0 vector++-- | Return the elements satisfying the predicate and their indices.+elemAndIndices :: (a -> Bool) -> Vector a -> [(a, Int)]+elemAndIndices p vector =+  let helper i vec+        | V.null vec     = []+        | p $ V.head vec = (V.head vec, i) : (helper (i + 1) $ V.tail vec)+        | otherwise      = helper (i + 1) $ V.tail vec+  in helper 0 vector++-- | 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 len = V.length vector1++      calc i result =+        let a = vector1 ! i+        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++-- | Orders a list of vectors from greatest to least length.+sortVecs :: [Vector a] -> [Vector a]+sortVecs []     = []+sortVecs (v:vs) =+  let len  = V.length v+      less = sortVecs $ L.filter (\u -> V.length u < len) vs+      more = sortVecs $ L.filter (\u -> V.length u >= len) vs+  in more L.++ [v] L.++ less++-- | Replace the element at the given index with the given element.+replaceElem :: Int -> a -> Vector a -> Vector a+replaceElem i e v = (V.take i v) V.++ (e `cons` (V.drop (i + 1) v))++-- | Replace the element at the given index with the given element.+replaceElemList :: Int -> a -> [a] -> [a]+replaceElemList i e l = (L.take i l) L.++ (e:(L.drop (i + 1) l))++-- | Quicksort treating the given predicate as the <= operator. Sorts in decreasing order.+quickSort :: (a -> a -> Bool) -> Vector a -> Vector a+quickSort rel vector = --rel is the <= operator+  if V.null vector then empty+  else+    let x       = V.head vector+        xs      = V.tail vector+        greater = V.filter (rel x) xs+        lesser  = V.filter (not . (rel x)) xs+    in (quickSort rel greater) V.++ (x `cons` (quickSort rel lesser))++-- | If the predicate is the <= operator and the vector is sorted in decreasing order, this inserts the value in the correct position.+orderedInsert :: (a -> a -> Bool) -> a -> Vector a -> Vector a+orderedInsert rel x vector = --rel is the <= operator+  case V.findIndex (\y -> y `rel` x) vector of+    Just i  ->+      case V.findIndex (\y -> x `rel` y) $ V.drop i vector of+        Just j  -> (V.take (i + j) vector) V.++ (x `cons` (V.drop (i + j) vector))+        Nothing -> (V.take i vector) V.++ (x `cons` (V.drop i vector))+    Nothing -> vector `snoc` x++-- | Takes the union of all of the vectors.+bigU :: Eq a => Vector (Vector a) -> Vector a+bigU =+  let exists x v+        | V.null v      = False+        | V.head v == x = True+        | otherwise     = exists x (V.tail v)+      union v1 v2 =+        if V.null v1 then v2+        else+          let x = V.head v1+          in+            if exists x v2 then union (V.tail v1) v2+            else union (V.tail v1) (x `cons` v2)+  in V.foldl1 union++-- | The element being searched for, the vector being searched, and the lower and upper bounds on the indices.+binarySearch :: Ord a => a -> Vector a -> Int -> Int -> Maybe Int+binarySearch value xs low high+  | high < low        = Nothing+  | xs ! mid > value  = binarySearch value xs low (mid - 1)+  | xs ! mid < value  = binarySearch value xs (mid + 1) high+  | otherwise         = Just mid+  where mid = low + ((high - low) `div` 2)++-- | Intersection of SORTED vectors.+(|^|) :: Ord a => Vector a -> Vector a -> Vector a+vector1 |^| vector2 =+  let len          = V.length vector2 - 1+      calc acc v =+        if V.null v then acc+        else+          let x = V.head v; xs = V.tail v+          in case binarySearch x vector2 0 len of+            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++-- | If the relation were the "greater than" operator, this would find the minimum element of the vector.+foldRelation :: (a -> a -> Bool) -> Vector a -> a+foldRelation rel vec =+  let calc w v+        | V.null v  = w+        | rel w x   = calc x xs+        | otherwise = calc w xs+        where x = V.head v; xs = V.tail v+  in calc (V.head vec) (V.tail vec)++-- | Unsafe index finding.+elemIndexUnsafe :: Eq a => a -> Vector a -> Int+elemIndexUnsafe elem vector =+  let find i v+        | V.null v         = error "Persistence.Util.elemIndexUnsafe. This is a bug. Please email the Persistence maintainers."+        | V.head v == elem = i+        | otherwise        = find (i + 1) $ V.tail v+  in find 0 vector++{- |+  Spark the first argument for parallel evaluation and force evaluation of the second argument,+  then return the first argument concatenated to the second. This is useful especially if the second+  argument is a recursive call that calls evalPar again, so that every elemment of the list will be+  sparked for parallelism.+-}+evalPar :: a -> [a] -> [a]+evalPar c r = runEval $ rpar c >> rseq r >> return (c:r)++-- | Union minus intersection of ordered vectors.+uin :: Ord a => Vector a -> Vector a -> Vector a+u `uin` v =+  let len = V.length v+      findAndInsert i elem vec+        | i == len  = vec `snoc` elem+        | elem == x = (V.take i vec) V.++ (V.drop i1 vec)+        | elem >  x = (V.take i vec) V.++ (elem `cons` (V.drop i vec))+        | otherwise = findAndInsert i1 elem vec+        where x = vec ! i; i1 = i + 1+  in+    if V.null u then v+    else (V.tail u) `uin` (findAndInsert 0 (V.head u) v)

@@ -0,0 +1,598 @@+{-# LANGUAGE FlexibleInstances #-}++import Persistence.Util+import Persistence.Matrix+import Persistence.SimplicialComplex+import Persistence.HasseDiagram+import Persistence.Filtration++import Data.Vector as V+import Data.List as L++class ToString a where+  toString :: a -> String++instance ToString Double where+  toString = show++instance ToString Int where+  toString = show++instance ToString Bool where+  toString b = if b then "1" else "0"++instance ToString IMatrix where+  toString 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++instance ToString [Int] where+  toString = show++instance ToString [[Int]] where+  toString = show++instance Show a => ToString (Vector (Vector (BarCode a))) where+  toString = unlines . V.toList . (V.map show)++instance ToString Landscape where+  toString = show++instance ToString a => ToString (Maybe a) where+  toString (Just a) = "Just " L.++ (toString a)+  toString Nothing  = "Nothing"++data TestResult = Success | Failure String deriving Eq++checkPass :: (Eq a, ToString a) => String -> a -> Maybe a -> TestResult+checkPass msg actual expect =+  case expect of+    Just e  ->+      if actual == e then Success+      else+        Failure $ msg L.++ ":\nActual:\n"+          L.++ (toString actual) L.++ "\nExpected:\n" L.++ (toString e)+    Nothing ->+      Failure $ msg L.++ ":\nActual:\n" L.++ (toString actual) L.++ "\nExpected: ?"+++formatTestResults :: [TestResult] -> String+formatTestResults results =+  let numTests  = L.length results+      numPassed = L.length $ L.filter (\r -> r == Success) results+      summary   = (show numPassed) L.++ "/" L.++ (show numTests) L.++ " Tests Passed.\n"+      accum a r =+        case r of+          Success   -> a+          Failure m -> a L.++ [m L.++ "\n"]+      details   = L.foldl accum [] results+  in '\n':(unlines $ summary:details)++matrix1 :: IMatrix+matrix1 = V.fromList $ L.map V.fromList+  [ [ 2, 3, 5]+  , [-4, 2, 3]+  ]++kernel1 :: IMatrix+kernel1 = V.fromList $ L.map V.fromList [[-1, -26, 16]]++matrix2 :: IMatrix+matrix2 = V.fromList $ L.map V.fromList+  [ [1, -1,  3]+  , [5,  6, -4]+  , [7,  4,  2]+  ]++kernel2 :: IMatrix+kernel2 = V.fromList $ L.map V.fromList [[14, -19, -11]]++matrix3 :: IMatrix+matrix3 = V.fromList $ L.map V.fromList $+  [ [1, 0, -3, 0,  2, -8]+  , [0, 1,  5, 0, -1,  4]+  , [0, 0,  0, 1,  7, -9]+  , [0, 0,  0, 0,  0,  0]+  ]++kernel3 :: IMatrix+kernel3 = V.fromList $ L.map V.fromList $+  [ [ 3, -5, 1,  0, 0, 0]+  , [-2,  1, 0, -7, 1, 0]+  , [ 8, -4, 0,  9, 0, 1]+  ]++matrix4 :: IMatrix+matrix4 = V.fromList $ L.map V.fromList+  [ [ 2,  4,   4]+  , [-6,  6,  12]+  , [10, -4, -16]+  ]++snf4 :: IMatrix+snf4 = V.fromList $ L.map V.fromList+  [ [2, 0,  0]+  , [0, 6,  0]+  , [0, 0, 12]+  ]++matrix5 :: IMatrix+matrix5 = V.fromList $ L.map V.fromList $+  [ [  9,  -36,   30]+  , [-36,  192, -180]+  , [ 30, -180,  180]+  ]++snf5 :: IMatrix+snf5 = V.fromList $ L.map V.fromList $+  [ [3, 0, 0]+  , [0, 12, 0]+  , [0, 0, 60]+  ]++matrix6 :: IMatrix+matrix6 = V.fromList $ L.map V.fromList $+  [ [  9, 0, -36,  30]+  , [  0, 0,   0,   0]+  , [-36, 0,  19, -18]+  , [ 30, 0, -18,  18]+  ]++snf6 :: IMatrix+snf6 = V.fromList $ L.map V.fromList $+  [ [ 1, 0,   0, 0]+  , [ 0, 3,   0, 0]+  , [ 0, 0, 462, 0]+  , [ 0, 0,   0, 0]+  ]++pointCloud1 = Right+  [ ( -7, -19)+  , ( -6, -16)+  , (  2, -16)+  , ( 11, -16)+  , ( 11, -11)+  , ( -6,  -9)+  , (  2,  -9)+  , ( 11,  -6)+  , (-11,  -5)+  , (  6,  -3)+  , (  0,   0)+  , ( -2,   3)+  , ( 11,   3)+  , (-14,   8)+  , (  1,   8)+  , ( -8,  14)+  , ( 10,  19)+  , ( 17,  20)+  ]++metric2 :: Floating a => (a, a) -> (a, a) -> a+metric2 (a, b) (c, d) =+  let x = a - c; y = b - d+  in sqrt $ x*x + y*y++scales1 = Right [10.0, 8.0, 6.0]+testFiltration1 = makeRipsFiltrationFast scales1 metric2 pointCloud1+nsTestFiltration1 = simple2Filtr testFiltration1++pointCloud2 = Right+  [ (  1, -22)+  , (  4, -21)+  , ( -3, -20)+  , (  4, -19)+  , ( -7, -18)+  , (  5, -17)+  , ( -7, -12)+  , ( -1, -12)+  , (  3, -12)+  , (  6, -12)+  , (  8, -12)+  , ( -9, -10)+  , ( -4,  -9)+  , (-12,  -8)+  , (  7,  -8)+  , ( -8,  -7)+  , ( -5,  -7)+  , (-11,  -6)+  , (  9,  -5)+  , ( -6,  -4)+  , (-14,  -3)+  , (-12,  -3)+  , (  7,  -3)+  , (  5,  -1)+  , ( -7,   0)+  , (  0,   0)+  , (  8,   0)+  , (-12,   1)+  , (  2,   1)+  , (  0,   2)+  , (-10,   3)+  , ( -7,   3)+  , (  6,   3)+  , ( -2,   5)+  , ( 10,   5)+  , ( -8,   9)+  , (  8,   9)+  , ( -9,  11)+  , (  6,  11)+  , ( -6,  12)+  , ( -4,  13)+  , (  2,  13)+  , (  8,  13)+  , ( -2,  14)+  , (  5,  14)+  , ( -4,  16)+  , (  1,  16)+  ]++--8 connected components at index 0, 2 connected components at index 1, 1 connected component at index 2+--1 loop lasting from 0 to 1, 1 loop lasting from 1 to 2, 1 loop starting at 1, 1 loop starting at 2+scales2 = Right [10.0, 8.0, 6.0]+testFiltration2 = makeRipsFiltrationFast scales2 metric2 pointCloud2+nsTestFiltration2 = simple2Filtr testFiltration2++octahedralCloud = Right+  [ ( 1, 0, 0)+  , (-1, 0, 0)+  , ( 0, 1, 0)+  , ( 0,-1, 0)+  , ( 0, 0, 1)+  , ( 0, 0,-1)+  ]++metric3 :: Floating a => (a, a, a) -> (a, a, a) -> a+metric3 (x0, y0, z0) (x1, y1, z1) =+  let dx = x1 - x0; dy = y1 - y0; dz = z1 - z0 in+  sqrt (dx*dx + dy*dy + dz*dz)++octaScales = Right [2.1, 1.6, 1.1, 0.5]++--5 connected components from 0 to 2, 1 void from 2 to 3+octahedron = makeRipsFiltrationFast octaScales metric3 octahedralCloud+nsoctahedron = simple2Filtr octahedron++sqrCloud = Right+  [ (0, 0)+  , (1, 0)+  , (2, 0)+  , (3, 0)+  , (3, 1)+  , (3, 2)+  , (3, 3)+  , (2, 3)+  , (1, 3)+  , (0, 3)+  , (0, 2)+  , (0, 1)+  ]++sqrScales = Right [2.0, 1.5, 1.1, 0.5]++--11 connected components from 0 to 1, 1 loop starting at 1+square = makeRipsFiltrationFast sqrScales metric2 sqrCloud+nssquare = simple2Filtr square++dGraph1 =+  [ ( 0,  1)+  , ( 0,  2)+  , ( 1,  2)+  , ( 1,  3)+  , ( 1,  4)+  , ( 1,  6)+  , ( 1,  7)+  , ( 2,  4)+  , ( 2,  6)+  , ( 2,  7)+  , ( 3,  2)+  , ( 3,  4)+  , ( 3,  6)+  , ( 3,  7)+  , ( 4,  5)+  , ( 4,  6)+  , ( 4,  7)+  , ( 5,  9)+  , ( 5, 11)+  , ( 6,  7)+  , ( 6, 10)+  , ( 6, 12)+  , ( 7,  8)+  , ( 7, 12)+  , ( 8,  4)+  , ( 8,  5)+  , ( 8,  9)+  , (11,  8)+  , (11,  9)+  , (11, 14)+  , (12, 10)+  , (12, 13)+  , (13, 15)+  , (13, 17)+  , (14, 13)+  , (14, 15)+  , (14, 16)+  , (14, 17)+  , (16, 13)+  , (16, 15)+  , (16, 17)+  , (17, 15)+  ]++directedGraph = encodeDirectedGraph 18 dGraph1+dCliqueComplex = hDiagram2sc $ directedFlagComplex directedGraph++--Mobius strip created in blender+mobius = Right+  [ ( 1.250000,  0.000000,  0.000000)+  , ( 0.750000,  0.000000,  0.000000)+  , (-1.250000,  0.000000,  0.000000)+  , (-1.199179,  0.064705,  0.321319)+  , (-1.053526,  0.125000,  0.608253)+  , (-0.832107,  0.176777,  0.832107)+  , (-0.562500,  0.216506,  0.974279)+  , (-0.275566,  0.241481,  1.028426)+  , ( 0.000000,  0.250000,  1.000000)+  , ( 0.242072,  0.241481,  0.903426)+  , ( 0.437500,  0.216506,  0.757772)+  , ( 0.582106,  0.176777,  0.582107)+  , ( 0.678525,  0.125000,  0.391747)+  , ( 0.732673,  0.064705,  0.196319)+  , (-1.207407, -0.000000, -0.323524)+  , (-1.082532, -0.000000, -0.625000)+  , (-0.883884, -0.000000, -0.883884)+  , (-0.625000, -0.000000, -1.082532)+  , (-0.323524, -0.000000, -1.207408)+  , ( 0.000000, -0.000000, -1.250001)+  , ( 0.323524, -0.000000, -1.207408)+  , ( 0.625001, -0.000000, -1.082533)+  , ( 0.883884, -0.000000, -0.883884)+  , ( 1.082533, -0.000000, -0.625000)+  , ( 1.207409, -0.000000, -0.323524)+  , (-0.258819,  0.000000,  0.965926)+  , (-0.562438, -0.099169,  0.732983)+  , ( 0.562437,  0.099169,  0.732983)+  , (-0.707107,  0.000000,  0.707107)+  , ( 0.000000,  0.000000,  1.000000)+  , (-1.030574,  0.047835,  0.426877)+  , ( 1.000000,  0.000000,  0.000000)+  , ( 0.868575,  0.016316,  0.114350)+  , (-1.114315,  0.016316,  0.146703)+  , ( 0.965926,  0.000000,  0.258819)+  , ( 1.030574, -0.047835,  0.426878)+  , (-0.872030,  0.076095,  0.669132)+  , ( 0.872030, -0.076095,  0.669132)+  , ( 0.707106,  0.000000,  0.707107)+  , (-0.400989,  0.115485,  0.968074)+  , ( 0.128396,  0.123931,  0.975269)+  , ( 0.364378,  0.115485,  0.879685)+  , (-0.128397, -0.123931,  0.975268)+  , (-1.115376, -0.000000, -0.146842)+  , (-1.039365, -0.000000, -0.430519)+  , (-0.866026, -0.000000, -0.500000)+  , (-0.694184,  0.000000, -0.532666)+  , (-0.684857, -0.000000, -0.892523)+  , (-0.500000, -0.000000, -0.866026)+  , (-0.334848,  0.000000, -0.808395)+  , (-0.146842, -0.000000, -1.115376)+  , ( 0.000000, -0.000000, -1.000000)+  , ( 0.114211,  0.000000, -0.867515)+  , ( 0.430519, -0.000000, -1.039365)+  , ( 0.500000, -0.000000, -0.866026)+  , ( 0.532667,  0.000000, -0.694185)+  , ( 0.892523, -0.000000, -0.684857)+  , ( 0.866026, -0.000000, -0.500000)+  , ( 0.808395,  0.000000, -0.334848)+  , ( 1.115377, -0.000000, -0.146842)+  ]++mobScales = Right [0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15]+mobiusStrip = makeRipsFiltrationFast mobScales metric3 mobius+nsmobius = simple2Filtr mobiusStrip++tests =+  [ let actual = kernelInt matrix1+        expect = Just $ kernel1+    in checkPass "kernelInt matrix1" actual expect++  , let actual = kernelIntPar matrix1+        expect = Just $ kernel1+    in checkPass "kernelIntPar matrix1" actual expect++  , let actual = kernelInt matrix2+        expect = Just $ kernel2+    in checkPass "kernelInt matrix2" actual expect++  , let actual = kernelIntPar matrix2+        expect = Just $ kernel2+    in checkPass "kernelIntPar matrix2" actual expect++  , let actual = kernelInt matrix3+        expect = Just $ kernel3+    in checkPass "kernelInt matrix3" actual expect++  , let actual = kernelIntPar matrix3+        expect = Just $ kernel3+    in checkPass "kernelIntPar matrix3" actual expect++  , let actual = normalFormInt matrix4+        expect = Just $ snf4+    in checkPass "normalFormInt matrix4" actual expect++  , let actual = normalFormIntPar matrix4+        expect = Just $ snf4+    in checkPass "normalFormIntPar matrix4" actual expect++  , let actual = normalFormInt matrix5+        expect = Just $ snf5+    in checkPass "normalFormInt matrix5" actual expect++  , let actual = normalFormIntPar matrix5+        expect = Just $ snf5+    in checkPass "normalFormIntPar matrix5" actual expect++  , let actual = normalFormInt matrix6+        expect = Just $ snf6+    in checkPass "normalFormInt matrix6" actual expect++  , let actual = normalFormIntPar matrix6+        expect = Just $ snf6+    in checkPass "normalFormIntPar matrix6" actual expect++  --test filtration 1++  , let actual = indexBarCodesSimple testFiltration1+        expect = Just $ V.fromList $ L.map V.fromList [[(0, Finite 2), (0, Finite 2), (0,Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Infinity), (0, Infinity), (0, Infinity)], [(2,Infinity),(2,Infinity)], []]+    in checkPass "indexBarCodesSimple testFiltration1" actual expect++  , let actual = scaleBarCodesSimple scales1 testFiltration1+        expect = Just $ V.fromList $ L.map V.fromList [[(6.0, Finite 10.0), (6.0, Finite 10.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Infinity), (6.0, Infinity), (6.0, Infinity)], [(10.0, Infinity), (10.0, Infinity)], []]+    in checkPass "scaleBarCodesSimple testFiltration1" actual expect++  , let actual = indexBarCodes nsTestFiltration1+        expect = Just $ V.fromList $ L.map V.fromList [[(0, Finite 2), (0, Finite 2), (0,Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Infinity), (0, Infinity), (0, Infinity)], [(2,Infinity),(2,Infinity)], [], []]+    in checkPass "indexBarCodes testFiltration1" actual expect++  , let actual = scaleBarCodes scales1 nsTestFiltration1+        expect = Just $ V.fromList $ L.map V.fromList [[(6.0, Finite 10.0), (6.0, Finite 10.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Finite 8.0), (6.0, Infinity), (6.0, Infinity), (6.0, Infinity)], [(10.0, Infinity), (10.0, Infinity)], [], []]+    in checkPass "scaleBarCodes testFiltration1" actual expect++  --octahedral filtration++  , let actual = indexBarCodesSimple octahedron+        expect = Just $ V.fromList $ L.map V.fromList [[(0, Finite 2), (0, Finite 2), (0, Finite 2), (0, Finite 2), (0,Finite 2),(0,Infinity)], [], [(2,Finite 3)], [], []]+    in checkPass "indexBarCodesSimple octahedron" actual expect++  , let actual = scaleBarCodesSimple octaScales octahedron+        expect = Just $ V.fromList $ L.map V.fromList [[(0.5, Finite 1.6),(0.5, Finite 1.6),(0.5, Finite 1.6),(0.5, Finite 1.6),(0.5, Finite 1.6), (0.5, Infinity)], [], [(1.6,Finite 2.1)], [], []]+    in checkPass "scaleBarCodesSimple octahedron" actual expect++  , let actual = indexBarCodes nsoctahedron+        expect = Just $ V.fromList $ L.map V.fromList [[(0, Finite 2), (0, Finite 2), (0, Finite 2), (0, Finite 2), (0,Finite 2),(0,Infinity)], [], [(2,Finite 3)], [], [], []]+    in checkPass "indexBarCodes octahedron" actual expect++  , let actual = scaleBarCodes octaScales nsoctahedron+        expect = Just $ V.fromList $ L.map V.fromList [[(0.5, Finite 1.6),(0.5, Finite 1.6),(0.5, Finite 1.6),(0.5, Finite 1.6),(0.5, Finite 1.6), (0.5, Infinity)], [], [(1.6,Finite 2.1)], [], [], []]+    in checkPass "scaleBarCodes octahedron" actual expect++  --square filtration++  , let actual = indexBarCodesSimple square+        expect = Just $ V.fromList $ L.map V.fromList [[(0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Infinity)], [(1,Infinity)]]+    in checkPass "indexBarCodesSimple square" actual expect++  , let actual = scaleBarCodesSimple sqrScales square+        expect = Just $ V.fromList $ L.map V.fromList [[(0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Infinity)], [(1.1,Infinity)]]+    in checkPass "scaleBarCodesSimple square" actual expect++  , let actual = indexBarCodes nssquare+        expect = Just $ V.fromList $ L.map V.fromList [[(0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Infinity)], [(1,Infinity)], []]+    in checkPass "indexBarCodes square" actual expect++  , let actual = scaleBarCodes sqrScales nssquare+        expect = Just $ V.fromList $ L.map V.fromList [[(0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Finite 1.1), (0.5, Infinity)], [(1.1,Infinity)], []]+    in checkPass "scaleBarCodes square" actual expect++  --mobius strip filtration++  , let actual = indexBarCodesSimple mobiusStrip+        expect = Just $ V.fromList $ L.map V.fromList [[(0, Finite 4), (0, Finite 4), (0, Finite 4), (0, Finite 4), (0, Finite 4), (0, Finite 2), (0, Finite 2), (0, Finite 2), (0, Finite 2), (0,Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0,Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1),(0,Finite 1),(0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Infinity)], [(4, Finite 5), (3, Finite 5), (3, Finite 5), (4, Infinity)], [], [], []]+    in checkPass "indexBarCodesSimple mobius" actual expect++  , let actual = scaleBarCodesSimple mobScales mobiusStrip+        expect = Just $ V.fromList $ L.map V.fromList $ [[(0.15, Finite 0.35), (0.15, Finite 0.35), (0.15, Finite 0.35), (0.15, Finite 0.35), (0.15, Finite 0.35), (0.15, Finite 0.25), (0.15, Finite 0.25), (0.15, Finite 0.25), (0.15, Finite 0.25), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Infinity)], [(0.35, Finite 0.4), (0.3, Finite 0.4), (0.3, Finite 0.4), (0.35, Infinity)], [], [], []]+    in checkPass "scaleBarCodesSimple mobius" actual expect++  , let actual = indexBarCodes nsmobius+        expect = Just $ V.fromList $ L.map V.fromList [[(0, Finite 4), (0, Finite 4), (0, Finite 4), (0, Finite 4), (0, Finite 4), (0, Finite 2), (0, Finite 2), (0, Finite 2), (0, Finite 2), (0,Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0,Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1),(0,Finite 1),(0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Finite 1), (0, Infinity)], [(4, Finite 5), (3, Finite 5), (3, Finite 5), (4, Infinity)], [], [], [], []]+    in checkPass "indexBarCodes mobius" actual expect++  , let actual = scaleBarCodes mobScales nsmobius+        expect = Just $ V.fromList $ L.map V.fromList $ [[(0.15, Finite 0.35), (0.15, Finite 0.35), (0.15, Finite 0.35), (0.15, Finite 0.35), (0.15, Finite 0.35), (0.15, Finite 0.25), (0.15, Finite 0.25), (0.15, Finite 0.25), (0.15, Finite 0.25), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Finite 0.2), (0.15, Infinity)], [(0.35, Finite 0.4), (0.3, Finite 0.4), (0.3, Finite 0.4), (0.35, Infinity)], [], [], [], []]+    in checkPass "scaleBarCodes mobius" actual expect++  , let complex = getComplex 1 $ Right nssquare+        actual  = bettiNumbers complex+        expect  = Just [1,1]+    in checkPass "bettiNumbers nssquare 1" actual expect++  , let complex = getComplex 1 $ Right nssquare+        actual  = bettiNumbersPar complex+        expect  = Just [1,1]+    in checkPass "bettiNumbersPar nssquare 1" actual expect++  , let complex = getComplex 2 $ Left octahedron+        actual  = bettiNumbers complex+        expect  = Just [1,0,1]+    in checkPass "bettiNumbers octahedron 2" actual expect++  , let complex = getComplex 2 $ Left octahedron+        actual  = bettiNumbersPar complex+        expect  = Just [1,0,1]+    in checkPass "bettiNumbersPar octahedron 2" actual expect++  , let complex = getComplex 2 $ Left testFiltration1+        actual  = bettiNumbers complex+        expect  = Just [3,2,0,0]+    in checkPass "bettiNumbers testFiltration1 2" actual expect++  , let complex = getComplex 2 $ Left testFiltration1+        actual  = bettiNumbersPar complex+        expect  = Just [3,2,0,0]+    in checkPass "bettiNumbersPar testFiltration1 2" actual expect++  --for now we only check if the landscape functions don't throw errors++  , let ibarcodes = indexBarCodesSimple testFiltration1+        actual    = calcLandscape $ ibarcodes ! 0+        expect    = Nothing+    in checkPass "calcLandscape testFiltration1 0-dimensional" actual expect++  , let ibarcodes = indexBarCodesSimple mobiusStrip+        actual    = calcLandscape $ ibarcodes ! 1+        expect    = Nothing+    in checkPass "calcLandscape mobiusStrip 1-dimensional" actual expect++  , let ibarcodestest = indexBarCodesSimple testFiltration1+        ibarcodesmob  = indexBarCodesSimple mobiusStrip+        testscape     = calcLandscape $ ibarcodestest ! 0+        mobscape      = calcLandscape $ ibarcodesmob ! 1+        actual        = avgLandscapes [testscape, mobscape]+        expect        = Nothing+    in checkPass "avgLandscapes testFiltration1 mobiuStrip" actual expect++  , let ibarcodestest = indexBarCodesSimple testFiltration1+        ibarcodesmob  = indexBarCodesSimple mobiusStrip+        testscape     = calcLandscape $ ibarcodestest ! 0+        mobscape      = calcLandscape $ ibarcodesmob ! 1+        actual        = diffLandscapes testscape mobscape+        expect        = Nothing+    in checkPass "diffLandscapes testFiltration1 mobiusStrip" actual expect++  , let ibarcodestest = indexBarCodesSimple testFiltration1+        ibarcodesmob  = indexBarCodesSimple mobiusStrip+        testscape     = calcLandscape $ ibarcodestest ! 0+        mobscape      = calcLandscape $ ibarcodesmob ! 1+        actual        = metricLp (Finite 3.0) (0.0, 5.0) 0.5 testscape mobscape+        expect        = Nothing+    in checkPass "metricLp 3.0 (0.0, 5.0) testFiltration1 mobiusStrip" actual expect++  --Hasse diagram tests++  , let complex = dGraph2sc 18 dGraph1+        actual  = bettiNumbers complex+        expect  = Just [1,25]+    in checkPass "bettiNumbers dGraph2sc dGraph1" actual expect++  , let actual = bettiNumbers dCliqueComplex+        expect = Just [1,8,2,0,0,0,0]+    in checkPass "bettiNumbers dCliquecComplex" actual expect+  ]++main = putStrLn $ formatTestResults tests