IndexedMatrix.hs

   1 module IndexedMatrix where
   2 import Control.Monad (liftM)
   3 import Data.Maybe (fromMaybe, isJust)
   4 import Data.Set (Set)
   5 import qualified Data.Set as Set
   6 import Numeric.Matrix (Matrix, MatrixElement)
   7 import qualified Numeric.Matrix as Matrix
   8 import Util.Function ((...))
   9 import qualified Util.Set as Set
  10
  11 data IndexedMatrix r c e = IndexedMatrix
  12   { rowIndices    :: Set r
  13   , columnIndices :: Set c
  14   , stripIndices  :: Matrix e
  15   } deriving (Show, Eq)
  16
  17 indexedMatrix :: MatrixElement e =>
  18   Set r -> Set c -> (r -> c -> e) -> IndexedMatrix r c e
  19 indexedMatrix rs cs f = IndexedMatrix
  20   { rowIndices    = rs
  21   , columnIndices = cs
  22   , stripIndices  = Matrix.matrix (Set.size rs, Set.size cs) $ \(i, j) ->
  23       f (Set.elemAt (i - 1) rs) $ Set.elemAt (j - 1) cs
  24   }
  25
  26 maybeInRange :: (Ord r, Ord c) =>
  27   (IndexedMatrix r c e -> r -> c -> a) ->
  28     IndexedMatrix r c e -> r -> c -> Maybe a
  29 maybeInRange f m r c =
  30   if (Set.member r (rowIndices m) && Set.member c (columnIndices m))
  31   then Just $ f m r c
  32   else Nothing
  33
  34 at :: (Ord r, Ord c, MatrixElement e) =>
  35   IndexedMatrix r c e -> r -> c -> Maybe e
  36 at = maybeInRange $ \m r c -> Matrix.at (stripIndices m)
  37   (Set.findIndex r (rowIndices m) + 1, Set.findIndex c (columnIndices m) + 1)
  38
  39 at' :: (Ord r, Ord c, MatrixElement e) => IndexedMatrix r c e -> r -> c -> e
  40 at' = fromMaybe 0 ... at
  41
  42 firstMinorMatrix :: (Ord r, Ord c, MatrixElement e) =>
  43   IndexedMatrix r c e -> r -> c -> Maybe (IndexedMatrix r c e)
  44 firstMinorMatrix = maybeInRange $ \m r c ->
  45   let
  46     rs = rowIndices m
  47     cs = columnIndices m
  48   in
  49   IndexedMatrix
  50     { rowIndices    = Set.delete r rs
  51     , columnIndices = Set.delete c cs
  52     , stripIndices  = flip Matrix.minorMatrix (stripIndices m)
  53         (Set.findIndex r rs + 1, Set.findIndex c cs + 1)
  54     }
  55
  56 inverse :: MatrixElement e => IndexedMatrix r c e -> Maybe (IndexedMatrix c r e)
  57 inverse m = flip liftM (Matrix.inv $ stripIndices m) $ \m' -> IndexedMatrix
  58   { rowIndices = columnIndices m
  59   , columnIndices = rowIndices m
  60   , stripIndices = m'
  61   }
  62
  63 invertible :: MatrixElement e => IndexedMatrix r c e -> Bool
  64 invertible = isJust . inverse
  65
  66 identity :: MatrixElement e => Set i -> IndexedMatrix i i e
  67 identity i = IndexedMatrix
  68   { rowIndices = i
  69   , columnIndices = i
  70   , stripIndices = Matrix.unit $ Set.size i
  71   }
  72
  73 times :: (Eq j, MatrixElement e) =>
  74   IndexedMatrix i j e -> IndexedMatrix j k e -> Maybe (IndexedMatrix i k e)
  75 times m n =
  76   if columnIndices m == rowIndices n
  77   then Just $ IndexedMatrix
  78     { rowIndices = rowIndices m
  79     , columnIndices = columnIndices n
  80     , stripIndices = stripIndices m * stripIndices n
  81     }
  82   else Nothing
  83
  84 numRows :: IndexedMatrix r c e -> Integer
  85 numRows = toInteger . Set.size . rowIndices
  86
  87 numColumns :: IndexedMatrix r c e -> Integer
  88 numColumns = toInteger . Set.size . columnIndices
  89