From: Michael Orlitzky Date: Thu, 21 Feb 2013 03:34:29 +0000 (-0500) Subject: Convert scalars to 1x1 vectors. X-Git-Url: http://gitweb.michael.orlitzky.com/?p=numerical-analysis.git;a=commitdiff_plain;h=2d0ecad8695e129443e53311d6494b2465f1a672 Convert scalars to 1x1 vectors. --- diff --git a/src/Linear/Matrix.hs b/src/Linear/Matrix.hs index 39576dc..ef0e9b6 100644 --- 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 NumericPrelude hiding (abs) +import NumericPrelude hiding ((*), abs) +import qualified NumericPrelude as NP ((*)) 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.Module as Module 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)) --- >>> (transpose (cholesky m1)) `mult` (cholesky m1) +-- >>> (transpose (cholesky m1)) * (cholesky m1) -- ((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 = - (((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 @@ -249,12 +252,13 @@ cholesky m = construct r -- -- 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)) -- -mult :: (Ring.C a, +infixl 7 * +(*) :: (Ring.C 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 -mult m1 m2 = construct lambda +(*) m1 m2 = construct lambda 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 - 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, @@ -299,26 +313,60 @@ instance (Algebraic.C a, 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 --- 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)) @@ -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)) +-- 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, + Dim v ~ S n, 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 -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: -- --- >>> 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, + Dim v ~ S n, Vector w (w a), Vector v [a], Vector v a, @@ -365,5 +420,5 @@ angle :: (Transcendental.C a, 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) diff --git a/src/Linear/Vector.hs b/src/Linear/Vector.hs index 9774dcd..d728334 100644 --- a/src/Linear/Vector.hs +++ b/src/Linear/Vector.hs @@ -63,7 +63,7 @@ instance Vector D4 a where -- -- Examples: -- --- >>> let v1 = Vec2D 1 2 +-- >>> let v1 = D2 1 2 -- >>> v1 ! 1 -- 2 -- @@ -74,7 +74,7 @@ instance Vector D4 a where -- -- Examples: -- --- >>> let v1 = Vec3D 1 2 3 +-- >>> let v1 = D3 1 2 3 -- >>> v1 !? 2 -- Just 3 -- >>> v1 !? 3 @@ -84,54 +84,3 @@ instance Vector D4 a where (!?) v1 idx | idx < 0 || idx >= V.length v1 = Nothing | otherwise = Just $ v1 ! idx - - - - ---instance (RealFloat a, Ord a, Vector v a) => Normed (Vn v a) where - -- | 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 = make3d (1,5,2) - -- >>> norm_infty v1 - -- 5 - -- --- norm_infty (Vn v1) = realToFrac $ V.foldl max 0 v1 - - -- | Generic p-norms. The usual norm in R^n is (norm_p 2). - -- - -- Examples: - -- - -- >>> let v1 = make2d (3,4) - -- >>> norm_p 1 v1 - -- 7.0 - -- >>> norm_p 2 v1 - -- 5.0 - -- --- norm_p p (Vn v1) = --- realToFrac $ root $ V.sum $ V.map (exponentiate . abs) v1 --- where --- exponentiate = (** (fromIntegral p)) --- root = (** (recip (fromIntegral p))) - - - - - --- | Convenient constructor for 2D vectors. --- --- Examples: --- --- >>> import Roots.Simple --- >>> let h = 0.5 :: Double --- >>> let g1 (Vn (Vec2D x y)) = 1.0 + h*exp(-(x^2))/(1.0 + y^2) --- >>> let g2 (Vn (Vec2D x y)) = 0.5 + h*atan(x^2 + y^2) --- >>> let g u = make2d ((g1 u), (g2 u)) --- >>> let u0 = make2d (1.0, 1.0) --- >>> let eps = 1/(10^9) --- >>> fixed_point g eps u0 --- (1.0728549599342185,1.0820591495686167) ---