{- |
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