diff options
Diffstat (limited to 'lib')
-rw-r--r-- | lib/Data/Packed/Matrix.hs | 40 | ||||
-rw-r--r-- | lib/Numeric/GSL/Fitting.hs | 2 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 132 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Interface.hs | 55 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Linear.hs | 63 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests.hs | 4 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests/Instances.hs | 20 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests/Properties.hs | 8 |
8 files changed, 188 insertions, 136 deletions
diff --git a/lib/Data/Packed/Matrix.hs b/lib/Data/Packed/Matrix.hs index 9059723..2855a90 100644 --- a/lib/Data/Packed/Matrix.hs +++ b/lib/Data/Packed/Matrix.hs | |||
@@ -21,7 +21,7 @@ | |||
21 | module Data.Packed.Matrix ( | 21 | module Data.Packed.Matrix ( |
22 | Element, RealElement, Container(..), | 22 | Element, RealElement, Container(..), |
23 | Convert(..), RealOf, ComplexOf, SingleOf, DoubleOf, ElementOf, | 23 | Convert(..), RealOf, ComplexOf, SingleOf, DoubleOf, ElementOf, |
24 | Precision(..), comp, real'', complex'', | 24 | Precision(..), comp, |
25 | AutoReal(..), | 25 | AutoReal(..), |
26 | Matrix,rows,cols, | 26 | Matrix,rows,cols, |
27 | (><), | 27 | (><), |
@@ -485,7 +485,10 @@ instance Precision (Complex Float) (Complex Double) where | |||
485 | float2DoubleG = asComplex . float2DoubleV . asReal | 485 | float2DoubleG = asComplex . float2DoubleV . asReal |
486 | 486 | ||
487 | -- | Supported real types | 487 | -- | Supported real types |
488 | class (Element t, Element (Complex t), RealFloat t) => RealElement t | 488 | class (Element t, Element (Complex t), RealFloat t |
489 | -- , RealOf t ~ t, RealOf (Complex t) ~ t | ||
490 | ) | ||
491 | => RealElement t | ||
489 | 492 | ||
490 | instance RealElement Double | 493 | instance RealElement Double |
491 | 494 | ||
@@ -501,6 +504,8 @@ class Container c where | |||
501 | single' :: Precision a b => c b -> c a | 504 | single' :: Precision a b => c b -> c a |
502 | double' :: Precision a b => c a -> c b | 505 | double' :: Precision a b => c a -> c b |
503 | 506 | ||
507 | comp x = complex' x | ||
508 | |||
504 | instance Container Vector where | 509 | instance Container Vector where |
505 | toComplex = toComplexV | 510 | toComplex = toComplexV |
506 | fromComplex = fromComplexV | 511 | fromComplex = fromComplexV |
@@ -592,34 +597,23 @@ instance Convert (Complex Float) where | |||
592 | 597 | ||
593 | ------------------------------------------------------------------- | 598 | ------------------------------------------------------------------- |
594 | 599 | ||
595 | |||
596 | -- | to be replaced by Convert | 600 | -- | to be replaced by Convert |
597 | class Convert t => AutoReal t where | 601 | class Convert t => AutoReal t where |
598 | real''' :: Container c => c Double -> c t | 602 | real'' :: Container c => c Double -> c t |
599 | complex''' :: Container c => c t -> c (Complex Double) | 603 | complex'' :: Container c => c t -> c (Complex Double) |
600 | 604 | ||
601 | instance AutoReal Double where | 605 | instance AutoReal Double where |
602 | real''' = real | 606 | real'' = real |
603 | complex''' = complex | 607 | complex'' = complex |
604 | 608 | ||
605 | instance AutoReal (Complex Double) where | 609 | instance AutoReal (Complex Double) where |
606 | real''' = real | 610 | real'' = real |
607 | complex''' = complex | 611 | complex'' = complex |
608 | 612 | ||
609 | instance AutoReal Float where | 613 | instance AutoReal Float where |
610 | real''' = real . single | 614 | real'' = real . single |
611 | complex''' = double . complex | 615 | complex'' = double . complex |
612 | 616 | ||
613 | instance AutoReal (Complex Float) where | 617 | instance AutoReal (Complex Float) where |
614 | real''' = real . single | 618 | real'' = real . single |
615 | complex''' = double . complex | 619 | complex'' = double . complex |
616 | |||
617 | |||
618 | comp x = complex' x | ||
619 | |||
620 | -- complex'' x = double (complex x) | ||
621 | -- real'' x = real (single x) | ||
622 | |||
623 | real'' x = real''' x | ||
624 | complex'' x = complex''' x | ||
625 | |||
diff --git a/lib/Numeric/GSL/Fitting.hs b/lib/Numeric/GSL/Fitting.hs index b1f3a32..e5aa528 100644 --- a/lib/Numeric/GSL/Fitting.hs +++ b/lib/Numeric/GSL/Fitting.hs | |||
@@ -109,7 +109,7 @@ err (model,deriv) dat vsol = zip sol errs where | |||
109 | sol = toList vsol | 109 | sol = toList vsol |
110 | c = max 1 (chi/sqrt (fromIntegral dof)) | 110 | c = max 1 (chi/sqrt (fromIntegral dof)) |
111 | dof = length dat - (rows cov) | 111 | dof = length dat - (rows cov) |
112 | chi = pnorm PNorm2 (fromList $ cost (resMs model) dat sol) | 112 | chi = norm2 (fromList $ cost (resMs model) dat sol) |
113 | js = fromLists $ jacobian (resDs deriv) dat sol | 113 | js = fromLists $ jacobian (resDs deriv) dat sol |
114 | cov = inv $ trans js <> js | 114 | cov = inv $ trans js <> js |
115 | errs = toList $ scalar c * sqrt (takeDiag cov) | 115 | errs = toList $ scalar c * sqrt (takeDiag cov) |
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index 7e258de..0f3ff0e 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs | |||
@@ -1,5 +1,8 @@ | |||
1 | {-# LANGUAGE FlexibleContexts, FlexibleInstances #-} | 1 | {-# LANGUAGE FlexibleContexts, FlexibleInstances #-} |
2 | {-# LANGUAGE CPP #-} | 2 | {-# LANGUAGE CPP #-} |
3 | {-# LANGUAGE MultiParamTypeClasses #-} | ||
4 | {-# LANGUAGE UndecidableInstances #-} | ||
5 | {-# LANGUAGE TypeFamilies #-} | ||
3 | ----------------------------------------------------------------------------- | 6 | ----------------------------------------------------------------------------- |
4 | {- | | 7 | {- | |
5 | Module : Numeric.LinearAlgebra.Algorithms | 8 | Module : Numeric.LinearAlgebra.Algorithms |
@@ -58,10 +61,9 @@ module Numeric.LinearAlgebra.Algorithms ( | |||
58 | nullspacePrec, | 61 | nullspacePrec, |
59 | nullVector, | 62 | nullVector, |
60 | nullspaceSVD, | 63 | nullspaceSVD, |
61 | -- * Norms | ||
62 | Normed(..), NormType(..), | ||
63 | -- * Misc | 64 | -- * Misc |
64 | eps, i, | 65 | eps, peps, i, |
66 | Normed(..), NormType(..), | ||
65 | -- * Util | 67 | -- * Util |
66 | haussholder, | 68 | haussholder, |
67 | unpackQR, unpackHess, | 69 | unpackQR, unpackHess, |
@@ -72,17 +74,15 @@ module Numeric.LinearAlgebra.Algorithms ( | |||
72 | 74 | ||
73 | 75 | ||
74 | import Data.Packed.Internal hiding ((//)) | 76 | import Data.Packed.Internal hiding ((//)) |
75 | --import Data.Packed.Vector | ||
76 | import Data.Packed.Matrix | 77 | import Data.Packed.Matrix |
77 | import Data.Complex | 78 | import Data.Complex |
78 | import Numeric.GSL.Vector | ||
79 | import Numeric.LinearAlgebra.LAPACK as LAPACK | 79 | import Numeric.LinearAlgebra.LAPACK as LAPACK |
80 | import Numeric.LinearAlgebra.Linear | 80 | import Numeric.LinearAlgebra.Linear |
81 | import Data.List(foldl1') | 81 | import Data.List(foldl1') |
82 | import Data.Array | 82 | import Data.Array |
83 | 83 | ||
84 | -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. | 84 | -- | Auxiliary typeclass used to define generic computations for both real and complex matrices. |
85 | class (AutoReal t, Prod t, Linear Vector t, Linear Matrix t) => Field t where | 85 | class (Product t, Linear Vector t, Linear Matrix t) => Field t where |
86 | svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) | 86 | svd' :: Matrix t -> (Matrix t, Vector Double, Matrix t) |
87 | thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) | 87 | thinSVD' :: Matrix t -> (Matrix t, Vector Double, Matrix t) |
88 | sv' :: Matrix t -> Vector Double | 88 | sv' :: Matrix t -> Vector Double |
@@ -172,6 +172,9 @@ thinSVD = {-# SCC "thinSVD" #-} thinSVD' | |||
172 | singularValues :: Field t => Matrix t -> Vector Double | 172 | singularValues :: Field t => Matrix t -> Vector Double |
173 | singularValues = {-# SCC "singularValues" #-} sv' | 173 | singularValues = {-# SCC "singularValues" #-} sv' |
174 | 174 | ||
175 | instance (Field t, RealOf t ~ Double) => Norm2 Matrix t where | ||
176 | norm2 m = singularValues m @> 0 | ||
177 | |||
175 | -- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. | 178 | -- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values. |
176 | -- | 179 | -- |
177 | -- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@. | 180 | -- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@. |
@@ -379,80 +382,16 @@ ranksv teps maxdim s = k where | |||
379 | eps :: Double | 382 | eps :: Double |
380 | eps = 2.22044604925031e-16 | 383 | eps = 2.22044604925031e-16 |
381 | 384 | ||
382 | peps :: RealFloat x => x -> x | 385 | |
383 | peps x = 2.0**(fromIntegral $ 1-floatDigits x) | 386 | -- | 1 + 0.5*peps == 1, 1 + 0.6*peps /= 1 |
387 | peps :: RealFloat x => x | ||
388 | peps = x where x = 2.0**(fromIntegral $ 1-floatDigits x) | ||
384 | 389 | ||
385 | 390 | ||
386 | -- | The imaginary unit: @i = 0.0 :+ 1.0@ | 391 | -- | The imaginary unit: @i = 0.0 :+ 1.0@ |
387 | i :: Complex Double | 392 | i :: Complex Double |
388 | i = 0:+1 | 393 | i = 0:+1 |
389 | 394 | ||
390 | --------------------------------------------------------------------------- | ||
391 | |||
392 | norm2 :: Vector Double -> Double | ||
393 | norm2 = toScalarR Norm2 | ||
394 | |||
395 | norm1 :: Vector Double -> Double | ||
396 | norm1 = toScalarR AbsSum | ||
397 | |||
398 | data NormType = Infinity | PNorm1 | PNorm2 -- PNorm Int | ||
399 | |||
400 | pnormRV PNorm2 = norm2 | ||
401 | pnormRV PNorm1 = norm1 | ||
402 | pnormRV Infinity = vectorMax . vectorMapR Abs | ||
403 | --pnormRV _ = error "pnormRV not yet defined" | ||
404 | |||
405 | pnormCV PNorm2 = norm2 . asReal | ||
406 | pnormCV PNorm1 = norm1 . mapVector magnitude | ||
407 | pnormCV Infinity = vectorMax . mapVector magnitude | ||
408 | --pnormCV _ = error "pnormCV not yet defined" | ||
409 | |||
410 | pnormRM PNorm2 m = singularValues m @> 0 | ||
411 | pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m | ||
412 | pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m) | ||
413 | --pnormRM _ _ = error "p norm not yet defined" | ||
414 | |||
415 | pnormCM PNorm2 m = singularValues m @> 0 | ||
416 | pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m | ||
417 | pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m) | ||
418 | --pnormCM _ _ = error "p norm not yet defined" | ||
419 | |||
420 | -- | Objects which have a p-norm. | ||
421 | -- Using it you can define convenient shortcuts: | ||
422 | -- | ||
423 | -- @norm2 x = pnorm PNorm2 x@ | ||
424 | -- | ||
425 | -- @frobenius m = norm2 . flatten $ m@ | ||
426 | class Normed t where | ||
427 | pnorm :: NormType -> t -> Double | ||
428 | |||
429 | instance Normed (Vector Double) where | ||
430 | pnorm = pnormRV | ||
431 | |||
432 | instance Normed (Vector (Complex Double)) where | ||
433 | pnorm = pnormCV | ||
434 | |||
435 | instance Normed (Matrix Double) where | ||
436 | pnorm = pnormRM | ||
437 | |||
438 | instance Normed (Matrix (Complex Double)) where | ||
439 | pnorm = pnormCM | ||
440 | |||
441 | ----------------------------------------------------------------------- | ||
442 | -- to be optimized | ||
443 | |||
444 | instance Normed (Vector Float) where | ||
445 | pnorm t = pnorm t . double | ||
446 | |||
447 | instance Normed (Vector (Complex Float)) where | ||
448 | pnorm t = pnorm t . double | ||
449 | |||
450 | instance Normed (Matrix Float) where | ||
451 | pnorm t = pnorm t . double | ||
452 | |||
453 | instance Normed (Matrix (Complex Float)) where | ||
454 | pnorm t = pnorm t . double | ||
455 | |||
456 | ----------------------------------------------------------------------- | 395 | ----------------------------------------------------------------------- |
457 | 396 | ||
458 | -- | The nullspace of a matrix from its SVD decomposition. | 397 | -- | The nullspace of a matrix from its SVD decomposition. |
@@ -588,8 +527,8 @@ diagonalize m = if rank v == n | |||
588 | -- | 527 | -- |
589 | -- @logm = matFunc log@ | 528 | -- @logm = matFunc log@ |
590 | -- | 529 | -- |
591 | matFunc :: (Field t) => (Complex Double -> Complex Double) -> Matrix t -> Matrix (Complex Double) | 530 | matFunc :: (Complex Double -> Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double) |
592 | matFunc f m = case diagonalize (complex'' m) of | 531 | matFunc f m = case diagonalize m of |
593 | Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v | 532 | Just (l,v) -> v `mXm` diag (mapVector f l) `mXm` inv v |
594 | Nothing -> error "Sorry, matFunc requires a diagonalizable matrix" | 533 | Nothing -> error "Sorry, matFunc requires a diagonalizable matrix" |
595 | 534 | ||
@@ -607,7 +546,7 @@ epslist = [ (fromIntegral k, golubeps k k) | k <- [1..]] | |||
607 | geps delta = head [ k | (k,g) <- epslist, g<delta] | 546 | geps delta = head [ k | (k,g) <- epslist, g<delta] |
608 | 547 | ||
609 | expGolub m = iterate msq f !! j | 548 | expGolub m = iterate msq f !! j |
610 | where j = max 0 $ floor $ log2 $ pnorm Infinity m | 549 | where j = max 0 $ floor $ log2 $ normInf m |
611 | log2 x = log x / log 2 | 550 | log2 x = log x / log 2 |
612 | a = m */ fromIntegral ((2::Int)^j) | 551 | a = m */ fromIntegral ((2::Int)^j) |
613 | q = geps eps -- 7 steps | 552 | q = geps eps -- 7 steps |
@@ -630,7 +569,7 @@ expGolub m = iterate msq f !! j | |||
630 | {- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan, | 569 | {- | Matrix exponential. It uses a direct translation of Algorithm 11.3.1 in Golub & Van Loan, |
631 | based on a scaled Pade approximation. | 570 | based on a scaled Pade approximation. |
632 | -} | 571 | -} |
633 | expm :: (Normed (Matrix t), Field t) => Matrix t -> Matrix t | 572 | expm :: (Norm Vector t, Field t) => Matrix t -> Matrix t |
634 | expm = expGolub | 573 | expm = expGolub |
635 | 574 | ||
636 | -------------------------------------------------------------- | 575 | -------------------------------------------------------------- |
@@ -646,11 +585,11 @@ It only works with invertible matrices that have a real solution. For diagonaliz | |||
646 | [ 2.0, 2.25 | 585 | [ 2.0, 2.25 |
647 | , 0.0, 2.0 ]@ | 586 | , 0.0, 2.0 ]@ |
648 | -} | 587 | -} |
649 | sqrtm :: (Normed (Matrix t), Field t) => Matrix t -> Matrix t | 588 | sqrtm :: (Norm Vector t, Field t) => Matrix t -> Matrix t |
650 | sqrtm = sqrtmInv | 589 | sqrtm = sqrtmInv |
651 | 590 | ||
652 | sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x)) | 591 | sqrtmInv x = fst $ fixedPoint $ iterate f (x, ident (rows x)) |
653 | where fixedPoint (a:b:rest) | pnorm PNorm1 (fst a |-| fst b) < eps = a | 592 | where fixedPoint (a:b:rest) | norm1 (fst a |-| fst b) < peps = a |
654 | | otherwise = fixedPoint (b:rest) | 593 | | otherwise = fixedPoint (b:rest) |
655 | fixedPoint _ = error "fixedpoint with impossible inputs" | 594 | fixedPoint _ = error "fixedpoint with impossible inputs" |
656 | f (y,z) = (0.5 .* (y |+| inv z), | 595 | f (y,z) = (0.5 .* (y |+| inv z), |
@@ -691,4 +630,35 @@ luFact (l_u,perm) | r <= c = (l ,u ,p, s) | |||
691 | (|+|) = add | 630 | (|+|) = add |
692 | (|*|) = mul | 631 | (|*|) = mul |
693 | 632 | ||
694 | -------------------------------------------------- | 633 | --------------------------------------------------------------------------- |
634 | |||
635 | data NormType = Infinity | PNorm1 | PNorm2 | ||
636 | |||
637 | -- | Old class | ||
638 | class Normed t where | ||
639 | pnorm :: NormType -> t -> Double | ||
640 | |||
641 | instance Norm Vector t => Normed (Vector t) where | ||
642 | pnorm PNorm1 = realToFrac . norm1 | ||
643 | pnorm PNorm2 = realToFrac . normFrob | ||
644 | pnorm Infinity = realToFrac . normInf | ||
645 | |||
646 | instance Normed (Matrix Double) where | ||
647 | pnorm PNorm1 = norm1 | ||
648 | pnorm PNorm2 = norm2 | ||
649 | pnorm Infinity = normInf | ||
650 | |||
651 | instance Normed (Matrix (Complex Double)) where | ||
652 | pnorm PNorm1 = norm1 | ||
653 | pnorm PNorm2 = norm2 | ||
654 | pnorm Infinity = normInf | ||
655 | |||
656 | instance Normed (Matrix Float) where | ||
657 | pnorm PNorm1 = realToFrac . norm1 | ||
658 | pnorm PNorm2 = norm2 . double -- not yet optimized for Float | ||
659 | pnorm Infinity = realToFrac . normInf | ||
660 | |||
661 | instance Normed (Matrix (Complex Float)) where | ||
662 | pnorm PNorm1 = realToFrac . norm1 | ||
663 | pnorm PNorm2 = norm2 . double -- not yet optimized for Float | ||
664 | pnorm Infinity = realToFrac . normInf | ||
diff --git a/lib/Numeric/LinearAlgebra/Interface.hs b/lib/Numeric/LinearAlgebra/Interface.hs index 542d76e..ec08694 100644 --- a/lib/Numeric/LinearAlgebra/Interface.hs +++ b/lib/Numeric/LinearAlgebra/Interface.hs | |||
@@ -1,4 +1,5 @@ | |||
1 | {-# OPTIONS_GHC -fglasgow-exts #-} | 1 | {-# OPTIONS_GHC -fglasgow-exts #-} |
2 | {-# LANGUAGE UndecidableInstances #-} | ||
2 | ----------------------------------------------------------------------------- | 3 | ----------------------------------------------------------------------------- |
3 | {- | | 4 | {- | |
4 | Module : Numeric.LinearAlgebra.Interface | 5 | Module : Numeric.LinearAlgebra.Interface |
@@ -18,7 +19,7 @@ In the context of the standard numeric operators, one-component vectors and matr | |||
18 | ----------------------------------------------------------------------------- | 19 | ----------------------------------------------------------------------------- |
19 | 20 | ||
20 | module Numeric.LinearAlgebra.Interface( | 21 | module Numeric.LinearAlgebra.Interface( |
21 | (<>),(<.>), | 22 | (<>),(<.>),mulG, Adapt, adaptElements, |
22 | (<\>), | 23 | (<\>), |
23 | (.*),(*/), | 24 | (.*),(*/), |
24 | (<|>),(<->), | 25 | (<|>),(<->), |
@@ -28,22 +29,28 @@ import Data.Packed.Vector | |||
28 | import Data.Packed.Matrix | 29 | import Data.Packed.Matrix |
29 | import Numeric.LinearAlgebra.Algorithms | 30 | import Numeric.LinearAlgebra.Algorithms |
30 | import Numeric.LinearAlgebra.Linear | 31 | import Numeric.LinearAlgebra.Linear |
32 | import Data.Complex | ||
33 | import Control.Arrow((***)) | ||
31 | 34 | ||
32 | --import Numeric.GSL.Vector | 35 | --import Numeric.GSL.Vector |
33 | 36 | ||
34 | class Mul a b c | a b -> c where | 37 | class Mul a b c | a b -> c where |
35 | infixl 7 <> | 38 | infixl 7 <> |
36 | -- | Matrix-matrix, matrix-vector, and vector-matrix products. | 39 | -- | Matrix-matrix, matrix-vector, and vector-matrix products. |
37 | (<>) :: Prod t => a t -> b t -> c t | 40 | (<>) :: Product t => a t -> b t -> c t |
41 | mulG :: (Element r, Element s, Adapt r s t t, Product t) => a r -> b s -> c t | ||
38 | 42 | ||
39 | instance Mul Matrix Matrix Matrix where | 43 | instance Mul Matrix Matrix Matrix where |
40 | (<>) = multiply | 44 | (<>) = mXm |
45 | mulG a b = uncurry mXm (curry adapt a b) | ||
41 | 46 | ||
42 | instance Mul Matrix Vector Vector where | 47 | instance Mul Matrix Vector Vector where |
43 | (<>) m v = flatten $ m <> (asColumn v) | 48 | (<>) m v = flatten $ m <> (asColumn v) |
49 | mulG m v = flatten $ m `mulG` (asColumn v) | ||
44 | 50 | ||
45 | instance Mul Vector Matrix Vector where | 51 | instance Mul Vector Matrix Vector where |
46 | (<>) v m = flatten $ (asRow v) <> m | 52 | (<>) v m = flatten $ (asRow v) <> m |
53 | mulG v m = flatten $ (asRow v) `mulG` m | ||
47 | 54 | ||
48 | --------------------------------------------------- | 55 | --------------------------------------------------- |
49 | 56 | ||
@@ -120,3 +127,45 @@ a <-> b = joinV a b | |||
120 | 127 | ||
121 | ---------------------------------------------------- | 128 | ---------------------------------------------------- |
122 | 129 | ||
130 | class Adapt a b c d | a b -> c, a b -> d where | ||
131 | adapt :: Container k => (k a, k b) -> (k c, k d) | ||
132 | |||
133 | --instance Adapt a a a a where | ||
134 | -- adapt = id *** id | ||
135 | |||
136 | instance Adapt Float Float Float Float where | ||
137 | adapt = id *** id | ||
138 | |||
139 | instance Adapt Double Double Double Double where | ||
140 | adapt = id *** id | ||
141 | |||
142 | instance Adapt (Complex Float) (Complex Float) (Complex Float) (Complex Float) where | ||
143 | adapt = id *** id | ||
144 | |||
145 | instance Adapt Float Double Double Double where | ||
146 | adapt = double *** id | ||
147 | |||
148 | instance Adapt Double Float Double Double where | ||
149 | adapt = id *** double | ||
150 | |||
151 | instance Adapt Float (Complex Float) (Complex Float) (Complex Float) where | ||
152 | adapt = complex *** id | ||
153 | |||
154 | instance Adapt (Complex Float) Float (Complex Float) (Complex Float) where | ||
155 | adapt = id *** complex | ||
156 | |||
157 | instance (Convert a, Convert (DoubleOf a), ComplexOf (DoubleOf a) ~ Complex Double) => Adapt a (Complex Double) (Complex Double) (Complex Double) where | ||
158 | adapt = complex.double *** id | ||
159 | |||
160 | instance (Convert a, Convert (DoubleOf a), ComplexOf (DoubleOf a) ~ Complex Double) => Adapt (Complex Double) a (Complex Double) (Complex Double) where | ||
161 | adapt = id *** complex.double | ||
162 | |||
163 | instance Adapt Double (Complex Float) (Complex Double) (Complex Double) where | ||
164 | adapt = complex *** double | ||
165 | |||
166 | instance Adapt (Complex Float) Double (Complex Double) (Complex Double) where | ||
167 | adapt = double *** complex | ||
168 | |||
169 | adaptElements:: (Adapt a b c d, Container k) => (k a, k b) -> (k c, k d) | ||
170 | adaptElements p = adapt p | ||
171 | |||
diff --git a/lib/Numeric/LinearAlgebra/Linear.hs b/lib/Numeric/LinearAlgebra/Linear.hs index 67921d8..6c21a16 100644 --- a/lib/Numeric/LinearAlgebra/Linear.hs +++ b/lib/Numeric/LinearAlgebra/Linear.hs | |||
@@ -1,5 +1,6 @@ | |||
1 | {-# LANGUAGE UndecidableInstances, MultiParamTypeClasses, FlexibleInstances #-} | 1 | {-# LANGUAGE UndecidableInstances, MultiParamTypeClasses, FlexibleInstances #-} |
2 | {-# LANGUAGE FlexibleContexts #-} | 2 | {-# LANGUAGE FlexibleContexts #-} |
3 | {-# LANGUAGE TypeFamilies #-} | ||
3 | ----------------------------------------------------------------------------- | 4 | ----------------------------------------------------------------------------- |
4 | {- | | 5 | {- | |
5 | Module : Numeric.LinearAlgebra.Linear | 6 | Module : Numeric.LinearAlgebra.Linear |
@@ -20,9 +21,11 @@ module Numeric.LinearAlgebra.Linear ( | |||
20 | Vectors(..), | 21 | Vectors(..), |
21 | Linear(..), | 22 | Linear(..), |
22 | -- * Products | 23 | -- * Products |
23 | Prod(..), | 24 | Product(..), |
24 | mXm,mXv,vXm, | 25 | mXm,mXv,vXm, |
25 | outer, kronecker, | 26 | outer, kronecker, |
27 | -- * Norms | ||
28 | Norm(..), Norm2(..), | ||
26 | -- * Creation of numeric vectors | 29 | -- * Creation of numeric vectors |
27 | constant, linspace | 30 | constant, linspace |
28 | ) where | 31 | ) where |
@@ -190,38 +193,38 @@ linspace n (a,b) = addConstant a $ scale s $ fromList [0 .. fromIntegral n-1] | |||
190 | 193 | ||
191 | ---------------------------------------------------- | 194 | ---------------------------------------------------- |
192 | 195 | ||
193 | class Element t => Prod t where | 196 | class Element t => Product t where |
194 | multiply :: Matrix t -> Matrix t -> Matrix t | 197 | multiply :: Matrix t -> Matrix t -> Matrix t |
195 | ctrans :: Matrix t -> Matrix t | 198 | ctrans :: Matrix t -> Matrix t |
196 | 199 | ||
197 | instance Prod Double where | 200 | instance Product Double where |
198 | multiply = multiplyR | 201 | multiply = multiplyR |
199 | ctrans = trans | 202 | ctrans = trans |
200 | 203 | ||
201 | instance Prod (Complex Double) where | 204 | instance Product (Complex Double) where |
202 | multiply = multiplyC | 205 | multiply = multiplyC |
203 | ctrans = conj . trans | 206 | ctrans = conj . trans |
204 | 207 | ||
205 | instance Prod Float where | 208 | instance Product Float where |
206 | multiply = multiplyF | 209 | multiply = multiplyF |
207 | ctrans = trans | 210 | ctrans = trans |
208 | 211 | ||
209 | instance Prod (Complex Float) where | 212 | instance Product (Complex Float) where |
210 | multiply = multiplyQ | 213 | multiply = multiplyQ |
211 | ctrans = conj . trans | 214 | ctrans = conj . trans |
212 | 215 | ||
213 | ---------------------------------------------------------- | 216 | ---------------------------------------------------------- |
214 | 217 | ||
215 | -- synonym for matrix product | 218 | -- synonym for matrix product |
216 | mXm :: Prod t => Matrix t -> Matrix t -> Matrix t | 219 | mXm :: Product t => Matrix t -> Matrix t -> Matrix t |
217 | mXm = multiply | 220 | mXm = multiply |
218 | 221 | ||
219 | -- matrix - vector product | 222 | -- matrix - vector product |
220 | mXv :: Prod t => Matrix t -> Vector t -> Vector t | 223 | mXv :: Product t => Matrix t -> Vector t -> Vector t |
221 | mXv m v = flatten $ m `mXm` (asColumn v) | 224 | mXv m v = flatten $ m `mXm` (asColumn v) |
222 | 225 | ||
223 | -- vector - matrix product | 226 | -- vector - matrix product |
224 | vXm :: Prod t => Vector t -> Matrix t -> Vector t | 227 | vXm :: Product t => Vector t -> Matrix t -> Vector t |
225 | vXm v m = flatten $ (asRow v) `mXm` m | 228 | vXm v m = flatten $ (asRow v) `mXm` m |
226 | 229 | ||
227 | {- | Outer product of two vectors. | 230 | {- | Outer product of two vectors. |
@@ -232,7 +235,7 @@ vXm v m = flatten $ (asRow v) `mXm` m | |||
232 | , 10.0, 4.0, 6.0 | 235 | , 10.0, 4.0, 6.0 |
233 | , 15.0, 6.0, 9.0 ]@ | 236 | , 15.0, 6.0, 9.0 ]@ |
234 | -} | 237 | -} |
235 | outer :: (Prod t) => Vector t -> Vector t -> Matrix t | 238 | outer :: (Product t) => Vector t -> Vector t -> Matrix t |
236 | outer u v = asColumn u `multiply` asRow v | 239 | outer u v = asColumn u `multiply` asRow v |
237 | 240 | ||
238 | {- | Kronecker product of two matrices. | 241 | {- | Kronecker product of two matrices. |
@@ -257,9 +260,47 @@ m2=(4><3) | |||
257 | , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 | 260 | , 0.0, 0.0, 0.0, -7.0, -8.0, -9.0, 21.0, 24.0, 27.0 |
258 | , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ | 261 | , 0.0, 0.0, 0.0, -10.0, -11.0, -12.0, 30.0, 33.0, 36.0 ]@ |
259 | -} | 262 | -} |
260 | kronecker :: (Prod t) => Matrix t -> Matrix t -> Matrix t | 263 | kronecker :: (Product t) => Matrix t -> Matrix t -> Matrix t |
261 | kronecker a b = fromBlocks | 264 | kronecker a b = fromBlocks |
262 | . splitEvery (cols a) | 265 | . splitEvery (cols a) |
263 | . map (reshape (cols b)) | 266 | . map (reshape (cols b)) |
264 | . toRows | 267 | . toRows |
265 | $ flatten a `outer` flatten b | 268 | $ flatten a `outer` flatten b |
269 | |||
270 | -------------------------------------------------- | ||
271 | |||
272 | -- | simple norms | ||
273 | class (Element t, RealFloat (RealOf t)) => Norm c t where | ||
274 | norm1 :: c t -> RealOf t | ||
275 | normInf :: c t -> RealOf t | ||
276 | normFrob :: c t -> RealOf t | ||
277 | |||
278 | instance Norm Vector Double where | ||
279 | normFrob = toScalarR Norm2 | ||
280 | norm1 = toScalarR AbsSum | ||
281 | normInf = vectorMax . vectorMapR Abs | ||
282 | |||
283 | instance Norm Vector Float where | ||
284 | normFrob = toScalarF Norm2 | ||
285 | norm1 = toScalarF AbsSum | ||
286 | normInf = vectorMax . vectorMapF Abs | ||
287 | |||
288 | instance (Norm Vector t, Vectors Vector t, RealElement t | ||
289 | , RealOf t ~ t, RealOf (Complex t) ~ t | ||
290 | ) => Norm Vector (Complex t) where | ||
291 | normFrob = normFrob . asReal | ||
292 | norm1 = norm1 . mapVector magnitude | ||
293 | normInf = vectorMax . mapVector magnitude | ||
294 | |||
295 | instance Norm Vector t => Norm Matrix t where | ||
296 | normFrob = normFrob . flatten | ||
297 | norm1 = maximum . map norm1 . toColumns | ||
298 | normInf = norm1 . trans | ||
299 | |||
300 | class Norm2 c t where | ||
301 | norm2 :: c t -> RealOf t | ||
302 | |||
303 | instance Norm Vector t => Norm2 Vector t where | ||
304 | norm2 = normFrob | ||
305 | |||
306 | -- (the instance Norm2 Matrix t requires singular values and is defined later) | ||
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs index 91f6742..77bb2f7 100644 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ b/lib/Numeric/LinearAlgebra/Tests.hs | |||
@@ -168,7 +168,7 @@ fittingTest = utest "levmar" (ok1 && ok2) | |||
168 | sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] | 168 | sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] |
169 | 169 | ||
170 | ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d | 170 | ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d |
171 | ok2 = pnorm PNorm2 (fromList (map fst sols) - fromList sol) < 1E-5 | 171 | ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5 |
172 | 172 | ||
173 | ----------------------------------------------------- | 173 | ----------------------------------------------------- |
174 | 174 | ||
@@ -318,7 +318,7 @@ runTests n = do | |||
318 | test (cholProp . rPosDef) | 318 | test (cholProp . rPosDef) |
319 | test (cholProp . cPosDef) | 319 | test (cholProp . cPosDef) |
320 | putStrLn "------ expm" | 320 | putStrLn "------ expm" |
321 | test (expmDiagProp . rSqWC) | 321 | test (expmDiagProp . complex. rSqWC) |
322 | test (expmDiagProp . cSqWC) | 322 | test (expmDiagProp . cSqWC) |
323 | putStrLn "------ fft" | 323 | putStrLn "------ fft" |
324 | test (\v -> ifft (fft v) |~| v) | 324 | test (\v -> ifft (fft v) |~| v) |
diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs index aaaff28..21a6f88 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs | |||
@@ -27,13 +27,10 @@ module Numeric.LinearAlgebra.Tests.Instances( | |||
27 | ) where | 27 | ) where |
28 | 28 | ||
29 | 29 | ||
30 | import Numeric.LinearAlgebra hiding (real,complex) | 30 | import Numeric.LinearAlgebra |
31 | import Control.Monad(replicateM) | 31 | import Control.Monad(replicateM) |
32 | #include "quickCheckCompat.h" | 32 | #include "quickCheckCompat.h" |
33 | 33 | ||
34 | real x = real'' x | ||
35 | complex x = complex'' x | ||
36 | |||
37 | #if MIN_VERSION_QuickCheck(2,0,0) | 34 | #if MIN_VERSION_QuickCheck(2,0,0) |
38 | shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] | 35 | shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] |
39 | shrinkListElementwise [] = [] | 36 | shrinkListElementwise [] = [] |
@@ -72,7 +69,7 @@ instance (Field a, Arbitrary a) => Arbitrary (Vector a) where | |||
72 | #if MIN_VERSION_QuickCheck(2,0,0) | 69 | #if MIN_VERSION_QuickCheck(2,0,0) |
73 | -- shrink any one of the components | 70 | -- shrink any one of the components |
74 | shrink = map fromList . shrinkListElementwise . toList | 71 | shrink = map fromList . shrinkListElementwise . toList |
75 | 72 | ||
76 | #else | 73 | #else |
77 | coarbitrary = undefined | 74 | coarbitrary = undefined |
78 | #endif | 75 | #endif |
@@ -140,7 +137,7 @@ instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where | |||
140 | 137 | ||
141 | -- a well-conditioned general matrix (the singular values are between 1 and 100) | 138 | -- a well-conditioned general matrix (the singular values are between 1 and 100) |
142 | newtype (WC a) = WC (Matrix a) deriving Show | 139 | newtype (WC a) = WC (Matrix a) deriving Show |
143 | instance (Field a, Arbitrary a) => Arbitrary (WC a) where | 140 | instance (AutoReal a, Field a, Arbitrary a) => Arbitrary (WC a) where |
144 | arbitrary = do | 141 | arbitrary = do |
145 | m <- arbitrary | 142 | m <- arbitrary |
146 | let (u,_,v) = svd m | 143 | let (u,_,v) = svd m |
@@ -149,7 +146,7 @@ instance (Field a, Arbitrary a) => Arbitrary (WC a) where | |||
149 | n = min r c | 146 | n = min r c |
150 | sv' <- replicateM n (choose (1,100)) | 147 | sv' <- replicateM n (choose (1,100)) |
151 | let s = diagRect (fromList sv') r c | 148 | let s = diagRect (fromList sv') r c |
152 | return $ WC (u <> real s <> trans v) | 149 | return $ WC (u <> real'' s <> trans v) |
153 | 150 | ||
154 | #if MIN_VERSION_QuickCheck(2,0,0) | 151 | #if MIN_VERSION_QuickCheck(2,0,0) |
155 | #else | 152 | #else |
@@ -159,14 +156,14 @@ instance (Field a, Arbitrary a) => Arbitrary (WC a) where | |||
159 | 156 | ||
160 | -- a well-conditioned square matrix (the singular values are between 1 and 100) | 157 | -- a well-conditioned square matrix (the singular values are between 1 and 100) |
161 | newtype (SqWC a) = SqWC (Matrix a) deriving Show | 158 | newtype (SqWC a) = SqWC (Matrix a) deriving Show |
162 | instance (Field a, Arbitrary a) => Arbitrary (SqWC a) where | 159 | instance (AutoReal a, Field a, Arbitrary a) => Arbitrary (SqWC a) where |
163 | arbitrary = do | 160 | arbitrary = do |
164 | Sq m <- arbitrary | 161 | Sq m <- arbitrary |
165 | let (u,_,v) = svd m | 162 | let (u,_,v) = svd m |
166 | n = rows m | 163 | n = rows m |
167 | sv' <- replicateM n (choose (1,100)) | 164 | sv' <- replicateM n (choose (1,100)) |
168 | let s = diag (fromList sv') | 165 | let s = diag (fromList sv') |
169 | return $ SqWC (u <> real s <> trans v) | 166 | return $ SqWC (u <> real'' s <> trans v) |
170 | 167 | ||
171 | #if MIN_VERSION_QuickCheck(2,0,0) | 168 | #if MIN_VERSION_QuickCheck(2,0,0) |
172 | #else | 169 | #else |
@@ -176,14 +173,14 @@ instance (Field a, Arbitrary a) => Arbitrary (SqWC a) where | |||
176 | 173 | ||
177 | -- a positive definite square matrix (the eigenvalues are between 0 and 100) | 174 | -- a positive definite square matrix (the eigenvalues are between 0 and 100) |
178 | newtype (PosDef a) = PosDef (Matrix a) deriving Show | 175 | newtype (PosDef a) = PosDef (Matrix a) deriving Show |
179 | instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (PosDef a) where | 176 | instance (AutoReal a, Field a, Arbitrary a, Num (Vector a)) => Arbitrary (PosDef a) where |
180 | arbitrary = do | 177 | arbitrary = do |
181 | Her m <- arbitrary | 178 | Her m <- arbitrary |
182 | let (_,v) = eigSH m | 179 | let (_,v) = eigSH m |
183 | n = rows m | 180 | n = rows m |
184 | l <- replicateM n (choose (0,100)) | 181 | l <- replicateM n (choose (0,100)) |
185 | let s = diag (fromList l) | 182 | let s = diag (fromList l) |
186 | p = v <> real s <> ctrans v | 183 | p = v <> real'' s <> ctrans v |
187 | return $ PosDef (0.5 * p + 0.5 * ctrans p) | 184 | return $ PosDef (0.5 * p + 0.5 * ctrans p) |
188 | 185 | ||
189 | #if MIN_VERSION_QuickCheck(2,0,0) | 186 | #if MIN_VERSION_QuickCheck(2,0,0) |
@@ -243,3 +240,4 @@ cPosDef (PosDef m) = m :: CM | |||
243 | 240 | ||
244 | rConsist (Consistent (a,b)) = (a,b::RM) | 241 | rConsist (Consistent (a,b)) = (a,b::RM) |
245 | cConsist (Consistent (a,b)) = (a,b::CM) | 242 | cConsist (Consistent (a,b)) = (a,b::CM) |
243 | |||
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs index 9891d8a..d6bb338 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs | |||
@@ -54,15 +54,15 @@ complex x = complex'' x | |||
54 | debug x = trace (show x) x | 54 | debug x = trace (show x) x |
55 | 55 | ||
56 | -- relative error | 56 | -- relative error |
57 | dist :: (Normed t, Num t) => t -> t -> Double | 57 | --dist :: (Normed t, Num t) => t -> t -> Double |
58 | dist a b = r | 58 | dist a b = r |
59 | where norm = pnorm Infinity | 59 | where norm = normInf |
60 | na = norm a | 60 | na = norm a |
61 | nb = norm b | 61 | nb = norm b |
62 | nab = norm (a-b) | 62 | nab = norm (a-b) |
63 | mx = max na nb | 63 | mx = max na nb |
64 | mn = min na nb | 64 | mn = min na nb |
65 | r = if mn < eps | 65 | r = if mn < peps |
66 | then mx | 66 | then mx |
67 | else nab/mx | 67 | else nab/mx |
68 | 68 | ||
@@ -71,7 +71,7 @@ a |~| b = a :~10~: b | |||
71 | --a |~| b = dist a b < 10^^(-10) | 71 | --a |~| b = dist a b < 10^^(-10) |
72 | 72 | ||
73 | data Aprox a = (:~) a Int | 73 | data Aprox a = (:~) a Int |
74 | (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool | 74 | -- (~:) :: (Normed a, Num a) => Aprox a -> a -> Bool |
75 | a :~n~: b = dist a b < 10^^(-n) | 75 | a :~n~: b = dist a b < 10^^(-n) |
76 | 76 | ||
77 | ------------------------------------------------------ | 77 | ------------------------------------------------------ |