module Numeric.ContainerBoot (
ident, diag, ctrans,
Container(..),
Product(..),
mXm,mXv,vXm,
outer, kronecker,
Convert(..),
Complexable(),
RealElement(),
RealOf, ComplexOf, SingleOf, DoubleOf,
IndexOf,
module Data.Complex,
build', konst',
(.*),(*/),(<|>),(<->),
vectorMax,vectorMin,
vectorMaxIndex, vectorMinIndex
) where
import Data.Packed
import Numeric.Conversion
import Data.Packed.Internal
import Numeric.GSL.Vector
import Data.Complex
import Control.Monad(ap)
import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ)
import System.IO.Unsafe
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
class (Complexable c, Fractional e, Element e) => Container c e where
scalar :: e -> c e
conj :: c e -> c e
scale :: e -> c e -> c e
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
mul :: c e -> c e -> c e
divide :: c e -> c e -> c e
equal :: c e -> c e -> Bool
arctan2 :: c e -> c e -> c e
cmap :: (Element a, Element b) => (a -> b) -> c a -> c b
konst :: e -> IndexOf c -> c e
build :: IndexOf c -> (ArgOf c e) -> c e
atIndex :: c e -> IndexOf c -> e
minIndex :: c e -> IndexOf c
maxIndex :: c e -> IndexOf c
minElement :: c e -> e
maxElement :: c e -> e
sumElements :: c e -> e
prodElements :: c e -> e
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 = constantD
build = buildV
conj = id
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
arctan2 = vectorZipR ATan2
scalar x = fromList [x]
konst = constantD
build = buildV
conj = id
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
arctan2 = vectorZipC ATan2
scalar x = fromList [x]
konst = constantD
build = buildV
conj = conjugateC
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
arctan2 = vectorZipQ ATan2
scalar x = fromList [x]
konst = constantD
build = buildV
conj = conjugateQ
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
arctan2 = liftMatrix2 arctan2
scalar x = (1><1) [x]
konst v (r,c) = reshape c (konst v (r*c))
build = buildM
conj = liftMatrix conj
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)
maxIndex m = let (r,c) = (rows m,cols m)
i = (maxIndex $ flatten m)
in (i `div` c,i `mod` c)
minElement = ap (@@>) minIndex
maxElement = ap (@@>) maxIndex
sumElements = sumElements . flatten
prodElements = prodElements . flatten
class Element e => Product e where
multiply :: Matrix e -> Matrix e -> Matrix e
dot :: Vector e -> Vector e -> e
absSum :: Vector e -> RealOf e
norm1 :: Vector e -> RealOf e
norm2 :: Vector e -> RealOf e
normInf :: Vector e -> RealOf e
instance Product Float where
norm2 = toScalarF Norm2
absSum = toScalarF AbsSum
dot = dotF
norm1 = toScalarF AbsSum
normInf = maxElement . vectorMapF Abs
multiply = multiplyF
instance Product Double where
norm2 = toScalarR Norm2
absSum = toScalarR AbsSum
dot = dotR
norm1 = toScalarR AbsSum
normInf = maxElement . vectorMapR Abs
multiply = multiplyR
instance Product (Complex Float) where
norm2 = toScalarQ Norm2
absSum = toScalarQ AbsSum
dot = dotQ
norm1 = sumElements . fst . fromComplex . vectorMapQ Abs
normInf = maxElement . fst . fromComplex . vectorMapQ Abs
multiply = multiplyQ
instance Product (Complex Double) where
norm2 = toScalarC Norm2
absSum = toScalarC AbsSum
dot = dotC
norm1 = sumElements . fst . fromComplex . vectorMapC Abs
normInf = maxElement . fst . fromComplex . vectorMapC Abs
multiply = multiplyC
mXm :: Product t => Matrix t -> Matrix t -> Matrix t
mXm = multiply
mXv :: Product t => Matrix t -> Vector t -> Vector t
mXv m v = flatten $ m `mXm` (asColumn v)
vXm :: Product t => Vector t -> Matrix t -> Vector t
vXm v m = flatten $ (asRow v) `mXm` m
outer :: (Product t) => Vector t -> Vector t -> Matrix t
outer u v = asColumn u `multiply` asRow v
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
conjugateAux fun x = unsafePerformIO $ do
v <- createVector (dim x)
app2 fun vec x vec v "conjugateAux"
return v
conjugateQ :: Vector (Complex Float) -> Vector (Complex Float)
conjugateQ = conjugateAux c_conjugateQ
foreign import ccall "conjugateQ" c_conjugateQ :: TQVQV
conjugateC :: Vector (Complex Double) -> Vector (Complex Double)
conjugateC = conjugateAux c_conjugateC
foreign import ccall "conjugateC" c_conjugateC :: TCVCV
infixl 7 .*
a .* x = scale a x
infixl 7 */
v */ x = scale (recip x) v
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 <->
a <|> b = joinH a b
a <-> b = joinV a b
vectorMin :: (Container Vector t, Element t) => Vector t -> t
vectorMin = minElement
vectorMax :: (Container Vector t, Element t) => Vector t -> t
vectorMax = maxElement
vectorMaxIndex :: Vector Double -> Int
vectorMaxIndex = round . toScalarR MaxIdx
vectorMinIndex :: Vector Double -> Int
vectorMinIndex = round . toScalarR MinIdx
class Build f where
build' :: BoundsOf f -> f -> ContainerOf f
type family BoundsOf x
type instance BoundsOf (a->a) = Int
type instance BoundsOf (a->a->a) = (Int,Int)
type family ContainerOf x
type instance ContainerOf (a->a) = Vector a
type instance ContainerOf (a->a->a) = Matrix a
instance (Element a, Num a) => Build (a->a) where
build' = buildV
instance (Element a, Num a) => Build (a->a->a) where
build' = buildM
buildM (rc,cc) f = fromLists [ [f r c | c <- cs] | r <- rs ]
where rs = map fromIntegral [0 .. (rc1)]
cs = map fromIntegral [0 .. (cc1)]
buildV n f = fromList [f k | k <- ks]
where ks = map fromIntegral [0 .. (n1)]
class Konst s where
konst' :: Element e => e -> s -> ContainerOf' s e
type family ContainerOf' x y
type instance ContainerOf' Int a = Vector a
type instance ContainerOf' (Int,Int) a = Matrix a
instance Konst Int where
konst' = constantD
instance Konst (Int,Int) where
konst' k (r,c) = reshape c $ konst' k (r*c)
ctrans :: (Container Vector e, Element e) => Matrix e -> Matrix e
ctrans = liftMatrix conj . trans
diag :: (Num a, Element a) => Vector a -> Matrix a
diag v = diagRect 0 v n n where n = dim v
ident :: (Num a, Element a) => Int -> Matrix a
ident n = diag (constantD 1 n)