{- |
<https://en.wikipedia.org/wiki/Row_echelon_form>
Duplicate of @htam@ package.
-}
module RowEchelon (
Matrix(Matrix), matrixRows, matrixWidth, matrixHeight,
shortBesidesTall, (|||),
matrixValid, matrixFromRows,
identity, matrixProduct,
nullspace, reducedRowEchelon, rowReduction, layoutEchelonBlocks,
scatter,
Zipper0(Zipper0), Zipper1(Zipper1), altOuter,
) where
import Orthogonals (matrix_matrix)
import qualified Data.Foldable as Fold
import qualified Data.NonEmpty as NonEmpty
import qualified Data.List.HT as ListHT
import qualified Data.List as List
import Data.NonEmpty ((!:))
import Data.Maybe (fromMaybe)
{- $setup
>>> import RowEchelon as Matrix
>>> import qualified Test.QuickCheck as QC
>>> import Test.QuickCheck ((===))
>>> import Control.Monad (replicateM)
>>> import qualified Data.Foldable as Fold
>>> import qualified Data.NonEmpty.Class as NonEmptyC
>>> import qualified Data.NonEmpty as NonEmpty
>>> import qualified Data.List.Match as Match
>>> import qualified Data.List.HT as ListHT
>>> import Data.NonEmpty ((!:))
>>> import Data.Ratio ((%))
>>>
>>> genElementUniform, genElementNearZero :: QC.Gen Integer
>>> genElementUniform = QC.choose (-10,10)
>>> genElementNearZero = fmap (flip rem 11) $ QC.arbitrary
>>>
>>> genMatrixForSize :: Int -> Int -> QC.Gen (Matrix Integer)
>>> genMatrixForSize m n = do
>>> fmap (Matrix n) $ replicateM m $ replicateM n genElementNearZero
>>>
>>> genMatrix :: QC.Gen (Matrix Integer)
>>> genMatrix = do
>>> m <- QC.choose (0,10)
>>> n <- QC.choose (0,10)
>>> genMatrixForSize m n
>>>
>>> shrinkMatrix :: (Eq a) => Matrix a -> [Matrix a]
>>> shrinkMatrix matrix@(Matrix width rows) =
>>> filter (matrix/=) $
>>> map (Matrix width . snd) (ListHT.removeEach rows)
>>> ++
>>> [Matrix (width-1) $ map (drop 1) rows]
>>>
>>> forMatrix :: (QC.Testable test) => (Matrix Rational -> test) -> QC.Property
>>> forMatrix prop =
>>> QC.forAllShrink genMatrix shrinkMatrix (prop . rationalMatrix)
>>>
>>> rationalMatrix :: Matrix Integer -> Matrix Rational
>>> rationalMatrix = fmap (%1)
-}
data Zipper0 a = Zipper0 [a] [a]
deriving (Eq, Show)
data Zipper1 a = Zipper1 [a] a [a]
deriving (Eq, Show)
instance Functor Zipper0 where
fmap f (Zipper0 xs ys) = Zipper0 (fmap f xs) (fmap f ys)
instance Functor Zipper1 where
fmap f (Zipper1 xs y zs) = Zipper1 (fmap f xs) (f y) (fmap f zs)
zipper0FromList :: [a] -> Zipper0 a
zipper0FromList = Zipper0 []
type List0 = []
type List1 = NonEmpty.T List0
-- data Matrix a = Matrix {matrixWidth :: Int, matrixRows :: [[a]]}
data Matrix a = Matrix Int [[a]]
deriving (Eq, Show)
instance Functor Matrix where
fmap f (Matrix width rows) = Matrix width $ map (map f) rows
instance Foldable Matrix where
foldMap f (Matrix _width rows) = foldMap (foldMap f) rows
matrixValid :: Matrix a -> Bool
matrixValid (Matrix width rows) = all ((width==) . length) rows
matrixWidth :: Matrix a -> Int
matrixWidth (Matrix width _) = width
matrixHeight :: Matrix a -> Int
matrixHeight = length . matrixRows
matrixRows :: Matrix a -> [[a]]
matrixRows (Matrix _ rows) = rows
matrixFromRows :: [[a]] -> Matrix a
matrixFromRows rows =
Matrix (ListHT.switchL 0 (const . length) rows) rows
infixr 3 |||
-- | requires that both matrices have the same height
(|||) :: Matrix a -> Matrix a -> Matrix a
Matrix widthA a ||| Matrix widthB b =
Matrix (widthA+widthB) (zipWith (++) a b)
{- |
The expression @shortBesidesTall a b@ means:
> /A B_upper\
> | |
> \0 B_lower/
Matrix @a@ must be at most as tall as @b@.
-}
shortBesidesTall :: (Num a) => Matrix a -> Matrix a -> Matrix a
shortBesidesTall a b =
Matrix (matrixWidth a + matrixWidth b) $
zipWith (++)
(matrixRows a ++ repeat (replicate (matrixWidth a) 0))
(matrixRows b)
matrixProduct :: (Num a) => Matrix a -> Matrix a -> Matrix a
matrixProduct (Matrix _ x) (Matrix width y) =
Matrix width $ matrix_matrix x $ List.transpose y
{- |
>>> nullspace (Matrix 2 []) :: Matrix Rational
Matrix 2 [[1 % 1,0 % 1],[0 % 1,1 % 1]]
prop> forMatrix $ matrixValid . nullspace
prop> :{
forMatrix $ \matrix ->
matrixWidth matrix == matrixHeight (nullspace matrix)
:}
prop> :{
forMatrix $ \matrix ->
matrixWidth (nullspace matrix) <= matrixWidth matrix
:}
max 0 (width matrix - height matrix) <= width nullspace
prop> :{
forMatrix $ \matrix ->
matrixWidth matrix <= matrixWidth (nullspace matrix) + matrixHeight matrix
:}
prop> :{
forMatrix $ \matrix ->
Fold.all (0==) $ matrixProduct matrix (nullspace matrix)
:}
-}
nullspace :: (RealFrac a) => Matrix a -> Matrix a
nullspace matrix =
let echelon = reducedRowEchelon matrix in
let flat = layoutEchelonBlocks echelon in
let nullDim = matrixWidth flat in
Matrix nullDim $
scatter
(fmap matrixWidth $ altOuter echelon)
(matrixRows $ fmap negate flat) $
matrixRows $ identity nullDim
identity :: (Num a) => Int -> Matrix a
identity n = Matrix n $ take n $ map (take n) $ iterate (0:) $ 1 : repeat 0
-- cf. event-list
data AlternatingList a b = AlternatingList a [(b,a)]
deriving (Eq)
instance (Show a, Show b) => Show (AlternatingList a b) where
showsPrec p xs =
showParen (p>=5) $
flip
(altFoldr
(\a -> showsPrec 5 a . showString " ./ ")
(\b -> showsPrec 5 b . showString " /. "))
xs
.
showString "[]"
altSingleton :: a -> AlternatingList a b
altSingleton a = AlternatingList a []
infixr 5 /. , ./
(./) = altConsOuter
(/.) = altConsInner
altConsOuter, (./) :: a -> List0 (b,a) -> AlternatingList a b
altConsOuter = AlternatingList
altConsInner, (/.) :: b -> AlternatingList a b -> List0 (b,a)
altConsInner b ~(AlternatingList a bas) = (b,a) : bas
altOuter :: AlternatingList a b -> List1 a
altOuter (AlternatingList a bas) = a !: map snd bas
altFoldr :: (a -> c -> d) -> (b -> d -> c) -> c -> AlternatingList a b -> d
altFoldr f g c (AlternatingList a0 bas) =
f a0 $ foldr (\(b,a) -> g b . f a) c bas
{- |
>>> scatter (3 !: 1 : 4 : 1 : []) "012" ['a'..'k']
"abc0d1efgh2ijk"
-}
scatter :: List1 Int -> [a] -> [a] -> [a]
scatter (NonEmpty.Cons k_ ks_) =
let go _k [] [] dst = dst
go k0 (k1:ks) (a:as) dst =
case splitAt k0 dst of
(prefix,suffix) -> prefix ++ a : go k1 ks as suffix
go _ _ _ _ = error "scatter: inconsistent lengths"
in go k_ ks_
{- |
> mapM_ print $ layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[1,3,0,5],[0,0,1,7::Rational]]
>>> layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[0::Rational]]
Matrix 1 []
>>> layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[1,3,0,5::Rational]]
Matrix 3 [[3 % 1,0 % 1,5 % 1]]
>>> layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[0,1,3,5::Rational]]
Matrix 3 [[0 % 1,3 % 1,5 % 1]]
>>> layoutEchelonBlocks $ reducedRowEchelon $ matrixFromRows [[1,3,0,5],[0,0,1,7::Rational]]
Matrix 2 [[3 % 1,5 % 1],[0 % 1,7 % 1]]
prop> :{
forMatrix $ \matrix ->
(layoutEchelonBlocks $ reducedRowEchelon $
identity (matrixHeight matrix) ||| matrix)
===
matrix
:}
Construct a matrix that is already in Echelon form
and check that it is preserved by Echelon decomposition.
prop> :{
QC.forAll (fmap (flip mod 10) . NonEmpty.mapTail (take 9) <$> QC.arbitrary) $
\blockWidths ->
QC.forAll (traverse (uncurry genMatrixForSize) $
NonEmptyC.zip (0!:[1..]) blockWidths) $
\blocksInt ->
let blocks :: NonEmpty.T [] (Matrix Rational)
blocks = fmap (fmap (%1)) blocksInt
unitBesidesBlock block =
Matrix (matrixWidth block + 1) $
zipWith (:)
(drop 1 $ Match.replicate (matrixRows block) 0 ++ [1])
(matrixRows block)
echelonWithInterleavedEye =
Fold.foldr1 shortBesidesTall $
NonEmpty.mapTail (map unitBesidesBlock) blocks
in
layoutEchelonBlocks (reducedRowEchelon echelonWithInterleavedEye)
===
Fold.foldr1 shortBesidesTall blocks
:}
-}
layoutEchelonBlocks :: (Num a) => AlternatingList (Matrix a) a -> Matrix a
layoutEchelonBlocks = Fold.foldr1 shortBesidesTall . altOuter
{-
/I B\ /0\
| | * P * x = | |
\0 0/ \0/
B is a rectangular upper block triangular matrix with non-uniform block size.
The solution is
/-B\
P*X = | |
\ I/
such that
/-B\
(I B) * | | = I*(-B) + B*I = 0
\ I/
This is the nullspace:
/-B\
P^T * | |
\ I/
The null-space vectors have the form
/* * *\
|1 0 0|
|* * *|
|* * *|
|* * *|
|0 1 0|
|* * *|
|0 0 1|
|* * *|
\* * */
-}
{- |
>>> reducedRowEchelon $ matrixFromRows [[2,1,3::Rational]]
Matrix 0 [] ./ 2 % 1 /. Matrix 2 [[1 % 2,3 % 2]] ./ []
>>> reducedRowEchelon $ matrixFromRows [[2],[1],[3::Rational]]
Matrix 0 [] ./ 3 % 1 /. Matrix 0 [[]] ./ []
>>> reducedRowEchelon $ matrixFromRows [[1,0,2],[0,1,3::Rational]]
Matrix 0 [] ./ 1 % 1 /. Matrix 0 [[]] ./ 1 % 1 /. Matrix 1 [[2 % 1],[3 % 1]] ./ []
>>> reducedRowEchelon $ matrixFromRows [[1,3,0,5],[0,0,1,7::Rational]]
Matrix 0 [] ./ 1 % 1 /. Matrix 1 [[3 % 1]] ./ 1 % 1 /. Matrix 1 [[5 % 1],[7 % 1]] ./ []
prop> :{
forMatrix $ \matrix ->
Fold.all matrixValid $ altOuter $ reducedRowEchelon matrix
:}
-}
reducedRowEchelon :: (RealFrac a) => Matrix a -> AlternatingList (Matrix a) a
reducedRowEchelon matrix@(Matrix _ []) = altSingleton matrix
reducedRowEchelon (Matrix _ rows) =
let go matrix@(Zipper0 upper _lower) =
case rowReduction matrix of
(n, Nothing) -> altSingleton (Matrix n $ reverse upper)
(n, Just (hd, block, remaining)) ->
Matrix n block ./ hd /. go remaining
in go $ zipper0FromList rows
{- |
>>> rowReduction $ Zipper0 [[1,2,3,4,5],[6,7,8,9,10::Rational]] [[0,2,0,0,0],[3,0,0,0,0]]
(0,Just (3 % 1,[[],[]],Zipper0 [[0 % 1,0 % 1,0 % 1,0 % 1],[2 % 1,3 % 1,4 % 1,5 % 1],[7 % 1,8 % 1,9 % 1,10 % 1]] [[2 % 1,0 % 1,0 % 1,0 % 1]]))
-}
rowReduction ::
(RealFrac a) => Zipper0 [a] -> (Int, Maybe (a, [[a]], Zipper0 [a]))
rowReduction (Zipper0 upper []) =
(ListHT.switchL 0 (\row _rows -> length row) upper, Nothing)
rowReduction (Zipper0 upper (focus:lower)) =
case shiftUntilNonZeroColumn 0 $ Zipper1 upper focus lower of
(n, Nothing) -> (n, Nothing)
(n, Just (Zipper1 upperNE (NonEmpty.Cons maxRowHead maxRow) lowerNE)) ->
(n, Just $
let maxRowNormalized = map (/maxRowHead) maxRow in
(maxRowHead,
map (take n) $ reverse upper,
Zipper0
(maxRowNormalized :
map (cancelRow maxRowNormalized) upperNE)
(map (cancelRow maxRowNormalized) lowerNE)))
cancelRow :: (Num a) => List0 a -> List1 a -> List0 a
cancelRow xs (NonEmpty.Cons y0 ys) = zipWith (-) ys $ map (y0*) xs
shiftUntilNonZeroColumn ::
(Real a) => Int -> Zipper1 [a] -> (Int, Maybe (Zipper1 (List1 a)))
shiftUntilNonZeroColumn n matrix =
case checkHeads matrix of
Nothing -> (n, Nothing)
Just matrixNE@(Zipper1 upperNE foucsNE lowerNE) ->
case maxHead (foucsNE!:lowerNE) of
(maxRow, lowerWithoutMax) ->
if NonEmpty.head maxRow == 0
then (shiftUntilNonZeroColumn $! (n+1)) $
fmap NonEmpty.tail matrixNE
else (n, Just (Zipper1 upperNE maxRow lowerWithoutMax))
maxHead :: (Real a) => List1 (List1 a) -> (List1 a, [List1 a])
maxHead = NonEmpty.maximumKey (abs . NonEmpty.head . fst) . NonEmpty.removeEach
checkHeads :: Zipper1 [a] -> Maybe (Zipper1 (List1 a))
checkHeads (Zipper1 upper focus lower) =
case NonEmpty.fetch focus of
Nothing ->
if all null (upper++lower) then Nothing else error "row lengths differ"
Just xs -> Just $
Zipper1
(fromMaybe (error "upper row lengths differ") $ allNonEmpty upper)
xs
(fromMaybe (error "lower row lengths differ") $ allNonEmpty lower)
allNonEmpty :: [List0 a] -> Maybe [List1 a]
allNonEmpty xs =
case ListHT.partitionMaybe NonEmpty.fetch xs of
(ne,[]) -> Just ne
_ -> Nothing