From 5b6de561f131d75049fdb999e98a07939ec2e8e7 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sat, 24 May 2014 13:32:58 +0200 Subject: backward compatibility --- packages/base/src/Data/Packed/Development.hs | 3 +- packages/base/src/Data/Packed/IO.hs | 2 +- packages/base/src/Data/Packed/Internal/Numeric.hs | 39 ++- packages/base/src/Data/Packed/Matrix.hs | 13 +- packages/base/src/Data/Packed/Numeric.hs | 288 ++++++++++++++++++++++ packages/base/src/Data/Packed/Vector.hs | 16 +- 6 files changed, 326 insertions(+), 35 deletions(-) create mode 100644 packages/base/src/Data/Packed/Numeric.hs (limited to 'packages/base/src/Data') diff --git a/packages/base/src/Data/Packed/Development.hs b/packages/base/src/Data/Packed/Development.hs index 1efedc9..72eb16b 100644 --- a/packages/base/src/Data/Packed/Development.hs +++ b/packages/base/src/Data/Packed/Development.hs @@ -25,8 +25,7 @@ module Data.Packed.Development ( unsafeFromForeignPtr, unsafeToForeignPtr, check, (//), - at', atM', fi, table, - conformMs, conformVs, shSize, splitEvery + at', atM', fi ) where import Data.Packed.Internal diff --git a/packages/base/src/Data/Packed/IO.hs b/packages/base/src/Data/Packed/IO.hs index 94fb30a..f7afa80 100644 --- a/packages/base/src/Data/Packed/IO.hs +++ b/packages/base/src/Data/Packed/IO.hs @@ -18,12 +18,12 @@ module Data.Packed.IO ( ) where import Data.Packed -import Data.Packed.Development import Text.Printf(printf) import Data.List(intersperse) import Data.Complex import Numeric.Vectorized(vectorScan,saveMatrix) import Control.Applicative((<$>)) +import Data.Packed.Internal {- | Creates a string from a matrix given a separator and a function to show each entry. Using this function the user can easily define any desired display function: diff --git a/packages/base/src/Data/Packed/Internal/Numeric.hs b/packages/base/src/Data/Packed/Internal/Numeric.hs index 9cd18df..3c1c1d0 100644 --- a/packages/base/src/Data/Packed/Internal/Numeric.hs +++ b/packages/base/src/Data/Packed/Internal/Numeric.hs @@ -49,6 +49,7 @@ import Data.Complex import Control.Applicative((<*>)) import Numeric.LinearAlgebra.LAPACK(multiplyR,multiplyC,multiplyF,multiplyQ) +import Data.Packed.Internal ------------------------------------------------------------------- @@ -65,7 +66,9 @@ 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 +class (Complexable c, Fractional e, Element e) => Container c e + where + size' :: c e -> IndexOf c scalar' :: e -> c e conj' :: c e -> c e scale' :: e -> c e -> c e @@ -114,7 +117,9 @@ class (Complexable c, Fractional e, Element e) => Container c e where -------------------------------------------------------------------------- -instance Container Vector Float where +instance Container Vector Float + where + size' = dim scale' = vectorMapValF Scale scaleRecip = vectorMapValF Recip addConstant = vectorMapValF AddConstant @@ -125,7 +130,7 @@ instance Container Vector Float where equal u v = dim u == dim v && maxElement (vectorMapF Abs (sub u v)) == 0.0 arctan2' = vectorZipF ATan2 scalar' x = fromList [x] - konst' = constant + konst' = constantD build' = buildV conj' = id cmap' = mapVector @@ -142,7 +147,9 @@ instance Container Vector Float where accum' = accumV cond' = condV condF -instance Container Vector Double where +instance Container Vector Double + where + size' = dim scale' = vectorMapValR Scale scaleRecip = vectorMapValR Recip addConstant = vectorMapValR AddConstant @@ -153,7 +160,7 @@ instance Container Vector Double where equal u v = dim u == dim v && maxElement (vectorMapR Abs (sub u v)) == 0.0 arctan2' = vectorZipR ATan2 scalar' x = fromList [x] - konst' = constant + konst' = constantD build' = buildV conj' = id cmap' = mapVector @@ -170,7 +177,9 @@ instance Container Vector Double where accum' = accumV cond' = condV condD -instance Container Vector (Complex Double) where +instance Container Vector (Complex Double) + where + size' = dim scale' = vectorMapValC Scale scaleRecip = vectorMapValC Recip addConstant = vectorMapValC AddConstant @@ -181,7 +190,7 @@ instance Container Vector (Complex Double) where equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 arctan2' = vectorZipC ATan2 scalar' x = fromList [x] - konst' = constant + konst' = constantD build' = buildV conj' = conjugateC cmap' = mapVector @@ -198,7 +207,9 @@ instance Container Vector (Complex Double) where accum' = accumV cond' = undefined -- cannot match -instance Container Vector (Complex Float) where +instance Container Vector (Complex Float) + where + size' = dim scale' = vectorMapValQ Scale scaleRecip = vectorMapValQ Recip addConstant = vectorMapValQ AddConstant @@ -209,7 +220,7 @@ instance Container Vector (Complex Float) where equal u v = dim u == dim v && maxElement (mapVector magnitude (sub u v)) == 0.0 arctan2' = vectorZipQ ATan2 scalar' x = fromList [x] - konst' = constant + konst' = constantD build' = buildV conj' = conjugateQ cmap' = mapVector @@ -228,7 +239,9 @@ instance Container Vector (Complex Float) where --------------------------------------------------------------- -instance (Container Vector a) => Container Matrix a where +instance (Container Vector a) => Container Matrix a + where + size' = size scale' x = liftMatrix (scale' x) scaleRecip x = liftMatrix (scaleRecip x) addConstant x = liftMatrix (addConstant x) @@ -637,7 +650,7 @@ 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) +ident n = diag (constantD 1 n) -------------------------------------------------------- @@ -681,8 +694,12 @@ condV f a b l e t = f a' b' l' e' t' class Transposable t where + -- | (conjugate) transpose tr :: t -> t +instance (Container Vector t) => Transposable (Matrix t) + where + tr = ctrans class Linear t v where diff --git a/packages/base/src/Data/Packed/Matrix.hs b/packages/base/src/Data/Packed/Matrix.hs index b3be823..2acb31a 100644 --- a/packages/base/src/Data/Packed/Matrix.hs +++ b/packages/base/src/Data/Packed/Matrix.hs @@ -226,15 +226,12 @@ takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (co ------------------------------------------------------------ -{- | An easy way to create a matrix: +{- | create a general matrix ->>> (2><3)[2,4,7,-3,11,0] +>>> (2><3) [2, 4, 7+2*𝑖, -3, 11, 0] (2><3) - [ 2.0, 4.0, 7.0 - , -3.0, 11.0, 0.0 ] - -This is the format produced by the instances of Show (Matrix a), which -can also be used for input. + [ 2.0 :+ 0.0, 4.0 :+ 0.0, 7.0 :+ 2.0 + , (-3.0) :+ (-0.0), 11.0 :+ 0.0, 0.0 :+ 0.0 ] The input list is explicitly truncated, so that it can safely be used with lists that are too long (like infinite lists). @@ -244,6 +241,8 @@ safely be used with lists that are too long (like infinite lists). [ 1.0, 2.0, 3.0 , 4.0, 5.0, 6.0 ] +This is the format produced by the instances of Show (Matrix a), which +can also be used for input. -} (><) :: (Storable a) => Int -> Int -> [a] -> Matrix a diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs new file mode 100644 index 0000000..01cf6c5 --- /dev/null +++ b/packages/base/src/Data/Packed/Numeric.hs @@ -0,0 +1,288 @@ +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE UndecidableInstances #-} + +----------------------------------------------------------------------------- +-- | +-- Module : Data.Packed.Numeric +-- 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 Data.Packed.Numeric ( + -- * Basic functions + module Data.Packed, + Konst(..), Build(..), + linspace, + diag, ident, + ctrans, + -- * Generic operations + Container(..), + -- add, mul, sub, divide, equal, scaleRecip, addConstant, + scalar, conj, scale, arctan2, cmap, + atIndex, minIndex, maxIndex, minElement, maxElement, + sumElements, prodElements, + step, cond, find, assoc, accum, + Transposable(..), Linear(..), + -- * Matrix product + Product(..), udot, dot, (◇), + Mul(..), + Contraction(..),(<.>), + optimiseMult, + mXm,mXv,vXm,LSDiv,(<\>), + outer, kronecker, + -- * Random numbers + RandDist(..), + randomVector, + gaussianSample, + uniformSample, + meanCov, + -- * Element conversion + Convert(..), + Complexable(), + RealElement(), + + RealOf, ComplexOf, SingleOf, DoubleOf, + + IndexOf, + module Data.Complex, + -- * IO + module Data.Packed.IO, + -- * Misc + Testable(..) +) where + +import Data.Packed +import Data.Packed.Internal.Numeric +import Data.Complex +import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD) +import Data.Monoid(Monoid(mconcat)) +import Data.Packed.IO +import Numeric.LinearAlgebra.Random + +------------------------------------------------------------------ + +{- | 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) + +-------------------------------------------------------- + +{- | Matrix product, matrix - vector product, and dot product (equivalent to 'contraction') + +(This operator can also be written using the unicode symbol ◇ (25c7).) + +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] +-} +infixl 7 <.> +(<.>) :: Contraction a b c => a -> b -> c +(<.>) = contraction + + +class Contraction a b c | a b -> c + where + -- | Matrix product, matrix - vector product, and dot product + contraction :: a -> b -> c + +instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where + u `contraction` v = conj u `udot` v + +instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where + contraction = mXv + +instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where + contraction v m = (conj v) `vXm` m + +instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where + contraction = 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 + +-------------------------------------------------------------------------------- + +-- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD) +infixl 7 <\> +(<\>) :: (LSDiv c, Field t) => Matrix t -> c t -> c t +(<\>) = linSolve + +class LSDiv c + where + linSolve :: Field t => Matrix t -> c t -> c t + +instance LSDiv Vector + where + linSolve m v = flatten (linearSolveSVD m (reshape 1 v)) + +instance LSDiv Matrix + where + linSolve = 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 unicode symbol (25c7) for 'contraction' +(◇) :: Contraction a b c => a -> b -> c +infixl 7 ◇ +(◇) = contraction + +-- | 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 + +-------------------------------------------------------------------------------- + + +{- | Compute mean vector and covariance matrix of the rows of a matrix. + +>>> meanCov $ gaussianSample 666 1000 (fromList[4,5]) (diagl[2,3]) +(fromList [4.010341078059521,5.0197204699640405], +(2><2) + [ 1.9862461923890056, -1.0127225830525157e-2 + , -1.0127225830525157e-2, 3.0373954915729318 ]) + +-} +meanCov :: Matrix Double -> (Vector Double, Matrix Double) +meanCov x = (med,cov) where + r = rows x + k = 1 / fromIntegral r + med = konst k r `vXm` x + meds = konst 1 r `outer` med + xc = x `sub` meds + cov = scale (recip (fromIntegral (r-1))) (trans xc `mXm` xc) + +-------------------------------------------------------------------------------- + diff --git a/packages/base/src/Data/Packed/Vector.hs b/packages/base/src/Data/Packed/Vector.hs index 53fe563..31dcf47 100644 --- a/packages/base/src/Data/Packed/Vector.hs +++ b/packages/base/src/Data/Packed/Vector.hs @@ -17,17 +17,15 @@ module Data.Packed.Vector ( Vector, - fromList, (|>), toList, buildVector, constant, + fromList, (|>), toList, buildVector, dim, (@>), subVector, takesV, vjoin, join, mapVector, mapVectorWithIndex, zipVector, zipVectorWith, unzipVector, unzipVectorWith, mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_, - foldLoop, foldVector, foldVectorG, foldVectorWithIndex, - stepD, stepF, condD, condF, conjugateC, conjugateQ + foldLoop, foldVector, foldVectorG, foldVectorWithIndex ) where import Data.Packed.Internal.Vector -import Data.Packed.Internal.Matrix import Foreign.Storable ------------------------------------------------------------------- @@ -95,13 +93,3 @@ unzipVector = unzipVectorWith id join :: Storable t => [Vector t] -> Vector t join = vjoin -{- | creates a vector with a given number of equal components: - -@> constant 2 7 -7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ --} -constant :: Element a => a -> Int -> Vector a --- constant x n = runSTVector (newVector x n) -constant = constantD-- about 2x faster - - -- cgit v1.2.3