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