1 files changed, 109 insertions, 48 deletions

diff --git a/lib/Numeric/Matrix.hs b/lib/Numeric/Matrix.hs
index fa3f94a..6ba870f 100644
--- a/lib/Numeric/Matrix.hs
+++ b/lib/Numeric/Matrix.hs
@@ -3,6 +3,7 @@
 {-# LANGUAGE FlexibleInstances #-}
 {-# LANGUAGE UndecidableInstances #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FunctionalDependencies #-}
 
 -----------------------------------------------------------------------------
 -- |
@@ -14,12 +15,22 @@
 -- Stability   :  provisional
 -- Portability :  portable
 --
--- Numeric instances and functions for 'Data.Packed.Matrix's
+-- Numeric instances and functions for 'Matrix'.
+-- In the context of the standard numeric operators, one-component
+-- vectors and matrices automatically expand to match the dimensions of the other operand.
 --
 -----------------------------------------------------------------------------
 
 module Numeric.Matrix (
-                       module Data.Packed.Matrix,
+    -- * Basic functions
+    module Data.Packed.Matrix,
+    module Numeric.Vector,
+    --module Numeric.Container,
+    -- * Operators
+    (<>), (<\>),
+    -- * Deprecated
+    (.*),(*/),(<|>),(<->),
+    vectorMax,vectorMin
                       ) where
 
 -------------------------------------------------------------------
@@ -28,18 +39,18 @@ import Data.Packed.Vector
 import Data.Packed.Matrix
 import Numeric.Container
 --import Numeric.LinearAlgebra.Linear
-import Numeric.Vector()
+import Numeric.Vector
+import Numeric.LinearAlgebra.Algorithms
+--import Control.Monad(ap)
 
-import Control.Monad(ap)
-
-import Control.Arrow((***))
+--import Control.Arrow((***))
 
 -------------------------------------------------------------------
 
-instance Linear Matrix a => Eq (Matrix a) where
+instance Container Matrix a => Eq (Matrix a) where
     (==) = equal
 
-instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
+instance (Container Matrix a, Num (Vector a)) => Num (Matrix a) where
     (+) = liftMatrix2Auto (+)
     (-) = liftMatrix2Auto (-)
     negate = liftMatrix negate
@@ -50,13 +61,13 @@ instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where
 
 ---------------------------------------------------
 
-instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
+instance (Container Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
     fromRational n = (1><1) [fromRational n]
     (/) = liftMatrix2Auto (/)
 
 ---------------------------------------------------------
 
-instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where
+instance (Container Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where
     sin   = liftMatrix sin
     cos   = liftMatrix cos
     tan   = liftMatrix tan
@@ -75,43 +86,93 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati
     sqrt  = liftMatrix sqrt
     pi    = (1><1) [pi]
 
----------------------------------------------------------------
-
-instance NumericContainer Matrix where
-    toComplex = uncurry $ liftMatrix2 $ curry toComplex
-    fromComplex z = (reshape c *** reshape c) . fromComplex . flatten $ z
-        where c = cols z
-    complex' = liftMatrix complex'
-    conj = liftMatrix conj
---    cmap f = liftMatrix (cmap f)
-    single' = liftMatrix single'
-    double' = liftMatrix double'
-
----------------------------------------------------------------
-
-instance (Linear Vector a, Container Matrix a) => Linear Matrix a where
-    scale x = liftMatrix (scale x)
-    scaleRecip x = liftMatrix (scaleRecip x)
-    addConstant x = liftMatrix (addConstant x)
-    add = liftMatrix2 add
-    sub = liftMatrix2 sub
-    mul = liftMatrix2 mul
-    divide = liftMatrix2 divide
-    equal a b = cols a == cols b && flatten a `equal` flatten b
-    scalar x = (1><1) [x]
-    --
-instance (Container Vector a) => Container Matrix a where
-    cmap f = liftMatrix (mapVector f)
-    atIndex = (@@>)
-    minIndex m = let (r,c) = (rows m,cols m)
-                     i = (minIndex $ flatten m)
-                 in (i `div` c,(i `mod` c) + 1)
-    maxIndex m = let (r,c) = (rows m,cols m)
-                     i = (maxIndex $ flatten m)
-                 in (i `div` c,(i `mod` c) + 1)
-    minElement = ap (@@>) minIndex
-    maxElement = ap (@@>) maxIndex
-    sumElements = sumElements . flatten
-    prodElements = prodElements . flatten
+--------------------------------------------------------
+
+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)
+
+instance Mul Vector Matrix Vector where
+    (<>) v m = flatten $ (asRow v) <> m
 
 ----------------------------------------------------
+
+{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-}
+
+-- -- | @x .* a = scale x a@
+-- (.*) :: (Linear c a) => a -> c a -> c a
+infixl 7 .*
+a .* x = scale a x
+
+----------------------------------------------------
+
+{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-}
+
+-- -- | @a *\/ x = scale (recip x) a@
+-- (*/) :: (Linear c a) => c a -> a -> c a
+infixl 7 */
+v */ x = scale (recip x) v
+
+-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD).
+(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a
+infixl 7 <\>
+m <\> v = flatten (linearSolveSVD m (reshape 1 v))
+
+------------------------------------------------
+
+{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-}
+{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-}
+
+class Joinable a b where
+    joinH :: Element t => a t -> b t -> Matrix t
+    joinV :: Element t => a t -> b t -> Matrix t
+
+instance Joinable Matrix Matrix where
+    joinH m1 m2 = fromBlocks [[m1,m2]]
+    joinV m1 m2 = fromBlocks [[m1],[m2]]
+
+instance Joinable Matrix Vector where
+    joinH m v = joinH m (asColumn v)
+    joinV m v = joinV m (asRow v)
+
+instance Joinable Vector Matrix where
+    joinH v m = joinH (asColumn v) m
+    joinV v m = joinV (asRow v) m
+
+infixl 4 <|>
+infixl 3 <->
+
+{-- - | Horizontal concatenation of matrices and vectors:
+
+@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0]
+(6><4)
+ [ 1.0, 0.0, 0.0, 1.0
+ , 0.0, 1.0, 0.0, 2.0
+ , 0.0, 0.0, 1.0, 3.0
+ , 3.0, 0.0, 0.0, 4.0
+ , 0.0, 3.0, 0.0, 5.0
+ , 0.0, 0.0, 3.0, 6.0 ]@
+-}
+-- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
+a <|> b = joinH a b
+
+-- -- | Vertical concatenation of matrices and vectors.
+-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t
+a <-> b = joinV a b
+
+-------------------------------------------------------------------
+
+{-# DEPRECATED vectorMin "use minElement" #-}
+vectorMin :: (Container Vector t, Element t) => Vector t -> t
+vectorMin = minElement
+
+{-# DEPRECATED vectorMax "use maxElement" #-}
+vectorMax :: (Container Vector t, Element t) => Vector t -> t
+vectorMax = maxElement