diff options
Diffstat (limited to 'packages/base/src/Numeric')
-rw-r--r-- | packages/base/src/Numeric/HMatrix.hs | 41 |
1 files changed, 32 insertions, 9 deletions
diff --git a/packages/base/src/Numeric/HMatrix.hs b/packages/base/src/Numeric/HMatrix.hs index 1b6c1ea..c88918f 100644 --- a/packages/base/src/Numeric/HMatrix.hs +++ b/packages/base/src/Numeric/HMatrix.hs | |||
@@ -47,21 +47,22 @@ module Numeric.HMatrix( | |||
47 | -- * Linear Systems | 47 | -- * Linear Systems |
48 | linSolve, (<\>), | 48 | linSolve, (<\>), |
49 | -- * Factorizations | 49 | -- * Factorizations |
50 | svd, svdTall, svdFlat, Eigen(..), | 50 | svd, withCompactSVD, svdTall, svdFlat, Eigen(..), |
51 | withNullspace, | 51 | withNullspace, qr, |
52 | -- * Misc | 52 | -- * Misc |
53 | mean, | 53 | mean, |
54 | Disp(..), Domain(..), | 54 | Disp(..), Domain(..), |
55 | withVector, withMatrix, | 55 | withVector, withMatrix, |
56 | toRows, toColumns, | 56 | toRows, toColumns, |
57 | Sized(..), Diag(..), Sym, sym | 57 | Sized(..), Diag(..), Sym, sym, mTm, unSym |
58 | ) where | 58 | ) where |
59 | 59 | ||
60 | 60 | ||
61 | import GHC.TypeLits | 61 | import GHC.TypeLits |
62 | import Numeric.LinearAlgebra.HMatrix hiding ( | 62 | import Numeric.LinearAlgebra.HMatrix hiding ( |
63 | (<>),(#>),(<·>),Konst(..),diag, disp,(¦),(——),row,col,vector,matrix,linspace,toRows,toColumns, | 63 | (<>),(#>),(<·>),Konst(..),diag, disp,(¦),(——),row,col,vector,matrix,linspace,toRows,toColumns, |
64 | (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',eigenvalues,eigenvaluesSH,eigenvaluesSH',build) | 64 | (<\>),fromList,takeDiag,svd,eig,eigSH,eigSH',eigenvalues,eigenvaluesSH,eigenvaluesSH',build, |
65 | qr) | ||
65 | import qualified Numeric.LinearAlgebra.HMatrix as LA | 66 | import qualified Numeric.LinearAlgebra.HMatrix as LA |
66 | import Data.Proxy(Proxy) | 67 | import Data.Proxy(Proxy) |
67 | import Numeric.LinearAlgebra.Static | 68 | import Numeric.LinearAlgebra.Static |
@@ -280,6 +281,12 @@ newtype Sym n = Sym (Sq n) deriving Show | |||
280 | sym :: KnownNat n => Sq n -> Sym n | 281 | sym :: KnownNat n => Sq n -> Sym n |
281 | sym m = Sym $ (m + tr m)/2 | 282 | sym m = Sym $ (m + tr m)/2 |
282 | 283 | ||
284 | mTm :: (KnownNat m, KnownNat n) => L m n -> Sym n | ||
285 | mTm x = Sym (tr x <> x) | ||
286 | |||
287 | unSym :: Sym n -> Sq n | ||
288 | unSym (Sym x) = x | ||
289 | |||
283 | 290 | ||
284 | 𝑖 :: Sized ℂ s c => s | 291 | 𝑖 :: Sized ℂ s c => s |
285 | 𝑖 = konst iC | 292 | 𝑖 = konst iC |
@@ -307,7 +314,6 @@ instance KnownNat n => Eigen (Sq n) (C n) (M n n) | |||
307 | 314 | ||
308 | -------------------------------------------------------------------------------- | 315 | -------------------------------------------------------------------------------- |
309 | 316 | ||
310 | |||
311 | withNullspace | 317 | withNullspace |
312 | :: forall m n z . (KnownNat m, KnownNat n) | 318 | :: forall m n z . (KnownNat m, KnownNat n) |
313 | => L m n | 319 | => L m n |
@@ -318,6 +324,26 @@ withNullspace (LA.nullspace . extract -> a) f = | |||
318 | Nothing -> error "static/dynamic mismatch" | 324 | Nothing -> error "static/dynamic mismatch" |
319 | Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k) | 325 | Just (SomeNat (_ :: Proxy k)) -> f (mkL a :: L n k) |
320 | 326 | ||
327 | |||
328 | withCompactSVD | ||
329 | :: forall m n z . (KnownNat m, KnownNat n) | ||
330 | => L m n | ||
331 | -> (forall k . (KnownNat k) => (L m k, R k, L n k) -> z) | ||
332 | -> z | ||
333 | withCompactSVD (LA.compactSVD . extract -> (u,s,v)) f = | ||
334 | case someNatVal $ fromIntegral $ size s of | ||
335 | Nothing -> error "static/dynamic mismatch" | ||
336 | Just (SomeNat (_ :: Proxy k)) -> f (mkL u :: L m k, mkR s :: R k, mkL v :: L n k) | ||
337 | |||
338 | -------------------------------------------------------------------------------- | ||
339 | |||
340 | qr :: (KnownNat m, KnownNat n) => L m n -> (L m m, L m n) | ||
341 | qr (extract -> x) = (mkL q, mkL r) | ||
342 | where | ||
343 | (q,r) = LA.qr x | ||
344 | |||
345 | -- use qrRaw? | ||
346 | |||
321 | -------------------------------------------------------------------------------- | 347 | -------------------------------------------------------------------------------- |
322 | 348 | ||
323 | split :: forall p n . (KnownNat p, KnownNat n, p<=n) => R n -> (R p, R (n-p)) | 349 | split :: forall p n . (KnownNat p, KnownNat n, p<=n) => R n -> (R p, R (n-p)) |
@@ -518,7 +544,7 @@ test :: (Bool, IO ()) | |||
518 | test = (ok,info) | 544 | test = (ok,info) |
519 | where | 545 | where |
520 | ok = extract (eye :: Sq 5) == ident 5 | 546 | ok = extract (eye :: Sq 5) == ident 5 |
521 | && unwrap (mTm sm :: Sq 3) == tr ((3><3)[1..]) LA.<> (3><3)[1..] | 547 | && (unwrap .unSym) (mTm sm :: Sym 3) == tr ((3><3)[1..]) LA.<> (3><3)[1..] |
522 | && unwrap (tm :: L 3 5) == LA.matrix 5 [1..15] | 548 | && unwrap (tm :: L 3 5) == LA.matrix 5 [1..15] |
523 | && thingS == thingD | 549 | && thingS == thingD |
524 | && precS == precD | 550 | && precS == precD |
@@ -546,9 +572,6 @@ test = (ok,info) | |||
546 | 572 | ||
547 | v = 𝕧 2 & 4 & 7 | 573 | v = 𝕧 2 & 4 & 7 |
548 | 574 | ||
549 | -- mTm :: L n m -> Sq m | ||
550 | mTm a = tr a <> a | ||
551 | |||
552 | tm :: GL | 575 | tm :: GL |
553 | tm = lmat 0 [1..] | 576 | tm = lmat 0 [1..] |
554 | 577 | ||