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
-- >>> 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.
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
--
-- 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,
=> 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] ]
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,
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))
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),
=> 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,
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)
--
-- Examples:
--
--- >>> let v1 = Vec2D 1 2
+-- >>> let v1 = D2 1 2
-- >>> v1 ! 1
-- 2
--
--
-- Examples:
--
--- >>> let v1 = Vec3D 1 2 3
+-- >>> let v1 = D3 1 2 3
-- >>> v1 !? 2
-- Just 3
-- >>> v1 !? 3
(!?) 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)
---
projects / numerical-analysis.git / commitdiff