mapMatrixWithIndex, Mul, multiply.hs

5 files changed, 164 insertions, 2 deletions

diff --git a/CHANGES b/CHANGES
index 9a07e00..b2d63bf 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,3 +1,9 @@
+0.11.1.0
+========
+
+- exported Mul
+- mapMatrixWithIndex{,M,M_}
+
 0.11.0.0
 ========
 
diff --git a/examples/multiply.hs b/examples/multiply.hs
new file mode 100644
index 0000000..d7c74ee
--- /dev/null
+++ b/examples/multiply.hs
@@ -0,0 +1,100 @@
+{-# LANGUAGE UnicodeSyntax
+           , MultiParamTypeClasses
+           , FunctionalDependencies
+           , FlexibleInstances
+           , FlexibleContexts
+--           , OverlappingInstances
+           , UndecidableInstances #-}
+
+import Numeric.LinearAlgebra
+
+class Scaling a b c | a b -> c where
+ -- ^ 0x22C5    8901    DOT OPERATOR, scaling
+ infixl 7 ⋅
+ (⋅) :: a -> b -> c
+
+class Contraction a b c | a b -> c where
+ -- ^ 0x00D7    215     MULTIPLICATION SIGN     ×, contraction
+ infixl 7 ×
+ (×) :: a -> b -> c
+
+class Outer a b c | a b -> c where
+ -- ^ 0x2297    8855    CIRCLED TIMES   ⊗, outer product (not associative)
+ infixl 7 ⊗
+ (⊗) :: a -> b -> c
+
+
+-------
+
+instance (Num t) => Scaling t t t where
+    (⋅) = (*)
+
+instance Container Vector t => Scaling t (Vector t) (Vector t) where
+    (⋅) = scale
+
+instance Container Vector t => Scaling (Vector t) t (Vector t) where
+    (⋅) = flip scale
+
+instance Container Vector t => Scaling t (Matrix t) (Matrix t) where
+    (⋅) = scale
+
+instance Container Vector t => Scaling (Matrix t) t (Matrix t) where
+    (⋅) = flip scale
+
+
+instance Product t => Contraction (Vector t) (Vector t) t where
+    (×) = dot
+
+instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
+    (×) = mXv
+
+instance Product t => Contraction (Vector t) (Matrix t) (Vector t) where
+    (×) = vXm
+
+instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
+    (×) = mXm
+
+
+--instance Scaling a b c => Contraction a b c where
+--    (×) = (⋅)
+
+-----
+
+instance Product t => Outer (Vector t) (Vector t) (Matrix t) where
+    (⊗) = outer
+
+instance Product t => Outer (Vector t) (Matrix t) (Matrix t) where
+    v ⊗ m = kronecker (asColumn v) m
+
+instance Product t => Outer (Matrix t) (Vector t) (Matrix t) where
+    m ⊗ v = kronecker m (asRow v)
+
+instance Product t => Outer (Matrix t) (Matrix t) (Matrix t) where
+    (⊗) = kronecker
+
+-----
+
+
+v = 3 |> [1..] :: Vector Double
+
+m = (3 >< 3) [1..] :: Matrix Double
+
+s = 3 :: Double
+
+a = s ⋅ v × m × m × v ⋅ s
+
+b = (v ⊗ m) ⊗ (v ⊗ m)
+
+c = v ⊗ m ⊗ v ⊗ m
+
+d = s ⋅ (3 |> [10,20..] :: Vector Double)
+
+main = do
+    print $ scale s v <> m <.> v 
+    print $ scale s v <.> (m <> v)
+    print $ s * (v <> m <.> v)
+    print $ s ⋅ v × m × v
+    print a
+    print (b == c)
+    print d
+
diff --git a/hmatrix.cabal b/hmatrix.cabal
index a8dcb43..a5563aa 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.11.0.5
+Version:            0.11.1.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
@@ -54,6 +54,7 @@ extra-source-files: examples/tests.hs
                     examples/vector.hs
                     examples/monadic.hs
                     examples/bool.hs
+                    examples/multiply.hs
 
 extra-source-files: lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.h,
                     lib/Numeric/LinearAlgebra/LAPACK/clapack.h
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index ab68618..0b23b2f 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -36,6 +36,7 @@ module Data.Packed.Matrix (
     subMatrix, takeRows, dropRows, takeColumns, dropColumns,
     extractRows,
     diagRect, takeDiag,
+    mapMatrix, mapMatrixWithIndex, mapMatrixWithIndexM, mapMatrixWithIndexM_,
     liftMatrix, liftMatrix2, liftMatrix2Auto,fromArray2D
 ) where
 
@@ -44,6 +45,7 @@ import qualified Data.Packed.ST as ST
 import Data.List(transpose,intersperse)
 import Data.Array
 import Foreign.Storable
+import Control.Arrow((***))
 
 -------------------------------------------------------------------
 
@@ -348,3 +350,56 @@ toBlocksEvery r c m = toBlocks rs cs m where
     cs = replicate qc c ++ if rc > 0 then [rc] else []
 
 -------------------------------------------------------------------
+
+mk c g = \k v -> g ((fromIntegral *** fromIntegral) (divMod k c)) v 
+
+{- | 
+
+@ghci> mapMatrixWithIndexM_ (\\(i,j) v -> printf \"m[%.0f,%.0f] = %.f\\n\" i j v :: IO()) ((2><3)[1 :: Double ..])
+m[0,0] = 1
+m[0,1] = 2
+m[0,2] = 3
+m[1,0] = 4
+m[1,1] = 5
+m[1,2] = 6@
+-}
+mapMatrixWithIndexM_
+  :: (Element a, Num a,
+      Functor f, Monad f) =>
+      ((a, a) -> a -> f ()) -> Matrix a -> f ()
+mapMatrixWithIndexM_ g m = mapVectorWithIndexM_ (mk c g) . flatten $ m 
+  where
+    c = cols m
+
+{- |
+
+@ghci> mapMatrixWithIndexM (\\(i,j) v -> Just $ 100*v + 10*i + j) (ident 3:: Matrix Double)
+Just (3><3)
+ [ 100.0,   1.0,   2.0
+ ,  10.0, 111.0,  12.0
+ ,  20.0,  21.0, 122.0 ]@
+-}
+mapMatrixWithIndexM
+  :: (Foreign.Storable.Storable t, 
+      Element a, Num a,
+      Functor f, Monad f) =>
+      ((a, a) -> a -> f t) -> Matrix a -> f (Matrix t)
+mapMatrixWithIndexM g m = fmap (reshape c) . mapVectorWithIndexM (mk c g) . flatten $ m 
+    where
+      c = cols m
+
+{- |
+@ghci> mapMatrixWithIndex (\\(i,j) v -> 100*v + 10*i + j) (ident 3:: Matrix Double)
+(3><3)
+ [ 100.0,   1.0,   2.0
+ ,  10.0, 111.0,  12.0
+ ,  20.0,  21.0, 122.0 ]@
+ -}
+mapMatrixWithIndex :: (Foreign.Storable.Storable t, 
+      Element a, Num a) =>
+      ((a, a) -> a -> t) -> Matrix a -> Matrix t
+mapMatrixWithIndex g = head . mapMatrixWithIndexM (\a b -> [g a b])
+
+mapMatrix :: (Storable a, Storable b) => (a -> b) -> Matrix a -> Matrix b
+mapMatrix f = liftMatrix (mapVector f)
+
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
index 621574e..90155fe 100644
--- a/lib/Numeric/Container.hs
+++ b/lib/Numeric/Container.hs
@@ -36,7 +36,7 @@ module Numeric.Container (
     -- * Matrix product
     Product(..),
     optimiseMult,
-    mXm,mXv,vXm,(<.>),(<>),(<\>),
+    mXm,mXv,vXm,(<.>),Mul(..),(<\>),
     outer, kronecker,
     -- * Random numbers
     RandDist(..),