instances moved to base

4 files changed, 405 insertions, 0 deletions

diff --git a/packages/base/hmatrix-base.cabal b/packages/base/hmatrix-base.cabal
index 8571e7a..8eeb97e 100644
--- a/packages/base/hmatrix-base.cabal
+++ b/packages/base/hmatrix-base.cabal
@@ -39,12 +39,15 @@ library
                         Numeric.LinearAlgebra.LAPACK
                         Data.Packed.Numeric
                         Numeric.Vectorized
+                        Numeric.Vector
+                        Numeric.Matrix
                         
     other-modules:      Data.Packed.Internal,
                         Data.Packed.Internal.Common,
                         Data.Packed.Internal.Signatures,
                         Data.Packed.Internal.Vector,
                         Data.Packed.Internal.Matrix
+                        Numeric.Chain
 
     C-sources:          src/C/lapack-aux.c
                         src/C/vector-aux.c
diff --git a/packages/base/src/Numeric/Chain.hs b/packages/base/src/Numeric/Chain.hs
new file mode 100644
index 0000000..fbdb01b
--- /dev/null
+++ b/packages/base/src/Numeric/Chain.hs
@@ -0,0 +1,143 @@
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Chain
+-- Copyright   :  (c) Vivian McPhail 2010
+-- License     :  GPL-style
+--
+-- Maintainer  :  Vivian McPhail <haskell.vivian.mcphail <at> gmail.com>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- optimisation of association order for chains of matrix multiplication
+--
+-----------------------------------------------------------------------------
+
+module Numeric.Chain (
+                      optimiseMult,
+                     ) where
+
+import Data.Maybe
+
+import Data.Packed.Matrix
+import Data.Packed.Numeric
+
+import qualified Data.Array.IArray as A
+
+-----------------------------------------------------------------------------
+{- | 
+     Provide optimal association order for a chain of matrix multiplications 
+     and apply the multiplications.
+
+     The algorithm is the well-known O(n\^3) dynamic programming algorithm
+     that builds a pyramid of optimal associations.
+
+> m1, m2, m3, m4 :: Matrix Double
+> m1 = (10><15) [1..]
+> m2 = (15><20) [1..]
+> m3 = (20><5) [1..]
+> m4 = (5><10) [1..]
+
+> >>> optimiseMult [m1,m2,m3,m4]
+
+will perform @((m1 `multiply` (m2 `multiply` m3)) `multiply` m4)@
+
+The naive left-to-right multiplication would take @4500@ scalar multiplications
+whereas the optimised version performs @2750@ scalar multiplications.  The complexity
+in this case is 32 (= 4^3/2) * (2 comparisons, 3 scalar multiplications, 3 scalar additions,
+5 lookups, 2 updates) + a constant (= three table allocations)
+-}
+optimiseMult :: Product t => [Matrix t] -> Matrix t
+optimiseMult = chain
+
+-----------------------------------------------------------------------------
+
+type Matrices a = A.Array Int (Matrix a)
+type Sizes      = A.Array Int (Int,Int)
+type Cost       = A.Array Int (A.Array Int (Maybe Int))
+type Indexes    = A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int))))
+
+update :: A.Array Int (A.Array Int a) -> (Int,Int) -> a -> A.Array Int (A.Array Int a)
+update a (r,c) e = a A.// [(r,(a A.! r) A.// [(c,e)])]
+
+newWorkSpaceCost :: Int -> A.Array Int (A.Array Int (Maybe Int))
+newWorkSpaceCost n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n]
+   where subArray i = A.listArray (1,i) (repeat Nothing)
+
+newWorkSpaceIndexes :: Int -> A.Array Int (A.Array Int (Maybe ((Int,Int),(Int,Int))))
+newWorkSpaceIndexes n = A.array (1,n) $ map (\i -> (i, subArray i)) [1..n]
+   where subArray i = A.listArray (1,i) (repeat Nothing)
+
+matricesToSizes :: [Matrix a] -> Sizes
+matricesToSizes ms = A.listArray (1,length ms) $ map (\m -> (rows m,cols m)) ms
+
+chain :: Product a => [Matrix a] -> Matrix a
+chain []  = error "chain: zero matrices to multiply"
+chain [m] = m
+chain [ml,mr] = ml `multiply` mr
+chain ms = let ln = length ms
+               ma = A.listArray (1,ln) ms
+               mz = matricesToSizes ms
+               i = chain_cost mz
+           in chain_paren (ln,ln) i ma
+
+chain_cost :: Sizes -> Indexes
+chain_cost mz = let (_,u) = A.bounds mz
+                    cost = newWorkSpaceCost u
+                    ixes = newWorkSpaceIndexes u
+                    (_,_,i) =  foldl chain_cost' (mz,cost,ixes) (order u)
+                in i
+
+chain_cost' :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes)
+chain_cost' sci@(mz,cost,ixes) (r,c) 
+    | c == 1                     = let cost' = update cost (r,c) (Just 0)
+                                       ixes' = update ixes (r,c) (Just ((r,c),(r,c)))
+                                       in (mz,cost',ixes')
+    | otherwise                  = minimum_cost sci (r,c)
+
+minimum_cost :: (Sizes,Cost,Indexes) -> (Int,Int) -> (Sizes,Cost,Indexes)
+minimum_cost sci fu = foldl (smaller_cost fu) sci (fulcrum_order fu)
+
+smaller_cost :: (Int,Int) -> (Sizes,Cost,Indexes) -> ((Int,Int),(Int,Int)) -> (Sizes,Cost,Indexes)
+smaller_cost (r,c) (mz,cost,ixes) ix@((lr,lc),(rr,rc)) = let op_cost =   fromJust ((cost A.! lr) A.! lc)
+                                                                       + fromJust ((cost A.! rr) A.! rc)
+                                                                       + fst (mz A.! (lr-lc+1))
+                                                                         * snd (mz A.! lc)
+                                                                         * snd (mz A.! rr)
+                                                             cost' = (cost A.! r) A.! c
+                                                         in case cost' of
+                                                                       Nothing -> let cost'' = update cost (r,c) (Just op_cost)
+                                                                                      ixes'' = update ixes (r,c) (Just ix)
+                                                                                  in (mz,cost'',ixes'')
+                                                                       Just ct -> if op_cost < ct then
+                                                                                  let cost'' = update cost (r,c) (Just op_cost)
+                                                                                      ixes'' = update ixes (r,c) (Just ix)
+                                                                                  in (mz,cost'',ixes'')
+                                                                                  else (mz,cost,ixes)
+                                                                         
+
+fulcrum_order (r,c) = let fs' = zip (repeat r) [1..(c-1)]
+                      in map (partner (r,c)) fs'
+
+partner (r,c) (a,b) = ((r-b, c-b), (a,b))
+
+order 0 = []
+order n = order (n-1) ++ zip (repeat n) [1..n]
+
+chain_paren :: Product a => (Int,Int) -> Indexes -> Matrices a -> Matrix a
+chain_paren (r,c) ixes ma = let ((lr,lc),(rr,rc)) = fromJust $ (ixes A.! r) A.! c
+                            in if lr == rr && lc == rc then (ma A.! lr)
+                               else (chain_paren (lr,lc) ixes ma) `multiply` (chain_paren (rr,rc) ixes ma) 
+
+--------------------------------------------------------------------------
+
+{- TESTS
+
+-- optimal association is ((m1*(m2*m3))*m4)
+m1, m2, m3, m4 :: Matrix Double
+m1 = (10><15) [1..]
+m2 = (15><20) [1..]
+m3 = (20><5) [1..]
+m4 = (5><10) [1..]
+
+-}
+
diff --git a/packages/base/src/Numeric/Matrix.hs b/packages/base/src/Numeric/Matrix.hs
new file mode 100644
index 0000000..3478aae
--- /dev/null
+++ b/packages/base/src/Numeric/Matrix.hs
@@ -0,0 +1,100 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Matrix
+-- Copyright   :  (c) Alberto Ruiz 2010
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Provides instances of standard classes 'Show', 'Read', 'Eq',
+-- 'Num', 'Fractional', and 'Floating' for 'Matrix'.
+--
+-- In arithmetic operations one-component
+-- vectors and matrices automatically expand to match the dimensions of the other operand.
+
+-----------------------------------------------------------------------------
+
+module Numeric.Matrix (
+                      ) where
+
+-------------------------------------------------------------------
+
+import Data.Packed
+import Data.Packed.Numeric
+import qualified Data.Monoid as M
+import Data.List(partition)
+import Numeric.Chain
+
+-------------------------------------------------------------------
+
+instance Container Matrix a => Eq (Matrix a) where
+    (==) = equal
+
+instance (Container Matrix a, Num (Vector a)) => Num (Matrix a) where
+    (+) = liftMatrix2Auto (+)
+    (-) = liftMatrix2Auto (-)
+    negate = liftMatrix negate
+    (*) = liftMatrix2Auto (*)
+    signum = liftMatrix signum
+    abs = liftMatrix abs
+    fromInteger = (1><1) . return . fromInteger
+
+---------------------------------------------------
+
+instance (Container Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where
+    fromRational n = (1><1) [fromRational n]
+    (/) = liftMatrix2Auto (/)
+
+---------------------------------------------------------
+
+instance (Floating a, Container Vector a, Floating (Vector a), Fractional (Matrix 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
+    (**)  = 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)
+
diff --git a/packages/base/src/Numeric/Vector.hs b/packages/base/src/Numeric/Vector.hs
new file mode 100644
index 0000000..2769cd9
--- /dev/null
+++ b/packages/base/src/Numeric/Vector.hs
@@ -0,0 +1,159 @@
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE UndecidableInstances #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+-----------------------------------------------------------------------------
+-- |
+-- Module      :  Numeric.Vector
+-- Copyright   :  (c) Alberto Ruiz 2011
+-- License     :  GPL-style
+--
+-- Maintainer  :  Alberto Ruiz <aruiz@um.es>
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Provides instances of standard classes 'Show', 'Read', 'Eq',
+-- 'Num', 'Fractional',  and 'Floating' for 'Vector'.
+-- 
+-----------------------------------------------------------------------------
+
+module Numeric.Vector () where
+
+import Numeric.Vectorized
+import Data.Packed.Vector
+import Data.Packed.Numeric
+
+-------------------------------------------------------------------
+
+adaptScalar f1 f2 f3 x y
+    | dim x == 1 = f1   (x@>0) y
+    | dim y == 1 = f3 x (y@>0)
+    | otherwise = f2 x y
+
+------------------------------------------------------------------
+
+instance Num (Vector Float) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapF Sign
+    abs = vectorMapF Abs
+    fromInteger = fromList . return . fromInteger
+
+instance Num (Vector Double) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapR Sign
+    abs = vectorMapR Abs
+    fromInteger = fromList . return . fromInteger
+
+instance Num (Vector (Complex Double)) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapC Sign
+    abs = vectorMapC Abs
+    fromInteger = fromList . return . fromInteger
+
+instance Num (Vector (Complex Float)) where
+    (+) = adaptScalar addConstant add (flip addConstant)
+    negate = scale (-1)
+    (*) = adaptScalar scale mul (flip scale)
+    signum = vectorMapQ Sign
+    abs = vectorMapQ Abs
+    fromInteger = fromList . return . fromInteger
+
+---------------------------------------------------
+
+instance (Container Vector a, Num (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 Floating (Vector Float) where
+    sin   = vectorMapF Sin
+    cos   = vectorMapF Cos
+    tan   = vectorMapF Tan
+    asin  = vectorMapF ASin
+    acos  = vectorMapF ACos
+    atan  = vectorMapF ATan
+    sinh  = vectorMapF Sinh
+    cosh  = vectorMapF Cosh
+    tanh  = vectorMapF Tanh
+    asinh = vectorMapF ASinh
+    acosh = vectorMapF ACosh
+    atanh = vectorMapF ATanh
+    exp   = vectorMapF Exp
+    log   = vectorMapF Log
+    sqrt  = vectorMapF Sqrt
+    (**)  = adaptScalar (vectorMapValF PowSV) (vectorZipF Pow) (flip (vectorMapValF PowVS))
+    pi    = fromList [pi]
+
+-------------------------------------------------------------
+
+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 Floating (Vector (Complex Float)) where
+    sin   = vectorMapQ Sin
+    cos   = vectorMapQ Cos
+    tan   = vectorMapQ Tan
+    asin  = vectorMapQ ASin
+    acos  = vectorMapQ ACos
+    atan  = vectorMapQ ATan
+    sinh  = vectorMapQ Sinh
+    cosh  = vectorMapQ Cosh
+    tanh  = vectorMapQ Tanh
+    asinh = vectorMapQ ASinh
+    acosh = vectorMapQ ACosh
+    atanh = vectorMapQ ATanh
+    exp   = vectorMapQ Exp
+    log   = vectorMapQ Log
+    sqrt  = vectorMapQ Sqrt
+    (**)  = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS))
+    pi    = fromList [pi]
+