diff options
author | Alberto Ruiz <aruiz@um.es> | 2007-09-21 18:28:08 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2007-09-21 18:28:08 +0000 |
commit | 0198366bba7a5f2d67338633f9eb90889ffc31b2 (patch) | |
tree | 4897d90233b333ee2092e63a4b74c7bcb2d22577 | |
parent | d4cb2692f9dae748da23371057a983deca4b2f80 (diff) |
add examples
-rw-r--r-- | examples/data.txt | 10 | ||||
-rw-r--r-- | examples/deriv.hs | 8 | ||||
-rw-r--r-- | examples/error.hs | 16 | ||||
-rw-r--r-- | examples/integrate.hs | 17 | ||||
-rw-r--r-- | examples/kalman.hs | 56 | ||||
-rw-r--r-- | examples/listlike.hs | 32 | ||||
-rw-r--r-- | examples/minimize.hs | 43 | ||||
-rw-r--r-- | examples/pca1.hs | 45 | ||||
-rw-r--r-- | examples/pca2.hs | 67 | ||||
-rw-r--r-- | examples/pinv1.hs | 18 | ||||
-rw-r--r-- | examples/pinv2.hs | 26 | ||||
-rw-r--r-- | examples/plot.hs | 20 | ||||
-rw-r--r-- | examples/testmnist.m | 29 |
13 files changed, 387 insertions, 0 deletions
diff --git a/examples/data.txt b/examples/data.txt new file mode 100644 index 0000000..7031462 --- /dev/null +++ b/examples/data.txt | |||
@@ -0,0 +1,10 @@ | |||
1 | 1 1.1 | ||
2 | 2 3.9 | ||
3 | 3 9.2 | ||
4 | 4 15.8 | ||
5 | 5 25.3 | ||
6 | 6 35.7 | ||
7 | 7 49.4 | ||
8 | 8 63.6 | ||
9 | 9 81.5 | ||
10 | 10 99.5 \ No newline at end of file | ||
diff --git a/examples/deriv.hs b/examples/deriv.hs new file mode 100644 index 0000000..472a284 --- /dev/null +++ b/examples/deriv.hs | |||
@@ -0,0 +1,8 @@ | |||
1 | -- Numerical differentiation | ||
2 | |||
3 | import GSL | ||
4 | |||
5 | d :: (Double -> Double) -> (Double -> Double) | ||
6 | d f x = fst $ derivCentral 0.01 f x | ||
7 | |||
8 | main = print $ d (\x-> x * d (\y-> x+y) 1) 1 | ||
diff --git a/examples/error.hs b/examples/error.hs new file mode 100644 index 0000000..bcfe224 --- /dev/null +++ b/examples/error.hs | |||
@@ -0,0 +1,16 @@ | |||
1 | import GSL | ||
2 | import Prelude hiding (catch) | ||
3 | import Control.Exception | ||
4 | |||
5 | test x = catch | ||
6 | (print x) | ||
7 | (\e -> putStrLn $ "captured ["++ show e++"]") | ||
8 | |||
9 | main = do | ||
10 | setErrorHandlerOff | ||
11 | |||
12 | test $ log_e (-1) | ||
13 | test $ 5 + (fst.exp_e) 1000 | ||
14 | test $ bessel_zero_Jnu_e (-0.3) 2 | ||
15 | |||
16 | putStrLn "Bye" \ No newline at end of file | ||
diff --git a/examples/integrate.hs b/examples/integrate.hs new file mode 100644 index 0000000..6da88ad --- /dev/null +++ b/examples/integrate.hs | |||
@@ -0,0 +1,17 @@ | |||
1 | -- Numerical integration | ||
2 | import GSL | ||
3 | |||
4 | quad f a b = fst $ integrateQAGS 1E-9 100 f a b | ||
5 | |||
6 | -- A multiple integral can be easily defined using partial application | ||
7 | quad2 f a b g1 g2 = quad h a b | ||
8 | where h x = quad (f x) (g1 x) (g2 x) | ||
9 | |||
10 | volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) | ||
11 | 0 r (const 0) (\x->sqrt (r*r-x*x)) | ||
12 | |||
13 | main = do | ||
14 | print $ quad (\x -> 4/(x^2+1)) 0 1 | ||
15 | print pi | ||
16 | print $ volSphere 2.5 | ||
17 | print $ 4/3*pi*2.5**3 | ||
diff --git a/examples/kalman.hs b/examples/kalman.hs new file mode 100644 index 0000000..ccb8083 --- /dev/null +++ b/examples/kalman.hs | |||
@@ -0,0 +1,56 @@ | |||
1 | import LinearAlgebra | ||
2 | import Graphics.Plot | ||
3 | import LinearAlgebra.Instances | ||
4 | |||
5 | --import GSLHaskell | ||
6 | |||
7 | vector l = fromList l :: Vector Double | ||
8 | matrix ls = fromLists ls :: Matrix Double | ||
9 | diagl = diag . vector | ||
10 | |||
11 | f = matrix [[1,0,0,0], | ||
12 | [1,1,0,0], | ||
13 | [0,0,1,0], | ||
14 | [0,0,0,1]] | ||
15 | |||
16 | h = matrix [[0,-1,1,0], | ||
17 | [0,-1,0,1]] | ||
18 | |||
19 | q = diagl [1,1,0,0] | ||
20 | |||
21 | r = diagl [2,2] | ||
22 | |||
23 | s0 = State (vector [0, 0, 10, -10]) (diagl [10,0, 100, 100]) | ||
24 | |||
25 | data System = System {kF, kH, kQ, kR :: Matrix Double} | ||
26 | data State = State {sX :: Vector Double , sP :: Matrix Double} | ||
27 | type Measurement = Vector Double | ||
28 | |||
29 | kalman :: System -> State -> Measurement -> State | ||
30 | kalman (System f h q r) (State x p) z = State x' p' where | ||
31 | px = f <> x -- prediction | ||
32 | pq = f <> p <> trans f -- its covariance | ||
33 | y = z - h <> px -- residue | ||
34 | cy = h <> pq <> trans h + r -- its covariance | ||
35 | k = pq <> trans h <> inv cy -- kalman gain | ||
36 | x' = px + k <> y -- new state | ||
37 | p' = (ident (dim x) - k <> h) <> pq -- its covariance | ||
38 | |||
39 | sys = System f h q r | ||
40 | |||
41 | zs = [vector [15-k,-20-k] | k <- [0..]] | ||
42 | |||
43 | xs = s0 : zipWith (kalman sys) xs zs | ||
44 | |||
45 | des = map (sqrt.takeDiag.sP) xs | ||
46 | |||
47 | evolution n (xs,des) = | ||
48 | vector [1.. fromIntegral n]:(toColumns $ fromRows $ take n (zipWith (-) (map sX xs) des)) ++ | ||
49 | (toColumns $ fromRows $ take n (zipWith (+) (map sX xs) des)) | ||
50 | |||
51 | main = do | ||
52 | print $ fromRows $ take 10 (map sX xs) | ||
53 | mapM_ (print . sP) $ take 10 xs | ||
54 | mplot (evolution 20 (xs,des)) | ||
55 | |||
56 | --(<>) = multiply \ No newline at end of file | ||
diff --git a/examples/listlike.hs b/examples/listlike.hs new file mode 100644 index 0000000..6c54f17 --- /dev/null +++ b/examples/listlike.hs | |||
@@ -0,0 +1,32 @@ | |||
1 | {-# OPTIONS_GHC -fglasgow-exts #-} | ||
2 | |||
3 | import qualified Data.ListLike as LL | ||
4 | import LinearAlgebra | ||
5 | import Data.Monoid | ||
6 | import Data.Packed.Internal.Vector | ||
7 | import Foreign | ||
8 | |||
9 | instance (Storable a) => Monoid (Vector a) where | ||
10 | mempty = V { dim = 0, fptr = undefined, ptr = undefined } | ||
11 | mappend a b = mconcat [a,b] | ||
12 | mconcat = j . filter ((>0).dim) | ||
13 | where j [] = mempty | ||
14 | j l = join l | ||
15 | |||
16 | instance Storable a => LL.FoldableLL (Vector a) a where | ||
17 | foldl f x v = foldl f x (toList v) | ||
18 | foldr f x v = foldr f x (toList v) | ||
19 | |||
20 | instance Storable a => LL.ListLike (Vector a) a where | ||
21 | singleton a = fromList [a] | ||
22 | head a = a @> 0 | ||
23 | tail a | dim a == 1 = mempty | ||
24 | | otherwise = subVector 1 (dim a -1) a | ||
25 | genericLength = fromIntegral.dim | ||
26 | |||
27 | |||
28 | v k = fromList [1..k] :: Vector Double | ||
29 | |||
30 | f k = k+(3::Double) | ||
31 | |||
32 | main = print $ (LL.map f [1..5] :: Vector Double) \ No newline at end of file | ||
diff --git a/examples/minimize.hs b/examples/minimize.hs new file mode 100644 index 0000000..0429a24 --- /dev/null +++ b/examples/minimize.hs | |||
@@ -0,0 +1,43 @@ | |||
1 | -- the multidimensional minimization example in the GSL manual | ||
2 | import GSL | ||
3 | import LinearAlgebra | ||
4 | import Graphics.Plot | ||
5 | |||
6 | -- the function to be minimized | ||
7 | f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30 | ||
8 | |||
9 | -- its gradient | ||
10 | df [x,y] = [20*(x-1), 40*(y-2)] | ||
11 | |||
12 | -- the conjugate gradient method | ||
13 | minimizeCG = minimizeConjugateGradient 1E-2 1E-4 1E-3 30 | ||
14 | |||
15 | -- a minimization algorithm which does not require the gradient | ||
16 | minimizeS f xi = minimizeNMSimplex f xi (replicate (length xi) 1) 1E-2 100 | ||
17 | |||
18 | -- Numerical estimation of the gradient | ||
19 | gradient f v = [partialDerivative k f v | k <- [0 .. length v -1]] | ||
20 | |||
21 | partialDerivative n f v = fst (derivCentral 0.01 g (v!!n)) where | ||
22 | g x = f (concat [a,x:b]) | ||
23 | (a,_:b) = splitAt n v | ||
24 | |||
25 | main = do | ||
26 | -- conjugate gradient with true gradient | ||
27 | let (s,p) = minimizeCG f df [5,7] | ||
28 | print s -- solution | ||
29 | dispR 2 p -- evolution of the algorithm | ||
30 | let [x,y] = drop 2 (toColumns p) | ||
31 | mplot [x,y] -- path from the starting point to the solution | ||
32 | |||
33 | -- conjugate gradient with estimated gradient | ||
34 | let (s,p) = minimizeCG f (gradient f) [5,7] | ||
35 | print s | ||
36 | dispR 2 p | ||
37 | mplot $ drop 2 (toColumns p) | ||
38 | |||
39 | -- without gradient, using the NM Simplex method | ||
40 | let (s,p) = minimizeS f [5,7] | ||
41 | print s | ||
42 | dispR 2 p | ||
43 | mplot $ drop 3 (toColumns p) | ||
diff --git a/examples/pca1.hs b/examples/pca1.hs new file mode 100644 index 0000000..2c5074d --- /dev/null +++ b/examples/pca1.hs | |||
@@ -0,0 +1,45 @@ | |||
1 | -- Principal component analysis | ||
2 | |||
3 | import LinearAlgebra | ||
4 | import System.Directory(doesFileExist) | ||
5 | import System(system) | ||
6 | import Control.Monad(when) | ||
7 | |||
8 | type Vec = Vector Double | ||
9 | type Mat = Matrix Double | ||
10 | |||
11 | sumColumns m = constant 1 (rows m) <> m | ||
12 | |||
13 | -- Vec with the mean value of the columns of a Mat | ||
14 | mean x = sumColumns x / fromIntegral (rows x) | ||
15 | |||
16 | -- covariance Mat of a list of observations as rows of a Mat | ||
17 | cov x = (trans xc <> xc) / fromIntegral (rows x -1) | ||
18 | where xc = center x | ||
19 | center m = m - constant 1 (rows m) `outer` mean m | ||
20 | |||
21 | -- creates the compression and decompression functions from the desired number of components | ||
22 | pca :: Int -> Mat -> (Vec -> Vec , Vec -> Vec) | ||
23 | pca n dataSet = (encode,decode) | ||
24 | where | ||
25 | encode x = vp <> (x - m) | ||
26 | decode x = x <> vp + m | ||
27 | m = mean dataSet | ||
28 | c = cov dataSet | ||
29 | (_,v) = eigS c | ||
30 | vp = takeRows n (trans v) | ||
31 | |||
32 | main = do | ||
33 | ok <- doesFileExist ("mnist.txt") | ||
34 | when (not ok) $ do | ||
35 | putStrLn "\nTrying to download test datafile..." | ||
36 | system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz") | ||
37 | system("gunzip mnist.txt.gz") | ||
38 | return () | ||
39 | m <- fromFile "mnist.txt" (5000,785) | ||
40 | let xs = takeColumns (cols m -1) m -- the last column is the digit type (class label) | ||
41 | let x = toRows xs !! 4 -- an arbitrary test Vec | ||
42 | let (pe,pd) = pca 10 xs | ||
43 | let y = pe x | ||
44 | print y -- compressed version | ||
45 | print $ norm (x - pd y) / norm x --reconstruction quality | ||
diff --git a/examples/pca2.hs b/examples/pca2.hs new file mode 100644 index 0000000..bd498e4 --- /dev/null +++ b/examples/pca2.hs | |||
@@ -0,0 +1,67 @@ | |||
1 | -- Improved PCA, including illustrative graphics | ||
2 | |||
3 | import LinearAlgebra | ||
4 | import Graphics.Plot | ||
5 | import System.Directory(doesFileExist) | ||
6 | import System(system) | ||
7 | import Control.Monad(when) | ||
8 | |||
9 | type Vec = Vector Double | ||
10 | type Mat = Matrix Double | ||
11 | |||
12 | sumColumns m = constant 1 (rows m) <> m | ||
13 | |||
14 | -- Vector with the mean value of the columns of a Mat | ||
15 | mean x = sumColumns x / fromIntegral (rows x) | ||
16 | |||
17 | -- covariance matrix of a list of observations as rows of a matrix | ||
18 | cov x = (trans xc <> xc) / fromIntegral (rows x -1) | ||
19 | where xc = center x | ||
20 | center m = m - constant 1 (rows m) `outer` mean m | ||
21 | |||
22 | type Stat = (Vec, [Double], Mat) | ||
23 | -- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov) | ||
24 | stat :: Mat -> Stat | ||
25 | stat x = (m, toList s, trans v) where | ||
26 | m = mean x | ||
27 | (s,v) = eigS (cov x) | ||
28 | |||
29 | -- creates the compression and decompression functions from the desired reconstruction | ||
30 | -- quality and the statistics of a data set | ||
31 | pca :: Double -> Stat -> (Vec -> Vec , Vec -> Vec) | ||
32 | pca prec (m,s,v) = (encode,decode) | ||
33 | where | ||
34 | encode x = vp <> (x - m) | ||
35 | decode x = x <> vp + m | ||
36 | vp = takeRows n v | ||
37 | n = 1 + (length $ fst $ span (< (prec'*sum s)) $ cumSum s) | ||
38 | cumSum = tail . scanl (+) 0.0 | ||
39 | prec' = if prec <=0.0 || prec >= 1.0 | ||
40 | then error "the precision in pca must be 0<prec<1" | ||
41 | else prec | ||
42 | |||
43 | shdigit :: Vec -> IO () | ||
44 | shdigit v = imshow (reshape 28 (-v)) | ||
45 | |||
46 | -- shows the effect of a given reconstruction quality on a test vector | ||
47 | test :: Stat -> Double -> Vec -> IO () | ||
48 | test st prec x = do | ||
49 | let (pe,pd) = pca prec st | ||
50 | let y = pe x | ||
51 | print $ dim y | ||
52 | shdigit (pd y) | ||
53 | |||
54 | main = do | ||
55 | ok <- doesFileExist ("mnist.txt") | ||
56 | when (not ok) $ do | ||
57 | putStrLn "\nTrying to download test datafile..." | ||
58 | system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz") | ||
59 | system("gunzip mnist.txt.gz") | ||
60 | return () | ||
61 | m <- fromFile "mnist.txt" (5000,785) | ||
62 | let xs = takeColumns (cols m -1) m | ||
63 | let x = toRows xs !! 4 -- an arbitrary test vector | ||
64 | shdigit x | ||
65 | let st = stat xs | ||
66 | test st 0.90 x | ||
67 | test st 0.50 x | ||
diff --git a/examples/pinv1.hs b/examples/pinv1.hs new file mode 100644 index 0000000..76fa0a9 --- /dev/null +++ b/examples/pinv1.hs | |||
@@ -0,0 +1,18 @@ | |||
1 | -- initial check for the polynomial model example | ||
2 | import LinearAlgebra | ||
3 | |||
4 | |||
5 | prepSyst :: Int -> Matrix Double -> (Matrix Double, Vector Double) | ||
6 | prepSyst n d = (a,b) where | ||
7 | [x,b] = toColumns d | ||
8 | a = fromColumns $ 1+0*x : map (x^) [1 .. n] | ||
9 | |||
10 | main = do | ||
11 | dat <- readMatrix `fmap` readFile "data.txt" | ||
12 | let (a,b) = prepSyst 3 dat | ||
13 | putStr "Coefficient matrix:\n" | ||
14 | dispR 3 a | ||
15 | putStr "Desired values:\n" | ||
16 | print b | ||
17 | putStr "\nLeast Squares Solution:\n" | ||
18 | print $ pinv a <> b -- equivalent to: linearSolveLSR a (fromColumns [b]) | ||
diff --git a/examples/pinv2.hs b/examples/pinv2.hs new file mode 100644 index 0000000..c1038d1 --- /dev/null +++ b/examples/pinv2.hs | |||
@@ -0,0 +1,26 @@ | |||
1 | -- MSE polynomial model using the pseudoinverse | ||
2 | import LinearAlgebra | ||
3 | import Graphics.Plot | ||
4 | |||
5 | expand :: Int -> Vector Double -> Matrix Double | ||
6 | expand n x = fromColumns $ constant 1 (dim x): map (x^) [1 .. n] | ||
7 | |||
8 | polynomialModel :: Matrix Double -> Int -> (Vector Double -> Vector Double) | ||
9 | polynomialModel d n = f where | ||
10 | f z = expand n z <> ws | ||
11 | ws = pinv a <> b | ||
12 | [x,b] = toColumns d | ||
13 | a = expand n x | ||
14 | |||
15 | mk fv = f where | ||
16 | f x = fv (fromList [x]) @> 0 | ||
17 | |||
18 | main = do | ||
19 | d <- readMatrix `fmap` readFile "data.txt" | ||
20 | let [x,y] = toColumns d | ||
21 | let pol = polynomialModel d | ||
22 | let view = [x, y, pol 1 x, pol 2 x, pol 3 x] | ||
23 | dispR 2 $ fromColumns view | ||
24 | mplot view | ||
25 | let f = mk (pol 2) | ||
26 | print (f 2.5) | ||
diff --git a/examples/plot.hs b/examples/plot.hs new file mode 100644 index 0000000..1177c11 --- /dev/null +++ b/examples/plot.hs | |||
@@ -0,0 +1,20 @@ | |||
1 | import LinearAlgebra | ||
2 | import Graphics.Plot | ||
3 | import GSL(erf_Z, erf) | ||
4 | |||
5 | sombrero n = f x y where | ||
6 | (x,y) = meshdom range range | ||
7 | range = linspace n (-2,2) | ||
8 | f x y = exp (-r2) * cos (2*r2) where | ||
9 | r2 = x*x+y*y | ||
10 | |||
11 | f x = sin x + 0.5 * sin (5*x) | ||
12 | |||
13 | gaussianPDF = erf_Z | ||
14 | cumdist x = 0.5 * (1+ erf (x/sqrt 2)) | ||
15 | |||
16 | main = do | ||
17 | let x = linspace 1000 (-4,4) | ||
18 | mplot [f x] | ||
19 | mplot [x, liftVector cumdist x, liftVector gaussianPDF x] | ||
20 | mesh (sombrero 40) \ No newline at end of file | ||
diff --git a/examples/testmnist.m b/examples/testmnist.m new file mode 100644 index 0000000..38625a7 --- /dev/null +++ b/examples/testmnist.m | |||
@@ -0,0 +1,29 @@ | |||
1 | #! /usr/bin/octave -qf | ||
2 | % measuring Octave computing times | ||
3 | |||
4 | t0=time(); | ||
5 | load mnist.txt | ||
6 | disp("load"); | ||
7 | disp(time()-t0) | ||
8 | |||
9 | |||
10 | x = mnist(:,1:784); | ||
11 | d = mnist(:,785); | ||
12 | |||
13 | |||
14 | t0=time(); | ||
15 | xc = x - repmat(mean(x),rows(x),1); | ||
16 | disp("x - repmat(mean(x),rows(x),1)"); | ||
17 | disp(time()-t0) | ||
18 | |||
19 | t0=time(); | ||
20 | mc = (xc'*xc)/rows(x); | ||
21 | disp("(xc'*xc)/rows(x)"); | ||
22 | disp(time()-t0) | ||
23 | |||
24 | t0=time(); | ||
25 | [v,l]=eig(mc); | ||
26 | disp("eig"); | ||
27 | disp(time()-t0) | ||
28 | |||
29 | disp(flipud(diag(l))(1:10)); \ No newline at end of file | ||