{-# LANGUAGE RebindableSyntax #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FlexibleInstances #-} {- | Copyright : (c) Henning Thielemann 2009, Mikael Johansson 2006 Maintainer : numericprelude@henning-thielemann.de Stability : provisional Portability : requires multi-parameter type classes Routines and abstractions for Matrices and basic linear algebra over fields or rings. We stick to simple Int indices. Although advanced indices would be nice e.g. for matrices with sub-matrices, this is not easily implemented since arrays do only support a lower and an upper bound but no additional parameters. ToDo: - Matrix inverse, determinant -} module MathObj.Matrix ( T, Dimension, format, transpose, rows, columns, index, fromRows, fromColumns, fromList, dimension, numRows, numColumns, zipWith, zero, one, diagonal, scale, random, randomR, ) where import qualified Algebra.Module as Module import qualified Algebra.Vector as Vector import qualified Algebra.Ring as Ring import qualified Algebra.Additive as Additive import Algebra.Module((*>), ) import Algebra.Ring((*), fromInteger, scalarProduct, ) import Algebra.Additive((+), (-), subtract, ) import qualified System.Random as Rnd import Data.Array (Array, array, listArray, accumArray, elems, bounds, (!), ixmap, range, ) import qualified Data.List as List import Control.Monad (liftM2, ) import Control.Exception (assert, ) import Data.Function.HT (powerAssociative, ) import Data.Tuple.HT (swap, mapFst, ) import Data.List.HT (outerProduct, ) import Text.Show.HT (concatS, ) import NumericPrelude.Numeric (Int, ) import NumericPrelude.Base hiding (zipWith, ) {- | A matrix is a twodimensional array, indexed by integers. -} data T a = Cons (Array (Dimension, Dimension) a) deriving (Eq,Ord,Read) type Dimension = Int {- | Transposition of matrices is just transposition in the sense of Data.List. -} transpose :: T a -> T a transpose (Cons m) = let (lower,upper) = bounds m in Cons (ixmap (swap lower, swap upper) swap m) rows :: T a -> [[a]] rows mM@(Cons m) = let ((lr,lc), (ur,uc)) = bounds m in outerProduct (index mM) (range (lr,ur)) (range (lc,uc)) columns :: T a -> [[a]] columns mM@(Cons m) = let ((lr,lc), (ur,uc)) = bounds m in outerProduct (flip (index mM)) (range (lc,uc)) (range (lr,ur)) index :: T a -> Dimension -> Dimension -> a index (Cons m) i j = m ! (i,j) fromRows :: Dimension -> Dimension -> [[a]] -> T a fromRows m n = Cons . array (indexBounds m n) . concat . List.zipWith (\r -> map (\(c,x) -> ((r,c),x))) allIndices . map (zip allIndices) fromColumns :: Dimension -> Dimension -> [[a]] -> T a fromColumns m n = Cons . array (indexBounds m n) . concat . List.zipWith (\r -> map (\(c,x) -> ((c,r),x))) allIndices . map (zip allIndices) fromList :: Dimension -> Dimension -> [a] -> T a fromList m n xs = Cons (listArray (indexBounds m n) xs) appPrec :: Int appPrec = 10 instance (Show a) => Show (T a) where showsPrec p m = showParen (p >= appPrec) (showString "Matrix.fromRows " . showsPrec appPrec (rows m)) format :: (Show a) => T a -> String format m = formatS m "" formatS :: (Show a) => T a -> ShowS formatS = concatS . map (\r -> showString "(" . concatS r . showString ")\n") . map (List.intersperse (' ':) . map (showsPrec 11)) . rows dimension :: T a -> (Dimension,Dimension) dimension (Cons m) = uncurry subtract (bounds m) + (1,1) numRows :: T a -> Dimension numRows = fst . dimension numColumns :: T a -> Dimension numColumns = snd . dimension -- These implementations may benefit from a better exception than -- just assertions to validate dimensionalities instance (Additive.C a) => Additive.C (T a) where (+) = zipWith (+) (-) = zipWith (-) zero = zero 1 1 zipWith :: (a -> b -> c) -> T a -> T b -> T c zipWith op mM@(Cons m) nM@(Cons n) = let d = dimension mM em = elems m en = elems n in assert (d == dimension nM) $ uncurry fromList d (List.zipWith op em en) zero :: (Additive.C a) => Dimension -> Dimension -> T a zero m n = fromList m n $ List.repeat Additive.zero -- List.replicate (fromInteger (m*n)) zero one :: (Ring.C a) => Dimension -> T a one n = Cons $ accumArray (flip const) Additive.zero (indexBounds n n) (map (\i -> ((i,i), Ring.one)) (indexRange n)) diagonal :: (Additive.C a) => [a] -> T a diagonal xs = let n = List.length xs in Cons $ accumArray (flip const) Additive.zero (indexBounds n n) (zip (map (\i -> (i,i)) allIndices) xs) scale :: (Ring.C a) => a -> T a -> T a scale s = Vector.functorScale s instance (Ring.C a) => Ring.C (T a) where mM * nM = assert (numColumns mM == numRows nM) $ fromList (numRows mM) (numColumns nM) $ liftM2 scalarProduct (rows mM) (columns nM) fromInteger n = fromList 1 1 [fromInteger n] mM ^ n = assert (numColumns mM == numRows mM) $ assert (n >= Additive.zero) $ powerAssociative (*) (one (numColumns mM)) mM n instance Functor T where fmap f (Cons m) = Cons (fmap f m) instance Vector.C T where zero = Additive.zero (<+>) = (+) (*>) = scale instance Module.C a b => Module.C a (T b) where x *> m = fmap (x*>) m random :: (Rnd.RandomGen g, Rnd.Random a) => Dimension -> Dimension -> g -> (T a, g) random = randomAux Rnd.random randomR :: (Rnd.RandomGen g, Rnd.Random a) => Dimension -> Dimension -> (a,a) -> g -> (T a, g) randomR m n rng = randomAux (Rnd.randomR rng) m n {- could be made nicer with the State monad, but I like to keep dependencies minimal -} randomAux :: (Rnd.RandomGen g, Rnd.Random a) => (g -> (a,g)) -> Dimension -> Dimension -> g -> (T a, g) randomAux rnd m n g0 = mapFst (fromList m n) $ swap $ List.mapAccumL (\g _i -> swap $ rnd g) g0 (indexRange (m*n)) {- What more do we need from our matrix type? We have addition, subtraction and multiplication, and thus composition of generic free-module-maps. We're going to want to solve linear equations with or without fields underneath, so we're going to want an implementation of the Gaussian algorithm as well as most probably Smith normal form. Determinants are cool, and these are to be calculated either with the Gaussian algorithm or some other goodish method. -} {- {- | We'll want generic linear equation solving, returning one solution, any solution really, or nothing. Basically, this is asking for the preimage of a given vector over the given map, so a_11 x_1 + .. + a_1n x_n = y_1 ... a_m1 x_1 + .. + a_mn a_n = y_m has really x_1,...,x_n as a preimage of the vector y_1,..,y_m under the map (a_ij), since obviously y_1,..,y_m = (a_ij) x_1,..,x_n So, generic linear equation solving boils down to the function -} preimage :: (Ring.C a) => T a -> T a -> Maybe (T a) preimage a y = assert (numRows a == numRows y && -- they match numColumns y == 1) -- and y is a column vector Nothing -} {- Cf. /usr/lib/hugs/demos/Matrix.hs -} -- these functions control whether we use 0 or 1 based indices indexRange :: Dimension -> [Dimension] indexRange n = [0..(n-1)] indexBounds :: Dimension -> Dimension -> ((Dimension,Dimension), (Dimension,Dimension)) indexBounds m n = ((0,0), (m-1,n-1)) allIndices :: [Dimension] allIndices = [0..]