diff options
author | Alberto Ruiz <aruiz@um.es> | 2014-04-30 16:10:20 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2014-04-30 16:10:20 +0200 |
commit | cf1379fed234cf99b2ccce6d9311bfc5c58ef4a3 (patch) | |
tree | 2532837500eb98af048404e33018fa88d8a7f535 /lib | |
parent | b8b7b47450a1007d4d8c77bb7d04e2744c3bfcd4 (diff) |
improved documentation
Diffstat (limited to 'lib')
-rw-r--r-- | lib/Data/Packed/Internal/Matrix.hs | 23 | ||||
-rw-r--r-- | lib/Data/Packed/Internal/Vector.hs | 26 | ||||
-rw-r--r-- | lib/Data/Packed/Matrix.hs | 81 | ||||
-rw-r--r-- | lib/Numeric/Container.hs | 65 | ||||
-rw-r--r-- | lib/Numeric/ContainerBoot.hs | 51 | ||||
-rw-r--r-- | lib/Numeric/GSL/Differentiation.hs | 10 | ||||
-rw-r--r-- | lib/Numeric/GSL/Fitting.hs | 10 | ||||
-rw-r--r-- | lib/Numeric/GSL/Fourier.hs | 4 | ||||
-rw-r--r-- | lib/Numeric/GSL/Integration.hs | 74 | ||||
-rw-r--r-- | lib/Numeric/GSL/Minimization.hs | 5 | ||||
-rw-r--r-- | lib/Numeric/GSL/ODE.hs | 8 | ||||
-rw-r--r-- | lib/Numeric/GSL/Polynomials.hs | 24 | ||||
-rw-r--r-- | lib/Numeric/GSL/Root.hs | 34 | ||||
-rw-r--r-- | lib/Numeric/IO.hs | 35 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 23 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Data.hs | 9 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Data/Devel.hs | 4 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Util.hs | 49 |
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) |
139 | fromList [1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0] | ||
140 | |||
141 | -} | ||
139 | flatten :: Element t => Matrix t -> Vector t | 142 | flatten :: Element t => Matrix t -> Vector t |
140 | flatten = xdat . cmat | 143 | flatten = xdat . cmat |
141 | 144 | ||
145 | {- | ||
142 | type Mt t s = Int -> Int -> Ptr t -> s | 146 | type Mt t s = Int -> Int -> Ptr t -> s |
143 | -- not yet admitted by my haddock version | 147 | |
144 | -- infixr 6 ::> | 148 | infixr 6 ::> |
145 | -- type t ::> s = Mt t s | 149 | type t ::> s = Mt t s |
150 | -} | ||
146 | 151 | ||
147 | -- | the inverse of 'Data.Packed.Matrix.fromLists' | 152 | -- | the inverse of 'Data.Packed.Matrix.fromLists' |
148 | toLists :: (Element t) => Matrix t -> [[t]] | 153 | toLists :: (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@ |
208 | where r is the desired number of rows.) | 213 | where 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 | -} |
217 | reshape :: Storable t => Int -> Vector t -> Matrix t | 222 | reshape :: 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 | -} |
120 | toList :: Storable a => Vector a -> [a] | 120 | toList :: Storable a => Vector a -> [a] |
121 | toList v = safeRead v $ peekArray (dim v) | 121 | toList 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..] | ||
128 | fromList [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 |
128 | infixl 9 |> | 132 | infixl 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]) |
163 | 3 |> [3.0,4.0,5.0]@ | 167 | fromList [3.0,4.0,5.0] |
164 | 168 | ||
165 | -} | 169 | -} |
166 | subVector :: Storable t => Int -- ^ index of the starting element | 170 | subVector :: 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 |
176 | 7.0@ | 180 | 7.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] |
187 | 8 |> [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0]@ | 191 | fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0] |
188 | 192 | ||
189 | -} | 193 | -} |
190 | vjoin :: Storable t => [Vector t] -> Vector t | 194 | vjoin :: 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 | -} |
212 | takesV :: Storable t => [Int] -> Vector t -> [Vector t] | 216 | takesV :: 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 | |||
116 | corresponding common row and column: | 116 | corresponding common row and column: |
117 | 117 | ||
118 | @ | 118 | @ |
119 | > let disp = putStr . dispf 2 | 119 | disp = 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]] | ||
126 | 8x10 | 123 | 8x10 |
127 | 1 0 0 0 0 7 7 7 10 20 | 124 | 1 0 0 0 0 7 7 7 10 20 |
128 | 0 1 0 0 0 7 7 7 10 20 | 125 | 0 1 0 0 0 7 7 7 10 20 |
@@ -132,7 +129,6 @@ corresponding common row and column: | |||
132 | 3 3 3 3 3 1 0 0 0 0 | 129 | 3 3 3 3 3 1 0 0 0 0 |
133 | 3 3 3 3 3 0 2 0 0 0 | 130 | 3 3 3 3 3 0 2 0 0 0 |
134 | 3 3 3 3 3 0 0 3 0 0 | 131 | 3 3 3 3 3 0 0 3 0 0 |
135 | @ | ||
136 | 132 | ||
137 | -} | 133 | -} |
138 | fromBlocks :: Element t => [[Matrix t]] -> Matrix t | 134 | fromBlocks :: 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]] | ||
165 | 7x8 | ||
166 | 1 1 0 0 0 0 0 0 | ||
167 | 1 1 0 0 0 0 0 0 | ||
168 | 0 0 2 2 2 2 2 0 | ||
169 | 0 0 2 2 2 2 2 0 | ||
170 | 0 0 2 2 2 2 2 0 | ||
171 | 0 0 0 0 0 0 0 5 | ||
172 | 0 0 0 0 0 0 0 7 | ||
173 | |||
174 | -} | ||
167 | diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t | 175 | diagBlock :: (Element t, Num t) => [Matrix t] -> Matrix t |
168 | diagBlock ms = fromBlocks $ zipWith f ms [0..] | 176 | diagBlock 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 | -} |
196 | diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t | 205 | diagRect :: (Storable t) => t -> Vector t -> Int -> Int -> Matrix t |
197 | diagRect z v r c = ST.runSTMatrix $ do | 206 | diagRect 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 | ||
216 | This is the format produced by the instances of Show (Matrix a), which | 225 | This is the format produced by the instances of Show (Matrix a), which |
217 | can also be used for input. | 226 | can also be used for input. |
@@ -219,12 +228,11 @@ can also be used for input. | |||
219 | The input list is explicitly truncated, so that it can | 228 | The input list is explicitly truncated, so that it can |
220 | safely be used with lists that are too long (like infinite lists). | 229 | safely be used with lists that are too long (like infinite lists). |
221 | 230 | ||
222 | Example: | 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 | -} |
262 | fromLists :: Element t => [[t]] -> Matrix t | 271 | fromLists :: Element t => [[t]] -> Matrix t |
263 | fromLists = fromRows . map fromList | 272 | fromLists = 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 | -- | ||
266 | asRow :: Storable a => Vector a -> Matrix a | 280 | asRow :: Storable a => Vector a -> Matrix a |
267 | asRow v = reshape (dim v) v | 281 | asRow 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 | -- | ||
270 | asColumn :: Storable a => Vector a -> Matrix a | 293 | asColumn :: Storable a => Vector a -> Matrix a |
271 | asColumn v = reshape 1 v | 294 | asColumn 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 | -} |
318 | repmat :: (Element t) => Matrix t -> Int -> Int -> Matrix t | 341 | repmat :: (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 ..]) |
380 | m[0,0] = 1 | 403 | m[0,0] = 1 |
381 | m[0,1] = 2 | 404 | m[0,1] = 2 |
382 | m[0,2] = 3 | 405 | m[0,2] = 3 |
383 | m[1,0] = 4 | 406 | m[1,0] = 4 |
384 | m[1,1] = 5 | 407 | m[1,1] = 5 |
385 | m[1,2] = 6@ | 408 | m[1,2] = 6 |
409 | |||
386 | -} | 410 | -} |
387 | mapMatrixWithIndexM_ | 411 | mapMatrixWithIndexM_ |
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) |
397 | Just (3><3) | 421 | Just (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 | -} |
402 | mapMatrixWithIndexM | 427 | mapMatrixWithIndexM |
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 | -} |
416 | mapMatrixWithIndex | 443 | mapMatrixWithIndex |
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) |
89 | 5 |> [-3.0,-0.5,2.0,4.5,7.0]@ | 89 | fromList [-3.0,-0.5,2.0,4.5,7.0]@ |
90 | 90 | ||
91 | Logarithmic spacing can be defined as follows: | 91 | Logarithmic spacing can be defined as follows: |
92 | 92 | ||
@@ -105,7 +105,32 @@ cdot u v = udot (conj u) v | |||
105 | class Contraction a b c | a b -> c, a c -> b, b c -> a | 105 | class 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) | ||
119 | 3x3 | ||
120 | 30 70 110 | ||
121 | 70 174 278 | ||
122 | 110 278 446 | ||
123 | |||
124 | >>> a <> fromList [1,0,2,-1::Double] | ||
125 | fromList [3.0,11.0,19.0] | ||
126 | |||
127 | >>> fromList [1,2,3::Double] <> a | ||
128 | fromList [38.0,44.0,50.0,56.0] | ||
129 | |||
130 | >>> fromList [1,i] <> fromList[2*i+1,3] | ||
131 | 1.0 :+ 5.0 | ||
132 | |||
133 | -} | ||
109 | (<>) :: a -> b -> c | 134 | (<>) :: a -> b -> c |
110 | 135 | ||
111 | instance Product t => Contraction (Vector t) (Vector t) t where | 136 | instance 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] | ||
168 | 1.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 |
142 | infixl 7 · | 172 | infixl 7 · |
143 | u · v = cdot u v | 173 | u · v = cdot u v |
@@ -147,6 +177,16 @@ u · v = cdot u v | |||
147 | -- bidirectional type inference | 177 | -- bidirectional type inference |
148 | class Konst e d c | d -> c, c -> d | 178 | class 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 | ||
152 | instance Container Vector e => Konst e Int Vector | 192 | instance Container Vector e => Konst e Int Vector |
@@ -161,6 +201,19 @@ instance Container Vector e => Konst e (Int,Int) Matrix | |||
161 | 201 | ||
162 | class Build d f c e | d -> c, c -> d, f -> e, f -> d, f -> c, c e -> f, d e -> f | 202 | class 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 | ||
166 | instance Container Vector e => Build Int (e -> e) Vector e | 219 | instance 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 |
67 | class (Complexable c, Fractional e, Element e) => Container c e where | 67 | class (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 | -} |
391 | outer :: (Product t) => Vector t -> Vector t -> Matrix t | 411 | outer :: (Product t) => Vector t -> Vector t -> Matrix t |
392 | outer u v = asColumn u `multiply` asRow v | 412 | outer 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 | -} |
416 | kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t | 437 | kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t |
417 | kronecker a b = fromBlocks | 438 | kronecker 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 | 62 | 0.7071067811865476 |
63 | 63 | ||
64 | -} | 64 | -} |
65 | derivCentral :: Double -- ^ initial step size | 65 | derivCentral :: 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 | ||
14 | The example program in the GSL manual (see examples/fitting.hs): | 14 | The example program in the GSL manual (see examples/fitting.hs): |
15 | 15 | ||
16 | @dat = [ | 16 | @ |
17 | dat = [ | ||
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] | |||
25 | expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]] | 26 | expModelDer [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]) |
39 | vector (4) [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]@ | 39 | fromList [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)] |
40 | 40 | ||
41 | -} | 41 | -} |
42 | fft :: Vector (Complex Double) -> Vector (Complex Double) | 42 | fft :: 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 | -} |
38 | foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) | 38 | foreign 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 | ||
50 | integrateQAGS :: Double -- ^ precision (e.g. 1E-9) | 50 | integrateQAGS :: Double -- ^ precision (e.g. 1E-9) |
@@ -56,7 +56,7 @@ integrateQAGS :: Double -- ^ precision (e.g. 1E-9) | |||
56 | integrateQAGS prec n f a b = unsafePerformIO $ do | 56 | integrateQAGS 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 | ||
82 | integrateQNG :: Double -- ^ precision (e.g. 1E-9) | 82 | integrateQNG :: Double -- ^ precision (e.g. 1E-9) |
@@ -87,7 +87,7 @@ integrateQNG :: Double -- ^ precision (e.g. 1E-9) | |||
87 | integrateQNG prec f a b = unsafePerformIO $ do | 87 | integrateQNG 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). |
106 | For example: | 106 | For 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 | ||
115 | integrateQAGI :: Double -- ^ precision (e.g. 1E-9) | 115 | integrateQAGI :: Double -- ^ precision (e.g. 1E-9) |
@@ -119,7 +119,7 @@ integrateQAGI :: Double -- ^ precision (e.g. 1E-9) | |||
119 | integrateQAGI prec n f = unsafePerformIO $ do | 119 | integrateQAGI 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). |
138 | For example: | 138 | For 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 | ||
147 | integrateQAGIU :: Double -- ^ precision (e.g. 1E-9) | 147 | integrateQAGIU :: Double -- ^ precision (e.g. 1E-9) |
@@ -152,7 +152,7 @@ integrateQAGIU :: Double -- ^ precision (e.g. 1E-9) | |||
152 | integrateQAGIU prec n f a = unsafePerformIO $ do | 152 | integrateQAGIU 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). |
171 | For example: | 171 | For 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 | ||
180 | integrateQAGIL :: Double -- ^ precision (e.g. 1E-9) | 180 | integrateQAGIL :: Double -- ^ precision (e.g. 1E-9) |
@@ -185,7 +185,7 @@ integrateQAGIL :: Double -- ^ precision (e.g. 1E-9) | |||
185 | integrateQAGIL prec n f b = unsafePerformIO $ do | 185 | integrateQAGIL 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 | ||
213 | For example: | 213 | For 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 | ||
220 | Unlike other quadrature methods, integrateCQUAD also returns the | 220 | Unlike other quadrature methods, integrateCQUAD also returns the |
221 | number of function evaluations required. | 221 | number 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 | |||
16 | The example in the GSL manual: | 16 | The example in the GSL manual: |
17 | 17 | ||
18 | @ | 18 | @ |
19 | |||
20 | f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30 | 19 | f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30 |
21 | 20 | ||
22 | main = do | 21 | main = 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 | ... |
34 | 22.000 30.001 0.013 0.992 1.997 | 34 | 22.000 30.001 0.013 0.992 1.997 |
35 | 23.000 30.001 0.008 0.992 1.997 | 35 | 23.000 30.001 0.008 0.992 1.997 |
36 | @ | ||
37 | 36 | ||
38 | The path to the solution can be graphically shown by means of: | 37 | The 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 | ||
14 | A simple example: | 14 | A simple example: |
15 | 15 | ||
16 | @import Numeric.GSL | 16 | @ |
17 | import Numeric.GSL.ODE | ||
17 | import Numeric.LinearAlgebra | 18 | import Numeric.LinearAlgebra |
18 | import Graphics.Plot | 19 | import Numeric.LinearAlgebra.Util(mplot) |
19 | 20 | ||
20 | xdot t [x,v] = [v, -0.95*x - 0.1*v] | 21 | xdot t [x,v] = [v, -0.95*x - 0.1*v] |
21 | 22 | ||
@@ -23,7 +24,8 @@ ts = linspace 100 (0,20 :: Double) | |||
23 | 24 | ||
24 | sol = odeSolve xdot [10,0] ts | 25 | sol = odeSolve xdot [10,0] ts |
25 | 26 | ||
26 | main = mplot (ts : toColumns sol)@ | 27 | main = 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) | |||
27 | import Foreign.C.Types (CInt(..)) | 27 | import 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 | |
32 | For 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 | ||
38 | The example in the GSL manual: To find the roots of x^5 -1 = 0: | 38 | The 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), |
43 | 0.30901699437494734 :+ 0.9510565162951536, | 43 | 0.30901699437494756 :+ 0.9510565162951535, |
44 | 0.30901699437494734 :+ (-0.9510565162951536), | 44 | 0.30901699437494756 :+ (-0.9510565162951535), |
45 | 1.0 :+ 0.0]@ | 45 | 1.0000000000000002 :+ 0.0] |
46 | 46 | ||
47 | -} | 47 | -} |
48 | polySolve :: [Double] -> [Complex Double] | 48 | polySolve :: [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 | ||
14 | The example in the GSL manual: | 14 | The 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) ] |
17 | import Numeric.LinearAlgebra(format) | 17 | >>> let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5] |
18 | import Text.Printf(printf) | 18 | >>> sol |
19 | |||
20 | rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ] | ||
21 | |||
22 | disp = putStrLn . format \" \" (printf \"%.3f\") | ||
23 | |||
24 | main = 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 | 21 | 11x5 |
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 |
41 | 10.000 1.000 1.000 0.000 0.000 | 31 | 10.000 1.000 0.989 -0.000 -0.108 |
42 | @ | 32 | 11.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..] | ||
45 | 3x4 E3 | 43 | 3x4 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 | -} |
51 | disps :: Int -> Matrix Double -> String | 49 | disps :: Int -> Matrix Double -> String |
52 | disps d x = sdims x ++ " " ++ formatScaled d x | 50 | disps 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..] | ||
58 | 3x4 | ||
59 | 1.00 1.50 2.00 2.50 | ||
60 | 3.00 3.50 4.00 4.50 | ||
61 | 5.00 5.50 6.00 6.50 | ||
57 | 62 | ||
58 | \> disp (1/3 + ident 4) | ||
59 | 4x4 | ||
60 | 1.333 0.333 0.333 0.333 | ||
61 | 0.333 1.333 0.333 0.333 | ||
62 | 0.333 0.333 1.333 0.333 | ||
63 | 0.333 0.333 0.333 1.333 | ||
64 | @ | ||
65 | -} | 63 | -} |
66 | dispf :: Int -> Matrix Double -> String | 64 | dispf :: Int -> Matrix Double -> String |
67 | dispf d x = sdims x ++ "\n" ++ formatFixed (if isInt x then 0 else d) x | 65 | dispf 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)) | ||
87 | 10 |> 0.00 0.11 0.22 0.33 0.44 0.56 0.67 0.78 0.89 1.00 | 83 | 10 |> 0.00 0.11 0.22 0.33 0.44 0.56 0.67 0.78 0.89 1.00 |
88 | @ | 84 | |
89 | -} | 85 | -} |
90 | vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String | 86 | vecdisp :: (Element t) => (Matrix t -> String) -> Vector t -> String |
91 | vecdisp f v | 87 | vecdisp 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 | -} | ||
98 | latexFormat :: String -- ^ type of braces: \"matrix\", \"bmatrix\", \"pmatrix\", etc. | 99 | latexFormat :: 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] | 375 | m = (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 | ||
379 | 1. 0. 0. | 381 | 1. 0. 0. |
380 | 0. 1. 0. | 382 | 0. 1. 0. |
381 | 0. 0. 10000000000. | 383 | 0. 0. 10000000000. |
382 | \ -- | 384 | |
383 | \> pinvTol 1E8 m | 385 | >>> pinvTol 1E8 m |
384 | 1. 0. 0. | 386 | 1. 0. 0. |
385 | 0. 1. 0. | 387 | 0. 1. 0. |
386 | 0. 0. 1.@ | 388 | 0. 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 | -} |
606 | sqrtm :: Field t => Matrix t -> Matrix t | 609 | sqrtm :: Field t => Matrix t -> Matrix t |
607 | sqrtm = sqrtmInv | 610 | sqrtm = 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 | ||
13 | module Numeric.LinearAlgebra.Data( | 13 | module 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 | ||
57 | import Data.Packed.Foreign | 57 | import Data.Packed.Foreign |
58 | import Data.Packed.Development | 58 | import Data.Packed.Development |
59 | import Data.Packed.ST | 59 | import Data.Packed.ST |
60 | import Numeric.Container(Container) | 60 | import Numeric.Container(Container,Contraction,LSDiv,Product) |
61 | import Data.Packed | 61 | import 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 | ||
15 | module Numeric.LinearAlgebra.Util( | 15 | module 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 | |||
65 | import Graphics.Plot | 65 | import 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 | ||
71 | 2x2 | ||
72 | 0.33333 0.00000 | ||
73 | 0.00000 0.33333 | ||
74 | |||
75 | -} | ||
68 | disp :: Int -> Matrix Double -> IO () | 76 | disp :: Int -> Matrix Double -> IO () |
69 | -- ^ show a matrix with given number of digits after the decimal point | 77 | |
70 | disp n = putStrLn . dispf n | 78 | disp 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 | |||
82 | rand :: Int -> Int -> IO (Matrix Double) | 90 | rand :: Int -> Int -> IO (Matrix Double) |
83 | rand = randm Uniform | 91 | rand = 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 | ||
97 | 3x5 | ||
98 | 0.386 -1.141 0.491 -0.510 1.512 | ||
99 | 0.069 -0.919 1.022 -0.181 0.745 | ||
100 | 0.313 -0.670 -0.097 -1.575 -0.583 | ||
101 | |||
102 | -} | ||
86 | randn :: Int -> Int -> IO (Matrix Double) | 103 | randn :: Int -> Int -> IO (Matrix Double) |
87 | randn = randm Gaussian | 104 | randn = randm Gaussian |
88 | 105 | ||
@@ -115,9 +132,17 @@ infixl 3 & | |||
115 | (&) :: Vector Double -> Vector Double -> Vector Double | 132 | (&) :: Vector Double -> Vector Double -> Vector Double |
116 | a & b = vjoin [a,b] | 133 | a & 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 | -} | ||
121 | infixl 3 ¦ | 146 | infixl 3 ¦ |
122 | (¦) :: Matrix Double -> Matrix Double -> Matrix Double | 147 | (¦) :: Matrix Double -> Matrix Double -> Matrix Double |
123 | a ¦ b = fromBlocks [[a,b]] | 148 | a ¦ b = fromBlocks [[a,b]] |
@@ -141,7 +166,15 @@ row = asRow . fromList | |||
141 | col :: [Double] -> Matrix Double | 166 | col :: [Double] -> Matrix Double |
142 | col = asColumn . fromList | 167 | col = 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 | -} | ||
145 | infixl 9 ? | 178 | infixl 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 | |||
172 | unitary :: Vector Double -> Vector Double | 205 | unitary :: Vector Double -> Vector Double |
173 | unitary v = v / scalar (norm v) | 206 | unitary v = v / scalar (norm v) |
174 | 207 | ||
175 | -- | (rows &&& cols) | 208 | -- | ('rows' &&& 'cols') |
176 | size :: Matrix t -> (Int, Int) | 209 | size :: Matrix t -> (Int, Int) |
177 | size m = (rows m, cols m) | 210 | size m = (rows m, cols m) |
178 | 211 | ||