From 1925c123d7d8184a1d2ddc0a413e0fd2776e1083 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Thu, 8 May 2014 08:48:12 +0200 Subject: empty hmatrix-base --- examples/bool.hs | 54 --------------- examples/data.txt | 10 --- examples/deriv.hs | 8 --- examples/devel/ej1/functions.c | 35 ---------- examples/devel/ej1/wrappers.hs | 44 ------------ examples/devel/ej2/functions.c | 24 ------- examples/devel/ej2/wrappers.hs | 32 --------- examples/error.hs | 21 ------ examples/fitting.hs | 24 ------- examples/inplace.hs | 152 ----------------------------------------- examples/integrate.hs | 24 ------- examples/kalman.hs | 51 -------------- examples/lie.hs | 65 ------------------ examples/minimize.hs | 50 -------------- examples/monadic.hs | 118 -------------------------------- examples/multiply.hs | 104 ---------------------------- examples/ode.hs | 50 -------------- examples/parallel.hs | 28 -------- examples/pca1.hs | 46 ------------- examples/pca2.hs | 65 ------------------ examples/pinv.hs | 20 ------ examples/plot.hs | 20 ------ examples/root.hs | 31 --------- examples/vector.hs | 31 --------- 24 files changed, 1107 deletions(-) delete mode 100644 examples/bool.hs delete mode 100644 examples/data.txt delete mode 100644 examples/deriv.hs delete mode 100644 examples/devel/ej1/functions.c delete mode 100644 examples/devel/ej1/wrappers.hs delete mode 100644 examples/devel/ej2/functions.c delete mode 100644 examples/devel/ej2/wrappers.hs delete mode 100644 examples/error.hs delete mode 100644 examples/fitting.hs delete mode 100644 examples/inplace.hs delete mode 100644 examples/integrate.hs delete mode 100644 examples/kalman.hs delete mode 100644 examples/lie.hs delete mode 100644 examples/minimize.hs delete mode 100644 examples/monadic.hs delete mode 100644 examples/multiply.hs delete mode 100644 examples/ode.hs delete mode 100644 examples/parallel.hs delete mode 100644 examples/pca1.hs delete mode 100644 examples/pca2.hs delete mode 100644 examples/pinv.hs delete mode 100644 examples/plot.hs delete mode 100644 examples/root.hs delete mode 100644 examples/vector.hs (limited to 'examples') diff --git a/examples/bool.hs b/examples/bool.hs deleted file mode 100644 index 679b8bf..0000000 --- a/examples/bool.hs +++ /dev/null @@ -1,54 +0,0 @@ --- vectorized boolean operations defined in terms of step or cond - -import Numeric.LinearAlgebra - -infix 4 .==., ./=., .<., .<=., .>=., .>. -infixr 3 .&&. -infixr 2 .||. - -a .<. b = step (b-a) -a .<=. b = cond a b 1 1 0 -a .==. b = cond a b 0 1 0 -a ./=. b = cond a b 1 0 1 -a .>=. b = cond a b 0 1 1 -a .>. b = step (a-b) - -a .&&. b = step (a*b) -a .||. b = step (a+b) -no a = 1-a -xor a b = a ./=. b -equiv a b = a .==. b -imp a b = no a .||. b - -taut x = minElement x == 1 - -minEvery a b = cond a b a a b -maxEvery a b = cond a b b b a - --- examples - -clip a b x = cond y b y y b where y = cond x a a x x - -disp = putStr . dispf 3 - -eye n = ident n :: Matrix Double -row = asRow . fromList :: [Double] -> Matrix Double -col = asColumn . fromList :: [Double] -> Matrix Double - -m = (3><4) [1..] :: Matrix Double - -p = row [0,0,1,1] -q = row [0,1,0,1] - -main = do - print $ find (>6) m - disp $ assoc (6,8) 7 $ zip (find (/=0) (eye 5)) [10..] - disp $ accum (eye 5) (+) [((0,2),3), ((3,1),7), ((1,1),1)] - disp $ m .>=. 10 .||. m .<. 4 - (disp . fromColumns . map flatten) [p, q, p.&&.q, p .||.q, p `xor` q, p `equiv` q, p `imp` q] - print $ taut $ (p `imp` q ) `equiv` (no q `imp` no p) - print $ taut $ (xor p q) `equiv` (p .&&. no q .||. no p .&&. q) - disp $ clip 3 8 m - disp $ col [1..7] .<=. row [1..5] - disp $ cond (col [1..3]) (row [1..4]) m 50 (3*m) - diff --git a/examples/data.txt b/examples/data.txt deleted file mode 100644 index 2b9bfea..0000000 --- a/examples/data.txt +++ /dev/null @@ -1,10 +0,0 @@ - 0.9 1.1 - 2.1 3.9 - 3.1 9.2 - 4.0 51.8 - 4.9 25.3 - 6.1 35.7 - 7.0 49.4 - 7.9 3.6 - 9.1 81.5 -10.2 99.5 \ No newline at end of file diff --git a/examples/deriv.hs b/examples/deriv.hs deleted file mode 100644 index c9456d1..0000000 --- a/examples/deriv.hs +++ /dev/null @@ -1,8 +0,0 @@ --- Numerical differentiation - -import Numeric.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/devel/ej1/functions.c b/examples/devel/ej1/functions.c deleted file mode 100644 index 02a4cdd..0000000 --- a/examples/devel/ej1/functions.c +++ /dev/null @@ -1,35 +0,0 @@ -/* assuming row order */ - -typedef struct { double r, i; } doublecomplex; - -#define DVEC(A) int A##n, double*A##p -#define CVEC(A) int A##n, doublecomplex*A##p -#define DMAT(A) int A##r, int A##c, double*A##p -#define CMAT(A) int A##r, int A##c, doublecomplex*A##p - -#define AT(M,row,col) (M##p[(row)*M##c + (col)]) - -/*-----------------------------------------------------*/ - -int c_scale_vector(double s, DVEC(x), DVEC(y)) { - int k; - for (k=0; k<=yn; k++) { - yp[k] = s*xp[k]; - } - return 0; -} - -/*-----------------------------------------------------*/ - -int c_diag(DMAT(m),DVEC(y),DMAT(z)) { - int i,j; - for (j=0; j<5) [1..] - ------------------------------------------------------ - -foreign import ccall unsafe "c_scale_vector" - cScaleVector :: Double -- scale - -> CInt -> Ptr Double -- argument - -> CInt -> Ptr Double -- result - -> IO CInt -- exit code - -myScale s x = unsafePerformIO $ do - y <- createVector (dim x) - app2 (cScaleVector s) vec x vec y "cScaleVector" - return y - ------------------------------------------------------ --- forcing row order - -foreign import ccall unsafe "c_diag" - cDiag :: CInt -> CInt -> Ptr Double -- argument - -> CInt -> Ptr Double -- result1 - -> CInt -> CInt -> Ptr Double -- result2 - -> IO CInt -- exit code - -myDiag m = unsafePerformIO $ do - y <- createVector (min r c) - z <- createMatrix RowMajor r c - app3 cDiag mat (cmat m) vec y mat z "cDiag" - return (y,z) - where r = rows m - c = cols m diff --git a/examples/devel/ej2/functions.c b/examples/devel/ej2/functions.c deleted file mode 100644 index 4dcd377..0000000 --- a/examples/devel/ej2/functions.c +++ /dev/null @@ -1,24 +0,0 @@ -/* general element order */ - -typedef struct { double r, i; } doublecomplex; - -#define DVEC(A) int A##n, double*A##p -#define CVEC(A) int A##n, doublecomplex*A##p -#define DMAT(A) int A##r, int A##c, double*A##p -#define CMAT(A) int A##r, int A##c, doublecomplex*A##p - -#define AT(M,r,c) (M##p[(r)*sr+(c)*sc]) - -int c_diag(int ro, DMAT(m),DVEC(y),DMAT(z)) { - int i,j,sr,sc; - if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;} - for (j=0; j<5) [1..] - ------------------------------------------------------ --- arbitrary data order - -foreign import ccall unsafe "c_diag" - cDiag :: CInt -- matrix order - -> CInt -> CInt -> Ptr Double -- argument - -> CInt -> Ptr Double -- result1 - -> CInt -> CInt -> Ptr Double -- result2 - -> IO CInt -- exit code - -myDiag m = unsafePerformIO $ do - y <- createVector (min r c) - z <- createMatrix (orderOf m) r c - app3 (cDiag o) mat m vec y mat z "cDiag" - return (y,z) - where r = rows m - c = cols m - o = if orderOf m == RowMajor then 1 else 0 diff --git a/examples/error.hs b/examples/error.hs deleted file mode 100644 index 5efae7c..0000000 --- a/examples/error.hs +++ /dev/null @@ -1,21 +0,0 @@ -import Numeric.GSL -import Numeric.GSL.Special -import Numeric.LinearAlgebra -import Prelude hiding (catch) -import Control.Exception - -test x = catch - (print x) - (\e -> putStrLn $ "captured ["++ show (e :: SomeException) ++"]") - -main = do - setErrorHandlerOff - - test $ log_e (-1) - 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) - - putStrLn "Bye" \ No newline at end of file diff --git a/examples/fitting.hs b/examples/fitting.hs deleted file mode 100644 index a8f6b1c..0000000 --- a/examples/fitting.hs +++ /dev/null @@ -1,24 +0,0 @@ --- nonlinear least-squares fitting - -import Numeric.GSL.Fitting -import Numeric.LinearAlgebra - -xs = map return [0 .. 39] -sigma = 0.1 -ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs) - + scalar sigma * (randomVector 0 Gaussian 40) - -dat :: [([Double],([Double],Double))] - -dat = zip xs (zip ys (repeat sigma)) - -expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b] - -expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]] - -(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] - -main = do - print dat - print path - print sol diff --git a/examples/inplace.hs b/examples/inplace.hs deleted file mode 100644 index dcfff56..0000000 --- a/examples/inplace.hs +++ /dev/null @@ -1,152 +0,0 @@ --- some tests of the interface for pure --- computations with inplace updates - -import Numeric.LinearAlgebra -import Data.Packed.ST -import Data.Packed.Convert - -import Data.Array.Unboxed -import Data.Array.ST -import Control.Monad.ST -import Control.Monad - -main = sequence_[ - print test1, - print test2, - print test3, - print test4, - test5, - test6, - print test7, - test8, - test0] - --- helper functions -vector l = fromList l :: Vector Double -norm v = pnorm PNorm2 v - --- hmatrix vector and matrix -v = vector [1..10] -m = (5><10) [1..50::Double] - ----------------------------------------------------------------------- - --- vector creation by in-place updates on a copy of the argument -test1 = fun v - -fun :: Element t => Vector t -> Vector t -fun x = runSTVector $ do - a <- thawVector x - mapM_ (flip (modifyVector a) (+57)) [0 .. dim x `div` 2 - 1] - return a - --- another example: creation of an antidiagonal matrix from a list -test2 = antiDiag 5 8 [1..] :: Matrix Double - -antiDiag :: (Element b) => Int -> Int -> [b] -> Matrix b -antiDiag r c l = runSTMatrix $ do - m <- newMatrix 0 r c - let d = min r c - 1 - sequence_ $ zipWith (\i v -> writeMatrix m i (c-1-i) v) [0..d] l - return m - --- using vector or matrix functions on mutable objects requires freezing: -test3 = g1 v - -g1 x = runST $ do - a <- thawVector x - writeVector a (dim x -1) 0 - b <- freezeVector a - return (norm b) - --- another possibility: -test4 = g2 v - -g2 x = runST $ do - a <- thawVector x - writeVector a (dim x -1) 0 - t <- liftSTVector norm a - return t - --------------------------------------------------------------- - --- haskell arrays -hv = listArray (0,9) [1..10::Double] -hm = listArray ((0,0),(4,9)) [1..50::Double] - - - --- conversion from standard Haskell arrays -test5 = do - print $ norm (vectorFromArray hv) - print $ norm v - print $ rcond (matrixFromArray hm) - print $ rcond m - - --- conversion to mutable ST arrays -test6 = do - let y = clearColumn m 1 - print y - print (matrixFromArray y) - -clearColumn x c = runSTUArray $ do - a <- mArrayFromMatrix x - forM_ [0..rows x-1] $ \i-> - writeArray a (i,c) (0::Double) - return a - --- hmatrix functions applied to mutable ST arrays -test7 = unitary (listArray (1,4) [3,5,7,2] :: UArray Int Double) - -unitary v = runSTUArray $ do - a <- thaw v - n <- norm `fmap` vectorFromMArray a - b <- mapArray (/n) a - return b - -------------------------------------------------- - --- (just to check that they are not affected) -test0 = do - print v - print m - --print hv - --print hm - -------------------------------------------------- - -histogram n ds = runSTVector $ do - h <- newVector (0::Double) n -- number of bins - let inc k = modifyVector h k (+1) - mapM_ inc ds - return h - --- check that newVector is really called with a fresh new array -histoCheck ds = runSTVector $ do - h <- newVector (0::Double) 15 -- > constant for this test - let inc k = modifyVector h k (+1) - mapM_ inc ds - return h - -hc = fromList [1 .. 15::Double] - --- check that thawVector creates a new array -histoCheck2 ds = runSTVector $ do - h <- thawVector hc - let inc k = modifyVector h k (+1) - mapM_ inc ds - return h - -test8 = do - let ds = [0..14] - print $ histogram 15 ds - print $ histogram 15 ds - print $ histogram 15 ds - print $ histoCheck ds - print $ histoCheck ds - print $ histoCheck ds - print $ histoCheck2 ds - print $ histoCheck2 ds - print $ histoCheck2 ds - putStrLn "----------------------" diff --git a/examples/integrate.hs b/examples/integrate.hs deleted file mode 100644 index 3a03da6..0000000 --- a/examples/integrate.hs +++ /dev/null @@ -1,24 +0,0 @@ --- Numerical integration -import Numeric.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 y1 y2 g1 g2 = quad h y1 y2 - where - h y = quad (flip f y) (g1 y) (g2 y) - -volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y)) - 0 r (const 0) (\x->sqrt (r*r-x*x)) - --- wikipedia example -exw = quad2 f 7 10 (const 11) (const 14) - where - f x y = x**2 + 4*y - -main = do - print $ quad (\x -> 4/(x^2+1)) 0 1 - print pi - print $ volSphere 2.5 - print $ 4/3*pi*2.5**3 - print $ exw diff --git a/examples/kalman.hs b/examples/kalman.hs deleted file mode 100644 index 7fbe3d2..0000000 --- a/examples/kalman.hs +++ /dev/null @@ -1,51 +0,0 @@ -import Numeric.LinearAlgebra -import Graphics.Plot - -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 + 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 - -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)) diff --git a/examples/lie.hs b/examples/lie.hs deleted file mode 100644 index db21ea8..0000000 --- a/examples/lie.hs +++ /dev/null @@ -1,65 +0,0 @@ --- The magic of Lie Algebra - -import Numeric.LinearAlgebra - -disp = putStrLn . dispf 5 - -rot1 :: Double -> Matrix Double -rot1 a = (3><3) - [ 1, 0, 0 - , 0, c, s - , 0,-s, c ] - where c = cos a - s = sin a - -g1,g2,g3 :: Matrix Double - -g1 = (3><3) [0, 0,0 - ,0, 0,1 - ,0,-1,0] - -rot2 :: Double -> Matrix Double -rot2 a = (3><3) - [ c, 0, s - , 0, 1, 0 - ,-s, 0, c ] - where c = cos a - s = sin a - -g2 = (3><3) [ 0,0,1 - , 0,0,0 - ,-1,0,0] - -rot3 :: Double -> Matrix Double -rot3 a = (3><3) - [ c, s, 0 - ,-s, c, 0 - , 0, 0, 1 ] - where c = cos a - s = sin a - -g3 = (3><3) [ 0,1,0 - ,-1,0,0 - , 0,0,0] - -deg=pi/180 - --- commutator -infix 8 & -a & b = a <> b - b <> a - -infixl 6 |+| -a |+| b = a + b + a&b /2 + (a-b)&(a & b) /12 - -main = do - let a = 45*deg - b = 50*deg - c = -30*deg - 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 - putStrLn "exponential map back to the group (2 terms):" - disp (expm lie) - putStrLn "exact position:" - disp exact diff --git a/examples/minimize.hs b/examples/minimize.hs deleted file mode 100644 index 19b2cb3..0000000 --- a/examples/minimize.hs +++ /dev/null @@ -1,50 +0,0 @@ --- the multidimensional minimization example in the GSL manual -import Numeric.GSL -import Numeric.LinearAlgebra -import Graphics.Plot -import Text.Printf(printf) - --- the function to be minimized -f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30 - --- exact gradient -df [x,y] = [20*(x-1), 40*(y-2)] - --- a minimization algorithm which does not require the gradient -minimizeS f xi = minimize NMSimplex2 1E-2 100 (replicate (length xi) 1) f xi - --- 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 - -disp = putStrLn . format " " (printf "%.3f") - -allMethods :: (Enum a, Bounded a) => [a] -allMethods = [minBound .. maxBound] - -test method = do - print method - let (s,p) = minimize method 1E-2 30 [1,1] f [5,7] - print s - 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 - -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 - -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 deleted file mode 100644 index 7c6e0f4..0000000 --- a/examples/monadic.hs +++ /dev/null @@ -1,118 +0,0 @@ --- monadic computations --- (contributed by Vivian McPhail) - -import Numeric.LinearAlgebra -import Control.Monad.State.Strict -import Control.Monad.Maybe -import Foreign.Storable(Storable) -import System.Random(randomIO) - -------------------------------------------- - --- an instance of MonadIO, a monad transformer -type VectorMonadT = StateT Int IO - -test1 :: Vector Int -> IO (Vector Int) -test1 = mapVectorM $ \x -> do - putStr $ (show x) ++ " " - return (x + 1) - --- we can have an arbitrary monad AND do IO -addInitialM :: Vector Int -> 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 = 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 - put $ v + x - v' <- get - liftIO $ putStr $ (show v') ++ " " - - -indexPlusSum :: Vector Int -> VectorMonadT () -indexPlusSum v' = do - let f i x = do - s <- get - let inc = x+s - liftIO $ putStr $ show (i,inc) ++ " " - put inc - return inc - v <- mapVectorWithIndexM f v' - liftIO $ do - putStrLn "" - putStrLn $ show v - -------------------------------------------- - --- short circuit -monoStep :: Double -> MaybeT (State Double) () -monoStep d = do - dp <- get - when (d < dp) (fail "negative difference") - put d -{-# INLINE monoStep #-} - -isMonotoneIncreasing :: Vector Double -> Bool -isMonotoneIncreasing v = - let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v @> 0) - in case res of - Nothing -> False - Just _ -> True - - -------------------------------------------- - --- | 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) - where step e = do - ep <- lift $ get - if t e ep - then lift $ put e - 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) - where step e = do - ep <- get - put e - return $ f ep e - -------------------------------------------- - -v :: Vector Int -v = 10 |> [0..] - -w = fromList ([1..10]++[10,9..1]) :: Vector Double - - -main = do - v' <- test1 v - putStrLn "" - putStrLn $ show v' - evalStateT (addInitialM v) 0 - putStrLn "" - putStrLn $ show (sumEvens v) - evalStateT (sumEvensAndPrint v) 0 - putStrLn "" - evalStateT (indexPlusSum v) 0 - putStrLn "-----------------------" - mapVectorM_ print v - print =<< (mapVectorM (const randomIO) v :: IO (Vector Double)) - print =<< (mapVectorM (\a -> fmap (+a) randomIO) (5|>[0,100..1000]) :: IO (Vector Double)) - putStrLn "-----------------------" - print $ isMonotoneIncreasing w - print $ isMonotoneIncreasing (subVector 0 7 w) - print $ successive_ (>) v - print $ successive_ (>) w - print $ successive (+) v diff --git a/examples/multiply.hs b/examples/multiply.hs deleted file mode 100644 index 572961c..0000000 --- a/examples/multiply.hs +++ /dev/null @@ -1,104 +0,0 @@ -{-# LANGUAGE UnicodeSyntax - , MultiParamTypeClasses - , FunctionalDependencies - , FlexibleInstances - , FlexibleContexts --- , OverlappingInstances - , UndecidableInstances #-} - -import Numeric.LinearAlgebra - -class Scaling a b c | a b -> c where - -- ^ 0x22C5 8901 DOT OPERATOR, scaling - infixl 7 ⋅ - (⋅) :: a -> b -> c - -instance (Num t) => Scaling t t t where - (⋅) = (*) - -instance Container Vector t => Scaling t (Vector t) (Vector t) where - (⋅) = scale - -instance Container Vector t => Scaling (Vector t) t (Vector t) where - (⋅) = flip scale - -instance Container Vector t => Scaling t (Matrix t) (Matrix t) where - (⋅) = scale - -instance Container Vector t => Scaling (Matrix t) t (Matrix t) where - (⋅) = flip scale - - -class Mul a b c | a b -> c, a c -> b, b c -> a where - -- ^ 0x00D7 215 MULTIPLICATION SIGN ×, contraction - infixl 7 × - (×) :: a -> b -> c - - -------- - - - -instance Product t => Mul (Vector t) (Vector t) t where - (×) = udot - -instance Product t => Mul (Matrix t) (Vector t) (Vector t) where - (×) = mXv - -instance Product t => Mul (Vector t) (Matrix t) (Vector t) where - (×) = vXm - -instance Product t => Mul (Matrix t) (Matrix t) (Matrix t) where - (×) = mXm - - ---instance Scaling a b c => Contraction a b c where --- (×) = (⋅) - --------------------------------------------------------------------------------- - -class Outer a - where - infixl 7 ⊗ - -- | unicode 0x2297 8855 CIRCLED TIMES ⊗ - -- - -- vector outer product and matrix Kronecker product - (⊗) :: Product t => a t -> a t -> Matrix t - -instance Outer Vector where - (⊗) = outer - -instance Outer Matrix where - (⊗) = kronecker - --------------------------------------------------------------------------------- - - -v = 3 |> [1..] :: Vector Double - -m = (3 >< 3) [1..] :: Matrix Double - -s = 3 :: Double - -a = s ⋅ v × m × m × v ⋅ s - ---b = (v ⊗ m) ⊗ (v ⊗ m) - ---c = v ⊗ m ⊗ v ⊗ m - -d = s ⋅ (3 |> [10,20..] :: Vector Double) - -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 $ s ⋅ v × m × v - print a --- print (b == c) - print d - print $ asColumn u ⊗ w - print $ w ⊗ asColumn u - diff --git a/examples/ode.hs b/examples/ode.hs deleted file mode 100644 index dc6e0ec..0000000 --- a/examples/ode.hs +++ /dev/null @@ -1,50 +0,0 @@ -{-# LANGUAGE ViewPatterns #-} -import Numeric.GSL.ODE -import Numeric.LinearAlgebra -import Graphics.Plot -import Debug.Trace(trace) -debug x = trace (show x) x - -vanderpol mu = do - let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)] - ts = linspace 1000 (0,50) - sol = toColumns $ odeSolve (xdot mu) [1,0] ts - mplot (ts : sol) - mplot sol - - -harmonic w d = do - let xdot w d t [x,v] = [v, a*x + b*v] where a = -w^2; b = -2*d*w - ts = linspace 100 (0,20) - sol = odeSolve (xdot w d) [1,0] ts - mplot (ts : toColumns sol) - - -kepler v a = mplot (take 2 $ toColumns sol) where - xdot t [x,y,vx,vy] = [vx,vy,x*k,y*k] - where g=1 - k=(-g)*(x*x+y*y)**(-1.5) - ts = linspace 100 (0,30) - sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts - degree = pi/180 - - -main = do - vanderpol 2 - harmonic 1 0 - harmonic 1 0.1 - kepler 0.3 60 - kepler 0.4 70 - vanderpol' 2 - --- example of odeSolveV with jacobian -vanderpol' mu = do - let xdot mu t (toList->[x,v]) = fromList [v, -x + mu * v * (1-x^2)] - 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 - sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts - mplot sol - - diff --git a/examples/parallel.hs b/examples/parallel.hs deleted file mode 100644 index e875407..0000000 --- a/examples/parallel.hs +++ /dev/null @@ -1,28 +0,0 @@ --- $ ghc --make -O -rtsopts -threaded parallel.hs --- $ ./parallel 3000 +RTS -N4 -s -A200M - -import System.Environment(getArgs) -import Numeric.LinearAlgebra -import Control.Parallel.Strategies -import System.Time - -inParallel = parMap rwhnf id - --- matrix product decomposed into p parallel subtasks -parMul p x y = fromBlocks [ inParallel ( map (x <>) ys ) ] - where [ys] = toBlocksEvery (rows y) (cols y `div` p) y - -main = do - n <- (read . head) `fmap` getArgs - let m = ident n :: Matrix Double - time $ print $ maxElement $ takeDiag $ m <> m - time $ print $ maxElement $ takeDiag $ parMul 2 m m - time $ print $ maxElement $ takeDiag $ parMul 4 m m - time $ print $ maxElement $ takeDiag $ parMul 8 m m - -time act = do - t0 <- getClockTime - act - t1 <- getClockTime - print $ tdSec $ normalizeTimeDiff $ diffClockTimes t1 t0 - diff --git a/examples/pca1.hs b/examples/pca1.hs deleted file mode 100644 index a11eba9..0000000 --- a/examples/pca1.hs +++ /dev/null @@ -1,46 +0,0 @@ --- Principal component analysis - -import Numeric.LinearAlgebra -import System.Directory(doesFileExist) -import System.Process(system) -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 - -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 <- loadMatrix "mnist.txt" -- 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 deleted file mode 100644 index e7ea95f..0000000 --- a/examples/pca2.hs +++ /dev/null @@ -1,65 +0,0 @@ --- Improved PCA, including illustrative graphics - -import Numeric.LinearAlgebra -import Graphics.Plot -import System.Directory(doesFileExist) -import System.Process(system) -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) - - -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) = 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 - 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 <- loadMatrix "mnist.txt" - 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/pinv.hs b/examples/pinv.hs deleted file mode 100644 index 7de50b8..0000000 --- a/examples/pinv.hs +++ /dev/null @@ -1,20 +0,0 @@ -import Numeric.LinearAlgebra -import Graphics.Plot -import Text.Printf(printf) - -expand :: Int -> Vector Double -> Matrix Double -expand n x = fromColumns $ map (x^) [0 .. n] - -polynomialModel :: Vector Double -> Vector Double -> Int - -> (Vector Double -> Vector Double) -polynomialModel x y n = f where - f z = expand n z <> ws - ws = expand n x <\> y - -main = do - [x,y] <- (toColumns . readMatrix) `fmap` readFile "data.txt" - let pol = polynomialModel x y - let view = [x, y, pol 1 x, pol 2 x, pol 3 x] - putStrLn $ " x y p 1 p 2 p 3" - putStrLn $ format " " (printf "%.2f") $ fromColumns view - mplot view diff --git a/examples/plot.hs b/examples/plot.hs deleted file mode 100644 index f950aa5..0000000 --- a/examples/plot.hs +++ /dev/null @@ -1,20 +0,0 @@ -import Numeric.LinearAlgebra -import Graphics.Plot -import Numeric.GSL.Special(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, mapVector cumdist x, mapVector gaussianPDF x] - mesh (sombrero 40) \ No newline at end of file diff --git a/examples/root.hs b/examples/root.hs deleted file mode 100644 index 8546ff5..0000000 --- a/examples/root.hs +++ /dev/null @@ -1,31 +0,0 @@ --- root finding examples -import Numeric.GSL -import Numeric.LinearAlgebra -import Text.Printf(printf) - -rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ] - -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 - -jacobian a b [x,y] = [ [-a , 0] - , [-2*b*x, b] ] - -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 = putStrLn . format " " (printf "%.3f") - -main = do - test Hybrids - test Hybrid - test DNewton - test Broyden - - mapM_ testJ [HybridsJ .. GNewton] 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