From a928a3a1713704cf3d5148bedc7ff8acb1347599 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Tue, 10 Jun 2014 13:08:17 +0200 Subject: module and function names --- packages/base/hmatrix.cabal | 11 +- packages/base/src/Data/Packed/Numeric.hs | 4 +- packages/base/src/Numeric/Complex.hs | 72 +++ packages/base/src/Numeric/HMatrix.hs | 634 +++++++++++++++------ packages/base/src/Numeric/LinearAlgebra/Complex.hs | 80 --- packages/base/src/Numeric/LinearAlgebra/Data.hs | 4 +- packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | 201 +++++++ packages/base/src/Numeric/LinearAlgebra/Real.hs | 480 ---------------- packages/base/src/Numeric/LinearAlgebra/Static.hs | 2 +- packages/base/src/Numeric/LinearAlgebra/Util.hs | 14 +- packages/base/src/Numeric/Sparse.hs | 2 +- 11 files changed, 757 insertions(+), 747 deletions(-) create mode 100644 packages/base/src/Numeric/Complex.hs delete mode 100644 packages/base/src/Numeric/LinearAlgebra/Complex.hs create mode 100644 packages/base/src/Numeric/LinearAlgebra/HMatrix.hs delete mode 100644 packages/base/src/Numeric/LinearAlgebra/Real.hs (limited to 'packages') diff --git a/packages/base/hmatrix.cabal b/packages/base/hmatrix.cabal index 30d88cd..3f16f5f 100644 --- a/packages/base/hmatrix.cabal +++ b/packages/base/hmatrix.cabal @@ -47,14 +47,14 @@ library Numeric.HMatrix Numeric.LinearAlgebra.Devel Numeric.LinearAlgebra.Data - Numeric.LinearAlgebra.Real - -- Numeric.LinearAlgebra.Complex + Numeric.LinearAlgebra.HMatrix + other-modules: Data.Packed.Internal, - Data.Packed.Internal.Common, - Data.Packed.Internal.Signatures, - Data.Packed.Internal.Vector, + Data.Packed.Internal.Common + Data.Packed.Internal.Signatures + Data.Packed.Internal.Vector Data.Packed.Internal.Matrix Data.Packed.IO Numeric.Chain @@ -68,6 +68,7 @@ library Numeric.LinearAlgebra.Random Numeric.Conversion Numeric.Sparse + Numeric.Complex Numeric.LinearAlgebra.Static C-sources: src/C/lapack-aux.c diff --git a/packages/base/src/Data/Packed/Numeric.hs b/packages/base/src/Data/Packed/Numeric.hs index d2a20be..e59c1cd 100644 --- a/packages/base/src/Data/Packed/Numeric.hs +++ b/packages/base/src/Data/Packed/Numeric.hs @@ -106,7 +106,7 @@ infixr 8 <ยท>, #> {- | dot product ->>> vect [1,2,3,4] <ยท> vect [-2,0,1,1] +>>> vector [1,2,3,4] <ยท> vector [-2,0,1,1] 5.0 >>> let ๐‘– = 0:+1 :: โ„‚ @@ -128,7 +128,7 @@ infixr 8 <ยท>, #> [ 1.0, 2.0, 3.0 , 4.0, 5.0, 6.0 ] ->>> let v = vect [10,20,30] +>>> let v = vector [10,20,30] >>> m #> v fromList [140.0,320.0] diff --git a/packages/base/src/Numeric/Complex.hs b/packages/base/src/Numeric/Complex.hs new file mode 100644 index 0000000..d194af3 --- /dev/null +++ b/packages/base/src/Numeric/Complex.hs @@ -0,0 +1,72 @@ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE EmptyDataDecls #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +{-# LANGUAGE GADTs #-} + + +{- | +Module : Numeric.HMatrix.Static.Complex +Copyright : (c) Alberto Ruiz 2006-14 +License : BSD3 +Stability : experimental + +-} + +module Numeric.Complex( + C, M, + vec2, vec3, vec4, (&), (#), + vect, + Her, her, ๐‘–, +) where + +import GHC.TypeLits +import Numeric.LinearAlgebra.Util(โ„‚,iC) +import qualified Numeric.LinearAlgebra.HMatrix as LA +import Numeric.LinearAlgebra.Static + + +๐‘– :: Sized โ„‚ s c => s +๐‘– = konst iC + +newtype Her n = Her (M n n) + +her :: KnownNat n => M n n -> Her n +her m = Her $ (m + LA.tr m)/2 + + + + +infixl 4 & +(&) :: forall n . KnownNat n + => C n -> โ„‚ -> C (n+1) +u & x = u # (mkC (LA.scalar x) :: C 1) + +infixl 4 # +(#) :: forall n m . (KnownNat n, KnownNat m) + => C n -> C m -> C (n+m) +(C u) # (C v) = C (vconcat u v) + + + +vec2 :: โ„‚ -> โ„‚ -> C 2 +vec2 a b = C (gvec2 a b) + +vec3 :: โ„‚ -> โ„‚ -> โ„‚ -> C 3 +vec3 a b c = C (gvec3 a b c) + + +vec4 :: โ„‚ -> โ„‚ -> โ„‚ -> โ„‚ -> C 4 +vec4 a b c d = C (gvec4 a b c d) + +vect :: forall n . KnownNat n => [โ„‚] -> C n +vect xs = C (gvect "C" xs) + diff --git a/packages/base/src/Numeric/HMatrix.hs b/packages/base/src/Numeric/HMatrix.hs index ec96bfc..421333a 100644 --- a/packages/base/src/Numeric/HMatrix.hs +++ b/packages/base/src/Numeric/HMatrix.hs @@ -1,201 +1,497 @@ ------------------------------------------------------------------------------ +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FunctionalDependencies #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE EmptyDataDecls #-} +{-# LANGUAGE Rank2Types #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE OverlappingInstances #-} +{-# LANGUAGE TypeFamilies #-} + + {- | Module : Numeric.HMatrix Copyright : (c) Alberto Ruiz 2006-14 License : BSD3 -Maintainer : Alberto Ruiz -Stability : provisional +Stability : experimental + +Experimental interface for real arrays with statically checked dimensions. -} ------------------------------------------------------------------------------ -module Numeric.HMatrix ( - - -- * Basic types and data processing - module Numeric.LinearAlgebra.Data, - - -- * Arithmetic and numeric classes - -- | - -- The standard numeric classes are defined elementwise: - -- - -- >>> vect [1,2,3] * vect [3,0,-2] - -- fromList [3.0,0.0,-6.0] - -- - -- >>> mat 3 [1..9] * ident 3 - -- (3><3) - -- [ 1.0, 0.0, 0.0 - -- , 0.0, 5.0, 0.0 - -- , 0.0, 0.0, 9.0 ] - -- - -- In arithmetic operations single-element vectors and matrices - -- (created from numeric literals or using 'scalar') automatically - -- expand to match the dimensions of the other operand: - -- - -- >>> 5 + 2*ident 3 :: Matrix Double - -- (3><3) - -- [ 7.0, 5.0, 5.0 - -- , 5.0, 7.0, 5.0 - -- , 5.0, 5.0, 7.0 ] - -- - -- >>> mat 3 [1..9] + mat 1 [10,20,30] - -- (3><3) - -- [ 11.0, 12.0, 13.0 - -- , 24.0, 25.0, 26.0 - -- , 37.0, 38.0, 39.0 ] - -- +module Numeric.HMatrix( + -- * Vector + R, + vec2, vec3, vec4, (&), (#), split, headTail, + vector, + linspace, range, dim, + -- * Matrix + L, Sq, build, + row, col, (ยฆ),(โ€”โ€”), splitRows, splitCols, + unrow, uncol, + + eye, + diagR, diag, + blockAt, + matrix, -- * Products - -- ** dot - (<ยท>), - -- ** matrix-vector - (#>), (!#>), - -- ** matrix-matrix - (<>), - -- | The matrix x matrix product is also implemented in the "Data.Monoid" instance, where - -- single-element matrices (created from numeric literals or using 'scalar') - -- are used for scaling. - -- - -- >>> import Data.Monoid as M - -- >>> let m = mat 3 [1..6] - -- >>> m M.<> 2 M.<> diagl[0.5,1,0] - -- (2><3) - -- [ 1.0, 4.0, 0.0 - -- , 4.0, 10.0, 0.0 ] - -- - -- 'mconcat' uses 'optimiseMult' to get the optimal association order. - - - -- ** other - outer, kronecker, cross, - scale, - sumElements, prodElements, - + (<>),(#>),(<ยท>), -- * Linear Systems - (<\>), - linearSolve, - linearSolveLS, - linearSolveSVD, - luSolve, - cholSolve, - cgSolve, - cgSolve', + linSolve, (<\>), + -- * Factorizations + svd, svdTall, svdFlat, Eigen(..), + withNullspace, + -- * Misc + Disp(..), + withVector, withMatrix, + toRows, toColumns, + Sized(..), Diag(..), Sym, sym, + module Numeric.LinearAlgebra.HMatrix +) where - -- * Inverse and pseudoinverse - inv, pinv, pinvTol, - -- * Determinant and rank - rcond, rank, - det, invlndet, +import GHC.TypeLits +import Numeric.LinearAlgebra.HMatrix hiding ( + (<>),(#>),(<ยท>),Konst(..),diag, disp,(ยฆ),(โ€”โ€”),row,col,vector,matrix,linspace,toRows,toColumns, + (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',eigenvalues,eigenvaluesSH,eigenvaluesSH',build) +import qualified Numeric.LinearAlgebra.HMatrix as LA +import Data.Proxy(Proxy) +import Numeric.LinearAlgebra.Static +import Control.Arrow((***)) - -- * Norms - Normed(..), - norm_Frob, norm_nuclear, - -- * Nullspace and range - orth, - nullspace, null1, null1sym, - -- * SVD - svd, - fullSVD, - thinSVD, - compactSVD, - singularValues, - leftSV, rightSV, - -- * Eigensystems - eig, eigSH, eigSH', - eigenvalues, eigenvaluesSH, eigenvaluesSH', - geigSH', - -- * QR - qr, rq, qrRaw, qrgr, +ud1 :: R n -> Vector โ„ +ud1 (R (Dim v)) = v - -- * Cholesky - chol, cholSH, mbCholSH, - -- * Hessenberg - hess, +infixl 4 & +(&) :: forall n . KnownNat n + => R n -> โ„ -> R (n+1) +u & x = u # (konst x :: R 1) - -- * Schur - schur, +infixl 4 # +(#) :: forall n m . (KnownNat n, KnownNat m) + => R n -> R m -> R (n+m) +(R u) # (R v) = R (vconcat u v) - -- * LU - lu, luPacked, - -- * Matrix functions - expm, - sqrtm, - matFunc, - -- * Correlation and convolution - corr, conv, corrMin, corr2, conv2, +vec2 :: โ„ -> โ„ -> R 2 +vec2 a b = R (gvec2 a b) - -- * Random arrays +vec3 :: โ„ -> โ„ -> โ„ -> R 3 +vec3 a b c = R (gvec3 a b c) + + +vec4 :: โ„ -> โ„ -> โ„ -> โ„ -> R 4 +vec4 a b c d = R (gvec4 a b c d) + +vector :: KnownNat n => [โ„] -> R n +vector = fromList + +matrix :: (KnownNat m, KnownNat n) => [โ„] -> L m n +matrix = fromList + +linspace :: forall n . KnownNat n => (โ„,โ„) -> R n +linspace (a,b) = mkR (LA.linspace d (a,b)) + where + d = fromIntegral . natVal $ (undefined :: Proxy n) + +range :: forall n . KnownNat n => R n +range = mkR (LA.linspace d (1,fromIntegral d)) + where + d = fromIntegral . natVal $ (undefined :: Proxy n) + +dim :: forall n . KnownNat n => R n +dim = mkR (scalar d) + where + d = fromIntegral . natVal $ (undefined :: Proxy n) + + +-------------------------------------------------------------------------------- + + +ud2 :: L m n -> Matrix โ„ +ud2 (L (Dim (Dim x))) = x + + +-------------------------------------------------------------------------------- +-------------------------------------------------------------------------------- + +diagR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => โ„ -> R k -> L m n +diagR x v = mkL (asRow (vjoin [scalar x, ev, zeros])) + where + ev = extract v + zeros = LA.konst x (max 0 ((min m' n') - size ev)) + m' = fromIntegral . natVal $ (undefined :: Proxy m) + n' = fromIntegral . natVal $ (undefined :: Proxy n) + +diag :: KnownNat n => R n -> Sq n +diag = diagR 0 + +eye :: KnownNat n => Sq n +eye = diag 1 + +-------------------------------------------------------------------------------- + +blockAt :: forall m n . (KnownNat m, KnownNat n) => โ„ -> Int -> Int -> Matrix Double -> L m n +blockAt x r c a = mkL res + where + z = scalar x + z1 = LA.konst x (r,c) + z2 = LA.konst x (max 0 (m'-(ra+r)), max 0 (n'-(ca+c))) + ra = min (rows a) . max 0 $ m'-r + ca = min (cols a) . max 0 $ n'-c + sa = subMatrix (0,0) (ra, ca) a + m' = fromIntegral . natVal $ (undefined :: Proxy m) + n' = fromIntegral . natVal $ (undefined :: Proxy n) + res = fromBlocks [[z1,z,z],[z,sa,z],[z,z,z2]] + + + + + +-------------------------------------------------------------------------------- + + +row :: R n -> L 1 n +row = mkL . asRow . ud1 + +--col :: R n -> L n 1 +col v = tr . row $ v + +unrow :: L 1 n -> R n +unrow = mkR . head . LA.toRows . ud2 + +--uncol :: L n 1 -> R n +uncol v = unrow . tr $ v + + +infixl 2 โ€”โ€” +(โ€”โ€”) :: (KnownNat r1, KnownNat r2, KnownNat c) => L r1 c -> L r2 c -> L (r1+r2) c +a โ€”โ€” b = mkL (extract a LA.โ€”โ€” extract b) + + +infixl 3 ยฆ +-- (ยฆ) :: (KnownNat r, KnownNat c1, KnownNat c2) => L r c1 -> L r c2 -> L r (c1+c2) +a ยฆ b = tr (tr a โ€”โ€” tr b) + + +type Sq n = L n n +--type CSq n = CL n n + +type GL = (KnownNat n, KnownNat m) => L m n +type GSq = KnownNat n => Sq n + +isKonst :: forall m n . (KnownNat m, KnownNat n) => L m n -> Maybe (โ„,(Int,Int)) +isKonst (unwrap -> x) + | singleM x = Just (x `atIndex` (0,0), (m',n')) + | otherwise = Nothing + where + m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int + n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int - Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, - -- * Misc - meanCov, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, - โ„,โ„‚,iC, - -- * Auxiliary classes - Element, Container, Product, Numeric, LSDiv, - Complexable, RealElement, - RealOf, ComplexOf, SingleOf, DoubleOf, - IndexOf, - Field, --- Normed, - Transposable, - CGState(..), - Testable(..) -) where -import Numeric.LinearAlgebra.Data - -import Numeric.Matrix() -import Numeric.Vector() -import Data.Packed.Numeric hiding ((<>)) -import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed,orth) -import qualified Numeric.LinearAlgebra.Algorithms as A -import Numeric.LinearAlgebra.Util -import Numeric.LinearAlgebra.Random -import Numeric.Sparse((!#>)) -import Numeric.LinearAlgebra.Util.CG - -{- | dense matrix product - ->>> let a = (3><5) [1..] ->>> a -(3><5) - [ 1.0, 2.0, 3.0, 4.0, 5.0 - , 6.0, 7.0, 8.0, 9.0, 10.0 - , 11.0, 12.0, 13.0, 14.0, 15.0 ] - ->>> let b = (5><2) [1,3, 0,2, -1,5, 7,7, 6,0] ->>> b -(5><2) - [ 1.0, 3.0 - , 0.0, 2.0 - , -1.0, 5.0 - , 7.0, 7.0 - , 6.0, 0.0 ] - ->>> a <> b -(3><2) - [ 56.0, 50.0 - , 121.0, 135.0 - , 186.0, 220.0 ] --} -(<>) :: Numeric t => Matrix t -> Matrix t -> Matrix t -(<>) = mXm infixr 8 <> +(<>) :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => L m k -> L k n -> L m n + +(isKonst -> Just (a,(_,k))) <> (isKonst -> Just (b,_)) = konst (a * b * fromIntegral k) + +(isDiag -> Just (0,a,_)) <> (isDiag -> Just (0,b,_)) = diagR 0 (mkR v :: R k) + where + v = a' * b' + n = min (size a) (size b) + a' = subVector 0 n a + b' = subVector 0 n b + +(isDiag -> Just (0,a,_)) <> (extract -> b) = mkL (asColumn a * takeRows (size a) b) + +(extract -> a) <> (isDiag -> Just (0,b,_)) = mkL (takeColumns (size b) a * asRow b) + +a <> b = mkL (extract a LA.<> extract b) + +infixr 8 #> +(#>) :: (KnownNat m, KnownNat n) => L m n -> R n -> R m +(isDiag -> Just (0, w, _)) #> v = mkR (w * subVector 0 (size w) (extract v)) +m #> v = mkR (extract m LA.#> extract v) + + +infixr 8 <ยท> +(<ยท>) :: R n -> R n -> โ„ +(ud1 -> u) <ยท> (ud1 -> v) + | singleV u || singleV v = sumElements (u * v) + | otherwise = udot u v + +-------------------------------------------------------------------------------- + +{- +class Minim (n :: Nat) (m :: Nat) + where + type Mini n m :: Nat + +instance forall (n :: Nat) . Minim n n + where + type Mini n n = n + + +instance forall (n :: Nat) (m :: Nat) . (n <= m+1) => Minim n m + where + type Mini n m = n + +instance forall (n :: Nat) (m :: Nat) . (m <= n+1) => Minim n m + where + type Mini n m = m +-} + +class Diag m d | m -> d + where + takeDiag :: m -> d + + + +instance forall n . (KnownNat n) => Diag (L n n) (R n) + where + takeDiag m = mkR (LA.takeDiag (extract m)) + + +instance forall m n . (KnownNat m, KnownNat n, m <= n+1) => Diag (L m n) (R m) + where + takeDiag m = mkR (LA.takeDiag (extract m)) + + +instance forall m n . (KnownNat m, KnownNat n, n <= m+1) => Diag (L m n) (R n) + where + takeDiag m = mkR (LA.takeDiag (extract m)) + + +-------------------------------------------------------------------------------- + +linSolve :: (KnownNat m, KnownNat n) => L m m -> L m n -> Maybe (L m n) +linSolve (extract -> a) (extract -> b) = fmap mkL (LA.linearSolve a b) + +(<\>) :: (KnownNat m, KnownNat n, KnownNat r) => L m n -> L m r -> L n r +(extract -> a) <\> (extract -> b) = mkL (a LA.<\> b) + +svd :: (KnownNat m, KnownNat n) => L m n -> (L m m, R n, L n n) +svd (extract -> m) = (mkL u, mkR s', mkL v) + where + (u,s,v) = LA.svd m + s' = vjoin [s, z] + z = LA.konst 0 (max 0 (cols m - size s)) + + +svdTall :: (KnownNat m, KnownNat n, n <= m) => L m n -> (L m n, R n, L n n) +svdTall (extract -> m) = (mkL u, mkR s, mkL v) + where + (u,s,v) = LA.thinSVD m + + +svdFlat :: (KnownNat m, KnownNat n, m <= n) => L m n -> (L m m, R m, L n m) +svdFlat (extract -> m) = (mkL u, mkR s, mkL v) + where + (u,s,v) = LA.thinSVD m + +-------------------------------------------------------------------------------- + +class Eigen m l v | m -> l, m -> v + where + eigensystem :: m -> (l,v) + eigenvalues :: m -> l + +newtype Sym n = Sym (Sq n) deriving Show + + +sym :: KnownNat n => Sq n -> Sym n +sym m = Sym $ (m + tr m)/2 + + + +instance KnownNat n => Eigen (Sym n) (R n) (L n n) + where + eigenvalues (Sym (extract -> m)) = mkR . LA.eigenvaluesSH' $ m + eigensystem (Sym (extract -> m)) = (mkR l, mkL v) + where + (l,v) = LA.eigSH' m + +instance KnownNat n => Eigen (Sq n) (C n) (M n n) + where + eigenvalues (extract -> m) = mkC . LA.eigenvalues $ m + eigensystem (extract -> m) = (mkC l, mkM v) + where + (l,v) = LA.eig m + +-------------------------------------------------------------------------------- + + +withNullspace + :: forall m n z . (KnownNat m, KnownNat n) + => L m n + -> (forall k . (KnownNat k) => L n k -> z) + -> z +withNullspace (LA.nullspace . extract -> a) f = + case someNatVal $ fromIntegral $ cols a of + Nothing -> error "static/dynamic mismatch" + Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k) + +-------------------------------------------------------------------------------- + +split :: forall p n . (KnownNat p, KnownNat n, p<=n) => R n -> (R p, R (n-p)) +split (extract -> v) = ( mkR (subVector 0 p' v) , + mkR (subVector p' (size v - p') v) ) + where + p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int + + +headTail :: (KnownNat n, 1<=n) => R n -> (โ„, R (n-1)) +headTail = ((!0) . extract *** id) . split + + +splitRows :: forall p m n . (KnownNat p, KnownNat m, KnownNat n, p<=m) => L m n -> (L p n, L (m-p) n) +splitRows (extract -> x) = ( mkL (takeRows p' x) , + mkL (dropRows p' x) ) + where + p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int + +splitCols :: forall p m n. (KnownNat p, KnownNat m, KnownNat n, KnownNat (n-p), p<=n) => L m n -> (L m p, L m (n-p)) +splitCols = (tr *** tr) . splitRows . tr + + +toRows :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R n] +toRows (LA.toRows . extract -> vs) = map mkR vs + + +toColumns :: forall m n . (KnownNat m, KnownNat n) => L m n -> [R m] +toColumns (LA.toColumns . extract -> vs) = map mkR vs + + +splittest + = do + let v = range :: R 7 + a = snd (split v) :: R 4 + print $ a + print $ snd . headTail . snd . headTail $ v + print $ first (vec3 1 2 3) + print $ second (vec3 1 2 3) + print $ third (vec3 1 2 3) + print $ (snd $ splitRows eye :: L 4 6) + where + first v = fst . headTail $ v + second v = first . snd . headTail $ v + third v = first . snd . headTail . snd . headTail $ v + +-------------------------------------------------------------------------------- + +build + :: forall m n . (KnownNat n, KnownNat m) + => (โ„ -> โ„ -> โ„) + -> L m n +build f = mkL $ LA.build (m',n') f + where + m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int + n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int --- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. -linearSolve m b = A.mbLinearSolve m b +-------------------------------------------------------------------------------- --- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'. -nullspace m = nullspaceSVD (Left (1*eps)) m (rightSV m) +withVector + :: forall z + . Vector โ„ + -> (forall n . (KnownNat n) => R n -> z) + -> z +withVector v f = + case someNatVal $ fromIntegral $ size v of + Nothing -> error "static/dynamic mismatch" + Just (SomeNat (_ :: Proxy m)) -> f (mkR v :: R m) + + +withMatrix + :: forall z + . Matrix โ„ + -> (forall m n . (KnownNat m, KnownNat n) => L m n -> z) + -> z +withMatrix a f = + case someNatVal $ fromIntegral $ rows a of + Nothing -> error "static/dynamic mismatch" + Just (SomeNat (_ :: Proxy m)) -> + case someNatVal $ fromIntegral $ cols a of + Nothing -> error "static/dynamic mismatch" + Just (SomeNat (_ :: Proxy n)) -> + f (mkL a :: L m n) + + +-------------------------------------------------------------------------------- + +test :: (Bool, IO ()) +test = (ok,info) + where + ok = extract (eye :: Sq 5) == ident 5 + && unwrap (mTm sm :: Sq 3) == tr ((3><3)[1..]) LA.<> (3><3)[1..] + && unwrap (tm :: L 3 5) == LA.matrix 5 [1..15] + && thingS == thingD + && precS == precD + && withVector (LA.vector [1..15]) sumV == sumElements (LA.fromList [1..15]) + + info = do + print $ u + print $ v + print (eye :: Sq 3) + print $ ((u & 5) + 1) <ยท> v + print (tm :: L 2 5) + print (tm <> sm :: L 2 3) + print thingS + print thingD + print precS + print precD + print $ withVector (LA.vector [1..15]) sumV + splittest + + sumV w = w <ยท> konst 1 + + u = vec2 3 5 + + ๐•ง x = vector [x] :: R 1 + + v = ๐•ง 2 & 4 & 7 + +-- mTm :: L n m -> Sq m + mTm a = tr a <> a + + tm :: GL + tm = lmat 0 [1..] + + lmat :: forall m n . (KnownNat m, KnownNat n) => โ„ -> [โ„] -> L m n + lmat z xs = mkL . reshape n' . LA.fromList . take (m'*n') $ xs ++ repeat z + where + m' = fromIntegral . natVal $ (undefined :: Proxy m) + n' = fromIntegral . natVal $ (undefined :: Proxy n) + + sm :: GSq + sm = lmat 0 [1..] + + thingS = (u & 1) <ยท> tr q #> q #> v + where + q = tm :: L 10 3 + + 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 (size (ud1 u) +1, size (ud1 v)) LA.#> ud1 v + + +instance (KnownNat n', KnownNat m') => Testable (L n' m') + where + checkT _ = test --- | return an orthonormal basis of the range space of a matrix. See also 'orthSVD'. -orth m = orthSVD (Left (1*eps)) m (leftSV m) diff --git a/packages/base/src/Numeric/LinearAlgebra/Complex.hs b/packages/base/src/Numeric/LinearAlgebra/Complex.hs deleted file mode 100644 index 17bc397..0000000 --- a/packages/base/src/Numeric/LinearAlgebra/Complex.hs +++ /dev/null @@ -1,80 +0,0 @@ -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE KindSignatures #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FunctionalDependencies #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE EmptyDataDecls #-} -{-# LANGUAGE Rank2Types #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE ViewPatterns #-} -{-# LANGUAGE GADTs #-} - - -{- | -Module : Numeric.LinearAlgebra.Complex -Copyright : (c) Alberto Ruiz 2006-14 -License : BSD3 -Stability : experimental - --} - -module Numeric.LinearAlgebra.Complex( - C, - vec2, vec3, vec4, (&), (#), - vect, - R -) where - -import GHC.TypeLits -import Numeric.HMatrix hiding ( - (<>),(#>),(<ยท>),Konst(..),diag, disp,(ยฆ),(โ€”โ€”),row,col,vect,mat,linspace) -import qualified Numeric.HMatrix as LA -import Data.Proxy(Proxy) -import Numeric.LinearAlgebra.Static - - - -instance forall n . KnownNat n => Show (C n) - where - show (ud1 -> v) - | size v == 1 = "("++show (v!0)++" :: C "++show d++")" - | otherwise = "(vect"++ drop 8 (show v)++" :: C "++show d++")" - where - d = fromIntegral . natVal $ (undefined :: Proxy n) :: Int - - -ud1 :: C n -> Vector โ„‚ -ud1 (C (Dim v)) = v - -mkC :: Vector โ„‚ -> C n -mkC = C . Dim - - -infixl 4 & -(&) :: forall n . KnownNat n - => C n -> โ„‚ -> C (n+1) -u & x = u # (mkC (LA.scalar x) :: C 1) - -infixl 4 # -(#) :: forall n m . (KnownNat n, KnownNat m) - => C n -> C m -> C (n+m) -(C u) # (C v) = C (vconcat u v) - - - -vec2 :: โ„‚ -> โ„‚ -> C 2 -vec2 a b = C (gvec2 a b) - -vec3 :: โ„‚ -> โ„‚ -> โ„‚ -> C 3 -vec3 a b c = C (gvec3 a b c) - - -vec4 :: โ„‚ -> โ„‚ -> โ„‚ -> โ„‚ -> C 4 -vec4 a b c d = C (gvec4 a b c d) - -vect :: forall n . KnownNat n => [โ„‚] -> C n -vect xs = C (gvect "C" xs) - diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs index 20f3b81..4e4868a 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Data.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs @@ -16,11 +16,11 @@ module Numeric.LinearAlgebra.Data( -- * Vector -- | 1D arrays are storable vectors from the vector package. - vect, (|>), + vector, (|>), -- * Matrix - mat, (><), tr, + matrix, (><), tr, -- * Indexing diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs new file mode 100644 index 0000000..54ddd68 --- /dev/null +++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs @@ -0,0 +1,201 @@ +----------------------------------------------------------------------------- +{- | +Module : Numeric.LinearAlgebra.HMatrix +Copyright : (c) Alberto Ruiz 2006-14 +License : BSD3 +Maintainer : Alberto Ruiz +Stability : provisional + +-} +----------------------------------------------------------------------------- +module Numeric.LinearAlgebra.HMatrix ( + + -- * Basic types and data processing + module Numeric.LinearAlgebra.Data, + + -- * Arithmetic and numeric classes + -- | + -- The standard numeric classes are defined elementwise: + -- + -- >>> vector [1,2,3] * vector [3,0,-2] + -- fromList [3.0,0.0,-6.0] + -- + -- >>> matrix 3 [1..9] * ident 3 + -- (3><3) + -- [ 1.0, 0.0, 0.0 + -- , 0.0, 5.0, 0.0 + -- , 0.0, 0.0, 9.0 ] + -- + -- In arithmetic operations single-element vectors and matrices + -- (created from numeric literals or using 'scalar') automatically + -- expand to match the dimensions of the other operand: + -- + -- >>> 5 + 2*ident 3 :: Matrix Double + -- (3><3) + -- [ 7.0, 5.0, 5.0 + -- , 5.0, 7.0, 5.0 + -- , 5.0, 5.0, 7.0 ] + -- + -- >>> matrix 3 [1..9] + matrix 1 [10,20,30] + -- (3><3) + -- [ 11.0, 12.0, 13.0 + -- , 24.0, 25.0, 26.0 + -- , 37.0, 38.0, 39.0 ] + -- + + -- * Products + -- ** dot + (<ยท>), + -- ** matrix-vector + (#>), (!#>), + -- ** matrix-matrix + (<>), + -- | The matrix x matrix product is also implemented in the "Data.Monoid" instance, where + -- single-element matrices (created from numeric literals or using 'scalar') + -- are used for scaling. + -- + -- >>> import Data.Monoid as M + -- >>> let m = matrix 3 [1..6] + -- >>> m M.<> 2 M.<> diagl[0.5,1,0] + -- (2><3) + -- [ 1.0, 4.0, 0.0 + -- , 4.0, 10.0, 0.0 ] + -- + -- 'mconcat' uses 'optimiseMult' to get the optimal association order. + + + -- ** other + outer, kronecker, cross, + scale, + sumElements, prodElements, + + -- * Linear Systems + (<\>), + linearSolve, + linearSolveLS, + linearSolveSVD, + luSolve, + cholSolve, + cgSolve, + cgSolve', + + -- * Inverse and pseudoinverse + inv, pinv, pinvTol, + + -- * Determinant and rank + rcond, rank, + det, invlndet, + + -- * Norms + Normed(..), + norm_Frob, norm_nuclear, + + -- * Nullspace and range + orth, + nullspace, null1, null1sym, + + -- * SVD + svd, + fullSVD, + thinSVD, + compactSVD, + singularValues, + leftSV, rightSV, + + -- * Eigensystems + eig, eigSH, eigSH', + eigenvalues, eigenvaluesSH, eigenvaluesSH', + geigSH', + + -- * QR + qr, rq, qrRaw, qrgr, + + -- * Cholesky + chol, cholSH, mbCholSH, + + -- * Hessenberg + hess, + + -- * Schur + schur, + + -- * LU + lu, luPacked, + + -- * Matrix functions + expm, + sqrtm, + matFunc, + + -- * Correlation and convolution + corr, conv, corrMin, corr2, conv2, + + -- * Random arrays + + Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, + + -- * Misc + meanCov, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, + โ„,โ„‚,iC, + -- * Auxiliary classes + Element, Container, Product, Numeric, LSDiv, + Complexable, RealElement, + RealOf, ComplexOf, SingleOf, DoubleOf, + IndexOf, + Field, +-- Normed, + Transposable, + CGState(..), + Testable(..) +) where + +import Numeric.LinearAlgebra.Data + +import Numeric.Matrix() +import Numeric.Vector() +import Data.Packed.Numeric hiding ((<>)) +import Numeric.LinearAlgebra.Algorithms hiding (linearSolve,Normed,orth) +import qualified Numeric.LinearAlgebra.Algorithms as A +import Numeric.LinearAlgebra.Util +import Numeric.LinearAlgebra.Random +import Numeric.Sparse((!#>)) +import Numeric.LinearAlgebra.Util.CG + +{- | dense matrix product + +>>> let a = (3><5) [1..] +>>> a +(3><5) + [ 1.0, 2.0, 3.0, 4.0, 5.0 + , 6.0, 7.0, 8.0, 9.0, 10.0 + , 11.0, 12.0, 13.0, 14.0, 15.0 ] + +>>> let b = (5><2) [1,3, 0,2, -1,5, 7,7, 6,0] +>>> b +(5><2) + [ 1.0, 3.0 + , 0.0, 2.0 + , -1.0, 5.0 + , 7.0, 7.0 + , 6.0, 0.0 ] + +>>> a <> b +(3><2) + [ 56.0, 50.0 + , 121.0, 135.0 + , 186.0, 220.0 ] + +-} +(<>) :: Numeric t => Matrix t -> Matrix t -> Matrix t +(<>) = mXm +infixr 8 <> + +-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. +linearSolve m b = A.mbLinearSolve m b + +-- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'. +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) + diff --git a/packages/base/src/Numeric/LinearAlgebra/Real.hs b/packages/base/src/Numeric/LinearAlgebra/Real.hs deleted file mode 100644 index 97c462e..0000000 --- a/packages/base/src/Numeric/LinearAlgebra/Real.hs +++ /dev/null @@ -1,480 +0,0 @@ -{-# LANGUAGE DataKinds #-} -{-# LANGUAGE KindSignatures #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE FunctionalDependencies #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE EmptyDataDecls #-} -{-# LANGUAGE Rank2Types #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE ViewPatterns #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE OverlappingInstances #-} -{-# LANGUAGE TypeFamilies #-} - - -{- | -Module : Numeric.LinearAlgebra.Real -Copyright : (c) Alberto Ruiz 2006-14 -License : BSD3 -Stability : experimental - -Experimental interface for real arrays with statically checked dimensions. - --} - -module Numeric.LinearAlgebra.Real( - -- * Vector - R, C, - vec2, vec3, vec4, (&), (#), split, headTail, - vector, - linspace, range, dim, - -- * Matrix - L, Sq, M, def, - row, col, (ยฆ),(โ€”โ€”), splitRows, splitCols, - unrow, uncol, - - eye, - diagR, diag, - blockAt, - matrix, - -- * Products - (<>),(#>),(<ยท>), - -- * Linear Systems - linSolve, (<\>), - -- * Factorizations - svd, svdTall, svdFlat, Eigen(..), - -- * Pretty printing - Disp(..), - -- * Misc - withVector, withMatrix, - Sized(..), Diag(..), Sym, sym, Her, her, ๐‘–, - module Numeric.HMatrix -) where - - -import GHC.TypeLits -import Numeric.HMatrix hiding ( - (<>),(#>),(<ยท>),Konst(..),diag, disp,(ยฆ),(โ€”โ€”),row,col,vect,mat,linspace, - (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',eigenvalues,eigenvaluesSH,eigenvaluesSH',build) -import qualified Numeric.HMatrix as LA -import Data.Proxy(Proxy) -import Numeric.LinearAlgebra.Static -import Control.Arrow((***)) - - -๐‘– :: Sized โ„‚ s c => s -๐‘– = konst iC - - - - - -ud1 :: R n -> Vector โ„ -ud1 (R (Dim v)) = v - - -infixl 4 & -(&) :: forall n . KnownNat n - => R n -> โ„ -> R (n+1) -u & x = u # (konst x :: R 1) - -infixl 4 # -(#) :: forall n m . (KnownNat n, KnownNat m) - => R n -> R m -> R (n+m) -(R u) # (R v) = R (vconcat u v) - - - -vec2 :: โ„ -> โ„ -> R 2 -vec2 a b = R (gvec2 a b) - -vec3 :: โ„ -> โ„ -> โ„ -> R 3 -vec3 a b c = R (gvec3 a b c) - - -vec4 :: โ„ -> โ„ -> โ„ -> โ„ -> R 4 -vec4 a b c d = R (gvec4 a b c d) - -vector :: KnownNat n => [โ„] -> R n -vector = fromList - -matrix :: (KnownNat m, KnownNat n) => [โ„] -> L m n -matrix = fromList - -linspace :: forall n . KnownNat n => (โ„,โ„) -> R n -linspace (a,b) = mkR (LA.linspace d (a,b)) - where - d = fromIntegral . natVal $ (undefined :: Proxy n) - -range :: forall n . KnownNat n => R n -range = mkR (LA.linspace d (1,fromIntegral d)) - where - d = fromIntegral . natVal $ (undefined :: Proxy n) - -dim :: forall n . KnownNat n => R n -dim = mkR (scalar d) - where - d = fromIntegral . natVal $ (undefined :: Proxy n) - - --------------------------------------------------------------------------------- - - -ud2 :: L m n -> Matrix โ„ -ud2 (L (Dim (Dim x))) = x - - --------------------------------------------------------------------------------- --------------------------------------------------------------------------------- - -diagR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => โ„ -> R k -> L m n -diagR x v = mkL (asRow (vjoin [scalar x, ev, zeros])) - where - ev = extract v - zeros = LA.konst x (max 0 ((min m' n') - size ev)) - m' = fromIntegral . natVal $ (undefined :: Proxy m) - n' = fromIntegral . natVal $ (undefined :: Proxy n) - -diag :: KnownNat n => R n -> Sq n -diag = diagR 0 - -eye :: KnownNat n => Sq n -eye = diag 1 - --------------------------------------------------------------------------------- - -blockAt :: forall m n . (KnownNat m, KnownNat n) => โ„ -> Int -> Int -> Matrix Double -> L m n -blockAt x r c a = mkL res - where - z = scalar x - z1 = LA.konst x (r,c) - z2 = LA.konst x (max 0 (m'-(ra+r)), max 0 (n'-(ca+c))) - ra = min (rows a) . max 0 $ m'-r - ca = min (cols a) . max 0 $ n'-c - sa = subMatrix (0,0) (ra, ca) a - m' = fromIntegral . natVal $ (undefined :: Proxy m) - n' = fromIntegral . natVal $ (undefined :: Proxy n) - res = fromBlocks [[z1,z,z],[z,sa,z],[z,z,z2]] - - - - - --------------------------------------------------------------------------------- - - -row :: R n -> L 1 n -row = mkL . asRow . ud1 - ---col :: R n -> L n 1 -col v = tr . row $ v - -unrow :: L 1 n -> R n -unrow = mkR . head . toRows . ud2 - ---uncol :: L n 1 -> R n -uncol v = unrow . tr $ v - - -infixl 2 โ€”โ€” -(โ€”โ€”) :: (KnownNat r1, KnownNat r2, KnownNat c) => L r1 c -> L r2 c -> L (r1+r2) c -a โ€”โ€” b = mkL (extract a LA.โ€”โ€” extract b) - - -infixl 3 ยฆ --- (ยฆ) :: (KnownNat r, KnownNat c1, KnownNat c2) => L r c1 -> L r c2 -> L r (c1+c2) -a ยฆ b = tr (tr a โ€”โ€” tr b) - - -type Sq n = L n n ---type CSq n = CL n n - -type GL = (KnownNat n, KnownNat m) => L m n -type GSq = KnownNat n => Sq n - -isKonst :: forall m n . (KnownNat m, KnownNat n) => L m n -> Maybe (โ„,(Int,Int)) -isKonst (unwrap -> x) - | singleM x = Just (x `atIndex` (0,0), (m',n')) - | otherwise = Nothing - where - m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int - n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int - - - - -infixr 8 <> -(<>) :: forall m k n. (KnownNat m, KnownNat k, KnownNat n) => L m k -> L k n -> L m n - -(isKonst -> Just (a,(_,k))) <> (isKonst -> Just (b,_)) = konst (a * b * fromIntegral k) - -(isDiag -> Just (0,a,_)) <> (isDiag -> Just (0,b,_)) = diagR 0 (mkR v :: R k) - where - v = a' * b' - n = min (size a) (size b) - a' = subVector 0 n a - b' = subVector 0 n b - -(isDiag -> Just (0,a,_)) <> (extract -> b) = mkL (asColumn a * takeRows (size a) b) - -(extract -> a) <> (isDiag -> Just (0,b,_)) = mkL (takeColumns (size b) a * asRow b) - -a <> b = mkL (extract a LA.<> extract b) - -infixr 8 #> -(#>) :: (KnownNat m, KnownNat n) => L m n -> R n -> R m -(isDiag -> Just (0, w, _)) #> v = mkR (w * subVector 0 (size w) (extract v)) -m #> v = mkR (extract m LA.#> extract v) - - -infixr 8 <ยท> -(<ยท>) :: R n -> R n -> โ„ -(ud1 -> u) <ยท> (ud1 -> v) - | singleV u || singleV v = sumElements (u * v) - | otherwise = udot u v - --------------------------------------------------------------------------------- - -{- -class Minim (n :: Nat) (m :: Nat) - where - type Mini n m :: Nat - -instance forall (n :: Nat) . Minim n n - where - type Mini n n = n - - -instance forall (n :: Nat) (m :: Nat) . (n <= m+1) => Minim n m - where - type Mini n m = n - -instance forall (n :: Nat) (m :: Nat) . (m <= n+1) => Minim n m - where - type Mini n m = m --} - -class Diag m d | m -> d - where - takeDiag :: m -> d - - - -instance forall n . (KnownNat n) => Diag (L n n) (R n) - where - takeDiag m = mkR (LA.takeDiag (extract m)) - - -instance forall m n . (KnownNat m, KnownNat n, m <= n+1) => Diag (L m n) (R m) - where - takeDiag m = mkR (LA.takeDiag (extract m)) - - -instance forall m n . (KnownNat m, KnownNat n, n <= m+1) => Diag (L m n) (R n) - where - takeDiag m = mkR (LA.takeDiag (extract m)) - - --------------------------------------------------------------------------------- - -linSolve :: (KnownNat m, KnownNat n) => L m m -> L m n -> Maybe (L m n) -linSolve (extract -> a) (extract -> b) = fmap mkL (LA.linearSolve a b) - -(<\>) :: (KnownNat m, KnownNat n, KnownNat r) => L m n -> L m r -> L n r -(extract -> a) <\> (extract -> b) = mkL (a LA.<\> b) - -svd :: (KnownNat m, KnownNat n) => L m n -> (L m m, R n, L n n) -svd (extract -> m) = (mkL u, mkR s', mkL v) - where - (u,s,v) = LA.svd m - s' = vjoin [s, z] - z = LA.konst 0 (max 0 (cols m - size s)) - - -svdTall :: (KnownNat m, KnownNat n, n <= m) => L m n -> (L m n, R n, L n n) -svdTall (extract -> m) = (mkL u, mkR s, mkL v) - where - (u,s,v) = LA.thinSVD m - - -svdFlat :: (KnownNat m, KnownNat n, m <= n) => L m n -> (L m m, R m, L m n) -svdFlat (extract -> m) = (mkL u, mkR s, mkL v) - where - (u,s,v) = LA.thinSVD m - --------------------------------------------------------------------------------- - -class Eigen m l v | m -> l, m -> v - where - eigensystem :: m -> (l,v) - eigenvalues :: m -> l - -newtype Sym n = Sym (Sq n) deriving Show - -newtype Her n = Her (M n n) - -sym :: KnownNat n => Sq n -> Sym n -sym m = Sym $ (m + tr m)/2 - -her :: KnownNat n => M n n -> Her n -her m = Her $ (m + tr m)/2 - -instance KnownNat n => Eigen (Sym n) (R n) (L n n) - where - eigenvalues (Sym (extract -> m)) = mkR . LA.eigenvaluesSH' $ m - eigensystem (Sym (extract -> m)) = (mkR l, mkL v) - where - (l,v) = LA.eigSH' m - -instance KnownNat n => Eigen (Sq n) (C n) (M n n) - where - eigenvalues (extract -> m) = mkC . LA.eigenvalues $ m - eigensystem (extract -> m) = (mkC l, mkM v) - where - (l,v) = LA.eig m - --------------------------------------------------------------------------------- - -split :: forall p n . (KnownNat p, KnownNat n, p<=n) => R n -> (R p, R (n-p)) -split (extract -> v) = ( mkR (subVector 0 p' v) , - mkR (subVector p' (size v - p') v) ) - where - p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int - - -headTail :: (KnownNat n, 1<=n) => R n -> (โ„, R (n-1)) -headTail = ((!0) . extract *** id) . split - - -splitRows :: forall p m n. (KnownNat p, KnownNat m, KnownNat n, p<=m) => L m n -> (L p n, L (m-p) n) -splitRows (extract -> x) = ( mkL (takeRows p' x) , - mkL (dropRows p' x) ) - where - p' = fromIntegral . natVal $ (undefined :: Proxy p) :: Int - -splitCols :: forall p m n. (KnownNat p, KnownNat m, KnownNat n, KnownNat (n-p), p<=n) => L m n -> (L m p, L m (n-p)) -splitCols = (tr *** tr) . splitRows . tr - - -splittest - = do - let v = range :: R 7 - a = snd (split v) :: R 4 - print $ a - print $ snd . headTail . snd . headTail $ v - print $ first (vec3 1 2 3) - print $ second (vec3 1 2 3) - print $ third (vec3 1 2 3) - print $ (snd $ splitRows eye :: L 4 6) - where - first v = fst . headTail $ v - second v = first . snd . headTail $ v - third v = first . snd . headTail . snd . headTail $ v - --------------------------------------------------------------------------------- - -def - :: forall m n . (KnownNat n, KnownNat m) - => (โ„ -> โ„ -> โ„) - -> L m n -def f = mkL $ LA.build (m',n') f - where - m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int - n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int - --------------------------------------------------------------------------------- - -withVector - :: forall z - . Vector โ„ - -> (forall n . (KnownNat n) => R n -> z) - -> z -withVector v f = - case someNatVal $ fromIntegral $ size v of - Nothing -> error "static/dynamic mismatch" - Just (SomeNat (_ :: Proxy m)) -> f (mkR v :: R m) - - -withMatrix - :: forall z - . Matrix โ„ - -> (forall m n . (KnownNat m, KnownNat n) => L m n -> z) - -> z -withMatrix a f = - case someNatVal $ fromIntegral $ rows a of - Nothing -> error "static/dynamic mismatch" - Just (SomeNat (_ :: Proxy m)) -> - case someNatVal $ fromIntegral $ cols a of - Nothing -> error "static/dynamic mismatch" - Just (SomeNat (_ :: Proxy n)) -> - f (mkL a :: L n m) - --------------------------------------------------------------------------------- - -test :: (Bool, IO ()) -test = (ok,info) - where - ok = extract (eye :: Sq 5) == ident 5 - && unwrap (mTm sm :: Sq 3) == tr ((3><3)[1..]) LA.<> (3><3)[1..] - && unwrap (tm :: L 3 5) == LA.mat 5 [1..15] - && thingS == thingD - && precS == precD - && withVector (LA.vect [1..15]) sumV == sumElements (LA.fromList [1..15]) - - info = do - print $ u - print $ v - print (eye :: Sq 3) - print $ ((u & 5) + 1) <ยท> v - print (tm :: L 2 5) - print (tm <> sm :: L 2 3) - print thingS - print thingD - print precS - print precD - print $ withVector (LA.vect [1..15]) sumV - splittest - - sumV w = w <ยท> konst 1 - - u = vec2 3 5 - - ๐•ง x = vector [x] :: R 1 - - v = ๐•ง 2 & 4 & 7 - --- mTm :: L n m -> Sq m - mTm a = tr a <> a - - tm :: GL - tm = lmat 0 [1..] - - lmat :: forall m n . (KnownNat m, KnownNat n) => โ„ -> [โ„] -> L m n - lmat z xs = mkL . reshape n' . LA.fromList . take (m'*n') $ xs ++ repeat z - where - m' = fromIntegral . natVal $ (undefined :: Proxy m) - n' = fromIntegral . natVal $ (undefined :: Proxy n) - - sm :: GSq - sm = lmat 0 [1..] - - thingS = (u & 1) <ยท> tr q #> q #> v - where - q = tm :: L 10 3 - - thingD = vjoin [ud1 u, 1] LA.<ยท> tr m LA.#> m LA.#> ud1 v - where - m = LA.mat 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 (size (ud1 u) +1, size (ud1 v)) LA.#> ud1 v - - -instance (KnownNat n', KnownNat m') => Testable (L n' m') - where - checkT _ = test - - diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs index 2647f20..daf8d80 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs @@ -25,7 +25,7 @@ module Numeric.LinearAlgebra.Static where import GHC.TypeLits -import Numeric.HMatrix as LA +import Numeric.LinearAlgebra.HMatrix as LA import Data.Packed as D import Data.Packed.ST import Data.Proxy(Proxy) diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs index 4824af4..0ac4634 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs @@ -19,7 +19,7 @@ Stability : provisional module Numeric.LinearAlgebra.Util( -- * Convenience functions - vect, mat, + vector, matrix, disp, zeros, ones, diagl, @@ -79,27 +79,27 @@ iC = 0:+1 {- | create a real vector ->>> vect [1..5] +>>> vector [1..5] fromList [1.0,2.0,3.0,4.0,5.0] -} -vect :: [โ„] -> Vector โ„ -vect = fromList +vector :: [โ„] -> Vector โ„ +vector = fromList {- | create a real matrix ->>> mat 5 [1..15] +>>> matrix 5 [1..15] (3><5) [ 1.0, 2.0, 3.0, 4.0, 5.0 , 6.0, 7.0, 8.0, 9.0, 10.0 , 11.0, 12.0, 13.0, 14.0, 15.0 ] -} -mat +matrix :: Int -- ^ columns -> [โ„] -- ^ elements -> Matrix โ„ -mat c = reshape c . fromList +matrix c = reshape c . fromList {- | print a real matrix with given number of digits after the decimal point diff --git a/packages/base/src/Numeric/Sparse.hs b/packages/base/src/Numeric/Sparse.hs index f495e3a..f1516ec 100644 --- a/packages/base/src/Numeric/Sparse.hs +++ b/packages/base/src/Numeric/Sparse.hs @@ -168,7 +168,7 @@ gmXv Dense{..} v {- | general matrix - vector product >>> let m = mkSparse [((0,999),1.0),((1,1999),2.0)] ->>> m !#> vect[1..2000] +>>> m !#> vector [1..2000] fromList [1000.0,4000.0] -} -- cgit v1.2.3