From fd94ecb3c3032beccdca4a4dee38bb306f57cd8b Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 16 May 2014 20:57:13 +0200 Subject: Numeric.Container compatible --- packages/base/src/Data/Packed/Internal/Numeric.hs | 607 ++++++++++++++++++ packages/base/src/Data/Packed/Numeric.hs | 680 +++++---------------- packages/base/src/Numeric/Chain.hs | 2 +- packages/base/src/Numeric/Container.hs | 240 -------- .../base/src/Numeric/LinearAlgebra/Algorithms.hs | 2 +- packages/base/src/Numeric/LinearAlgebra/Base.hs | 2 +- packages/base/src/Numeric/LinearAlgebra/Data.hs | 3 +- packages/base/src/Numeric/LinearAlgebra/Devel.hs | 2 +- packages/base/src/Numeric/LinearAlgebra/Util.hs | 14 +- .../src/Numeric/LinearAlgebra/Util/Convolution.hs | 2 +- packages/base/src/Numeric/Matrix.hs | 2 +- packages/base/src/Numeric/Vector.hs | 2 +- 12 files changed, 775 insertions(+), 783 deletions(-) create mode 100644 packages/base/src/Data/Packed/Internal/Numeric.hs delete mode 100644 packages/base/src/Numeric/Container.hs (limited to 'packages/base/src') diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs new file mode 100644 index 0000000..81a8083 --- /dev/null +++ b/packages/base/src/Data/Packed/Internal/Numeric.hs @@ -0,0 +1,607 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Data.Packed.Internal.Numeric +-- Copyright : (c) Alberto Ruiz 2010-14 +-- License : BSD3 +-- Maintainer : Alberto Ruiz +-- Stability : provisional +-- +----------------------------------------------------------------------------- + +module Data.Packed.Internal.Numeric ( + -- * Basic functions + ident, diag, ctrans, + -- * Generic operations + Container(..), + -- * Matrix product and related functions + Product(..), udot, + mXm,mXv,vXm, + outer, kronecker, + -- * Element conversion + Convert(..), + Complexable(), + RealElement(), + + RealOf, ComplexOf, SingleOf, DoubleOf, + + IndexOf, + module Data.Complex +) where + +import Data.Packed +import Data.Packed.ST as ST +import Numeric.Conversion +import Data.Packed.Development +import Numeric.Vectorized +import Data.Complex +import Control.Applicative((<*>)) + +import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) + +------------------------------------------------------------------- + +type family IndexOf (c :: * -> *) + +type instance IndexOf Vector = Int +type instance IndexOf Matrix = (Int,Int) + +type family ArgOf (c :: * -> *) a + +type instance ArgOf Vector a = a -> a +type instance ArgOf Matrix a = a -> a -> a + +------------------------------------------------------------------- + +-- | Basic element-by-element functions for numeric containers +class (Complexable c, Fractional e, Element e) => Container c e where + -- | create a structure with a single element + -- + -- >>> let v = fromList [1..3::Double] + -- >>> v / scalar (norm2 v) + -- fromList [0.2672612419124244,0.5345224838248488,0.8017837257372732] + -- + scalar :: e -> c e + -- | complex conjugate + conj :: c 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 + -- + -- element by element inverse tangent + arctan2 :: c e -> c e -> c e + -- + -- | cannot implement instance Functor because of Element class constraint + cmap :: (Element b) => (e -> b) -> c e -> c b + -- | constant structure of given size + konst' :: e -> IndexOf c -> c e + -- | create a structure using a function + -- + -- Hilbert matrix of order N: + -- + -- @hilb n = build' (n,n) (\\i j -> 1/(i+j+1))@ + build' :: IndexOf c -> (ArgOf c e) -> c e + -- | indexing function + atIndex :: c e -> IndexOf c -> e + -- | index of min element + minIndex :: c e -> IndexOf c + -- | index of max element + maxIndex :: c e -> IndexOf c + -- | value of min element + minElement :: c e -> e + -- | value of max element + maxElement :: c e -> e + -- the C functions sumX/prodX are twice as fast as using foldVector + -- | the sum of elements (faster than using @fold@) + sumElements :: c e -> e + -- | the product of elements (faster than using @fold@) + prodElements :: c e -> e + + -- | A more efficient implementation of @cmap (\\x -> if x>0 then 1 else 0)@ + -- + -- >>> step $ linspace 5 (-1,1::Double) + -- 5 |> [0.0,0.0,0.0,1.0,1.0] + -- + + step :: RealElement e => c e -> c e + + -- | Element by element version of @case compare a b of {LT -> l; EQ -> e; GT -> g}@. + -- + -- Arguments with any dimension = 1 are automatically expanded: + -- + -- >>> cond ((1><4)[1..]) ((3><1)[1..]) 0 100 ((3><4)[1..]) :: Matrix Double + -- (3><4) + -- [ 100.0, 2.0, 3.0, 4.0 + -- , 0.0, 100.0, 7.0, 8.0 + -- , 0.0, 0.0, 100.0, 12.0 ] + -- + + cond :: RealElement e + => c e -- ^ a + -> c e -- ^ b + -> c e -- ^ l + -> c e -- ^ e + -> c e -- ^ g + -> c e -- ^ result + + -- | Find index of elements which satisfy a predicate + -- + -- >>> find (>0) (ident 3 :: Matrix Double) + -- [(0,0),(1,1),(2,2)] + -- + + find :: (e -> Bool) -> c e -> [IndexOf c] + + -- | Create a structure from an association list + -- + -- >>> assoc 5 0 [(3,7),(1,4)] :: Vector Double + -- fromList [0.0,4.0,0.0,7.0,0.0] + -- + -- >>> assoc (2,3) 0 [((0,2),7),((1,0),2*i-3)] :: Matrix (Complex Double) + -- (2><3) + -- [ 0.0 :+ 0.0, 0.0 :+ 0.0, 7.0 :+ 0.0 + -- , (-3.0) :+ 2.0, 0.0 :+ 0.0, 0.0 :+ 0.0 ] + -- + assoc :: IndexOf c -- ^ size + -> e -- ^ default value + -> [(IndexOf c, e)] -- ^ association list + -> c e -- ^ result + + -- | Modify a structure using an update function + -- + -- >>> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double + -- (5><5) + -- [ 1.0, 0.0, 0.0, 3.0, 0.0 + -- , 0.0, 6.0, 0.0, 0.0, 0.0 + -- , 0.0, 0.0, 1.0, 0.0, 0.0 + -- , 0.0, 0.0, 0.0, 1.0, 0.0 + -- , 0.0, 0.0, 0.0, 0.0, 1.0 ] + -- + -- computation of histogram: + -- + -- >>> accum (konst 0 7) (+) (map (flip (,) 1) [4,5,4,1,5,2,5]) :: Vector Double + -- fromList [0.0,1.0,1.0,0.0,2.0,3.0,0.0] + -- + + accum :: c e -- ^ initial structure + -> (e -> e -> e) -- ^ update function + -> [(IndexOf c, e)] -- ^ association list + -> c e -- ^ result + +-------------------------------------------------------------------------- + +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 + arctan2 = vectorZipF ATan2 + scalar x = fromList [x] + konst' = constant + build' = buildV + conj = id + cmap = mapVector + atIndex = (@>) + minIndex = emptyErrorV "minIndex" (round . toScalarF MinIdx) + maxIndex = emptyErrorV "maxIndex" (round . toScalarF MaxIdx) + minElement = emptyErrorV "minElement" (toScalarF Min) + maxElement = emptyErrorV "maxElement" (toScalarF Max) + sumElements = sumF + prodElements = prodF + step = stepF + find = findV + assoc = assocV + accum = accumV + cond = condV condF + +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 + arctan2 = vectorZipR ATan2 + scalar x = fromList [x] + konst' = constant + build' = buildV + conj = id + cmap = mapVector + atIndex = (@>) + minIndex = emptyErrorV "minIndex" (round . toScalarR MinIdx) + maxIndex = emptyErrorV "maxIndex" (round . toScalarR MaxIdx) + minElement = emptyErrorV "minElement" (toScalarR Min) + maxElement = emptyErrorV "maxElement" (toScalarR Max) + sumElements = sumR + prodElements = prodR + step = stepD + find = findV + assoc = assocV + accum = accumV + cond = condV condD + +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 + arctan2 = vectorZipC ATan2 + scalar x = fromList [x] + konst' = constant + build' = buildV + conj = conjugateC + cmap = mapVector + atIndex = (@>) + minIndex = emptyErrorV "minIndex" (minIndex . fst . fromComplex . (mul <*> conj)) + maxIndex = emptyErrorV "maxIndex" (maxIndex . fst . fromComplex . (mul <*> conj)) + minElement = emptyErrorV "minElement" (atIndex <*> minIndex) + maxElement = emptyErrorV "maxElement" (atIndex <*> maxIndex) + sumElements = sumC + prodElements = prodC + step = undefined -- cannot match + find = findV + assoc = assocV + accum = accumV + cond = undefined -- cannot match + +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 + arctan2 = vectorZipQ ATan2 + scalar x = fromList [x] + konst' = constant + build' = buildV + conj = conjugateQ + cmap = mapVector + atIndex = (@>) + minIndex = emptyErrorV "minIndex" (minIndex . fst . fromComplex . (mul <*> conj)) + maxIndex = emptyErrorV "maxIndex" (maxIndex . fst . fromComplex . (mul <*> conj)) + minElement = emptyErrorV "minElement" (atIndex <*> minIndex) + maxElement = emptyErrorV "maxElement" (atIndex <*> maxIndex) + sumElements = sumQ + prodElements = prodQ + step = undefined -- cannot match + find = findV + assoc = assocV + accum = accumV + cond = undefined -- cannot match + +--------------------------------------------------------------- + +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 + arctan2 = liftMatrix2 arctan2 + scalar x = (1><1) [x] + konst' v (r,c) = matrixFromVector RowMajor r c (konst' v (r*c)) + build' = buildM + conj = liftMatrix conj + cmap f = liftMatrix (mapVector f) + atIndex = (@@>) + minIndex = emptyErrorM "minIndex of Matrix" $ + \m -> divMod (minIndex $ flatten m) (cols m) + maxIndex = emptyErrorM "maxIndex of Matrix" $ + \m -> divMod (maxIndex $ flatten m) (cols m) + minElement = emptyErrorM "minElement of Matrix" (atIndex <*> minIndex) + maxElement = emptyErrorM "maxElement of Matrix" (atIndex <*> maxIndex) + sumElements = sumElements . flatten + prodElements = prodElements . flatten + step = liftMatrix step + find = findM + assoc = assocM + accum = accumM + cond = condM + + +emptyErrorV msg f v = + if dim v > 0 + then f v + else error $ msg ++ " of Vector with dim = 0" + +emptyErrorM msg f m = + if rows m > 0 && cols m > 0 + then f m + else error $ msg++" "++shSize m + +---------------------------------------------------- + +-- | Matrix product and related functions +class (Num e, Element e) => Product e where + -- | matrix product + multiply :: Matrix e -> Matrix e -> Matrix e + -- | sum of absolute value of elements (differs in complex case from @norm1@) + absSum :: Vector e -> RealOf e + -- | sum of absolute value of elements + norm1 :: Vector e -> RealOf e + -- | euclidean norm + norm2 :: Vector e -> RealOf e + -- | element of maximum magnitude + normInf :: Vector e -> RealOf e + +instance Product Float where + norm2 = emptyVal (toScalarF Norm2) + absSum = emptyVal (toScalarF AbsSum) + norm1 = emptyVal (toScalarF AbsSum) + normInf = emptyVal (maxElement . vectorMapF Abs) + multiply = emptyMul multiplyF + +instance Product Double where + norm2 = emptyVal (toScalarR Norm2) + absSum = emptyVal (toScalarR AbsSum) + norm1 = emptyVal (toScalarR AbsSum) + normInf = emptyVal (maxElement . vectorMapR Abs) + multiply = emptyMul multiplyR + +instance Product (Complex Float) where + norm2 = emptyVal (toScalarQ Norm2) + absSum = emptyVal (toScalarQ AbsSum) + norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapQ Abs) + normInf = emptyVal (maxElement . fst . fromComplex . vectorMapQ Abs) + multiply = emptyMul multiplyQ + +instance Product (Complex Double) where + norm2 = emptyVal (toScalarC Norm2) + absSum = emptyVal (toScalarC AbsSum) + norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapC Abs) + normInf = emptyVal (maxElement . fst . fromComplex . vectorMapC Abs) + multiply = emptyMul multiplyC + +emptyMul m a b + | x1 == 0 && x2 == 0 || r == 0 || c == 0 = konst' 0 (r,c) + | otherwise = m a b + where + r = rows a + x1 = cols a + x2 = rows b + c = cols b + +emptyVal f v = + if dim v > 0 + then f v + else 0 + +-- FIXME remove unused C wrappers +-- | unconjugated dot product +udot :: Product e => Vector e -> Vector e -> e +udot u v + | dim u == dim v = val (asRow u `multiply` asColumn v) + | otherwise = error $ "different dimensions "++show (dim u)++" and "++show (dim v)++" in dot product" + where + val m | dim u > 0 = m@@>(0,0) + | otherwise = 0 + +---------------------------------------------------------- + +-- 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 + +------------------------------------------------------------------- + + +class Convert t where + real :: Container c t => c (RealOf t) -> c t + complex :: Container c t => c t -> c (ComplexOf t) + single :: Container c t => c t -> c (SingleOf t) + double :: Container c t => c t -> c (DoubleOf t) + toComplex :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t) + fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t) + + +instance Convert Double where + real = id + complex = comp' + single = single' + double = id + toComplex = toComplex' + fromComplex = fromComplex' + +instance Convert Float where + real = id + complex = comp' + single = id + double = double' + toComplex = toComplex' + fromComplex = fromComplex' + +instance Convert (Complex Double) where + real = comp' + complex = id + single = single' + double = id + toComplex = toComplex' + fromComplex = fromComplex' + +instance Convert (Complex Float) where + real = comp' + complex = id + single = id + double = double' + toComplex = toComplex' + fromComplex = fromComplex' + +------------------------------------------------------------------- + +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 + +------------------------------------------------------------ + +buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ] + where rs = map fromIntegral [0 .. (rc-1)] + cs = map fromIntegral [0 .. (cc-1)] + +buildV n f = fromList [f k | k <- ks] + where ks = map fromIntegral [0 .. (n-1)] + +-------------------------------------------------------- +-- | conjugate transpose +ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e +ctrans = liftMatrix conj . trans + +-- | Creates a square matrix with a given diagonal. +diag :: (Num a, Element a) => Vector a -> Matrix a +diag v = diagRect 0 v n n where n = dim v + +-- | creates the identity matrix of given dimension +ident :: (Num a, Element a) => Int -> Matrix a +ident n = diag (constant 1 n) + +-------------------------------------------------------- + +findV p x = foldVectorWithIndex g [] x where + g k z l = if p z then k:l else l + +findM p x = map ((`divMod` cols x)) $ findV p (flatten x) + +assocV n z xs = ST.runSTVector $ do + v <- ST.newVector z n + mapM_ (\(k,x) -> ST.writeVector v k x) xs + return v + +assocM (r,c) z xs = ST.runSTMatrix $ do + m <- ST.newMatrix z r c + mapM_ (\((i,j),x) -> ST.writeMatrix m i j x) xs + return m + +accumV v0 f xs = ST.runSTVector $ do + v <- ST.thawVector v0 + mapM_ (\(k,x) -> ST.modifyVector v k (f x)) xs + return v + +accumM m0 f xs = ST.runSTMatrix $ do + m <- ST.thawMatrix m0 + mapM_ (\((i,j),x) -> ST.modifyMatrix m i j (f x)) xs + return m + +---------------------------------------------------------------------- + +condM a b l e t = matrixFromVector RowMajor (rows a'') (cols a'') $ cond a' b' l' e' t' + where + args@(a'':_) = conformMs [a,b,l,e,t] + [a', b', l', e', t'] = map flatten args + +condV f a b l e t = f a' b' l' e' t' + where + [a', b', l', e', t'] = conformVs [a,b,l,e,t] + diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs index c13e91d..6036e8c 100644 --- a/packages/base/src/Data/Packed/Numeric.hs +++ b/packages/base/src/Data/Packed/Numeric.hs @@ -1,8 +1,8 @@ -{-# LANGUAGE CPP #-} {-# LANGUAGE TypeFamilies #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE UndecidableInstances #-} ----------------------------------------------------------------------------- @@ -13,16 +13,32 @@ -- Maintainer : Alberto Ruiz -- Stability : provisional -- +-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines. +-- +-- The 'Container' class is used to define optimized generic functions which work +-- on 'Vector' and 'Matrix' with real or complex elements. +-- +-- Some of these functions are also available in the instances of the standard +-- numeric Haskell classes provided by "Numeric.LinearAlgebra". +-- ----------------------------------------------------------------------------- +{-# OPTIONS_HADDOCK hide #-} module Data.Packed.Numeric ( -- * Basic functions - ident, diag, ctrans, + module Data.Packed, + konst, build, + linspace, + diag, ident, + ctrans, -- * Generic operations Container(..), - -- * Matrix product and related functions - Product(..), udot, - mXm,mXv,vXm, + -- * Matrix product + Product(..), udot, dot, (◇), + Mul(..), + Contraction(..), + optimiseMult, + mXm,mXv,vXm,LSDiv(..), outer, kronecker, -- * Element conversion Convert(..), @@ -32,576 +48,196 @@ module Data.Packed.Numeric ( RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf, - module Data.Complex + module Data.Complex, + -- * IO + module Data.Packed.IO ) where -import Data.Packed -import Data.Packed.ST as ST -import Numeric.Conversion -import Data.Packed.Development -import Numeric.Vectorized +import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ) +import Data.Packed.Internal.Numeric import Data.Complex -import Control.Applicative((<*>)) - -import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) +import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD) +import Data.Monoid(Monoid(mconcat)) +import Data.Packed.IO -------------------------------------------------------------------- +------------------------------------------------------------------ -type family IndexOf (c :: * -> *) +{- | Creates a real vector containing a range of values: -type instance IndexOf Vector = Int -type instance IndexOf Matrix = (Int,Int) +>>> linspace 5 (-3,7::Double) +fromList [-3.0,-0.5,2.0,4.5,7.0]@ -type family ArgOf (c :: * -> *) a +>>> linspace 5 (8,2+i) :: Vector (Complex Double) +fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0] -type instance ArgOf Vector a = a -> a -type instance ArgOf Matrix a = a -> a -> a +Logarithmic spacing can be defined as follows: -------------------------------------------------------------------- - --- | Basic element-by-element functions for numeric containers -class (Complexable c, Fractional e, Element e) => Container c e where - -- | create a structure with a single element - -- - -- >>> let v = fromList [1..3::Double] - -- >>> v / scalar (norm2 v) - -- fromList [0.2672612419124244,0.5345224838248488,0.8017837257372732] - -- - scalar :: e -> c e - -- | complex conjugate - conj :: c 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 - -- - -- element by element inverse tangent - arctan2 :: c e -> c e -> c e - -- - -- | cannot implement instance Functor because of Element class constraint - cmap :: (Element b) => (e -> b) -> c e -> c b - -- | constant structure of given size - konst' :: e -> IndexOf c -> c e - -- | create a structure using a function - -- - -- Hilbert matrix of order N: - -- - -- @hilb n = build' (n,n) (\\i j -> 1/(i+j+1))@ - build' :: IndexOf c -> (ArgOf c e) -> c e - -- | indexing function - atIndex :: c e -> IndexOf c -> e - -- | index of min element - minIndex :: c e -> IndexOf c - -- | index of max element - maxIndex :: c e -> IndexOf c - -- | value of min element - minElement :: c e -> e - -- | value of max element - maxElement :: c e -> e - -- the C functions sumX/prodX are twice as fast as using foldVector - -- | the sum of elements (faster than using @fold@) - sumElements :: c e -> e - -- | the product of elements (faster than using @fold@) - prodElements :: c e -> e - - -- | A more efficient implementation of @cmap (\\x -> if x>0 then 1 else 0)@ - -- - -- >>> step $ linspace 5 (-1,1::Double) - -- 5 |> [0.0,0.0,0.0,1.0,1.0] - -- - - step :: RealElement e => c e -> c e - - -- | Element by element version of @case compare a b of {LT -> l; EQ -> e; GT -> g}@. - -- - -- Arguments with any dimension = 1 are automatically expanded: - -- - -- >>> cond ((1><4)[1..]) ((3><1)[1..]) 0 100 ((3><4)[1..]) :: Matrix Double - -- (3><4) - -- [ 100.0, 2.0, 3.0, 4.0 - -- , 0.0, 100.0, 7.0, 8.0 - -- , 0.0, 0.0, 100.0, 12.0 ] - -- - - cond :: RealElement e - => c e -- ^ a - -> c e -- ^ b - -> c e -- ^ l - -> c e -- ^ e - -> c e -- ^ g - -> c e -- ^ result - - -- | Find index of elements which satisfy a predicate - -- - -- >>> find (>0) (ident 3 :: Matrix Double) - -- [(0,0),(1,1),(2,2)] - -- - - find :: (e -> Bool) -> c e -> [IndexOf c] +@logspace n (a,b) = 10 ** linspace n (a,b)@ +-} +linspace :: (Container Vector e) => Int -> (e, e) -> Vector e +linspace 0 (a,b) = fromList[(a+b)/2] +linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1] + where s = (b-a)/fromIntegral (n-1) - -- | Create a structure from an association list - -- - -- >>> assoc 5 0 [(3,7),(1,4)] :: Vector Double - -- fromList [0.0,4.0,0.0,7.0,0.0] - -- - -- >>> assoc (2,3) 0 [((0,2),7),((1,0),2*i-3)] :: Matrix (Complex Double) - -- (2><3) - -- [ 0.0 :+ 0.0, 0.0 :+ 0.0, 7.0 :+ 0.0 - -- , (-3.0) :+ 2.0, 0.0 :+ 0.0, 0.0 :+ 0.0 ] - -- - assoc :: IndexOf c -- ^ size - -> e -- ^ default value - -> [(IndexOf c, e)] -- ^ association list - -> c e -- ^ result +-------------------------------------------------------- - -- | Modify a structure using an update function - -- - -- >>> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double - -- (5><5) - -- [ 1.0, 0.0, 0.0, 3.0, 0.0 - -- , 0.0, 6.0, 0.0, 0.0, 0.0 - -- , 0.0, 0.0, 1.0, 0.0, 0.0 - -- , 0.0, 0.0, 0.0, 1.0, 0.0 - -- , 0.0, 0.0, 0.0, 0.0, 1.0 ] - -- - -- computation of histogram: - -- - -- >>> accum (konst 0 7) (+) (map (flip (,) 1) [4,5,4,1,5,2,5]) :: Vector Double - -- fromList [0.0,1.0,1.0,0.0,2.0,3.0,0.0] - -- - - accum :: c e -- ^ initial structure - -> (e -> e -> e) -- ^ update function - -> [(IndexOf c, e)] -- ^ association list - -> c e -- ^ result - --------------------------------------------------------------------------- - -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 - arctan2 = vectorZipF ATan2 - scalar x = fromList [x] - konst' = constant - build' = buildV - conj = id - cmap = mapVector - atIndex = (@>) - minIndex = emptyErrorV "minIndex" (round . toScalarF MinIdx) - maxIndex = emptyErrorV "maxIndex" (round . toScalarF MaxIdx) - minElement = emptyErrorV "minElement" (toScalarF Min) - maxElement = emptyErrorV "maxElement" (toScalarF Max) - sumElements = sumF - prodElements = prodF - step = stepF - find = findV - assoc = assocV - accum = accumV - cond = condV condF - -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 - arctan2 = vectorZipR ATan2 - scalar x = fromList [x] - konst' = constant - build' = buildV - conj = id - cmap = mapVector - atIndex = (@>) - minIndex = emptyErrorV "minIndex" (round . toScalarR MinIdx) - maxIndex = emptyErrorV "maxIndex" (round . toScalarR MaxIdx) - minElement = emptyErrorV "minElement" (toScalarR Min) - maxElement = emptyErrorV "maxElement" (toScalarR Max) - sumElements = sumR - prodElements = prodR - step = stepD - find = findV - assoc = assocV - accum = accumV - cond = condV condD - -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 - arctan2 = vectorZipC ATan2 - scalar x = fromList [x] - konst' = constant - build' = buildV - conj = conjugateC - cmap = mapVector - atIndex = (@>) - minIndex = emptyErrorV "minIndex" (minIndex . fst . fromComplex . (mul <*> conj)) - maxIndex = emptyErrorV "maxIndex" (maxIndex . fst . fromComplex . (mul <*> conj)) - minElement = emptyErrorV "minElement" (atIndex <*> minIndex) - maxElement = emptyErrorV "maxElement" (atIndex <*> maxIndex) - sumElements = sumC - prodElements = prodC - step = undefined -- cannot match - find = findV - assoc = assocV - accum = accumV - cond = undefined -- cannot match - -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 - arctan2 = vectorZipQ ATan2 - scalar x = fromList [x] - konst' = constant - build' = buildV - conj = conjugateQ - cmap = mapVector - atIndex = (@>) - minIndex = emptyErrorV "minIndex" (minIndex . fst . fromComplex . (mul <*> conj)) - maxIndex = emptyErrorV "maxIndex" (maxIndex . fst . fromComplex . (mul <*> conj)) - minElement = emptyErrorV "minElement" (atIndex <*> minIndex) - maxElement = emptyErrorV "maxElement" (atIndex <*> maxIndex) - sumElements = sumQ - prodElements = prodQ - step = undefined -- cannot match - find = findV - assoc = assocV - accum = accumV - cond = undefined -- cannot match - ---------------------------------------------------------------- - -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 - arctan2 = liftMatrix2 arctan2 - scalar x = (1><1) [x] - konst' v (r,c) = matrixFromVector RowMajor r c (konst' v (r*c)) - build' = buildM - conj = liftMatrix conj - cmap f = liftMatrix (mapVector f) - atIndex = (@@>) - minIndex = emptyErrorM "minIndex of Matrix" $ - \m -> divMod (minIndex $ flatten m) (cols m) - maxIndex = emptyErrorM "maxIndex of Matrix" $ - \m -> divMod (maxIndex $ flatten m) (cols m) - minElement = emptyErrorM "minElement of Matrix" (atIndex <*> minIndex) - maxElement = emptyErrorM "maxElement of Matrix" (atIndex <*> maxIndex) - sumElements = sumElements . flatten - prodElements = prodElements . flatten - step = liftMatrix step - find = findM - assoc = assocM - accum = accumM - cond = condM - - -emptyErrorV msg f v = - if dim v > 0 - then f v - else error $ msg ++ " of Vector with dim = 0" - -emptyErrorM msg f m = - if rows m > 0 && cols m > 0 - then f m - else error $ msg++" "++shSize m - ----------------------------------------------------- - --- | Matrix product and related functions -class (Num e, Element e) => Product e where - -- | matrix product - multiply :: Matrix e -> Matrix e -> Matrix e - -- | sum of absolute value of elements (differs in complex case from @norm1@) - absSum :: Vector e -> RealOf e - -- | sum of absolute value of elements - norm1 :: Vector e -> RealOf e - -- | euclidean norm - norm2 :: Vector e -> RealOf e - -- | element of maximum magnitude - normInf :: Vector e -> RealOf e - -instance Product Float where - norm2 = emptyVal (toScalarF Norm2) - absSum = emptyVal (toScalarF AbsSum) - norm1 = emptyVal (toScalarF AbsSum) - normInf = emptyVal (maxElement . vectorMapF Abs) - multiply = emptyMul multiplyF - -instance Product Double where - norm2 = emptyVal (toScalarR Norm2) - absSum = emptyVal (toScalarR AbsSum) - norm1 = emptyVal (toScalarR AbsSum) - normInf = emptyVal (maxElement . vectorMapR Abs) - multiply = emptyMul multiplyR - -instance Product (Complex Float) where - norm2 = emptyVal (toScalarQ Norm2) - absSum = emptyVal (toScalarQ AbsSum) - norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapQ Abs) - normInf = emptyVal (maxElement . fst . fromComplex . vectorMapQ Abs) - multiply = emptyMul multiplyQ - -instance Product (Complex Double) where - norm2 = emptyVal (toScalarC Norm2) - absSum = emptyVal (toScalarC AbsSum) - norm1 = emptyVal (sumElements . fst . fromComplex . vectorMapC Abs) - normInf = emptyVal (maxElement . fst . fromComplex . vectorMapC Abs) - multiply = emptyMul multiplyC - -emptyMul m a b - | x1 == 0 && x2 == 0 || r == 0 || c == 0 = konst' 0 (r,c) - | otherwise = m a b - where - r = rows a - x1 = cols a - x2 = rows b - c = cols b - -emptyVal f v = - if dim v > 0 - then f v - else 0 - --- FIXME remove unused C wrappers --- | unconjugated dot product -udot :: Product e => Vector e -> Vector e -> e -udot u v - | dim u == dim v = val (asRow u `multiply` asColumn v) - | otherwise = error $ "different dimensions "++show (dim u)++" and "++show (dim v)++" in dot product" +class Contraction a b c | a b -> c where - val m | dim u > 0 = m@@>(0,0) - | otherwise = 0 + infixl 7 <.> + {- | Matrix product, matrix vector product, and dot product ----------------------------------------------------------- +Examples: --- synonym for matrix product -mXm :: Product t => Matrix t -> Matrix t -> Matrix t -mXm = multiply +>>> let a = (3><4) [1..] :: Matrix Double +>>> let v = fromList [1,0,2,-1] :: Vector Double +>>> let u = fromList [1,2,3] :: Vector Double --- matrix - vector product -mXv :: Product t => Matrix t -> Vector t -> Vector t -mXv m v = flatten $ m `mXm` (asColumn v) +>>> a +(3><4) + [ 1.0, 2.0, 3.0, 4.0 + , 5.0, 6.0, 7.0, 8.0 + , 9.0, 10.0, 11.0, 12.0 ] --- vector - matrix product -vXm :: Product t => Vector t -> Matrix t -> Vector t -vXm v m = flatten $ (asRow v) `mXm` m +matrix × matrix: -{- | Outer product of two vectors. +>>> disp 2 (a <.> trans a) +3x3 + 30 70 110 + 70 174 278 +110 278 446 ->>> 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 +matrix × vector: -------------------------------------------------------------------- +>>> a <.> v +fromList [3.0,11.0,19.0] +dot product: -class Convert t where - real :: Container c t => c (RealOf t) -> c t - complex :: Container c t => c t -> c (ComplexOf t) - single :: Container c t => c t -> c (SingleOf t) - double :: Container c t => c t -> c (DoubleOf t) - toComplex :: (Container c t, RealElement t) => (c t, c t) -> c (Complex t) - fromComplex :: (Container c t, RealElement t) => c (Complex t) -> (c t, c t) +>>> u <.> fromList[3,2,1::Double] +10 +For complex vectors the first argument is conjugated: -instance Convert Double where - real = id - complex = comp' - single = single' - double = id - toComplex = toComplex' - fromComplex = fromComplex' +>>> fromList [1,i] <.> fromList[2*i+1,3] +1.0 :+ (-1.0) -instance Convert Float where - real = id - complex = comp' - single = id - double = double' - toComplex = toComplex' - fromComplex = fromComplex' +>>> fromList [1,i,1-i] <.> complex a +fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0] -instance Convert (Complex Double) where - real = comp' - complex = id - single = single' - double = id - toComplex = toComplex' - fromComplex = fromComplex' - -instance Convert (Complex Float) where - real = comp' - complex = id - single = id - double = double' - toComplex = toComplex' - fromComplex = fromComplex' - -------------------------------------------------------------------- +-} + (<.>) :: a -> b -> c -type family RealOf x -type instance RealOf Double = Double -type instance RealOf (Complex Double) = Double +instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where + u <.> v = conj u `udot` v -type instance RealOf Float = Float -type instance RealOf (Complex Float) = Float +instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where + (<.>) = mXv -type family ComplexOf x +instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where + (<.>) v m = (conj v) `vXm` m -type instance ComplexOf Double = Complex Double -type instance ComplexOf (Complex Double) = Complex Double +instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where + (<.>) = mXm -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 +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 -type instance SingleOf (Complex a) = Complex (SingleOf a) +instance Mul Matrix Matrix Matrix where + (<>) = mXm -type family DoubleOf x +instance Mul Matrix Vector Vector where + (<>) m v = flatten $ m <> asColumn v -type instance DoubleOf Double = Double -type instance DoubleOf Float = Double +instance Mul Vector Matrix Vector where + (<>) v m = flatten $ asRow v <> m -type instance DoubleOf (Complex a) = Complex (DoubleOf a) +-------------------------------------------------------------------------------- -type family ElementOf c +class LSDiv c where + infixl 7 <\> + -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) + (<\>) :: Field t => Matrix t -> c t -> c t -type instance ElementOf (Vector a) = a -type instance ElementOf (Matrix a) = a +instance LSDiv Vector where + m <\> v = flatten (linearSolveSVD m (reshape 1 v)) ------------------------------------------------------------- +instance LSDiv Matrix where + (<\>) = linearSolveSVD -buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ] - where rs = map fromIntegral [0 .. (rc-1)] - cs = map fromIntegral [0 .. (cc-1)] +-------------------------------------------------------------------------------- -buildV n f = fromList [f k | k <- ks] - where ks = map fromIntegral [0 .. (n-1)] +class Konst e d c | d -> c, c -> d + where + -- | + -- >>> konst 7 3 :: Vector Float + -- fromList [7.0,7.0,7.0] + -- + -- >>> konst i (3::Int,4::Int) + -- (3><4) + -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 + -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 + -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ] + -- + konst :: e -> d -> c e --------------------------------------------------------- --- | conjugate transpose -ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e -ctrans = liftMatrix conj . trans +instance Container Vector e => Konst e Int Vector + where + konst = konst' --- | Creates a square matrix with a given diagonal. -diag :: (Num a, Element a) => Vector a -> Matrix a -diag v = diagRect 0 v n n where n = dim v +instance Container Vector e => Konst e (Int,Int) Matrix + where + konst = konst' --- | creates the identity matrix of given dimension -ident :: (Num a, Element a) => Int -> Matrix a -ident n = diag (constant 1 n) +-------------------------------------------------------------------------------- --------------------------------------------------------- +class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f + where + -- | + -- >>> build 5 (**2) :: Vector Double + -- fromList [0.0,1.0,4.0,9.0,16.0] + -- + -- Hilbert matrix of order N: + -- + -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double + -- >>> putStr . dispf 2 $ hilb 3 + -- 3x3 + -- 1.00 0.50 0.33 + -- 0.50 0.33 0.25 + -- 0.33 0.25 0.20 + -- + build :: d -> f -> c e -findV p x = foldVectorWithIndex g [] x where - g k z l = if p z then k:l else l +instance Container Vector e => Build Int (e -> e) Vector e + where + build = build' -findM p x = map ((`divMod` cols x)) $ findV p (flatten x) +instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e + where + build = build' -assocV n z xs = ST.runSTVector $ do - v <- ST.newVector z n - mapM_ (\(k,x) -> ST.writeVector v k x) xs - return v +-------------------------------------------------------------------------------- -assocM (r,c) z xs = ST.runSTMatrix $ do - m <- ST.newMatrix z r c - mapM_ (\((i,j),x) -> ST.writeMatrix m i j x) xs - return m +{- | alternative operator for '(\<.\>)' -accumV v0 f xs = ST.runSTVector $ do - v <- ST.thawVector v0 - mapM_ (\(k,x) -> ST.modifyVector v k (f x)) xs - return v +x25c7, white diamond -accumM m0 f xs = ST.runSTMatrix $ do - m <- ST.thawMatrix m0 - mapM_ (\((i,j),x) -> ST.modifyMatrix m i j (f x)) xs - return m +-} +(◇) :: Contraction a b c => a -> b -> c +infixl 7 ◇ +(◇) = (<.>) ----------------------------------------------------------------------- +-- | dot product: @cdot u v = 'udot' ('conj' u) v@ +dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t +dot u v = udot (conj u) v -condM a b l e t = matrixFromVector RowMajor (rows a'') (cols a'') $ cond a' b' l' e' t' - where - args@(a'':_) = conformMs [a,b,l,e,t] - [a', b', l', e', t'] = map flatten args +-------------------------------------------------------------------------------- -condV f a b l e t = f a' b' l' e' t' - where - [a', b', l', e', t'] = conformVs [a,b,l,e,t] +optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t +optimiseMult = mconcat diff --git a/packages/base/src/Numeric/Chain.hs b/packages/base/src/Numeric/Chain.hs index fbdb01b..c6160e9 100644 --- a/packages/base/src/Numeric/Chain.hs +++ b/packages/base/src/Numeric/Chain.hs @@ -19,7 +19,7 @@ module Numeric.Chain ( import Data.Maybe import Data.Packed.Matrix -import Data.Packed.Numeric +import Data.Packed.Internal.Numeric import qualified Data.Array.IArray as A diff --git a/packages/base/src/Numeric/Container.hs b/packages/base/src/Numeric/Container.hs deleted file mode 100644 index 240e5f5..0000000 --- a/packages/base/src/Numeric/Container.hs +++ /dev/null @@ -1,240 +0,0 @@ -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FunctionalDependencies #-} -{-# LANGUAGE UndecidableInstances #-} - ------------------------------------------------------------------------------ --- | --- Module : Numeric.Container --- Copyright : (c) Alberto Ruiz 2010-14 --- License : BSD3 --- Maintainer : Alberto Ruiz --- Stability : provisional --- --- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines. --- --- The 'Container' class is used to define optimized generic functions which work --- on 'Vector' and 'Matrix' with real or complex elements. --- --- Some of these functions are also available in the instances of the standard --- numeric Haskell classes provided by "Numeric.LinearAlgebra". --- ------------------------------------------------------------------------------ -{-# OPTIONS_HADDOCK hide #-} - -module Numeric.Container ( - -- * Basic functions - module Data.Packed, - konst, build, - linspace, - diag, ident, - ctrans, - -- * Generic operations - Container(..), - -- * Matrix product - Product(..), udot, dot, (◇), - Mul(..), - Contraction(..), - optimiseMult, - mXm,mXv,vXm,LSDiv(..), - outer, kronecker, - -- * Element conversion - Convert(..), - Complexable(), - RealElement(), - - RealOf, ComplexOf, SingleOf, DoubleOf, - - IndexOf, - module Data.Complex -) where - -import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ) -import Data.Packed.Numeric -import Data.Complex -import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD) -import Data.Monoid(Monoid(mconcat)) - ------------------------------------------------------------------- - -{- | Creates a real vector containing a range of values: - ->>> linspace 5 (-3,7::Double) -fromList [-3.0,-0.5,2.0,4.5,7.0]@ - ->>> linspace 5 (8,2+i) :: Vector (Complex Double) -fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0] - -Logarithmic spacing can be defined as follows: - -@logspace n (a,b) = 10 ** linspace n (a,b)@ --} -linspace :: (Container Vector e) => Int -> (e, e) -> Vector e -linspace 0 (a,b) = fromList[(a+b)/2] -linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1] - where s = (b-a)/fromIntegral (n-1) - --------------------------------------------------------- - -class Contraction a b c | a b -> c - where - infixl 7 <.> - {- | Matrix product, matrix vector product, and dot product - -Examples: - ->>> let a = (3><4) [1..] :: Matrix Double ->>> let v = fromList [1,0,2,-1] :: Vector Double ->>> let u = fromList [1,2,3] :: Vector Double - ->>> a -(3><4) - [ 1.0, 2.0, 3.0, 4.0 - , 5.0, 6.0, 7.0, 8.0 - , 9.0, 10.0, 11.0, 12.0 ] - -matrix × matrix: - ->>> disp 2 (a <.> trans a) -3x3 - 30 70 110 - 70 174 278 -110 278 446 - -matrix × vector: - ->>> a <.> v -fromList [3.0,11.0,19.0] - -dot product: - ->>> u <.> fromList[3,2,1::Double] -10 - -For complex vectors the first argument is conjugated: - ->>> fromList [1,i] <.> fromList[2*i+1,3] -1.0 :+ (-1.0) - ->>> fromList [1,i,1-i] <.> complex a -fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0] - --} - (<.>) :: a -> b -> c - - -instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where - u <.> v = conj u `udot` v - -instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where - (<.>) = mXv - -instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where - (<.>) v m = (conj v) `vXm` m - -instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where - (<.>) = mXm - - --------------------------------------------------------------------------------- - -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 - --------------------------------------------------------------------------------- - -class LSDiv c where - infixl 7 <\> - -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) - (<\>) :: Field t => Matrix t -> c t -> c t - -instance LSDiv Vector where - m <\> v = flatten (linearSolveSVD m (reshape 1 v)) - -instance LSDiv Matrix where - (<\>) = linearSolveSVD - --------------------------------------------------------------------------------- - -class Konst e d c | d -> c, c -> d - where - -- | - -- >>> konst 7 3 :: Vector Float - -- fromList [7.0,7.0,7.0] - -- - -- >>> konst i (3::Int,4::Int) - -- (3><4) - -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 - -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 - -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ] - -- - konst :: e -> d -> c e - -instance Container Vector e => Konst e Int Vector - where - konst = konst' - -instance Container Vector e => Konst e (Int,Int) Matrix - where - konst = konst' - --------------------------------------------------------------------------------- - -class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f - where - -- | - -- >>> build 5 (**2) :: Vector Double - -- fromList [0.0,1.0,4.0,9.0,16.0] - -- - -- Hilbert matrix of order N: - -- - -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double - -- >>> putStr . dispf 2 $ hilb 3 - -- 3x3 - -- 1.00 0.50 0.33 - -- 0.50 0.33 0.25 - -- 0.33 0.25 0.20 - -- - build :: d -> f -> c e - -instance Container Vector e => Build Int (e -> e) Vector e - where - build = build' - -instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e - where - build = build' - --------------------------------------------------------------------------------- - -{- | alternative operator for '(\<.\>)' - -x25c7, white diamond - --} -(◇) :: Contraction a b c => a -> b -> c -infixl 7 ◇ -(◇) = (<.>) - --- | dot product: @cdot u v = 'udot' ('conj' u) v@ -dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t -dot u v = udot (conj u) v - --------------------------------------------------------------------------------- - -optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t -optimiseMult = mconcat - diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs index 92761be..063bfc9 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs @@ -81,7 +81,7 @@ import Data.Packed import Numeric.LinearAlgebra.LAPACK as LAPACK import Data.List(foldl1') import Data.Array -import Data.Packed.Numeric +import Data.Packed.Internal.Numeric {- | Generic linear algebra functions for double precision real and complex matrices. diff --git a/packages/base/src/Numeric/LinearAlgebra/Base.hs b/packages/base/src/Numeric/LinearAlgebra/Base.hs index 1af4711..395c84a 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Base.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Base.hs @@ -129,7 +129,7 @@ import Numeric.LinearAlgebra.Data import Numeric.Matrix() import Numeric.Vector() -import Numeric.Container +import Data.Packed.Numeric import Numeric.LinearAlgebra.Algorithms import Numeric.LinearAlgebra.Util diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs index 3bc88f9..2754576 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Data.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs @@ -61,8 +61,7 @@ module Numeric.LinearAlgebra.Data( import Data.Packed.Vector import Data.Packed.Matrix -import Numeric.Container -import Data.Packed.IO +import Data.Packed.Numeric import Numeric.LinearAlgebra.Util import Data.Complex diff --git a/packages/base/src/Numeric/LinearAlgebra/Devel.hs b/packages/base/src/Numeric/LinearAlgebra/Devel.hs index c41db2d..b5ef60d 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Devel.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Devel.hs @@ -60,7 +60,7 @@ module Numeric.LinearAlgebra.Devel( import Data.Packed.Foreign import Data.Packed.Development import Data.Packed.ST -import Numeric.Container(Container,Contraction,LSDiv,Product, +import Data.Packed.Numeric(Container,Contraction,LSDiv,Product, Complexable(),RealElement(), RealOf, ComplexOf, SingleOf, DoubleOf, IndexOf) import Data.Packed diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs index f0470ab..440f6d1 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs @@ -45,25 +45,15 @@ module Numeric.LinearAlgebra.Util( vec, vech, dup, - vtrans, -{- -- * Plot - mplot, - plot, parametricPlot, - splot, mesh, meshdom, - matrixToPGM, imshow, - gnuplotX, gnuplotpdf, gnuplotWin --} + vtrans ) where -import Numeric.Container -import Data.Packed.IO +import Data.Packed.Numeric import Numeric.LinearAlgebra.Algorithms hiding (i) import Numeric.Matrix() import Numeric.Vector() import Numeric.LinearAlgebra.Util.Convolution ---import Graphics.Plot - {- | print a real matrix with given number of digits after the decimal point diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs index 1775f14..3cad8d7 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs @@ -16,7 +16,7 @@ module Numeric.LinearAlgebra.Util.Convolution( corr2, conv2, separable ) where -import Numeric.Container +import Data.Packed.Numeric vectSS :: Element t => Int -> Vector t -> Matrix t diff --git a/packages/base/src/Numeric/Matrix.hs b/packages/base/src/Numeric/Matrix.hs index 962ee84..719b591 100644 --- a/packages/base/src/Numeric/Matrix.hs +++ b/packages/base/src/Numeric/Matrix.hs @@ -27,7 +27,7 @@ module Numeric.Matrix ( ------------------------------------------------------------------- import Data.Packed -import Data.Packed.Numeric +import Data.Packed.Internal.Numeric import qualified Data.Monoid as M import Data.List(partition) import Numeric.Chain diff --git a/packages/base/src/Numeric/Vector.hs b/packages/base/src/Numeric/Vector.hs index 3a425f5..28b453f 100644 --- a/packages/base/src/Numeric/Vector.hs +++ b/packages/base/src/Numeric/Vector.hs @@ -21,7 +21,7 @@ module Numeric.Vector () where import Numeric.Vectorized import Data.Packed.Vector -import Data.Packed.Numeric +import Data.Packed.Internal.Numeric ------------------------------------------------------------------- -- cgit v1.2.3