{-# LANGUAGE TypeOperators #-}
module Math.LinearCircuit (resistance) where

import qualified Data.Graph.Comfort as Graph
import Data.Graph.Comfort (Graph)

import qualified Numeric.LAPACK.Matrix.Shape as MatrixShape
import qualified Numeric.LAPACK.Matrix.Symmetric as Symmetric
import qualified Numeric.LAPACK.Matrix as Matrix
import qualified Numeric.LAPACK.Vector as Vector
import qualified Numeric.Netlib.Class as Class
import Numeric.LAPACK.Matrix.Symmetric ((#%%%#))
import Numeric.LAPACK.Matrix ((#\|))

import qualified Data.Array.Comfort.Boxed as BoxedArray
import qualified Data.Array.Comfort.Storable as Array
import Data.Array.Comfort.Storable ((!))
import Data.Array.Comfort.Shape ((::+)((::+)))

import Control.Monad.Trans.Identity (IdentityT(IdentityT))

import qualified Data.Map as Map
import Data.Set (Set)


type Wrap = IdentityT

voltageMatrix ::
   (Graph.Edge edge, Ord node, Class.Floating a) =>
   Set node -> Set (Wrap edge node) ->
   Matrix.General (Set (Wrap edge node)) (Set node) a
voltageMatrix :: Set node
-> Set (Wrap edge node)
-> General (Set (Wrap edge node)) (Set node) a
voltageMatrix Set node
nodes =
   Set node
-> Array (Set (Wrap edge node)) (Vector (Set node) a)
-> General (Set (Wrap edge node)) (Set node) a
forall height width a.
(C height, C width, Eq width, Storable a) =>
width -> Array height (Vector width a) -> General height width a
Matrix.fromRowArray Set node
nodes (Array (Set (Wrap edge node)) (Vector (Set node) a)
 -> General (Set (Wrap edge node)) (Set node) a)
-> (Set (Wrap edge node)
    -> Array (Set (Wrap edge node)) (Vector (Set node) a))
-> Set (Wrap edge node)
-> General (Set (Wrap edge node)) (Set node) a
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   (Wrap edge node -> Vector (Set node) a)
-> Array (Set (Wrap edge node)) (Wrap edge node)
-> Array (Set (Wrap edge node)) (Vector (Set node) a)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap
      (\(IdentityT edge node
e) ->
         a -> Set node -> [(Index (Set node), a)] -> Vector (Set node) a
forall sh a.
(Indexed sh, Storable a) =>
a -> sh -> [(Index sh, a)] -> Array sh a
Array.fromAssociations a
0 Set node
nodes
            [(edge node -> node
forall (edge :: * -> *) node. Edge edge => edge node -> node
Graph.from edge node
e, a
1), (edge node -> node
forall (edge :: * -> *) node. Edge edge => edge node -> node
Graph.to edge node
e, -a
1)]) (Array (Set (Wrap edge node)) (Wrap edge node)
 -> Array (Set (Wrap edge node)) (Vector (Set node) a))
-> (Set (Wrap edge node)
    -> Array (Set (Wrap edge node)) (Wrap edge node))
-> Set (Wrap edge node)
-> Array (Set (Wrap edge node)) (Vector (Set node) a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
   Set (Wrap edge node)
-> Array (Set (Wrap edge node)) (Wrap edge node)
forall sh. Indexed sh => sh -> Array sh (Index sh)
BoxedArray.indices

fullMatrix ::
   (Graph.Edge edge, Ord node, Class.Floating a) =>
   Graph edge node a nodeLabel -> node ->
   Matrix.Symmetric (() ::+ Set (Wrap edge node) ::+ Set node) a
fullMatrix :: Graph edge node a nodeLabel
-> node -> Symmetric (() ::+ (Set (Wrap edge node) ::+ Set node)) a
fullMatrix Graph edge node a nodeLabel
gr node
src =
   let edges :: Map (Wrap edge node) a
edges = (edge node -> Wrap edge node)
-> Map (edge node) a -> Map (Wrap edge node) a
forall k1 k2 a. (k1 -> k2) -> Map k1 a -> Map k2 a
Map.mapKeysMonotonic edge node -> Wrap edge node
forall k (f :: k -> *) (a :: k). f a -> IdentityT f a
IdentityT (Map (edge node) a -> Map (Wrap edge node) a)
-> Map (edge node) a -> Map (Wrap edge node) a
forall a b. (a -> b) -> a -> b
$ Graph edge node a nodeLabel -> Map (edge node) a
forall (edge :: * -> *) node edgeLabel nodeLabel.
(Edge edge, Ord node) =>
Graph edge node edgeLabel nodeLabel -> Map (edge node) edgeLabel
Graph.edgeLabels Graph edge node a nodeLabel
gr
       nodes :: Set node
nodes = Graph edge node a nodeLabel -> Set node
forall (e :: * -> *) n el nl. Graph e n el nl -> Set n
Graph.nodeSet Graph edge node a nodeLabel
gr
       order :: Order
order = Order
MatrixShape.RowMajor
       symmetricZero :: width
-> ArrayMatrix
     Packed Symmetric Filled Filled Shape Small Small width width a
symmetricZero width
sh = Omni Packed Symmetric Filled Filled Shape Small Small width width
-> ArrayMatrix
     Packed Symmetric Filled Filled Shape Small Small width width a
forall property meas vert horiz height width a pack lower upper.
(Homogeneous property, Measure meas, C vert, C horiz, C height,
 C width, Floating a) =>
Omni pack property lower upper meas vert horiz height width
-> ArrayMatrix
     pack property lower upper meas vert horiz height width a
Matrix.zero (Omni Packed Symmetric Filled Filled Shape Small Small width width
 -> ArrayMatrix
      Packed Symmetric Filled Filled Shape Small Small width width a)
-> Omni
     Packed Symmetric Filled Filled Shape Small Small width width
-> ArrayMatrix
     Packed Symmetric Filled Filled Shape Small Small width width a
forall a b. (a -> b) -> a -> b
$ Order
-> width
-> Omni
     Packed Symmetric Filled Filled Shape Small Small width width
forall size. Order -> size -> Symmetric size
MatrixShape.symmetric Order
order width
sh
       unit :: Vector (Set (Wrap edge node) ::+ Set node) a
unit = (Set (Wrap edge node) ::+ Set node)
-> Index (Set (Wrap edge node) ::+ Set node)
-> Vector (Set (Wrap edge node) ::+ Set node) a
forall sh a.
(Indexed sh, Floating a) =>
sh -> Index sh -> Vector sh a
Vector.unit (Map (Wrap edge node) a -> Set (Wrap edge node)
forall k a. Map k a -> Set k
Map.keysSet Map (Wrap edge node) a
edges Set (Wrap edge node)
-> Set node -> Set (Wrap edge node) ::+ Set node
forall sh0 sh1. sh0 -> sh1 -> sh0 ::+ sh1
::+ Set node
nodes) (node -> Either (Wrap edge node) node
forall a b. b -> Either a b
Right node
src)
   in (()
-> ArrayMatrix
     Packed Symmetric Filled Filled Shape Small Small () () a
forall width a.
(C width, Floating a) =>
width
-> ArrayMatrix
     Packed Symmetric Filled Filled Shape Small Small width width a
symmetricZero (), Order
-> Vector (Set (Wrap edge node) ::+ Set node) a
-> General () (Set (Wrap edge node) ::+ Set node) a
forall width a. Order -> Vector width a -> General () width a
Matrix.singleRow Order
order Vector (Set (Wrap edge node) ::+ Set node) a
unit)
      #%%%#
      (Order
-> Vector (Set (Wrap edge node)) a
-> Symmetric (Set (Wrap edge node)) a
forall sh a.
(C sh, Floating a) =>
Order -> Vector sh a -> Symmetric sh a
Symmetric.diagonal Order
order (Vector (Set (Wrap edge node)) a
 -> Symmetric (Set (Wrap edge node)) a)
-> Vector (Set (Wrap edge node)) a
-> Symmetric (Set (Wrap edge node)) a
forall a b. (a -> b) -> a -> b
$ Map (Wrap edge node) a -> Vector (Set (Wrap edge node)) a
forall k a. (Ord k, Storable a) => Map k a -> Array (Set k) a
Array.fromMap Map (Wrap edge node) a
edges,
            Set node
-> Set (Wrap edge node)
-> General (Set (Wrap edge node)) (Set node) a
forall (edge :: * -> *) node a.
(Edge edge, Ord node, Floating a) =>
Set node
-> Set (Wrap edge node)
-> General (Set (Wrap edge node)) (Set node) a
voltageMatrix Set node
nodes (Set (Wrap edge node)
 -> General (Set (Wrap edge node)) (Set node) a)
-> Set (Wrap edge node)
-> General (Set (Wrap edge node)) (Set node) a
forall a b. (a -> b) -> a -> b
$ Map (Wrap edge node) a -> Set (Wrap edge node)
forall k a. Map k a -> Set k
Map.keysSet Map (Wrap edge node) a
edges)
      #%%%#
      Set node -> SymmetricP Packed (Set node) a
forall width a.
(C width, Floating a) =>
width
-> ArrayMatrix
     Packed Symmetric Filled Filled Shape Small Small width width a
symmetricZero Set node
nodes

_resistance ::
   (Graph.Edge edge, Ord node, Class.Floating a) =>
   Graph edge node a nodeLabel -> node -> node -> a
_resistance :: Graph edge node a nodeLabel -> node -> node -> a
_resistance Graph edge node a nodeLabel
gr node
src node
dst =
   let m :: Symmetric (() ::+ (Set (Wrap edge node) ::+ Set node)) a
m = Graph edge node a nodeLabel
-> node -> Symmetric (() ::+ (Set (Wrap edge node) ::+ Set node)) a
forall (edge :: * -> *) node a nodeLabel.
(Edge edge, Ord node, Floating a) =>
Graph edge node a nodeLabel
-> node -> Symmetric (() ::+ (Set (Wrap edge node) ::+ Set node)) a
fullMatrix Graph edge node a nodeLabel
gr node
src
       ix :: Either a (Either a node)
ix = Either a node -> Either a (Either a node)
forall a b. b -> Either a b
Right (node -> Either a node
forall a b. b -> Either a b
Right node
dst)
   in - (Symmetric (() ::+ (Set (Wrap edge node) ::+ Set node)) a
m Symmetric (() ::+ (Set (Wrap edge node) ::+ Set node)) a
-> Vector (() ::+ (Set (Wrap edge node) ::+ Set node)) a
-> Vector (() ::+ (Set (Wrap edge node) ::+ Set node)) a
forall typ xl xu lower upper meas height width a.
(Solve typ, ToQuadratic typ, SolveExtra typ xl, SolveExtra typ xu,
 BoxExtra typ xl, BoxExtra typ xu, Strip lower, Strip upper,
 Measure meas, C height, C width, Eq height, Floating a) =>
QuadraticMeas typ xl xu lower upper meas height width a
-> Vector height a -> Vector width a
#\| (() ::+ (Set (Wrap edge node) ::+ Set node))
-> Index (() ::+ (Set (Wrap edge node) ::+ Set node))
-> Vector (() ::+ (Set (Wrap edge node) ::+ Set node)) a
forall sh a.
(Indexed sh, Floating a) =>
sh -> Index sh -> Vector sh a
Vector.unit (Symmetric (() ::+ (Set (Wrap edge node) ::+ Set node)) a
-> () ::+ (Set (Wrap edge node) ::+ Set node)
forall pack sh a. SymmetricP pack sh a -> sh
Symmetric.size Symmetric (() ::+ (Set (Wrap edge node) ::+ Set node)) a
m) Index (() ::+ (Set (Wrap edge node) ::+ Set node))
forall a a. Either a (Either a node)
ix) Vector (() ::+ (Set (Wrap edge node) ::+ Set node)) a
-> Index (() ::+ (Set (Wrap edge node) ::+ Set node)) -> a
forall sh a.
(Indexed sh, Storable a) =>
Array sh a -> Index sh -> a
! Index (() ::+ (Set (Wrap edge node) ::+ Set node))
forall a a. Either a (Either a node)
ix

{-
The above computation with blockwise inversion.
-}
resistance ::
   (Graph.Edge edge, Ord node, Class.Floating a) =>
   Graph edge node a nodeLabel -> node -> node -> a
resistance :: Graph edge node a nodeLabel -> node -> node -> a
resistance Graph edge node a nodeLabel
gr node
src node
dst =
   let edges :: Map (IdentityT edge node) a
edges = (edge node -> IdentityT edge node)
-> Map (edge node) a -> Map (IdentityT edge node) a
forall k1 k2 a. (k1 -> k2) -> Map k1 a -> Map k2 a
Map.mapKeysMonotonic edge node -> IdentityT edge node
forall k (f :: k -> *) (a :: k). f a -> IdentityT f a
IdentityT (Map (edge node) a -> Map (IdentityT edge node) a)
-> Map (edge node) a -> Map (IdentityT edge node) a
forall a b. (a -> b) -> a -> b
$ Graph edge node a nodeLabel -> Map (edge node) a
forall (edge :: * -> *) node edgeLabel nodeLabel.
(Edge edge, Ord node) =>
Graph edge node edgeLabel nodeLabel -> Map (edge node) edgeLabel
Graph.edgeLabels Graph edge node a nodeLabel
gr
       nodes :: Set node
nodes = Graph edge node a nodeLabel -> Set node
forall (e :: * -> *) n el nl. Graph e n el nl -> Set n
Graph.nodeSet Graph edge node a nodeLabel
gr
       order :: Order
order = Order
MatrixShape.RowMajor
       schurComplement :: SymmetricP Packed (() ::+ Set node) a
schurComplement =
         (Omni Packed Symmetric Filled Filled Shape Small Small () ()
-> ArrayMatrix
     Packed Symmetric Filled Filled Shape Small Small () () a
forall property meas vert horiz height width a pack lower upper.
(Homogeneous property, Measure meas, C vert, C horiz, C height,
 C width, Floating a) =>
Omni pack property lower upper meas vert horiz height width
-> ArrayMatrix
     pack property lower upper meas vert horiz height width a
Matrix.zero (Omni Packed Symmetric Filled Filled Shape Small Small () ()
 -> ArrayMatrix
      Packed Symmetric Filled Filled Shape Small Small () () a)
-> Omni Packed Symmetric Filled Filled Shape Small Small () ()
-> ArrayMatrix
     Packed Symmetric Filled Filled Shape Small Small () () a
forall a b. (a -> b) -> a -> b
$ Order
-> ()
-> Omni Packed Symmetric Filled Filled Shape Small Small () ()
forall size. Order -> size -> Symmetric size
MatrixShape.symmetric Order
order (),
            Order -> Vector (Set node) a -> General () (Set node) a
forall width a. Order -> Vector width a -> General () width a
Matrix.singleRow Order
order (Vector (Set node) a -> General () (Set node) a)
-> Vector (Set node) a -> General () (Set node) a
forall a b. (a -> b) -> a -> b
$ Vector (Set node) a -> Vector (Set node) a
forall sh a. (C sh, Floating a) => Vector sh a -> Vector sh a
Vector.negate (Vector (Set node) a -> Vector (Set node) a)
-> Vector (Set node) a -> Vector (Set node) a
forall a b. (a -> b) -> a -> b
$ Set node -> Index (Set node) -> Vector (Set node) a
forall sh a.
(Indexed sh, Floating a) =>
sh -> Index sh -> Vector sh a
Vector.unit Set node
nodes node
Index (Set node)
src)
         #%%%#
         Vector (Set (IdentityT edge node)) a
-> General (Set (IdentityT edge node)) (Set node) a
-> SymmetricP Packed (Set node) a
forall pack height width a.
(Packing pack, C height, Eq height, C width, Floating a) =>
Vector height a
-> General height width a -> SymmetricP pack width a
Symmetric.congruenceDiagonal
            (Vector (Set (IdentityT edge node)) a
-> Vector (Set (IdentityT edge node)) a
forall sh a. (C sh, Floating a) => Vector sh a -> Vector sh a
Vector.recip (Vector (Set (IdentityT edge node)) a
 -> Vector (Set (IdentityT edge node)) a)
-> Vector (Set (IdentityT edge node)) a
-> Vector (Set (IdentityT edge node)) a
forall a b. (a -> b) -> a -> b
$ Map (IdentityT edge node) a -> Vector (Set (IdentityT edge node)) a
forall k a. (Ord k, Storable a) => Map k a -> Array (Set k) a
Array.fromMap Map (IdentityT edge node) a
edges)
            (Set node
-> Set (IdentityT edge node)
-> General (Set (IdentityT edge node)) (Set node) a
forall (edge :: * -> *) node a.
(Edge edge, Ord node, Floating a) =>
Set node
-> Set (Wrap edge node)
-> General (Set (Wrap edge node)) (Set node) a
voltageMatrix Set node
nodes (Set (IdentityT edge node)
 -> General (Set (IdentityT edge node)) (Set node) a)
-> Set (IdentityT edge node)
-> General (Set (IdentityT edge node)) (Set node) a
forall a b. (a -> b) -> a -> b
$ Map (IdentityT edge node) a -> Set (IdentityT edge node)
forall k a. Map k a -> Set k
Map.keysSet Map (IdentityT edge node) a
edges)
       ix :: Either a node
ix = node -> Either a node
forall a b. b -> Either a b
Right node
dst
   in (SymmetricP Packed (() ::+ Set node) a
schurComplement SymmetricP Packed (() ::+ Set node) a
-> Vector (() ::+ Set node) a -> Vector (() ::+ Set node) a
forall typ xl xu lower upper meas height width a.
(Solve typ, ToQuadratic typ, SolveExtra typ xl, SolveExtra typ xu,
 BoxExtra typ xl, BoxExtra typ xu, Strip lower, Strip upper,
 Measure meas, C height, C width, Eq height, Floating a) =>
QuadraticMeas typ xl xu lower upper meas height width a
-> Vector height a -> Vector width a
#\| (() ::+ Set node)
-> Index (() ::+ Set node) -> Vector (() ::+ Set node) a
forall sh a.
(Indexed sh, Floating a) =>
sh -> Index sh -> Vector sh a
Vector.unit (() () -> Set node -> () ::+ Set node
forall sh0 sh1. sh0 -> sh1 -> sh0 ::+ sh1
::+ Set node
nodes) Index (() ::+ Set node)
forall a. Either a node
ix) Vector (() ::+ Set node) a -> Index (() ::+ Set node) -> a
forall sh a.
(Indexed sh, Storable a) =>
Array sh a -> Index sh -> a
! Index (() ::+ Set node)
forall a. Either a node
ix