From 7376d022b12a27db5a396f89806a709555c1c522 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 15 Jun 2015 12:52:46 +0200 Subject: documentation, more general cond, remove some unicode, minor changes --- packages/base/src/Internal/CG.hs | 8 +- packages/base/src/Internal/Container.hs | 85 ++++++++++++++++++---- packages/base/src/Internal/Element.hs | 38 ++++++++-- packages/base/src/Internal/Matrix.hs | 1 - packages/base/src/Internal/Modular.hs | 9 ++- packages/base/src/Internal/Numeric.hs | 12 +-- packages/base/src/Internal/Util.hs | 38 ++++++---- packages/base/src/Numeric/LinearAlgebra.hs | 48 +++++++++++- packages/base/src/Numeric/LinearAlgebra/Data.hs | 19 +++-- packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | 5 +- packages/base/src/Numeric/LinearAlgebra/Static.hs | 16 ++-- 11 files changed, 216 insertions(+), 63 deletions(-) (limited to 'packages/base/src') diff --git a/packages/base/src/Internal/CG.hs b/packages/base/src/Internal/CG.hs index fd14212..758d130 100644 --- a/packages/base/src/Internal/CG.hs +++ b/packages/base/src/Internal/CG.hs @@ -45,13 +45,13 @@ cg sym at a (CGState p r r2 x _) = CGState p' r' r'2 x' rdx ap1 = a p ap | sym = ap1 | otherwise = at ap1 - pap | sym = p <·> ap1 + pap | sym = p <.> ap1 | otherwise = norm2 ap1 ** 2 alpha = r2 / pap dx = scale alpha p x' = x + dx r' = r - scale alpha ap - r'2 = r' <·> r' + r'2 = r' <.> r' beta = r'2 / r2 p' = r' + scale beta p @@ -75,9 +75,9 @@ solveG mat ma meth rawb x0' ϵb ϵx b = mat rawb x0 = if x0' == 0 then konst 0 (dim b) else x0' r0 = b - a x0 - r20 = r0 <·> r0 + r20 = r0 <.> r0 p0 = r0 - nb2 = b <·> b + nb2 = b <.> b ok CGState {..} = cgr2 --- | An infix synonym for 'dot' -(<.>) :: Numeric t => Vector t -> Vector t -> t -(<.>) = dot +infixr 8 <.> +{- | An infix synonym for 'dot' +>>> vector [1,2,3,4] <.> vector [-2,0,1,1] +5.0 -infixr 8 <·>, #> +>>> let 𝑖 = 0:+1 :: C +>>> fromList [1+𝑖,1] <.> fromList [1,1+𝑖] +2.0 :+ 0.0 -{- | infix synonym for 'dot' +-} ->>> vector [1,2,3,4] <·> vector [-2,0,1,1] -5.0 +(<.>) :: Numeric t => Vector t -> Vector t -> t +(<.>) = dot ->>> let 𝑖 = 0:+1 :: ℂ ->>> fromList [1+𝑖,1] <·> fromList [1,1+𝑖] -2.0 :+ 0.0 -(the dot symbol "·" is obtained by Alt-Gr .) --} -(<·>) :: Numeric t => Vector t -> Vector t -> t -(<·>) = dot {- | infix synonym for 'app' @@ -91,6 +86,7 @@ infixr 8 <·>, #> fromList [140.0,320.0] -} +infixr 8 #> (#>) :: Numeric t => Matrix t -> Vector t -> Vector t (#>) = mXv @@ -264,6 +260,24 @@ instance Numeric Z sortVector :: (Ord t, Element t) => Vector t -> Vector t sortVector = sortV +{- | + +>>> m <- randn 4 10 +>>> disp 2 m +4x10 +-0.31 0.41 0.43 -0.19 -0.17 -0.23 -0.17 -1.04 -0.07 -1.24 + 0.26 0.19 0.14 0.83 -1.54 -0.09 0.37 -0.63 0.71 -0.50 +-0.11 -0.10 -1.29 -1.40 -1.04 -0.89 -0.68 0.35 -1.46 1.86 + 1.04 -0.29 0.19 -0.75 -2.20 -0.01 1.06 0.11 -2.09 -1.58 + +>>> disp 2 $ m ?? (All, Pos $ sortIndex (m!1)) +4x10 +-0.17 -1.04 -1.24 -0.23 0.43 0.41 -0.31 -0.17 -0.07 -0.19 +-1.54 -0.63 -0.50 -0.09 0.14 0.19 0.26 0.37 0.71 0.83 +-1.04 0.35 1.86 -0.89 -1.29 -0.10 -0.11 -0.68 -1.46 -1.40 +-2.20 0.11 -1.58 -0.01 0.19 -0.29 1.04 1.06 -2.09 -0.75 + +-} sortIndex :: (Ord t, Element t) => Vector t -> Vector I sortIndex = sortI @@ -273,11 +287,50 @@ ccompare = ccompare' cselect :: (Container c t) => c I -> c t -> c t -> c t -> c t cselect = cselect' +{- | Extract elements from positions given in matrices of rows and columns. + +>>> r +(3><3) + [ 1, 1, 1 + , 1, 2, 2 + , 1, 2, 3 ] +>>> c +(3><3) + [ 0, 1, 5 + , 2, 2, 1 + , 4, 4, 1 ] +>>> m +(4><6) + [ 0, 1, 2, 3, 4, 5 + , 6, 7, 8, 9, 10, 11 + , 12, 13, 14, 15, 16, 17 + , 18, 19, 20, 21, 22, 23 ] + +>>> remap r c m +(3><3) + [ 6, 7, 11 + , 8, 14, 13 + , 10, 16, 19 ] + +The indexes are autoconformable. + +>>> c' +(3><1) + [ 1 + , 2 + , 4 ] +>>> remap r c' m +(3><3) + [ 7, 7, 7 + , 8, 14, 14 + , 10, 16, 22 ] + +-} remap :: Element t => Matrix I -> Matrix I -> Matrix t -> Matrix t remap i j m | minElement i >= 0 && maxElement i < fromIntegral (rows m) && minElement j >= 0 && maxElement j < fromIntegral (cols m) = remapM i' j' m - | otherwise = error $ "out of range index in rmap" + | otherwise = error $ "out of range index in remap" where [i',j'] = conformMs [i,j] diff --git a/packages/base/src/Internal/Element.hs b/packages/base/src/Internal/Element.hs index 4007491..51d5686 100644 --- a/packages/base/src/Internal/Element.hs +++ b/packages/base/src/Internal/Element.hs @@ -80,7 +80,7 @@ breakAt c l = (a++[c],tail b) where (a,b) = break (==c) l -------------------------------------------------------------------------------- - +-- | Specification of indexes for the operator '??'. data Extractor = All | Range Int Int Int @@ -102,7 +102,32 @@ ppext (Drop n) = printf "Drop %d" n ppext (TakeLast n) = printf "TakeLast %d" n ppext (DropLast n) = printf "DropLast %d" n +{- | General matrix slicing. + +>>> m +(4><5) + [ 0, 1, 2, 3, 4 + , 5, 6, 7, 8, 9 + , 10, 11, 12, 13, 14 + , 15, 16, 17, 18, 19 ] + +>>> m ?? (Take 3, DropLast 2) +(3><3) + [ 0, 1, 2 + , 5, 6, 7 + , 10, 11, 12 ] + +>>> m ?? (Pos (idxs[2,1]), All) +(2><5) + [ 10, 11, 12, 13, 14 + , 5, 6, 7, 8, 9 ] + +>>> m ?? (PosCyc (idxs[-7,80]), Range 4 (-2) 0) +(2><3) + [ 9, 7, 5 + , 4, 2, 0 ] +-} infixl 9 ?? (??) :: Element t => Matrix t -> (Extractor,Extractor) -> Matrix t @@ -328,27 +353,30 @@ r >< c = f where ---------------------------------------------------------------- --- | Creates a matrix with the first n rows of another matrix takeRows :: Element t => Int -> Matrix t -> Matrix t takeRows n mt = subMatrix (0,0) (n, cols mt) mt + -- | Creates a matrix with the last n rows of another matrix takeLastRows :: Element t => Int -> Matrix t -> Matrix t takeLastRows n mt = subMatrix (rows mt - n, 0) (n, cols mt) mt --- | Creates a copy of a matrix without the first n rows + dropRows :: Element t => Int -> Matrix t -> Matrix t dropRows n mt = subMatrix (n,0) (rows mt - n, cols mt) mt + -- | Creates a copy of a matrix without the last n rows dropLastRows :: Element t => Int -> Matrix t -> Matrix t dropLastRows n mt = subMatrix (0,0) (rows mt - n, cols mt) mt --- |Creates a matrix with the first n columns of another matrix + takeColumns :: Element t => Int -> Matrix t -> Matrix t takeColumns n mt = subMatrix (0,0) (rows mt, n) mt + -- |Creates a matrix with the last n columns of another matrix takeLastColumns :: Element t => Int -> Matrix t -> Matrix t takeLastColumns n mt = subMatrix (0, cols mt - n) (rows mt, n) mt --- | Creates a copy of a matrix without the first n columns + dropColumns :: Element t => Int -> Matrix t -> Matrix t dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt + -- | Creates a copy of a matrix without the last n columns dropLastColumns :: Element t => Int -> Matrix t -> Matrix t dropLastColumns n mt = subMatrix (0,0) (rows mt, cols mt - n) mt diff --git a/packages/base/src/Internal/Matrix.hs b/packages/base/src/Internal/Matrix.hs index e4b1226..75e92a5 100644 --- a/packages/base/src/Internal/Matrix.hs +++ b/packages/base/src/Internal/Matrix.hs @@ -378,7 +378,6 @@ foreign import ccall unsafe "transL" ctransL :: TMM Z ---------------------------------------------------------------------- --- | Extracts a submatrix from a matrix. subMatrix :: Element a => (Int,Int) -- ^ (r0,c0) starting position -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix diff --git a/packages/base/src/Internal/Modular.hs b/packages/base/src/Internal/Modular.hs index 3b27310..0d906bb 100644 --- a/packages/base/src/Internal/Modular.hs +++ b/packages/base/src/Internal/Modular.hs @@ -9,8 +9,8 @@ {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE UndecidableInstances #-} {-# LANGUAGE GADTs #-} -{-# LANGUAGE TypeFamilies #-} - +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} {- | Module : Internal.Modular @@ -23,7 +23,7 @@ Proof of concept of statically checked modular arithmetic. -} module Internal.Modular( - Mod + Mod, type (./.) ) where import Internal.Vector @@ -49,6 +49,9 @@ import Data.Complex newtype Mod (n :: Nat) t = Mod {unMod:: t} deriving (Storable) +infixr 5 ./. +type (./.) x n = Mod n x + instance (Integral t, Enum t, KnownNat m) => Enum (Mod m t) where toEnum = l0 (\m x -> fromIntegral $ x `mod` (fromIntegral m)) diff --git a/packages/base/src/Internal/Numeric.hs b/packages/base/src/Internal/Numeric.hs index ca17c23..efcde2c 100644 --- a/packages/base/src/Internal/Numeric.hs +++ b/packages/base/src/Internal/Numeric.hs @@ -481,14 +481,16 @@ step = step' -- , 0.0, 100.0, 7.0, 8.0 -- , 0.0, 0.0, 100.0, 12.0 ] -- +-- >>> let chop x = cond (abs x) 1E-6 0 0 x +-- cond - :: (Ord e, Container c e) + :: (Ord e, Container c e, Container c x) => c e -- ^ a -> c e -- ^ b - -> c e -- ^ l - -> c e -- ^ e - -> c e -- ^ g - -> c e -- ^ result + -> c x -- ^ l + -> c x -- ^ e + -> c x -- ^ g + -> c x -- ^ result cond a b l e g = cselect' (ccompare' a b) l e g diff --git a/packages/base/src/Internal/Util.hs b/packages/base/src/Internal/Util.hs index 09ba21c..f08f710 100644 --- a/packages/base/src/Internal/Util.hs +++ b/packages/base/src/Internal/Util.hs @@ -85,7 +85,7 @@ type ℤ = Int type ℂ = Complex Double -- | imaginary unit -iC :: ℂ +iC :: C iC = 0:+1 {- | Create a real vector. @@ -94,7 +94,7 @@ iC = 0:+1 fromList [1.0,2.0,3.0,4.0,5.0] -} -vector :: [ℝ] -> Vector ℝ +vector :: [R] -> Vector R vector = fromList {- | Create a real matrix. @@ -108,8 +108,8 @@ vector = fromList -} matrix :: Int -- ^ number of columns - -> [ℝ] -- ^ elements in row order - -> Matrix ℝ + -> [R] -- ^ elements in row order + -> Matrix R matrix c = reshape c . fromList @@ -260,34 +260,34 @@ norm = pnorm PNorm2 class Normed a where - norm_0 :: a -> ℝ - norm_1 :: a -> ℝ - norm_2 :: a -> ℝ - norm_Inf :: a -> ℝ + norm_0 :: a -> R + norm_1 :: a -> R + norm_2 :: a -> R + norm_Inf :: a -> R -instance Normed (Vector ℝ) +instance Normed (Vector R) where norm_0 v = sumElements (step (abs v - scalar (eps*normInf v))) norm_1 = pnorm PNorm1 norm_2 = pnorm PNorm2 norm_Inf = pnorm Infinity -instance Normed (Vector ℂ) +instance Normed (Vector C) where norm_0 v = sumElements (step (fst (fromComplex (abs v)) - scalar (eps*normInf v))) norm_1 = pnorm PNorm1 norm_2 = pnorm PNorm2 norm_Inf = pnorm Infinity -instance Normed (Matrix ℝ) +instance Normed (Matrix R) where norm_0 = norm_0 . flatten norm_1 = pnorm PNorm1 norm_2 = pnorm PNorm2 norm_Inf = pnorm Infinity -instance Normed (Matrix ℂ) +instance Normed (Matrix C) where norm_0 = norm_0 . flatten norm_1 = pnorm PNorm1 @@ -323,12 +323,22 @@ instance Normed (Vector (Complex Float)) norm_Inf = norm_Inf . double -norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> ℝ +norm_Frob :: (Normed (Vector t), Element t) => Matrix t -> R norm_Frob = norm_2 . flatten -norm_nuclear :: Field t => Matrix t -> ℝ +norm_nuclear :: Field t => Matrix t -> R norm_nuclear = sumElements . singularValues +{- | Check if the absolute value or complex magnitude is greater than a given threshold + +>>> magnit 1E-6 (1E-12 :: R) +False +>>> magnit 1E-6 (3+iC :: C) +True +>>> magnit 0 (3 :: I ./. 5) +True + +-} magnit :: (Element t, Normed (Vector t)) => R -> t -> Bool magnit e x = norm_1 (fromList [x]) > e diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs index fe524cc..2e6b8ca 100644 --- a/packages/base/src/Numeric/LinearAlgebra.hs +++ b/packages/base/src/Numeric/LinearAlgebra.hs @@ -47,7 +47,7 @@ module Numeric.LinearAlgebra ( -- * Products -- ** dot - dot, (<·>), + dot, (<.>), -- ** matrix-vector app, (#>), (<#), (!#>), -- ** matrix-matrix @@ -240,7 +240,53 @@ nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m) -- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'. orth m = orthSVD (Left (1*eps)) m (leftSV m) +{- | Experimental implementation of LU factorization + working on any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'. + +>>> let m = ident 5 + (5><5) [0..] :: Matrix (Z ./. 17) +(5><5) + [ 1, 1, 2, 3, 4 + , 5, 7, 7, 8, 9 + , 10, 11, 13, 13, 14 + , 15, 16, 0, 2, 2 + , 3, 4, 5, 6, 8 ] + +>>> let (l,u,p,s) = luFact $ luPacked' m +>>> l +(5><5) + [ 1, 0, 0, 0, 0 + , 6, 1, 0, 0, 0 + , 12, 7, 1, 0, 0 + , 7, 10, 7, 1, 0 + , 8, 2, 6, 11, 1 ] +>>> u +(5><5) + [ 15, 16, 0, 2, 2 + , 0, 13, 7, 13, 14 + , 0, 0, 15, 0, 11 + , 0, 0, 0, 15, 15 + , 0, 0, 0, 0, 1 ] + +-} luPacked' x = mutable (luST (magnit 0)) x +{- | Experimental implementation of gaussian elimination + working on any Fractional element type, including 'Mod' n 'I' and 'Mod' n 'Z'. + +>>> let a = (2><2) [1,2,3,5] :: Matrix (Z ./. 13) +(2><2) + [ 1, 2 + , 3, 5 ] +>>> b +(2><3) + [ 5, 1, 3 + , 8, 6, 3 ] + +>>> let x = linearSolve' a b +(2><3) + [ 4, 7, 4 + , 7, 10, 6 ] + +-} linearSolve' x y = gaussElim x y diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs index fffc2bd..2a45771 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Data.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs @@ -1,3 +1,5 @@ +{-# LANGUAGE TypeOperators #-} + -------------------------------------------------------------------------------- {- | Module : Numeric.LinearAlgebra.Data @@ -17,11 +19,11 @@ module Numeric.LinearAlgebra.Data( -- | 1D arrays are storable vectors from the vector package. There is no distinction -- between row and column vectors. - fromList, toList, vector, (|>), + fromList, toList, (|>), vector, range, idxs, -- * Matrix - matrix, (><), tr, tr', + (><), matrix, tr, tr', -- * Dimensions @@ -43,7 +45,7 @@ module Numeric.LinearAlgebra.Data( Indexable(..), -- * Construction - scalar, Konst(..), Build(..), assoc, accum, linspace, range, idxs, -- ones, zeros, + scalar, Konst(..), Build(..), assoc, accum, linspace, -- ones, zeros, -- * Diagonal ident, diag, diagl, diagRect, takeDiag, @@ -53,16 +55,19 @@ module Numeric.LinearAlgebra.Data( -- * Matrix extraction Extractor(..), (??), + takeRows, dropRows, takeColumns, dropColumns, - subMatrix, (?), (¿), fliprl, flipud, remap, + subMatrix, (?), (¿), fliprl, flipud, + + remap, -- * Block matrix fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery, -- * Mapping functions conj, cmap, cmod, - - step, cond, ccompare, cselect, + + step, cond, -- * Find elements find, maxIndex, minIndex, maxElement, minElement, @@ -87,7 +92,7 @@ module Numeric.LinearAlgebra.Data( separable, fromArray2D, module Data.Complex, - R,C,I,Z,Mod, + R,C,I,Z,Mod, type(./.), Vector, Matrix, GMatrix, nRows, nCols ) where diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs index 11c2487..8adaaaf 100644 --- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs +++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs @@ -13,10 +13,13 @@ compatibility with previous version, to be removed module Numeric.LinearAlgebra.HMatrix ( module Numeric.LinearAlgebra, - (¦),(——),ℝ,ℂ, + (¦),(——),ℝ,ℂ,(<·>) ) where import Numeric.LinearAlgebra import Internal.Util +infixr 8 <·> +(<·>) :: Numeric t => Vector t -> Vector t -> t +(<·>) = dot diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs index a657bd0..d0a790d 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs @@ -44,7 +44,7 @@ module Numeric.LinearAlgebra.Static( -- * Complex C, M, Her, her, 𝑖, -- * Products - (<>),(#>),(<·>), + (<>),(#>),(<.>), -- * Linear Systems linSolve, (<\>), -- * Factorizations @@ -55,13 +55,13 @@ module Numeric.LinearAlgebra.Static( Disp(..), Domain(..), withVector, withMatrix, toRows, toColumns, - Sized(..), Diag(..), Sym, sym, mTm, unSym + Sized(..), Diag(..), Sym, sym, mTm, unSym, (<·>) ) where import GHC.TypeLits import Numeric.LinearAlgebra hiding ( - (<>),(#>),(<·>),Konst(..),diag, disp,(===),(|||), + (<>),(#>),(<.>),Konst(..),diag, disp,(===),(|||), row,col,vector,matrix,linspace,toRows,toColumns, (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH', eigenvalues,eigenvaluesSH,eigenvaluesSH',build, @@ -207,6 +207,10 @@ infixr 8 <·> (<·>) :: R n -> R n -> ℝ (<·>) = dotR +infixr 8 <.> +(<.>) :: R n -> R n -> ℝ +(<.>) = dotR + -------------------------------------------------------------------------------- class Diag m d | m -> d @@ -496,7 +500,7 @@ appC m v = mkC (extract m LA.#> extract v) dotC :: KnownNat n => C n -> C n -> ℂ dotC (unwrap -> u) (unwrap -> v) | singleV u || singleV v = sumElements (conj u * v) - | otherwise = u LA.<·> v + | otherwise = u LA.<.> v crossC :: C 3 -> C 3 -> C 3 @@ -584,12 +588,12 @@ test = (ok,info) where q = tm :: L 10 3 - thingD = vjoin [ud1 u, 1] LA.<·> tr m LA.#> m LA.#> ud1 v + thingD = vjoin [ud1 u, 1] LA.<.> tr m LA.#> m LA.#> ud1 v where m = LA.matrix 3 [1..30] precS = (1::Double) + (2::Double) * ((1 :: R 3) * (u & 6)) <·> konst 2 #> v - precD = 1 + 2 * vjoin[ud1 u, 6] LA.<·> LA.konst 2 (LA.size (ud1 u) +1, LA.size (ud1 v)) LA.#> ud1 v + precD = 1 + 2 * vjoin[ud1 u, 6] LA.<.> LA.konst 2 (LA.size (ud1 u) +1, LA.size (ud1 v)) LA.#> ud1 v splittest -- cgit v1.2.3