code refactoring

2 files changed, 168 insertions, 1 deletions

diff --git a/lib/LinearAlgebra/Instances.hs b/lib/LinearAlgebra/Instances.hs
new file mode 100644
index 0000000..3dbe5a7
--- /dev/null
+++ b/lib/LinearAlgebra/Instances.hs
@@ -0,0 +1,140 @@
+{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
+-----------------------------------------------------------------------------
+{- |
+Module      :  LinearAlgebra.Instances
+Copyright   :  (c) Alberto Ruiz 2006
+License     :  GPL-style
+
+Maintainer  :  Alberto Ruiz (aruiz at um dot es)
+Stability   :  provisional
+Portability :  portable
+
+Numeric instances for Vector and 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 LinearAlgebra.Instances(
+) where
+
+import LinearAlgebra.Linear
+import GSL.Vector
+import Data.Packed.Matrix
+import Data.Packed.Vector
+import Complex
+
+adaptScalar f1 f2 f3 x y
+    | dim x == 1 = f1   (x@>0) y
+    | dim y == 1 = f3 x (y@>0)
+    | otherwise = f2 x y
+
+liftMatrix2' :: (Field t, Field a, Field b) => (Vector a -> Vector b -> Vector t) -> Matrix a -> Matrix b -> Matrix t
+liftMatrix2' f m1 m2 | compat' m1 m2 = reshape (max (cols m1) (cols m2)) (f (flatten m1) (flatten m2))
+                     | otherwise    = error "nonconformant matrices in liftMatrix2'"
+
+compat' :: Matrix a -> Matrix b -> Bool
+compat' m1 m2 = rows m1 == 1 && cols m1 == 1
+             || rows m2 == 1 && cols m2 == 1
+             || rows m1 == rows m2 && cols m1 == cols m2
+
+instance (Eq a, Field a) => Eq (Vector a) where
+    a == b = dim a == dim b && toList a == toList b
+
+instance (Linear Vector a) => Num (Vector a) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = liftVector signum
+    abs = liftVector abs
+    fromInteger = fromList . return . fromInteger
+
+instance (Eq a, Field a) => Eq (Matrix a) where
+    a == b = cols a == cols b && flatten a == flatten b
+
+instance (Linear Vector a) => Num (Matrix a) where
+    (+) = liftMatrix2' (+)
+    (-) = liftMatrix2' (-)
+    negate = liftMatrix negate
+    (*) = liftMatrix2' (*)
+    signum = liftMatrix signum
+    abs = liftMatrix abs
+    fromInteger = (1><1) . return . fromInteger
+
+---------------------------------------------------
+
+instance (Linear Vector a) => Fractional (Vector a) where
+    fromRational n = fromList [fromRational n]
+    (/) = adaptScalar f divide g where
+        r `f` v = scaleRecip r v
+        v `g` r = scale (recip r) v
+
+-------------------------------------------------------
+
+instance (Linear Vector a, Fractional (Vector a)) => Fractional (Matrix a) where
+    fromRational n = (1><1) [fromRational n]
+    (/) = liftMatrix2' (/)
+
+---------------------------------------------------------
+
+instance Floating (Vector Double) where
+    sin   = vectorMapR Sin
+    cos   = vectorMapR Cos
+    tan   = vectorMapR Tan
+    asin  = vectorMapR ASin
+    acos  = vectorMapR ACos
+    atan  = vectorMapR ATan
+    sinh  = vectorMapR Sinh
+    cosh  = vectorMapR Cosh
+    tanh  = vectorMapR Tanh
+    asinh = vectorMapR ASinh
+    acosh = vectorMapR ACosh
+    atanh = vectorMapR ATanh
+    exp   = vectorMapR Exp
+    log   = vectorMapR Log
+    sqrt  = vectorMapR Sqrt
+    (**)  = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS))
+    pi    = fromList [pi]
+
+-------------------------------------------------------------
+
+instance Floating (Vector (Complex Double)) where
+    sin   = vectorMapC Sin
+    cos   = vectorMapC Cos
+    tan   = vectorMapC Tan
+    asin  = vectorMapC ASin
+    acos  = vectorMapC ACos
+    atan  = vectorMapC ATan
+    sinh  = vectorMapC Sinh
+    cosh  = vectorMapC Cosh
+    tanh  = vectorMapC Tanh
+    asinh = vectorMapC ASinh
+    acosh = vectorMapC ACosh
+    atanh = vectorMapC ATanh
+    exp   = vectorMapC Exp
+    log   = vectorMapC Log
+    sqrt  = vectorMapC Sqrt
+    (**)  = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS))
+    pi    = fromList [pi]
+
+-----------------------------------------------------------
+
+instance (Linear Vector a, Floating (Vector a)) => Floating (Matrix a) where
+    sin   = liftMatrix sin
+    cos   = liftMatrix cos
+    tan   = liftMatrix tan
+    asin  = liftMatrix asin
+    acos  = liftMatrix acos
+    atan  = liftMatrix atan
+    sinh  = liftMatrix sinh
+    cosh  = liftMatrix cosh
+    tanh  = liftMatrix tanh
+    asinh = liftMatrix asinh
+    acosh = liftMatrix acosh
+    atanh = liftMatrix atanh
+    exp   = liftMatrix exp
+    log   = liftMatrix log
+    (**)  = liftMatrix2' (**)
+    sqrt  = liftMatrix sqrt
+    pi    = (1><1) [pi]
diff --git a/lib/LinearAlgebra/Linear.hs b/lib/LinearAlgebra/Linear.hs
index a2071ed..148bbd3 100644
--- a/lib/LinearAlgebra/Linear.hs
+++ b/lib/LinearAlgebra/Linear.hs
@@ -15,7 +15,8 @@ Portability :  uses ffi
 
 module LinearAlgebra.Linear (
     Linear(..),
-    multiply, dot, outer
+    dot, outer,
+    Mul(..)
 ) where
 
 
@@ -27,10 +28,12 @@ import Complex
 
 class (Field e) => Linear c e where
     scale       :: e -> c e -> c e
+    scaleRecip  :: e -> c e -> c e
     addConstant :: e -> c e -> c e
     add         :: c e -> c e -> c e
     sub         :: c e -> c e -> c e
     mul         :: c e -> c e -> c e
+    divide      :: c e -> c e -> c e
     toComplex   :: RealFloat e => (c e, c e) -> c (Complex e)
     fromComplex :: RealFloat e => c (Complex e) -> (c e, c e)
     comp        :: RealFloat e => c e -> c (Complex e)
@@ -38,10 +41,12 @@ class (Field e) => Linear c e where
 
 instance Linear Vector Double where
     scale = vectorMapValR Scale
+    scaleRecip = vectorMapValR Recip
     addConstant = vectorMapValR AddConstant
     add = vectorZipR Add
     sub = vectorZipR Sub
     mul = vectorZipR Mul
+    divide = vectorZipR Div
     toComplex = Data.Packed.Internal.toComplex
     fromComplex = Data.Packed.Internal.fromComplex
     comp = Data.Packed.Internal.comp
@@ -49,10 +54,12 @@ instance Linear Vector Double where
 
 instance Linear Vector (Complex Double) where
     scale = vectorMapValC Scale
+    scaleRecip = vectorMapValC Recip
     addConstant = vectorMapValC AddConstant
     add = vectorZipC Add
     sub = vectorZipC Sub
     mul = vectorZipC Mul
+    divide = vectorZipC Div
     toComplex = undefined -- can't match
     fromComplex = undefined
     comp = undefined
@@ -60,10 +67,12 @@ instance Linear Vector (Complex Double) where
 
 instance Linear Matrix Double 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
     toComplex = uncurry $ liftMatrix2 $ curry LinearAlgebra.Linear.toComplex
     fromComplex z = (reshape c r, reshape c i)
         where (r,i) = LinearAlgebra.Linear.fromComplex (cdat z)
@@ -73,10 +82,12 @@ instance Linear Matrix Double where
 
 instance Linear Matrix (Complex Double) 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
     toComplex = undefined
     fromComplex = undefined
     comp = undefined
@@ -102,3 +113,19 @@ dot u v = dat (multiply r c) `at` 0
 -}
 outer :: (Field t) => Vector t -> Vector t -> Matrix t
 outer u v = asColumn u `multiply` asRow v
+
+
+class Mul a b c | a b -> c where
+ infixl 7 <>
+ -- | matrix product
+ (<>) :: Field t => a t -> b t -> c t
+
+instance Mul Matrix Matrix Matrix where
+    (<>) = multiply
+
+instance Mul Matrix Vector Vector where
+    (<>) m v = flatten $ m <> (asColumn v)
+
+instance Mul Vector Matrix Vector where
+    (<>) v m = flatten $ (asRow v) <> m
+