From ae104ebd5891c84f9c8b4a40501fefdeeb1280c4 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 2 May 2014 14:12:22 +0200 Subject: (<>) and (×) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit (<>) back for compat, added (×), LSDiv, +documentation, classes to dev, Monoid using optimiseMult, --- lib/Numeric/Container.hs | 99 ++++++++++++++++++++++++++++----------- lib/Numeric/HMatrix.hs | 43 ++++++++++++----- lib/Numeric/HMatrix/Data.hs | 6 --- lib/Numeric/HMatrix/Devel.hs | 13 +++-- lib/Numeric/LinearAlgebra/Util.hs | 8 ++-- lib/Numeric/Matrix.hs | 27 +++++++++++ 6 files changed, 143 insertions(+), 53 deletions(-) (limited to 'lib') diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs index d1ce588..a71fdfe 100644 --- a/lib/Numeric/Container.hs +++ b/lib/Numeric/Container.hs @@ -37,7 +37,8 @@ module Numeric.Container ( Container(..), -- * Matrix product Product(..), - Contraction(..), + Mul(..), + Contraction(..), mmul, optimiseMult, mXm,mXv,vXm,LSDiv(..), cdot, (·), dot, (<.>), outer, kronecker, @@ -102,12 +103,19 @@ cdot u v = udot (conj u) v -------------------------------------------------------- -class Contraction a b c | a b -> c, a c -> b, b c -> a +class Contraction a b c | a b -> c, c -> a b where - infixl 7 <> + infixr 7 × {- | Matrix-matrix product, matrix-vector product, and unconjugated dot product ->>> let a = (3><4) [1..] :: Matrix Double +(unicode 0x00d7, multiplication sign) + +Examples: + +>>> let a = (3><4) [1..] :: Matrix Double +>>> let v = fromList [1,0,2,-1] :: Vector Double +>>> let u = fromList [1,2,3] :: Vector Double + >>> a (3><4) [ 1.0, 2.0, 3.0, 4.0 @@ -116,7 +124,7 @@ class Contraction a b c | a b -> c, a c -> b, b c -> a matrix × matrix: ->>> disp 2 (a <> trans a) +>>> disp 2 (a × trans a) 3x3 30 70 110 70 174 278 @@ -124,52 +132,79 @@ matrix × matrix: matrix × vector: ->>> a <> fromList [1,0,2,-1::Double] +>>> a × v fromList [3.0,11.0,19.0] -vector × matrix: - ->>> fromList [1,2,3::Double] <> a -fromList [38.0,44.0,50.0,56.0] - unconjugated dot product: ->>> fromList [1,i] <> fromList[2*i+1,3] +>>> fromList [1,i] × fromList[2*i+1,3] 1.0 :+ 5.0 --} - (<>) :: a -> b -> c +(×) is right associative, so we can write: -instance Product t => Contraction (Vector t) (Vector t) t where - (<>) = udot +>>> u × a × v +82.0 :: Double -instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where - (<>) = mXv +-} + (×) :: a -> b -> c -instance Product t => Contraction (Vector t) (Matrix t) (Vector t) where - (<>) = vXm +instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where + (×) = mXv instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where + (×) = mXm + +instance Contraction (Vector Double) (Vector Double) Double where + (×) = udot + +instance Contraction (Vector Float) (Vector Float) Float where + (×) = udot + +instance Contraction (Vector (Complex Double)) (Vector (Complex Double)) (Complex Double) where + (×) = udot + +instance Contraction (Vector (Complex Float)) (Vector (Complex Float)) (Complex Float) where + (×) = udot + + +-- | alternative function for the matrix product (×) +mmul :: Contraction a b c => a -> b -> c +mmul = (×) + +-------------------------------------------------------------------------------- + +class Mul a b c | a b -> c where + infixl 7 <> + -- | Matrix-matrix, matrix-vector, and vector-matrix products. + (<>) :: Product t => a t -> b t -> c t + +instance Mul Matrix Matrix Matrix where (<>) = mXm --------------------------------------------------------- +instance Mul Matrix Vector Vector where + (<>) m v = flatten $ m <> asColumn v -class LSDiv b c | b -> c, c->b where +instance Mul Vector Matrix Vector where + (<>) v m = flatten $ asRow v <> m + +-------------------------------------------------------------------------------- + +class LSDiv c where infixl 7 <\> -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) - (<\>) :: Field t => Matrix t -> b t -> c t + (<\>) :: Field t => Matrix t -> c t -> c t -instance LSDiv Vector Vector where +instance LSDiv Vector where m <\> v = flatten (linearSolveSVD m (reshape 1 v)) -instance LSDiv Matrix Matrix where +instance LSDiv Matrix where (<\>) = linearSolveSVD -------------------------------------------------------- {- | Dot product : @u · v = 'cdot' u v@ - (unicode 0x00b7, Alt-Gr .) + (unicode 0x00b7, middle dot, Alt-Gr .) >>> fromList [1,i] · fromList[2*i+1,3] 1.0 :+ (-1.0) @@ -233,7 +268,15 @@ instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e -------------------------------------------------------------------------------- --- | Compute mean vector and covariance matrix of the rows of a matrix. +{- | Compute mean vector and covariance matrix of the rows of a matrix. + +>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3]) +(fromList [4.010341078059521,5.0197204699640405], +(2><2) + [ 1.9862461923890056, -1.0127225830525157e-2 + , -1.0127225830525157e-2, 3.0373954915729318 ]) + +-} meanCov :: Matrix Double -> (Vector Double, Matrix Double) meanCov x = (med,cov) where r = rows x @@ -249,7 +292,7 @@ meanCov x = (med,cov) where dot :: Product e => Vector e -> Vector e -> e dot = udot -{-# DEPRECATED (<.>) "use udot or (<>)" #-} +{-# DEPRECATED (<.>) "use udot or (×)" #-} infixl 7 <.> (<.>) :: Product e => Vector e -> Vector e -> e (<.>) = udot diff --git a/lib/Numeric/HMatrix.hs b/lib/Numeric/HMatrix.hs index 8e0b4a2..a2f09df 100644 --- a/lib/Numeric/HMatrix.hs +++ b/lib/Numeric/HMatrix.hs @@ -16,25 +16,45 @@ module Numeric.HMatrix ( -- * Basic types and data processing module Numeric.HMatrix.Data, - -- | The standard numeric classes are defined elementwise. + -- | The standard numeric classes are defined elementwise: -- -- >>> fromList [1,2,3] * fromList [3,0,-2 :: Double] -- fromList [3.0,0.0,-6.0] -- - -- In arithmetic operations single-element vectors and matrices automatically - -- expand to match the dimensions of the other operand. + -- >>> (3><3) [1..9] * ident 3 :: Matrix Double + -- (3><3) + -- [ 1.0, 0.0, 0.0 + -- , 0.0, 5.0, 0.0 + -- , 0.0, 0.0, 9.0 ] + -- + -- In arithmetic operations single-element vectors and matrices + -- (created from numeric literals or using 'scalar') automatically + -- expand to match the dimensions of the other operand: -- - -- >>> 2 * ident 3 - -- 2 * ident 3 :: Matrix Double + -- >>> 5 + 2*ident 3 :: Matrix Double -- (3><3) - -- [ 2.0, 0.0, 0.0 - -- , 0.0, 2.0, 0.0 - -- , 0.0, 0.0, 2.0 ] + -- [ 7.0, 5.0, 5.0 + -- , 5.0, 7.0, 5.0 + -- , 5.0, 5.0, 7.0 ] -- -- * Products - (<>), (·), outer, kronecker, cross, - optimiseMult, scale, + (×), + + -- | The matrix product is also implemented in the "Data.Monoid" instance for Matrix, where + -- single-element matrices (created from numeric literals or using 'scalar') + -- are used for scaling. + -- + -- >>> let m = (2><3)[1..] :: Matrix Double + -- >>> m <> 2 <> diagl[0.5,1,0] + -- (2><3) + -- [ 1.0, 4.0, 0.0 + -- , 4.0, 10.0, 0.0 ] + -- + -- mconcat uses 'optimiseMult' to get the optimal association order. + + (·), outer, kronecker, cross, + scale, sumElements, prodElements, absSum, -- * Linear Systems @@ -103,7 +123,7 @@ module Numeric.HMatrix ( rand, randn, RandDist(..), randomVector, gaussianSample, uniformSample, -- * Misc - meanCov, peps, relativeError, haussholder + meanCov, peps, relativeError, haussholder, optimiseMult, udot, cdot, mmul ) where import Numeric.HMatrix.Data @@ -114,4 +134,3 @@ import Numeric.Container import Numeric.LinearAlgebra.Algorithms import Numeric.LinearAlgebra.Util - diff --git a/lib/Numeric/HMatrix/Data.hs b/lib/Numeric/HMatrix/Data.hs index 49dad10..288b0af 100644 --- a/lib/Numeric/HMatrix/Data.hs +++ b/lib/Numeric/HMatrix/Data.hs @@ -51,12 +51,6 @@ module Numeric.HMatrix.Data( -- * Conversion Convert(..), - Complexable(), - RealElement(), - - RealOf, ComplexOf, SingleOf, DoubleOf, - - IndexOf, -- * Misc arctan2, diff --git a/lib/Numeric/HMatrix/Devel.hs b/lib/Numeric/HMatrix/Devel.hs index 37bf826..7363477 100644 --- a/lib/Numeric/HMatrix/Devel.hs +++ b/lib/Numeric/HMatrix/Devel.hs @@ -50,15 +50,20 @@ module Numeric.HMatrix.Devel( mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_, liftMatrix, liftMatrix2, liftMatrix2Auto, - -- * Misc - Element, Container, Product, Contraction, LSDiv, Field + -- * Auxiliary classes + Element, Container, Product, Contraction, LSDiv, + Complexable(), RealElement(), + RealOf, ComplexOf, SingleOf, DoubleOf, + IndexOf, + Field, ) where import Data.Packed.Foreign import Data.Packed.Development import Data.Packed.ST -import Numeric.Container(Container,Contraction,LSDiv,Product) +import Numeric.Container(Container,Contraction,LSDiv,Product, + Complexable(),RealElement(), + RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf) import Data.Packed import Numeric.LinearAlgebra.Algorithms(Field) - diff --git a/lib/Numeric/LinearAlgebra/Util.hs b/lib/Numeric/LinearAlgebra/Util.hs index 21b6188..7164827 100644 --- a/lib/Numeric/LinearAlgebra/Util.hs +++ b/lib/Numeric/LinearAlgebra/Util.hs @@ -134,7 +134,7 @@ a & b = vjoin [a,b] {- | horizontal concatenation of real matrices - (0x00a6 broken bar) + (unicode 0x00a6, broken bar) >>> ident 3 ¦ konst 7 (3,4) (3><7) @@ -149,7 +149,7 @@ a ¦ b = fromBlocks [[a,b]] -- | vertical concatenation of real matrices -- --- (0x2014, em dash) +-- (unicode 0x2014, em dash) (——) :: Matrix Double -> Matrix Double -> Matrix Double infixl 2 —— a —— b = fromBlocks [[a],[b]] @@ -179,7 +179,9 @@ infixl 9 ? (?) :: Element t => Matrix t -> [Int] -> Matrix t (?) = flip extractRows --- | (00BF) extract selected columns +-- | extract selected columns +-- +-- (unicode 0x00bf, inverted question mark) infixl 9 ¿ (¿) :: Element t => Matrix t -> [Int] -> Matrix t m ¿ ks = trans . extractRows ks . trans $ m diff --git a/lib/Numeric/Matrix.hs b/lib/Numeric/Matrix.hs index 8397911..e285ff2 100644 --- a/lib/Numeric/Matrix.hs +++ b/lib/Numeric/Matrix.hs @@ -28,6 +28,8 @@ module Numeric.Matrix ( ------------------------------------------------------------------- import Numeric.Container +import qualified Data.Monoid as M +import Data.List(partition) ------------------------------------------------------------------- @@ -69,3 +71,28 @@ instance (Floating a, Container Vector a, Floating (Vector a), Fractional (Matri (**) = liftMatrix2Auto (**) sqrt = liftMatrix sqrt pi = (1><1) [pi] + +-------------------------------------------------------------------------------- + +isScalar m = rows m == 1 && cols m == 1 + +adaptScalarM f1 f2 f3 x y + | isScalar x = f1 (x @@>(0,0) ) y + | isScalar y = f3 x (y @@>(0,0) ) + | otherwise = f2 x y + +instance (Container Vector t, Eq t, Num (Vector t), Product t) => M.Monoid (Matrix t) + where + mempty = 1 + mappend = adaptScalarM scale mXm (flip scale) + + mconcat xs = work (partition isScalar xs) + where + work (ss,[]) = product ss + work (ss,ms) = scale' (product ss) (optimiseMult ms) + scale' x m + | isScalar x && x00 == 1 = m + | otherwise = scale x00 m + where + x00 = x @@> (0,0) + -- cgit v1.2.3