diff options
Diffstat (limited to 'packages/base/src/Numeric')
7 files changed, 89 insertions, 17 deletions
diff --git a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs index 25700bc..02ac6a0 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Algorithms.hs | |||
@@ -809,7 +809,7 @@ expGolub m = iterate msq f !! j | |||
809 | -------------------------------------------------------------- | 809 | -------------------------------------------------------------- |
810 | 810 | ||
811 | {- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia. | 811 | {- | Matrix square root. Currently it uses a simple iterative algorithm described in Wikipedia. |
812 | It only works with invertible matrices that have a real solution. For diagonalizable matrices you can try @matFunc sqrt@. | 812 | It only works with invertible matrices that have a real solution. |
813 | 813 | ||
814 | @m = (2><2) [4,9 | 814 | @m = (2><2) [4,9 |
815 | ,0,4] :: Matrix Double@ | 815 | ,0,4] :: Matrix Double@ |
@@ -819,6 +819,13 @@ It only works with invertible matrices that have a real solution. For diagonaliz | |||
819 | [ 2.0, 2.25 | 819 | [ 2.0, 2.25 |
820 | , 0.0, 2.0 ] | 820 | , 0.0, 2.0 ] |
821 | 821 | ||
822 | For diagonalizable matrices you can try 'matFunc' @sqrt@: | ||
823 | |||
824 | >>> matFunc sqrt ((2><2) [1,0,0,-1]) | ||
825 | (2><2) | ||
826 | [ 1.0 :+ 0.0, 0.0 :+ 0.0 | ||
827 | , 0.0 :+ 0.0, 0.0 :+ 1.0 ] | ||
828 | |||
822 | -} | 829 | -} |
823 | sqrtm :: Field t => Matrix t -> Matrix t | 830 | sqrtm :: Field t => Matrix t -> Matrix t |
824 | sqrtm = sqrtmInv | 831 | sqrtm = sqrtmInv |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Data.hs b/packages/base/src/Numeric/LinearAlgebra/Data.hs index b1a31fc..6dea407 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Data.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Data.hs | |||
@@ -40,7 +40,7 @@ module Numeric.LinearAlgebra.Data( | |||
40 | takeRows, dropRows, takeColumns, dropColumns, subMatrix, (?), (¿), fliprl, flipud, | 40 | takeRows, dropRows, takeColumns, dropColumns, subMatrix, (?), (¿), fliprl, flipud, |
41 | 41 | ||
42 | -- * Block matrix | 42 | -- * Block matrix |
43 | fromBlocks, (¦), (——), diagBlock, repmat, toBlocks, toBlocksEvery, | 43 | fromBlocks, (|||), (===), diagBlock, repmat, toBlocks, toBlocksEvery, |
44 | 44 | ||
45 | -- * Mapping functions | 45 | -- * Mapping functions |
46 | conj, cmap, step, cond, | 46 | conj, cmap, step, cond, |
@@ -66,7 +66,7 @@ module Numeric.LinearAlgebra.Data( | |||
66 | arctan2, | 66 | arctan2, |
67 | rows, cols, | 67 | rows, cols, |
68 | separable, | 68 | separable, |
69 | 69 | (¦),(——), | |
70 | module Data.Complex, | 70 | module Data.Complex, |
71 | 71 | ||
72 | Vector, Matrix, GMatrix, nRows, nCols | 72 | Vector, Matrix, GMatrix, nRows, nCols |
diff --git a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs index d2cae6c..677f9ee 100644 --- a/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs +++ b/packages/base/src/Numeric/LinearAlgebra/HMatrix.hs | |||
@@ -134,7 +134,7 @@ module Numeric.LinearAlgebra.HMatrix ( | |||
134 | Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, | 134 | Seed, RandDist(..), randomVector, rand, randn, gaussianSample, uniformSample, |
135 | 135 | ||
136 | -- * Misc | 136 | -- * Misc |
137 | meanCov, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, | 137 | meanCov, rowOuters, peps, relativeError, haussholder, optimiseMult, udot, nullspaceSVD, orthSVD, ranksv, |
138 | ℝ,ℂ,iC, | 138 | ℝ,ℂ,iC, |
139 | -- * Auxiliary classes | 139 | -- * Auxiliary classes |
140 | Element, Container, Product, Numeric, LSDiv, | 140 | Element, Container, Product, Numeric, LSDiv, |
@@ -194,7 +194,37 @@ mul :: Numeric t => Matrix t -> Matrix t -> Matrix t | |||
194 | mul = mXm | 194 | mul = mXm |
195 | 195 | ||
196 | 196 | ||
197 | -- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. | 197 | {- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, returning Nothing for a singular system. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'. |
198 | |||
199 | @ | ||
200 | a = (2><2) | ||
201 | [ 1.0, 2.0 | ||
202 | , 3.0, 5.0 ] | ||
203 | @ | ||
204 | |||
205 | @ | ||
206 | b = (2><3) | ||
207 | [ 6.0, 1.0, 10.0 | ||
208 | , 15.0, 3.0, 26.0 ] | ||
209 | @ | ||
210 | |||
211 | >>> linearSolve a b | ||
212 | Just (2><3) | ||
213 | [ -1.4802973661668753e-15, 0.9999999999999997, 1.999999999999997 | ||
214 | , 3.000000000000001, 1.6653345369377348e-16, 4.000000000000002 ] | ||
215 | |||
216 | >>> let Just x = it | ||
217 | >>> disp 5 x | ||
218 | 2x3 | ||
219 | -0.00000 1.00000 2.00000 | ||
220 | 3.00000 0.00000 4.00000 | ||
221 | |||
222 | >>> a <> x | ||
223 | (2><3) | ||
224 | [ 6.0, 1.0, 10.0 | ||
225 | , 15.0, 3.0, 26.0 ] | ||
226 | |||
227 | -} | ||
198 | linearSolve m b = A.mbLinearSolve m b | 228 | linearSolve m b = A.mbLinearSolve m b |
199 | 229 | ||
200 | -- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'. | 230 | -- | return an orthonormal basis of the null space of a matrix. See also 'nullspaceSVD'. |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static.hs b/packages/base/src/Numeric/LinearAlgebra/Static.hs index 5749c40..037396d 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static.hs | |||
@@ -512,7 +512,10 @@ crossC (extract -> x) (extract -> y) = mkC (LA.fromList [z1, z2, z3]) | |||
512 | -------------------------------------------------------------------------------- | 512 | -------------------------------------------------------------------------------- |
513 | 513 | ||
514 | diagRectR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℝ -> R k -> L m n | 514 | diagRectR :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℝ -> R k -> L m n |
515 | diagRectR x v = r | 515 | diagRectR x v |
516 | | m' == 1 = mkL (LA.diagRect x ev m' n') | ||
517 | | m'*n' > 0 = r | ||
518 | | otherwise = matrix [] | ||
516 | where | 519 | where |
517 | r = mkL (asRow (vjoin [scalar x, ev, zeros])) | 520 | r = mkL (asRow (vjoin [scalar x, ev, zeros])) |
518 | ev = extract v | 521 | ev = extract v |
@@ -521,7 +524,10 @@ diagRectR x v = r | |||
521 | 524 | ||
522 | 525 | ||
523 | diagRectC :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℂ -> C k -> M m n | 526 | diagRectC :: forall m n k . (KnownNat m, KnownNat n, KnownNat k) => ℂ -> C k -> M m n |
524 | diagRectC x v = r | 527 | diagRectC x v |
528 | | m' == 1 = mkM (LA.diagRect x ev m' n') | ||
529 | | m'*n' > 0 = r | ||
530 | | otherwise = fromList [] | ||
525 | where | 531 | where |
526 | r = mkM (asRow (vjoin [scalar x, ev, zeros])) | 532 | r = mkM (asRow (vjoin [scalar x, ev, zeros])) |
527 | ev = extract v | 533 | ev = extract v |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs b/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs index 339ef7d..ec02cf6 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Static/Internal.hs | |||
@@ -150,7 +150,7 @@ gmat st xs' | |||
150 | (xs,rest) = splitAt (m'*n') xs' | 150 | (xs,rest) = splitAt (m'*n') xs' |
151 | v = LA.fromList xs | 151 | v = LA.fromList xs |
152 | x = reshape n' v | 152 | x = reshape n' v |
153 | ok = rem (LA.size v) n' == 0 && LA.size x == (m',n') && null rest | 153 | ok = null rest && ((n' == 0 && dim v == 0) || n'> 0 && (rem (LA.size v) n' == 0) && LA.size x == (m',n')) |
154 | m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int | 154 | m' = fromIntegral . natVal $ (undefined :: Proxy m) :: Int |
155 | n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int | 155 | n' = fromIntegral . natVal $ (undefined :: Proxy n) :: Int |
156 | abort info = error $ st ++" "++show m' ++ " " ++ show n'++" can't be created from elements " ++ info | 156 | abort info = error $ st ++" "++show m' ++ " " ++ show n'++" can't be created from elements " ++ info |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util.hs b/packages/base/src/Numeric/LinearAlgebra/Util.hs index 6bb9d15..043aa21 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util.hs | |||
@@ -33,7 +33,7 @@ module Numeric.LinearAlgebra.Util( | |||
33 | diagl, | 33 | diagl, |
34 | row, | 34 | row, |
35 | col, | 35 | col, |
36 | (&), (¦), (——), (#), | 36 | (&), (¦), (|||), (——), (===), (#), |
37 | (?), (¿), | 37 | (?), (¿), |
38 | Indexable(..), size, | 38 | Indexable(..), size, |
39 | Numeric, | 39 | Numeric, |
@@ -157,25 +157,34 @@ a & b = vjoin [a,b] | |||
157 | 157 | ||
158 | {- | horizontal concatenation of real matrices | 158 | {- | horizontal concatenation of real matrices |
159 | 159 | ||
160 | (unicode 0x00a6, broken bar) | 160 | >>> ident 3 ||| konst 7 (3,4) |
161 | |||
162 | >>> ident 3 ¦ konst 7 (3,4) | ||
163 | (3><7) | 161 | (3><7) |
164 | [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0 | 162 | [ 1.0, 0.0, 0.0, 7.0, 7.0, 7.0, 7.0 |
165 | , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0 | 163 | , 0.0, 1.0, 0.0, 7.0, 7.0, 7.0, 7.0 |
166 | , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ] | 164 | , 0.0, 0.0, 1.0, 7.0, 7.0, 7.0, 7.0 ] |
167 | 165 | ||
168 | -} | 166 | -} |
167 | infixl 3 ||| | ||
168 | (|||) :: Matrix Double -> Matrix Double -> Matrix Double | ||
169 | a ||| b = fromBlocks [[a,b]] | ||
170 | |||
171 | -- | a synonym for ('|||') (unicode 0x00a6, broken bar) | ||
169 | infixl 3 ¦ | 172 | infixl 3 ¦ |
170 | (¦) :: Matrix Double -> Matrix Double -> Matrix Double | 173 | (¦) :: Matrix Double -> Matrix Double -> Matrix Double |
171 | a ¦ b = fromBlocks [[a,b]] | 174 | (¦) = (|||) |
175 | |||
172 | 176 | ||
173 | -- | vertical concatenation of real matrices | 177 | -- | vertical concatenation of real matrices |
174 | -- | 178 | -- |
175 | -- (unicode 0x2014, em dash) | 179 | (===) :: Matrix Double -> Matrix Double -> Matrix Double |
180 | infixl 2 === | ||
181 | a === b = fromBlocks [[a],[b]] | ||
182 | |||
183 | -- | a synonym for ('===') (unicode 0x2014, em dash) | ||
176 | (——) :: Matrix Double -> Matrix Double -> Matrix Double | 184 | (——) :: Matrix Double -> Matrix Double -> Matrix Double |
177 | infixl 2 —— | 185 | infixl 2 —— |
178 | a —— b = fromBlocks [[a],[b]] | 186 | (——) = (===) |
187 | |||
179 | 188 | ||
180 | (#) :: Matrix Double -> Matrix Double -> Matrix Double | 189 | (#) :: Matrix Double -> Matrix Double -> Matrix Double |
181 | infixl 2 # | 190 | infixl 2 # |
@@ -356,7 +365,26 @@ pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y | |||
356 | 365 | ||
357 | -------------------------------------------------------------------------------- | 366 | -------------------------------------------------------------------------------- |
358 | 367 | ||
359 | -- | outer products of rows | 368 | {- | outer products of rows |
369 | |||
370 | >>> a | ||
371 | (3><2) | ||
372 | [ 1.0, 2.0 | ||
373 | , 10.0, 20.0 | ||
374 | , 100.0, 200.0 ] | ||
375 | >>> b | ||
376 | (3><3) | ||
377 | [ 1.0, 2.0, 3.0 | ||
378 | , 4.0, 5.0, 6.0 | ||
379 | , 7.0, 8.0, 9.0 ] | ||
380 | |||
381 | >>> rowOuters a (b ||| 1) | ||
382 | (3><8) | ||
383 | [ 1.0, 2.0, 3.0, 1.0, 2.0, 4.0, 6.0, 2.0 | ||
384 | , 40.0, 50.0, 60.0, 10.0, 80.0, 100.0, 120.0, 20.0 | ||
385 | , 700.0, 800.0, 900.0, 100.0, 1400.0, 1600.0, 1800.0, 200.0 ] | ||
386 | |||
387 | -} | ||
360 | rowOuters :: Matrix Double -> Matrix Double -> Matrix Double | 388 | rowOuters :: Matrix Double -> Matrix Double -> Matrix Double |
361 | rowOuters a b = a' * b' | 389 | rowOuters a b = a' * b' |
362 | where | 390 | where |
diff --git a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs index c8c7536..c9e75de 100644 --- a/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs +++ b/packages/base/src/Numeric/LinearAlgebra/Util/Convolution.hs | |||
@@ -16,6 +16,7 @@ module Numeric.LinearAlgebra.Util.Convolution( | |||
16 | corr2, conv2, separable | 16 | corr2, conv2, separable |
17 | ) where | 17 | ) where |
18 | 18 | ||
19 | import qualified Data.Vector.Storable as SV | ||
19 | import Data.Packed.Numeric | 20 | import Data.Packed.Numeric |
20 | 21 | ||
21 | 22 | ||
@@ -51,7 +52,7 @@ conv ker v | |||
51 | | dim ker == 0 = konst 0 (dim v) | 52 | | dim ker == 0 = konst 0 (dim v) |
52 | | otherwise = corr ker' v' | 53 | | otherwise = corr ker' v' |
53 | where | 54 | where |
54 | ker' = (flatten.fliprl.asRow) ker | 55 | ker' = SV.reverse ker |
55 | v' = vjoin [z,v,z] | 56 | v' = vjoin [z,v,z] |
56 | z = konst 0 (dim ker -1) | 57 | z = konst 0 (dim ker -1) |
57 | 58 | ||