6 files changed, 143 insertions, 53 deletions

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)
+