add examples

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
+ 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<prec<1"
+                else prec
+
+shdigit :: Vec -> 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