From 0500032a1d954058b94cf9a0fa2a662e5666a526 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sun, 4 Oct 2015 14:16:57 +0200 Subject: update examples --- examples/error.hs | 8 +++++--- examples/inplace.hs | 4 +++- examples/kalman.hs | 33 ++++++++++++++++----------------- examples/lie.hs | 10 +++++----- examples/minimize.hs | 9 +++++---- examples/monadic.hs | 31 ++++++++++++++++++------------- examples/multiply.hs | 22 +++++++++++----------- examples/ode.hs | 2 +- examples/pca1.hs | 18 ++++++++---------- examples/pca2.hs | 17 ++++++++++------- examples/plot.hs | 4 ++-- examples/root.hs | 6 +++--- examples/vector.hs | 31 ------------------------------- 13 files changed, 87 insertions(+), 108 deletions(-) delete mode 100644 examples/vector.hs (limited to 'examples') 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 (print x) (\e -> putStrLn $ "captured ["++ show (e :: SomeException) ++"]") + main = do setErrorHandlerOff @@ -15,7 +16,8 @@ main = do test $ 5 + (fst.exp_e) 1000 test $ bessel_zero_Jnu_e (-0.3) 2 - test $ (linearSolve 0 4 :: Matrix Double) - test $ (linearSolve 5 (sqrt (-1)) :: Matrix Double) + test $ (inv 0 :: Matrix Double) + test $ (linearSolveLS 5 (sqrt (-1)) :: Matrix Double) + + putStrLn "Bye" - 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 @@ -- some tests of the interface for pure -- computations with inplace updates -import Numeric.LinearAlgebra.HMatrix +{-# LANGUAGE FlexibleContexts #-} + +import Numeric.LinearAlgebra import Numeric.LinearAlgebra.Devel 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 @@ import Numeric.LinearAlgebra import Graphics.Plot -vector l = fromList l :: Vector Double -matrix ls = fromLists ls :: Matrix Double -diagl = diag . vector +f = fromLists + [[1,0,0,0], + [1,1,0,0], + [0,0,1,0], + [0,0,0,1]] -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]] +h = fromLists + [[0,-1,1,0], + [0,-1,0,1]] q = diagl [1,1,0,0] @@ -25,13 +23,13 @@ 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 + q -- 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 + px = f #> x -- prediction + pq = f <> p <> tr f + q -- its covariance + y = z - h #> px -- residue + cy = h <> pq <> tr h + r -- its covariance + k = pq <> tr h <> inv cy -- kalman gain + x' = px + k #> y -- new state + p' = (ident (size x) - k <> h) <> pq -- its covariance sys = System f h q r @@ -49,3 +47,4 @@ main = do print $ fromRows $ take 10 (map sX xs) mapM_ (print . sP) $ take 10 xs mplot (evolution 20 (xs,des)) + 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 @@ -- The magic of Lie Algebra -import Numeric.LinearAlgebra +{-# LANGUAGE FlexibleContexts #-} -disp = putStrLn . dispf 5 +import Numeric.LinearAlgebra rot1 :: Double -> Matrix Double rot1 a = (3><3) @@ -58,8 +58,8 @@ main = do exact = rot3 a <> rot1 b <> rot2 c lie = scalar a * g3 |+| scalar b * g1 |+| scalar c * g2 putStrLn "position in the tangent space:" - disp lie + disp 5 lie putStrLn "exponential map back to the group (2 terms):" - disp (expm lie) + disp 5 (expm lie) putStrLn "exact position:" - disp exact + 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 g x = f (concat [a,x:b]) (a,_:b) = splitAt n v -disp = putStrLn . format " " (printf "%.3f") +disp' = putStrLn . format " " (printf "%.3f") allMethods :: (Enum a, Bounded a) => [a] allMethods = [minBound .. maxBound] @@ -29,22 +29,23 @@ test method = do print method let (s,p) = minimize method 1E-2 30 [1,1] f [5,7] print s - disp p + disp' p testD method = do print method let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f df [5,7] print s - disp p + disp' p testD' method = do putStrLn $ show method ++ " with estimated gradient" let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f (gradient f) [5,7] print s - disp p + disp' p main = do mapM_ test [NMSimplex, NMSimplex2] mapM_ testD allMethods testD' ConjugateFR mplot $ drop 3 . toColumns . snd $ minimizeS f [5,7] + 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 @@ -- monadic computations -- (contributed by Vivian McPhail) +{-# LANGUAGE FlexibleContexts #-} + import Numeric.LinearAlgebra +import Numeric.LinearAlgebra.Devel import Control.Monad.State.Strict -import Control.Monad.Maybe +import Control.Monad.Trans.Maybe import Foreign.Storable(Storable) import System.Random(randomIO) ------------------------------------------- -- an instance of MonadIO, a monad transformer -type VectorMonadT = StateT Int IO +type VectorMonadT = StateT I IO -test1 :: Vector Int -> IO (Vector Int) +test1 :: Vector I -> IO (Vector I) test1 = mapVectorM $ \x -> do putStr $ (show x) ++ " " return (x + 1) -- we can have an arbitrary monad AND do IO -addInitialM :: Vector Int -> VectorMonadT () +addInitialM :: Vector I -> VectorMonadT () addInitialM = mapVectorM_ $ \x -> do i <- get liftIO $ putStr $ (show $ x + i) ++ " " put $ x + i -- sum the values of the even indiced elements -sumEvens :: Vector Int -> Int +sumEvens :: Vector I -> I sumEvens = foldVectorWithIndex (\x a b -> if x `mod` 2 == 0 then a + b else b) 0 -- sum and print running total of evens -sumEvensAndPrint :: Vector Int -> VectorMonadT () sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do when (i `mod` 2 == 0) $ do v <- get @@ -38,7 +40,7 @@ sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do liftIO $ putStr $ (show v') ++ " " -indexPlusSum :: Vector Int -> VectorMonadT () +--indexPlusSum :: Vector I -> VectorMonadT () indexPlusSum v' = do let f i x = do s <- get @@ -63,7 +65,7 @@ monoStep d = do isMonotoneIncreasing :: Vector Double -> Bool isMonotoneIncreasing v = - let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v @> 0) + let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v ! 0) in case res of Nothing -> False Just _ -> True @@ -72,8 +74,8 @@ isMonotoneIncreasing v = ------------------------------------------- -- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs -successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool -successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0) +successive_ :: (Container Vector a, Indexable (Vector a) a) => (a -> a -> Bool) -> Vector a -> Bool +successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (size v - 1) v))) (v ! 0) where step e = do ep <- lift $ get if t e ep @@ -81,8 +83,10 @@ successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ s else (fail "successive_ test failed") -- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input -successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b -successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0) +successive + :: (Storable b, Container Vector s, Indexable (Vector s) s) + => (s -> s -> b) -> Vector s -> Vector b +successive f v = evalState (mapVectorM step (subVector 1 (size v - 1) v)) (v ! 0) where step e = do ep <- get put e @@ -90,7 +94,7 @@ successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0 ------------------------------------------- -v :: Vector Int +v :: Vector I v = 10 |> [0..] w = fromList ([1..10]++[10,9..1]) :: Vector Double @@ -116,3 +120,4 @@ main = do print $ successive_ (>) v print $ successive_ (>) w print $ successive (+) v + 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 instance Container Vector t => Scaling (Vector t) t (Vector t) where (⋅) = flip scale -instance Container Vector t => Scaling t (Matrix t) (Matrix t) where +instance (Num t, Container Vector t) => Scaling t (Matrix t) (Matrix t) where (⋅) = scale -instance Container Vector t => Scaling (Matrix t) t (Matrix t) where +instance (Num t, Container Vector t) => Scaling (Matrix t) t (Matrix t) where (⋅) = flip scale @@ -42,14 +42,14 @@ class Mul a b c | a b -> c, a c -> b, b c -> a where instance Product t => Mul (Vector t) (Vector t) t where (×) = udot -instance Product t => Mul (Matrix t) (Vector t) (Vector t) where - (×) = mXv +instance (Numeric t, Product t) => Mul (Matrix t) (Vector t) (Vector t) where + (×) = (#>) -instance Product t => Mul (Vector t) (Matrix t) (Vector t) where - (×) = vXm +instance (Numeric t, Product t) => Mul (Vector t) (Matrix t) (Vector t) where + (×) = (<#) -instance Product t => Mul (Matrix t) (Matrix t) (Matrix t) where - (×) = mXm +instance (Numeric t, Product t) => Mul (Matrix t) (Matrix t) (Matrix t) where + (×) = (<>) --instance Scaling a b c => Contraction a b c where @@ -92,9 +92,9 @@ u = fromList [3,0,5] w = konst 1 (2,3) :: Matrix Double main = do - print $ (scale s v <> m) `udot` v - print $ scale s v `udot` (m <> v) - print $ s * ((v <> m) `udot` v) + print $ (scale s v <# m) `udot` v + print $ scale s v `udot` (m #> v) + print $ s * ((v <# m) `udot` v) print $ s ⋅ v × m × v print a -- 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 jac t (toList->[x,v]) = (2><2) [ 0 , 1 , -1-2*x*v*mu, mu*(1-x**2) ] ts = linspace 1000 (0,50) - hi = (ts@>1 - ts@>0)/100 + hi = (ts!1 - ts!0)/100 sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts mplot sol 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) type Vec = Vector Double type Mat = Matrix Double - +{- -- Vector with the mean value of the columns of a matrix mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a -- covariance matrix of a list of observations stored as rows cov x = (trans xc <> xc) / fromIntegral (rows x - 1) where xc = x - asRow (mean x) +-} -- 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) = eigSH' c - vp = takeRows n (trans v) - -norm = pnorm PNorm2 + encode x = vp #> (x - m) + decode x = x <# vp + m + (m,c) = meanCov dataSet + (_,v) = eigSH (trustSym c) + vp = tr $ takeColumns n v main = do ok <- doesFileExist ("mnist.txt") @@ -43,4 +41,4 @@ main = do let (pe,pd) = pca 10 xs let y = pe x print y -- compressed version - print $ norm (x - pd y) / norm x --reconstruction quality + 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 @@ -- Improved PCA, including illustrative graphics +{-# LANGUAGE FlexibleContexts #-} + import Numeric.LinearAlgebra import Graphics.Plot import System.Directory(doesFileExist) @@ -10,27 +12,27 @@ type Vec = Vector Double type Mat = Matrix Double -- Vector with the mean value of the columns of a matrix -mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a +mean a = konst (recip . fromIntegral . rows $ a) (rows a) <# a -- covariance matrix of a list of observations stored as rows -cov x = (trans xc <> xc) / fromIntegral (rows x - 1) +cov x = (mTm xc) -- / fromIntegral (rows x - 1) where xc = x - asRow (mean x) 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 +stat x = (m, toList s, tr v) where m = mean x - (s,v) = eigSH' (cov x) + (s,v) = eigSH (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 + 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 @@ -46,7 +48,7 @@ test :: Stat -> Double -> Vec -> IO () test st prec x = do let (pe,pd) = pca prec st let y = pe x - print $ dim y + print $ size y shdigit (pd y) main = do @@ -63,3 +65,4 @@ main = do let st = stat xs test st 0.90 x test st 0.50 x + 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)) main = do let x = linspace 1000 (-4,4) mplot [f x] - mplot [x, mapVector cumdist x, mapVector gaussianPDF x] - mesh (sombrero 40) \ No newline at end of file + mplot [x, cmap cumdist x, cmap gaussianPDF x] + 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 print method let (s,p) = root method 1E-7 30 (rosenbrock 1 10) [-10,-5] print s -- solution - disp p -- evolution of the algorithm + disp' p -- evolution of the algorithm jacobian a b [x,y] = [ [-a , 0] , [-2*b*x, b] ] @@ -18,9 +18,9 @@ testJ method = do print method let (s,p) = rootJ method 1E-7 30 (rosenbrock 1 10) (jacobian 1 10) [-10,-5] print s - disp p + disp' p -disp = putStrLn . format " " (printf "%.3f") +disp' = putStrLn . format " " (printf "%.3f") main = do 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 @@ --- conversion to/from Data.Vector.Storable --- from Roman Leshchinskiy "vector" package --- --- In the future Data.Packed.Vector will be replaced by Data.Vector.Storable - -------------------------------------------- - -import Numeric.LinearAlgebra as H -import Data.Packed.Development(unsafeFromForeignPtr, unsafeToForeignPtr) -import Foreign.Storable -import qualified Data.Vector.Storable as V - -fromVector :: Storable t => V.Vector t -> H.Vector t -fromVector v = unsafeFromForeignPtr p i n where - (p,i,n) = V.unsafeToForeignPtr v - -toVector :: Storable t => H.Vector t -> V.Vector t -toVector v = V.unsafeFromForeignPtr p i n where - (p,i,n) = unsafeToForeignPtr v - -------------------------------------------- - -v = V.slice 5 10 (V.fromList [1 .. 10::Double] V.++ V.replicate 10 7) - -w = subVector 2 3 (linspace 5 (0,1)) :: Vector Double - -main = do - print v - print $ fromVector v - print w - print $ toVector w -- cgit v1.2.3