src/Linear/Vector.hs

   1 {-# LANGUAGE FlexibleContexts #-}
   2 {-# LANGUAGE FlexibleInstances #-}
   3 {-# LANGUAGE MultiParamTypeClasses #-}
   4 {-# LANGUAGE ScopedTypeVariables #-}
   5 {-# LANGUAGE TypeFamilies #-}
   6
   7 module Linear.Vector
   8 where
   9
  10 import Data.List (intercalate)
  11 import Data.Vector.Fixed (
  12   Dim,
  13   Fun(..),
  14   N1,
  15   N2,
  16   N3,
  17   N4,
  18   Vector(..),
  19   construct,
  20   inspect,
  21   toList,
  22   )
  23 import qualified Data.Vector.Fixed as V (
  24   length,
  25   )
  26
  27 import Normed
  28
  29
  30 -- * Low-dimension vector wrappers.
  31 --
  32 -- These wrappers are instances of 'Vector', so they inherit all of
  33 -- the userful instances defined above. But, they use fixed
  34 -- constructors, so you can pattern match out the individual
  35 -- components.
  36
  37 data D1 a = D1 a
  38 type instance Dim D1 = N1
  39 instance Vector D1 a where
  40   inspect (D1 x) (Fun f) = f x
  41   construct = Fun D1
  42
  43 data D2 a = D2 a a
  44 type instance Dim D2 = N2
  45 instance Vector D2 a where
  46   inspect (D2 x y) (Fun f) = f x y
  47   construct = Fun D2
  48
  49 data D3 a = D3 a a a
  50 type instance Dim D3 = N3
  51 instance Vector D3 a where
  52   inspect (D3 x y z) (Fun f) = f x y z
  53   construct = Fun D3
  54
  55 data D4 a = D4 a a a a
  56 type instance Dim D4 = N4
  57 instance Vector D4 a where
  58   inspect (D4 w x y z) (Fun f) = f w x y z
  59   construct = Fun D4
  60
  61
  62 -- | Unsafe indexing.
  63 --
  64 --   Examples:
  65 --
  66 --   >>> let v1 = Vec2D 1 2
  67 --   >>> v1 ! 1
  68 --   2
  69 --
  70 (!) :: (Vector v a) => v a -> Int -> a
  71 (!) v1 idx = (toList v1) !! idx
  72
  73 -- | Safe indexing.
  74 --
  75 --   Examples:
  76 --
  77 --   >>> let v1 = Vec3D 1 2 3
  78 --   >>> v1 !? 2
  79 --   Just 3
  80 --   >>> v1 !? 3
  81 --   Nothing
  82 --
  83 (!?) :: (Vector v a) => v a -> Int -> Maybe a
  84 (!?) v1 idx
  85   | idx < 0 || idx >= V.length v1 = Nothing
  86   | otherwise                     = Just $ v1 ! idx
  87
  88
  89
  90
  91 --instance (RealFloat a, Ord a, Vector v a) => Normed (Vn v a) where
  92   -- | The infinity norm. We don't use V.maximum here because it
  93   --   relies on a type constraint that the vector be non-empty and I
  94   --   don't know how to pattern match it away.
  95   --
  96   --   Examples:
  97   --
  98   --   >>> let v1 = make3d (1,5,2)
  99   --   >>> norm_infty v1
 100   --   5
 101   --
 102 --  norm_infty (Vn v1) = realToFrac $ V.foldl max 0 v1
 103
 104   -- | Generic p-norms. The usual norm in R^n is (norm_p 2).
 105   --
 106   --   Examples:
 107   --
 108   --   >>> let v1 = make2d (3,4)
 109   --   >>> norm_p 1 v1
 110   --   7.0
 111   --   >>> norm_p 2 v1
 112   --   5.0
 113   --
 114 --  norm_p p (Vn v1) =
 115 --    realToFrac $ root $ V.sum $ V.map (exponentiate . abs) v1
 116 --    where
 117 --      exponentiate = (** (fromIntegral p))
 118 --      root = (** (recip (fromIntegral p)))
 119
 120
 121
 122
 123
 124 -- | Convenient constructor for 2D vectors.
 125 --
 126 --   Examples:
 127 --
 128 --   >>> import Roots.Simple
 129 --   >>> let h = 0.5 :: Double
 130 --   >>> let g1 (Vn (Vec2D x y)) = 1.0 + h*exp(-(x^2))/(1.0 + y^2)
 131 --   >>> let g2 (Vn (Vec2D x y)) = 0.5 + h*atan(x^2 + y^2)
 132 --   >>> let g u = make2d ((g1 u), (g2 u))
 133 --   >>> let u0 = make2d (1.0, 1.0)
 134 --   >>> let eps = 1/(10^9)
 135 --   >>> fixed_point g eps u0
 136 --   (1.0728549599342185,1.0820591495686167)
 137 --