Convert scalars to 1x1 vectors.

[numerical-analysis.git] / src / Linear / Matrix.hs

diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs
index 39576dc41343c2bbad77b745a33bc720d44f666f..ef0e9b6ff3ccba65c310b81126ea915d543e0411 100644 (file)
--- a/src/Linear/Matrix.hs
+++ b/src/Linear/Matrix.hs
@@ -33,13 +33,16 @@ import Data.Vector.Fixed.Internal (Arity, arity, S, Dim)
 import Linear.Vector
 import Normed
 
 import Linear.Vector
 import Normed
 

-import NumericPrelude hiding (abs)
+import NumericPrelude hiding ((*), abs)
+import qualified NumericPrelude as NP ((*))

 import qualified Algebra.Algebraic as Algebraic
 import qualified Algebra.Algebraic as Algebraic

+import Algebra.Algebraic (root)

 import qualified Algebra.Absolute as Absolute
 import qualified Algebra.Additive as Additive
 import qualified Algebra.Ring as Ring
 import Algebra.Absolute (abs)
 import qualified Algebra.Field as Field
 import qualified Algebra.Absolute as Absolute
 import qualified Algebra.Additive as Additive
 import qualified Algebra.Ring as Ring
 import Algebra.Absolute (abs)
 import qualified Algebra.Field as Field

+import qualified Algebra.Module as Module

 import qualified Algebra.RealField as RealField
 import qualified Algebra.RealRing as RealRing
 import qualified Algebra.ToRational as ToRational
 import qualified Algebra.RealField as RealField
 import qualified Algebra.RealRing as RealRing
 import qualified Algebra.ToRational as ToRational


@@ -225,7 +228,7 @@ construct lambda = Mat rows
 --   >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
 --   >>> cholesky m1
 --   ((4.47213595499958,-0.22360679774997896),(0.0,4.466542286825459))
 --   >>> let m1 = fromList [[20,-1], [-1,20]] :: Mat2 Double
 --   >>> cholesky m1
 --   ((4.47213595499958,-0.22360679774997896),(0.0,4.466542286825459))

---   >>> (transpose (cholesky m1)) `mult` (cholesky m1)
+--   >>> (transpose (cholesky m1)) * (cholesky m1)

 --   ((20.000000000000004,-1.0),(-1.0,20.0))
 --
 cholesky :: forall a v w.
 --   ((20.000000000000004,-1.0),(-1.0,20.0))
 --
 cholesky :: forall a v w.


@@ -240,7 +243,7 @@ cholesky m = construct r
     r :: Int -> Int -> a
     r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]])
           | i < j =
     r :: Int -> Int -> a
     r i j | i == j = sqrt(m !!! (i,j) - sum [(r k i)^2 | k <- [0..i-1]])
           | i < j =

-              (((m !!! (i,j)) - sum [(r k i)*(r k j) | k <- [0..i-1]]))/(r i i)
+              (((m !!! (i,j)) - sum [(r k i) NP.* (r k j) | k <- [0..i-1]]))/(r i i)

           | otherwise = 0
 
 
           | otherwise = 0
 
 


@@ -249,12 +252,13 @@ cholesky m = construct r
 --
 --   Examples:
 --
 --
 --   Examples:
 --

---   >>> let m1 = fromList [[1,2,3], [4,5,6]]  :: Mat Vec2D Vec3D Int
---   >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat Vec3D Vec2D Int
---   >>> m1 `mult` m2
+--   >>> let m1 = fromList [[1,2,3], [4,5,6]]  :: Mat D2 D3 Int
+--   >>> let m2 = fromList [[1,2],[3,4],[5,6]] :: Mat D3 D2 Int
+--   >>> m1 * m2

 --   ((22,28),(49,64))
 --
 --   ((22,28),(49,64))
 --

-mult :: (Ring.C a,
+infixl 7 *
+(*) :: (Ring.C a,

          Vector v a,
          Vector w a,
          Vector z a,
          Vector v a,
          Vector w a,
          Vector z a,


@@ -262,10 +266,10 @@ mult :: (Ring.C a,
         => Mat v w a
         -> Mat w z a
         -> Mat v z a
         => Mat v w a
         -> Mat w z a
         -> Mat v z a

-mult m1 m2 = construct lambda
+(*) m1 m2 = construct lambda

   where
     lambda i j =
   where
     lambda i j =

-      sum [(m1 !!! (i,k)) * (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ]
+      sum [(m1 !!! (i,k)) NP.* (m2 !!! (k,j)) | k <- [0..(ncols m1)-1] ]

 
 
 
 
 
 


@@ -288,8 +292,18 @@ instance (Ring.C a,
           Vector w a,
           v ~ w)
          => Ring.C (Mat v w a) where
           Vector w a,
           v ~ w)
          => Ring.C (Mat v w a) where

-    one = Mat (V.replicate $ V.replicate (fromInteger 1))
-    m1 * m2 = m1 `mult` m2
+  -- The first * is ring multiplication, the second is matrix
+  -- multiplication.
+  m1 * m2 = m1 * m1
+
+
+instance (Ring.C a,
+          Vector v (w a),
+          Vector w a)
+         => Module.C a (Mat v w a) where
+  -- We can multiply a matrix by a scalar of the same type as its
+  -- elements.
+  x *> (Mat rows) = Mat $ V.map (V.map (NP.* x)) rows

 
 
 instance (Algebraic.C a,
 
 
 instance (Algebraic.C a,


@@ -299,26 +313,60 @@ instance (Algebraic.C a,
           Vector v a,
           Vector v [a])
          => Normed (Mat v w a) where
           Vector v a,
           Vector v [a])
          => Normed (Mat v w a) where

-   -- Treat the matrix as a big vector.
-   norm_p p (Mat rows) =
-     sqrt $ sum [(fromRational' $ toRational x)^2 | x <- xs]
-     where
-       xs = concat $ V.toList $ V.map V.toList rows
+  -- | Generic p-norms. The usual norm in R^n is (norm_p 2). We treat
+  --   all matrices as big vectors.
+  --
+  --   Examples:
+  --
+  --   >>> let v1 = vec2d (3,4)
+  --   >>> norm_p 1 v1
+  --   7.0
+  --   >>> norm_p 2 v1
+  --   5.0
+  --
+  norm_p p (Mat rows) =
+    (root p') $ sum [(fromRational' $ toRational x)^p' | x <- xs]
+    where
+      p' = toInteger p
+      xs = concat $ V.toList $ V.map V.toList rows

 
 

-   norm_infty m@(Mat rows)
-     | nrows m == 0 || ncols m == 0 = 0
-     | otherwise =
-         fromRational' $ toRational $
-         P.maximum $ V.toList $ V.map (P.maximum . V.toList) rows
+  -- | The infinity norm. We don't use V.maximum here because it
+  --   relies on a type constraint that the vector be non-empty and I
+  --   don't know how to pattern match it away.
+  --
+  --   Examples:
+  --
+  --   >>> let v1 = vec3d (1,5,2)
+  --   >>> norm_infty v1
+  --   5
+  --
+  norm_infty m@(Mat rows)
+    | nrows m == 0 || ncols m == 0 = 0
+    | otherwise =
+        fromRational' $ toRational $
+        P.maximum $ V.toList $ V.map (P.maximum . V.toList) rows

 
 
 
 
 
 -- Vector helpers. We want it to be easy to create low-dimension
 
 
 
 
 
 -- Vector helpers. We want it to be easy to create low-dimension

--- column vectors.
-type Vec a b = Mat a D1 b
+-- column vectors, which are nx1 matrices.

 
 

+-- | Convenient constructor for 2D vectors.
+--
+--   Examples:
+--
+--   >>> import Roots.Simple
+--   >>> let h = 0.5 :: Double
+--   >>> let g1 (Mat (D2 (D1 x) (D1 y))) = 1.0 + h*exp(-(x^2))/(1.0 + y^2)
+--   >>> let g2 (Mat (D2 (D1 x) (D1 y))) = 0.5 + h*atan(x^2 + y^2)
+--   >>> let g u = vec2d ((g1 u), (g2 u))
+--   >>> let u0 = vec2d (1.0, 1.0)
+--   >>> let eps = 1/(10^9)
+--   >>> fixed_point g eps u0
+--   (1.0728549599342185,1.0820591495686167)
+--

 vec2d :: (a,a) -> Mat D2 D1 a
 vec2d (x,y) = Mat (D2 (D1 x) (D1 y))
 
 vec2d :: (a,a) -> Mat D2 D1 a
 vec2d (x,y) = Mat (D2 (D1 x) (D1 y))
 


@@ -328,8 +376,14 @@ vec3d (x,y,z) = Mat (D3 (D1 x) (D1 y) (D1 z))
 vec4d :: (a,a,a,a) -> Mat D4 D1 a
 vec4d (w,x,y,z) = Mat (D4 (D1 w) (D1 x) (D1 y) (D1 z))
 
 vec4d :: (a,a,a,a) -> Mat D4 D1 a
 vec4d (w,x,y,z) = Mat (D4 (D1 w) (D1 x) (D1 y) (D1 z))
 

+-- Since we commandeered multiplication, we need to create 1x1
+-- matrices in order to multiply things.
+scalar :: a -> Mat D1 D1 a
+scalar x = Mat (D1 (D1 x))
+

 dot :: (RealRing.C a,
         Dim w ~ V.N1,
 dot :: (RealRing.C a,
         Dim w ~ V.N1,

+        Dim v ~ S n,

         Vector v a,
         Vector w a,
         Vector w (v a),
         Vector v a,
         Vector w a,
         Vector w (v a),


@@ -337,21 +391,22 @@ dot :: (RealRing.C a,
        => Mat v w a
        -> Mat v w a
        -> a
        => Mat v w a
        -> Mat v w a
        -> a

-v1 `dot` v2 = ((transpose v1) `mult` v2) !!! (0, 0)
+v1 `dot` v2 = ((transpose v1) * v2) !!! (0, 0)

 
 
 -- | The angle between @v1@ and @v2@ in Euclidean space.
 --
 --   Examples:
 --
 
 
 -- | The angle between @v1@ and @v2@ in Euclidean space.
 --
 --   Examples:
 --

---   >>> let v1 = make2d (1.0, 0.0)
---   >>> let v2 = make2d (0.0, 1.0)
+--   >>> let v1 = vec2d (1.0, 0.0)
+--   >>> let v2 = vec2d (0.0, 1.0)

 --   >>> angle v1 v2 == pi/2.0
 --   True
 --
 angle :: (Transcendental.C a,
           RealRing.C a,
           Dim w ~ V.N1,
 --   >>> angle v1 v2 == pi/2.0
 --   True
 --
 angle :: (Transcendental.C a,
           RealRing.C a,
           Dim w ~ V.N1,

+          Dim v ~ S n,

           Vector w (w a),
           Vector v [a],
           Vector v a,
           Vector w (w a),
           Vector v [a],
           Vector v a,


@@ -365,5 +420,5 @@ angle :: (Transcendental.C a,
 angle v1 v2 =
   acos theta
   where
 angle v1 v2 =
   acos theta
   where

-   theta = (recip norms) * (v1 `dot` v2)
-   norms = (norm v1) * (norm v2)
+   theta = (recip norms) NP.* (v1 `dot` v2)
+   norms = (norm v1) NP.* (norm v2)