summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/Data/Packed/Internal/Matrix.hs23
-rw-r--r--lib/Data/Packed/Internal/Vector.hs26
-rw-r--r--lib/Data/Packed/Matrix.hs81
-rw-r--r--lib/Numeric/Container.hs65
-rw-r--r--lib/Numeric/ContainerBoot.hs51
-rw-r--r--lib/Numeric/GSL/Differentiation.hs10
-rw-r--r--lib/Numeric/GSL/Fitting.hs10
-rw-r--r--lib/Numeric/GSL/Fourier.hs4
-rw-r--r--lib/Numeric/GSL/Integration.hs74
-rw-r--r--lib/Numeric/GSL/Minimization.hs5
-rw-r--r--lib/Numeric/GSL/ODE.hs8
-rw-r--r--lib/Numeric/GSL/Polynomials.hs24
-rw-r--r--lib/Numeric/GSL/Root.hs34
-rw-r--r--lib/Numeric/IO.hs35
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs23
-rw-r--r--lib/Numeric/LinearAlgebra/Data.hs9
-rw-r--r--lib/Numeric/LinearAlgebra/Data/Devel.hs4
-rw-r--r--lib/Numeric/LinearAlgebra/Util.hs49
18 files changed, 339 insertions, 196 deletions
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 5e6e649..8709a00 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -132,17 +132,22 @@ mat a f =
132 let m g = do 132 let m g = do
133 g (fi (rows a)) (fi (cols a)) p 133 g (fi (rows a)) (fi (cols a)) p
134 f m 134 f m
135-- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose. 135
136-- 136{- | Creates a vector by concatenation of rows. If the matrix is ColumnMajor, this operation requires a transpose.
137-- @\> flatten ('ident' 3) 137
138-- 9 |> [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]@ 138>>> flatten (ident 3)
139fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0]
140
141-}
139flatten :: Element t => Matrix t -> Vector t 142flatten :: Element t => Matrix t -> Vector t
140flatten = xdat . cmat 143flatten = xdat . cmat
141 144
145{-
142type Mt t s = Int -> Int -> Ptr t -> s 146type Mt t s = Int -> Int -> Ptr t -> s
143-- not yet admitted by my haddock version 147
144-- infixr 6 ::> 148infixr 6 ::>
145-- type t ::> s = Mt t s 149type t ::> s = Mt t s
150-}
146 151
147-- | the inverse of 'Data.Packed.Matrix.fromLists' 152-- | the inverse of 'Data.Packed.Matrix.fromLists'
148toLists :: (Element t) => Matrix t -> [[t]] 153toLists :: (Element t) => Matrix t -> [[t]]
@@ -207,11 +212,11 @@ createMatrix ord r c = do
207{- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns. (GNU-Octave groups by columns. To do it you can define @reshapeF r = trans . reshape r@ 212{- | Creates a matrix from a vector by grouping the elements in rows with the desired number of columns. (GNU-Octave groups by columns. To do it you can define @reshapeF r = trans . reshape r@
208where r is the desired number of rows.) 213where r is the desired number of rows.)
209 214
210@\> reshape 4 ('fromList' [1..12]) 215>>> reshape 4 (fromList [1..12])
211(3><4) 216(3><4)
212 [ 1.0, 2.0, 3.0, 4.0 217 [ 1.0, 2.0, 3.0, 4.0
213 , 5.0, 6.0, 7.0, 8.0 218 , 5.0, 6.0, 7.0, 8.0
214 , 9.0, 10.0, 11.0, 12.0 ]@ 219 , 9.0, 10.0, 11.0, 12.0 ]
215 220
216-} 221-}
217reshape :: Storable t => Int -> Vector t -> Matrix t 222reshape :: Storable t => Int -> Vector t -> Matrix t
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index f0bd9aa..415c972 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -113,16 +113,20 @@ inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r
113 113
114{- | extracts the Vector elements to a list 114{- | extracts the Vector elements to a list
115 115
116@> toList (linspace 5 (1,10)) 116>>> toList (linspace 5 (1,10))
117[1.0,3.25,5.5,7.75,10.0]@ 117[1.0,3.25,5.5,7.75,10.0]
118 118
119-} 119-}
120toList :: Storable a => Vector a -> [a] 120toList :: Storable a => Vector a -> [a]
121toList v = safeRead v $ peekArray (dim v) 121toList v = safeRead v $ peekArray (dim v)
122 122
123{- | An alternative to 'fromList' with explicit dimension. The input 123{- | Create a vector from a list of elements and explicit dimension. The input
124 list is explicitly truncated if it is too long, so it may safely 124 list is explicitly truncated if it is too long, so it may safely
125 be used, for instance, with infinite lists. 125 be used, for instance, with infinite lists.
126
127>>> 5 |> [1..]
128fromList [1.0,2.0,3.0,4.0,5.0]
129
126-} 130-}
127(|>) :: (Storable a) => Int -> [a] -> Vector a 131(|>) :: (Storable a) => Int -> [a] -> Vector a
128infixl 9 |> 132infixl 9 |>
@@ -159,8 +163,8 @@ at v n
159 163
160{- | takes a number of consecutive elements from a Vector 164{- | takes a number of consecutive elements from a Vector
161 165
162@> subVector 2 3 (fromList [1..10]) 166>>> subVector 2 3 (fromList [1..10])
1633 |> [3.0,4.0,5.0]@ 167fromList [3.0,4.0,5.0]
164 168
165-} 169-}
166subVector :: Storable t => Int -- ^ index of the starting element 170subVector :: Storable t => Int -- ^ index of the starting element
@@ -172,8 +176,8 @@ subVector = Vector.slice
172 176
173{- | Reads a vector position: 177{- | Reads a vector position:
174 178
175@> fromList [0..9] \@\> 7 179>>> fromList [0..9] @> 7
1767.0@ 1807.0
177 181
178-} 182-}
179(@>) :: Storable t => Vector t -> Int -> t 183(@>) :: Storable t => Vector t -> Int -> t
@@ -183,8 +187,8 @@ infixl 9 @>
183 187
184{- | concatenate a list of vectors 188{- | concatenate a list of vectors
185 189
186@> vjoin [fromList [1..5], constant 1 3] 190>>> vjoin [fromList [1..5::Double], konst 1 3]
1878 |> [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]@ 191fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]
188 192
189-} 193-}
190vjoin :: Storable t => [Vector t] -> Vector t 194vjoin :: Storable t => [Vector t] -> Vector t
@@ -205,8 +209,8 @@ vjoin as = unsafePerformIO $ do
205 209
206{- | Extract consecutive subvectors of the given sizes. 210{- | Extract consecutive subvectors of the given sizes.
207 211
208@> takesV [3,4] (linspace 10 (1,10)) 212>>> takesV [3,4] (linspace 10 (1,10::Double))
209[3 |> [1.0,2.0,3.0],4 |> [4.0,5.0,6.0,7.0]]@ 213[fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]]
210 214
211-} 215-}
212takesV :: Storable t => [Int] -> Vector t -> [Vector t] 216takesV :: Storable t => [Int] -> Vector t -> [Vector t]
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs
index 5365c1c..f72bd15 100644
--- a/lib/Data/Packed/Matrix.hs
+++ b/lib/Data/Packed/Matrix.hs
@@ -116,13 +116,10 @@ Single row-column components are automatically expanded to match the
116corresponding common row and column: 116corresponding common row and column:
117 117
118@ 118@
119> let disp = putStr . dispf 2 119disp = putStr . dispf 2
120> let vector xs = fromList xs :: Vector Double 120@
121> let diagl = diag . vector
122> let rowm = asRow . vector
123
124> disp $ fromBlocks [[ident 5, 7, rowm[10,20]], [3, diagl[1,2,3], 0]]
125 121
122>>> disp $ fromBlocks [[ident 5, 7, row[10,20]], [3, diagl[1,2,3], 0]]
1268x10 1238x10
1271 0 0 0 0 7 7 7 10 20 1241 0 0 0 0 7 7 7 10 20
1280 1 0 0 0 7 7 7 10 20 1250 1 0 0 0 7 7 7 10 20
@@ -132,7 +129,6 @@ corresponding common row and column:
1323 3 3 3 3 1 0 0 0 0 1293 3 3 3 3 1 0 0 0 0
1333 3 3 3 3 0 2 0 0 0 1303 3 3 3 3 0 2 0 0 0
1343 3 3 3 3 0 0 3 0 0 1313 3 3 3 3 0 0 3 0 0
135@
136 132
137-} 133-}
138fromBlocks :: Element t => [[Matrix t]] -> Matrix t 134fromBlocks :: Element t => [[Matrix t]] -> Matrix t
@@ -163,7 +159,19 @@ adaptBlocks ms = ms' where
163 159
164-------------------------------------------------------------------------------- 160--------------------------------------------------------------------------------
165 161
166-- | create a block diagonal matrix 162{- | create a block diagonal matrix
163
164>>> disp 2 $ diagBlock [konst 1 (2,2), konst 2 (3,5), col [5,7]]
1657x8
1661 1 0 0 0 0 0 0
1671 1 0 0 0 0 0 0
1680 0 2 2 2 2 2 0
1690 0 2 2 2 2 2 0
1700 0 2 2 2 2 2 0
1710 0 0 0 0 0 0 5
1720 0 0 0 0 0 0 7
173
174-}
167diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t 175diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t
168diagBlock ms = fromBlocks $ zipWith f ms [0..] 176diagBlock ms = fromBlocks $ zipWith f ms [0..]
169 where 177 where
@@ -186,12 +194,13 @@ fliprl m = fromColumns . reverse . toColumns $ m
186 194
187{- | creates a rectangular diagonal matrix: 195{- | creates a rectangular diagonal matrix:
188 196
189@> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double 197>>> diagRect 7 (fromList [10,20,30]) 4 5 :: Matrix Double
190(4><5) 198(4><5)
191 [ 10.0, 7.0, 7.0, 7.0, 7.0 199 [ 10.0, 7.0, 7.0, 7.0, 7.0
192 , 7.0, 20.0, 7.0, 7.0, 7.0 200 , 7.0, 20.0, 7.0, 7.0, 7.0
193 , 7.0, 7.0, 30.0, 7.0, 7.0 201 , 7.0, 7.0, 30.0, 7.0, 7.0
194 , 7.0, 7.0, 7.0, 7.0, 7.0 ]@ 202 , 7.0, 7.0, 7.0, 7.0, 7.0 ]
203
195-} 204-}
196diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t 205diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t
197diagRect z v r c = ST.runSTMatrix $ do 206diagRect z v r c = ST.runSTMatrix $ do
@@ -208,10 +217,10 @@ takeDiag m = fromList [flatten m `at` (k*cols m+k) | k <- [0 .. min (rows m) (co
208 217
209{- | An easy way to create a matrix: 218{- | An easy way to create a matrix:
210 219
211@\> (2><3)[1..6] 220>>> (2><3)[2,4,7,-3,11,0]
212(2><3) 221(2><3)
213 [ 1.0, 2.0, 3.0 222 [ 2.0, 4.0, 7.0
214 , 4.0, 5.0, 6.0 ]@ 223 , -3.0, 11.0, 0.0 ]
215 224
216This is the format produced by the instances of Show (Matrix a), which 225This is the format produced by the instances of Show (Matrix a), which
217can also be used for input. 226can also be used for input.
@@ -219,12 +228,11 @@ can also be used for input.
219The input list is explicitly truncated, so that it can 228The input list is explicitly truncated, so that it can
220safely be used with lists that are too long (like infinite lists). 229safely be used with lists that are too long (like infinite lists).
221 230
222Example: 231>>> (2><3)[1..]
223
224@\> (2><3)[1..]
225(2><3) 232(2><3)
226 [ 1.0, 2.0, 3.0 233 [ 1.0, 2.0, 3.0
227 , 4.0, 5.0, 6.0 ]@ 234 , 4.0, 5.0, 6.0 ]
235
228 236
229-} 237-}
230(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a 238(><) :: (Storable a) => Int -> Int -> [a] -> Matrix a
@@ -253,20 +261,35 @@ dropColumns n mt = subMatrix (0,n) (rows mt, cols mt - n) mt
253 261
254{- | Creates a 'Matrix' from a list of lists (considered as rows). 262{- | Creates a 'Matrix' from a list of lists (considered as rows).
255 263
256@\> fromLists [[1,2],[3,4],[5,6]] 264>>> fromLists [[1,2],[3,4],[5,6]]
257(3><2) 265(3><2)
258 [ 1.0, 2.0 266 [ 1.0, 2.0
259 , 3.0, 4.0 267 , 3.0, 4.0
260 , 5.0, 6.0 ]@ 268 , 5.0, 6.0 ]
269
261-} 270-}
262fromLists :: Element t => [[t]] -> Matrix t 271fromLists :: Element t => [[t]] -> Matrix t
263fromLists = fromRows . map fromList 272fromLists = fromRows . map fromList
264 273
265-- | creates a 1-row matrix from a vector 274-- | creates a 1-row matrix from a vector
275--
276-- >>> asRow (fromList [1..5])
277-- (1><5)
278-- [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
279--
266asRow :: Storable a => Vector a -> Matrix a 280asRow :: Storable a => Vector a -> Matrix a
267asRow v = reshape (dim v) v 281asRow v = reshape (dim v) v
268 282
269-- | creates a 1-column matrix from a vector 283-- | creates a 1-column matrix from a vector
284--
285-- >>> asColumn (fromList [1..5])
286-- (5><1)
287-- [ 1.0
288-- , 2.0
289-- , 3.0
290-- , 4.0
291-- , 5.0 ]
292--
270asColumn :: Storable a => Vector a -> Matrix a 293asColumn :: Storable a => Vector a -> Matrix a
271asColumn v = reshape 1 v 294asColumn v = reshape 1 v
272 295
@@ -307,12 +330,12 @@ extractRows l m = fromRows $ extract (toRows m) l
307 330
308{- | creates matrix by repetition of a matrix a given number of rows and columns 331{- | creates matrix by repetition of a matrix a given number of rows and columns
309 332
310@> repmat (ident 2) 2 3 :: Matrix Double 333>>> repmat (ident 2) 2 3
311(4><6) 334(4><6)
312 [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0 335 [ 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
313 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 336 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0
314 , 1.0, 0.0, 1.0, 0.0, 1.0, 0.0 337 , 1.0, 0.0, 1.0, 0.0, 1.0, 0.0
315 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]@ 338 , 0.0, 1.0, 0.0, 1.0, 0.0, 1.0 ]
316 339
317-} 340-}
318repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t 341repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t
@@ -376,13 +399,14 @@ mk c g = \k -> g (divMod k c)
376 399
377{- | 400{- |
378 401
379@ghci> mapMatrixWithIndexM_ (\\(i,j) v -> printf \"m[%.0f,%.0f] = %.f\\n\" i j v :: IO()) ((2><3)[1 :: Double ..]) 402>>> mapMatrixWithIndexM_ (\(i,j) v -> printf "m[%d,%d] = %.f\n" i j v :: IO()) ((2><3)[1 :: Double ..])
380m[0,0] = 1 403m[0,0] = 1
381m[0,1] = 2 404m[0,1] = 2
382m[0,2] = 3 405m[0,2] = 3
383m[1,0] = 4 406m[1,0] = 4
384m[1,1] = 5 407m[1,1] = 5
385m[1,2] = 6@ 408m[1,2] = 6
409
386-} 410-}
387mapMatrixWithIndexM_ 411mapMatrixWithIndexM_
388 :: (Element a, Num a, Monad m) => 412 :: (Element a, Num a, Monad m) =>
@@ -393,11 +417,12 @@ mapMatrixWithIndexM_ g m = mapVectorWithIndexM_ (mk c g) . flatten $ m
393 417
394{- | 418{- |
395 419
396@ghci> mapMatrixWithIndexM (\\(i,j) v -> Just $ 100*v + 10*i + j) (ident 3:: Matrix Double) 420>>> mapMatrixWithIndexM (\(i,j) v -> Just $ 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
397Just (3><3) 421Just (3><3)
398 [ 100.0, 1.0, 2.0 422 [ 100.0, 1.0, 2.0
399 , 10.0, 111.0, 12.0 423 , 10.0, 111.0, 12.0
400 , 20.0, 21.0, 122.0 ]@ 424 , 20.0, 21.0, 122.0 ]
425
401-} 426-}
402mapMatrixWithIndexM 427mapMatrixWithIndexM
403 :: (Element a, Storable b, Monad m) => 428 :: (Element a, Storable b, Monad m) =>
@@ -407,11 +432,13 @@ mapMatrixWithIndexM g m = liftM (reshape c) . mapVectorWithIndexM (mk c g) . fla
407 c = cols m 432 c = cols m
408 433
409{- | 434{- |
410@ghci> mapMatrixWithIndex (\\(i,j) v -> 100*v + 10*i + j) (ident 3:: Matrix Double) 435
436>>> mapMatrixWithIndex (\\(i,j) v -> 100*v + 10*fromIntegral i + fromIntegral j) (ident 3:: Matrix Double)
411(3><3) 437(3><3)
412 [ 100.0, 1.0, 2.0 438 [ 100.0, 1.0, 2.0
413 , 10.0, 111.0, 12.0 439 , 10.0, 111.0, 12.0
414 , 20.0, 21.0, 122.0 ]@ 440 , 20.0, 21.0, 122.0 ]
441
415 -} 442 -}
416mapMatrixWithIndex 443mapMatrixWithIndex
417 :: (Element a, Storable b) => 444 :: (Element a, Storable b) =>
diff --git a/lib/Numeric/Container.hs b/lib/Numeric/Container.hs
index d2fd351..5339c7e 100644
--- a/lib/Numeric/Container.hs
+++ b/lib/Numeric/Container.hs
@@ -85,8 +85,8 @@ constant = constantD-- about 2x faster
85 85
86{- | Creates a real vector containing a range of values: 86{- | Creates a real vector containing a range of values:
87 87
88@\> linspace 5 (-3,7) 88>>> linspace 5 (-3,7)
895 |> [-3.0,-0.5,2.0,4.5,7.0]@ 89fromList [-3.0,-0.5,2.0,4.5,7.0]@
90 90
91Logarithmic spacing can be defined as follows: 91Logarithmic spacing can be defined as follows:
92 92
@@ -105,7 +105,32 @@ cdot u v = udot (conj u) v
105class Contraction a b c | a b -> c, a c -> b, b c -> a 105class Contraction a b c | a b -> c, a c -> b, b c -> a
106 where 106 where
107 infixl 7 <> 107 infixl 7 <>
108 -- | matrix-matrix product, matrix-vector product, unconjugated dot product 108 {- | matrix-matrix product, matrix-vector product, unconjugated dot product
109
110>>> let a = (3><4) [1..] :: Matrix Double
111
112>>> a
113(3><4)
114 [ 1.0, 2.0, 3.0, 4.0
115 , 5.0, 6.0, 7.0, 8.0
116 , 9.0, 10.0, 11.0, 12.0 ]
117
118>>> disp 2 (a <> trans a)
1193x3
120 30 70 110
121 70 174 278
122110 278 446
123
124>>> a <> fromList [1,0,2,-1::Double]
125fromList [3.0,11.0,19.0]
126
127>>> fromList [1,2,3::Double] <> a
128fromList [38.0,44.0,50.0,56.0]
129
130>>> fromList [1,i] <> fromList[2*i+1,3]
1311.0 :+ 5.0
132
133-}
109 (<>) :: a -> b -> c 134 (<>) :: a -> b -> c
110 135
111instance Product t => Contraction (Vector t) (Vector t) t where 136instance Product t => Contraction (Vector t) (Vector t) t where
@@ -135,9 +160,14 @@ instance LSDiv Matrix Matrix where
135 160
136-------------------------------------------------------- 161--------------------------------------------------------
137 162
138-- | dot product : @u · v = 'cdot' u v@ 163{- | dot product : @u · v = 'cdot' u v@
139-- 164
140-- unicode 0x00b7, Alt-Gr . 165 unicode 0x00b7, Alt-Gr .
166
167>>> fromList [1,i] · fromList[2*i+1,3]
1681.0 :+ (-1.0)
169
170-}
141(·) :: (Container Vector t, Product t) => Vector t -> Vector t -> t 171(·) :: (Container Vector t, Product t) => Vector t -> Vector t -> t
142infixl 7 · 172infixl 7 ·
143u · v = cdot u v 173u · v = cdot u v
@@ -147,6 +177,16 @@ u · v = cdot u v
147-- bidirectional type inference 177-- bidirectional type inference
148class Konst e d c | d -> c, c -> d 178class Konst e d c | d -> c, c -> d
149 where 179 where
180 -- |
181 -- >>> konst 7 3 :: Vector Float
182 -- fromList [7.0,7.0,7.0]
183 --
184 -- >>> konst i (3::Int,4::Int)
185 -- (3><4)
186 -- [ 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
187 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0
188 -- , 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0, 0.0 :+ 1.0 ]
189 --
150 konst :: e -> d -> c e 190 konst :: e -> d -> c e
151 191
152instance Container Vector e => Konst e Int Vector 192instance Container Vector e => Konst e Int Vector
@@ -161,6 +201,19 @@ instance Container Vector e => Konst e (Int,Int) Matrix
161 201
162class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f 202class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f
163 where 203 where
204 -- |
205 -- >>> build 5 (**2) :: Vector Double
206 -- fromList [0.0,1.0,4.0,9.0,16.0]
207 --
208 -- Hilbert matrix of order N:
209 --
210 -- >>> let hilb n = build (n,n) (\i j -> 1/(i+j+1)) :: Matrix Double
211 -- >>> putStr . dispf 2 $ hilb 3
212 -- 3x3
213 -- 1.00 0.50 0.33
214 -- 0.50 0.33 0.25
215 -- 0.33 0.25 0.20
216 --
164 build :: d -> f -> c e 217 build :: d -> f -> c e
165 218
166instance Container Vector e => Build Int (e -> e) Vector e 219instance Container Vector e => Build Int (e -> e) Vector e
diff --git a/lib/Numeric/ContainerBoot.hs b/lib/Numeric/ContainerBoot.hs
index d61633e..a333489 100644
--- a/lib/Numeric/ContainerBoot.hs
+++ b/lib/Numeric/ContainerBoot.hs
@@ -66,6 +66,11 @@ type instance ArgOf Matrix a = a -> a -> a
66-- | Basic element-by-element functions for numeric containers 66-- | Basic element-by-element functions for numeric containers
67class (Complexable c, Fractional e, Element e) => Container c e where 67class (Complexable c, Fractional e, Element e) => Container c e where
68 -- | create a structure with a single element 68 -- | create a structure with a single element
69 --
70 -- >>> let v = fromList [1..3::Double]
71 -- >>> v / scalar (norm2 v)
72 -- fromList [0.2672612419124244,0.5345224838248488,0.8017837257372732]
73 --
69 scalar :: e -> c e 74 scalar :: e -> c e
70 -- | complex conjugate 75 -- | complex conjugate
71 conj :: c e -> c e 76 conj :: c e -> c e
@@ -114,8 +119,9 @@ class (Complexable c, Fractional e, Element e) => Container c e where
114 119
115 -- | A more efficient implementation of @cmap (\\x -> if x>0 then 1 else 0)@ 120 -- | A more efficient implementation of @cmap (\\x -> if x>0 then 1 else 0)@
116 -- 121 --
117 -- @> step $ linspace 5 (-1,1::Double) 122 -- >>> step $ linspace 5 (-1,1::Double)
118 -- 5 |> [0.0,0.0,0.0,1.0,1.0]@ 123 -- 5 |> [0.0,0.0,0.0,1.0,1.0]
124 --
119 125
120 step :: RealElement e => c e -> c e 126 step :: RealElement e => c e -> c e
121 127
@@ -123,11 +129,12 @@ class (Complexable c, Fractional e, Element e) => Container c e where
123 -- 129 --
124 -- Arguments with any dimension = 1 are automatically expanded: 130 -- Arguments with any dimension = 1 are automatically expanded:
125 -- 131 --
126 -- @> cond ((1>\<4)[1..]) ((3>\<1)[1..]) 0 100 ((3>\<4)[1..]) :: Matrix Double 132 -- >>> cond ((1><4)[1..]) ((3><1)[1..]) 0 100 ((3><4)[1..]) :: Matrix Double
127 -- (3><4) 133 -- (3><4)
128 -- [ 100.0, 2.0, 3.0, 4.0 134 -- [ 100.0, 2.0, 3.0, 4.0
129 -- , 0.0, 100.0, 7.0, 8.0 135 -- , 0.0, 100.0, 7.0, 8.0
130 -- , 0.0, 0.0, 100.0, 12.0 ]@ 136 -- , 0.0, 0.0, 100.0, 12.0 ]
137 --
131 138
132 cond :: RealElement e 139 cond :: RealElement e
133 => c e -- ^ a 140 => c e -- ^ a
@@ -139,16 +146,22 @@ class (Complexable c, Fractional e, Element e) => Container c e where
139 146
140 -- | Find index of elements which satisfy a predicate 147 -- | Find index of elements which satisfy a predicate
141 -- 148 --
142 -- @> find (>0) (ident 3 :: Matrix Double) 149 -- >>> find (>0) (ident 3 :: Matrix Double)
143 -- [(0,0),(1,1),(2,2)]@ 150 -- [(0,0),(1,1),(2,2)]
151 --
144 152
145 find :: (e -> Bool) -> c e -> [IndexOf c] 153 find :: (e -> Bool) -> c e -> [IndexOf c]
146 154
147 -- | Create a structure from an association list 155 -- | Create a structure from an association list
148 -- 156 --
149 -- @> assoc 5 0 [(2,7),(1,3)] :: Vector Double 157 -- >>> assoc 5 0 [(3,7),(1,4)] :: Vector Double
150 -- 5 |> [0.0,3.0,7.0,0.0,0.0]@ 158 -- fromList [0.0,4.0,0.0,7.0,0.0]
151 159 --
160 -- >>> assoc (2,3) 0 [((0,2),7),((1,0),2*i-3)] :: Matrix (Complex Double)
161 -- (2><3)
162 -- [ 0.0 :+ 0.0, 0.0 :+ 0.0, 7.0 :+ 0.0
163 -- , (-3.0) :+ 2.0, 0.0 :+ 0.0, 0.0 :+ 0.0 ]
164 --
152 assoc :: IndexOf c -- ^ size 165 assoc :: IndexOf c -- ^ size
153 -> e -- ^ default value 166 -> e -- ^ default value
154 -> [(IndexOf c, e)] -- ^ association list 167 -> [(IndexOf c, e)] -- ^ association list
@@ -156,13 +169,19 @@ class (Complexable c, Fractional e, Element e) => Container c e where
156 169
157 -- | Modify a structure using an update function 170 -- | Modify a structure using an update function
158 -- 171 --
159 -- @> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double 172 -- >>> accum (ident 5) (+) [((1,1),5),((0,3),3)] :: Matrix Double
160 -- (5><5) 173 -- (5><5)
161 -- [ 1.0, 0.0, 0.0, 3.0, 0.0 174 -- [ 1.0, 0.0, 0.0, 3.0, 0.0
162 -- , 0.0, 6.0, 0.0, 0.0, 0.0 175 -- , 0.0, 6.0, 0.0, 0.0, 0.0
163 -- , 0.0, 0.0, 1.0, 0.0, 0.0 176 -- , 0.0, 0.0, 1.0, 0.0, 0.0
164 -- , 0.0, 0.0, 0.0, 1.0, 0.0 177 -- , 0.0, 0.0, 0.0, 1.0, 0.0
165 -- , 0.0, 0.0, 0.0, 0.0, 1.0 ]@ 178 -- , 0.0, 0.0, 0.0, 0.0, 1.0 ]
179 --
180 -- computation of histogram:
181 --
182 -- >>> accum (konst 0 7) (+) (map (flip (,) 1) [4,5,4,1,5,2,5]) :: Vector Double
183 -- fromList [0.0,1.0,1.0,0.0,2.0,3.0,0.0]
184 --
166 185
167 accum :: c e -- ^ initial structure 186 accum :: c e -- ^ initial structure
168 -> (e -> e -> e) -- ^ update function 187 -> (e -> e -> e) -- ^ update function
@@ -382,11 +401,12 @@ vXm v m = flatten $ (asRow v) `mXm` m
382 401
383{- | Outer product of two vectors. 402{- | Outer product of two vectors.
384 403
385@\> 'fromList' [1,2,3] \`outer\` 'fromList' [5,2,3] 404>>> fromList [1,2,3] `outer` fromList [5,2,3]
386(3><3) 405(3><3)
387 [ 5.0, 2.0, 3.0 406 [ 5.0, 2.0, 3.0
388 , 10.0, 4.0, 6.0 407 , 10.0, 4.0, 6.0
389 , 15.0, 6.0, 9.0 ]@ 408 , 15.0, 6.0, 9.0 ]
409
390-} 410-}
391outer :: (Product t) => Vector t -> Vector t -> Matrix t 411outer :: (Product t) => Vector t -> Vector t -> Matrix t
392outer u v = asColumn u `multiply` asRow v 412outer u v = asColumn u `multiply` asRow v
@@ -402,7 +422,7 @@ m2=(4><3)
402 , 7.0, 8.0, 9.0 422 , 7.0, 8.0, 9.0
403 , 10.0, 11.0, 12.0 ]@ 423 , 10.0, 11.0, 12.0 ]@
404 424
405@\> kronecker m1 m2 425>>> kronecker m1 m2
406(8><9) 426(8><9)
407 [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0 427 [ 1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 0.0, 0.0, 0.0
408 , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0 428 , 4.0, 5.0, 6.0, 8.0, 10.0, 12.0, 0.0, 0.0, 0.0
@@ -411,7 +431,8 @@ m2=(4><3)
411 , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0 431 , 0.0, 0.0, 0.0, -1.0, -2.0, -3.0, 3.0, 6.0, 9.0
412 , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0 432 , 0.0, 0.0, 0.0, -4.0, -5.0, -6.0, 12.0, 15.0, 18.0
413 , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 433 , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0
414 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ 434 , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]
435
415-} 436-}
416kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t 437kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t
417kronecker a b = fromBlocks 438kronecker a b = fromBlocks
diff --git a/lib/Numeric/GSL/Differentiation.hs b/lib/Numeric/GSL/Differentiation.hs
index d2a332c..93c5007 100644
--- a/lib/Numeric/GSL/Differentiation.hs
+++ b/lib/Numeric/GSL/Differentiation.hs
@@ -55,11 +55,11 @@ foreign import ccall safe "gsl-aux.h deriv"
55 55
56{- | Adaptive central difference algorithm, /gsl_deriv_central/. For example: 56{- | Adaptive central difference algorithm, /gsl_deriv_central/. For example:
57 57
58> > let deriv = derivCentral 0.01 58>>> let deriv = derivCentral 0.01
59> > deriv sin (pi/4) 59>>> deriv sin (pi/4)
60>(0.7071067812000676,1.0600063101654055e-10) 60(0.7071067812000676,1.0600063101654055e-10)
61> > cos (pi/4) 61>>> cos (pi/4)
62>0.7071067811865476 620.7071067811865476
63 63
64-} 64-}
65derivCentral :: Double -- ^ initial step size 65derivCentral :: Double -- ^ initial step size
diff --git a/lib/Numeric/GSL/Fitting.hs b/lib/Numeric/GSL/Fitting.hs
index 6343b76..c4f3a91 100644
--- a/lib/Numeric/GSL/Fitting.hs
+++ b/lib/Numeric/GSL/Fitting.hs
@@ -13,7 +13,8 @@ Nonlinear Least-Squares Fitting
13 13
14The example program in the GSL manual (see examples/fitting.hs): 14The example program in the GSL manual (see examples/fitting.hs):
15 15
16@dat = [ 16@
17dat = [
17 ([0.0],([6.0133918608118675],0.1)), 18 ([0.0],([6.0133918608118675],0.1)),
18 ([1.0],([5.5153769909966535],0.1)), 19 ([1.0],([5.5153769909966535],0.1)),
19 ([2.0],([5.261094606015287],0.1)), 20 ([2.0],([5.261094606015287],0.1)),
@@ -25,8 +26,9 @@ expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]
25expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]] 26expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
26 27
27(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] 28(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
29@
28 30
29\> path 31>>> path
30(6><5) 32(6><5)
31 [ 1.0, 76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797 33 [ 1.0, 76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797
32 , 2.0, 37.683816318260355, 2.858760367632973, 8.092094813253975e-2, 1.4479636296208662 34 , 2.0, 37.683816318260355, 2.858760367632973, 8.092094813253975e-2, 1.4479636296208662
@@ -34,10 +36,10 @@ expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]
34 , 4.0, 5.630494933603935, 5.021755718065913, 0.10287787128056883, 1.0338835440862608 36 , 4.0, 5.630494933603935, 5.021755718065913, 0.10287787128056883, 1.0338835440862608
35 , 5.0, 5.443976278682909, 5.045204331329302, 0.10405523433131504, 1.019416067207375 37 , 5.0, 5.443976278682909, 5.045204331329302, 0.10405523433131504, 1.019416067207375
36 , 6.0, 5.4439736648994685, 5.045357818922331, 0.10404905846029407, 1.0192487112786812 ] 38 , 6.0, 5.4439736648994685, 5.045357818922331, 0.10404905846029407, 1.0192487112786812 ]
37\> sol 39>>> sol
38[(5.045357818922331,6.027976702418132e-2), 40[(5.045357818922331,6.027976702418132e-2),
39(0.10404905846029407,3.157045047172834e-3), 41(0.10404905846029407,3.157045047172834e-3),
40(1.0192487112786812,3.782067731353722e-2)]@ 42(1.0192487112786812,3.782067731353722e-2)]
41 43
42-} 44-}
43----------------------------------------------------------------------------- 45-----------------------------------------------------------------------------
diff --git a/lib/Numeric/GSL/Fourier.hs b/lib/Numeric/GSL/Fourier.hs
index 4ef19b3..86aedd6 100644
--- a/lib/Numeric/GSL/Fourier.hs
+++ b/lib/Numeric/GSL/Fourier.hs
@@ -35,8 +35,8 @@ foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCVCV
35 35
36{- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave. 36{- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave.
37 37
38@> fft ('fromList' [1,2,3,4]) 38>>> fft (fromList [1,2,3,4])
39vector (4) [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]@ 39fromList [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]
40 40
41-} 41-}
42fft :: Vector (Complex Double) -> Vector (Complex Double) 42fft :: Vector (Complex Double) -> Vector (Complex Double)
diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs
index 7c1cdb9..5f0a415 100644
--- a/lib/Numeric/GSL/Integration.hs
+++ b/lib/Numeric/GSL/Integration.hs
@@ -35,16 +35,16 @@ eps = 1e-12
35 35
36{- | conversion of Haskell functions into function pointers that can be used in the C side 36{- | conversion of Haskell functions into function pointers that can be used in the C side
37-} 37-}
38foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) 38foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double))
39 39
40-------------------------------------------------------------------- 40--------------------------------------------------------------------
41{- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example: 41{- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example:
42 42
43@\> let quad = integrateQAGS 1E-9 1000 43>>> let quad = integrateQAGS 1E-9 1000
44\> let f a x = x**(-0.5) * log (a*x) 44>>> let f a x = x**(-0.5) * log (a*x)
45\> quad (f 1) 0 1 45>>> quad (f 1) 0 1
46(-3.999999999999974,4.871658632055187e-13)@ 46(-3.999999999999974,4.871658632055187e-13)
47 47
48-} 48-}
49 49
50integrateQAGS :: Double -- ^ precision (e.g. 1E-9) 50integrateQAGS :: Double -- ^ precision (e.g. 1E-9)
@@ -56,7 +56,7 @@ integrateQAGS :: Double -- ^ precision (e.g. 1E-9)
56integrateQAGS prec n f a b = unsafePerformIO $ do 56integrateQAGS prec n f a b = unsafePerformIO $ do
57 r <- malloc 57 r <- malloc
58 e <- malloc 58 e <- malloc
59 fp <- mkfun (\x _ -> f x) 59 fp <- mkfun (\x _ -> f x)
60 c_integrate_qags fp a b eps prec (fromIntegral n) r e // check "integrate_qags" 60 c_integrate_qags fp a b eps prec (fromIntegral n) r e // check "integrate_qags"
61 vr <- peek r 61 vr <- peek r
62 ve <- peek e 62 ve <- peek e
@@ -73,10 +73,10 @@ foreign import ccall safe "integrate_qags" c_integrate_qags
73----------------------------------------------------------------- 73-----------------------------------------------------------------
74{- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example: 74{- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example:
75 75
76@\> let quad = integrateQNG 1E-6 76>>> let quad = integrateQNG 1E-6
77\> quad (\\x -> 4\/(1+x*x)) 0 1 77>>> quad (\x -> 4/(1+x*x)) 0 1
78(3.141592653589793,3.487868498008632e-14)@ 78(3.141592653589793,3.487868498008632e-14)
79 79
80-} 80-}
81 81
82integrateQNG :: Double -- ^ precision (e.g. 1E-9) 82integrateQNG :: Double -- ^ precision (e.g. 1E-9)
@@ -87,7 +87,7 @@ integrateQNG :: Double -- ^ precision (e.g. 1E-9)
87integrateQNG prec f a b = unsafePerformIO $ do 87integrateQNG prec f a b = unsafePerformIO $ do
88 r <- malloc 88 r <- malloc
89 e <- malloc 89 e <- malloc
90 fp <- mkfun (\x _ -> f x) 90 fp <- mkfun (\x _ -> f x)
91 c_integrate_qng fp a b eps prec r e // check "integrate_qng" 91 c_integrate_qng fp a b eps prec r e // check "integrate_qng"
92 vr <- peek r 92 vr <- peek r
93 ve <- peek e 93 ve <- peek e
@@ -102,14 +102,14 @@ foreign import ccall safe "integrate_qng" c_integrate_qng
102 -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt 102 -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt
103 103
104-------------------------------------------------------------------- 104--------------------------------------------------------------------
105{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS). 105{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS).
106For example: 106For example:
107 107
108@\> let quad = integrateQAGI 1E-9 1000 108>>> let quad = integrateQAGI 1E-9 1000
109\> let f a x = exp(-a * x^2) 109>>> let f a x = exp(-a * x^2)
110\> quad (f 0.5) 110>>> quad (f 0.5)
111(2.5066282746310002,6.229215880648858e-11)@ 111(2.5066282746310002,6.229215880648858e-11)
112 112
113-} 113-}
114 114
115integrateQAGI :: Double -- ^ precision (e.g. 1E-9) 115integrateQAGI :: Double -- ^ precision (e.g. 1E-9)
@@ -119,7 +119,7 @@ integrateQAGI :: Double -- ^ precision (e.g. 1E-9)
119integrateQAGI prec n f = unsafePerformIO $ do 119integrateQAGI prec n f = unsafePerformIO $ do
120 r <- malloc 120 r <- malloc
121 e <- malloc 121 e <- malloc
122 fp <- mkfun (\x _ -> f x) 122 fp <- mkfun (\x _ -> f x)
123 c_integrate_qagi fp eps prec (fromIntegral n) r e // check "integrate_qagi" 123 c_integrate_qagi fp eps prec (fromIntegral n) r e // check "integrate_qagi"
124 vr <- peek r 124 vr <- peek r
125 ve <- peek e 125 ve <- peek e
@@ -134,14 +134,14 @@ foreign import ccall safe "integrate_qagi" c_integrate_qagi
134 -> CInt -> Ptr Double -> Ptr Double -> IO CInt 134 -> CInt -> Ptr Double -> Ptr Double -> IO CInt
135 135
136-------------------------------------------------------------------- 136--------------------------------------------------------------------
137{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf). 137{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf).
138For example: 138For example:
139 139
140@\> let quad = integrateQAGIU 1E-9 1000 140>>> let quad = integrateQAGIU 1E-9 1000
141\> let f a x = exp(-a * x^2) 141>>> let f a x = exp(-a * x^2)
142\> quad (f 0.5) 0 142>>> quad (f 0.5) 0
143(1.2533141373155001,3.114607940324429e-11)@ 143(1.2533141373155001,3.114607940324429e-11)
144 144
145-} 145-}
146 146
147integrateQAGIU :: Double -- ^ precision (e.g. 1E-9) 147integrateQAGIU :: Double -- ^ precision (e.g. 1E-9)
@@ -152,7 +152,7 @@ integrateQAGIU :: Double -- ^ precision (e.g. 1E-9)
152integrateQAGIU prec n f a = unsafePerformIO $ do 152integrateQAGIU prec n f a = unsafePerformIO $ do
153 r <- malloc 153 r <- malloc
154 e <- malloc 154 e <- malloc
155 fp <- mkfun (\x _ -> f x) 155 fp <- mkfun (\x _ -> f x)
156 c_integrate_qagiu fp a eps prec (fromIntegral n) r e // check "integrate_qagiu" 156 c_integrate_qagiu fp a eps prec (fromIntegral n) r e // check "integrate_qagiu"
157 vr <- peek r 157 vr <- peek r
158 ve <- peek e 158 ve <- peek e
@@ -167,14 +167,14 @@ foreign import ccall safe "integrate_qagiu" c_integrate_qagiu
167 -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt 167 -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
168 168
169-------------------------------------------------------------------- 169--------------------------------------------------------------------
170{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b). 170{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b).
171For example: 171For example:
172 172
173@\> let quad = integrateQAGIL 1E-9 1000 173>>> let quad = integrateQAGIL 1E-9 1000
174\> let f a x = exp(-a * x^2) 174>>> let f a x = exp(-a * x^2)
175\> quad (f 0.5) 0 175>>> quad (f 0.5) 0
176(1.2533141373155001,3.114607940324429e-11)@ 176(1.2533141373155001,3.114607940324429e-11)
177 177
178-} 178-}
179 179
180integrateQAGIL :: Double -- ^ precision (e.g. 1E-9) 180integrateQAGIL :: Double -- ^ precision (e.g. 1E-9)
@@ -185,7 +185,7 @@ integrateQAGIL :: Double -- ^ precision (e.g. 1E-9)
185integrateQAGIL prec n f b = unsafePerformIO $ do 185integrateQAGIL prec n f b = unsafePerformIO $ do
186 r <- malloc 186 r <- malloc
187 e <- malloc 187 e <- malloc
188 fp <- mkfun (\x _ -> f x) 188 fp <- mkfun (\x _ -> f x)
189 c_integrate_qagil fp b eps prec (fromIntegral n) r e // check "integrate_qagil" 189 c_integrate_qagil fp b eps prec (fromIntegral n) r e // check "integrate_qagil"
190 vr <- peek r 190 vr <- peek r
191 ve <- peek e 191 ve <- peek e
@@ -212,10 +212,10 @@ routines in QUADPACK, yet fails less often for difficult integrands.@
212 212
213For example: 213For example:
214 214
215@\> let quad = integrateCQUAD 1E-12 1000 215>>> let quad = integrateCQUAD 1E-12 1000
216\> let f a x = exp(-a * x^2) 216>>> let f a x = exp(-a * x^2)
217\> quad (f 0.5) 2 5 217>>> quad (f 0.5) 2 5
218(5.7025405463957006e-2,9.678874441303705e-16,95)@ 218(5.7025405463957006e-2,9.678874441303705e-16,95)
219 219
220Unlike other quadrature methods, integrateCQUAD also returns the 220Unlike other quadrature methods, integrateCQUAD also returns the
221number of function evaluations required. 221number of function evaluations required.
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
index 7ccca81..1879dab 100644
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -16,15 +16,15 @@ Minimization of a multidimensional function using some of the algorithms describ
16The example in the GSL manual: 16The example in the GSL manual:
17 17
18@ 18@
19
20f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30 19f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
21 20
22main = do 21main = do
23 let (s,p) = minimize NMSimplex2 1E-2 30 [1,1] f [5,7] 22 let (s,p) = minimize NMSimplex2 1E-2 30 [1,1] f [5,7]
24 print s 23 print s
25 print p 24 print p
25@
26 26
27\> main 27>>> main
28[0.9920430849306288,1.9969168063253182] 28[0.9920430849306288,1.9969168063253182]
29 0.000 512.500 1.130 6.500 5.000 29 0.000 512.500 1.130 6.500 5.000
30 1.000 290.625 1.409 5.250 4.000 30 1.000 290.625 1.409 5.250 4.000
@@ -33,7 +33,6 @@ main = do
33 ... 33 ...
3422.000 30.001 0.013 0.992 1.997 3422.000 30.001 0.013 0.992 1.997
3523.000 30.001 0.008 0.992 1.997 3523.000 30.001 0.008 0.992 1.997
36@
37 36
38The path to the solution can be graphically shown by means of: 37The path to the solution can be graphically shown by means of:
39 38
diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs
index ab037bd..9a29085 100644
--- a/lib/Numeric/GSL/ODE.hs
+++ b/lib/Numeric/GSL/ODE.hs
@@ -13,9 +13,10 @@ Solution of ordinary differential equation (ODE) initial value problems.
13 13
14A simple example: 14A simple example:
15 15
16@import Numeric.GSL 16@
17import Numeric.GSL.ODE
17import Numeric.LinearAlgebra 18import Numeric.LinearAlgebra
18import Graphics.Plot 19import Numeric.LinearAlgebra.Util(mplot)
19 20
20xdot t [x,v] = [v, -0.95*x - 0.1*v] 21xdot t [x,v] = [v, -0.95*x - 0.1*v]
21 22
@@ -23,7 +24,8 @@ ts = linspace 100 (0,20 :: Double)
23 24
24sol = odeSolve xdot [10,0] ts 25sol = odeSolve xdot [10,0] ts
25 26
26main = mplot (ts : toColumns sol)@ 27main = mplot (ts : toColumns sol)
28@
27 29
28-} 30-}
29----------------------------------------------------------------------------- 31-----------------------------------------------------------------------------
diff --git a/lib/Numeric/GSL/Polynomials.hs b/lib/Numeric/GSL/Polynomials.hs
index 694c003..290c615 100644
--- a/lib/Numeric/GSL/Polynomials.hs
+++ b/lib/Numeric/GSL/Polynomials.hs
@@ -27,22 +27,22 @@ import System.IO.Unsafe (unsafePerformIO)
27import Foreign.C.Types (CInt(..)) 27import Foreign.C.Types (CInt(..))
28#endif 28#endif
29 29
30{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/. For example, 30{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/.
31 the three solutions of x^3 + 8 = 0 31
32For example, the three solutions of x^3 + 8 = 0
33
34>>> polySolve [8,0,0,1]
35[(-2.0) :+ 0.0,1.0 :+ 1.7320508075688776,1.0 :+ (-1.7320508075688776)]
32 36
33@\> polySolve [8,0,0,1]
34[(-1.9999999999999998) :+ 0.0,
35 1.0 :+ 1.732050807568877,
36 1.0 :+ (-1.732050807568877)]@
37 37
38The example in the GSL manual: To find the roots of x^5 -1 = 0: 38The example in the GSL manual: To find the roots of x^5 -1 = 0:
39 39
40@\> polySolve [-1, 0, 0, 0, 0, 1] 40>>> polySolve [-1, 0, 0, 0, 0, 1]
41[(-0.8090169943749475) :+ 0.5877852522924731, 41[(-0.8090169943749472) :+ 0.5877852522924731,
42(-0.8090169943749475) :+ (-0.5877852522924731), 42(-0.8090169943749472) :+ (-0.5877852522924731),
430.30901699437494734 :+ 0.9510565162951536, 430.30901699437494756 :+ 0.9510565162951535,
440.30901699437494734 :+ (-0.9510565162951536), 440.30901699437494756 :+ (-0.9510565162951535),
451.0 :+ 0.0]@ 451.0000000000000002 :+ 0.0]
46 46
47-} 47-}
48polySolve :: [Double] -> [Complex Double] 48polySolve :: [Double] -> [Complex Double]
diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs
index 6da15e5..9d561c4 100644
--- a/lib/Numeric/GSL/Root.hs
+++ b/lib/Numeric/GSL/Root.hs
@@ -13,33 +13,23 @@ Multidimensional root finding.
13 13
14The example in the GSL manual: 14The example in the GSL manual:
15 15
16@import Numeric.GSL 16>>> let rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
17import Numeric.LinearAlgebra(format) 17>>> let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
18import Text.Printf(printf) 18>>> sol
19
20rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
21
22disp = putStrLn . format \" \" (printf \"%.3f\")
23
24main = do
25 let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
26 print sol
27 disp path
28
29\> main
30[1.0,1.0] 19[1.0,1.0]
31 0.000 -10.000 -5.000 11.000 -1050.000 20>>> disp 3 path
32 1.000 -3.976 24.827 4.976 90.203 2111x5
22 1.000 -10.000 -5.000 11.000 -1050.000
33 2.000 -3.976 24.827 4.976 90.203 23 2.000 -3.976 24.827 4.976 90.203
34 3.000 -3.976 24.827 4.976 90.203 24 3.000 -3.976 24.827 4.976 90.203
35 4.000 -1.274 -5.680 2.274 -73.018 25 4.000 -3.976 24.827 4.976 90.203
36 5.000 -1.274 -5.680 2.274 -73.018 26 5.000 -1.274 -5.680 2.274 -73.018
37 6.000 0.249 0.298 0.751 2.359 27 6.000 -1.274 -5.680 2.274 -73.018
38 7.000 0.249 0.298 0.751 2.359 28 7.000 0.249 0.298 0.751 2.359
39 8.000 1.000 0.878 -0.000 -1.218 29 8.000 0.249 0.298 0.751 2.359
40 9.000 1.000 0.989 -0.000 -0.108 30 9.000 1.000 0.878 -0.000 -1.218
4110.000 1.000 1.000 0.000 0.000 3110.000 1.000 0.989 -0.000 -0.108
42@ 3211.000 1.000 1.000 0.000 0.000
43 33
44-} 34-}
45----------------------------------------------------------------------------- 35-----------------------------------------------------------------------------
diff --git a/lib/Numeric/IO.hs b/lib/Numeric/IO.hs
index dacfa8b..57275ac 100644
--- a/lib/Numeric/IO.hs
+++ b/lib/Numeric/IO.hs
@@ -39,29 +39,27 @@ format sep f m = table sep . map (map f) . toLists $ m
39 39
40{- | Show a matrix with \"autoscaling\" and a given number of decimal places. 40{- | Show a matrix with \"autoscaling\" and a given number of decimal places.
41 41
42@disp = putStr . disps 2 42>>> putStr . disps 2 $ 120 * (3><4) [1..]
43
44\> disp $ 120 * (3><4) [1..]
453x4 E3 433x4 E3
46 0.12 0.24 0.36 0.48 44 0.12 0.24 0.36 0.48
47 0.60 0.72 0.84 0.96 45 0.60 0.72 0.84 0.96
48 1.08 1.20 1.32 1.44 46 1.08 1.20 1.32 1.44
49@ 47
50-} 48-}
51disps :: Int -> Matrix Double -> String 49disps :: Int -> Matrix Double -> String
52disps d x = sdims x ++ " " ++ formatScaled d x 50disps d x = sdims x ++ " " ++ formatScaled d x
53 51
54{- | Show a matrix with a given number of decimal places. 52{- | Show a matrix with a given number of decimal places.
55 53
56@disp = putStr . dispf 3 54>>> dispf 2 (1/3 + ident 3)
55"3x3\n1.33 0.33 0.33\n0.33 1.33 0.33\n0.33 0.33 1.33\n"
56
57>>> putStr . dispf 2 $ (3><4)[1,1.5..]
583x4
591.00 1.50 2.00 2.50
603.00 3.50 4.00 4.50
615.00 5.50 6.00 6.50
57 62
58\> disp (1/3 + ident 4)
594x4
601.333 0.333 0.333 0.333
610.333 1.333 0.333 0.333
620.333 0.333 1.333 0.333
630.333 0.333 0.333 1.333
64@
65-} 63-}
66dispf :: Int -> Matrix Double -> String 64dispf :: Int -> Matrix Double -> String
67dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x 65dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x
@@ -81,11 +79,9 @@ formatScaled dec t = "E"++show o++"\n" ++ ss
81 79
82{- | Show a vector using a function for showing matrices. 80{- | Show a vector using a function for showing matrices.
83 81
84@disp = putStr . vecdisp ('dispf' 2) 82>>> putStr . vecdisp (dispf 2) $ linspace 10 (0,1)
85
86\> disp ('linspace' 10 (0,1))
8710 |> 0.00 0.11 0.22 0.33 0.44 0.56 0.67 0.78 0.89 1.00 8310 |> 0.00 0.11 0.22 0.33 0.44 0.56 0.67 0.78 0.89 1.00
88@ 84
89-} 85-}
90vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String 86vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String
91vecdisp f v 87vecdisp f v
@@ -94,7 +90,12 @@ vecdisp f v
94 . f . trans . reshape 1 90 . f . trans . reshape 1
95 $ v 91 $ v
96 92
97-- | Tool to display matrices with latex syntax. 93{- | Tool to display matrices with latex syntax.
94
95>>> latexFormat "bmatrix" (dispf 2 $ ident 2)
96"\\begin{bmatrix}\n1 & 0\n\\\\\n0 & 1\n\\end{bmatrix}"
97
98-}
98latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc. 99latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc.
99 -> String -- ^ Formatted matrix, with elements separated by spaces and newlines 100 -> String -- ^ Formatted matrix, with elements separated by spaces and newlines
100 -> String 101 -> String
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 7223cd9..7c1c032 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -371,19 +371,21 @@ pinv = pinvTol 1
371 371
372{- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value. 372{- | @pinvTol r@ computes the pseudoinverse of a matrix with tolerance @tol=r*g*eps*(max rows cols)@, where g is the greatest singular value.
373 373
374@\> let m = 'fromLists' [[1,0, 0] 374@
375 ,[0,1, 0] 375m = (3><3) [ 1, 0, 0
376 ,[0,0,1e-10]] 376 , 0, 1, 0
377\ -- 377 , 0, 0, 1e-10] :: Matrix Double
378\> 'pinv' m 378@
379
380>>> pinv m
3791. 0. 0. 3811. 0. 0.
3800. 1. 0. 3820. 1. 0.
3810. 0. 10000000000. 3830. 0. 10000000000.
382\ -- 384
383\> pinvTol 1E8 m 385>>> pinvTol 1E8 m
3841. 0. 0. 3861. 0. 0.
3850. 1. 0. 3870. 1. 0.
3860. 0. 1.@ 3880. 0. 1.
387 389
388-} 390-}
389 391
@@ -598,10 +600,11 @@ It only works with invertible matrices that have a real solution. For diagonaliz
598@m = (2><2) [4,9 600@m = (2><2) [4,9
599 ,0,4] :: Matrix Double@ 601 ,0,4] :: Matrix Double@
600 602
601@\>sqrtm m 603>>> sqrtm m
602(2><2) 604(2><2)
603 [ 2.0, 2.25 605 [ 2.0, 2.25
604 , 0.0, 2.0 ]@ 606 , 0.0, 2.0 ]
607
605-} 608-}
606sqrtm :: Field t => Matrix t -> Matrix t 609sqrtm :: Field t => Matrix t -> Matrix t
607sqrtm = sqrtmInv 610sqrtm = sqrtmInv
diff --git a/lib/Numeric/LinearAlgebra/Data.hs b/lib/Numeric/LinearAlgebra/Data.hs
index 6e5dcef..a3639d5 100644
--- a/lib/Numeric/LinearAlgebra/Data.hs
+++ b/lib/Numeric/LinearAlgebra/Data.hs
@@ -12,6 +12,8 @@ Stability : provisional
12 12
13module Numeric.LinearAlgebra.Data( 13module Numeric.LinearAlgebra.Data(
14 -- * Vector 14 -- * Vector
15 -- | 1D arrays are storable vectors from the vector package.
16
15 Vector, (|>), dim, (@>), 17 Vector, (|>), dim, (@>),
16 18
17 -- * Matrix 19 -- * Matrix
@@ -49,13 +51,13 @@ module Numeric.LinearAlgebra.Data(
49 51
50 -- * Products 52 -- * Products
51 53
52 (<>), (·), scale, outer, kronecker, cross, 54 (<>), (·), outer, kronecker, cross,
53 sumElements, prodElements, absSum, 55 sumElements, prodElements, absSum,
54 optimiseMult, 56 optimiseMult,
55 57
56 corr, conv, corrMin, corr2, conv2, 58 corr, conv, corrMin, corr2, conv2,
57 59
58 LSDiv(..), 60 (<\>),
59 61
60 -- * Random arrays 62 -- * Random arrays
61 63
@@ -80,7 +82,8 @@ module Numeric.LinearAlgebra.Data(
80 module Data.Complex, 82 module Data.Complex,
81 83
82 -- * Misc 84 -- * Misc
83 meanCov, arctan2, 85 scale, meanCov, arctan2,
86 rows, cols,
84 separable, 87 separable,
85 fromArray2D 88 fromArray2D
86 89
diff --git a/lib/Numeric/LinearAlgebra/Data/Devel.hs b/lib/Numeric/LinearAlgebra/Data/Devel.hs
index 6358d1f..88c980c 100644
--- a/lib/Numeric/LinearAlgebra/Data/Devel.hs
+++ b/lib/Numeric/LinearAlgebra/Data/Devel.hs
@@ -51,13 +51,13 @@ module Numeric.LinearAlgebra.Data.Devel(
51 liftMatrix, liftMatrix2, liftMatrix2Auto, 51 liftMatrix, liftMatrix2, liftMatrix2Auto,
52 52
53 -- * Misc 53 -- * Misc
54 Element, Container 54 Element, Container, Product, Contraction, LSDiv
55) where 55) where
56 56
57import Data.Packed.Foreign 57import Data.Packed.Foreign
58import Data.Packed.Development 58import Data.Packed.Development
59import Data.Packed.ST 59import Data.Packed.ST
60import Numeric.Container(Container) 60import Numeric.Container(Container,Contraction,LSDiv,Product)
61import Data.Packed 61import Data.Packed
62 62
63 63
diff --git a/lib/Numeric/LinearAlgebra/Util.hs b/lib/Numeric/LinearAlgebra/Util.hs
index be01bc1..21b6188 100644
--- a/lib/Numeric/LinearAlgebra/Util.hs
+++ b/lib/Numeric/LinearAlgebra/Util.hs
@@ -13,7 +13,7 @@ Stability : provisional
13{-# OPTIONS_HADDOCK hide #-} 13{-# OPTIONS_HADDOCK hide #-}
14 14
15module Numeric.LinearAlgebra.Util( 15module Numeric.LinearAlgebra.Util(
16 16
17 -- * Convenience functions 17 -- * Convenience functions
18 size, disp, 18 size, disp,
19 zeros, ones, 19 zeros, ones,
@@ -65,8 +65,16 @@ import Numeric.LinearAlgebra.Util.Convolution
65import Graphics.Plot 65import Graphics.Plot
66 66
67 67
68{- | print a real matrix with given number of digits after the decimal point
69
70>>> disp 5 $ ident 2 / 3
712x2
720.33333 0.00000
730.00000 0.33333
74
75-}
68disp :: Int -> Matrix Double -> IO () 76disp :: Int -> Matrix Double -> IO ()
69-- ^ show a matrix with given number of digits after the decimal point 77
70disp n = putStrLn . dispf n 78disp n = putStrLn . dispf n
71 79
72-- | pseudorandom matrix with uniform elements between 0 and 1 80-- | pseudorandom matrix with uniform elements between 0 and 1
@@ -82,7 +90,16 @@ randm d r c = do
82rand :: Int -> Int -> IO (Matrix Double) 90rand :: Int -> Int -> IO (Matrix Double)
83rand = randm Uniform 91rand = randm Uniform
84 92
85-- | pseudorandom matrix with normal elements 93{- | pseudorandom matrix with normal elements
94
95>>> x <- randn 3 5
96>>> disp 3 x
973x5
980.386 -1.141 0.491 -0.510 1.512
990.069 -0.919 1.022 -0.181 0.745
1000.313 -0.670 -0.097 -1.575 -0.583
101
102-}
86randn :: Int -> Int -> IO (Matrix Double) 103randn :: Int -> Int -> IO (Matrix Double)
87randn = randm Gaussian 104randn = randm Gaussian
88 105
@@ -115,9 +132,17 @@ infixl 3 &
115(&) :: Vector Double -> Vector Double -> Vector Double 132(&) :: Vector Double -> Vector Double -> Vector Double
116a & b = vjoin [a,b] 133a & b = vjoin [a,b]
117 134
118-- | horizontal concatenation of real matrices 135{- | horizontal concatenation of real matrices
119-- 136
120-- (0x00a6 broken bar) 137 (0x00a6 broken bar)
138
139>>> ident 3 ¦ konst 7 (3,4)
140(3><7)
141 [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0
142 , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0
143 , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ]
144
145-}
121infixl 3 ¦ 146infixl 3 ¦
122(¦) :: Matrix Double -> Matrix Double -> Matrix Double 147(¦) :: Matrix Double -> Matrix Double -> Matrix Double
123a ¦ b = fromBlocks [[a,b]] 148a ¦ b = fromBlocks [[a,b]]
@@ -141,7 +166,15 @@ row = asRow . fromList
141col :: [Double] -> Matrix Double 166col :: [Double] -> Matrix Double
142col = asColumn . fromList 167col = asColumn . fromList
143 168
144-- | extract selected rows 169{- | extract selected rows
170
171>>> (20><4) [1..] ? [2,1,1]
172(3><4)
173 [ 9.0, 10.0, 11.0, 12.0
174 , 5.0, 6.0, 7.0, 8.0
175 , 5.0, 6.0, 7.0, 8.0 ]
176
177-}
145infixl 9 ? 178infixl 9 ?
146(?) :: Element t => Matrix t -> [Int] -> Matrix t 179(?) :: Element t => Matrix t -> [Int] -> Matrix t
147(?) = flip extractRows 180(?) = flip extractRows
@@ -172,7 +205,7 @@ norm = pnorm PNorm2
172unitary :: Vector Double -> Vector Double 205unitary :: Vector Double -> Vector Double
173unitary v = v / scalar (norm v) 206unitary v = v / scalar (norm v)
174 207
175-- | (rows &&& cols) 208-- | ('rows' &&& 'cols')
176size :: Matrix t -> (Int, Int) 209size :: Matrix t -> (Int, Int)
177size m = (rows m, cols m) 210size m = (rows m, cols m)
178 211