summaryrefslogtreecommitdiff
path: root/packages/hmatrix/examples
diff options
context:
space:
mode:
authorAlberto Ruiz <aruiz@um.es>2014-05-21 10:30:55 +0200
committerAlberto Ruiz <aruiz@um.es>2014-05-21 10:30:55 +0200
commit197e88c3b56d28840217010a2871c6ea3a4dd1a4 (patch)
tree825be9d6c9d87d23f7e5497c0133d11d52c63535 /packages/hmatrix/examples
parente07c3dee7235496b71a89233106d93f6cc94ada1 (diff)
update dependencies, move examples etc
Diffstat (limited to 'packages/hmatrix/examples')
-rw-r--r--packages/hmatrix/examples/bool.hs54
-rw-r--r--packages/hmatrix/examples/data.txt10
-rw-r--r--packages/hmatrix/examples/deriv.hs8
-rw-r--r--packages/hmatrix/examples/devel/ej1/functions.c35
-rw-r--r--packages/hmatrix/examples/devel/ej1/wrappers.hs44
-rw-r--r--packages/hmatrix/examples/devel/ej2/functions.c24
-rw-r--r--packages/hmatrix/examples/devel/ej2/wrappers.hs32
-rw-r--r--packages/hmatrix/examples/error.hs21
-rw-r--r--packages/hmatrix/examples/fitting.hs24
-rw-r--r--packages/hmatrix/examples/inplace.hs152
-rw-r--r--packages/hmatrix/examples/integrate.hs24
-rw-r--r--packages/hmatrix/examples/kalman.hs51
-rw-r--r--packages/hmatrix/examples/lie.hs65
-rw-r--r--packages/hmatrix/examples/minimize.hs50
-rw-r--r--packages/hmatrix/examples/monadic.hs118
-rw-r--r--packages/hmatrix/examples/multiply.hs104
-rw-r--r--packages/hmatrix/examples/ode.hs50
-rw-r--r--packages/hmatrix/examples/parallel.hs28
-rw-r--r--packages/hmatrix/examples/pca1.hs46
-rw-r--r--packages/hmatrix/examples/pca2.hs65
-rw-r--r--packages/hmatrix/examples/pinv.hs20
-rw-r--r--packages/hmatrix/examples/plot.hs20
-rw-r--r--packages/hmatrix/examples/root.hs31
-rw-r--r--packages/hmatrix/examples/vector.hs31
24 files changed, 0 insertions, 1107 deletions
diff --git a/packages/hmatrix/examples/bool.hs b/packages/hmatrix/examples/bool.hs
deleted file mode 100644
index 679b8bf..0000000
--- a/packages/hmatrix/examples/bool.hs
+++ /dev/null
@@ -1,54 +0,0 @@
1-- vectorized boolean operations defined in terms of step or cond
2
3import Numeric.LinearAlgebra
4
5infix 4 .==., ./=., .<., .<=., .>=., .>.
6infixr 3 .&&.
7infixr 2 .||.
8
9a .<. b = step (b-a)
10a .<=. b = cond a b 1 1 0
11a .==. b = cond a b 0 1 0
12a ./=. b = cond a b 1 0 1
13a .>=. b = cond a b 0 1 1
14a .>. b = step (a-b)
15
16a .&&. b = step (a*b)
17a .||. b = step (a+b)
18no a = 1-a
19xor a b = a ./=. b
20equiv a b = a .==. b
21imp a b = no a .||. b
22
23taut x = minElement x == 1
24
25minEvery a b = cond a b a a b
26maxEvery a b = cond a b b b a
27
28-- examples
29
30clip a b x = cond y b y y b where y = cond x a a x x
31
32disp = putStr . dispf 3
33
34eye n = ident n :: Matrix Double
35row = asRow . fromList :: [Double] -> Matrix Double
36col = asColumn . fromList :: [Double] -> Matrix Double
37
38m = (3><4) [1..] :: Matrix Double
39
40p = row [0,0,1,1]
41q = row [0,1,0,1]
42
43main = do
44 print $ find (>6) m
45 disp $ assoc (6,8) 7 $ zip (find (/=0) (eye 5)) [10..]
46 disp $ accum (eye 5) (+) [((0,2),3), ((3,1),7), ((1,1),1)]
47 disp $ m .>=. 10 .||. m .<. 4
48 (disp . fromColumns . map flatten) [p, q, p.&&.q, p .||.q, p `xor` q, p `equiv` q, p `imp` q]
49 print $ taut $ (p `imp` q ) `equiv` (no q `imp` no p)
50 print $ taut $ (xor p q) `equiv` (p .&&. no q .||. no p .&&. q)
51 disp $ clip 3 8 m
52 disp $ col [1..7] .<=. row [1..5]
53 disp $ cond (col [1..3]) (row [1..4]) m 50 (3*m)
54
diff --git a/packages/hmatrix/examples/data.txt b/packages/hmatrix/examples/data.txt
deleted file mode 100644
index 2b9bfea..0000000
--- a/packages/hmatrix/examples/data.txt
+++ /dev/null
@@ -1,10 +0,0 @@
1 0.9 1.1
2 2.1 3.9
3 3.1 9.2
4 4.0 51.8
5 4.9 25.3
6 6.1 35.7
7 7.0 49.4
8 7.9 3.6
9 9.1 81.5
1010.2 99.5 \ No newline at end of file
diff --git a/packages/hmatrix/examples/deriv.hs b/packages/hmatrix/examples/deriv.hs
deleted file mode 100644
index c9456d1..0000000
--- a/packages/hmatrix/examples/deriv.hs
+++ /dev/null
@@ -1,8 +0,0 @@
1-- Numerical differentiation
2
3import Numeric.GSL
4
5d :: (Double -> Double) -> (Double -> Double)
6d f x = fst $ derivCentral 0.01 f x
7
8main = print $ d (\x-> x * d (\y-> x+y) 1) 1
diff --git a/packages/hmatrix/examples/devel/ej1/functions.c b/packages/hmatrix/examples/devel/ej1/functions.c
deleted file mode 100644
index 02a4cdd..0000000
--- a/packages/hmatrix/examples/devel/ej1/functions.c
+++ /dev/null
@@ -1,35 +0,0 @@
1/* assuming row order */
2
3typedef struct { double r, i; } doublecomplex;
4
5#define DVEC(A) int A##n, double*A##p
6#define CVEC(A) int A##n, doublecomplex*A##p
7#define DMAT(A) int A##r, int A##c, double*A##p
8#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
9
10#define AT(M,row,col) (M##p[(row)*M##c + (col)])
11
12/*-----------------------------------------------------*/
13
14int c_scale_vector(double s, DVEC(x), DVEC(y)) {
15 int k;
16 for (k=0; k<=yn; k++) {
17 yp[k] = s*xp[k];
18 }
19 return 0;
20}
21
22/*-----------------------------------------------------*/
23
24int c_diag(DMAT(m),DVEC(y),DMAT(z)) {
25 int i,j;
26 for (j=0; j<yn; j++) {
27 yp[j] = AT(m,j,j);
28 }
29 for (i=0; i<mr; i++) {
30 for (j=0; j<mc; j++) {
31 AT(z,i,j) = i==j?yp[i]:0;
32 }
33 }
34 return 0;
35}
diff --git a/packages/hmatrix/examples/devel/ej1/wrappers.hs b/packages/hmatrix/examples/devel/ej1/wrappers.hs
deleted file mode 100644
index a88f74b..0000000
--- a/packages/hmatrix/examples/devel/ej1/wrappers.hs
+++ /dev/null
@@ -1,44 +0,0 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2
3-- $ ghc -O2 --make wrappers.hs functions.c
4
5import Numeric.LinearAlgebra
6import Data.Packed.Development
7import Foreign(Ptr,unsafePerformIO)
8import Foreign.C.Types(CInt)
9
10-----------------------------------------------------
11
12main = do
13 print $ myScale 3.0 (fromList [1..10])
14 print $ myDiag $ (3><5) [1..]
15
16-----------------------------------------------------
17
18foreign import ccall unsafe "c_scale_vector"
19 cScaleVector :: Double -- scale
20 -> CInt -> Ptr Double -- argument
21 -> CInt -> Ptr Double -- result
22 -> IO CInt -- exit code
23
24myScale s x = unsafePerformIO $ do
25 y <- createVector (dim x)
26 app2 (cScaleVector s) vec x vec y "cScaleVector"
27 return y
28
29-----------------------------------------------------
30-- forcing row order
31
32foreign import ccall unsafe "c_diag"
33 cDiag :: CInt -> CInt -> Ptr Double -- argument
34 -> CInt -> Ptr Double -- result1
35 -> CInt -> CInt -> Ptr Double -- result2
36 -> IO CInt -- exit code
37
38myDiag m = unsafePerformIO $ do
39 y <- createVector (min r c)
40 z <- createMatrix RowMajor r c
41 app3 cDiag mat (cmat m) vec y mat z "cDiag"
42 return (y,z)
43 where r = rows m
44 c = cols m
diff --git a/packages/hmatrix/examples/devel/ej2/functions.c b/packages/hmatrix/examples/devel/ej2/functions.c
deleted file mode 100644
index 4dcd377..0000000
--- a/packages/hmatrix/examples/devel/ej2/functions.c
+++ /dev/null
@@ -1,24 +0,0 @@
1/* general element order */
2
3typedef struct { double r, i; } doublecomplex;
4
5#define DVEC(A) int A##n, double*A##p
6#define CVEC(A) int A##n, doublecomplex*A##p
7#define DMAT(A) int A##r, int A##c, double*A##p
8#define CMAT(A) int A##r, int A##c, doublecomplex*A##p
9
10#define AT(M,r,c) (M##p[(r)*sr+(c)*sc])
11
12int c_diag(int ro, DMAT(m),DVEC(y),DMAT(z)) {
13 int i,j,sr,sc;
14 if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;}
15 for (j=0; j<yn; j++) {
16 yp[j] = AT(m,j,j);
17 }
18 for (i=0; i<mr; i++) {
19 for (j=0; j<mc; j++) {
20 AT(z,i,j) = i==j?yp[i]:0;
21 }
22 }
23 return 0;
24}
diff --git a/packages/hmatrix/examples/devel/ej2/wrappers.hs b/packages/hmatrix/examples/devel/ej2/wrappers.hs
deleted file mode 100644
index 1c02a24..0000000
--- a/packages/hmatrix/examples/devel/ej2/wrappers.hs
+++ /dev/null
@@ -1,32 +0,0 @@
1{-# LANGUAGE ForeignFunctionInterface #-}
2
3-- $ ghc -O2 --make wrappers.hs functions.c
4
5import Numeric.LinearAlgebra
6import Data.Packed.Development
7import Foreign(Ptr,unsafePerformIO)
8import Foreign.C.Types(CInt)
9
10-----------------------------------------------------
11
12main = do
13 print $ myDiag $ (3><5) [1..]
14
15-----------------------------------------------------
16-- arbitrary data order
17
18foreign import ccall unsafe "c_diag"
19 cDiag :: CInt -- matrix order
20 -> CInt -> CInt -> Ptr Double -- argument
21 -> CInt -> Ptr Double -- result1
22 -> CInt -> CInt -> Ptr Double -- result2
23 -> IO CInt -- exit code
24
25myDiag m = unsafePerformIO $ do
26 y <- createVector (min r c)
27 z <- createMatrix (orderOf m) r c
28 app3 (cDiag o) mat m vec y mat z "cDiag"
29 return (y,z)
30 where r = rows m
31 c = cols m
32 o = if orderOf m == RowMajor then 1 else 0
diff --git a/packages/hmatrix/examples/error.hs b/packages/hmatrix/examples/error.hs
deleted file mode 100644
index 5efae7c..0000000
--- a/packages/hmatrix/examples/error.hs
+++ /dev/null
@@ -1,21 +0,0 @@
1import Numeric.GSL
2import Numeric.GSL.Special
3import Numeric.LinearAlgebra
4import Prelude hiding (catch)
5import Control.Exception
6
7test x = catch
8 (print x)
9 (\e -> putStrLn $ "captured ["++ show (e :: SomeException) ++"]")
10
11main = do
12 setErrorHandlerOff
13
14 test $ log_e (-1)
15 test $ 5 + (fst.exp_e) 1000
16 test $ bessel_zero_Jnu_e (-0.3) 2
17
18 test $ (linearSolve 0 4 :: Matrix Double)
19 test $ (linearSolve 5 (sqrt (-1)) :: Matrix Double)
20
21 putStrLn "Bye" \ No newline at end of file
diff --git a/packages/hmatrix/examples/fitting.hs b/packages/hmatrix/examples/fitting.hs
deleted file mode 100644
index a8f6b1c..0000000
--- a/packages/hmatrix/examples/fitting.hs
+++ /dev/null
@@ -1,24 +0,0 @@
1-- nonlinear least-squares fitting
2
3import Numeric.GSL.Fitting
4import Numeric.LinearAlgebra
5
6xs = map return [0 .. 39]
7sigma = 0.1
8ys = map return $ toList $ fromList (map (head . expModel [5,0.1,1]) xs)
9 + scalar sigma * (randomVector 0 Gaussian 40)
10
11dat :: [([Double],([Double],Double))]
12
13dat = zip xs (zip ys (repeat sigma))
14
15expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]
16
17expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
18
19(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
20
21main = do
22 print dat
23 print path
24 print sol
diff --git a/packages/hmatrix/examples/inplace.hs b/packages/hmatrix/examples/inplace.hs
deleted file mode 100644
index dcfff56..0000000
--- a/packages/hmatrix/examples/inplace.hs
+++ /dev/null
@@ -1,152 +0,0 @@
1-- some tests of the interface for pure
2-- computations with inplace updates
3
4import Numeric.LinearAlgebra
5import Data.Packed.ST
6import Data.Packed.Convert
7
8import Data.Array.Unboxed
9import Data.Array.ST
10import Control.Monad.ST
11import Control.Monad
12
13main = sequence_[
14 print test1,
15 print test2,
16 print test3,
17 print test4,
18 test5,
19 test6,
20 print test7,
21 test8,
22 test0]
23
24-- helper functions
25vector l = fromList l :: Vector Double
26norm v = pnorm PNorm2 v
27
28-- hmatrix vector and matrix
29v = vector [1..10]
30m = (5><10) [1..50::Double]
31
32----------------------------------------------------------------------
33
34-- vector creation by in-place updates on a copy of the argument
35test1 = fun v
36
37fun :: Element t => Vector t -> Vector t
38fun x = runSTVector $ do
39 a <- thawVector x
40 mapM_ (flip (modifyVector a) (+57)) [0 .. dim x `div` 2 - 1]
41 return a
42
43-- another example: creation of an antidiagonal matrix from a list
44test2 = antiDiag 5 8 [1..] :: Matrix Double
45
46antiDiag :: (Element b) => Int -> Int -> [b] -> Matrix b
47antiDiag r c l = runSTMatrix $ do
48 m <- newMatrix 0 r c
49 let d = min r c - 1
50 sequence_ $ zipWith (\i v -> writeMatrix m i (c-1-i) v) [0..d] l
51 return m
52
53-- using vector or matrix functions on mutable objects requires freezing:
54test3 = g1 v
55
56g1 x = runST $ do
57 a <- thawVector x
58 writeVector a (dim x -1) 0
59 b <- freezeVector a
60 return (norm b)
61
62-- another possibility:
63test4 = g2 v
64
65g2 x = runST $ do
66 a <- thawVector x
67 writeVector a (dim x -1) 0
68 t <- liftSTVector norm a
69 return t
70
71--------------------------------------------------------------
72
73-- haskell arrays
74hv = listArray (0,9) [1..10::Double]
75hm = listArray ((0,0),(4,9)) [1..50::Double]
76
77
78
79-- conversion from standard Haskell arrays
80test5 = do
81 print $ norm (vectorFromArray hv)
82 print $ norm v
83 print $ rcond (matrixFromArray hm)
84 print $ rcond m
85
86
87-- conversion to mutable ST arrays
88test6 = do
89 let y = clearColumn m 1
90 print y
91 print (matrixFromArray y)
92
93clearColumn x c = runSTUArray $ do
94 a <- mArrayFromMatrix x
95 forM_ [0..rows x-1] $ \i->
96 writeArray a (i,c) (0::Double)
97 return a
98
99-- hmatrix functions applied to mutable ST arrays
100test7 = unitary (listArray (1,4) [3,5,7,2] :: UArray Int Double)
101
102unitary v = runSTUArray $ do
103 a <- thaw v
104 n <- norm `fmap` vectorFromMArray a
105 b <- mapArray (/n) a
106 return b
107
108-------------------------------------------------
109
110-- (just to check that they are not affected)
111test0 = do
112 print v
113 print m
114 --print hv
115 --print hm
116
117-------------------------------------------------
118
119histogram n ds = runSTVector $ do
120 h <- newVector (0::Double) n -- number of bins
121 let inc k = modifyVector h k (+1)
122 mapM_ inc ds
123 return h
124
125-- check that newVector is really called with a fresh new array
126histoCheck ds = runSTVector $ do
127 h <- newVector (0::Double) 15 -- > constant for this test
128 let inc k = modifyVector h k (+1)
129 mapM_ inc ds
130 return h
131
132hc = fromList [1 .. 15::Double]
133
134-- check that thawVector creates a new array
135histoCheck2 ds = runSTVector $ do
136 h <- thawVector hc
137 let inc k = modifyVector h k (+1)
138 mapM_ inc ds
139 return h
140
141test8 = do
142 let ds = [0..14]
143 print $ histogram 15 ds
144 print $ histogram 15 ds
145 print $ histogram 15 ds
146 print $ histoCheck ds
147 print $ histoCheck ds
148 print $ histoCheck ds
149 print $ histoCheck2 ds
150 print $ histoCheck2 ds
151 print $ histoCheck2 ds
152 putStrLn "----------------------"
diff --git a/packages/hmatrix/examples/integrate.hs b/packages/hmatrix/examples/integrate.hs
deleted file mode 100644
index 3a03da6..0000000
--- a/packages/hmatrix/examples/integrate.hs
+++ /dev/null
@@ -1,24 +0,0 @@
1-- Numerical integration
2import Numeric.GSL
3
4quad f a b = fst $ integrateQAGS 1E-9 100 f a b
5
6-- A multiple integral can be easily defined using partial application
7quad2 f y1 y2 g1 g2 = quad h y1 y2
8 where
9 h y = quad (flip f y) (g1 y) (g2 y)
10
11volSphere r = 8 * quad2 (\x y -> sqrt (r*r-x*x-y*y))
12 0 r (const 0) (\x->sqrt (r*r-x*x))
13
14-- wikipedia example
15exw = quad2 f 7 10 (const 11) (const 14)
16 where
17 f x y = x**2 + 4*y
18
19main = do
20 print $ quad (\x -> 4/(x^2+1)) 0 1
21 print pi
22 print $ volSphere 2.5
23 print $ 4/3*pi*2.5**3
24 print $ exw
diff --git a/packages/hmatrix/examples/kalman.hs b/packages/hmatrix/examples/kalman.hs
deleted file mode 100644
index 7fbe3d2..0000000
--- a/packages/hmatrix/examples/kalman.hs
+++ /dev/null
@@ -1,51 +0,0 @@
1import Numeric.LinearAlgebra
2import Graphics.Plot
3
4vector l = fromList l :: Vector Double
5matrix ls = fromLists ls :: Matrix Double
6diagl = diag . vector
7
8f = matrix [[1,0,0,0],
9 [1,1,0,0],
10 [0,0,1,0],
11 [0,0,0,1]]
12
13h = matrix [[0,-1,1,0],
14 [0,-1,0,1]]
15
16q = diagl [1,1,0,0]
17
18r = diagl [2,2]
19
20s0 = State (vector [0, 0, 10, -10]) (diagl [10,0, 100, 100])
21
22data System = System {kF, kH, kQ, kR :: Matrix Double}
23data State = State {sX :: Vector Double , sP :: Matrix Double}
24type Measurement = Vector Double
25
26kalman :: System -> State -> Measurement -> State
27kalman (System f h q r) (State x p) z = State x' p' where
28 px = f <> x -- prediction
29 pq = f <> p <> trans f + q -- its covariance
30 y = z - h <> px -- residue
31 cy = h <> pq <> trans h + r -- its covariance
32 k = pq <> trans h <> inv cy -- kalman gain
33 x' = px + k <> y -- new state
34 p' = (ident (dim x) - k <> h) <> pq -- its covariance
35
36sys = System f h q r
37
38zs = [vector [15-k,-20-k] | k <- [0..]]
39
40xs = s0 : zipWith (kalman sys) xs zs
41
42des = map (sqrt.takeDiag.sP) xs
43
44evolution n (xs,des) =
45 vector [1.. fromIntegral n]:(toColumns $ fromRows $ take n (zipWith (-) (map sX xs) des)) ++
46 (toColumns $ fromRows $ take n (zipWith (+) (map sX xs) des))
47
48main = do
49 print $ fromRows $ take 10 (map sX xs)
50 mapM_ (print . sP) $ take 10 xs
51 mplot (evolution 20 (xs,des))
diff --git a/packages/hmatrix/examples/lie.hs b/packages/hmatrix/examples/lie.hs
deleted file mode 100644
index db21ea8..0000000
--- a/packages/hmatrix/examples/lie.hs
+++ /dev/null
@@ -1,65 +0,0 @@
1-- The magic of Lie Algebra
2
3import Numeric.LinearAlgebra
4
5disp = putStrLn . dispf 5
6
7rot1 :: Double -> Matrix Double
8rot1 a = (3><3)
9 [ 1, 0, 0
10 , 0, c, s
11 , 0,-s, c ]
12 where c = cos a
13 s = sin a
14
15g1,g2,g3 :: Matrix Double
16
17g1 = (3><3) [0, 0,0
18 ,0, 0,1
19 ,0,-1,0]
20
21rot2 :: Double -> Matrix Double
22rot2 a = (3><3)
23 [ c, 0, s
24 , 0, 1, 0
25 ,-s, 0, c ]
26 where c = cos a
27 s = sin a
28
29g2 = (3><3) [ 0,0,1
30 , 0,0,0
31 ,-1,0,0]
32
33rot3 :: Double -> Matrix Double
34rot3 a = (3><3)
35 [ c, s, 0
36 ,-s, c, 0
37 , 0, 0, 1 ]
38 where c = cos a
39 s = sin a
40
41g3 = (3><3) [ 0,1,0
42 ,-1,0,0
43 , 0,0,0]
44
45deg=pi/180
46
47-- commutator
48infix 8 &
49a & b = a <> b - b <> a
50
51infixl 6 |+|
52a |+| b = a + b + a&b /2 + (a-b)&(a & b) /12
53
54main = do
55 let a = 45*deg
56 b = 50*deg
57 c = -30*deg
58 exact = rot3 a <> rot1 b <> rot2 c
59 lie = scalar a * g3 |+| scalar b * g1 |+| scalar c * g2
60 putStrLn "position in the tangent space:"
61 disp lie
62 putStrLn "exponential map back to the group (2 terms):"
63 disp (expm lie)
64 putStrLn "exact position:"
65 disp exact
diff --git a/packages/hmatrix/examples/minimize.hs b/packages/hmatrix/examples/minimize.hs
deleted file mode 100644
index 19b2cb3..0000000
--- a/packages/hmatrix/examples/minimize.hs
+++ /dev/null
@@ -1,50 +0,0 @@
1-- the multidimensional minimization example in the GSL manual
2import Numeric.GSL
3import Numeric.LinearAlgebra
4import Graphics.Plot
5import Text.Printf(printf)
6
7-- the function to be minimized
8f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
9
10-- exact gradient
11df [x,y] = [20*(x-1), 40*(y-2)]
12
13-- a minimization algorithm which does not require the gradient
14minimizeS f xi = minimize NMSimplex2 1E-2 100 (replicate (length xi) 1) f xi
15
16-- Numerical estimation of the gradient
17gradient f v = [partialDerivative k f v | k <- [0 .. length v -1]]
18
19partialDerivative n f v = fst (derivCentral 0.01 g (v!!n)) where
20 g x = f (concat [a,x:b])
21 (a,_:b) = splitAt n v
22
23disp = putStrLn . format " " (printf "%.3f")
24
25allMethods :: (Enum a, Bounded a) => [a]
26allMethods = [minBound .. maxBound]
27
28test method = do
29 print method
30 let (s,p) = minimize method 1E-2 30 [1,1] f [5,7]
31 print s
32 disp p
33
34testD method = do
35 print method
36 let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f df [5,7]
37 print s
38 disp p
39
40testD' method = do
41 putStrLn $ show method ++ " with estimated gradient"
42 let (s,p) = minimizeD method 1E-3 30 1E-2 1E-4 f (gradient f) [5,7]
43 print s
44 disp p
45
46main = do
47 mapM_ test [NMSimplex, NMSimplex2]
48 mapM_ testD allMethods
49 testD' ConjugateFR
50 mplot $ drop 3 . toColumns . snd $ minimizeS f [5,7]
diff --git a/packages/hmatrix/examples/monadic.hs b/packages/hmatrix/examples/monadic.hs
deleted file mode 100644
index 7c6e0f4..0000000
--- a/packages/hmatrix/examples/monadic.hs
+++ /dev/null
@@ -1,118 +0,0 @@
1-- monadic computations
2-- (contributed by Vivian McPhail)
3
4import Numeric.LinearAlgebra
5import Control.Monad.State.Strict
6import Control.Monad.Maybe
7import Foreign.Storable(Storable)
8import System.Random(randomIO)
9
10-------------------------------------------
11
12-- an instance of MonadIO, a monad transformer
13type VectorMonadT = StateT Int IO
14
15test1 :: Vector Int -> IO (Vector Int)
16test1 = mapVectorM $ \x -> do
17 putStr $ (show x) ++ " "
18 return (x + 1)
19
20-- we can have an arbitrary monad AND do IO
21addInitialM :: Vector Int -> VectorMonadT ()
22addInitialM = mapVectorM_ $ \x -> do
23 i <- get
24 liftIO $ putStr $ (show $ x + i) ++ " "
25 put $ x + i
26
27-- sum the values of the even indiced elements
28sumEvens :: Vector Int -> Int
29sumEvens = foldVectorWithIndex (\x a b -> if x `mod` 2 == 0 then a + b else b) 0
30
31-- sum and print running total of evens
32sumEvensAndPrint :: Vector Int -> VectorMonadT ()
33sumEvensAndPrint = mapVectorWithIndexM_ $ \ i x -> do
34 when (i `mod` 2 == 0) $ do
35 v <- get
36 put $ v + x
37 v' <- get
38 liftIO $ putStr $ (show v') ++ " "
39
40
41indexPlusSum :: Vector Int -> VectorMonadT ()
42indexPlusSum v' = do
43 let f i x = do
44 s <- get
45 let inc = x+s
46 liftIO $ putStr $ show (i,inc) ++ " "
47 put inc
48 return inc
49 v <- mapVectorWithIndexM f v'
50 liftIO $ do
51 putStrLn ""
52 putStrLn $ show v
53
54-------------------------------------------
55
56-- short circuit
57monoStep :: Double -> MaybeT (State Double) ()
58monoStep d = do
59 dp <- get
60 when (d < dp) (fail "negative difference")
61 put d
62{-# INLINE monoStep #-}
63
64isMonotoneIncreasing :: Vector Double -> Bool
65isMonotoneIncreasing v =
66 let res = evalState (runMaybeT $ (mapVectorM_ monoStep v)) (v @> 0)
67 in case res of
68 Nothing -> False
69 Just _ -> True
70
71
72-------------------------------------------
73
74-- | apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs
75successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool
76successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ step (subVector 1 (dim v - 1) v))) (v @> 0)
77 where step e = do
78 ep <- lift $ get
79 if t e ep
80 then lift $ put e
81 else (fail "successive_ test failed")
82
83-- | operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input
84successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b
85successive f v = evalState (mapVectorM step (subVector 1 (dim v - 1) v)) (v @> 0)
86 where step e = do
87 ep <- get
88 put e
89 return $ f ep e
90
91-------------------------------------------
92
93v :: Vector Int
94v = 10 |> [0..]
95
96w = fromList ([1..10]++[10,9..1]) :: Vector Double
97
98
99main = do
100 v' <- test1 v
101 putStrLn ""
102 putStrLn $ show v'
103 evalStateT (addInitialM v) 0
104 putStrLn ""
105 putStrLn $ show (sumEvens v)
106 evalStateT (sumEvensAndPrint v) 0
107 putStrLn ""
108 evalStateT (indexPlusSum v) 0
109 putStrLn "-----------------------"
110 mapVectorM_ print v
111 print =<< (mapVectorM (const randomIO) v :: IO (Vector Double))
112 print =<< (mapVectorM (\a -> fmap (+a) randomIO) (5|>[0,100..1000]) :: IO (Vector Double))
113 putStrLn "-----------------------"
114 print $ isMonotoneIncreasing w
115 print $ isMonotoneIncreasing (subVector 0 7 w)
116 print $ successive_ (>) v
117 print $ successive_ (>) w
118 print $ successive (+) v
diff --git a/packages/hmatrix/examples/multiply.hs b/packages/hmatrix/examples/multiply.hs
deleted file mode 100644
index 572961c..0000000
--- a/packages/hmatrix/examples/multiply.hs
+++ /dev/null
@@ -1,104 +0,0 @@
1{-# LANGUAGE UnicodeSyntax
2 , MultiParamTypeClasses
3 , FunctionalDependencies
4 , FlexibleInstances
5 , FlexibleContexts
6-- , OverlappingInstances
7 , UndecidableInstances #-}
8
9import Numeric.LinearAlgebra
10
11class Scaling a b c | a b -> c where
12 -- ^ 0x22C5 8901 DOT OPERATOR, scaling
13 infixl 7 â‹…
14 (â‹…) :: a -> b -> c
15
16instance (Num t) => Scaling t t t where
17 (â‹…) = (*)
18
19instance Container Vector t => Scaling t (Vector t) (Vector t) where
20 (â‹…) = scale
21
22instance Container Vector t => Scaling (Vector t) t (Vector t) where
23 (â‹…) = flip scale
24
25instance Container Vector t => Scaling t (Matrix t) (Matrix t) where
26 (â‹…) = scale
27
28instance Container Vector t => Scaling (Matrix t) t (Matrix t) where
29 (â‹…) = flip scale
30
31
32class Mul a b c | a b -> c, a c -> b, b c -> a where
33 -- ^ 0x00D7 215 MULTIPLICATION SIGN ×, contraction
34 infixl 7 ×
35 (×) :: a -> b -> c
36
37
38-------
39
40
41
42instance Product t => Mul (Vector t) (Vector t) t where
43 (×) = udot
44
45instance Product t => Mul (Matrix t) (Vector t) (Vector t) where
46 (×) = mXv
47
48instance Product t => Mul (Vector t) (Matrix t) (Vector t) where
49 (×) = vXm
50
51instance Product t => Mul (Matrix t) (Matrix t) (Matrix t) where
52 (×) = mXm
53
54
55--instance Scaling a b c => Contraction a b c where
56-- (×) = (⋅)
57
58--------------------------------------------------------------------------------
59
60class Outer a
61 where
62 infixl 7 ⊗
63 -- | unicode 0x2297 8855 CIRCLED TIMES ⊗
64 --
65 -- vector outer product and matrix Kronecker product
66 (⊗) :: Product t => a t -> a t -> Matrix t
67
68instance Outer Vector where
69 (⊗) = outer
70
71instance Outer Matrix where
72 (⊗) = kronecker
73
74--------------------------------------------------------------------------------
75
76
77v = 3 |> [1..] :: Vector Double
78
79m = (3 >< 3) [1..] :: Matrix Double
80
81s = 3 :: Double
82
83a = s ⋅ v × m × m × v ⋅ s
84
85--b = (v ⊗ m) ⊗ (v ⊗ m)
86
87--c = v ⊗ m ⊗ v ⊗ m
88
89d = s â‹… (3 |> [10,20..] :: Vector Double)
90
91u = fromList [3,0,5]
92w = konst 1 (2,3) :: Matrix Double
93
94main = do
95 print $ (scale s v <> m) `udot` v
96 print $ scale s v `udot` (m <> v)
97 print $ s * ((v <> m) `udot` v)
98 print $ s ⋅ v × m × v
99 print a
100-- print (b == c)
101 print d
102 print $ asColumn u ⊗ w
103 print $ w ⊗ asColumn u
104
diff --git a/packages/hmatrix/examples/ode.hs b/packages/hmatrix/examples/ode.hs
deleted file mode 100644
index dc6e0ec..0000000
--- a/packages/hmatrix/examples/ode.hs
+++ /dev/null
@@ -1,50 +0,0 @@
1{-# LANGUAGE ViewPatterns #-}
2import Numeric.GSL.ODE
3import Numeric.LinearAlgebra
4import Graphics.Plot
5import Debug.Trace(trace)
6debug x = trace (show x) x
7
8vanderpol mu = do
9 let xdot mu t [x,v] = [v, -x + mu * v * (1-x^2)]
10 ts = linspace 1000 (0,50)
11 sol = toColumns $ odeSolve (xdot mu) [1,0] ts
12 mplot (ts : sol)
13 mplot sol
14
15
16harmonic w d = do
17 let xdot w d t [x,v] = [v, a*x + b*v] where a = -w^2; b = -2*d*w
18 ts = linspace 100 (0,20)
19 sol = odeSolve (xdot w d) [1,0] ts
20 mplot (ts : toColumns sol)
21
22
23kepler v a = mplot (take 2 $ toColumns sol) where
24 xdot t [x,y,vx,vy] = [vx,vy,x*k,y*k]
25 where g=1
26 k=(-g)*(x*x+y*y)**(-1.5)
27 ts = linspace 100 (0,30)
28 sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts
29 degree = pi/180
30
31
32main = do
33 vanderpol 2
34 harmonic 1 0
35 harmonic 1 0.1
36 kepler 0.3 60
37 kepler 0.4 70
38 vanderpol' 2
39
40-- example of odeSolveV with jacobian
41vanderpol' mu = do
42 let xdot mu t (toList->[x,v]) = fromList [v, -x + mu * v * (1-x^2)]
43 jac t (toList->[x,v]) = (2><2) [ 0 , 1
44 , -1-2*x*v*mu, mu*(1-x**2) ]
45 ts = linspace 1000 (0,50)
46 hi = (ts@>1 - ts@>0)/100
47 sol = toColumns $ odeSolveV (MSBDF jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts
48 mplot sol
49
50
diff --git a/packages/hmatrix/examples/parallel.hs b/packages/hmatrix/examples/parallel.hs
deleted file mode 100644
index e875407..0000000
--- a/packages/hmatrix/examples/parallel.hs
+++ /dev/null
@@ -1,28 +0,0 @@
1-- $ ghc --make -O -rtsopts -threaded parallel.hs
2-- $ ./parallel 3000 +RTS -N4 -s -A200M
3
4import System.Environment(getArgs)
5import Numeric.LinearAlgebra
6import Control.Parallel.Strategies
7import System.Time
8
9inParallel = parMap rwhnf id
10
11-- matrix product decomposed into p parallel subtasks
12parMul p x y = fromBlocks [ inParallel ( map (x <>) ys ) ]
13 where [ys] = toBlocksEvery (rows y) (cols y `div` p) y
14
15main = do
16 n <- (read . head) `fmap` getArgs
17 let m = ident n :: Matrix Double
18 time $ print $ maxElement $ takeDiag $ m <> m
19 time $ print $ maxElement $ takeDiag $ parMul 2 m m
20 time $ print $ maxElement $ takeDiag $ parMul 4 m m
21 time $ print $ maxElement $ takeDiag $ parMul 8 m m
22
23time act = do
24 t0 <- getClockTime
25 act
26 t1 <- getClockTime
27 print $ tdSec $ normalizeTimeDiff $ diffClockTimes t1 t0
28
diff --git a/packages/hmatrix/examples/pca1.hs b/packages/hmatrix/examples/pca1.hs
deleted file mode 100644
index a11eba9..0000000
--- a/packages/hmatrix/examples/pca1.hs
+++ /dev/null
@@ -1,46 +0,0 @@
1-- Principal component analysis
2
3import Numeric.LinearAlgebra
4import System.Directory(doesFileExist)
5import System.Process(system)
6import Control.Monad(when)
7
8type Vec = Vector Double
9type Mat = Matrix Double
10
11
12-- Vector with the mean value of the columns of a matrix
13mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
14
15-- covariance matrix of a list of observations stored as rows
16cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
17 where xc = x - asRow (mean x)
18
19
20-- creates the compression and decompression functions from the desired number of components
21pca :: Int -> Mat -> (Vec -> Vec , Vec -> Vec)
22pca n dataSet = (encode,decode)
23 where
24 encode x = vp <> (x - m)
25 decode x = x <> vp + m
26 m = mean dataSet
27 c = cov dataSet
28 (_,v) = eigSH' c
29 vp = takeRows n (trans v)
30
31norm = pnorm PNorm2
32
33main = do
34 ok <- doesFileExist ("mnist.txt")
35 when (not ok) $ do
36 putStrLn "\nTrying to download test datafile..."
37 system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
38 system("gunzip mnist.txt.gz")
39 return ()
40 m <- loadMatrix "mnist.txt" -- fromFile "mnist.txt" (5000,785)
41 let xs = takeColumns (cols m -1) m -- the last column is the digit type (class label)
42 let x = toRows xs !! 4 -- an arbitrary test Vec
43 let (pe,pd) = pca 10 xs
44 let y = pe x
45 print y -- compressed version
46 print $ norm (x - pd y) / norm x --reconstruction quality
diff --git a/packages/hmatrix/examples/pca2.hs b/packages/hmatrix/examples/pca2.hs
deleted file mode 100644
index e7ea95f..0000000
--- a/packages/hmatrix/examples/pca2.hs
+++ /dev/null
@@ -1,65 +0,0 @@
1-- Improved PCA, including illustrative graphics
2
3import Numeric.LinearAlgebra
4import Graphics.Plot
5import System.Directory(doesFileExist)
6import System.Process(system)
7import Control.Monad(when)
8
9type Vec = Vector Double
10type Mat = Matrix Double
11
12-- Vector with the mean value of the columns of a matrix
13mean a = constant (recip . fromIntegral . rows $ a) (rows a) <> a
14
15-- covariance matrix of a list of observations stored as rows
16cov x = (trans xc <> xc) / fromIntegral (rows x - 1)
17 where xc = x - asRow (mean x)
18
19
20type Stat = (Vec, [Double], Mat)
21-- 1st and 2nd order statistics of a dataset (mean, eigenvalues and eigenvectors of cov)
22stat :: Mat -> Stat
23stat x = (m, toList s, trans v) where
24 m = mean x
25 (s,v) = eigSH' (cov x)
26
27-- creates the compression and decompression functions from the desired reconstruction
28-- quality and the statistics of a data set
29pca :: Double -> Stat -> (Vec -> Vec , Vec -> Vec)
30pca prec (m,s,v) = (encode,decode)
31 where
32 encode x = vp <> (x - m)
33 decode x = x <> vp + m
34 vp = takeRows n v
35 n = 1 + (length $ fst $ span (< (prec'*sum s)) $ cumSum s)
36 cumSum = tail . scanl (+) 0.0
37 prec' = if prec <=0.0 || prec >= 1.0
38 then error "the precision in pca must be 0<prec<1"
39 else prec
40
41shdigit :: Vec -> IO ()
42shdigit v = imshow (reshape 28 (-v))
43
44-- shows the effect of a given reconstruction quality on a test vector
45test :: Stat -> Double -> Vec -> IO ()
46test st prec x = do
47 let (pe,pd) = pca prec st
48 let y = pe x
49 print $ dim y
50 shdigit (pd y)
51
52main = do
53 ok <- doesFileExist ("mnist.txt")
54 when (not ok) $ do
55 putStrLn "\nTrying to download test datafile..."
56 system("wget -nv http://dis.um.es/~alberto/material/sp/mnist.txt.gz")
57 system("gunzip mnist.txt.gz")
58 return ()
59 m <- loadMatrix "mnist.txt"
60 let xs = takeColumns (cols m -1) m
61 let x = toRows xs !! 4 -- an arbitrary test vector
62 shdigit x
63 let st = stat xs
64 test st 0.90 x
65 test st 0.50 x
diff --git a/packages/hmatrix/examples/pinv.hs b/packages/hmatrix/examples/pinv.hs
deleted file mode 100644
index 7de50b8..0000000
--- a/packages/hmatrix/examples/pinv.hs
+++ /dev/null
@@ -1,20 +0,0 @@
1import Numeric.LinearAlgebra
2import Graphics.Plot
3import Text.Printf(printf)
4
5expand :: Int -> Vector Double -> Matrix Double
6expand n x = fromColumns $ map (x^) [0 .. n]
7
8polynomialModel :: Vector Double -> Vector Double -> Int
9 -> (Vector Double -> Vector Double)
10polynomialModel x y n = f where
11 f z = expand n z <> ws
12 ws = expand n x <\> y
13
14main = do
15 [x,y] <- (toColumns . readMatrix) `fmap` readFile "data.txt"
16 let pol = polynomialModel x y
17 let view = [x, y, pol 1 x, pol 2 x, pol 3 x]
18 putStrLn $ " x y p 1 p 2 p 3"
19 putStrLn $ format " " (printf "%.2f") $ fromColumns view
20 mplot view
diff --git a/packages/hmatrix/examples/plot.hs b/packages/hmatrix/examples/plot.hs
deleted file mode 100644
index f950aa5..0000000
--- a/packages/hmatrix/examples/plot.hs
+++ /dev/null
@@ -1,20 +0,0 @@
1import Numeric.LinearAlgebra
2import Graphics.Plot
3import Numeric.GSL.Special(erf_Z, erf)
4
5sombrero 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
11f x = sin x + 0.5 * sin (5*x)
12
13gaussianPDF = erf_Z
14cumdist x = 0.5 * (1+ erf (x/sqrt 2))
15
16main = do
17 let x = linspace 1000 (-4,4)
18 mplot [f x]
19 mplot [x, mapVector cumdist x, mapVector gaussianPDF x]
20 mesh (sombrero 40) \ No newline at end of file
diff --git a/packages/hmatrix/examples/root.hs b/packages/hmatrix/examples/root.hs
deleted file mode 100644
index 8546ff5..0000000
--- a/packages/hmatrix/examples/root.hs
+++ /dev/null
@@ -1,31 +0,0 @@
1-- root finding examples
2import Numeric.GSL
3import Numeric.LinearAlgebra
4import Text.Printf(printf)
5
6rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
7
8test method = do
9 print method
10 let (s,p) = root method 1E-7 30 (rosenbrock 1 10) [-10,-5]
11 print s -- solution
12 disp p -- evolution of the algorithm
13
14jacobian a b [x,y] = [ [-a , 0]
15 , [-2*b*x, b] ]
16
17testJ method = do
18 print method
19 let (s,p) = rootJ method 1E-7 30 (rosenbrock 1 10) (jacobian 1 10) [-10,-5]
20 print s
21 disp p
22
23disp = putStrLn . format " " (printf "%.3f")
24
25main = do
26 test Hybrids
27 test Hybrid
28 test DNewton
29 test Broyden
30
31 mapM_ testJ [HybridsJ .. GNewton]
diff --git a/packages/hmatrix/examples/vector.hs b/packages/hmatrix/examples/vector.hs
deleted file mode 100644
index f531cbd..0000000
--- a/packages/hmatrix/examples/vector.hs
+++ /dev/null
@@ -1,31 +0,0 @@
1-- conversion to/from Data.Vector.Storable
2-- from Roman Leshchinskiy "vector" package
3--
4-- In the future Data.Packed.Vector will be replaced by Data.Vector.Storable
5
6-------------------------------------------
7
8import Numeric.LinearAlgebra as H
9import Data.Packed.Development(unsafeFromForeignPtr, unsafeToForeignPtr)
10import Foreign.Storable
11import qualified Data.Vector.Storable as V
12
13fromVector :: Storable t => V.Vector t -> H.Vector t
14fromVector v = unsafeFromForeignPtr p i n where
15 (p,i,n) = V.unsafeToForeignPtr v
16
17toVector :: Storable t => H.Vector t -> V.Vector t
18toVector v = V.unsafeFromForeignPtr p i n where
19 (p,i,n) = unsafeToForeignPtr v
20
21-------------------------------------------
22
23v = V.slice 5 10 (V.fromList [1 .. 10::Double] V.++ V.replicate 10 7)
24
25w = subVector 2 3 (linspace 5 (0,1)) :: Vector Double
26
27main = do
28 print v
29 print $ fromVector v
30 print w
31 print $ toVector w
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-11 15:30:29 +0000