From a858bf910291b63603a226c3190ecb36de01b5ba Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 8 Sep 2010 08:14:27 +0000 Subject: re-export changes --- lib/Data/Packed.hs | 10 +- lib/Data/Packed/Random.hs | 7 +- lib/Graphics/Plot.hs | 4 +- lib/Numeric/Container.hs | 418 +++++++++++++++++++------------- lib/Numeric/LinearAlgebra.hs | 10 +- lib/Numeric/LinearAlgebra/Algorithms.hs | 9 +- lib/Numeric/LinearAlgebra/Instances.hs | 276 --------------------- lib/Numeric/LinearAlgebra/Interface.hs | 118 --------- lib/Numeric/LinearAlgebra/LAPACK.hs | 3 +- lib/Numeric/LinearAlgebra/Linear.hs | 160 ------------ lib/Numeric/Matrix.hs | 157 ++++++++---- lib/Numeric/Vector.hs | 141 ++--------- 12 files changed, 403 insertions(+), 910 deletions(-) delete mode 100644 lib/Numeric/LinearAlgebra/Instances.hs delete mode 100644 lib/Numeric/LinearAlgebra/Interface.hs delete mode 100644 lib/Numeric/LinearAlgebra/Linear.hs (limited to 'lib') diff --git a/lib/Data/Packed.hs b/lib/Data/Packed.hs index 8128667..bfc2d8b 100644 --- a/lib/Data/Packed.hs +++ b/lib/Data/Packed.hs @@ -16,11 +16,13 @@ The Vector and Matrix types and some utilities. module Data.Packed ( module Data.Packed.Vector, module Data.Packed.Matrix, - module Data.Packed.Random, - module Data.Complex +-- module Numeric.Conversion, +-- module Data.Packed.Random, +-- module Data.Complex ) where import Data.Packed.Vector import Data.Packed.Matrix -import Data.Packed.Random -import Data.Complex +--import Data.Packed.Random +--import Data.Complex +--import Numeric.Conversion \ No newline at end of file diff --git a/lib/Data/Packed/Random.hs b/lib/Data/Packed/Random.hs index 0bc5350..b30b299 100644 --- a/lib/Data/Packed/Random.hs +++ b/lib/Data/Packed/Random.hs @@ -20,11 +20,12 @@ module Data.Packed.Random ( ) where import Numeric.GSL.Vector -import Data.Packed.Matrix -import Numeric.Vector -import Numeric.LinearAlgebra.Linear +import Data.Packed +import Numeric.Container +import Data.Packed.Internal(constantD) import Numeric.LinearAlgebra.Algorithms +constant k v = constantD k v -- | Obtains a matrix whose rows are pseudorandom samples from a multivariate -- Gaussian distribution. diff --git a/lib/Graphics/Plot.hs b/lib/Graphics/Plot.hs index ae88a92..c5b5a4c 100644 --- a/lib/Graphics/Plot.hs +++ b/lib/Graphics/Plot.hs @@ -28,9 +28,7 @@ module Graphics.Plot( ) where -import Data.Packed -import Numeric.Vector -import Numeric.LinearAlgebra.Linear +import Numeric.Matrix import Data.List(intersperse) import System.Process (system) diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs index d6b1342..fc4fb0d 100644 --- a/lib/Numeric/Container.hs +++ b/lib/Numeric/Container.hs @@ -3,6 +3,7 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- -- | @@ -19,103 +20,35 @@ ----------------------------------------------------------------------------- module Numeric.Container ( - Linear(..), - Container(..), RealElement, Precision(..), NumericContainer(..), comp, - Convert(..), --AutoReal(..), - RealOf, ComplexOf, SingleOf, DoubleOf, - --- ElementOf, + --Linear(..), + Container(..), + Vectors(..), + Product(..), + mXm,mXv,vXm, + outer, kronecker, + + RealElement, --Precision, + ComplexContainer(toComplex,fromComplex,comp,conj), + Convert(..), --AutoReal(..), + RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf, - module Data.Complex ) where -import Data.Packed.Vector -import Data.Packed.Matrix -import Data.Packed.Internal.Vector ---import Data.Packed.Internal.Matrix ---import qualified Data.Packed.ST as ST +import Data.Packed +import Numeric.Conversion +import Data.Packed.Internal +import Numeric.GSL.Vector ---import Control.Arrow((***)) import Data.Complex +import Control.Monad(ap) +--import Control.Arrow((***)) -------------------------------------------------------------------- - --- | Supported single-double precision type pairs -class (Element s, Element d) => Precision s d | s -> d, d -> s where - double2FloatG :: Vector d -> Vector s - float2DoubleG :: Vector s -> Vector d - -instance Precision Float Double where - double2FloatG = double2FloatV - float2DoubleG = float2DoubleV - -instance Precision (Complex Float) (Complex Double) where - double2FloatG = asComplex . double2FloatV . asReal - float2DoubleG = asComplex . float2DoubleV . asReal - --- | Supported real types -class (Element t, Element (Complex t), RealFloat t --- , RealOf t ~ t, RealOf (Complex t) ~ t - ) - => RealElement t - -instance RealElement Double - -instance RealElement Float - --- | Conversion utilities -class NumericContainer c where - toComplex :: (RealElement e) => (c e, c e) -> c (Complex e) - fromComplex :: (RealElement e) => c (Complex e) -> (c e, c e) - complex' :: (RealElement e) => c e -> c (Complex e) - conj :: (RealElement e) => c (Complex e) -> c (Complex e) --- cmap :: (Element a, Element b) => (a -> b) -> c a -> c b - single' :: Precision a b => c b -> c a - double' :: Precision a b => c a -> c b - --- | a synonym for "complex'" -comp :: (NumericContainer c, RealElement e) => c e -> c (Complex e) -comp x = complex' x +import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) ------------------------------------------------------------------- -type family RealOf x - -type instance RealOf Double = Double -type instance RealOf (Complex Double) = Double - -type instance RealOf Float = Float -type instance RealOf (Complex Float) = Float - -type family ComplexOf x - -type instance ComplexOf Double = Complex Double -type instance ComplexOf (Complex Double) = Complex Double - -type instance ComplexOf Float = Complex Float -type instance ComplexOf (Complex Float) = Complex Float - -type family SingleOf x - -type instance SingleOf Double = Float -type instance SingleOf Float = Float - -type instance SingleOf (Complex a) = Complex (SingleOf a) - -type family DoubleOf x - -type instance DoubleOf Double = Double -type instance DoubleOf Float = Double - -type instance DoubleOf (Complex a) = Complex (DoubleOf a) - -type family ElementOf c - -type instance ElementOf (Vector a) = a -type instance ElementOf (Matrix a) = a - type family IndexOf c type instance IndexOf Vector = Int @@ -123,75 +56,16 @@ type instance IndexOf Matrix = (Int,Int) ------------------------------------------------------------------- -class (Element t, Element (RealOf t)) => Convert t where - real :: NumericContainer c => c (RealOf t) -> c t - complex :: NumericContainer c => c t -> c (ComplexOf t) - single :: NumericContainer c => c t -> c (SingleOf t) - double :: NumericContainer c => c t -> c (DoubleOf t) - - -instance Convert Double where - real = id - complex = complex' - single = single' - double = id - -instance Convert Float where - real = id - complex = complex' - single = id - double = double' - -instance Convert (Complex Double) where - real = complex' - complex = id - single = single' - double = id - -instance Convert (Complex Float) where - real = complex' - complex = id - single = id - double = double' - -------------------------------------------------------------------- - --- | to be replaced by Convert -class Convert t => AutoReal t where - real'' :: NumericContainer c => c Double -> c t - complex'' :: NumericContainer c => c t -> c (Complex Double) - - -instance AutoReal Double where - real'' = real - complex'' = complex - -instance AutoReal (Complex Double) where - real'' = real - complex'' = complex - -instance AutoReal Float where - real'' = real . single - complex'' = double . complex - -instance AutoReal (Complex Float) where - real'' = real . single - complex'' = double . complex - -------------------------------------------------------------------- - -- | Basic element-by-element functions for numeric containers class (Element e) => Container c e where -{- + -- | create a structure with a single element scalar :: e -> c e - -- | multiply every element by a scalar scale :: e -> c e -> c e -- | scale the element by element reciprocal of the object: -- -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ scaleRecip :: e -> c e -> c e - -- | add a constant to each element addConstant :: e -> c e -> c e add :: c e -> c e -> c e sub :: c e -> c e -> c e @@ -200,7 +74,7 @@ class (Element e) => Container c e where -- | element by element division divide :: c e -> c e -> c e equal :: c e -> c e -> Bool --} + -- | cannot implement instance Functor because of Element class constraint cmap :: (Element a, Element b) => (a -> b) -> c a -> c b -- @@ -217,23 +91,239 @@ class (Element e) => Container c e where sumElements :: c e -> e prodElements :: c e -> e --- | Basic element-by-element functions. -class (Element e, Container c e) => Linear c e where - -- | create a structure with a single element - scalar :: e -> c e - scale :: e -> c e -> c e - -- | scale the element by element reciprocal of the object: - -- - -- @scaleRecip 2 (fromList [5,i]) == 2 |> [0.4 :+ 0.0,0.0 :+ (-2.0)]@ - scaleRecip :: e -> c e -> c e - addConstant :: e -> c e -> c e - add :: c e -> c e -> c e - sub :: c e -> c e -> c e - -- | element by element multiplication - mul :: c e -> c e -> c e - -- | element by element division - divide :: c e -> c e -> c e - equal :: c e -> c e -> Bool +-- -- | Basic element-by-element functions. +-- class (Element e, Container c e) => Linear c e where +-------------------------------------------------------------------------- +instance Container Vector Float where + scale = vectorMapValF Scale + scaleRecip = vectorMapValF Recip + addConstant = vectorMapValF AddConstant + add = vectorZipF Add + sub = vectorZipF Sub + mul = vectorZipF Mul + divide = vectorZipF Div + equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 + scalar x = fromList [x] + -- +--instance Container Vector Float where + cmap = mapVector + atIndex = (@>) + minIndex = round . toScalarF MinIdx + maxIndex = round . toScalarF MaxIdx + minElement = toScalarF Min + maxElement = toScalarF Max + sumElements = sumF + prodElements = prodF + +instance Container Vector Double where + scale = vectorMapValR Scale + scaleRecip = vectorMapValR Recip + addConstant = vectorMapValR AddConstant + add = vectorZipR Add + sub = vectorZipR Sub + mul = vectorZipR Mul + divide = vectorZipR Div + equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 + scalar x = fromList [x] + -- +--instance Container Vector Double where + cmap = mapVector + atIndex = (@>) + minIndex = round . toScalarR MinIdx + maxIndex = round . toScalarR MaxIdx + minElement = toScalarR Min + maxElement = toScalarR Max + sumElements = sumR + prodElements = prodR + +instance Container Vector (Complex Double) where + scale = vectorMapValC Scale + scaleRecip = vectorMapValC Recip + addConstant = vectorMapValC AddConstant + add = vectorZipC Add + sub = vectorZipC Sub + mul = vectorZipC Mul + divide = vectorZipC Div + equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 + scalar x = fromList [x] + -- +--instance Container Vector (Complex Double) where + cmap = mapVector + atIndex = (@>) + minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) + maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) + minElement = ap (@>) minIndex + maxElement = ap (@>) maxIndex + sumElements = sumC + prodElements = prodC + +instance Container Vector (Complex Float) where + scale = vectorMapValQ Scale + scaleRecip = vectorMapValQ Recip + addConstant = vectorMapValQ AddConstant + add = vectorZipQ Add + sub = vectorZipQ Sub + mul = vectorZipQ Mul + divide = vectorZipQ Div + equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 + scalar x = fromList [x] + -- +--instance Container Vector (Complex Float) where + cmap = mapVector + atIndex = (@>) + minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) + maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) + minElement = ap (@>) minIndex + maxElement = ap (@>) maxIndex + sumElements = sumQ + prodElements = prodQ + +--------------------------------------------------------------- + +instance (Container Vector a) => Container Matrix a where + scale x = liftMatrix (scale x) + scaleRecip x = liftMatrix (scaleRecip x) + addConstant x = liftMatrix (addConstant x) + add = liftMatrix2 add + sub = liftMatrix2 sub + mul = liftMatrix2 mul + divide = liftMatrix2 divide + equal a b = cols a == cols b && flatten a `equal` flatten b + scalar x = (1><1) [x] + -- +--instance (Container Vector a) => Container Matrix a where + cmap f = liftMatrix (mapVector f) + atIndex = (@@>) + minIndex m = let (r,c) = (rows m,cols m) + i = (minIndex $ flatten m) + in (i `div` c,(i `mod` c) + 1) + maxIndex m = let (r,c) = (rows m,cols m) + i = (maxIndex $ flatten m) + in (i `div` c,(i `mod` c) + 1) + minElement = ap (@@>) minIndex + maxElement = ap (@@>) maxIndex + sumElements = sumElements . flatten + prodElements = prodElements . flatten + +---------------------------------------------------- + + +-- | Linear algebraic properties of objects +class Num e => Vectors a e where + -- | dot (inner) product + dot :: a e -> a e -> e + -- | sum of absolute value of elements (differs in complex case from @norm1@ + absSum :: a e -> e + -- | sum of absolute value of elements + norm1 :: a e -> e + -- | euclidean norm + norm2 :: a e -> e + -- | element of maximum magnitude + normInf :: a e -> e + +instance Vectors Vector Float where + norm2 = toScalarF Norm2 + absSum = toScalarF AbsSum + dot = dotF + norm1 = toScalarF AbsSum + normInf = maxElement . vectorMapF Abs + +instance Vectors Vector Double where + norm2 = toScalarR Norm2 + absSum = toScalarR AbsSum + dot = dotR + norm1 = toScalarR AbsSum + normInf = maxElement . vectorMapR Abs + +instance Vectors Vector (Complex Float) where + norm2 = (:+ 0) . toScalarQ Norm2 + absSum = (:+ 0) . toScalarQ AbsSum + dot = dotQ + norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapQ Abs + normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapQ Abs + +instance Vectors Vector (Complex Double) where + norm2 = (:+ 0) . toScalarC Norm2 + absSum = (:+ 0) . toScalarC AbsSum + dot = dotC + norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapC Abs + normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapC Abs + +---------------------------------------------------- + +class Element t => Product t where + multiply :: Matrix t -> Matrix t -> Matrix t + ctrans :: Matrix t -> Matrix t + +instance Product Double where + multiply = multiplyR + ctrans = trans + +instance Product (Complex Double) where + multiply = multiplyC + ctrans = conj . trans + +instance Product Float where + multiply = multiplyF + ctrans = trans + +instance Product (Complex Float) where + multiply = multiplyQ + ctrans = conj . trans + +---------------------------------------------------------- + +-- synonym for matrix product +mXm :: Product t => Matrix t -> Matrix t -> Matrix t +mXm = multiply + +-- matrix - vector product +mXv :: Product t => Matrix t -> Vector t -> Vector t +mXv m v = flatten $ m `mXm` (asColumn v) + +-- vector - matrix product +vXm :: Product t => Vector t -> Matrix t -> Vector t +vXm v m = flatten $ (asRow v) `mXm` m + +{- | Outer product of two vectors. + +@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3] +(3><3) + [ 5.0, 2.0, 3.0 + , 10.0, 4.0, 6.0 + , 15.0, 6.0, 9.0 ]@ +-} +outer :: (Product t) => Vector t -> Vector t -> Matrix t +outer u v = asColumn u `multiply` asRow v + +{- | Kronecker product of two matrices. + +@m1=(2><3) + [ 1.0, 2.0, 0.0 + , 0.0, -1.0, 3.0 ] +m2=(4><3) + [ 1.0, 2.0, 3.0 + , 4.0, 5.0, 6.0 + , 7.0, 8.0, 9.0 + , 10.0, 11.0, 12.0 ]@ + +@\> kronecker m1 m2 +(8><9) + [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0 + , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0 + , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0 + , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0 + , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0 + , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0 + , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 + , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ +-} +kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t +kronecker a b = fromBlocks + . splitEvery (cols a) + . map (reshape (cols b)) + . toRows + $ flatten a `outer` flatten b diff --git a/lib/Numeric/LinearAlgebra.hs b/lib/Numeric/LinearAlgebra.hs index da4a4b5..edfd87e 100644 --- a/lib/Numeric/LinearAlgebra.hs +++ b/lib/Numeric/LinearAlgebra.hs @@ -14,14 +14,12 @@ This module reexports all normally required functions for Linear Algebra applica ----------------------------------------------------------------------------- module Numeric.LinearAlgebra ( module Numeric.LinearAlgebra.Algorithms, - module Numeric.LinearAlgebra.Interface, - module Numeric.LinearAlgebra.Linear, +-- module Numeric.LinearAlgebra.Interface, module Numeric.Matrix, - module Numeric.Vector + module Data.Packed.Random ) where import Numeric.LinearAlgebra.Algorithms -import Numeric.LinearAlgebra.Interface -import Numeric.LinearAlgebra.Linear +--import Numeric.LinearAlgebra.Interface import Numeric.Matrix -import Numeric.Vector +import Data.Packed.Random diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index 8306961..fa6c475 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs @@ -77,15 +77,16 @@ module Numeric.LinearAlgebra.Algorithms ( import Data.Packed.Internal hiding ((//)) import Data.Packed.Matrix import Data.Complex -import Numeric.LinearAlgebra.Linear +--import Numeric.LinearAlgebra.Linear import Numeric.LinearAlgebra.LAPACK as LAPACK import Data.List(foldl1') import Data.Array -import Numeric.Vector -import Numeric.Matrix() +import Numeric.Container + +constant x = constantD x -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. -class (Product t, Linear Vector t, Container Vector t, Container Matrix t) => Field t where +class (Product t, Container Vector t, Container Matrix t) => Field t where svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) sv' :: Matrix t -> Vector Double diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs deleted file mode 100644 index 04a9d88..0000000 --- a/lib/Numeric/LinearAlgebra/Instances.hs +++ /dev/null @@ -1,276 +0,0 @@ -{-# LANGUAGE UndecidableInstances, FlexibleInstances #-} ------------------------------------------------------------------------------ -{- | -Module : Numeric.LinearAlgebra.Instances -Copyright : (c) Alberto Ruiz 2006 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) -Stability : provisional -Portability : portable - -This module exports Show, Read, Eq, Num, Fractional, and Floating instances for Vector and Matrix. - -In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand. - --} ------------------------------------------------------------------------------ - -module Numeric.LinearAlgebra.Instances( -) where - -import Numeric.LinearAlgebra.Linear -import Numeric.GSL.Vector -import Data.Packed.Matrix -import Data.Complex -import Data.List(transpose,intersperse) -import Data.Packed.Internal.Vector - -#ifndef VECTOR -import Foreign(Storable) -#endif - ------------------------------------------------------------------- - -instance (Show a, Element a) => (Show (Matrix a)) where - show m = (sizes++) . dsp . map (map show) . toLists $ m - where sizes = "("++show (rows m)++"><"++show (cols m)++")\n" - -dsp as = (++" ]") . (" ["++) . init . drop 2 . unlines . map (" , "++) . map unwords' $ transpose mtp - where - mt = transpose as - longs = map (maximum . map length) mt - mtp = zipWith (\a b -> map (pad a) b) longs mt - pad n str = replicate (n - length str) ' ' ++ str - unwords' = concat . intersperse ", " - -#ifndef VECTOR - -instance (Show a, Storable a) => (Show (Vector a)) where - show v = (show (dim v))++" |> " ++ show (toList v) - -#endif - ------------------------------------------------------------------- - -instance (Element a, Read a) => Read (Matrix a) where - readsPrec _ s = [((rs>' $ dims - -#ifdef VECTOR - -instance (Element a, Read a) => Read (Vector a) where - readsPrec _ s = [(fromList . read $ listnums, rest)] - where (thing,trest) = breakAt ']' s - (dims,listnums) = breakAt ' ' (dropWhile (==' ') thing) - rest = drop 31 trest -#else - -instance (Element a, Read a) => Read (Vector a) where - readsPrec _ s = [((d |>) . read $ listnums, rest)] - where (thing,rest) = breakAt ']' s - (dims,listnums) = breakAt '>' thing - d = read . init . fst . breakAt '|' $ dims - -#endif - -breakAt c l = (a++[c],tail b) where - (a,b) = break (==c) l - ------------------------------------------------------------------- - -adaptScalar f1 f2 f3 x y - | dim x == 1 = f1 (x@>0) y - | dim y == 1 = f3 x (y@>0) - | otherwise = f2 x y - -#ifndef VECTOR - -instance Linear Vector a => Eq (Vector a) where - (==) = equal - -#endif - -instance Num (Vector Float) where - (+) = adaptScalar addConstant add (flip addConstant) - negate = scale (-1) - (*) = adaptScalar scale mul (flip scale) - signum = vectorMapF Sign - abs = vectorMapF Abs - fromInteger = fromList . return . fromInteger - -instance Num (Vector Double) where - (+) = adaptScalar addConstant add (flip addConstant) - negate = scale (-1) - (*) = adaptScalar scale mul (flip scale) - signum = vectorMapR Sign - abs = vectorMapR Abs - fromInteger = fromList . return . fromInteger - -instance Num (Vector (Complex Double)) where - (+) = adaptScalar addConstant add (flip addConstant) - negate = scale (-1) - (*) = adaptScalar scale mul (flip scale) - signum = vectorMapC Sign - abs = vectorMapC Abs - fromInteger = fromList . return . fromInteger - -instance Num (Vector (Complex Float)) where - (+) = adaptScalar addConstant add (flip addConstant) - negate = scale (-1) - (*) = adaptScalar scale mul (flip scale) - signum = vectorMapQ Sign - abs = vectorMapQ Abs - fromInteger = fromList . return . fromInteger - -instance Linear Matrix a => Eq (Matrix a) where - (==) = equal - -instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where - (+) = liftMatrix2Auto (+) - (-) = liftMatrix2Auto (-) - negate = liftMatrix negate - (*) = liftMatrix2Auto (*) - signum = liftMatrix signum - abs = liftMatrix abs - fromInteger = (1><1) . return . fromInteger - ---------------------------------------------------- - -instance (Linear Vector a, Num (Vector a)) => Fractional (Vector a) where - fromRational n = fromList [fromRational n] - (/) = adaptScalar f divide g where - r `f` v = scaleRecip r v - v `g` r = scale (recip r) v - -------------------------------------------------------- - -instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where - fromRational n = (1><1) [fromRational n] - (/) = liftMatrix2Auto (/) - ---------------------------------------------------------- - -instance Floating (Vector Float) where - sin = vectorMapF Sin - cos = vectorMapF Cos - tan = vectorMapF Tan - asin = vectorMapF ASin - acos = vectorMapF ACos - atan = vectorMapF ATan - sinh = vectorMapF Sinh - cosh = vectorMapF Cosh - tanh = vectorMapF Tanh - asinh = vectorMapF ASinh - acosh = vectorMapF ACosh - atanh = vectorMapF ATanh - exp = vectorMapF Exp - log = vectorMapF Log - sqrt = vectorMapF Sqrt - (**) = adaptScalar (vectorMapValF PowSV) (vectorZipF Pow) (flip (vectorMapValF PowVS)) - pi = fromList [pi] - -------------------------------------------------------------- - -instance Floating (Vector Double) where - sin = vectorMapR Sin - cos = vectorMapR Cos - tan = vectorMapR Tan - asin = vectorMapR ASin - acos = vectorMapR ACos - atan = vectorMapR ATan - sinh = vectorMapR Sinh - cosh = vectorMapR Cosh - tanh = vectorMapR Tanh - asinh = vectorMapR ASinh - acosh = vectorMapR ACosh - atanh = vectorMapR ATanh - exp = vectorMapR Exp - log = vectorMapR Log - sqrt = vectorMapR Sqrt - (**) = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS)) - pi = fromList [pi] - -------------------------------------------------------------- - -instance Floating (Vector (Complex Double)) where - sin = vectorMapC Sin - cos = vectorMapC Cos - tan = vectorMapC Tan - asin = vectorMapC ASin - acos = vectorMapC ACos - atan = vectorMapC ATan - sinh = vectorMapC Sinh - cosh = vectorMapC Cosh - tanh = vectorMapC Tanh - asinh = vectorMapC ASinh - acosh = vectorMapC ACosh - atanh = vectorMapC ATanh - exp = vectorMapC Exp - log = vectorMapC Log - sqrt = vectorMapC Sqrt - (**) = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS)) - pi = fromList [pi] - ------------------------------------------------------------ - -instance Floating (Vector (Complex Float)) where - sin = vectorMapQ Sin - cos = vectorMapQ Cos - tan = vectorMapQ Tan - asin = vectorMapQ ASin - acos = vectorMapQ ACos - atan = vectorMapQ ATan - sinh = vectorMapQ Sinh - cosh = vectorMapQ Cosh - tanh = vectorMapQ Tanh - asinh = vectorMapQ ASinh - acosh = vectorMapQ ACosh - atanh = vectorMapQ ATanh - exp = vectorMapQ Exp - log = vectorMapQ Log - sqrt = vectorMapQ Sqrt - (**) = adaptScalar (vectorMapValQ PowSV) (vectorZipQ Pow) (flip (vectorMapValQ PowVS)) - pi = fromList [pi] - ------------------------------------------------------------ - -instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where - sin = liftMatrix sin - cos = liftMatrix cos - tan = liftMatrix tan - asin = liftMatrix asin - acos = liftMatrix acos - atan = liftMatrix atan - sinh = liftMatrix sinh - cosh = liftMatrix cosh - tanh = liftMatrix tanh - asinh = liftMatrix asinh - acosh = liftMatrix acosh - atanh = liftMatrix atanh - exp = liftMatrix exp - log = liftMatrix log - (**) = liftMatrix2Auto (**) - sqrt = liftMatrix sqrt - pi = (1><1) [pi] - ---------------------------------------------------------------- - --- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where --- mempty = 0 { idim = 0 } --- mappend a b = mconcat [a,b] --- mconcat = j . filter ((>0).dim) --- where j [] = mempty --- j l = join l - ---------------------------------------------------------------- - --- instance (NFData a, Storable a) => NFData (Vector a) where --- rnf = rnf . (@>0) --- --- instance (NFData a, Element a) => NFData (Matrix a) where --- rnf = rnf . flatten - diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs deleted file mode 100644 index fa3e209..0000000 --- a/lib/Numeric/LinearAlgebra/Interface.hs +++ /dev/null @@ -1,118 +0,0 @@ -{-# OPTIONS_GHC -fglasgow-exts #-} -{-# LANGUAGE UndecidableInstances #-} ------------------------------------------------------------------------------ -{- | -Module : Numeric.LinearAlgebra.Interface -Copyright : (c) Alberto Ruiz 2007 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) -Stability : provisional -Portability : portable - -Some useful operators, and Show, Read, Eq, Num, Fractional, and Floating instances for Vector and Matrix. - -In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand. - - --} ------------------------------------------------------------------------------ - -module Numeric.LinearAlgebra.Interface( - (<>),(<.>), - (<\>), - (.*),(*/), - (<|>),(<->), -) where - -import Numeric.Vector -import Numeric.Matrix -import Numeric.LinearAlgebra.Algorithms -import Numeric.LinearAlgebra.Linear - -class Mul a b c | a b -> c where - infixl 7 <> - -- | Matrix-matrix, matrix-vector, and vector-matrix products. - (<>) :: Product t => a t -> b t -> c t - -instance Mul Matrix Matrix Matrix where - (<>) = mXm - -instance Mul Matrix Vector Vector where - (<>) m v = flatten $ m <> (asColumn v) - -instance Mul Vector Matrix Vector where - (<>) v m = flatten $ (asRow v) <> m - ---------------------------------------------------- - --- | Dot product: @u \<.\> v = dot u v@ ---(<.>) :: (Field t) => Vector t -> Vector t -> t -(<.>) :: Vectors Vector t => Vector t -> Vector t -> t -infixl 7 <.> -(<.>) = dot - ----------------------------------------------------- - -{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-} - --- -- | @x .* a = scale x a@ --- (.*) :: (Linear c a) => a -> c a -> c a -infixl 7 .* -a .* x = scale a x - ----------------------------------------------------- - -{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-} - --- -- | @a *\/ x = scale (recip x) a@ --- (*/) :: (Linear c a) => c a -> a -> c a -infixl 7 */ -v */ x = scale (recip x) v - --- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD). -(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a -infixl 7 <\> -m <\> v = flatten (linearSolveSVD m (reshape 1 v)) - ------------------------------------------------- - -{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-} -{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-} - -class Joinable a b where - joinH :: Element t => a t -> b t -> Matrix t - joinV :: Element t => a t -> b t -> Matrix t - -instance Joinable Matrix Matrix where - joinH m1 m2 = fromBlocks [[m1,m2]] - joinV m1 m2 = fromBlocks [[m1],[m2]] - -instance Joinable Matrix Vector where - joinH m v = joinH m (asColumn v) - joinV m v = joinV m (asRow v) - -instance Joinable Vector Matrix where - joinH v m = joinH (asColumn v) m - joinV v m = joinV (asRow v) m - -infixl 4 <|> -infixl 3 <-> - -{-- - | Horizontal concatenation of matrices and vectors: - -@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0] -(6><4) - [ 1.0, 0.0, 0.0, 1.0 - , 0.0, 1.0, 0.0, 2.0 - , 0.0, 0.0, 1.0, 3.0 - , 3.0, 0.0, 0.0, 4.0 - , 0.0, 3.0, 0.0, 5.0 - , 0.0, 0.0, 3.0, 6.0 ]@ --} --- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t -a <|> b = joinH a b - --- -- | Vertical concatenation of matrices and vectors. --- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t -a <-> b = joinV a b diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs index 5d0154d..288439f 100644 --- a/lib/Numeric/LinearAlgebra/LAPACK.hs +++ b/lib/Numeric/LinearAlgebra/LAPACK.hs @@ -44,8 +44,7 @@ module Numeric.LinearAlgebra.LAPACK ( import Data.Packed.Internal import Data.Packed.Matrix import Data.Complex -import Numeric.Vector() -import Numeric.Container +import Numeric.Conversion import Numeric.GSL.Vector(vectorMapValR, FunCodeSV(Scale)) import Foreign import Foreign.C.Types (CInt) diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs deleted file mode 100644 index 775060e..0000000 --- a/lib/Numeric/LinearAlgebra/Linear.hs +++ /dev/null @@ -1,160 +0,0 @@ -{-# LANGUAGE UndecidableInstances, MultiParamTypeClasses, FlexibleInstances #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE TypeFamilies #-} ------------------------------------------------------------------------------ -{- | -Module : Numeric.LinearAlgebra.Linear -Copyright : (c) Alberto Ruiz 2006-7 -License : GPL-style - -Maintainer : Alberto Ruiz (aruiz at um dot es) -Stability : provisional -Portability : uses ffi - -Basic optimized operations on vectors and matrices. - --} ------------------------------------------------------------------------------ - -module Numeric.LinearAlgebra.Linear ( - -- * Linear Algebra Typeclasses - Vectors(..), - -- * Products - Product(..), - mXm,mXv,vXm, - outer, kronecker, - -- * Modules - --module Numeric.Vector, - --module Numeric.Matrix, - module Numeric.Container -) where - -import Data.Packed.Internal.Common -import Data.Packed.Matrix -import Data.Packed.Vector -import Data.Complex -import Numeric.Container -import Numeric.Vector() -import Numeric.Matrix() -import Numeric.GSL.Vector -import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) - --- | Linear algebraic properties of objects -class Num e => Vectors a e where - -- | dot (inner) product - dot :: a e -> a e -> e - -- | sum of absolute value of elements (differs in complex case from @norm1@ - absSum :: a e -> e - -- | sum of absolute value of elements - norm1 :: a e -> e - -- | euclidean norm - norm2 :: a e -> e - -- | element of maximum magnitude - normInf :: a e -> e - -instance Vectors Vector Float where - norm2 = toScalarF Norm2 - absSum = toScalarF AbsSum - dot = dotF - norm1 = toScalarF AbsSum - normInf = maxElement . vectorMapF Abs - -instance Vectors Vector Double where - norm2 = toScalarR Norm2 - absSum = toScalarR AbsSum - dot = dotR - norm1 = toScalarR AbsSum - normInf = maxElement . vectorMapR Abs - -instance Vectors Vector (Complex Float) where - norm2 = (:+ 0) . toScalarQ Norm2 - absSum = (:+ 0) . toScalarQ AbsSum - dot = dotQ - norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapQ Abs - normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapQ Abs - -instance Vectors Vector (Complex Double) where - norm2 = (:+ 0) . toScalarC Norm2 - absSum = (:+ 0) . toScalarC AbsSum - dot = dotC - norm1 = (:+ 0) . sumElements . fst . fromComplex . vectorMapC Abs - normInf = (:+ 0) . maxElement . fst . fromComplex . vectorMapC Abs - ----------------------------------------------------- - -class Element t => Product t where - multiply :: Matrix t -> Matrix t -> Matrix t - ctrans :: Matrix t -> Matrix t - -instance Product Double where - multiply = multiplyR - ctrans = trans - -instance Product (Complex Double) where - multiply = multiplyC - ctrans = conj . trans - -instance Product Float where - multiply = multiplyF - ctrans = trans - -instance Product (Complex Float) where - multiply = multiplyQ - ctrans = conj . trans - ----------------------------------------------------------- - --- synonym for matrix product -mXm :: Product t => Matrix t -> Matrix t -> Matrix t -mXm = multiply - --- matrix - vector product -mXv :: Product t => Matrix t -> Vector t -> Vector t -mXv m v = flatten $ m `mXm` (asColumn v) - --- vector - matrix product -vXm :: Product t => Vector t -> Matrix t -> Vector t -vXm v m = flatten $ (asRow v) `mXm` m - -{- | Outer product of two vectors. - -@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3] -(3><3) - [ 5.0, 2.0, 3.0 - , 10.0, 4.0, 6.0 - , 15.0, 6.0, 9.0 ]@ --} -outer :: (Product t) => Vector t -> Vector t -> Matrix t -outer u v = asColumn u `multiply` asRow v - -{- | Kronecker product of two matrices. - -@m1=(2><3) - [ 1.0, 2.0, 0.0 - , 0.0, -1.0, 3.0 ] -m2=(4><3) - [ 1.0, 2.0, 3.0 - , 4.0, 5.0, 6.0 - , 7.0, 8.0, 9.0 - , 10.0, 11.0, 12.0 ]@ - -@\> kronecker m1 m2 -(8><9) - [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0 - , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0 - , 7.0, 8.0, 9.0, 14.0, 16.0, 18.0, 0.0, 0.0, 0.0 - , 10.0, 11.0, 12.0, 20.0, 22.0, 24.0, 0.0, 0.0, 0.0 - , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0 - , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0 - , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 - , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ --} -kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t -kronecker a b = fromBlocks - . splitEvery (cols a) - . map (reshape (cols b)) - . toRows - $ flatten a `outer` flatten b - - -------------------------------------------------------------------- diff --git a/lib/Numeric/Matrix.hs b/lib/Numeric/Matrix.hs index fa3f94a..6ba870f 100644 --- a/lib/Numeric/Matrix.hs +++ b/lib/Numeric/Matrix.hs @@ -3,6 +3,7 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} ----------------------------------------------------------------------------- -- | @@ -14,12 +15,22 @@ -- Stability : provisional -- Portability : portable -- --- Numeric instances and functions for 'Data.Packed.Matrix's +-- Numeric instances and functions for 'Matrix'. +-- In the context of the standard numeric operators, one-component +-- vectors and matrices automatically expand to match the dimensions of the other operand. -- ----------------------------------------------------------------------------- module Numeric.Matrix ( - module Data.Packed.Matrix, + -- * Basic functions + module Data.Packed.Matrix, + module Numeric.Vector, + --module Numeric.Container, + -- * Operators + (<>), (<\>), + -- * Deprecated + (.*),(*/),(<|>),(<->), + vectorMax,vectorMin ) where ------------------------------------------------------------------- @@ -28,18 +39,18 @@ import Data.Packed.Vector import Data.Packed.Matrix import Numeric.Container --import Numeric.LinearAlgebra.Linear -import Numeric.Vector() +import Numeric.Vector +import Numeric.LinearAlgebra.Algorithms +--import Control.Monad(ap) -import Control.Monad(ap) - -import Control.Arrow((***)) +--import Control.Arrow((***)) ------------------------------------------------------------------- -instance Linear Matrix a => Eq (Matrix a) where +instance Container Matrix a => Eq (Matrix a) where (==) = equal -instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where +instance (Container Matrix a, Num (Vector a)) => Num (Matrix a) where (+) = liftMatrix2Auto (+) (-) = liftMatrix2Auto (-) negate = liftMatrix negate @@ -50,13 +61,13 @@ instance (Linear Matrix a, Num (Vector a)) => Num (Matrix a) where --------------------------------------------------- -instance (Linear Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where +instance (Container Vector a, Fractional (Vector a), Num (Matrix a)) => Fractional (Matrix a) where fromRational n = (1><1) [fromRational n] (/) = liftMatrix2Auto (/) --------------------------------------------------------- -instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where +instance (Container Vector a, Floating (Vector a), Fractional (Matrix a)) => Floating (Matrix a) where sin = liftMatrix sin cos = liftMatrix cos tan = liftMatrix tan @@ -75,43 +86,93 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati sqrt = liftMatrix sqrt pi = (1><1) [pi] ---------------------------------------------------------------- - -instance NumericContainer Matrix where - toComplex = uncurry $ liftMatrix2 $ curry toComplex - fromComplex z = (reshape c *** reshape c) . fromComplex . flatten $ z - where c = cols z - complex' = liftMatrix complex' - conj = liftMatrix conj --- cmap f = liftMatrix (cmap f) - single' = liftMatrix single' - double' = liftMatrix double' - ---------------------------------------------------------------- - -instance (Linear Vector a, Container Matrix a) => Linear Matrix a where - scale x = liftMatrix (scale x) - scaleRecip x = liftMatrix (scaleRecip x) - addConstant x = liftMatrix (addConstant x) - add = liftMatrix2 add - sub = liftMatrix2 sub - mul = liftMatrix2 mul - divide = liftMatrix2 divide - equal a b = cols a == cols b && flatten a `equal` flatten b - scalar x = (1><1) [x] - -- -instance (Container Vector a) => Container Matrix a where - cmap f = liftMatrix (mapVector f) - atIndex = (@@>) - minIndex m = let (r,c) = (rows m,cols m) - i = (minIndex $ flatten m) - in (i `div` c,(i `mod` c) + 1) - maxIndex m = let (r,c) = (rows m,cols m) - i = (maxIndex $ flatten m) - in (i `div` c,(i `mod` c) + 1) - minElement = ap (@@>) minIndex - maxElement = ap (@@>) maxIndex - sumElements = sumElements . flatten - prodElements = prodElements . flatten +-------------------------------------------------------- + +class Mul a b c | a b -> c where + infixl 7 <> + -- | Matrix-matrix, matrix-vector, and vector-matrix products. + (<>) :: Product t => a t -> b t -> c t + +instance Mul Matrix Matrix Matrix where + (<>) = mXm + +instance Mul Matrix Vector Vector where + (<>) m v = flatten $ m <> (asColumn v) + +instance Mul Vector Matrix Vector where + (<>) v m = flatten $ (asRow v) <> m ---------------------------------------------------- + +{-# DEPRECATED (.*) "use scale a x or scalar a * x" #-} + +-- -- | @x .* a = scale x a@ +-- (.*) :: (Linear c a) => a -> c a -> c a +infixl 7 .* +a .* x = scale a x + +---------------------------------------------------- + +{-# DEPRECATED (*/) "use scale (recip a) x or x / scalar a" #-} + +-- -- | @a *\/ x = scale (recip x) a@ +-- (*/) :: (Linear c a) => c a -> a -> c a +infixl 7 */ +v */ x = scale (recip x) v + +-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD). +(<\>) :: (Field a) => Matrix a -> Vector a -> Vector a +infixl 7 <\> +m <\> v = flatten (linearSolveSVD m (reshape 1 v)) + +------------------------------------------------ + +{-# DEPRECATED (<|>) "define operator a & b = fromBlocks[[a,b]] and use asRow/asColumn to join vectors" #-} +{-# DEPRECATED (<->) "define operator a // b = fromBlocks[[a],[b]] and use asRow/asColumn to join vectors" #-} + +class Joinable a b where + joinH :: Element t => a t -> b t -> Matrix t + joinV :: Element t => a t -> b t -> Matrix t + +instance Joinable Matrix Matrix where + joinH m1 m2 = fromBlocks [[m1,m2]] + joinV m1 m2 = fromBlocks [[m1],[m2]] + +instance Joinable Matrix Vector where + joinH m v = joinH m (asColumn v) + joinV m v = joinV m (asRow v) + +instance Joinable Vector Matrix where + joinH v m = joinH (asColumn v) m + joinV v m = joinV (asRow v) m + +infixl 4 <|> +infixl 3 <-> + +{-- - | Horizontal concatenation of matrices and vectors: + +@> (ident 3 \<-\> 3 * ident 3) \<|\> fromList [1..6.0] +(6><4) + [ 1.0, 0.0, 0.0, 1.0 + , 0.0, 1.0, 0.0, 2.0 + , 0.0, 0.0, 1.0, 3.0 + , 3.0, 0.0, 0.0, 4.0 + , 0.0, 3.0, 0.0, 5.0 + , 0.0, 0.0, 3.0, 6.0 ]@ +-} +-- (<|>) :: (Element t, Joinable a b) => a t -> b t -> Matrix t +a <|> b = joinH a b + +-- -- | Vertical concatenation of matrices and vectors. +-- (<->) :: (Element t, Joinable a b) => a t -> b t -> Matrix t +a <-> b = joinV a b + +------------------------------------------------------------------- + +{-# DEPRECATED vectorMin "use minElement" #-} +vectorMin :: (Container Vector t, Element t) => Vector t -> t +vectorMin = minElement + +{-# DEPRECATED vectorMax "use maxElement" #-} +vectorMax :: (Container Vector t, Element t) => Vector t -> t +vectorMax = maxElement diff --git a/lib/Numeric/Vector.hs b/lib/Numeric/Vector.hs index 55d645a..9427243 100644 --- a/lib/Numeric/Vector.hs +++ b/lib/Numeric/Vector.hs @@ -14,24 +14,26 @@ -- Stability : provisional -- Portability : portable -- --- Numeric instances and functions for 'Data.Packed.Vector's +-- Numeric instances and functions for 'Vector'. -- ----------------------------------------------------------------------------- module Numeric.Vector ( - -- * Vector creation - constant, linspace, - module Data.Packed.Vector + -- * Basic vector functions + module Data.Packed.Vector, + module Numeric.Container, + -- * Vector creation + constant, linspace, + -- * Operators + (<.>) ) where import Data.Complex -import Control.Monad(ap) - import Data.Packed.Vector -import Data.Packed.Internal.Matrix(Element(..)) -import Data.Packed.Internal.Vector(asComplex,asReal) -import Data.Packed.Matrix(toColumns,fromColumns,flatten,reshape) +import Data.Packed.Internal.Matrix +--import Data.Packed.Internal.Vector(asComplex,asReal) +--import Data.Packed.Matrix(toColumns,fromColumns,flatten,reshape) import Numeric.GSL.Vector import Numeric.Container @@ -92,10 +94,15 @@ Logarithmic spacing can be defined as follows: @logspace n (a,b) = 10 ** linspace n (a,b)@ -} -linspace :: (Enum e, Linear Vector e) => Int -> (e, e) -> Vector e +linspace :: (Enum e, Container Vector e) => Int -> (e, e) -> Vector e linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1] where s = (b-a)/fromIntegral (n-1) +-- | Dot product: @u \<.\> v = dot u v@ +(<.>) :: Vectors Vector t => Vector t -> Vector t -> t +infixl 7 <.> +(<.>) = dot + ------------------------------------------------------------------ adaptScalar f1 f2 f3 x y @@ -107,7 +114,7 @@ adaptScalar f1 f2 f3 x y #ifndef VECTOR -instance Linear Vector a => Eq (Vector a) where +instance Container Vector a => Eq (Vector a) where (==) = equal #endif @@ -146,7 +153,7 @@ instance Num (Vector (Complex Float)) where --------------------------------------------------- -instance (Linear Vector a, Num (Vector a)) => Fractional (Vector a) where +instance (Container Vector a, Num (Vector a)) => Fractional (Vector a) where fromRational n = fromList [fromRational n] (/) = adaptScalar f divide g where r `f` v = scaleRecip r v @@ -254,116 +261,6 @@ instance Floating (Vector (Complex Float)) where -- instance (NFData a, Element a) => NFData (Matrix a) where -- rnf = rnf . flatten ---------------------------------------------------------------- - --- | obtains the complex conjugate of a complex vector -conjV :: (RealElement a) => Vector (Complex a) -> Vector (Complex a) -conjV = mapVector conjugate - --- | creates a complex vector from vectors with real and imaginary parts -toComplexV :: (RealElement a) => (Vector a, Vector a) -> Vector (Complex a) -toComplexV (r,i) = asComplex $ flatten $ fromColumns [r,i] - --- | the inverse of 'toComplex' -fromComplexV :: (RealElement a) => Vector (Complex a) -> (Vector a, Vector a) -fromComplexV z = (r,i) where - [r,i] = toColumns $ reshape 2 $ asReal z -------------------------------------------------------------------------- -instance NumericContainer Vector where - toComplex = toComplexV - fromComplex = fromComplexV - complex' v = toComplex (v,constant 0 (dim v)) - conj = conjV --- cmap = mapVector - single' = double2FloatG - double' = float2DoubleG - --------------------------------------------------------------------------- - -instance Linear Vector Float where - scale = vectorMapValF Scale - scaleRecip = vectorMapValF Recip - addConstant = vectorMapValF AddConstant - add = vectorZipF Add - sub = vectorZipF Sub - mul = vectorZipF Mul - divide = vectorZipF Div - equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 - scalar x = fromList [x] - -- -instance Container Vector Float where - cmap = mapVector - atIndex = (@>) - minIndex = round . toScalarF MinIdx - maxIndex = round . toScalarF MaxIdx - minElement = toScalarF Min - maxElement = toScalarF Max - sumElements = sumF - prodElements = prodF - -instance Linear Vector Double where - scale = vectorMapValR Scale - scaleRecip = vectorMapValR Recip - addConstant = vectorMapValR AddConstant - add = vectorZipR Add - sub = vectorZipR Sub - mul = vectorZipR Mul - divide = vectorZipR Div - equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 - scalar x = fromList [x] - -- -instance Container Vector Double where - cmap = mapVector - atIndex = (@>) - minIndex = round . toScalarR MinIdx - maxIndex = round . toScalarR MaxIdx - minElement = toScalarR Min - maxElement = toScalarR Max - sumElements = sumR - prodElements = prodR - -instance Linear Vector (Complex Double) where - scale = vectorMapValC Scale - scaleRecip = vectorMapValC Recip - addConstant = vectorMapValC AddConstant - add = vectorZipC Add - sub = vectorZipC Sub - mul = vectorZipC Mul - divide = vectorZipC Div - equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 - scalar x = fromList [x] - -- -instance Container Vector (Complex Double) where - cmap = mapVector - atIndex = (@>) - minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) - maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) - minElement = ap (@>) minIndex - maxElement = ap (@>) maxIndex - sumElements = sumC - prodElements = prodC - -instance Linear Vector (Complex Float) where - scale = vectorMapValQ Scale - scaleRecip = vectorMapValQ Recip - addConstant = vectorMapValQ AddConstant - add = vectorZipQ Add - sub = vectorZipQ Sub - mul = vectorZipQ Mul - divide = vectorZipQ Div - equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 - scalar x = fromList [x] - -- -instance Container Vector (Complex Float) where - cmap = mapVector - atIndex = (@>) - minIndex = minIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) - maxIndex = maxIndex . fst . fromComplex . (zipVectorWith (*) `ap` mapVector conjugate) - minElement = ap (@>) minIndex - maxElement = ap (@>) maxIndex - sumElements = sumQ - prodElements = prodQ - ---------------------------------------------------------------- -- cgit v1.2.3