Persistence (empty) → 1.0

raw patch · 8 files changed

+1171/−0 lines, 8 filesdep +basedep +containersdep +maximal-cliquessetup-changed

Dependencies added: base, containers, maximal-cliques, parallel, vector

Files

• ChangeLog.md +5/−0
• LICENSE +30/−0
• Matrix.hs +271/−0
• Persistence.cabal +24/−0
• Persistence.hs +253/−0
• Setup.hs +2/−0
• SimplicialComplex.hs +271/−0
• Util.hs +315/−0

@@ -0,0 +1,5 @@+# Revision history for Persistence++## 1.0  -- YYYY-mm-dd++* First version. Released on an unsuspecting world.

@@ -0,0 +1,30 @@+Copyright (c) 2018, Eben Cowley++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.++    * Redistributions in binary form must reproduce the above+      copyright notice, this list of conditions and the following+      disclaimer in the documentation and/or other materials provided+      with the distribution.++    * Neither the name of Eben Cowley nor the names of other+      contributors may be used to endorse or promote products derived+      from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

@@ -0,0 +1,271 @@+{- |+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+  , getDiagonal+  , getUnsignedDiagonal+  , transposeMat+  , transposePar+  , multiply+  , multiplyPar+  , 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 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)++-- | 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 --the following variable should be useful for quickly determining whether or not there are more entries to eleiminate+      indices = V.filter (\index -> index > colIndex) $ V.findIndices id row --but that method is not working for some reason+  in+    if not $ row ! colIndex then+      case indices of+        v | V.null v -> Nothing+        v            ->+          let j = V.head v+          in Just (V.length v > 0, V.map (switchElems colIndex j) mat, Just (colIndex, j))+    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

@@ -0,0 +1,24 @@+-- Initial Persistence.cabal generated by cabal init.  For further +-- documentation, see http://haskell.org/cabal/users-guide/++name:                Persistence+version:             1.0+synopsis:            Quickly detect clusters and holes in data.+description:         Persistence is a topological data analysis library motivated by flexibility when it comes to the type of data being analyzed. If you have data that comes with a meaningful function into something of the Ord typeclass, you can use Persistence to detect clusters and holes in the data. You can also use the library to analyze undirected graphs, and interesting features for directed graphs will be added soon.+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++library+  exposed-modules:     SimplicialComplex, Persistence+  other-modules:       Util, Matrix+  -- other-extensions:    +  build-depends:       base >=4.10 && <4.11, vector >=0.12 && <0.13, parallel >=3.2 && <3.3, containers >=0.5 && <0.6, maximal-cliques >=0.1 && <0.2+  -- hs-source-dirs:      +  default-language:    Haskell2010

@@ -0,0 +1,253 @@+{- |+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.++-}++module Persistence+  ( Filtration+  , BarCode+  , sim2String+  , filtr2String+  , getComplex+  , getNumSimplices+  , getDimension+  , filterByWeightsFast+  , makeVRFiltrationFast+  , filterByWeightsLight+  , makeVRFiltrationLight+  , persistentHomology+  ) 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)++-- | 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.+  Really, this is a helper function to be called by makeVRFiltrationFast, but I decided to expose it in case you have a simplicial complex and weighted graph lying around.+  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)+              (if begin == i then codes else (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) (if i == 0 then codes else (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)

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

@@ -0,0 +1,271 @@+{- |+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.++The clique complex of a graph is a simplicial complex whose simplices are the complete subgraphs of the given graph. The Vietoris-Rips complex is the clique complex of the neighborhood graph.++Betti numbers measure the number of n-dimensional holes in a given simplicial complex. The 0th Betti number is the number of connected components, or clusters; the 1st Betti number is the number of loops or tunnels; then 2nd Betti number measures voids or hollow volumes; and if you have high-dimensional data, you might have higher Betti numbers representing the number of high-dimensional holes.++-}++module SimplicialComplex+  ( SimplicialComplex+  , Graph+  , sc2String+  , getDim+  , 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.+  The dimension starts at 1 because there is no need to store individual vertices.+  A simplex is represented as a pair: the vector of its vertices in the original data set from which the complex was constructed,+  and the vector of the indices of the faces in the next lowest dimension. Edges are the exception - they do not store reference to their+  faces because it would be redundant with their vertices.+-}+type SimplicialComplex = (Int, Vector (Vector (Vector Int, Vector Int)))++-- | 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++{- |+  This represents the (symmetric) adjacency matrix of some undirected graph. The type `a` is whatever distance is in your data analysis regime.+  The reason graphs are represented like this is because their main function is too speed up the construction of simplicial complexes and filtrations.+  If the clique complex of this graph were to be constructed, only the adjacencies would matter. But if a filtration is constructed all the distances+  will be required over and over again - this allows them to be accessed in constant time.+-}+type Graph a = Vector (Vector (a, Bool))++{- |+  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

@@ -0,0 +1,315 @@+{- |+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 calc i result =+        let a = vector1 ! i+        in case V.find (\b -> rel a b) vector2 of+          Just b  -> calc (i + 1) $ result `snoc` (a, b)+          Nothing -> calc (i + 1) result+  in 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++-- | 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)