summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--HSSL.cabal2
-rw-r--r--examples/oldtests.hs121
-rw-r--r--examples/tests.hs2
-rw-r--r--lib/Data/Packed/Instances.hs391
-rw-r--r--lib/Data/Packed/Internal/Vector.hs3
-rw-r--r--lib/GSL.hs6
-rw-r--r--lib/GSL/Compat.hs370
-rw-r--r--lib/GSL/Minimization.hs8
8 files changed, 504 insertions, 399 deletions
diff --git a/HSSL.cabal b/HSSL.cabal
index 84ca5f5..1377232 100644
--- a/HSSL.cabal
+++ b/HSSL.cabal
@@ -24,7 +24,6 @@ Exposed-modules: Data.Packed.Internal,
24 Data.Packed.Internal.Matrix, Data.Packed.Matrix, 24 Data.Packed.Internal.Matrix, Data.Packed.Matrix,
25 Data.Packed.Internal.Tensor, Data.Packed.Tensor, 25 Data.Packed.Internal.Tensor, Data.Packed.Tensor,
26 Data.Packed.Plot 26 Data.Packed.Plot
27 Data.Packed.Instances
28 LAPACK, 27 LAPACK,
29 GSL.Vector, 28 GSL.Vector,
30 GSL.Matrix 29 GSL.Matrix
@@ -34,6 +33,7 @@ Exposed-modules: Data.Packed.Internal,
34 GSL.Polynomials, 33 GSL.Polynomials,
35 GSL.Minimization, 34 GSL.Minimization,
36 LinearAlgebra, LinearAlgebra.Algorithms, 35 LinearAlgebra, LinearAlgebra.Algorithms,
36 GSL.Compat,
37 GSL, HSSL 37 GSL, HSSL
38Other-modules: 38Other-modules:
39C-sources: lib/Data/Packed/Internal/aux.c, 39C-sources: lib/Data/Packed/Internal/aux.c,
diff --git a/examples/oldtests.hs b/examples/oldtests.hs
new file mode 100644
index 0000000..987ef98
--- /dev/null
+++ b/examples/oldtests.hs
@@ -0,0 +1,121 @@
1import Test.HUnit
2
3import GSL
4import GSL.Matrix
5import System.Random(randomRs,mkStdGen)
6
7realMatrix = fromLists :: [[Double]] -> Matrix Double
8realVector = fromList :: [Double] -> Vector Double
9
10toComplexM = uncurry $ liftMatrix2 (curry toComplex)
11
12infixl 2 =~=
13a =~= b = pnorm PNorm1 (flatten (a - b)) < 1E-6
14
15randomMatrix seed (n,m) = reshape m $ realVector $ take (n*m) $ randomRs (-100,100) $ mkStdGen seed
16
17randomMatrixC seed (n,m) = toComplexM (randomMatrix seed (n,m), randomMatrix (seed+1) (n,m))
18
19besselTest = do
20 let (r,e) = bessel_J0_e 5.0
21 let expected = -0.17759677131433830434739701
22 assertBool "bessel_J0_e" ( abs (r-expected) < e )
23
24exponentialTest = do
25 let (v,e,err) = exp_e10_e 30.0
26 let expected = exp 30.0
27 assertBool "exp_e10_e" ( abs (v*10^e - expected) < 4E-2 )
28
29disp m = putStrLn (format " " show m)
30
31ms = realMatrix [[1,2,3]
32 ,[-4,1,7]]
33
34ms' = randomMatrix 27 (50,100)
35
36ms'' = toComplexM (randomMatrix 100 (50,100),randomMatrix 101 (50,100))
37
38fullsvdTest method mat msg = do
39 let (u,s,vt) = method mat
40 assertBool msg (u <> s <> trans vt =~= mat)
41
42svdg' m = (u, diag s, v) where (u,s,v) = svdg m
43
44full_svd_Rd = svdRdd
45
46--------------------------------------------------------------------
47
48mcu = toComplexM (randomMatrix 33 (20,20),randomMatrix 34 (20,20))
49
50mcur = randomMatrix 35 (40,40)
51
52-- eigenvectors are columns
53eigTest method m msg = do
54 let (s,v) = method m
55 assertBool msg $ m <> v =~= v <> diag s
56
57bigmat = m + trans m where m = randomMatrix 18 (1000,1000)
58bigmatc = mc + conjTrans mc where mc = toComplexM(m,m)
59 m = randomMatrix 19 (1000,1000)
60
61--------------------------------------------------------------------
62
63invTest msg m = do
64 assertBool msg $ m <> inv m =~= ident (rows m)
65
66invComplexTest msg m = do
67 assertBool msg $ m <> invC m =~= ident (rows m)
68
69invC m = linearSolveC m (ident (rows m))
70
71--identC n = toComplexM(ident n, (0::Double) <>ident n)
72
73--------------------------------------------------------------------
74
75pinvTest f msg m = do
76 assertBool msg $ m <> f m <> m =~= m
77
78pinvC m = linearSolveLSC m (ident (rows m))
79
80pinvSVDR m = linearSolveSVDR Nothing m (ident (rows m))
81
82pinvSVDC m = linearSolveSVDC Nothing m (ident (rows m))
83
84--------------------------------------------------------------------
85
86
87tests = TestList [
88 TestCase $ besselTest
89 , TestCase $ exponentialTest
90 , TestCase $ invTest "inv 100x100" (randomMatrix 18 (100,100))
91 , TestCase $ invComplexTest "complex inv 100x100" (randomMatrixC 18 (100,100))
92 , TestCase $ pinvTest (pinvTolg 1) "pinvg 100x50" (randomMatrix 18 (100,50))
93 , TestCase $ pinvTest pinv "pinv 100x50" (randomMatrix 18 (100,50))
94 , TestCase $ pinvTest pinv "pinv 50x100" (randomMatrix 18 (50,100))
95 , TestCase $ pinvTest pinvSVDR "pinvSVDR 100x50" (randomMatrix 18 (100,50))
96 , TestCase $ pinvTest pinvSVDR "pinvSVDR 50x100" (randomMatrix 18 (50,100))
97 , TestCase $ pinvTest pinvC "pinvC 100x50" (randomMatrixC 18 (100,50))
98 , TestCase $ pinvTest pinvC "pinvC 50x100" (randomMatrixC 18 (50,100))
99 , TestCase $ pinvTest pinvSVDC "pinvSVDC 100x50" (randomMatrixC 18 (100,50))
100 , TestCase $ pinvTest pinvSVDC "pinvSVDC 50x100" (randomMatrixC 18 (50,100))
101 , TestCase $ eigTest eigC mcu "eigC"
102 , TestCase $ eigTest eigR mcur "eigR"
103 , TestCase $ eigTest eigS (mcur+trans mcur) "eigS"
104 , TestCase $ eigTest eigSg (mcur+trans mcur) "eigSg"
105 , TestCase $ eigTest eigH (mcu+ (conjTrans) mcu) "eigH"
106 , TestCase $ eigTest eigHg (mcu+ (conjTrans) mcu) "eigHg"
107 , TestCase $ fullsvdTest svdg' ms "GSL svd small"
108 , TestCase $ fullsvdTest svdR ms "fullsvdR small"
109 , TestCase $ fullsvdTest svdR (trans ms) "fullsvdR small"
110 , TestCase $ fullsvdTest svdR ms' "fullsvdR"
111 , TestCase $ fullsvdTest svdR (trans ms') "fullsvdR"
112 , TestCase $ fullsvdTest full_svd_Rd ms' "fullsvdRd"
113 , TestCase $ fullsvdTest full_svd_Rd (trans ms') "fullsvdRd"
114 , TestCase $ fullsvdTest svdC ms'' "fullsvdC"
115 , TestCase $ fullsvdTest svdC (trans ms'') "fullsvdC"
116 , TestCase $ eigTest eigS bigmat "big eigS"
117 , TestCase $ eigTest eigH bigmatc "big eigH"
118 , TestCase $ eigTest eigR bigmat "big eigR"
119 ]
120
121main = runTestTT tests
diff --git a/examples/tests.hs b/examples/tests.hs
index 9de9339..5acfa9d 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -20,7 +20,7 @@ import Test.HUnit hiding ((~:))
20import Complex 20import Complex
21import LinearAlgebra.Algorithms 21import LinearAlgebra.Algorithms
22import GSL.Matrix 22import GSL.Matrix
23import Data.Packed.Instances hiding ((<>)) 23import GSL.Compat hiding ((<>))
24 24
25dist :: (Normed t, Num t) => t -> t -> Double 25dist :: (Normed t, Num t) => t -> t -> Double
26dist a b = norm (a-b) 26dist a b = norm (a-b)
diff --git a/lib/Data/Packed/Instances.hs b/lib/Data/Packed/Instances.hs
deleted file mode 100644
index 4478469..0000000
--- a/lib/Data/Packed/Instances.hs
+++ /dev/null
@@ -1,391 +0,0 @@
1{-# OPTIONS_GHC -fglasgow-exts #-}
2-----------------------------------------------------------------------------
3{- |
4Module : Data.Packed.Instances
5Copyright : (c) Alberto Ruiz 2006
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : uses -fffi and -fglasgow-exts
11
12Creates reasonable numeric instances for Vectors and Matrices. In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.
13
14-}
15-----------------------------------------------------------------------------
16
17module Data.Packed.Instances(
18 Contractible(..)
19) where
20
21import Data.Packed.Internal
22import Data.Packed.Vector
23import Data.Packed.Matrix
24import GSL.Vector
25import GSL.Matrix
26import LinearAlgebra.Algorithms
27import Complex
28
29instance (Eq a, Field a) => Eq (Vector a) where
30 a == b = dim a == dim b && toList a == toList b
31
32instance (Num a, Field a) => Num (Vector a) where
33 (+) = add
34 (-) = sub
35 (*) = mul
36 signum = liftVector signum
37 abs = liftVector abs
38 fromInteger = fromList . return . fromInteger
39
40instance (Eq a, Field a) => Eq (Matrix a) where
41 a == b = rows a == rows b && cols a == cols b && cdat a == cdat b && fdat a == fdat b
42
43instance (Num a, Field a) => Num (Matrix a) where
44 (+) = liftMatrix2 add
45 (-) = liftMatrix2 sub
46 (*) = liftMatrix2 mul
47 signum = liftMatrix signum
48 abs = liftMatrix abs
49 fromInteger = (1><1) . return . fromInteger
50
51---------------------------------------------------
52
53adaptScalar f1 f2 f3 x y
54 | dim x == 1 = f1 (x@>0) y
55 | dim y == 1 = f3 x (y@>0)
56 | otherwise = f2 x y
57
58{-
59subvv = vectorZip 4
60subvc v c = addConstant (-c) v
61subcv c v = addConstant c (scale (-1) v)
62
63mul = vectorZip 1
64
65instance Num (Vector Double) where
66 (+) = adaptScalar addConstant add (flip addConstant)
67 (-) = adaptScalar subcv subvv subvc
68 (*) = adaptScalar scale mul (flip scale)
69 abs = vectorMap 3
70 signum = vectorMap 15
71 fromInteger n = fromList [fromInteger n]
72
73----------------------------------------------------
74
75--addConstantC a = gmap (+a)
76--subCvv u v = u `add` scale (-1) v
77subCvv = vectorZipComplex 4 -- faster?
78subCvc v c = addConstantC (-c) v
79subCcv c v = addConstantC c (scale (-1) v)
80
81
82instance Num (Vector (Complex Double)) where
83 (+) = adaptScalar addConstantC add (flip addConstantC)
84 (-) = adaptScalar subCcv subCvv subCvc
85 (*) = adaptScalar scale (vectorZipComplex 1) (flip scale)
86 abs = gmap abs
87 signum = gmap signum
88 fromInteger n = fromList [fromInteger n]
89
90
91-- | adapts a function on two vectors to work on all the elements of two matrices
92liftMatrix2' :: (Vector a -> Vector b -> Vector c) -> Matrix a -> Matrix b -> Matrix c
93liftMatrix2' f m1@(M r1 c1 _) m2@(M r2 c2 _)
94 | sameShape m1 m2 || r1*c1==1 || r2*c2==1
95 = reshape (max c1 c2) $ f (flatten m1) (flatten m2)
96 | otherwise = error "inconsistent matrix dimensions"
97
98---------------------------------------------------
99
100instance (Eq a, Field a) => Eq (Matrix a) where
101 a == b = rows a == rows b && cdat a == cdat b
102
103instance Num (Matrix Double) where
104 (+) = liftMatrix2' (+)
105 (-) = liftMatrix2' (-)
106 (*) = liftMatrix2' (*)
107 abs = liftMatrix abs
108 signum = liftMatrix signum
109 fromInteger n = fromLists [[fromInteger n]]
110
111----------------------------------------------------
112
113instance Num (Matrix (Complex Double)) where
114 (+) = liftMatrix2' (+)
115 (-) = liftMatrix2' (-)
116 (*) = liftMatrix2' (*)
117 abs = liftMatrix abs
118 signum = liftMatrix signum
119 fromInteger n = fromLists [[fromInteger n]]
120
121------------------------------------------------------
122
123instance Fractional (Vector Double) where
124 fromRational n = fromList [fromRational n]
125 (/) = adaptScalar f (vectorZip 2) g where
126 r `f` v = vectorZip 2 (constant r (dim v)) v
127 v `g` r = scale (recip r) v
128
129-------------------------------------------------------
130
131instance Fractional (Vector (Complex Double)) where
132 fromRational n = fromList [fromRational n]
133 (/) = adaptScalar f (vectorZipComplex 2) g where
134 r `f` v = gmap ((*r).recip) v
135 v `g` r = gmap (/r) v
136
137------------------------------------------------------
138
139instance Fractional (Matrix Double) where
140 fromRational n = fromLists [[fromRational n]]
141 (/) = liftMatrix2' (/)
142
143-------------------------------------------------------
144
145instance Fractional (Matrix (Complex Double)) where
146 fromRational n = fromLists [[fromRational n]]
147 (/) = liftMatrix2' (/)
148
149---------------------------------------------------------
150
151instance Floating (Vector Double) where
152 sin = vectorMap 0
153 cos = vectorMap 1
154 tan = vectorMap 2
155 asin = vectorMap 4
156 acos = vectorMap 5
157 atan = vectorMap 6
158 sinh = vectorMap 7
159 cosh = vectorMap 8
160 tanh = vectorMap 9
161 asinh = vectorMap 10
162 acosh = vectorMap 11
163 atanh = vectorMap 12
164 exp = vectorMap 13
165 log = vectorMap 14
166 sqrt = vectorMap 16
167 (**) = adaptScalar f (vectorZip 5) g where f s v = constant s (dim v) ** v
168 g v s = v ** constant s (dim v)
169 pi = fromList [pi]
170
171-----------------------------------------------------------
172
173instance Floating (Matrix Double) where
174 sin = liftMatrix sin
175 cos = liftMatrix cos
176 tan = liftMatrix tan
177 asin = liftMatrix asin
178 acos = liftMatrix acos
179 atan = liftMatrix atan
180 sinh = liftMatrix sinh
181 cosh = liftMatrix cosh
182 tanh = liftMatrix tanh
183 asinh = liftMatrix asinh
184 acosh = liftMatrix acosh
185 atanh = liftMatrix atanh
186 exp = liftMatrix exp
187 log = liftMatrix log
188 sqrt = liftMatrix sqrt
189 (**) = liftMatrix2 (**)
190 pi = fromLists [[pi]]
191
192-------------------------------------------------------------
193
194instance Floating (Vector (Complex Double)) where
195 sin = vectorMapComplex 0
196 cos = vectorMapComplex 1
197 tan = vectorMapComplex 2
198 asin = vectorMapComplex 4
199 acos = vectorMapComplex 5
200 atan = vectorMapComplex 6
201 sinh = vectorMapComplex 7
202 cosh = vectorMapComplex 8
203 tanh = vectorMapComplex 9
204 asinh = vectorMapComplex 10
205 acosh = vectorMapComplex 11
206 atanh = vectorMapComplex 12
207 exp = vectorMapComplex 13
208 log = vectorMapComplex 14
209 sqrt = vectorMapComplex 16
210 (**) = adaptScalar f (vectorZipComplex 5) g where f s v = constantC s (dim v) ** v
211 g v s = v ** constantC s (dim v)
212 pi = fromList [pi]
213
214---------------------------------------------------------------
215
216instance Floating (Matrix (Complex Double)) where
217 sin = liftMatrix sin
218 cos = liftMatrix cos
219 tan = liftMatrix tan
220 asin = liftMatrix asin
221 acos = liftMatrix acos
222 atan = liftMatrix atan
223 sinh = liftMatrix sinh
224 cosh = liftMatrix cosh
225 tanh = liftMatrix tanh
226 asinh = liftMatrix asinh
227 acosh = liftMatrix acosh
228 atanh = liftMatrix atanh
229 exp = liftMatrix exp
230 log = liftMatrix log
231 (**) = liftMatrix2 (**)
232 sqrt = liftMatrix sqrt
233 pi = fromLists [[pi]]
234
235---------------------------------------------------------------
236-}
237
238class Contractible a b c | a b -> c where
239 infixl 7 <>
240{- | An overloaded operator for matrix products, matrix-vector and vector-matrix products, dot products and scaling of vectors and matrices. Type consistency is statically checked. Alternatively, you can use the specific functions described below, but using this operator you can automatically combine real and complex objects.
241
242@v = 'fromList' [1,2,3] :: Vector Double
243cv = 'fromList' [1+'i',2]
244m = 'fromLists' [[1,2,3],
245 [4,5,7]] :: Matrix Double
246cm = 'fromLists' [[ 1, 2],
247 [3+'i',7*'i'],
248 [ 'i', 1]]
249\
250\> m \<\> v
25114. 35.
252\
253\> cv \<\> m
2549.+1.i 12.+2.i 17.+3.i
255\
256\> m \<\> cm
257 7.+5.i 5.+14.i
25819.+12.i 15.+35.i
259\
260\> v \<\> 'i'
2611.i 2.i 3.i
262\
263\> v \<\> v
26414.0
265\
266\> cv \<\> cv
2674.0 :+ 2.0@
268
269-}
270 (<>) :: a -> b -> c
271
272
273instance Contractible Double Double Double where
274 (<>) = (*)
275
276instance Contractible Double (Complex Double) (Complex Double) where
277 a <> b = (a:+0) * b
278
279instance Contractible (Complex Double) Double (Complex Double) where
280 a <> b = a * (b:+0)
281
282instance Contractible (Complex Double) (Complex Double) (Complex Double) where
283 (<>) = (*)
284
285--------------------------------- matrix matrix
286
287instance Contractible (Matrix Double) (Matrix Double) (Matrix Double) where
288 (<>) = mXm
289
290instance Contractible (Matrix (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where
291 (<>) = mXm
292
293instance Contractible (Matrix (Complex Double)) (Matrix Double) (Matrix (Complex Double)) where
294 c <> r = c <> liftMatrix comp r
295
296instance Contractible (Matrix Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where
297 r <> c = liftMatrix comp r <> c
298
299--------------------------------- (Matrix Double) (Vector Double)
300
301instance Contractible (Matrix Double) (Vector Double) (Vector Double) where
302 (<>) = mXv
303
304instance Contractible (Matrix (Complex Double)) (Vector (Complex Double)) (Vector (Complex Double)) where
305 (<>) = mXv
306
307instance Contractible (Matrix (Complex Double)) (Vector Double) (Vector (Complex Double)) where
308 m <> v = m <> comp v
309
310instance Contractible (Matrix Double) (Vector (Complex Double)) (Vector (Complex Double)) where
311 m <> v = liftMatrix comp m <> v
312
313--------------------------------- (Vector Double) (Matrix Double)
314
315instance Contractible (Vector Double) (Matrix Double) (Vector Double) where
316 (<>) = vXm
317
318instance Contractible (Vector (Complex Double)) (Matrix (Complex Double)) (Vector (Complex Double)) where
319 (<>) = vXm
320
321instance Contractible (Vector (Complex Double)) (Matrix Double) (Vector (Complex Double)) where
322 v <> m = v <> liftMatrix comp m
323
324instance Contractible (Vector Double) (Matrix (Complex Double)) (Vector (Complex Double)) where
325 v <> m = comp v <> m
326
327--------------------------------- dot product
328
329instance Contractible (Vector Double) (Vector Double) Double where
330 (<>) = dot
331
332instance Contractible (Vector (Complex Double)) (Vector (Complex Double)) (Complex Double) where
333 (<>) = dot
334
335instance Contractible (Vector Double) (Vector (Complex Double)) (Complex Double) where
336 a <> b = comp a <> b
337
338instance Contractible (Vector (Complex Double)) (Vector Double) (Complex Double) where
339 (<>) = flip (<>)
340
341--------------------------------- scaling vectors
342
343instance Contractible Double (Vector Double) (Vector Double) where
344 (<>) = scale
345
346instance Contractible (Vector Double) Double (Vector Double) where
347 (<>) = flip (<>)
348
349instance Contractible (Complex Double) (Vector (Complex Double)) (Vector (Complex Double)) where
350 (<>) = scale
351
352instance Contractible (Vector (Complex Double)) (Complex Double) (Vector (Complex Double)) where
353 (<>) = flip (<>)
354
355instance Contractible Double (Vector (Complex Double)) (Vector (Complex Double)) where
356 a <> v = (a:+0) <> v
357
358instance Contractible (Vector (Complex Double)) Double (Vector (Complex Double)) where
359 (<>) = flip (<>)
360
361instance Contractible (Complex Double) (Vector Double) (Vector (Complex Double)) where
362 a <> v = a <> comp v
363
364instance Contractible (Vector Double) (Complex Double) (Vector (Complex Double)) where
365 (<>) = flip (<>)
366
367--------------------------------- scaling matrices
368
369instance Contractible Double (Matrix Double) (Matrix Double) where
370 (<>) a = liftMatrix (a <>)
371
372instance Contractible (Matrix Double) Double (Matrix Double) where
373 (<>) = flip (<>)
374
375instance Contractible (Complex Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where
376 (<>) a = liftMatrix (a <>)
377
378instance Contractible (Matrix (Complex Double)) (Complex Double) (Matrix (Complex Double)) where
379 (<>) = flip (<>)
380
381instance Contractible Double (Matrix (Complex Double)) (Matrix (Complex Double)) where
382 a <> m = (a:+0) <> m
383
384instance Contractible (Matrix (Complex Double)) Double (Matrix (Complex Double)) where
385 (<>) = flip (<>)
386
387instance Contractible (Complex Double) (Matrix Double) (Matrix (Complex Double)) where
388 a <> m = a <> liftMatrix comp m
389
390instance Contractible (Matrix Double) (Complex Double) (Matrix (Complex Double)) where
391 (<>) = flip (<>)
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index f1addf4..ab93577 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -129,5 +129,8 @@ constant x n | isReal id x = scast $ constantR (scast x) n
129 | isComp id x = scast $ constantC (scast x) n 129 | isComp id x = scast $ constantC (scast x) n
130 | otherwise = constantG x n 130 | otherwise = constantG x n
131 131
132liftVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b
132liftVector f = fromList . map f . toList 133liftVector f = fromList . map f . toList
134
135liftVector2 :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c
133liftVector2 f u v = fromList $ zipWith f (toList u) (toList v) 136liftVector2 f u v = fromList $ zipWith f (toList u) (toList v)
diff --git a/lib/GSL.hs b/lib/GSL.hs
index 8e033c3..ac3a2e8 100644
--- a/lib/GSL.hs
+++ b/lib/GSL.hs
@@ -17,7 +17,6 @@ module GSL (
17module Data.Packed.Vector, 17module Data.Packed.Vector,
18module Data.Packed.Matrix, 18module Data.Packed.Matrix,
19module Data.Packed.Tensor, 19module Data.Packed.Tensor,
20module Data.Packed.Instances,
21module LinearAlgebra.Algorithms, 20module LinearAlgebra.Algorithms,
22module LAPACK, 21module LAPACK,
23module GSL.Integration, 22module GSL.Integration,
@@ -26,14 +25,14 @@ module GSL.Special,
26module GSL.Fourier, 25module GSL.Fourier,
27module GSL.Polynomials, 26module GSL.Polynomials,
28module GSL.Minimization, 27module GSL.Minimization,
29module Data.Packed.Plot 28module Data.Packed.Plot,
29module GSL.Compat
30 30
31) where 31) where
32 32
33import Data.Packed.Vector 33import Data.Packed.Vector
34import Data.Packed.Matrix 34import Data.Packed.Matrix
35import Data.Packed.Tensor 35import Data.Packed.Tensor
36import Data.Packed.Instances
37import LinearAlgebra.Algorithms 36import LinearAlgebra.Algorithms
38import LAPACK 37import LAPACK
39import GSL.Integration 38import GSL.Integration
@@ -43,3 +42,4 @@ import GSL.Fourier
43import GSL.Polynomials 42import GSL.Polynomials
44import GSL.Minimization 43import GSL.Minimization
45import Data.Packed.Plot 44import Data.Packed.Plot
45import GSL.Compat
diff --git a/lib/GSL/Compat.hs b/lib/GSL/Compat.hs
new file mode 100644
index 0000000..6a94191
--- /dev/null
+++ b/lib/GSL/Compat.hs
@@ -0,0 +1,370 @@
1{-# OPTIONS_GHC -fglasgow-exts #-}
2-----------------------------------------------------------------------------
3{- |
4Module : GSL.Compat
5Copyright : (c) Alberto Ruiz 2006
6License : GPL-style
7
8Maintainer : Alberto Ruiz (aruiz at um dot es)
9Stability : provisional
10Portability : uses -fffi and -fglasgow-exts
11
12Creates reasonable numeric instances for Vectors and Matrices. In the context of the standard numeric operators, one-component vectors and matrices automatically expand to match the dimensions of the other operand.
13
14-}
15-----------------------------------------------------------------------------
16
17module GSL.Compat(
18 Mul,(<>), fromFile, readMatrix, size, dispR, dispC, format, gmap
19) where
20
21import Data.Packed.Internal hiding (dsp)
22import Data.Packed.Vector
23import Data.Packed.Matrix
24import GSL.Vector
25import GSL.Matrix
26import LinearAlgebra.Algorithms
27import Complex
28import Numeric(showGFloat)
29import Data.List(transpose,intersperse)
30
31
32adaptScalar f1 f2 f3 x y
33 | dim x == 1 = f1 (x@>0) y
34 | dim y == 1 = f3 x (y@>0)
35 | otherwise = f2 x y
36
37instance (Eq a, Field a) => Eq (Vector a) where
38 a == b = dim a == dim b && toList a == toList b
39
40instance (Num a, Field a) => Num (Vector a) where
41 (+) = adaptScalar addConstant add (flip addConstant)
42 negate = scale (-1)
43 (*) = adaptScalar scale mul (flip scale)
44 signum = liftVector signum
45 abs = liftVector abs
46 fromInteger = fromList . return . fromInteger
47
48instance (Eq a, Field a) => Eq (Matrix a) where
49 a == b = rows a == rows b && cols a == cols b && cdat a == cdat b && fdat a == fdat b
50
51instance (Num a, Field a) => Num (Matrix a) where
52 (+) = liftMatrix2 (+)
53 negate = liftMatrix negate
54 (*) = liftMatrix2 (*)
55 signum = liftMatrix signum
56 abs = liftMatrix abs
57 fromInteger = (1><1) . return . fromInteger
58
59---------------------------------------------------
60
61instance Fractional (Vector Double) where
62 fromRational n = fromList [fromRational n]
63 (/) = adaptScalar f (vectorZipR Div) g where
64 r `f` v = vectorMapValR Recip r v
65 v `g` r = scale (recip r) v
66
67-------------------------------------------------------
68
69instance Fractional (Vector (Complex Double)) where
70 fromRational n = fromList [fromRational n]
71 (/) = adaptScalar f (vectorZipC Div) g where
72 r `f` v = vectorMapValC Recip r v
73 v `g` r = scale (recip r) v
74
75------------------------------------------------------
76
77instance Fractional (Matrix Double) where
78 fromRational n = (1><1) [fromRational n]
79 (/) = liftMatrix2 (/)
80
81-------------------------------------------------------
82
83instance Fractional (Matrix (Complex Double)) where
84 fromRational n = (1><1) [fromRational n]
85 (/) = liftMatrix2 (/)
86
87---------------------------------------------------------
88
89instance Floating (Vector Double) where
90 sin = vectorMapR Sin
91 cos = vectorMapR Cos
92 tan = vectorMapR Tan
93 asin = vectorMapR ASin
94 acos = vectorMapR ACos
95 atan = vectorMapR ATan
96 sinh = vectorMapR Sinh
97 cosh = vectorMapR Cosh
98 tanh = vectorMapR Tanh
99 asinh = vectorMapR ASinh
100 acosh = vectorMapR ACosh
101 atanh = vectorMapR ATanh
102 exp = vectorMapR Exp
103 log = vectorMapR Log
104 sqrt = vectorMapR Sqrt
105 (**) = adaptScalar (vectorMapValR PowSV) (vectorZipR Pow) (flip (vectorMapValR PowVS))
106 pi = fromList [pi]
107
108-----------------------------------------------------------
109
110instance Floating (Matrix Double) where
111 sin = liftMatrix sin
112 cos = liftMatrix cos
113 tan = liftMatrix tan
114 asin = liftMatrix asin
115 acos = liftMatrix acos
116 atan = liftMatrix atan
117 sinh = liftMatrix sinh
118 cosh = liftMatrix cosh
119 tanh = liftMatrix tanh
120 asinh = liftMatrix asinh
121 acosh = liftMatrix acosh
122 atanh = liftMatrix atanh
123 exp = liftMatrix exp
124 log = liftMatrix log
125 (**) = liftMatrix2 (**)
126 sqrt = liftMatrix sqrt
127 pi = (1><1) [pi]
128-------------------------------------------------------------
129
130instance Floating (Vector (Complex Double)) where
131 sin = vectorMapC Sin
132 cos = vectorMapC Cos
133 tan = vectorMapC Tan
134 asin = vectorMapC ASin
135 acos = vectorMapC ACos
136 atan = vectorMapC ATan
137 sinh = vectorMapC Sinh
138 cosh = vectorMapC Cosh
139 tanh = vectorMapC Tanh
140 asinh = vectorMapC ASinh
141 acosh = vectorMapC ACosh
142 atanh = vectorMapC ATanh
143 exp = vectorMapC Exp
144 log = vectorMapC Log
145 sqrt = vectorMapC Sqrt
146 (**) = adaptScalar (vectorMapValC PowSV) (vectorZipC Pow) (flip (vectorMapValC PowVS))
147 pi = fromList [pi]
148
149---------------------------------------------------------------
150
151instance Floating (Matrix (Complex Double)) where
152 sin = liftMatrix sin
153 cos = liftMatrix cos
154 tan = liftMatrix tan
155 asin = liftMatrix asin
156 acos = liftMatrix acos
157 atan = liftMatrix atan
158 sinh = liftMatrix sinh
159 cosh = liftMatrix cosh
160 tanh = liftMatrix tanh
161 asinh = liftMatrix asinh
162 acosh = liftMatrix acosh
163 atanh = liftMatrix atanh
164 exp = liftMatrix exp
165 log = liftMatrix log
166 (**) = liftMatrix2 (**)
167 sqrt = liftMatrix sqrt
168 pi = (1><1) [pi]
169
170---------------------------------------------------------------
171
172
173class Mul a b c | a b -> c where
174 infixl 7 <>
175{- | An overloaded operator for matrix products, matrix-vector and vector-matrix products, dot products and scaling of vectors and matrices. Type consistency is statically checked. Alternatively, you can use the specific functions described below, but using this operator you can automatically combine real and complex objects.
176
177@v = 'fromList' [1,2,3] :: Vector Double
178cv = 'fromList' [1+'i',2]
179m = 'fromLists' [[1,2,3],
180 [4,5,7]] :: Matrix Double
181cm = 'fromLists' [[ 1, 2],
182 [3+'i',7*'i'],
183 [ 'i', 1]]
184\
185\> m \<\> v
18614. 35.
187\
188\> cv \<\> m
1899.+1.i 12.+2.i 17.+3.i
190\
191\> m \<\> cm
192 7.+5.i 5.+14.i
19319.+12.i 15.+35.i
194\
195\> v \<\> 'i'
1961.i 2.i 3.i
197\
198\> v \<\> v
19914.0
200\
201\> cv \<\> cv
2024.0 :+ 2.0@
203
204-}
205 (<>) :: a -> b -> c
206
207
208instance Mul Double Double Double where
209 (<>) = (*)
210
211instance Mul Double (Complex Double) (Complex Double) where
212 a <> b = (a:+0) * b
213
214instance Mul (Complex Double) Double (Complex Double) where
215 a <> b = a * (b:+0)
216
217instance Mul (Complex Double) (Complex Double) (Complex Double) where
218 (<>) = (*)
219
220--------------------------------- matrix matrix
221
222instance Mul (Matrix Double) (Matrix Double) (Matrix Double) where
223 (<>) = mXm
224
225instance Mul (Matrix (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where
226 (<>) = mXm
227
228instance Mul (Matrix (Complex Double)) (Matrix Double) (Matrix (Complex Double)) where
229 c <> r = c <> liftMatrix comp r
230
231instance Mul (Matrix Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where
232 r <> c = liftMatrix comp r <> c
233
234--------------------------------- (Matrix Double) (Vector Double)
235
236instance Mul (Matrix Double) (Vector Double) (Vector Double) where
237 (<>) = mXv
238
239instance Mul (Matrix (Complex Double)) (Vector (Complex Double)) (Vector (Complex Double)) where
240 (<>) = mXv
241
242instance Mul (Matrix (Complex Double)) (Vector Double) (Vector (Complex Double)) where
243 m <> v = m <> comp v
244
245instance Mul (Matrix Double) (Vector (Complex Double)) (Vector (Complex Double)) where
246 m <> v = liftMatrix comp m <> v
247
248--------------------------------- (Vector Double) (Matrix Double)
249
250instance Mul (Vector Double) (Matrix Double) (Vector Double) where
251 (<>) = vXm
252
253instance Mul (Vector (Complex Double)) (Matrix (Complex Double)) (Vector (Complex Double)) where
254 (<>) = vXm
255
256instance Mul (Vector (Complex Double)) (Matrix Double) (Vector (Complex Double)) where
257 v <> m = v <> liftMatrix comp m
258
259instance Mul (Vector Double) (Matrix (Complex Double)) (Vector (Complex Double)) where
260 v <> m = comp v <> m
261
262--------------------------------- dot product
263
264instance Mul (Vector Double) (Vector Double) Double where
265 (<>) = dot
266
267instance Mul (Vector (Complex Double)) (Vector (Complex Double)) (Complex Double) where
268 (<>) = dot
269
270instance Mul (Vector Double) (Vector (Complex Double)) (Complex Double) where
271 a <> b = comp a <> b
272
273instance Mul (Vector (Complex Double)) (Vector Double) (Complex Double) where
274 (<>) = flip (<>)
275
276--------------------------------- scaling vectors
277
278instance Mul Double (Vector Double) (Vector Double) where
279 (<>) = scale
280
281instance Mul (Vector Double) Double (Vector Double) where
282 (<>) = flip (<>)
283
284instance Mul (Complex Double) (Vector (Complex Double)) (Vector (Complex Double)) where
285 (<>) = scale
286
287instance Mul (Vector (Complex Double)) (Complex Double) (Vector (Complex Double)) where
288 (<>) = flip (<>)
289
290instance Mul Double (Vector (Complex Double)) (Vector (Complex Double)) where
291 a <> v = (a:+0) <> v
292
293instance Mul (Vector (Complex Double)) Double (Vector (Complex Double)) where
294 (<>) = flip (<>)
295
296instance Mul (Complex Double) (Vector Double) (Vector (Complex Double)) where
297 a <> v = a <> comp v
298
299instance Mul (Vector Double) (Complex Double) (Vector (Complex Double)) where
300 (<>) = flip (<>)
301
302--------------------------------- scaling matrices
303
304instance Mul Double (Matrix Double) (Matrix Double) where
305 (<>) a = liftMatrix (a <>)
306
307instance Mul (Matrix Double) Double (Matrix Double) where
308 (<>) = flip (<>)
309
310instance Mul (Complex Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where
311 (<>) a = liftMatrix (a <>)
312
313instance Mul (Matrix (Complex Double)) (Complex Double) (Matrix (Complex Double)) where
314 (<>) = flip (<>)
315
316instance Mul Double (Matrix (Complex Double)) (Matrix (Complex Double)) where
317 a <> m = (a:+0) <> m
318
319instance Mul (Matrix (Complex Double)) Double (Matrix (Complex Double)) where
320 (<>) = flip (<>)
321
322instance Mul (Complex Double) (Matrix Double) (Matrix (Complex Double)) where
323 a <> m = a <> liftMatrix comp m
324
325instance Mul (Matrix Double) (Complex Double) (Matrix (Complex Double)) where
326 (<>) = flip (<>)
327
328-----------------------------------------------------------------------------------
329
330size :: Vector a -> Int
331size = dim
332
333gmap f v = liftVector f v
334
335
336-- shows a Double with n digits after the decimal point
337shf :: (RealFloat a) => Int -> a -> String
338shf dec n | abs n < 1e-10 = "0."
339 | abs (n - (fromIntegral.round $ n)) < 1e-10 = show (round n) ++"."
340 | otherwise = showGFloat (Just dec) n ""
341-- shows a Complex Double as a pair, with n digits after the decimal point
342shfc n z@ (a:+b)
343 | magnitude z <1e-10 = "0."
344 | abs b < 1e-10 = shf n a
345 | abs a < 1e-10 = shf n b ++"i"
346 | b > 0 = shf n a ++"+"++shf n b ++"i"
347 | otherwise = shf n a ++shf n b ++"i"
348
349dsp :: String -> [[String]] -> String
350dsp sep as = unlines . map unwords' $ transpose mtp where
351 mt = transpose as
352 longs = map (maximum . map length) mt
353 mtp = zipWith (\a b -> map (pad a) b) longs mt
354 pad n str = replicate (n - length str) ' ' ++ str
355 unwords' = concat . intersperse sep
356
357format :: (Field t) => String -> (t -> String) -> Matrix t -> String
358format sep f m = dsp sep . map (map f) . toLists $ m
359
360disp m f = putStrLn $ "matrix ("++show (rows m) ++"x"++ show (cols m) ++")\n"++format " | " f m
361
362dispR :: Int -> Matrix Double -> IO ()
363dispR d m = disp m (shf d)
364
365dispC :: Int -> Matrix (Complex Double) -> IO ()
366dispC d m = disp m (shfc d)
367
368-- | creates a matrix from a table of numbers.
369readMatrix :: String -> Matrix Double
370readMatrix = fromLists . map (map read). map words . filter (not.null) . lines \ No newline at end of file
diff --git a/lib/GSL/Minimization.hs b/lib/GSL/Minimization.hs
index a59977e..aa89475 100644
--- a/lib/GSL/Minimization.hs
+++ b/lib/GSL/Minimization.hs
@@ -150,9 +150,11 @@ minimizeConjugateGradient istep minimpar tol maxit f df xi = unsafePerformIO $ d
150 df' = (fromList . df . toList) 150 df' = (fromList . df . toList)
151 fp <- mkVecfun (iv f') 151 fp <- mkVecfun (iv f')
152 dfp <- mkVecVecfun (aux_vTov df') 152 dfp <- mkVecVecfun (aux_vTov df')
153 print "entro"
153 rawpath <- createMIO maxit (n+2) 154 rawpath <- createMIO maxit (n+2)
154 (c_minimizeConjugateGradient fp dfp istep minimpar tol maxit // vec xiv) 155 (c_minimizeConjugateGradient fp dfp istep minimpar tol maxit // vec xiv)
155 "minimizeDerivV" [xiv] 156 "minimizeDerivV" [xiv]
157 print "salgo"
156 let it = round (rawpath @@> (maxit-1,0)) 158 let it = round (rawpath @@> (maxit-1,0))
157 path = takeRows it rawpath 159 path = takeRows it rawpath
158 sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path 160 sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path
@@ -186,12 +188,12 @@ foreign import ccall "wrapper"
186 188
187aux_vTov :: (Vector Double -> Vector Double) -> (Int -> Ptr Double -> Ptr Double -> IO()) 189aux_vTov :: (Vector Double -> Vector Double) -> (Int -> Ptr Double -> Ptr Double -> IO())
188aux_vTov f n p r = g where 190aux_vTov f n p r = g where
189 V {fptr = pr, ptr = p} = f x 191 V {fptr = pr, ptr = t} = f x
190 x = createV n copy "aux_vTov" [] 192 x = createV n copy "aux_vTov" []
191 copy n q = do 193 copy n q = do
192 copyArray q p n 194 copyArray q p n
193 return 0 195 return 0
194 g = withForeignPtr pr $ \_ -> copyArray r p n 196 g = withForeignPtr pr $ \_ -> copyArray r t n
195 197
196-------------------------------------------------------------------- 198--------------------------------------------------------------------
197 199
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-10 06:05:21 +0000