summaryrefslogtreecommitdiff
path: root/lib/GSLHaskell.hs
diff options
context:
space:
mode:
Diffstat (limited to 'lib/GSLHaskell.hs')
-rw-r--r--lib/GSLHaskell.hs327
1 files changed, 0 insertions, 327 deletions
diff --git a/lib/GSLHaskell.hs b/lib/GSLHaskell.hs
deleted file mode 100644
index 254a957..0000000
--- a/lib/GSLHaskell.hs
+++ /dev/null
@@ -1,327 +0,0 @@
1{-# OPTIONS_GHC -fglasgow-exts -fallow-undecidable-instances #-}
2-----------------------------------------------------------------------------
3{- |
4Module : GSLHaskell
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
12Old GSLHaskell interface.
13
14-}
15-----------------------------------------------------------------------------
16
17module GSLHaskell(
18 module Data.Packed.Vector,
19 module Data.Packed.Matrix,
20 module LinearAlgebra.Algorithms,
21 module LinearAlgebra.Linear,
22 module LAPACK,
23 module GSL.Integration,
24 module GSL.Differentiation,
25 module GSL.Fourier,
26 module GSL.Polynomials,
27 module GSL.Minimization,
28 module GSL.Matrix,
29 module GSL.Special,
30 module Graphics.Plot,
31 module Complex,
32 Mul,(<>), readMatrix, size, dispR, dispC, format, gmap, Joinable, (<|>),(<->),
33 fromArray2D, GSLHaskell.pnorm,
34) where
35
36import GSL.Integration
37import GSL.Differentiation
38import GSL.Fourier
39import GSL.Polynomials
40import GSL.Minimization
41import Graphics.Plot
42import Complex
43import GSL.Special(setErrorHandlerOff,
44 erf,
45 erf_Z,
46 bessel_J0_e,
47 exp_e10_e,
48 gamma)
49import Data.Packed.Vector
50import Data.Packed.Matrix
51import Data.Packed.Matrix hiding ((><))
52import GSL.Vector
53import qualified LinearAlgebra.Algorithms
54import LAPACK
55import GSL.Matrix
56import LinearAlgebra.Algorithms hiding (pnorm)
57import LinearAlgebra.Linear hiding (Mul,(<>))
58import Data.Packed.Internal.Matrix(multiply)
59import Complex
60import Numeric(showGFloat)
61import Data.List(transpose,intersperse)
62import Foreign(Storable)
63import Data.Array
64import LinearAlgebra.Instances
65
66
67class Mul a b c | a b -> c where
68 infixl 7 <>
69{- | 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.
70
71@v = 'fromList' [1,2,3] :: Vector Double
72cv = 'fromList' [1+'i',2]
73m = 'fromLists' [[1,2,3],
74 [4,5,7]] :: Matrix Double
75cm = 'fromLists' [[ 1, 2],
76 [3+'i',7*'i'],
77 [ 'i', 1]]
78\
79\> m \<\> v
8014. 35.
81\
82\> cv \<\> m
839.+1.i 12.+2.i 17.+3.i
84\
85\> m \<\> cm
86 7.+5.i 5.+14.i
8719.+12.i 15.+35.i
88\
89\> v \<\> 'i'
901.i 2.i 3.i
91\
92\> v \<\> v
9314.0
94\
95\> cv \<\> cv
964.0 :+ 2.0@
97
98-}
99 (<>) :: a -> b -> c
100
101
102instance Mul Double Double Double where
103 (<>) = (*)
104
105instance Mul Double (Complex Double) (Complex Double) where
106 a <> b = (a:+0) * b
107
108instance Mul (Complex Double) Double (Complex Double) where
109 a <> b = a * (b:+0)
110
111instance Mul (Complex Double) (Complex Double) (Complex Double) where
112 (<>) = (*)
113
114--------------------------------- matrix matrix
115
116instance Mul (Matrix Double) (Matrix Double) (Matrix Double) where
117 (<>) = multiply
118
119instance Mul (Matrix (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where
120 (<>) = multiply
121
122instance Mul (Matrix (Complex Double)) (Matrix Double) (Matrix (Complex Double)) where
123 c <> r = c <> liftMatrix comp r
124
125instance Mul (Matrix Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where
126 r <> c = liftMatrix comp r <> c
127
128--------------------------------- (Matrix Double) (Vector Double)
129
130instance Mul (Matrix Double) (Vector Double) (Vector Double) where
131 (<>) = mXv
132
133instance Mul (Matrix (Complex Double)) (Vector (Complex Double)) (Vector (Complex Double)) where
134 (<>) = mXv
135
136instance Mul (Matrix (Complex Double)) (Vector Double) (Vector (Complex Double)) where
137 m <> v = m <> comp v
138
139instance Mul (Matrix Double) (Vector (Complex Double)) (Vector (Complex Double)) where
140 m <> v = liftMatrix comp m <> v
141
142--------------------------------- (Vector Double) (Matrix Double)
143
144instance Mul (Vector Double) (Matrix Double) (Vector Double) where
145 (<>) = vXm
146
147instance Mul (Vector (Complex Double)) (Matrix (Complex Double)) (Vector (Complex Double)) where
148 (<>) = vXm
149
150instance Mul (Vector (Complex Double)) (Matrix Double) (Vector (Complex Double)) where
151 v <> m = v <> liftMatrix comp m
152
153instance Mul (Vector Double) (Matrix (Complex Double)) (Vector (Complex Double)) where
154 v <> m = comp v <> m
155
156--------------------------------- dot product
157
158instance Mul (Vector Double) (Vector Double) Double where
159 (<>) = dot
160
161instance Mul (Vector (Complex Double)) (Vector (Complex Double)) (Complex Double) where
162 (<>) = dot
163
164instance Mul (Vector Double) (Vector (Complex Double)) (Complex Double) where
165 a <> b = comp a <> b
166
167instance Mul (Vector (Complex Double)) (Vector Double) (Complex Double) where
168 (<>) = flip (<>)
169
170--------------------------------- scaling vectors
171
172instance Mul Double (Vector Double) (Vector Double) where
173 (<>) = scale
174
175instance Mul (Vector Double) Double (Vector Double) where
176 (<>) = flip (<>)
177
178instance Mul (Complex Double) (Vector (Complex Double)) (Vector (Complex Double)) where
179 (<>) = scale
180
181instance Mul (Vector (Complex Double)) (Complex Double) (Vector (Complex Double)) where
182 (<>) = flip (<>)
183
184instance Mul Double (Vector (Complex Double)) (Vector (Complex Double)) where
185 a <> v = (a:+0) <> v
186
187instance Mul (Vector (Complex Double)) Double (Vector (Complex Double)) where
188 (<>) = flip (<>)
189
190instance Mul (Complex Double) (Vector Double) (Vector (Complex Double)) where
191 a <> v = a <> comp v
192
193instance Mul (Vector Double) (Complex Double) (Vector (Complex Double)) where
194 (<>) = flip (<>)
195
196--------------------------------- scaling matrices
197
198instance Mul Double (Matrix Double) (Matrix Double) where
199 (<>) a = liftMatrix (a <>)
200
201instance Mul (Matrix Double) Double (Matrix Double) where
202 (<>) = flip (<>)
203
204instance Mul (Complex Double) (Matrix (Complex Double)) (Matrix (Complex Double)) where
205 (<>) a = liftMatrix (a <>)
206
207instance Mul (Matrix (Complex Double)) (Complex Double) (Matrix (Complex Double)) where
208 (<>) = flip (<>)
209
210instance Mul Double (Matrix (Complex Double)) (Matrix (Complex Double)) where
211 a <> m = (a:+0) <> m
212
213instance Mul (Matrix (Complex Double)) Double (Matrix (Complex Double)) where
214 (<>) = flip (<>)
215
216instance Mul (Complex Double) (Matrix Double) (Matrix (Complex Double)) where
217 a <> m = a <> liftMatrix comp m
218
219instance Mul (Matrix Double) (Complex Double) (Matrix (Complex Double)) where
220 (<>) = flip (<>)
221
222-----------------------------------------------------------------------------------
223
224size :: Vector a -> Int
225size = dim
226
227gmap :: (Storable a, Storable b) => (a->b) -> Vector a -> Vector b
228gmap f v = liftVector f v
229
230-- shows a Double with n digits after the decimal point
231shf :: (RealFloat a) => Int -> a -> String
232shf dec n | abs n < 1e-10 = "0."
233 | abs (n - (fromIntegral.round $ n)) < 1e-10 = show (round n) ++"."
234 | otherwise = showGFloat (Just dec) n ""
235-- shows a Complex Double as a pair, with n digits after the decimal point
236shfc n z@ (a:+b)
237 | magnitude z <1e-10 = "0."
238 | abs b < 1e-10 = shf n a
239 | abs a < 1e-10 = shf n b ++"i"
240 | b > 0 = shf n a ++"+"++shf n b ++"i"
241 | otherwise = shf n a ++shf n b ++"i"
242
243dsp :: String -> [[String]] -> String
244dsp sep as = unlines . map unwords' $ transpose mtp where
245 mt = transpose as
246 longs = map (maximum . map length) mt
247 mtp = zipWith (\a b -> map (pad a) b) longs mt
248 pad n str = replicate (n - length str) ' ' ++ str
249 unwords' = concat . intersperse sep
250
251format :: (Field t) => String -> (t -> String) -> Matrix t -> String
252format sep f m = dsp sep . map (map f) . toLists $ m
253
254disp m f = putStrLn $ "matrix ("++show (rows m) ++"x"++ show (cols m) ++")\n"++format " | " f m
255
256dispR :: Int -> Matrix Double -> IO ()
257dispR d m = disp m (shf d)
258
259dispC :: Int -> Matrix (Complex Double) -> IO ()
260dispC d m = disp m (shfc d)
261
262-- | creates a matrix from a table of numbers.
263readMatrix :: String -> Matrix Double
264readMatrix = fromLists . map (map read). map words . filter (not.null) . lines
265
266-------------------------------------------------------------
267
268class Joinable a b c | a b -> c where
269 joinH :: a -> b -> c
270 joinV :: a -> b -> c
271
272instance Joinable (Matrix Double) (Vector Double) (Matrix Double) where
273 joinH m v = fromBlocks [[m,reshape 1 v]]
274 joinV m v = fromBlocks [[m],[reshape (size v) v]]
275
276instance Joinable (Vector Double) (Matrix Double) (Matrix Double) where
277 joinH v m = fromBlocks [[reshape 1 v,m]]
278 joinV v m = fromBlocks [[reshape (size v) v],[m]]
279
280instance Joinable (Matrix Double) (Matrix Double) (Matrix Double) where
281 joinH m1 m2 = fromBlocks [[m1,m2]]
282 joinV m1 m2 = fromBlocks [[m1],[m2]]
283
284instance Joinable (Matrix (Complex Double)) (Vector (Complex Double)) (Matrix (Complex Double)) where
285 joinH m v = fromBlocks [[m,reshape 1 v]]
286 joinV m v = fromBlocks [[m],[reshape (size v) v]]
287
288instance Joinable (Vector (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where
289 joinH v m = fromBlocks [[reshape 1 v,m]]
290 joinV v m = fromBlocks [[reshape (size v) v],[m]]
291
292instance Joinable (Matrix (Complex Double)) (Matrix (Complex Double)) (Matrix (Complex Double)) where
293 joinH m1 m2 = fromBlocks [[m1,m2]]
294 joinV m1 m2 = fromBlocks [[m1],[m2]]
295
296infixl 3 <|>, <->
297
298{- | Horizontal concatenation of matrices and vectors:
299
300@\> 'ident' 3 \<-\> i\<\>'ident' 3 \<|\> 'fromList' [1..6]
301 1. 0. 0. 1.
302 0. 1. 0. 2.
303 0. 0. 1. 3.
3041.i 0. 0. 4.
305 0. 1.i 0. 5.
306 0. 0. 1.i 6.@
307-}
308(<|>) :: (Joinable a b c) => a -> b -> c
309a <|> b = joinH a b
310
311-- | Vertical concatenation of matrices and vectors.
312(<->) :: (Joinable a b c) => a -> b -> c
313a <-> b = joinV a b
314
315----------------------------------------------------------
316
317fromArray2D :: (Field e) => Array (Int, Int) e -> Matrix e
318fromArray2D m = (r><c) (elems m)
319 where ((r0,c0),(r1,c1)) = bounds m
320 r = r1-r0+1
321 c = c1-c0+1
322
323
324pnorm :: (Normed t1, Num t) => t -> t1 -> Double
325pnorm 0 = LinearAlgebra.Algorithms.pnorm Infinity
326pnorm 1 = LinearAlgebra.Algorithms.pnorm PNorm1
327pnorm 2 = LinearAlgebra.Algorithms.pnorm PNorm2 \ No newline at end of file