diff options
Diffstat (limited to 'examples')
-rw-r--r-- | examples/error.hs | 8 | ||||
-rw-r--r-- | examples/inplace.hs | 4 | ||||
-rw-r--r-- | examples/kalman.hs | 33 | ||||
-rw-r--r-- | examples/lie.hs | 10 | ||||
-rw-r--r-- | examples/minimize.hs | 9 | ||||
-rw-r--r-- | examples/monadic.hs | 31 | ||||
-rw-r--r-- | examples/multiply.hs | 22 | ||||
-rw-r--r-- | examples/ode.hs | 2 | ||||
-rw-r--r-- | examples/pca1.hs | 18 | ||||
-rw-r--r-- | examples/pca2.hs | 17 | ||||
-rw-r--r-- | examples/plot.hs | 4 | ||||
-rw-r--r-- | examples/root.hs | 6 | ||||
-rw-r--r-- | examples/vector.hs | 31 |
13 files changed, 87 insertions, 108 deletions
diff --git a/examples/error.hs b/examples/error.hs index 5efae7c..77467df 100644 --- a/examples/error.hs +++ b/examples/error.hs | |||
@@ -8,6 +8,7 @@ test x = catch | |||
8 | (print x) | 8 | (print x) |
9 | (\e -> putStrLn $ "captured ["++ show (e :: SomeException) ++"]") | 9 | (\e -> putStrLn $ "captured ["++ show (e :: SomeException) ++"]") |
10 | 10 | ||
11 | |||
11 | main = do | 12 | main = do |
12 | setErrorHandlerOff | 13 | setErrorHandlerOff |
13 | 14 | ||
@@ -15,7 +16,8 @@ main = do | |||
15 | test $ 5 + (fst.exp_e) 1000 | 16 | test $ 5 + (fst.exp_e) 1000 |
16 | test $ bessel_zero_Jnu_e (-0.3) 2 | 17 | test $ bessel_zero_Jnu_e (-0.3) 2 |
17 | 18 | ||
18 | test $ (linearSolve 0 4 :: Matrix Double) | 19 | test $ (inv 0 :: Matrix Double) |
19 | test $ (linearSolve 5 (sqrt (-1)) :: Matrix Double) | 20 | test $ (linearSolveLS 5 (sqrt (-1)) :: Matrix Double) |
21 | |||
22 | putStrLn "Bye" | ||
20 | 23 | ||
21 | putStrLn "Bye" \ No newline at end of file | ||
diff --git a/examples/inplace.hs b/examples/inplace.hs index 574aa44..19f9bc9 100644 --- a/examples/inplace.hs +++ b/examples/inplace.hs | |||
@@ -1,7 +1,9 @@ | |||
1 | -- some tests of the interface for pure | 1 | -- some tests of the interface for pure |
2 | -- computations with inplace updates | 2 | -- computations with inplace updates |
3 | 3 | ||
4 | import Numeric.LinearAlgebra.HMatrix | 4 | {-# LANGUAGE FlexibleContexts #-} |
5 | |||
6 | import Numeric.LinearAlgebra | ||
5 | import Numeric.LinearAlgebra.Devel | 7 | import Numeric.LinearAlgebra.Devel |
6 | 8 | ||
7 | import Data.Array.Unboxed | 9 | import Data.Array.Unboxed |
diff --git a/examples/kalman.hs b/examples/kalman.hs index 7fbe3d2..9756aa0 100644 --- a/examples/kalman.hs +++ b/examples/kalman.hs | |||
@@ -1,17 +1,15 @@ | |||
1 | import Numeric.LinearAlgebra | 1 | import Numeric.LinearAlgebra |
2 | import Graphics.Plot | 2 | import Graphics.Plot |
3 | 3 | ||
4 | vector l = fromList l :: Vector Double | 4 | f = fromLists |
5 | matrix ls = fromLists ls :: Matrix Double | 5 | [[1,0,0,0], |
6 | diagl = diag . vector | 6 | [1,1,0,0], |
7 | [0,0,1,0], | ||
8 | [0,0,0,1]] | ||
7 | 9 | ||
8 | f = matrix [[1,0,0,0], | 10 | h = fromLists |
9 | [1,1,0,0], | 11 | [[0,-1,1,0], |
10 | [0,0,1,0], | 12 | [0,-1,0,1]] |
11 | [0,0,0,1]] | ||
12 | |||
13 | h = matrix [[0,-1,1,0], | ||
14 | [0,-1,0,1]] | ||
15 | 13 | ||
16 | q = diagl [1,1,0,0] | 14 | q = diagl [1,1,0,0] |
17 | 15 | ||
@@ -25,13 +23,13 @@ type Measurement = Vector Double | |||
25 | 23 | ||
26 | kalman :: System -> State -> Measurement -> State | 24 | kalman :: System -> State -> Measurement -> State |
27 | kalman (System f h q r) (State x p) z = State x' p' where | 25 | kalman (System f h q r) (State x p) z = State x' p' where |
28 | px = f <> x -- prediction | 26 | px = f #> x -- prediction |
29 | pq = f <> p <> trans f + q -- its covariance | 27 | pq = f <> p <> tr f + q -- its covariance |
30 | y = z - h <> px -- residue | 28 | y = z - h #> px -- residue |
31 | cy = h <> pq <> trans h + r -- its covariance | 29 | cy = h <> pq <> tr h + r -- its covariance |
32 | k = pq <> trans h <> inv cy -- kalman gain | 30 | k = pq <> tr h <> inv cy -- kalman gain |
33 | x' = px + k <> y -- new state | 31 | x' = px + k #> y -- new state |
34 | p' = (ident (dim x) - k <> h) <> pq -- its covariance | 32 | p' = (ident (size x) - k <> h) <> pq -- its covariance |
35 | 33 | ||
36 | sys = System f h q r | 34 | sys = System f h q r |
37 | 35 | ||
@@ -49,3 +47,4 @@ main = do | |||
49 | print $ fromRows $ take 10 (map sX xs) | 47 | print $ fromRows $ take 10 (map sX xs) |
50 | mapM_ (print . sP) $ take 10 xs | 48 | mapM_ (print . sP) $ take 10 xs |
51 | mplot (evolution 20 (xs,des)) | 49 | mplot (evolution 20 (xs,des)) |
50 | |||
diff --git a/examples/lie.hs b/examples/lie.hs index db21ea8..4933df6 100644 --- a/examples/lie.hs +++ b/examples/lie.hs | |||
@@ -1,8 +1,8 @@ | |||
1 | -- The magic of Lie Algebra | 1 | -- The magic of Lie Algebra |
2 | 2 | ||
3 | import Numeric.LinearAlgebra | 3 | {-# LANGUAGE FlexibleContexts #-} |
4 | 4 | ||
5 | disp = putStrLn . dispf 5 | 5 | import Numeric.LinearAlgebra |
6 | 6 | ||
7 | rot1 :: Double -> Matrix Double | 7 | rot1 :: Double -> Matrix Double |
8 | rot1 a = (3><3) | 8 | rot1 a = (3><3) |
@@ -58,8 +58,8 @@ main = do | |||
58 | exact = rot3 a <> rot1 b <> rot2 c | 58 | exact = rot3 a <> rot1 b <> rot2 c |
59 | lie = scalar a * g3 |+| scalar b * g1 |+| scalar c * g2 | 59 | lie = scalar a * g3 |+| scalar b * g1 |+| scalar c * g2 |
60 | putStrLn "position in the tangent space:" | 60 | putStrLn "position in the tangent space:" |
61 | disp lie | 61 | disp 5 lie |
62 | putStrLn "exponential map back to the group (2 terms):" | 62 | putStrLn "exponential map back to the group (2 terms):" |
63 | disp (expm lie) | 63 | disp 5 (expm lie) |
64 | putStrLn "exact position:" | 64 | putStrLn "exact position:" |
65 | disp exact | 65 | disp 5 exact |
diff --git a/examples/minimize.hs b/examples/minimize.hs index 19b2cb3..c27afc2 100644 --- a/examples/minimize.hs +++ b/examples/minimize.hs | |||
@@ -20,7 +20,7 @@ partialDerivative n f v = fst (derivCentral 0.01 g (v!!n)) where | |||
20 | g x = f (concat [a,x:b]) | 20 | g x = f (concat [a,x:b]) |
21 | (a,_:b) = splitAt n v | 21 | (a,_:b) = splitAt n v |
22 | 22 | ||
23 | disp = putStrLn . format " " (printf "%.3f") | 23 | disp' = putStrLn . format " " (printf "%.3f") |
24 | 24 | ||
25 | allMethods :: (Enum a, Bounded a) => [a] | 25 | allMethods :: (Enum a, Bounded a) => [a] |
26 | allMethods = [minBound .. maxBound] | 26 | allMethods = [minBound .. maxBound] |
@@ -29,22 +29,23 @@ test method = do | |||
29 | print method | 29 | print method |
30 | let (s,p) = minimize method 1E-2 30 [1,1] f [5,7] | 30 | let (s,p) = minimize method 1E-2 30 [1,1] f [5,7] |
31 | print s | 31 | print s |
32 | disp p | 32 | disp' p |
33 | 33 | ||
34 | testD method = do | 34 | testD method = do |
35 | print method | 35 | print method |
36 | let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f df [5,7] | 36 | let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f df [5,7] |
37 | print s | 37 | print s |
38 | disp p | 38 | disp' p |
39 | 39 | ||
40 | testD' method = do | 40 | testD' method = do |
41 | putStrLn $ show method ++ " with estimated gradient" | 41 | putStrLn $ show method ++ " with estimated gradient" |
42 | let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f (gradient f) [5,7] | 42 | let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f (gradient f) [5,7] |
43 | print s | 43 | print s |
44 | disp p | 44 | disp' p |
45 | 45 | ||
46 | main = do | 46 | main = do |
47 | mapM_ test [NMSimplex, NMSimplex2] | 47 | mapM_ test [NMSimplex, NMSimplex2] |
48 | mapM_ testD allMethods | 48 | mapM_ testD allMethods |
49 | testD' ConjugateFR | 49 | testD' ConjugateFR |
50 | mplot $ drop 3 . toColumns . snd $ minimizeS f [5,7] | 50 | mplot $ drop 3 . toColumns . snd $ minimizeS f [5,7] |
51 | |||
diff --git a/examples/monadic.hs b/examples/monadic.hs index 7c6e0f4..cf8aacc 100644 --- a/examples/monadic.hs +++ b/examples/monadic.hs | |||
@@ -1,35 +1,37 @@ | |||
1 | -- monadic computations | 1 | -- monadic computations |
2 | -- (contributed by Vivian McPhail) | 2 | -- (contributed by Vivian McPhail) |
3 | 3 | ||
4 | {-# LANGUAGE FlexibleContexts #-} | ||
5 | |||
4 | import Numeric.LinearAlgebra | 6 | import Numeric.LinearAlgebra |
7 | import Numeric.LinearAlgebra.Devel | ||
5 | import Control.Monad.State.Strict | 8 | import Control.Monad.State.Strict |
6 | import Control.Monad.Maybe | 9 | import Control.Monad.Trans.Maybe |
7 | import Foreign.Storable(Storable) | 10 | import Foreign.Storable(Storable) |
8 | import System.Random(randomIO) | 11 | import System.Random(randomIO) |
9 | 12 | ||
10 | ------------------------------------------- | 13 | ------------------------------------------- |
11 | 14 | ||
12 | -- an instance of MonadIO, a monad transformer | 15 | -- an instance of MonadIO, a monad transformer |
13 | type VectorMonadT = StateT Int IO | 16 | type VectorMonadT = StateT I IO |
14 | 17 | ||
15 | test1 :: Vector Int -> IO (Vector Int) | 18 | test1 :: Vector I -> IO (Vector I) |
16 | test1 = mapVectorM $ \x -> do | 19 | test1 = mapVectorM $ \x -> do |
17 | putStr $ (show x) ++ " " | 20 | putStr $ (show x) ++ " " |
18 | return (x + 1) | 21 | return (x + 1) |
19 | 22 | ||
20 | -- we can have an arbitrary monad AND do IO | 23 | -- we can have an arbitrary monad AND do IO |
21 | addInitialM :: Vector Int -> VectorMonadT () | 24 | addInitialM :: Vector I -> VectorMonadT () |
22 | addInitialM = mapVectorM_ $ \x -> do | 25 | addInitialM = mapVectorM_ $ \x -> do |
23 | i <- get | 26 | i <- get |
24 | liftIO $ putStr $ (show $ x + i) ++ " " | 27 | liftIO $ putStr $ (show $ x + i) ++ " " |
25 | put $ x + i | 28 | put $ x + i |
26 | 29 | ||
27 | -- sum the values of the even indiced elements | 30 | -- sum the values of the even indiced elements |
28 | sumEvens :: Vector Int -> Int | 31 | sumEvens :: Vector I -> I |
29 | sumEvens = foldVectorWithIndex (\x a b -> if x `mod` 2 == 0 then a + b else b) 0 | 32 | sumEvens = foldVectorWithIndex (\x a b -> if x `mod` 2 == 0 then a + b else b) 0 |
30 | 33 | ||
31 | -- sum and print running total of evens | 34 | -- sum and print running total of evens |
32 | sumEvensAndPrint :: Vector Int -> VectorMonadT () | ||
33 | sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do | 35 | sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do |
34 | when (i `mod` 2 == 0) $ do | 36 | when (i `mod` 2 == 0) $ do |
35 | v <- get | 37 | v <- get |
@@ -38,7 +40,7 @@ sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do | |||
38 | liftIO $ putStr $ (show v') ++ " " | 40 | liftIO $ putStr $ (show v') ++ " " |
39 | 41 | ||
40 | 42 | ||
41 | indexPlusSum :: Vector Int -> VectorMonadT () | 43 | --indexPlusSum :: Vector I -> VectorMonadT () |
42 | indexPlusSum v' = do | 44 | indexPlusSum v' = do |
43 | let f i x = do | 45 | let f i x = do |
44 | s <- get | 46 | s <- get |
@@ -63,7 +65,7 @@ monoStep d = do | |||
63 | 65 | ||
64 | isMonotoneIncreasing :: Vector Double -> Bool | 66 | isMonotoneIncreasing :: Vector Double -> Bool |
65 | isMonotoneIncreasing v = | 67 | isMonotoneIncreasing v = |
66 | let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v @> 0) | 68 | let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v ! 0) |
67 | in case res of | 69 | in case res of |
68 | Nothing -> False | 70 | Nothing -> False |
69 | Just _ -> True | 71 | Just _ -> True |
@@ -72,8 +74,8 @@ isMonotoneIncreasing v = | |||
72 | ------------------------------------------- | 74 | ------------------------------------------- |
73 | 75 | ||
74 | -- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs | 76 | -- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs |
75 | successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool | 77 | successive_ :: (Container Vector a, Indexable (Vector a) a) => (a -> a -> Bool) -> Vector a -> Bool |
76 | successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0) | 78 | successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (size v - 1) v))) (v ! 0) |
77 | where step e = do | 79 | where step e = do |
78 | ep <- lift $ get | 80 | ep <- lift $ get |
79 | if t e ep | 81 | if t e ep |
@@ -81,8 +83,10 @@ successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ s | |||
81 | else (fail "successive_ test failed") | 83 | else (fail "successive_ test failed") |
82 | 84 | ||
83 | -- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input | 85 | -- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input |
84 | successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b | 86 | successive |
85 | successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0) | 87 | :: (Storable b, Container Vector s, Indexable (Vector s) s) |
88 | => (s -> s -> b) -> Vector s -> Vector b | ||
89 | successive f v = evalState (mapVectorM step (subVector 1 (size v - 1) v)) (v ! 0) | ||
86 | where step e = do | 90 | where step e = do |
87 | ep <- get | 91 | ep <- get |
88 | put e | 92 | put e |
@@ -90,7 +94,7 @@ successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0 | |||
90 | 94 | ||
91 | ------------------------------------------- | 95 | ------------------------------------------- |
92 | 96 | ||
93 | v :: Vector Int | 97 | v :: Vector I |
94 | v = 10 |> [0..] | 98 | v = 10 |> [0..] |
95 | 99 | ||
96 | w = fromList ([1..10]++[10,9..1]) :: Vector Double | 100 | w = fromList ([1..10]++[10,9..1]) :: Vector Double |
@@ -116,3 +120,4 @@ main = do | |||
116 | print $ successive_ (>) v | 120 | print $ successive_ (>) v |
117 | print $ successive_ (>) w | 121 | print $ successive_ (>) w |
118 | print $ successive (+) v | 122 | print $ successive (+) v |
123 | |||
diff --git a/examples/multiply.hs b/examples/multiply.hs index 572961c..be8fa73 100644 --- a/examples/multiply.hs +++ b/examples/multiply.hs | |||
@@ -22,10 +22,10 @@ instance Container Vector t => Scaling t (Vector t) (Vector t) where | |||
22 | instance Container Vector t => Scaling (Vector t) t (Vector t) where | 22 | instance Container Vector t => Scaling (Vector t) t (Vector t) where |
23 | (⋅) = flip scale | 23 | (⋅) = flip scale |
24 | 24 | ||
25 | instance Container Vector t => Scaling t (Matrix t) (Matrix t) where | 25 | instance (Num t, Container Vector t) => Scaling t (Matrix t) (Matrix t) where |
26 | (⋅) = scale | 26 | (⋅) = scale |
27 | 27 | ||
28 | instance Container Vector t => Scaling (Matrix t) t (Matrix t) where | 28 | instance (Num t, Container Vector t) => Scaling (Matrix t) t (Matrix t) where |
29 | (⋅) = flip scale | 29 | (⋅) = flip scale |
30 | 30 | ||
31 | 31 | ||
@@ -42,14 +42,14 @@ class Mul a b c | a b -> c, a c -> b, b c -> a where | |||
42 | instance Product t => Mul (Vector t) (Vector t) t where | 42 | instance Product t => Mul (Vector t) (Vector t) t where |
43 | (×) = udot | 43 | (×) = udot |
44 | 44 | ||
45 | instance Product t => Mul (Matrix t) (Vector t) (Vector t) where | 45 | instance (Numeric t, Product t) => Mul (Matrix t) (Vector t) (Vector t) where |
46 | (×) = mXv | 46 | (×) = (#>) |
47 | 47 | ||
48 | instance Product t => Mul (Vector t) (Matrix t) (Vector t) where | 48 | instance (Numeric t, Product t) => Mul (Vector t) (Matrix t) (Vector t) where |
49 | (×) = vXm | 49 | (×) = (<#) |
50 | 50 | ||
51 | instance Product t => Mul (Matrix t) (Matrix t) (Matrix t) where | 51 | instance (Numeric t, Product t) => Mul (Matrix t) (Matrix t) (Matrix t) where |
52 | (×) = mXm | 52 | (×) = (<>) |
53 | 53 | ||
54 | 54 | ||
55 | --instance Scaling a b c => Contraction a b c where | 55 | --instance Scaling a b c => Contraction a b c where |
@@ -92,9 +92,9 @@ u = fromList [3,0,5] | |||
92 | w = konst 1 (2,3) :: Matrix Double | 92 | w = konst 1 (2,3) :: Matrix Double |
93 | 93 | ||
94 | main = do | 94 | main = do |
95 | print $ (scale s v <> m) `udot` v | 95 | print $ (scale s v <# m) `udot` v |
96 | print $ scale s v `udot` (m <> v) | 96 | print $ scale s v `udot` (m #> v) |
97 | print $ s * ((v <> m) `udot` v) | 97 | print $ s * ((v <# m) `udot` v) |
98 | print $ s ⋅ v × m × v | 98 | print $ s ⋅ v × m × v |
99 | print a | 99 | print a |
100 | -- print (b == c) | 100 | -- print (b == c) |
diff --git a/examples/ode.hs b/examples/ode.hs index dc6e0ec..4cf1673 100644 --- a/examples/ode.hs +++ b/examples/ode.hs | |||
@@ -43,7 +43,7 @@ vanderpol' mu = do | |||
43 | jac t (toList->[x,v]) = (2><2) [ 0 , 1 | 43 | jac t (toList->[x,v]) = (2><2) [ 0 , 1 |
44 | , -1-2*x*v*mu, mu*(1-x**2) ] | 44 | , -1-2*x*v*mu, mu*(1-x**2) ] |
45 | ts = linspace 1000 (0,50) | 45 | ts = linspace 1000 (0,50) |
46 | hi = (ts@>1 - ts@>0)/100 | 46 | hi = (ts!1 - ts!0)/100 |
47 | sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts | 47 | sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts |
48 | mplot sol | 48 | mplot sol |
49 | 49 | ||
diff --git a/examples/pca1.hs b/examples/pca1.hs index a11eba9..ad2214d 100644 --- a/examples/pca1.hs +++ b/examples/pca1.hs | |||
@@ -8,27 +8,25 @@ import Control.Monad(when) | |||
8 | type Vec = Vector Double | 8 | type Vec = Vector Double |
9 | type Mat = Matrix Double | 9 | type Mat = Matrix Double |
10 | 10 | ||
11 | 11 | {- | |
12 | -- Vector with the mean value of the columns of a matrix | 12 | -- Vector with the mean value of the columns of a matrix |
13 | mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a | 13 | mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a |
14 | 14 | ||
15 | -- covariance matrix of a list of observations stored as rows | 15 | -- covariance matrix of a list of observations stored as rows |
16 | cov x = (trans xc <> xc) / fromIntegral (rows x - 1) | 16 | cov x = (trans xc <> xc) / fromIntegral (rows x - 1) |
17 | where xc = x - asRow (mean x) | 17 | where xc = x - asRow (mean x) |
18 | -} | ||
18 | 19 | ||
19 | 20 | ||
20 | -- creates the compression and decompression functions from the desired number of components | 21 | -- creates the compression and decompression functions from the desired number of components |
21 | pca :: Int -> Mat -> (Vec -> Vec , Vec -> Vec) | 22 | pca :: Int -> Mat -> (Vec -> Vec , Vec -> Vec) |
22 | pca n dataSet = (encode,decode) | 23 | pca n dataSet = (encode,decode) |
23 | where | 24 | where |
24 | encode x = vp <> (x - m) | 25 | encode x = vp #> (x - m) |
25 | decode x = x <> vp + m | 26 | decode x = x <# vp + m |
26 | m = mean dataSet | 27 | (m,c) = meanCov dataSet |
27 | c = cov dataSet | 28 | (_,v) = eigSH (trustSym c) |
28 | (_,v) = eigSH' c | 29 | vp = tr $ takeColumns n v |
29 | vp = takeRows n (trans v) | ||
30 | |||
31 | norm = pnorm PNorm2 | ||
32 | 30 | ||
33 | main = do | 31 | main = do |
34 | ok <- doesFileExist ("mnist.txt") | 32 | ok <- doesFileExist ("mnist.txt") |
@@ -43,4 +41,4 @@ main = do | |||
43 | let (pe,pd) = pca 10 xs | 41 | let (pe,pd) = pca 10 xs |
44 | let y = pe x | 42 | let y = pe x |
45 | print y -- compressed version | 43 | print y -- compressed version |
46 | print $ norm (x - pd y) / norm x --reconstruction quality | 44 | print $ norm_2 (x - pd y) / norm_2 x --reconstruction quality |
diff --git a/examples/pca2.hs b/examples/pca2.hs index e7ea95f..892d382 100644 --- a/examples/pca2.hs +++ b/examples/pca2.hs | |||
@@ -1,5 +1,7 @@ | |||
1 | -- Improved PCA, including illustrative graphics | 1 | -- Improved PCA, including illustrative graphics |
2 | 2 | ||
3 | {-# LANGUAGE FlexibleContexts #-} | ||
4 | |||
3 | import Numeric.LinearAlgebra | 5 | import Numeric.LinearAlgebra |
4 | import Graphics.Plot | 6 | import Graphics.Plot |
5 | import System.Directory(doesFileExist) | 7 | import System.Directory(doesFileExist) |
@@ -10,27 +12,27 @@ type Vec = Vector Double | |||
10 | type Mat = Matrix Double | 12 | type Mat = Matrix Double |
11 | 13 | ||
12 | -- Vector with the mean value of the columns of a matrix | 14 | -- Vector with the mean value of the columns of a matrix |
13 | mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a | 15 | mean a = konst (recip . fromIntegral . rows $ a) (rows a) <# a |
14 | 16 | ||
15 | -- covariance matrix of a list of observations stored as rows | 17 | -- covariance matrix of a list of observations stored as rows |
16 | cov x = (trans xc <> xc) / fromIntegral (rows x - 1) | 18 | cov x = (mTm xc) -- / fromIntegral (rows x - 1) |
17 | where xc = x - asRow (mean x) | 19 | where xc = x - asRow (mean x) |
18 | 20 | ||
19 | 21 | ||
20 | type Stat = (Vec, [Double], Mat) | 22 | type Stat = (Vec, [Double], Mat) |
21 | -- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov) | 23 | -- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov) |
22 | stat :: Mat -> Stat | 24 | stat :: Mat -> Stat |
23 | stat x = (m, toList s, trans v) where | 25 | stat x = (m, toList s, tr v) where |
24 | m = mean x | 26 | m = mean x |
25 | (s,v) = eigSH' (cov x) | 27 | (s,v) = eigSH (cov x) |
26 | 28 | ||
27 | -- creates the compression and decompression functions from the desired reconstruction | 29 | -- creates the compression and decompression functions from the desired reconstruction |
28 | -- quality and the statistics of a data set | 30 | -- quality and the statistics of a data set |
29 | pca :: Double -> Stat -> (Vec -> Vec , Vec -> Vec) | 31 | pca :: Double -> Stat -> (Vec -> Vec , Vec -> Vec) |
30 | pca prec (m,s,v) = (encode,decode) | 32 | pca prec (m,s,v) = (encode,decode) |
31 | where | 33 | where |
32 | encode x = vp <> (x - m) | 34 | encode x = vp #> (x - m) |
33 | decode x = x <> vp + m | 35 | decode x = x <# vp + m |
34 | vp = takeRows n v | 36 | vp = takeRows n v |
35 | n = 1 + (length $ fst $ span (< (prec'*sum s)) $ cumSum s) | 37 | n = 1 + (length $ fst $ span (< (prec'*sum s)) $ cumSum s) |
36 | cumSum = tail . scanl (+) 0.0 | 38 | cumSum = tail . scanl (+) 0.0 |
@@ -46,7 +48,7 @@ test :: Stat -> Double -> Vec -> IO () | |||
46 | test st prec x = do | 48 | test st prec x = do |
47 | let (pe,pd) = pca prec st | 49 | let (pe,pd) = pca prec st |
48 | let y = pe x | 50 | let y = pe x |
49 | print $ dim y | 51 | print $ size y |
50 | shdigit (pd y) | 52 | shdigit (pd y) |
51 | 53 | ||
52 | main = do | 54 | main = do |
@@ -63,3 +65,4 @@ main = do | |||
63 | let st = stat xs | 65 | let st = stat xs |
64 | test st 0.90 x | 66 | test st 0.90 x |
65 | test st 0.50 x | 67 | test st 0.50 x |
68 | |||
diff --git a/examples/plot.hs b/examples/plot.hs index f950aa5..90643ed 100644 --- a/examples/plot.hs +++ b/examples/plot.hs | |||
@@ -16,5 +16,5 @@ cumdist x = 0.5 * (1+ erf (x/sqrt 2)) | |||
16 | main = do | 16 | main = do |
17 | let x = linspace 1000 (-4,4) | 17 | let x = linspace 1000 (-4,4) |
18 | mplot [f x] | 18 | mplot [f x] |
19 | mplot [x, mapVector cumdist x, mapVector gaussianPDF x] | 19 | mplot [x, cmap cumdist x, cmap gaussianPDF x] |
20 | mesh (sombrero 40) \ No newline at end of file | 20 | mesh (sombrero 40) |
diff --git a/examples/root.hs b/examples/root.hs index 8546ff5..fa6e77a 100644 --- a/examples/root.hs +++ b/examples/root.hs | |||
@@ -9,7 +9,7 @@ test method = do | |||
9 | print method | 9 | print method |
10 | let (s,p) = root method 1E-7 30 (rosenbrock 1 10) [-10,-5] | 10 | let (s,p) = root method 1E-7 30 (rosenbrock 1 10) [-10,-5] |
11 | print s -- solution | 11 | print s -- solution |
12 | disp p -- evolution of the algorithm | 12 | disp' p -- evolution of the algorithm |
13 | 13 | ||
14 | jacobian a b [x,y] = [ [-a , 0] | 14 | jacobian a b [x,y] = [ [-a , 0] |
15 | , [-2*b*x, b] ] | 15 | , [-2*b*x, b] ] |
@@ -18,9 +18,9 @@ testJ method = do | |||
18 | print method | 18 | print method |
19 | let (s,p) = rootJ method 1E-7 30 (rosenbrock 1 10) (jacobian 1 10) [-10,-5] | 19 | let (s,p) = rootJ method 1E-7 30 (rosenbrock 1 10) (jacobian 1 10) [-10,-5] |
20 | print s | 20 | print s |
21 | disp p | 21 | disp' p |
22 | 22 | ||
23 | disp = putStrLn . format " " (printf "%.3f") | 23 | disp' = putStrLn . format " " (printf "%.3f") |
24 | 24 | ||
25 | main = do | 25 | main = do |
26 | test Hybrids | 26 | test Hybrids |
diff --git a/examples/vector.hs b/examples/vector.hs deleted file mode 100644 index f531cbd..0000000 --- a/examples/vector.hs +++ /dev/null | |||
@@ -1,31 +0,0 @@ | |||
1 | -- conversion to/from Data.Vector.Storable | ||
2 | -- from Roman Leshchinskiy "vector" package | ||
3 | -- | ||
4 | -- In the future Data.Packed.Vector will be replaced by Data.Vector.Storable | ||
5 | |||
6 | ------------------------------------------- | ||
7 | |||
8 | import Numeric.LinearAlgebra as H | ||
9 | import Data.Packed.Development(unsafeFromForeignPtr, unsafeToForeignPtr) | ||
10 | import Foreign.Storable | ||
11 | import qualified Data.Vector.Storable as V | ||
12 | |||
13 | fromVector :: Storable t => V.Vector t -> H.Vector t | ||
14 | fromVector v = unsafeFromForeignPtr p i n where | ||
15 | (p,i,n) = V.unsafeToForeignPtr v | ||
16 | |||
17 | toVector :: Storable t => H.Vector t -> V.Vector t | ||
18 | toVector v = V.unsafeFromForeignPtr p i n where | ||
19 | (p,i,n) = unsafeToForeignPtr v | ||
20 | |||
21 | ------------------------------------------- | ||
22 | |||
23 | v = V.slice 5 10 (V.fromList [1 .. 10::Double] V.++ V.replicate 10 7) | ||
24 | |||
25 | w = subVector 2 3 (linspace 5 (0,1)) :: Vector Double | ||
26 | |||
27 | main = do | ||
28 | print v | ||
29 | print $ fromVector v | ||
30 | print w | ||
31 | print $ toVector w | ||