summaryrefslogtreecommitdiff
path: root/packages/hmatrix
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-05-16 20:10:57 +0200
committerAlberto Ruiz <aruiz@um.es>2014-05-16 20:10:57 +0200
commitd4d9082a8d7d3eed6cb5f188fc3b476847dcac27 (patch)
tree0a19bc18be429454fd59c16fd46e29e1e7bae722 /packages/hmatrix
parent51f4cc7b4b301142b8df73568ffaa448f9e6dd50 (diff)
GSL.LinearAlgebra reorganized
Diffstat (limited to 'packages/hmatrix')
-rw-r--r--packages/hmatrix/hmatrix.cabal16
-rw-r--r--packages/hmatrix/src/Numeric/Container.hs241
-rw-r--r--packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs (renamed from packages/hmatrix/src/Numeric/GSL/Vector.hs)8
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra.hs33
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/GSL.hs136
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/IO.hs (renamed from packages/hmatrix/src/Numeric/IO.hs)22
-rw-r--r--packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs (renamed from packages/hmatrix/src/Numeric/Random.hs)6
7 files changed, 190 insertions, 272 deletions
diff --git a/packages/hmatrix/hmatrix.cabal b/packages/hmatrix/hmatrix.cabal
index c79da8a..a509df4 100644
--- a/packages/hmatrix/hmatrix.cabal
+++ b/packages/hmatrix/hmatrix.cabal
@@ -90,15 +90,13 @@ library
90 Numeric.GSL.Fitting, 90 Numeric.GSL.Fitting,
91 Numeric.GSL.ODE, 91 Numeric.GSL.ODE,
92 Numeric.GSL, 92 Numeric.GSL,
93 Graphics.Plot, 93 Numeric.LinearAlgebra,
94 Numeric.HMatrix, 94 Graphics.Plot
95 Numeric.HMatrix.Data, 95 other-modules: Numeric.GSL.Internal,
96 Numeric.HMatrix.Devel 96 Numeric.GSL.LinearAlgebra,
97 other-modules: Numeric.Random, 97 Numeric.LinearAlgebra.IO,
98 Numeric.GSL.Internal, 98 Numeric.LinearAlgebra.Random
99 Numeric.GSL.Vector, 99
100 Numeric.IO
101
102 100
103 C-sources: src/Numeric/GSL/gsl-aux.c 101 C-sources: src/Numeric/GSL/gsl-aux.c
104 102
diff --git a/packages/hmatrix/src/Numeric/Container.hs b/packages/hmatrix/src/Numeric/Container.hs
deleted file mode 100644
index 146780d..0000000
--- a/packages/hmatrix/src/Numeric/Container.hs
+++ /dev/null
@@ -1,241 +0,0 @@
1{-# LANGUAGE TypeFamilies #-}
2{-# LANGUAGE FlexibleContexts #-}
3{-# LANGUAGE FlexibleInstances #-}
4{-# LANGUAGE MultiParamTypeClasses #-}
5{-# LANGUAGE FunctionalDependencies #-}
6{-# LANGUAGE UndecidableInstances #-}
7
8-----------------------------------------------------------------------------
9-- |
10-- Module : Numeric.Container
11-- Copyright : (c) Alberto Ruiz 2010-14
12-- License : BSD3
13-- Maintainer : Alberto Ruiz
14-- Stability : provisional
15-- Portability : portable
16--
17-- Basic numeric operations on 'Vector' and 'Matrix', including conversion routines.
18--
19-- The 'Container' class is used to define optimized generic functions which work
20-- on 'Vector' and 'Matrix' with real or complex elements.
21--
22-- Some of these functions are also available in the instances of the standard
23-- numeric Haskell classes provided by "Numeric.LinearAlgebra".
24--
25-----------------------------------------------------------------------------
26{-# OPTIONS_HADDOCK hide #-}
27
28module Numeric.Container (
29 -- * Basic functions
30 module Data.Packed,
31 konst, build,
32 linspace,
33 diag, ident,
34 ctrans,
35 -- * Generic operations
36 Container(..),
37 -- * Matrix product
38 Product(..), udot, dot, (◇),
39 Mul(..),
40 Contraction(..),
41 optimiseMult,
42 mXm,mXv,vXm,LSDiv(..),
43 outer, kronecker,
44 -- * Element conversion
45 Convert(..),
46 Complexable(),
47 RealElement(),
48
49 RealOf, ComplexOf, SingleOf, DoubleOf,
50
51 IndexOf,
52 module Data.Complex
53) where
54
55import Data.Packed hiding (stepD, stepF, condD, condF, conjugateC, conjugateQ)
56import Data.Packed.Numeric
57import Data.Complex
58import Numeric.LinearAlgebra.Algorithms(Field,linearSolveSVD)
59import Data.Monoid(Monoid(mconcat))
60
61------------------------------------------------------------------
62
63{- | Creates a real vector containing a range of values:
64
65>>> linspace 5 (-3,7::Double)
66fromList [-3.0,-0.5,2.0,4.5,7.0]@
67
68>>> linspace 5 (8,2+i) :: Vector (Complex Double)
69fromList [8.0 :+ 0.0,6.5 :+ 0.25,5.0 :+ 0.5,3.5 :+ 0.75,2.0 :+ 1.0]
70
71Logarithmic spacing can be defined as follows:
72
73@logspace n (a,b) = 10 ** linspace n (a,b)@
74-}
75linspace :: (Container Vector e) => Int -> (e, e) -> Vector e
76linspace 0 (a,b) = fromList[(a+b)/2]
77linspace n (a,b) = addConstant a $ scale s $ fromList $ map fromIntegral [0 .. n-1]
78 where s = (b-a)/fromIntegral (n-1)
79
80--------------------------------------------------------
81
82class Contraction a b c | a b -> c
83 where
84 infixl 7 <.>
85 {- | Matrix product, matrix vector product, and dot product
86
87Examples:
88
89>>> let a = (3><4) [1..] :: Matrix Double
90>>> let v = fromList [1,0,2,-1] :: Vector Double
91>>> let u = fromList [1,2,3] :: Vector Double
92
93>>> a
94(3><4)
95 [ 1.0, 2.0, 3.0, 4.0
96 , 5.0, 6.0, 7.0, 8.0
97 , 9.0, 10.0, 11.0, 12.0 ]
98
99matrix × matrix:
100
101>>> disp 2 (a <.> trans a)
1023x3
103 30 70 110
104 70 174 278
105110 278 446
106
107matrix × vector:
108
109>>> a <.> v
110fromList [3.0,11.0,19.0]
111
112dot product:
113
114>>> u <.> fromList[3,2,1::Double]
11510
116
117For complex vectors the first argument is conjugated:
118
119>>> fromList [1,i] <.> fromList[2*i+1,3]
1201.0 :+ (-1.0)
121
122>>> fromList [1,i,1-i] <.> complex a
123fromList [10.0 :+ 4.0,12.0 :+ 4.0,14.0 :+ 4.0,16.0 :+ 4.0]
124
125-}
126 (<.>) :: a -> b -> c
127
128
129instance (Product t, Container Vector t) => Contraction (Vector t) (Vector t) t where
130 u <.> v = conj u `udot` v
131
132instance Product t => Contraction (Matrix t) (Vector t) (Vector t) where
133 (<.>) = mXv
134
135instance (Container Vector t, Product t) => Contraction (Vector t) (Matrix t) (Vector t) where
136 (<.>) v m = (conj v) `vXm` m
137
138instance Product t => Contraction (Matrix t) (Matrix t) (Matrix t) where
139 (<.>) = mXm
140
141
142--------------------------------------------------------------------------------
143
144class Mul a b c | a b -> c where
145 infixl 7 <>
146 -- | Matrix-matrix, matrix-vector, and vector-matrix products.
147 (<>) :: Product t => a t -> b t -> c t
148
149instance Mul Matrix Matrix Matrix where
150 (<>) = mXm
151
152instance Mul Matrix Vector Vector where
153 (<>) m v = flatten $ m <> asColumn v
154
155instance Mul Vector Matrix Vector where
156 (<>) v m = flatten $ asRow v <> m
157
158--------------------------------------------------------------------------------
159
160class LSDiv c where
161 infixl 7 <\>
162 -- | least squares solution of a linear system, similar to the \\ operator of Matlab\/Octave (based on linearSolveSVD)
163 (<\>) :: Field t => Matrix t -> c t -> c t
164
165instance LSDiv Vector where
166 m <\> v = flatten (linearSolveSVD m (reshape 1 v))
167
168instance LSDiv Matrix where
169 (<\>) = linearSolveSVD
170
171--------------------------------------------------------------------------------
172
173class Konst e d c | d -> c, c -> d
174 where
175 -- |
176 -- >>> konst 7 3 :: Vector Float
177 -- fromList [7.0,7.0,7.0]
178 --
179 -- >>> konst i (3::Int,4::Int)
180 -- (3><4)
181 -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
182 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
183 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
184 --
185 konst :: e -> d -> c e
186
187instance Container Vector e => Konst e Int Vector
188 where
189 konst = konst'
190
191instance Container Vector e => Konst e (Int,Int) Matrix
192 where
193 konst = konst'
194
195--------------------------------------------------------------------------------
196
197class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
198 where
199 -- |
200 -- >>> build 5 (**2) :: Vector Double
201 -- fromList [0.0,1.0,4.0,9.0,16.0]
202 --
203 -- Hilbert matrix of order N:
204 --
205 -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
206 -- >>> putStr . dispf 2 $ hilb 3
207 -- 3x3
208 -- 1.00 0.50 0.33
209 -- 0.50 0.33 0.25
210 -- 0.33 0.25 0.20
211 --
212 build :: d -> f -> c e
213
214instance Container Vector e => Build Int (e -> e) Vector e
215 where
216 build = build'
217
218instance Container Matrix e => Build (Int,Int) (e -> e -> e) Matrix e
219 where
220 build = build'
221
222--------------------------------------------------------------------------------
223
224{- | alternative operator for '(\<.\>)'
225
226x25c7, white diamond
227
228-}
229(◇) :: Contraction a b c => a -> b -> c
230infixl 7 ◇
231(◇) = (<.>)
232
233-- | dot product: @cdot u v = 'udot' ('conj' u) v@
234dot :: (Container Vector t, Product t) => Vector t -> Vector t -> t
235dot u v = udot (conj u) v
236
237--------------------------------------------------------------------------------
238
239optimiseMult :: Monoid (Matrix t) => [Matrix t] -> Matrix t
240optimiseMult = mconcat
241
diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
index bca074f..43e46c5 100644
--- a/packages/hmatrix/src/Numeric/GSL/Vector.hs
+++ b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
@@ -1,16 +1,14 @@
1----------------------------------------------------------------------------- 1-----------------------------------------------------------------------------
2-- | 2-- |
3-- Module : Numeric.GSL.Vector 3-- Module : Numeric.GSL.LinearAlgebra
4-- Copyright : (c) Alberto Ruiz 2007 4-- Copyright : (c) Alberto Ruiz 2007-14
5-- License : GPL 5-- License : GPL
6-- Maintainer : Alberto Ruiz 6-- Maintainer : Alberto Ruiz
7-- Stability : provisional 7-- Stability : provisional
8-- 8--
9-- Low level interface to vector operations.
10--
11----------------------------------------------------------------------------- 9-----------------------------------------------------------------------------
12 10
13module Numeric.GSL.Vector ( 11module Numeric.GSL.LinearAlgebra (
14 RandDist(..), randomVector, 12 RandDist(..), randomVector,
15 saveMatrix, 13 saveMatrix,
16 fwriteVector, freadVector, fprintfVector, fscanfVector 14 fwriteVector, freadVector, fprintfVector, fscanfVector
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/LinearAlgebra.hs
new file mode 100644
index 0000000..1c2b1bb
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra.hs
@@ -0,0 +1,33 @@
1-----------------------------------------------------------------------------
2{- |
3Module : Numeric.LinearAlgebra
4Copyright : (c) Alberto Ruiz 2006-14
5License : GPL
6
7Maintainer : Alberto Ruiz
8Stability : provisional
9
10
11This module reexports all normally required functions for Linear Algebra applications.
12
13It also provides instances of standard classes 'Show', 'Read', 'Eq',
14'Num', 'Fractional', and 'Floating' for 'Vector' and 'Matrix'.
15In arithmetic operations one-component vectors and matrices automatically
16expand to match the dimensions of the other operand.
17
18-}
19-----------------------------------------------------------------------------
20
21module Numeric.LinearAlgebra (
22 module Numeric.Container,
23 module Numeric.LinearAlgebra.Algorithms,
24 module Numeric.LinearAlgebra.IO,
25 module Numeric.LinearAlgebra.Random
26) where
27
28import Numeric.Container
29import Numeric.LinearAlgebra.Algorithms
30import Numeric.LinearAlgebra.Util()
31import Numeric.LinearAlgebra.IO
32import Numeric.LinearAlgebra.Random
33
diff --git a/packages/hmatrix/src/Numeric/LinearAlgebra/GSL.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/GSL.hs
new file mode 100644
index 0000000..1fde621
--- /dev/null
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/GSL.hs
@@ -0,0 +1,136 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Numeric.LinearAlgebra.GSL
4-- Copyright : (c) Alberto Ruiz 2007-14
5-- License : GPL
6-- Maintainer : Alberto Ruiz
7-- Stability : provisional
8--
9-----------------------------------------------------------------------------
10
11module Numeric.LinearAlgebra.GSL (
12 RandDist(..), randomVector,
13 saveMatrix,
14 fwriteVector, freadVector, fprintfVector, fscanfVector,
15 fileDimensions, loadMatrix, fromFile
16) where
17
18import Data.Packed
19import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
20
21import Data.Complex
22import Foreign.Marshal.Alloc(free)
23import Foreign.Ptr(Ptr)
24import Foreign.C.Types
25import Foreign.C.String(newCString)
26import System.IO.Unsafe(unsafePerformIO)
27import System.Process(readProcess)
28
29fromei x = fromIntegral (fromEnum x) :: CInt
30
31-----------------------------------------------------------------------
32
33data RandDist = Uniform -- ^ uniform distribution in [0,1)
34 | Gaussian -- ^ normal distribution with mean zero and standard deviation one
35 deriving Enum
36
37-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
38randomVector :: Int -- ^ seed
39 -> RandDist -- ^ distribution
40 -> Int -- ^ vector size
41 -> Vector Double
42randomVector seed dist n = unsafePerformIO $ do
43 r <- createVector n
44 app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector"
45 return r
46
47foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV
48
49--------------------------------------------------------------------------------
50
51-- | Saves a matrix as 2D ASCII table.
52saveMatrix :: FilePath
53 -> String -- ^ format (%f, %g, %e)
54 -> Matrix Double
55 -> IO ()
56saveMatrix filename fmt m = do
57 charname <- newCString filename
58 charfmt <- newCString fmt
59 let o = if orderOf m == RowMajor then 1 else 0
60 app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
61 free charname
62 free charfmt
63
64foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
65
66--------------------------------------------------------------------------------
67
68-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
69fscanfVector :: FilePath -> Int -> IO (Vector Double)
70fscanfVector filename n = do
71 charname <- newCString filename
72 res <- createVector n
73 app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
74 free charname
75 return res
76
77foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
78
79-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
80fprintfVector :: FilePath -> String -> Vector Double -> IO ()
81fprintfVector filename fmt v = do
82 charname <- newCString filename
83 charfmt <- newCString fmt
84 app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
85 free charname
86 free charfmt
87
88foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
89
90-- | Loads a vector from a binary file (the number of elements must be known in advance).
91freadVector :: FilePath -> Int -> IO (Vector Double)
92freadVector filename n = do
93 charname <- newCString filename
94 res <- createVector n
95 app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
96 free charname
97 return res
98
99foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
100
101-- | Saves the elements of a vector to a binary file.
102fwriteVector :: FilePath -> Vector Double -> IO ()
103fwriteVector filename v = do
104 charname <- newCString filename
105 app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
106 free charname
107
108foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
109
110type PD = Ptr Double --
111type TV = CInt -> PD -> IO CInt --
112type TM = CInt -> CInt -> PD -> IO CInt --
113
114--------------------------------------------------------------------------------
115
116{- | obtains the number of rows and columns in an ASCII data file
117 (provisionally using unix's wc).
118-}
119fileDimensions :: FilePath -> IO (Int,Int)
120fileDimensions fname = do
121 wcres <- readProcess "wc" ["-w",fname] ""
122 contents <- readFile fname
123 let tot = read . head . words $ wcres
124 c = length . head . dropWhile null . map words . lines $ contents
125 if tot > 0
126 then return (tot `div` c, c)
127 else return (0,0)
128
129-- | Loads a matrix from an ASCII file formatted as a 2D table.
130loadMatrix :: FilePath -> IO (Matrix Double)
131loadMatrix file = fromFile file =<< fileDimensions file
132
133-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
134fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
135fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
136
diff --git a/packages/hmatrix/src/Numeric/IO.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/IO.hs
index 08d812b..d2278ee 100644
--- a/packages/hmatrix/src/Numeric/IO.hs
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/IO.hs
@@ -1,28 +1,22 @@
1----------------------------------------------------------------------------- 1-----------------------------------------------------------------------------
2-- | 2-- |
3-- Module : Numeric.IO 3-- Module : Numeric.LinearAlgebra.IO
4-- Copyright : (c) Alberto Ruiz 2010 4-- Copyright : (c) Alberto Ruiz 2007-14
5-- License : GPL 5-- License : GPL
6-- 6-- Maintainer : Alberto Ruiz
7-- Maintainer : Alberto Ruiz <aruiz@um.es>
8-- Stability : provisional 7-- Stability : provisional
9-- Portability : portable
10--
11-- Display, formatting and IO functions for numeric 'Vector' and 'Matrix'
12-- 8--
13----------------------------------------------------------------------------- 9-----------------------------------------------------------------------------
14 10
15module Numeric.IO ( 11module Numeric.LinearAlgebra.IO (
16 dispf, disps, dispcf, vecdisp, latexFormat, format, 12 saveMatrix,
17 loadMatrix, saveMatrix, fromFile, fileDimensions, 13 fwriteVector, freadVector, fprintfVector, fscanfVector,
18 readMatrix, fromArray2D, 14 fileDimensions, loadMatrix, fromFile
19 fscanfVector, fprintfVector, freadVector, fwriteVector
20) where 15) where
21 16
22import Data.Packed 17import Data.Packed
23import Data.Packed.IO 18import Numeric.GSL.LinearAlgebra
24import System.Process(readProcess) 19import System.Process(readProcess)
25import Numeric.GSL.Vector
26 20
27 21
28{- | obtains the number of rows and columns in an ASCII data file 22{- | obtains the number of rows and columns in an ASCII data file
diff --git a/packages/hmatrix/src/Numeric/Random.hs b/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
index 0e722fd..e9ab2b7 100644
--- a/packages/hmatrix/src/Numeric/Random.hs
+++ b/packages/hmatrix/src/Numeric/LinearAlgebra/Random.hs
@@ -1,6 +1,6 @@
1----------------------------------------------------------------------------- 1-----------------------------------------------------------------------------
2-- | 2-- |
3-- Module : Numeric.Random 3-- Module : Numeric.LinearAlgebra.Random
4-- Copyright : (c) Alberto Ruiz 2009-14 4-- Copyright : (c) Alberto Ruiz 2009-14
5-- License : GPL 5-- License : GPL
6-- 6--
@@ -11,7 +11,7 @@
11-- 11--
12----------------------------------------------------------------------------- 12-----------------------------------------------------------------------------
13 13
14module Numeric.Random ( 14module Numeric.LinearAlgebra.Random (
15 Seed, 15 Seed,
16 RandDist(..), 16 RandDist(..),
17 randomVector, 17 randomVector,
@@ -20,7 +20,7 @@ module Numeric.Random (
20 rand, randn 20 rand, randn
21) where 21) where
22 22
23import Numeric.GSL.Vector 23import Numeric.GSL.LinearAlgebra
24import Numeric.Container 24import Numeric.Container
25import Numeric.LinearAlgebra.Algorithms 25import Numeric.LinearAlgebra.Algorithms
26import System.Random(randomIO) 26import System.Random(randomIO)