{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DerivingStrategies #-}
{-# LANGUAGE ExplicitNamespaces #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}

-- |
-- Module: Data.DiGraph.FloydWarshall
-- Copyright: Copyright © 2018 Kadena LLC.
-- License: MIT
-- Maintainer: Lars Kuhtz <lars@kadena.io>
-- Stability: experimental
--
-- TODO
--
module Data.DiGraph.FloydWarshall
(
-- * Graph Representations
  DenseAdjMatrix
, AdjacencySets

-- * Conversions
, fromAdjacencySets
, toAdjacencySets

-- * FloydWarshall Algorithm
, ShortestPathMatrix(..)
, floydWarshall
, shortestPath
, distance
, diameter

-- * Legacy exports
, distMatrix_
, floydWarshall_
, diameter_
, shortestPaths_
) where

import Control.DeepSeq

import Data.Foldable
import qualified Data.HashMap.Strict as HM
import qualified Data.HashSet as HS
import Data.Massiv.Array as M
#if !MIN_VERSION_base(4,11,0)
import Data.Semigroup
#endif

import GHC.Generics

import Numeric.Natural

-- -------------------------------------------------------------------------- --
-- Graph Representations

-- | Adjacency matrix representation of a directed graph.
--
type DenseAdjMatrix = Array U Ix2 Int

-- | Adjacency set representation of a directed graph.
--
type AdjacencySets = HM.HashMap Int (HS.HashSet Int)

-- -------------------------------------------------------------------------- --
-- Adjacency Matrix for dense Graphs
--
-- (uses massiv)

-- | Assumes that the input is an undirected graph and that the vertex
-- set is a prefix of the natural numbers.
--
fromAdjacencySets :: AdjacencySets -> DenseAdjMatrix
fromAdjacencySets :: AdjacencySets -> DenseAdjMatrix
fromAdjacencySets AdjacencySets
g = Comp -> Sz Ix2 -> (Ix2 -> Int) -> DenseAdjMatrix
forall r ix e.
Construct r ix e =>
Comp -> Sz ix -> (ix -> e) -> Array r ix e
makeArray Comp
Seq (Ix2 -> Sz Ix2
forall ix. Index ix => ix -> Sz ix
Sz (Int
n Int -> Int -> Ix2
:. Int
n)) Ix2 -> Int
forall p. Num p => Ix2 -> p
go
  where
    n :: Int
n = AdjacencySets -> Int
forall k v. HashMap k v -> Int
HM.size AdjacencySets
g
    go :: Ix2 -> p
go (Int
i :. Int
j)
        | (Int, Int) -> Bool
isEdge (Int
i, Int
j) = p
1
        | (Int, Int) -> Bool
isEdge (Int
j, Int
i) = p
1
        | Bool
otherwise = p
0
    isEdge :: (Int, Int) -> Bool
isEdge (Int
a, Int
b) = Bool -> (HashSet Int -> Bool) -> Maybe (HashSet Int) -> Bool
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Bool
False (Int -> HashSet Int -> Bool
forall a. (Eq a, Hashable a) => a -> HashSet a -> Bool
HS.member Int
b) (Maybe (HashSet Int) -> Bool) -> Maybe (HashSet Int) -> Bool
forall a b. (a -> b) -> a -> b
$ Int -> AdjacencySets -> Maybe (HashSet Int)
forall k v. (Eq k, Hashable k) => k -> HashMap k v -> Maybe v
HM.lookup Int
a AdjacencySets
g

-- | Converts an adjacency matrix into a graph in adjacnency set representation.
--
toAdjacencySets :: DenseAdjMatrix -> AdjacencySets
toAdjacencySets :: DenseAdjMatrix -> AdjacencySets
toAdjacencySets = (AdjacencySets -> Ix2 -> Int -> AdjacencySets)
-> AdjacencySets -> DenseAdjMatrix -> AdjacencySets
forall r ix e a.
Source r ix e =>
(a -> ix -> e -> a) -> a -> Array r ix e -> a
ifoldlS AdjacencySets -> Ix2 -> Int -> AdjacencySets
forall a.
(Eq a, Num a) =>
AdjacencySets -> Ix2 -> a -> AdjacencySets
f AdjacencySets
forall a. Monoid a => a
mempty
  where
    f :: AdjacencySets -> Ix2 -> a -> AdjacencySets
f AdjacencySets
a (Int
i :. Int
j) a
x
        | a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = AdjacencySets
a
        | Bool
otherwise = (HashSet Int -> HashSet Int -> HashSet Int)
-> Int -> HashSet Int -> AdjacencySets -> AdjacencySets
forall k v.
(Eq k, Hashable k) =>
(v -> v -> v) -> k -> v -> HashMap k v -> HashMap k v
HM.insertWith HashSet Int -> HashSet Int -> HashSet Int
forall a. Semigroup a => a -> a -> a
(<>) Int
i (Int -> HashSet Int
forall a. Hashable a => a -> HashSet a
HS.singleton Int
j) AdjacencySets
a

-- -------------------------------------------------------------------------- --
-- Floyd Warshall with Paths

-- | Shortest path matrix of a graph.
--
newtype ShortestPathMatrix = ShortestPathMatrix (Array U Ix2 (Double, Int))
    deriving (Int -> ShortestPathMatrix -> ShowS
[ShortestPathMatrix] -> ShowS
ShortestPathMatrix -> String
(Int -> ShortestPathMatrix -> ShowS)
-> (ShortestPathMatrix -> String)
-> ([ShortestPathMatrix] -> ShowS)
-> Show ShortestPathMatrix
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [ShortestPathMatrix] -> ShowS
$cshowList :: [ShortestPathMatrix] -> ShowS
show :: ShortestPathMatrix -> String
$cshow :: ShortestPathMatrix -> String
showsPrec :: Int -> ShortestPathMatrix -> ShowS
$cshowsPrec :: Int -> ShortestPathMatrix -> ShowS
Show, ShortestPathMatrix -> ShortestPathMatrix -> Bool
(ShortestPathMatrix -> ShortestPathMatrix -> Bool)
-> (ShortestPathMatrix -> ShortestPathMatrix -> Bool)
-> Eq ShortestPathMatrix
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
$c/= :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
== :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
$c== :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
Eq, Eq ShortestPathMatrix
Eq ShortestPathMatrix
-> (ShortestPathMatrix -> ShortestPathMatrix -> Ordering)
-> (ShortestPathMatrix -> ShortestPathMatrix -> Bool)
-> (ShortestPathMatrix -> ShortestPathMatrix -> Bool)
-> (ShortestPathMatrix -> ShortestPathMatrix -> Bool)
-> (ShortestPathMatrix -> ShortestPathMatrix -> Bool)
-> (ShortestPathMatrix -> ShortestPathMatrix -> ShortestPathMatrix)
-> (ShortestPathMatrix -> ShortestPathMatrix -> ShortestPathMatrix)
-> Ord ShortestPathMatrix
ShortestPathMatrix -> ShortestPathMatrix -> Bool
ShortestPathMatrix -> ShortestPathMatrix -> Ordering
ShortestPathMatrix -> ShortestPathMatrix -> ShortestPathMatrix
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: ShortestPathMatrix -> ShortestPathMatrix -> ShortestPathMatrix
$cmin :: ShortestPathMatrix -> ShortestPathMatrix -> ShortestPathMatrix
max :: ShortestPathMatrix -> ShortestPathMatrix -> ShortestPathMatrix
$cmax :: ShortestPathMatrix -> ShortestPathMatrix -> ShortestPathMatrix
>= :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
$c>= :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
> :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
$c> :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
<= :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
$c<= :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
< :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
$c< :: ShortestPathMatrix -> ShortestPathMatrix -> Bool
compare :: ShortestPathMatrix -> ShortestPathMatrix -> Ordering
$ccompare :: ShortestPathMatrix -> ShortestPathMatrix -> Ordering
$cp1Ord :: Eq ShortestPathMatrix
Ord, (forall x. ShortestPathMatrix -> Rep ShortestPathMatrix x)
-> (forall x. Rep ShortestPathMatrix x -> ShortestPathMatrix)
-> Generic ShortestPathMatrix
forall x. Rep ShortestPathMatrix x -> ShortestPathMatrix
forall x. ShortestPathMatrix -> Rep ShortestPathMatrix x
forall a.
(forall x. a -> Rep a x) -> (forall x. Rep a x -> a) -> Generic a
$cto :: forall x. Rep ShortestPathMatrix x -> ShortestPathMatrix
$cfrom :: forall x. ShortestPathMatrix -> Rep ShortestPathMatrix x
Generic)
    deriving newtype (ShortestPathMatrix -> ()
(ShortestPathMatrix -> ()) -> NFData ShortestPathMatrix
forall a. (a -> ()) -> NFData a
rnf :: ShortestPathMatrix -> ()
$crnf :: ShortestPathMatrix -> ()
NFData)

-- | Shortest path computation for integral matrixes.
--
floydWarshall :: Unbox a => Real a => Array U Ix2 a -> ShortestPathMatrix
floydWarshall :: Array U Ix2 a -> ShortestPathMatrix
floydWarshall = Array U Ix2 (Double, Int) -> ShortestPathMatrix
ShortestPathMatrix
    (Array U Ix2 (Double, Int) -> ShortestPathMatrix)
-> (Array U Ix2 a -> Array U Ix2 (Double, Int))
-> Array U Ix2 a
-> ShortestPathMatrix
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Array U Ix2 (Double, Int) -> Array U Ix2 (Double, Int)
floydWarshallInternal
    (Array U Ix2 (Double, Int) -> Array U Ix2 (Double, Int))
-> (Array U Ix2 a -> Array U Ix2 (Double, Int))
-> Array U Ix2 a
-> Array U Ix2 (Double, Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. U -> Array D Ix2 (Double, Int) -> Array U Ix2 (Double, Int)
forall r ix e r'.
(Mutable r ix e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
computeAs U
U
    (Array D Ix2 (Double, Int) -> Array U Ix2 (Double, Int))
-> (Array U Ix2 a -> Array D Ix2 (Double, Int))
-> Array U Ix2 a
-> Array U Ix2 (Double, Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Array U Ix2 a -> Array D Ix2 (Double, Int)
forall a r.
(Real a, Source r Ix2 a) =>
Array r Ix2 a -> Array D Ix2 (Double, Int)
intDistMatrix

-- | Compute a shortest path between two vertices of a graph from the shortest
-- path matrix of the graph.
--
shortestPath
    :: ShortestPathMatrix
    -> Int
    -> Int
    -> Maybe [Int]
shortestPath :: ShortestPathMatrix -> Int -> Int -> Maybe [Int]
shortestPath (ShortestPathMatrix Array U Ix2 (Double, Int)
m) Int
src Int
trg
    | DenseAdjMatrix -> Bool
forall r ix e. Load r ix e => Array r ix e -> Bool
M.isEmpty DenseAdjMatrix
mat = Maybe [Int]
forall a. Maybe a
Nothing
    | Bool -> Bool
not (Sz Ix2 -> Ix2 -> Bool
forall ix. Index ix => Sz ix -> ix -> Bool
M.isSafeIndex (Array U Ix2 (Double, Int) -> Sz Ix2
forall r ix e. Load r ix e => Array r ix e -> Sz ix
size Array U Ix2 (Double, Int)
m) (Int
src Int -> Int -> Ix2
:. Int
trg)) = Maybe [Int]
forall a. Maybe a
Nothing
    | (DenseAdjMatrix
mat DenseAdjMatrix -> Ix2 -> Int
forall r ix e. Manifest r ix e => Array r ix e -> ix -> e
M.! (Int
src Int -> Int -> Ix2
:. Int
trg)) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== (-Int
1) = Maybe [Int]
forall a. Maybe a
Nothing
    | Bool
otherwise = Int -> Int -> Maybe [Int]
go Int
src Int
trg
  where
    mat :: DenseAdjMatrix
mat = U -> Array D Ix2 Int -> DenseAdjMatrix
forall r ix e r'.
(Mutable r ix e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
M.computeAs U
U (Array D Ix2 Int -> DenseAdjMatrix)
-> Array D Ix2 Int -> DenseAdjMatrix
forall a b. (a -> b) -> a -> b
$ ((Double, Int) -> Int)
-> Array U Ix2 (Double, Int) -> Array D Ix2 Int
forall r ix e' e.
Source r ix e' =>
(e' -> e) -> Array r ix e' -> Array D ix e
M.map (Double, Int) -> Int
forall a b. (a, b) -> b
snd Array U Ix2 (Double, Int)
m
    go :: Int -> Int -> Maybe [Int]
go Int
a Int
b
        | Int
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
b = [Int] -> Maybe [Int]
forall (m :: * -> *) a. Monad m => a -> m a
return []
        | Bool
otherwise = do
            Int
n <- DenseAdjMatrix -> Ix2 -> Maybe Int
forall r ix e. Manifest r ix e => Array r ix e -> ix -> Maybe e
M.index DenseAdjMatrix
mat (Int
a Int -> Int -> Ix2
:. Int
b)
            (:) Int
n ([Int] -> [Int]) -> Maybe [Int] -> Maybe [Int]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> Int -> Maybe [Int]
go Int
n Int
b

-- | Compute the distance between two vertices of a graph from the shortest path
-- matrix of the graph.
--
distance :: ShortestPathMatrix -> Int -> Int -> Maybe Double
distance :: ShortestPathMatrix -> Int -> Int -> Maybe Double
distance (ShortestPathMatrix Array U Ix2 (Double, Int)
m) Int
src Int
trg
    | Array U Ix2 (Double, Int) -> Bool
forall r ix e. Load r ix e => Array r ix e -> Bool
M.isEmpty Array U Ix2 (Double, Int)
m = Maybe Double
forall a. Maybe a
Nothing
    | Bool
otherwise = Double -> Maybe Double
forall a. RealFrac a => a -> Maybe a
toDistance (Double -> Maybe Double)
-> ((Double, Int) -> Double) -> (Double, Int) -> Maybe Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double, Int) -> Double
forall a b. (a, b) -> a
fst ((Double, Int) -> Maybe Double)
-> Maybe (Double, Int) -> Maybe Double
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<< Array U Ix2 (Double, Int) -> Ix2 -> Maybe (Double, Int)
forall r ix e. Manifest r ix e => Array r ix e -> ix -> Maybe e
M.index Array U Ix2 (Double, Int)
m (Int
src Int -> Int -> Ix2
:. Int
trg)

-- | Compute the diameter of a graph from the shortest path matrix of the graph.
--
diameter :: ShortestPathMatrix -> Maybe Double
diameter :: ShortestPathMatrix -> Maybe Double
diameter (ShortestPathMatrix Array U Ix2 (Double, Int)
m)
    | Array U Ix2 (Double, Int) -> Bool
forall r ix e. Load r ix e => Array r ix e -> Bool
M.isEmpty Array U Ix2 (Double, Int)
m = Double -> Maybe Double
forall a. a -> Maybe a
Just Double
0
    | Bool
otherwise = Double -> Maybe Double
forall a. RealFrac a => a -> Maybe a
toDistance (Double -> Maybe Double) -> Double -> Maybe Double
forall a b. (a -> b) -> a -> b
$ Array D Ix2 Double -> Double
forall r ix e. (Source r ix e, Ord e) => Array r ix e -> e
maximum' (Array D Ix2 Double -> Double) -> Array D Ix2 Double -> Double
forall a b. (a -> b) -> a -> b
$ ((Double, Int) -> Double)
-> Array U Ix2 (Double, Int) -> Array D Ix2 Double
forall r ix e' e.
Source r ix e' =>
(e' -> e) -> Array r ix e' -> Array D ix e
M.map (Double, Int) -> Double
forall a b. (a, b) -> a
fst Array U Ix2 (Double, Int)
m

-- -------------------------------------------------------------------------- --
-- Internal

toDistance :: RealFrac a => a -> Maybe a
toDistance :: a -> Maybe a
toDistance a
x
    | a
x a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
0 = Maybe a
forall a. Maybe a
Nothing
    | Bool
otherwise = a -> Maybe a
forall a. a -> Maybe a
Just a
x

-- | Distance matrix for int inputs.
--
intDistMatrix
    :: Real a
    => Source r Ix2 a
    => Array r Ix2 a
    -> Array M.D Ix2 (Double, Int)
intDistMatrix :: Array r Ix2 a -> Array D Ix2 (Double, Int)
intDistMatrix = (Ix2 -> a -> (Double, Int))
-> Array r Ix2 a -> Array D Ix2 (Double, Int)
forall r ix e' e.
Source r ix e' =>
(ix -> e' -> e) -> Array r ix e' -> Array D ix e
M.imap Ix2 -> a -> (Double, Int)
forall a a. (Real a, Fractional a) => Ix2 -> a -> (a, Int)
go
  where
    go :: Ix2 -> a -> (a, Int)
go (Int
x :. Int
y) a
e
        | Int
x Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
y = (a
0, Int
y)
        | a
e a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 = (a -> a
forall a b. (Real a, Fractional b) => a -> b
realToFrac a
e, Int
y)
        | Bool
otherwise = (a
1a -> a -> a
forall a. Fractional a => a -> a -> a
/a
0, -Int
1)

-- | Floyd-Warshall With Path Matrix
--
-- TODO: use a mutable array?
-- TODO: implement Dijkstra's algorithm for adj matrix representation.
--
floydWarshallInternal
    :: Array U Ix2 (Double, Int)
    -> Array U Ix2 (Double,Int)
floydWarshallInternal :: Array U Ix2 (Double, Int) -> Array U Ix2 (Double, Int)
floydWarshallInternal Array U Ix2 (Double, Int)
a = (Array U Ix2 (Double, Int) -> Int -> Array U Ix2 (Double, Int))
-> Array U Ix2 (Double, Int) -> [Int] -> Array U Ix2 (Double, Int)
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Array U Ix2 (Double, Int) -> Int -> Array U Ix2 (Double, Int)
go Array U Ix2 (Double, Int)
a [Int
0..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
  where
    Sz (Int
n :. Int
_) = Array U Ix2 (Double, Int) -> Sz Ix2
forall r ix e. Load r ix e => Array r ix e -> Sz ix
size Array U Ix2 (Double, Int)
a

    go :: Array U Ix2 (Double, Int) -> Int -> Array U Ix2 (Double,Int)
    go :: Array U Ix2 (Double, Int) -> Int -> Array U Ix2 (Double, Int)
go Array U Ix2 (Double, Int)
c Int
k = Comp
-> Sz Ix2 -> (Ix2 -> (Double, Int)) -> Array U Ix2 (Double, Int)
forall r ix e.
Construct r ix e =>
Comp -> Sz ix -> (ix -> e) -> Array r ix e
makeArray Comp
Seq (Ix2 -> Sz Ix2
forall ix. Index ix => ix -> Sz ix
Sz (Int
n Int -> Int -> Ix2
:. Int
n)) ((Ix2 -> (Double, Int)) -> Array U Ix2 (Double, Int))
-> (Ix2 -> (Double, Int)) -> Array U Ix2 (Double, Int)
forall a b. (a -> b) -> a -> b
$ \(Int
x :. Int
y) ->
        let
            !xy :: Double
xy = (Double, Int) -> Double
forall a b. (a, b) -> a
fst ((Double, Int) -> Double) -> (Double, Int) -> Double
forall a b. (a -> b) -> a -> b
$! Array U Ix2 (Double, Int)
c Array U Ix2 (Double, Int) -> Ix2 -> (Double, Int)
forall r ix e. Manifest r ix e => Array r ix e -> ix -> e
M.! (Int
x Int -> Int -> Ix2
:. Int
y)
            !xk :: Double
xk = (Double, Int) -> Double
forall a b. (a, b) -> a
fst ((Double, Int) -> Double) -> (Double, Int) -> Double
forall a b. (a -> b) -> a -> b
$! Array U Ix2 (Double, Int)
c Array U Ix2 (Double, Int) -> Ix2 -> (Double, Int)
forall r ix e. Manifest r ix e => Array r ix e -> ix -> e
M.! (Int
x Int -> Int -> Ix2
:. Int
k)
            !ky :: Double
ky = (Double, Int) -> Double
forall a b. (a, b) -> a
fst ((Double, Int) -> Double) -> (Double, Int) -> Double
forall a b. (a -> b) -> a -> b
$! Array U Ix2 (Double, Int)
c Array U Ix2 (Double, Int) -> Ix2 -> (Double, Int)
forall r ix e. Manifest r ix e => Array r ix e -> ix -> e
M.! (Int
k Int -> Int -> Ix2
:. Int
y)
            !nxy :: Int
nxy = (Double, Int) -> Int
forall a b. (a, b) -> b
snd ((Double, Int) -> Int) -> (Double, Int) -> Int
forall a b. (a -> b) -> a -> b
$! Array U Ix2 (Double, Int)
c Array U Ix2 (Double, Int) -> Ix2 -> (Double, Int)
forall r ix e. Manifest r ix e => Array r ix e -> ix -> e
M.! (Int
x Int -> Int -> Ix2
:. Int
y)
            !nxk :: Int
nxk = (Double, Int) -> Int
forall a b. (a, b) -> b
snd ((Double, Int) -> Int) -> (Double, Int) -> Int
forall a b. (a -> b) -> a -> b
$! Array U Ix2 (Double, Int)
c Array U Ix2 (Double, Int) -> Ix2 -> (Double, Int)
forall r ix e. Manifest r ix e => Array r ix e -> ix -> e
M.! (Int
x Int -> Int -> Ix2
:. Int
k)
        in if Double
xy Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
xk Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ky then (Double
xk Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ky, Int
nxk) else (Double
xy, Int
nxy)

-- -------------------------------------------------------------------------- --
-- Floyd Warshall Without Paths (more efficient, by factor of 2)

-- | Floyd Warshall Without Paths (more efficient, by factor of 2).
--
distMatrix_
    :: Source r Ix2 Int
    => Array r Ix2 Int
    -> Array M.D Ix2 Double
distMatrix_ :: Array r Ix2 Int -> Array D Ix2 Double
distMatrix_ = (Ix2 -> Int -> Double) -> Array r Ix2 Int -> Array D Ix2 Double
forall r ix e' e.
Source r ix e' =>
(ix -> e' -> e) -> Array r ix e' -> Array D ix e
M.imap Ix2 -> Int -> Double
forall a p. (Real a, Fractional p) => Ix2 -> a -> p
go
  where
    go :: Ix2 -> a -> p
go (Int
x :. Int
y) a
e
        | Int
x Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
y = p
0
        | a
e a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0 = a -> p
forall a b. (Real a, Fractional b) => a -> b
realToFrac a
e
        | Bool
otherwise = p
1p -> p -> p
forall a. Fractional a => a -> a -> a
/p
0

-- | Floyd Warshall Without Paths (more efficient, by factor of 2).
--
-- TODO: use a mutable array?
-- TODO: implement Dijkstra's algorithm for adj matrix representation.
--
floydWarshall_
    :: Array U Ix2 Double
    -> Array U Ix2 Double
floydWarshall_ :: Array U Ix2 Double -> Array U Ix2 Double
floydWarshall_ Array U Ix2 Double
a = (Array U Ix2 Double -> Int -> Array U Ix2 Double)
-> Array U Ix2 Double -> [Int] -> Array U Ix2 Double
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' Array U Ix2 Double -> Int -> Array U Ix2 Double
go Array U Ix2 Double
a [Int
0..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]
  where
    Sz (Int
n :. Int
_) = Array U Ix2 Double -> Sz Ix2
forall r ix e. Load r ix e => Array r ix e -> Sz ix
size Array U Ix2 Double
a

    go :: Array U Ix2 Double -> Int -> Array U Ix2 Double
    go :: Array U Ix2 Double -> Int -> Array U Ix2 Double
go Array U Ix2 Double
c Int
k = Comp -> Sz Ix2 -> (Ix2 -> Double) -> Array U Ix2 Double
forall r ix e.
Construct r ix e =>
Comp -> Sz ix -> (ix -> e) -> Array r ix e
makeArray Comp
Seq (Ix2 -> Sz Ix2
forall ix. Index ix => ix -> Sz ix
Sz (Int
n Int -> Int -> Ix2
:. Int
n)) ((Ix2 -> Double) -> Array U Ix2 Double)
-> (Ix2 -> Double) -> Array U Ix2 Double
forall a b. (a -> b) -> a -> b
$ \(Int
x :. Int
y) ->
        let
            !xy :: Double
xy = Array U Ix2 Double
c Array U Ix2 Double -> Ix2 -> Double
forall r ix e. Manifest r ix e => Array r ix e -> ix -> e
M.! (Int
x Int -> Int -> Ix2
:. Int
y)
            !xk :: Double
xk = Array U Ix2 Double
c Array U Ix2 Double -> Ix2 -> Double
forall r ix e. Manifest r ix e => Array r ix e -> ix -> e
M.! (Int
x Int -> Int -> Ix2
:. Int
k)
            !ky :: Double
ky = Array U Ix2 Double
c Array U Ix2 Double -> Ix2 -> Double
forall r ix e. Manifest r ix e => Array r ix e -> ix -> e
M.! (Int
k Int -> Int -> Ix2
:. Int
y)
        in if Double
xy Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
xk Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ky then Double
xk Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ky else Double
xy

-- | Shortest path matrix.
--
-- All entries of the result matrix are either whole numbers or @Infinity@.
--
shortestPaths_ :: Array U Ix2 Int -> Array U Ix2 Double
shortestPaths_ :: DenseAdjMatrix -> Array U Ix2 Double
shortestPaths_ = Array U Ix2 Double -> Array U Ix2 Double
floydWarshall_ (Array U Ix2 Double -> Array U Ix2 Double)
-> (DenseAdjMatrix -> Array U Ix2 Double)
-> DenseAdjMatrix
-> Array U Ix2 Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. U -> Array D Ix2 Double -> Array U Ix2 Double
forall r ix e r'.
(Mutable r ix e, Load r' ix e) =>
r -> Array r' ix e -> Array r ix e
computeAs U
U (Array D Ix2 Double -> Array U Ix2 Double)
-> (DenseAdjMatrix -> Array D Ix2 Double)
-> DenseAdjMatrix
-> Array U Ix2 Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. DenseAdjMatrix -> Array D Ix2 Double
forall r. Source r Ix2 Int => Array r Ix2 Int -> Array D Ix2 Double
distMatrix_

-- | Diameter of a graph.
--
diameter_ :: Array U Ix2 Int -> Maybe Natural
diameter_ :: DenseAdjMatrix -> Maybe Natural
diameter_ DenseAdjMatrix
g
    | DenseAdjMatrix -> Bool
forall r ix e. Load r ix e => Array r ix e -> Bool
M.isEmpty DenseAdjMatrix
g = Natural -> Maybe Natural
forall a. a -> Maybe a
Just Natural
0
    | Bool
otherwise = let x :: Natural
x = Double -> Natural
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Natural) -> Double -> Natural
forall a b. (a -> b) -> a -> b
$ Array U Ix2 Double -> Double
forall r ix e. (Source r ix e, Ord e) => Array r ix e -> e
maximum' (Array U Ix2 Double -> Double) -> Array U Ix2 Double -> Double
forall a b. (a -> b) -> a -> b
$ DenseAdjMatrix -> Array U Ix2 Double
shortestPaths_ DenseAdjMatrix
g
        in if Natural
x Natural -> Natural -> Bool
forall a. Eq a => a -> a -> Bool
== Double -> Natural
forall a b. (RealFrac a, Integral b) => a -> b
round (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0 :: Double) then Maybe Natural
forall a. Maybe a
Nothing else Natural -> Maybe Natural
forall a. a -> Maybe a
Just Natural
x