summaryrefslogtreecommitdiff
path: root/packages/hmatrix/src/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'packages/hmatrix/src/Numeric')
-rw-r--r--packages/hmatrix/src/Numeric/GSL.hs43
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Differentiation.hs85
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Fitting.hs180
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Fourier.hs44
-rw-r--r--packages/hmatrix/src/Numeric/GSL/IO.hs42
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Integration.hs246
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Internal.hs126
-rw-r--r--packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs135
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Minimization.hs222
-rw-r--r--packages/hmatrix/src/Numeric/GSL/ODE.hs140
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Polynomials.hs57
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Random.hs84
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Root.hs199
-rw-r--r--packages/hmatrix/src/Numeric/GSL/Vector.hs107
-rw-r--r--packages/hmatrix/src/Numeric/GSL/gsl-aux.c945
-rw-r--r--packages/hmatrix/src/Numeric/GSL/gsl-ode.c182
16 files changed, 0 insertions, 2837 deletions
diff --git a/packages/hmatrix/src/Numeric/GSL.hs b/packages/hmatrix/src/Numeric/GSL.hs
deleted file mode 100644
index 61b8d7b..0000000
--- a/packages/hmatrix/src/Numeric/GSL.hs
+++ /dev/null
@@ -1,43 +0,0 @@
1{- |
2
3Module : Numeric.GSL
4Copyright : (c) Alberto Ruiz 2006-14
5License : GPL
6
7Maintainer : Alberto Ruiz
8Stability : provisional
9
10This module reexports all available GSL functions.
11
12The GSL special functions are in the separate package hmatrix-special.
13
14-}
15
16module Numeric.GSL (
17 module Numeric.GSL.Integration
18, module Numeric.GSL.Differentiation
19, module Numeric.GSL.Fourier
20, module Numeric.GSL.Polynomials
21, module Numeric.GSL.Minimization
22, module Numeric.GSL.Root
23, module Numeric.GSL.ODE
24, module Numeric.GSL.Fitting
25, module Data.Complex
26, setErrorHandlerOff
27) where
28
29import Numeric.GSL.Integration
30import Numeric.GSL.Differentiation
31import Numeric.GSL.Fourier
32import Numeric.GSL.Polynomials
33import Numeric.GSL.Minimization
34import Numeric.GSL.Root
35import Numeric.GSL.ODE
36import Numeric.GSL.Fitting
37import Data.Complex
38
39
40-- | This action removes the GSL default error handler (which aborts the program), so that
41-- GSL errors can be handled by Haskell (using Control.Exception) and ghci doesn't abort.
42foreign import ccall unsafe "GSL/gsl-aux.h no_abort_on_error" setErrorHandlerOff :: IO ()
43
diff --git a/packages/hmatrix/src/Numeric/GSL/Differentiation.hs b/packages/hmatrix/src/Numeric/GSL/Differentiation.hs
deleted file mode 100644
index 0fb58ef..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Differentiation.hs
+++ /dev/null
@@ -1,85 +0,0 @@
1{- |
2Module : Numeric.GSL.Differentiation
3Copyright : (c) Alberto Ruiz 2006
4License : GPL
5
6Maintainer : Alberto Ruiz
7Stability : provisional
8
9Numerical differentiation.
10
11<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Differentiation.html#Numerical-Differentiation>
12
13From the GSL manual: \"The functions described in this chapter compute numerical derivatives by finite differencing. An adaptive algorithm is used to find the best choice of finite difference and to estimate the error in the derivative.\"
14-}
15
16
17module Numeric.GSL.Differentiation (
18 derivCentral,
19 derivForward,
20 derivBackward
21) where
22
23import Foreign.C.Types
24import Foreign.Marshal.Alloc(malloc, free)
25import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
26import Foreign.Storable(peek)
27import Numeric.GSL.Internal
28import System.IO.Unsafe(unsafePerformIO)
29
30derivGen ::
31 CInt -- ^ type: 0 central, 1 forward, 2 backward
32 -> Double -- ^ initial step size
33 -> (Double -> Double) -- ^ function
34 -> Double -- ^ point where the derivative is taken
35 -> (Double, Double) -- ^ result and error
36derivGen c h f x = unsafePerformIO $ do
37 r <- malloc
38 e <- malloc
39 fp <- mkfun (\y _ -> f y)
40 c_deriv c fp x h r e // check "deriv"
41 vr <- peek r
42 ve <- peek e
43 let result = (vr,ve)
44 free r
45 free e
46 freeHaskellFunPtr fp
47 return result
48
49foreign import ccall safe "gsl-aux.h deriv"
50 c_deriv :: CInt -> FunPtr (Double -> Ptr () -> Double) -> Double -> Double
51 -> Ptr Double -> Ptr Double -> IO CInt
52
53
54{- | Adaptive central difference algorithm, /gsl_deriv_central/. For example:
55
56>>> let deriv = derivCentral 0.01
57>>> deriv sin (pi/4)
58(0.7071067812000676,1.0600063101654055e-10)
59>>> cos (pi/4)
600.7071067811865476
61
62-}
63derivCentral :: Double -- ^ initial step size
64 -> (Double -> Double) -- ^ function
65 -> Double -- ^ point where the derivative is taken
66 -> (Double, Double) -- ^ result and absolute error
67derivCentral = derivGen 0
68
69-- | Adaptive forward difference algorithm, /gsl_deriv_forward/. The function is evaluated only at points greater than x, and never at x itself. The derivative is returned in result and an estimate of its absolute error is returned in abserr. This function should be used if f(x) has a discontinuity at x, or is undefined for values less than x. A backward derivative can be obtained using a negative step.
70derivForward :: Double -- ^ initial step size
71 -> (Double -> Double) -- ^ function
72 -> Double -- ^ point where the derivative is taken
73 -> (Double, Double) -- ^ result and absolute error
74derivForward = derivGen 1
75
76-- | Adaptive backward difference algorithm, /gsl_deriv_backward/.
77derivBackward ::Double -- ^ initial step size
78 -> (Double -> Double) -- ^ function
79 -> Double -- ^ point where the derivative is taken
80 -> (Double, Double) -- ^ result and absolute error
81derivBackward = derivGen 2
82
83{- | conversion of Haskell functions into function pointers that can be used in the C side
84-}
85foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double))
diff --git a/packages/hmatrix/src/Numeric/GSL/Fitting.hs b/packages/hmatrix/src/Numeric/GSL/Fitting.hs
deleted file mode 100644
index 93fb281..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Fitting.hs
+++ /dev/null
@@ -1,180 +0,0 @@
1{- |
2Module : Numeric.GSL.Fitting
3Copyright : (c) Alberto Ruiz 2010
4License : GPL
5
6Maintainer : Alberto Ruiz
7Stability : provisional
8
9Nonlinear Least-Squares Fitting
10
11<http://www.gnu.org/software/gsl/manual/html_node/Nonlinear-Least_002dSquares-Fitting.html>
12
13The example program in the GSL manual (see examples/fitting.hs):
14
15@
16dat = [
17 ([0.0],([6.0133918608118675],0.1)),
18 ([1.0],([5.5153769909966535],0.1)),
19 ([2.0],([5.261094606015287],0.1)),
20 ...
21 ([39.0],([1.0619821710802808],0.1))]
22
23expModel [a,lambda,b] [t] = [a * exp (-lambda * t) + b]
24
25expModelDer [a,lambda,b] [t] = [[exp (-lambda * t), -t * a * exp(-lambda*t) , 1]]
26
27(sol,path) = fitModelScaled 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
28@
29
30>>> path
31(6><5)
32 [ 1.0, 76.45780563978782, 1.6465931240727802, 1.8147715267618197e-2, 0.6465931240727797
33 , 2.0, 37.683816318260355, 2.858760367632973, 8.092094813253975e-2, 1.4479636296208662
34 , 3.0, 9.5807893736187, 4.948995119561291, 0.11942927999921617, 1.0945766509238248
35 , 4.0, 5.630494933603935, 5.021755718065913, 0.10287787128056883, 1.0338835440862608
36 , 5.0, 5.443976278682909, 5.045204331329302, 0.10405523433131504, 1.019416067207375
37 , 6.0, 5.4439736648994685, 5.045357818922331, 0.10404905846029407, 1.0192487112786812 ]
38>>> sol
39[(5.045357818922331,6.027976702418132e-2),
40(0.10404905846029407,3.157045047172834e-3),
41(1.0192487112786812,3.782067731353722e-2)]
42
43-}
44
45
46module Numeric.GSL.Fitting (
47 -- * Levenberg-Marquardt
48 nlFitting, FittingMethod(..),
49 -- * Utilities
50 fitModelScaled, fitModel
51) where
52
53import Numeric.LinearAlgebra
54import Numeric.GSL.Internal
55
56import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
57import Foreign.C.Types
58import System.IO.Unsafe(unsafePerformIO)
59
60-------------------------------------------------------------------------
61
62type TVV = TV (TV Res)
63type TVM = TV (TM Res)
64
65data FittingMethod = LevenbergMarquardtScaled -- ^ Interface to gsl_multifit_fdfsolver_lmsder. This is a robust and efficient version of the Levenberg-Marquardt algorithm as implemented in the scaled lmder routine in minpack. Minpack was written by Jorge J. More, Burton S. Garbow and Kenneth E. Hillstrom.
66 | LevenbergMarquardt -- ^ This is an unscaled version of the lmder algorithm. The elements of the diagonal scaling matrix D are set to 1. This algorithm may be useful in circumstances where the scaled version of lmder converges too slowly, or the function is already scaled appropriately.
67 deriving (Enum,Eq,Show,Bounded)
68
69
70-- | Nonlinear multidimensional least-squares fitting.
71nlFitting :: FittingMethod
72 -> Double -- ^ absolute tolerance
73 -> Double -- ^ relative tolerance
74 -> Int -- ^ maximum number of iterations allowed
75 -> (Vector Double -> Vector Double) -- ^ function to be minimized
76 -> (Vector Double -> Matrix Double) -- ^ Jacobian
77 -> Vector Double -- ^ starting point
78 -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
79
80nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit
81
82nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
83 let p = dim xiv
84 n = dim (f xiv)
85 fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
86 jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
87 rawpath <- createMatrix RowMajor maxit (2+p)
88 app2 (c_nlfit m fp jp epsabs epsrel (fi maxit) (fi n)) vec xiv mat rawpath "c_nlfit"
89 let it = round (rawpath @@> (maxit-1,0))
90 path = takeRows it rawpath
91 [sol] = toRows $ dropRows (it-1) path
92 freeHaskellFunPtr fp
93 freeHaskellFunPtr jp
94 return (subVector 2 p sol, path)
95
96foreign import ccall safe "nlfit"
97 c_nlfit:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> Double -> CInt -> CInt -> TVM
98
99-------------------------------------------------------
100
101checkdim1 n _p v
102 | dim v == n = v
103 | otherwise = error $ "Error: "++ show n
104 ++ " components expected in the result of the function supplied to nlFitting"
105
106checkdim2 n p m
107 | rows m == n && cols m == p = m
108 | otherwise = error $ "Error: "++ show n ++ "x" ++ show p
109 ++ " Jacobian expected in nlFitting"
110
111------------------------------------------------------------
112
113err (model,deriv) dat vsol = zip sol errs where
114 sol = toList vsol
115 c = max 1 (chi/sqrt (fromIntegral dof))
116 dof = length dat - (rows cov)
117 chi = norm2 (fromList $ cost (resMs model) dat sol)
118 js = fromLists $ jacobian (resDs deriv) dat sol
119 cov = inv $ trans js <.> js
120 errs = toList $ scalar c * sqrt (takeDiag cov)
121
122
123
124-- | Higher level interface to 'nlFitting' 'LevenbergMarquardtScaled'. The optimization function and
125-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
126-- instances (x, (y,sigma)) to be fitted.
127
128fitModelScaled
129 :: Double -- ^ absolute tolerance
130 -> Double -- ^ relative tolerance
131 -> Int -- ^ maximum number of iterations allowed
132 -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
133 -> [(x, ([Double], Double))] -- ^ instances
134 -> [Double] -- ^ starting point
135 -> ([(Double, Double)], Matrix Double) -- ^ (solution, error) and optimization path
136fitModelScaled epsabs epsrel maxit (model,deriv) dt xin = (err (model,deriv) dt sol, path) where
137 (sol,path) = nlFitting LevenbergMarquardtScaled epsabs epsrel maxit
138 (fromList . cost (resMs model) dt . toList)
139 (fromLists . jacobian (resDs deriv) dt . toList)
140 (fromList xin)
141
142
143
144-- | Higher level interface to 'nlFitting' 'LevenbergMarquardt'. The optimization function and
145-- Jacobian are automatically built from a model f vs x = y and its derivatives, and a list of
146-- instances (x,y) to be fitted.
147
148fitModel :: Double -- ^ absolute tolerance
149 -> Double -- ^ relative tolerance
150 -> Int -- ^ maximum number of iterations allowed
151 -> ([Double] -> x -> [Double], [Double] -> x -> [[Double]]) -- ^ (model, derivatives)
152 -> [(x, [Double])] -- ^ instances
153 -> [Double] -- ^ starting point
154 -> ([Double], Matrix Double) -- ^ solution and optimization path
155fitModel epsabs epsrel maxit (model,deriv) dt xin = (toList sol, path) where
156 (sol,path) = nlFitting LevenbergMarquardt epsabs epsrel maxit
157 (fromList . cost (resM model) dt . toList)
158 (fromLists . jacobian (resD deriv) dt . toList)
159 (fromList xin)
160
161cost model ds vs = concatMap (model vs) ds
162
163jacobian modelDer ds vs = concatMap (modelDer vs) ds
164
165-- | Model-to-residual for association pairs with sigma, to be used with 'fitModel'.
166resMs :: ([Double] -> x -> [Double]) -> [Double] -> (x, ([Double], Double)) -> [Double]
167resMs m v = \(x,(ys,s)) -> zipWith (g s) (m v x) ys where g s a b = (a-b)/s
168
169-- | Associated derivative for 'resMs'.
170resDs :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, ([Double], Double)) -> [[Double]]
171resDs m v = \(x,(_,s)) -> map (map (/s)) (m v x)
172
173-- | Model-to-residual for association pairs, to be used with 'fitModel'. It is equivalent
174-- to 'resMs' with all sigmas = 1.
175resM :: ([Double] -> x -> [Double]) -> [Double] -> (x, [Double]) -> [Double]
176resM m v = \(x,ys) -> zipWith g (m v x) ys where g a b = a-b
177
178-- | Associated derivative for 'resM'.
179resD :: ([Double] -> x -> [[Double]]) -> [Double] -> (x, [Double]) -> [[Double]]
180resD m v = \(x,_) -> m v x
diff --git a/packages/hmatrix/src/Numeric/GSL/Fourier.hs b/packages/hmatrix/src/Numeric/GSL/Fourier.hs
deleted file mode 100644
index 734325b..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Fourier.hs
+++ /dev/null
@@ -1,44 +0,0 @@
1{- |
2Module : Numeric.GSL.Fourier
3Copyright : (c) Alberto Ruiz 2006
4License : GPL
5Maintainer : Alberto Ruiz
6Stability : provisional
7
8Fourier Transform.
9
10<http://www.gnu.org/software/gsl/manual/html_node/Fast-Fourier-Transforms.html#Fast-Fourier-Transforms>
11
12-}
13
14module Numeric.GSL.Fourier (
15 fft,
16 ifft
17) where
18
19import Data.Packed
20import Numeric.GSL.Internal
21import Data.Complex
22import Foreign.C.Types
23import System.IO.Unsafe (unsafePerformIO)
24
25genfft code v = unsafePerformIO $ do
26 r <- createVector (dim v)
27 app2 (c_fft code) vec v vec r "fft"
28 return r
29
30foreign import ccall unsafe "gsl-aux.h fft" c_fft :: CInt -> TCV (TCV Res)
31
32
33{- | Fast 1D Fourier transform of a 'Vector' @(@'Complex' 'Double'@)@ using /gsl_fft_complex_forward/. It uses the same scaling conventions as GNU Octave.
34
35>>> fft (fromList [1,2,3,4])
36fromList [10.0 :+ 0.0,(-2.0) :+ 2.0,(-2.0) :+ 0.0,(-2.0) :+ (-2.0)]
37
38-}
39fft :: Vector (Complex Double) -> Vector (Complex Double)
40fft = genfft 0
41
42-- | The inverse of 'fft', using /gsl_fft_complex_inverse/.
43ifft :: Vector (Complex Double) -> Vector (Complex Double)
44ifft = genfft 1
diff --git a/packages/hmatrix/src/Numeric/GSL/IO.hs b/packages/hmatrix/src/Numeric/GSL/IO.hs
deleted file mode 100644
index 0d6031a..0000000
--- a/packages/hmatrix/src/Numeric/GSL/IO.hs
+++ /dev/null
@@ -1,42 +0,0 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Numeric.GSL.IO
4-- Copyright : (c) Alberto Ruiz 2007-14
5-- License : GPL
6-- Maintainer : Alberto Ruiz
7-- Stability : provisional
8--
9-----------------------------------------------------------------------------
10
11module Numeric.GSL.IO (
12 saveMatrix,
13 fwriteVector, freadVector, fprintfVector, fscanfVector,
14 fileDimensions, loadMatrix, fromFile
15) where
16
17import Data.Packed
18import Numeric.GSL.Vector
19import System.Process(readProcess)
20
21
22{- | obtains the number of rows and columns in an ASCII data file
23 (provisionally using unix's wc).
24-}
25fileDimensions :: FilePath -> IO (Int,Int)
26fileDimensions fname = do
27 wcres <- readProcess "wc" ["-w",fname] ""
28 contents <- readFile fname
29 let tot = read . head . words $ wcres
30 c = length . head . dropWhile null . map words . lines $ contents
31 if tot > 0
32 then return (tot `div` c, c)
33 else return (0,0)
34
35-- | Loads a matrix from an ASCII file formatted as a 2D table.
36loadMatrix :: FilePath -> IO (Matrix Double)
37loadMatrix file = fromFile file =<< fileDimensions file
38
39-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
40fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
41fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
42
diff --git a/packages/hmatrix/src/Numeric/GSL/Integration.hs b/packages/hmatrix/src/Numeric/GSL/Integration.hs
deleted file mode 100644
index 9c1d43a..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Integration.hs
+++ /dev/null
@@ -1,246 +0,0 @@
1{- |
2Module : Numeric.GSL.Integration
3Copyright : (c) Alberto Ruiz 2006
4License : GPL
5Maintainer : Alberto Ruiz
6Stability : provisional
7
8Numerical integration routines.
9
10<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html#Numerical-Integration>
11-}
12
13
14module Numeric.GSL.Integration (
15 integrateQNG,
16 integrateQAGS,
17 integrateQAGI,
18 integrateQAGIU,
19 integrateQAGIL,
20 integrateCQUAD
21) where
22
23import Foreign.C.Types
24import Foreign.Marshal.Alloc(malloc, free)
25import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
26import Foreign.Storable(peek)
27import Numeric.GSL.Internal
28import System.IO.Unsafe(unsafePerformIO)
29
30eps = 1e-12
31
32{- | conversion of Haskell functions into function pointers that can be used in the C side
33-}
34foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double))
35
36--------------------------------------------------------------------
37{- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example:
38
39>>> let quad = integrateQAGS 1E-9 1000
40>>> let f a x = x**(-0.5) * log (a*x)
41>>> quad (f 1) 0 1
42(-3.999999999999974,4.871658632055187e-13)
43
44-}
45
46integrateQAGS :: Double -- ^ precision (e.g. 1E-9)
47 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
48 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
49 -> Double -- ^ a
50 -> Double -- ^ b
51 -> (Double, Double) -- ^ result of the integration and error
52integrateQAGS prec n f a b = unsafePerformIO $ do
53 r <- malloc
54 e <- malloc
55 fp <- mkfun (\x _ -> f x)
56 c_integrate_qags fp a b eps prec (fromIntegral n) r e // check "integrate_qags"
57 vr <- peek r
58 ve <- peek e
59 let result = (vr,ve)
60 free r
61 free e
62 freeHaskellFunPtr fp
63 return result
64
65foreign import ccall safe "integrate_qags" c_integrate_qags
66 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
67 -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
68
69-----------------------------------------------------------------
70{- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example:
71
72>>> let quad = integrateQNG 1E-6
73>>> quad (\x -> 4/(1+x*x)) 0 1
74(3.141592653589793,3.487868498008632e-14)
75
76-}
77
78integrateQNG :: Double -- ^ precision (e.g. 1E-9)
79 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
80 -> Double -- ^ a
81 -> Double -- ^ b
82 -> (Double, Double) -- ^ result of the integration and error
83integrateQNG prec f a b = unsafePerformIO $ do
84 r <- malloc
85 e <- malloc
86 fp <- mkfun (\x _ -> f x)
87 c_integrate_qng fp a b eps prec r e // check "integrate_qng"
88 vr <- peek r
89 ve <- peek e
90 let result = (vr,ve)
91 free r
92 free e
93 freeHaskellFunPtr fp
94 return result
95
96foreign import ccall safe "integrate_qng" c_integrate_qng
97 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
98 -> Double -> Double -> Ptr Double -> Ptr Double -> IO CInt
99
100--------------------------------------------------------------------
101{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS).
102For example:
103
104>>> let quad = integrateQAGI 1E-9 1000
105>>> let f a x = exp(-a * x^2)
106>>> quad (f 0.5)
107(2.5066282746310002,6.229215880648858e-11)
108
109-}
110
111integrateQAGI :: Double -- ^ precision (e.g. 1E-9)
112 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
113 -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf)
114 -> (Double, Double) -- ^ result of the integration and error
115integrateQAGI prec n f = unsafePerformIO $ do
116 r <- malloc
117 e <- malloc
118 fp <- mkfun (\x _ -> f x)
119 c_integrate_qagi fp eps prec (fromIntegral n) r e // check "integrate_qagi"
120 vr <- peek r
121 ve <- peek e
122 let result = (vr,ve)
123 free r
124 free e
125 freeHaskellFunPtr fp
126 return result
127
128foreign import ccall safe "integrate_qagi" c_integrate_qagi
129 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
130 -> CInt -> Ptr Double -> Ptr Double -> IO CInt
131
132--------------------------------------------------------------------
133{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf).
134For example:
135
136>>> let quad = integrateQAGIU 1E-9 1000
137>>> let f a x = exp(-a * x^2)
138>>> quad (f 0.5) 0
139(1.2533141373155001,3.114607940324429e-11)
140
141-}
142
143integrateQAGIU :: Double -- ^ precision (e.g. 1E-9)
144 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
145 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
146 -> Double -- ^ a
147 -> (Double, Double) -- ^ result of the integration and error
148integrateQAGIU prec n f a = unsafePerformIO $ do
149 r <- malloc
150 e <- malloc
151 fp <- mkfun (\x _ -> f x)
152 c_integrate_qagiu fp a eps prec (fromIntegral n) r e // check "integrate_qagiu"
153 vr <- peek r
154 ve <- peek e
155 let result = (vr,ve)
156 free r
157 free e
158 freeHaskellFunPtr fp
159 return result
160
161foreign import ccall safe "integrate_qagiu" c_integrate_qagiu
162 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
163 -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
164
165--------------------------------------------------------------------
166{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b).
167For example:
168
169>>> let quad = integrateQAGIL 1E-9 1000
170>>> let f a x = exp(-a * x^2)
171>>> quad (f 0.5) 0
172(1.2533141373155001,3.114607940324429e-11)
173
174-}
175
176integrateQAGIL :: Double -- ^ precision (e.g. 1E-9)
177 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
178 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
179 -> Double -- ^ b
180 -> (Double, Double) -- ^ result of the integration and error
181integrateQAGIL prec n f b = unsafePerformIO $ do
182 r <- malloc
183 e <- malloc
184 fp <- mkfun (\x _ -> f x)
185 c_integrate_qagil fp b eps prec (fromIntegral n) r e // check "integrate_qagil"
186 vr <- peek r
187 ve <- peek e
188 let result = (vr,ve)
189 free r
190 free e
191 freeHaskellFunPtr fp
192 return result
193
194foreign import ccall safe "gsl-aux.h integrate_qagil" c_integrate_qagil
195 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
196 -> Double -> CInt -> Ptr Double -> Ptr Double -> IO CInt
197
198
199--------------------------------------------------------------------
200{- | Numerical integration using /gsl_integration_cquad/ (quadrature
201for general integrands). From the GSL manual:
202
203@CQUAD is a new doubly-adaptive general-purpose quadrature routine
204which can handle most types of singularities, non-numerical function
205values such as Inf or NaN, as well as some divergent integrals. It
206generally requires more function evaluations than the integration
207routines in QUADPACK, yet fails less often for difficult integrands.@
208
209For example:
210
211>>> let quad = integrateCQUAD 1E-12 1000
212>>> let f a x = exp(-a * x^2)
213>>> quad (f 0.5) 2 5
214(5.7025405463957006e-2,9.678874441303705e-16,95)
215
216Unlike other quadrature methods, integrateCQUAD also returns the
217number of function evaluations required.
218
219-}
220
221integrateCQUAD :: Double -- ^ precision (e.g. 1E-9)
222 -> Int -- ^ size of auxiliary workspace (e.g. 1000)
223 -> (Double -> Double) -- ^ function to be integrated on the interval (a, b)
224 -> Double -- ^ a
225 -> Double -- ^ b
226 -> (Double, Double, Int) -- ^ result of the integration, error and number of function evaluations performed
227integrateCQUAD prec n f a b = unsafePerformIO $ do
228 r <- malloc
229 e <- malloc
230 neval <- malloc
231 fp <- mkfun (\x _ -> f x)
232 c_integrate_cquad fp a b eps prec (fromIntegral n) r e neval // check "integrate_cquad"
233 vr <- peek r
234 ve <- peek e
235 vneval <- peek neval
236 let result = (vr,ve,vneval)
237 free r
238 free e
239 free neval
240 freeHaskellFunPtr fp
241 return result
242
243foreign import ccall safe "integrate_cquad" c_integrate_cquad
244 :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double
245 -> Double -> Double -> CInt -> Ptr Double -> Ptr Double -> Ptr Int -> IO CInt
246
diff --git a/packages/hmatrix/src/Numeric/GSL/Internal.hs b/packages/hmatrix/src/Numeric/GSL/Internal.hs
deleted file mode 100644
index a1c4e0c..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Internal.hs
+++ /dev/null
@@ -1,126 +0,0 @@
1-- |
2-- Module : Numeric.GSL.Internal
3-- Copyright : (c) Alberto Ruiz 2009
4-- License : GPL
5-- Maintainer : Alberto Ruiz
6-- Stability : provisional
7--
8--
9-- Auxiliary functions.
10--
11
12
13module Numeric.GSL.Internal(
14 iv,
15 mkVecfun,
16 mkVecVecfun,
17 mkDoubleVecVecfun,
18 mkDoublefun,
19 aux_vTov,
20 mkVecMatfun,
21 mkDoubleVecMatfun,
22 aux_vTom,
23 createV,
24 createMIO,
25 module Data.Packed.Development,
26 check,
27 Res,TV,TM,TCV,TCM
28) where
29
30import Data.Packed
31import Data.Packed.Development hiding (check)
32import Data.Complex
33
34import Foreign.Marshal.Array(copyArray)
35import Foreign.Ptr(Ptr, FunPtr)
36import Foreign.C.Types
37import Foreign.C.String(peekCString)
38import System.IO.Unsafe(unsafePerformIO)
39import Data.Vector.Storable(unsafeWith)
40import Control.Monad(when)
41
42iv :: (Vector Double -> Double) -> (CInt -> Ptr Double -> Double)
43iv f n p = f (createV (fromIntegral n) copy "iv") where
44 copy n' q = do
45 copyArray q p (fromIntegral n')
46 return 0
47
48-- | conversion of Haskell functions into function pointers that can be used in the C side
49foreign import ccall safe "wrapper"
50 mkVecfun :: (CInt -> Ptr Double -> Double)
51 -> IO( FunPtr (CInt -> Ptr Double -> Double))
52
53foreign import ccall safe "wrapper"
54 mkVecVecfun :: TVV -> IO (FunPtr TVV)
55
56foreign import ccall safe "wrapper"
57 mkDoubleVecVecfun :: (Double -> TVV) -> IO (FunPtr (Double -> TVV))
58
59foreign import ccall safe "wrapper"
60 mkDoublefun :: (Double -> Double) -> IO (FunPtr (Double -> Double))
61
62aux_vTov :: (Vector Double -> Vector Double) -> TVV
63aux_vTov f n p nr r = g where
64 v = f x
65 x = createV (fromIntegral n) copy "aux_vTov"
66 copy n' q = do
67 copyArray q p (fromIntegral n')
68 return 0
69 g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral nr)
70 return 0
71
72foreign import ccall safe "wrapper"
73 mkVecMatfun :: TVM -> IO (FunPtr TVM)
74
75foreign import ccall safe "wrapper"
76 mkDoubleVecMatfun :: (Double -> TVM) -> IO (FunPtr (Double -> TVM))
77
78aux_vTom :: (Vector Double -> Matrix Double) -> TVM
79aux_vTom f n p rr cr r = g where
80 v = flatten $ f x
81 x = createV (fromIntegral n) copy "aux_vTov"
82 copy n' q = do
83 copyArray q p (fromIntegral n')
84 return 0
85 g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
86 return 0
87
88createV n fun msg = unsafePerformIO $ do
89 r <- createVector n
90 app1 fun vec r msg
91 return r
92
93createMIO r c fun msg = do
94 res <- createMatrix RowMajor r c
95 app1 fun mat res msg
96 return res
97
98--------------------------------------------------------------------------------
99
100-- | check the error code
101check :: String -> IO CInt -> IO ()
102check msg f = do
103 err <- f
104 when (err/=0) $ do
105 ps <- gsl_strerror err
106 s <- peekCString ps
107 error (msg++": "++s)
108 return ()
109
110-- | description of GSL error codes
111foreign import ccall unsafe "gsl_strerror" gsl_strerror :: CInt -> IO (Ptr CChar)
112
113type PF = Ptr Float
114type PD = Ptr Double
115type PQ = Ptr (Complex Float)
116type PC = Ptr (Complex Double)
117
118type Res = IO CInt
119type TV x = CInt -> PD -> x
120type TM x = CInt -> CInt -> PD -> x
121type TCV x = CInt -> PC -> x
122type TCM x = CInt -> CInt -> PC -> x
123
124type TVV = TV (TV Res)
125type TVM = TV (TM Res)
126
diff --git a/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs b/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
deleted file mode 100644
index 17e2258..0000000
--- a/packages/hmatrix/src/Numeric/GSL/LinearAlgebra.hs
+++ /dev/null
@@ -1,135 +0,0 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Numeric.GSL.LinearAlgebra
4-- Copyright : (c) Alberto Ruiz 2007-14
5-- License : GPL
6-- Maintainer : Alberto Ruiz
7-- Stability : provisional
8--
9-----------------------------------------------------------------------------
10
11module Numeric.GSL.LinearAlgebra (
12 RandDist(..), randomVector,
13 saveMatrix,
14 fwriteVector, freadVector, fprintfVector, fscanfVector,
15 fileDimensions, loadMatrix, fromFile
16) where
17
18import Data.Packed
19import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
20
21import Foreign.Marshal.Alloc(free)
22import Foreign.Ptr(Ptr)
23import Foreign.C.Types
24import Foreign.C.String(newCString)
25import System.IO.Unsafe(unsafePerformIO)
26import System.Process(readProcess)
27
28fromei x = fromIntegral (fromEnum x) :: CInt
29
30-----------------------------------------------------------------------
31
32data RandDist = Uniform -- ^ uniform distribution in [0,1)
33 | Gaussian -- ^ normal distribution with mean zero and standard deviation one
34 deriving Enum
35
36-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
37randomVector :: Int -- ^ seed
38 -> RandDist -- ^ distribution
39 -> Int -- ^ vector size
40 -> Vector Double
41randomVector seed dist n = unsafePerformIO $ do
42 r <- createVector n
43 app1 (c_random_vector (fi seed) ((fi.fromEnum) dist)) vec r "randomVector"
44 return r
45
46foreign import ccall unsafe "random_vector" c_random_vector :: CInt -> CInt -> TV
47
48--------------------------------------------------------------------------------
49
50-- | Saves a matrix as 2D ASCII table.
51saveMatrix :: FilePath
52 -> String -- ^ format (%f, %g, %e)
53 -> Matrix Double
54 -> IO ()
55saveMatrix filename fmt m = do
56 charname <- newCString filename
57 charfmt <- newCString fmt
58 let o = if orderOf m == RowMajor then 1 else 0
59 app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
60 free charname
61 free charfmt
62
63foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
64
65--------------------------------------------------------------------------------
66
67-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
68fscanfVector :: FilePath -> Int -> IO (Vector Double)
69fscanfVector filename n = do
70 charname <- newCString filename
71 res <- createVector n
72 app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
73 free charname
74 return res
75
76foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
77
78-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
79fprintfVector :: FilePath -> String -> Vector Double -> IO ()
80fprintfVector filename fmt v = do
81 charname <- newCString filename
82 charfmt <- newCString fmt
83 app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
84 free charname
85 free charfmt
86
87foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
88
89-- | Loads a vector from a binary file (the number of elements must be known in advance).
90freadVector :: FilePath -> Int -> IO (Vector Double)
91freadVector filename n = do
92 charname <- newCString filename
93 res <- createVector n
94 app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
95 free charname
96 return res
97
98foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
99
100-- | Saves the elements of a vector to a binary file.
101fwriteVector :: FilePath -> Vector Double -> IO ()
102fwriteVector filename v = do
103 charname <- newCString filename
104 app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
105 free charname
106
107foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
108
109type PD = Ptr Double --
110type TV = CInt -> PD -> IO CInt --
111type TM = CInt -> CInt -> PD -> IO CInt --
112
113--------------------------------------------------------------------------------
114
115{- | obtains the number of rows and columns in an ASCII data file
116 (provisionally using unix's wc).
117-}
118fileDimensions :: FilePath -> IO (Int,Int)
119fileDimensions fname = do
120 wcres <- readProcess "wc" ["-w",fname] ""
121 contents <- readFile fname
122 let tot = read . head . words $ wcres
123 c = length . head . dropWhile null . map words . lines $ contents
124 if tot > 0
125 then return (tot `div` c, c)
126 else return (0,0)
127
128-- | Loads a matrix from an ASCII file formatted as a 2D table.
129loadMatrix :: FilePath -> IO (Matrix Double)
130loadMatrix file = fromFile file =<< fileDimensions file
131
132-- | Loads a matrix from an ASCII file (the number of rows and columns must be known in advance).
133fromFile :: FilePath -> (Int,Int) -> IO (Matrix Double)
134fromFile filename (r,c) = reshape c `fmap` fscanfVector filename (r*c)
135
diff --git a/packages/hmatrix/src/Numeric/GSL/Minimization.hs b/packages/hmatrix/src/Numeric/GSL/Minimization.hs
deleted file mode 100644
index 056d463..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Minimization.hs
+++ /dev/null
@@ -1,222 +0,0 @@
1{- |
2Module : Numeric.GSL.Minimization
3Copyright : (c) Alberto Ruiz 2006-9
4License : GPL
5Maintainer : Alberto Ruiz
6Stability : provisional
7
8Minimization of a multidimensional function using some of the algorithms described in:
9
10<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Minimization.html>
11
12The example in the GSL manual:
13
14@
15f [x,y] = 10*(x-1)^2 + 20*(y-2)^2 + 30
16
17main = do
18 let (s,p) = minimize NMSimplex2 1E-2 30 [1,1] f [5,7]
19 print s
20 print p
21@
22
23>>> main
24[0.9920430849306288,1.9969168063253182]
25 0.000 512.500 1.130 6.500 5.000
26 1.000 290.625 1.409 5.250 4.000
27 2.000 290.625 1.409 5.250 4.000
28 3.000 252.500 1.409 5.500 1.000
29 ...
3022.000 30.001 0.013 0.992 1.997
3123.000 30.001 0.008 0.992 1.997
32
33The path to the solution can be graphically shown by means of:
34
35@'Graphics.Plot.mplot' $ drop 3 ('toColumns' p)@
36
37Taken from the GSL manual:
38
39The vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm is a quasi-Newton method which builds up an approximation to the second derivatives of the function f using the difference between successive gradient vectors. By combining the first and second derivatives the algorithm is able to take Newton-type steps towards the function minimum, assuming quadratic behavior in that region.
40
41The bfgs2 version of this minimizer is the most efficient version available, and is a faithful implementation of the line minimization scheme described in Fletcher's Practical Methods of Optimization, Algorithms 2.6.2 and 2.6.4. It supercedes the original bfgs routine and requires substantially fewer function and gradient evaluations. The user-supplied tolerance tol corresponds to the parameter \sigma used by Fletcher. A value of 0.1 is recommended for typical use (larger values correspond to less accurate line searches).
42
43The nmsimplex2 version is a new O(N) implementation of the earlier O(N^2) nmsimplex minimiser. It calculates the size of simplex as the rms distance of each vertex from the center rather than the mean distance, which has the advantage of allowing a linear update.
44
45-}
46
47
48module Numeric.GSL.Minimization (
49 minimize, minimizeV, MinimizeMethod(..),
50 minimizeD, minimizeVD, MinimizeMethodD(..),
51 uniMinimize, UniMinimizeMethod(..),
52
53 minimizeNMSimplex,
54 minimizeConjugateGradient,
55 minimizeVectorBFGS2
56) where
57
58
59import Data.Packed
60import Numeric.GSL.Internal
61
62import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
63import Foreign.C.Types
64import System.IO.Unsafe(unsafePerformIO)
65
66------------------------------------------------------------------------
67
68{-# DEPRECATED minimizeNMSimplex "use minimize NMSimplex2 eps maxit sizes f xi" #-}
69minimizeNMSimplex f xi szs eps maxit = minimize NMSimplex eps maxit szs f xi
70
71{-# DEPRECATED minimizeConjugateGradient "use minimizeD ConjugateFR eps maxit step tol f g xi" #-}
72minimizeConjugateGradient step tol eps maxit f g xi = minimizeD ConjugateFR eps maxit step tol f g xi
73
74{-# DEPRECATED minimizeVectorBFGS2 "use minimizeD VectorBFGS2 eps maxit step tol f g xi" #-}
75minimizeVectorBFGS2 step tol eps maxit f g xi = minimizeD VectorBFGS2 eps maxit step tol f g xi
76
77-------------------------------------------------------------------------
78
79data UniMinimizeMethod = GoldenSection
80 | BrentMini
81 | QuadGolden
82 deriving (Enum, Eq, Show, Bounded)
83
84-- | Onedimensional minimization.
85
86uniMinimize :: UniMinimizeMethod -- ^ The method used.
87 -> Double -- ^ desired precision of the solution
88 -> Int -- ^ maximum number of iterations allowed
89 -> (Double -> Double) -- ^ function to minimize
90 -> Double -- ^ guess for the location of the minimum
91 -> Double -- ^ lower bound of search interval
92 -> Double -- ^ upper bound of search interval
93 -> (Double, Matrix Double) -- ^ solution and optimization path
94
95uniMinimize method epsrel maxit fun xmin xl xu = uniMinimizeGen (fi (fromEnum method)) fun xmin xl xu epsrel maxit
96
97uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do
98 fp <- mkDoublefun f
99 rawpath <- createMIO maxit 4
100 (c_uniMinize m fp epsrel (fi maxit) xmin xl xu)
101 "uniMinimize"
102 let it = round (rawpath @@> (maxit-1,0))
103 path = takeRows it rawpath
104 [sol] = toLists $ dropRows (it-1) path
105 freeHaskellFunPtr fp
106 return (sol !! 1, path)
107
108
109foreign import ccall safe "uniMinimize"
110 c_uniMinize:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> Double -> TM Res
111
112data MinimizeMethod = NMSimplex
113 | NMSimplex2
114 deriving (Enum,Eq,Show,Bounded)
115
116-- | Minimization without derivatives
117minimize :: MinimizeMethod
118 -> Double -- ^ desired precision of the solution (size test)
119 -> Int -- ^ maximum number of iterations allowed
120 -> [Double] -- ^ sizes of the initial search box
121 -> ([Double] -> Double) -- ^ function to minimize
122 -> [Double] -- ^ starting point
123 -> ([Double], Matrix Double) -- ^ solution vector and optimization path
124
125-- | Minimization without derivatives (vector version)
126minimizeV :: MinimizeMethod
127 -> Double -- ^ desired precision of the solution (size test)
128 -> Int -- ^ maximum number of iterations allowed
129 -> Vector Double -- ^ sizes of the initial search box
130 -> (Vector Double -> Double) -- ^ function to minimize
131 -> Vector Double -- ^ starting point
132 -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
133
134minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi)
135 where v2l (v,m) = (toList v, m)
136
137ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
138
139minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
140 let n = dim xiv
141 fp <- mkVecfun (iv f)
142 rawpath <- ww2 vec xiv vec szv $ \xiv' szv' ->
143 createMIO maxit (n+3)
144 (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
145 "minimize"
146 let it = round (rawpath @@> (maxit-1,0))
147 path = takeRows it rawpath
148 sol = flatten $ dropColumns 3 $ dropRows (it-1) path
149 freeHaskellFunPtr fp
150 return (sol, path)
151
152
153foreign import ccall safe "gsl-aux.h minimize"
154 c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TV(TV(TM Res))
155
156----------------------------------------------------------------------------------
157
158
159data MinimizeMethodD = ConjugateFR
160 | ConjugatePR
161 | VectorBFGS
162 | VectorBFGS2
163 | SteepestDescent
164 deriving (Enum,Eq,Show,Bounded)
165
166-- | Minimization with derivatives.
167minimizeD :: MinimizeMethodD
168 -> Double -- ^ desired precision of the solution (gradient test)
169 -> Int -- ^ maximum number of iterations allowed
170 -> Double -- ^ size of the first trial step
171 -> Double -- ^ tol (precise meaning depends on method)
172 -> ([Double] -> Double) -- ^ function to minimize
173 -> ([Double] -> [Double]) -- ^ gradient
174 -> [Double] -- ^ starting point
175 -> ([Double], Matrix Double) -- ^ solution vector and optimization path
176
177-- | Minimization with derivatives (vector version)
178minimizeVD :: MinimizeMethodD
179 -> Double -- ^ desired precision of the solution (gradient test)
180 -> Int -- ^ maximum number of iterations allowed
181 -> Double -- ^ size of the first trial step
182 -> Double -- ^ tol (precise meaning depends on method)
183 -> (Vector Double -> Double) -- ^ function to minimize
184 -> (Vector Double -> Vector Double) -- ^ gradient
185 -> Vector Double -- ^ starting point
186 -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
187
188minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD
189 method eps maxit istep tol (f.toList) (fromList.df.toList) (fromList xi)
190 where v2l (v,m) = (toList v, m)
191
192
193minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
194 let n = dim xiv
195 f' = f
196 df' = (checkdim1 n . df)
197 fp <- mkVecfun (iv f')
198 dfp <- mkVecVecfun (aux_vTov df')
199 rawpath <- vec xiv $ \xiv' ->
200 createMIO maxit (n+2)
201 (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
202 "minimizeD"
203 let it = round (rawpath @@> (maxit-1,0))
204 path = takeRows it rawpath
205 sol = flatten $ dropColumns 2 $ dropRows (it-1) path
206 freeHaskellFunPtr fp
207 freeHaskellFunPtr dfp
208 return (sol,path)
209
210foreign import ccall safe "gsl-aux.h minimizeD"
211 c_minimizeD :: CInt
212 -> FunPtr (CInt -> Ptr Double -> Double)
213 -> FunPtr (TV (TV Res))
214 -> Double -> Double -> Double -> CInt
215 -> TV (TM Res)
216
217---------------------------------------------------------------------
218
219checkdim1 n v
220 | dim v == n = v
221 | otherwise = error $ "Error: "++ show n
222 ++ " components expected in the result of the gradient supplied to minimizeD"
diff --git a/packages/hmatrix/src/Numeric/GSL/ODE.hs b/packages/hmatrix/src/Numeric/GSL/ODE.hs
deleted file mode 100644
index 7549a65..0000000
--- a/packages/hmatrix/src/Numeric/GSL/ODE.hs
+++ /dev/null
@@ -1,140 +0,0 @@
1{- |
2Module : Numeric.GSL.ODE
3Copyright : (c) Alberto Ruiz 2010
4License : GPL
5Maintainer : Alberto Ruiz
6Stability : provisional
7
8Solution of ordinary differential equation (ODE) initial value problems.
9
10<http://www.gnu.org/software/gsl/manual/html_node/Ordinary-Differential-Equations.html>
11
12A simple example:
13
14@
15import Numeric.GSL.ODE
16import Numeric.LinearAlgebra
17import Numeric.LinearAlgebra.Util(mplot)
18
19xdot t [x,v] = [v, -0.95*x - 0.1*v]
20
21ts = linspace 100 (0,20 :: Double)
22
23sol = odeSolve xdot [10,0] ts
24
25main = mplot (ts : toColumns sol)
26@
27
28-}
29-----------------------------------------------------------------------------
30
31module Numeric.GSL.ODE (
32 odeSolve, odeSolveV, ODEMethod(..), Jacobian
33) where
34
35import Data.Packed
36import Numeric.GSL.Internal
37
38import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)
39import Foreign.C.Types
40import System.IO.Unsafe(unsafePerformIO)
41
42-------------------------------------------------------------------------
43
44type TVV = TV (TV Res)
45type TVM = TV (TM Res)
46type TVVM = TV (TV (TM Res))
47
48type Jacobian = Double -> Vector Double -> Matrix Double
49
50-- | Stepping functions
51data ODEMethod = RK2 -- ^ Embedded Runge-Kutta (2, 3) method.
52 | RK4 -- ^ 4th order (classical) Runge-Kutta. The error estimate is obtained by halving the step-size. For more efficient estimate of the error, use the embedded methods.
53 | RKf45 -- ^ Embedded Runge-Kutta-Fehlberg (4, 5) method. This method is a good general-purpose integrator.
54 | RKck -- ^ Embedded Runge-Kutta Cash-Karp (4, 5) method.
55 | RK8pd -- ^ Embedded Runge-Kutta Prince-Dormand (8,9) method.
56 | RK2imp Jacobian -- ^ Implicit 2nd order Runge-Kutta at Gaussian points.
57 | RK4imp Jacobian -- ^ Implicit 4th order Runge-Kutta at Gaussian points.
58 | BSimp Jacobian -- ^ Implicit Bulirsch-Stoer method of Bader and Deuflhard. The method is generally suitable for stiff problems.
59 | RK1imp Jacobian -- ^ Implicit Gaussian first order Runge-Kutta. Also known as implicit Euler or backward Euler method. Error estimation is carried out by the step doubling method.
60 | MSAdams -- ^ A variable-coefficient linear multistep Adams method in Nordsieck form. This stepper uses explicit Adams-Bashforth (predictor) and implicit Adams-Moulton (corrector) methods in P(EC)^m functional iteration mode. Method order varies dynamically between 1 and 12.
61 | MSBDF Jacobian -- ^ A variable-coefficient linear multistep backward differentiation formula (BDF) method in Nordsieck form. This stepper uses the explicit BDF formula as predictor and implicit BDF formula as corrector. A modified Newton iteration method is used to solve the system of non-linear equations. Method order varies dynamically between 1 and 5. The method is generally suitable for stiff problems.
62
63
64-- | A version of 'odeSolveV' with reasonable default parameters and system of equations defined using lists.
65odeSolve
66 :: (Double -> [Double] -> [Double]) -- ^ xdot(t,x)
67 -> [Double] -- ^ initial conditions
68 -> Vector Double -- ^ desired solution times
69 -> Matrix Double -- ^ solution
70odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
71 where hi = (ts@>1 - ts@>0)/100
72 epsAbs = 1.49012e-08
73 epsRel = 1.49012e-08
74 l2v f = \t -> fromList . f t . toList
75
76-- | Evolution of the system with adaptive step-size control.
77odeSolveV
78 :: ODEMethod
79 -> Double -- ^ initial step size
80 -> Double -- ^ absolute tolerance for the state vector
81 -> Double -- ^ relative tolerance for the state vector
82 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
83 -> Vector Double -- ^ initial conditions
84 -> Vector Double -- ^ desired solution times
85 -> Matrix Double -- ^ solution
86odeSolveV RK2 = odeSolveV' 0 Nothing
87odeSolveV RK4 = odeSolveV' 1 Nothing
88odeSolveV RKf45 = odeSolveV' 2 Nothing
89odeSolveV RKck = odeSolveV' 3 Nothing
90odeSolveV RK8pd = odeSolveV' 4 Nothing
91odeSolveV (RK2imp jac) = odeSolveV' 5 (Just jac)
92odeSolveV (RK4imp jac) = odeSolveV' 6 (Just jac)
93odeSolveV (BSimp jac) = odeSolveV' 7 (Just jac)
94odeSolveV (RK1imp jac) = odeSolveV' 8 (Just jac)
95odeSolveV MSAdams = odeSolveV' 9 Nothing
96odeSolveV (MSBDF jac) = odeSolveV' 10 (Just jac)
97
98
99odeSolveV'
100 :: CInt
101 -> Maybe (Double -> Vector Double -> Matrix Double) -- ^ optional jacobian
102 -> Double -- ^ initial step size
103 -> Double -- ^ absolute tolerance for the state vector
104 -> Double -- ^ relative tolerance for the state vector
105 -> (Double -> Vector Double -> Vector Double) -- ^ xdot(t,x)
106 -> Vector Double -- ^ initial conditions
107 -> Vector Double -- ^ desired solution times
108 -> Matrix Double -- ^ solution
109odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do
110 let n = dim xiv
111 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
112 jp <- case mbjac of
113 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
114 Nothing -> return nullFunPtr
115 sol <- vec xiv $ \xiv' ->
116 vec (checkTimes ts) $ \ts' ->
117 createMIO (dim ts) n
118 (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
119 "ode"
120 freeHaskellFunPtr fp
121 return sol
122
123foreign import ccall safe "ode"
124 ode_c :: CInt -> Double -> Double -> Double -> FunPtr (Double -> TVV) -> FunPtr (Double -> TVM) -> TVVM
125
126-------------------------------------------------------
127
128checkdim1 n v
129 | dim v == n = v
130 | otherwise = error $ "Error: "++ show n
131 ++ " components expected in the result of the function supplied to odeSolve"
132
133checkdim2 n m
134 | rows m == n && cols m == n = m
135 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
136 ++ " Jacobian expected in odeSolve"
137
138checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
139 | otherwise = error "odeSolve requires increasing times"
140 where ts' = toList ts
diff --git a/packages/hmatrix/src/Numeric/GSL/Polynomials.hs b/packages/hmatrix/src/Numeric/GSL/Polynomials.hs
deleted file mode 100644
index b1be85d..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Polynomials.hs
+++ /dev/null
@@ -1,57 +0,0 @@
1{- |
2Module : Numeric.GSL.Polynomials
3Copyright : (c) Alberto Ruiz 2006
4License : GPL
5Maintainer : Alberto Ruiz
6Stability : provisional
7
8Polynomials.
9
10<http://www.gnu.org/software/gsl/manual/html_node/General-Polynomial-Equations.html#General-Polynomial-Equations>
11
12-}
13
14
15module Numeric.GSL.Polynomials (
16 polySolve
17) where
18
19import Data.Packed
20import Numeric.GSL.Internal
21import Data.Complex
22import System.IO.Unsafe (unsafePerformIO)
23
24#if __GLASGOW_HASKELL__ >= 704
25import Foreign.C.Types (CInt(..))
26#endif
27
28{- | Solution of general polynomial equations, using /gsl_poly_complex_solve/.
29
30For example, the three solutions of x^3 + 8 = 0
31
32>>> polySolve [8,0,0,1]
33[(-2.0) :+ 0.0,1.0 :+ 1.7320508075688776,1.0 :+ (-1.7320508075688776)]
34
35
36The example in the GSL manual: To find the roots of x^5 -1 = 0:
37
38>>> polySolve [-1, 0, 0, 0, 0, 1]
39[(-0.8090169943749472) :+ 0.5877852522924731,
40(-0.8090169943749472) :+ (-0.5877852522924731),
410.30901699437494756 :+ 0.9510565162951535,
420.30901699437494756 :+ (-0.9510565162951535),
431.0000000000000002 :+ 0.0]
44
45-}
46polySolve :: [Double] -> [Complex Double]
47polySolve = toList . polySolve' . fromList
48
49polySolve' :: Vector Double -> Vector (Complex Double)
50polySolve' v | dim v > 1 = unsafePerformIO $ do
51 r <- createVector (dim v-1)
52 app2 c_polySolve vec v vec r "polySolve"
53 return r
54 | otherwise = error "polySolve on a polynomial of degree zero"
55
56foreign import ccall unsafe "gsl-aux.h polySolve" c_polySolve:: TV (TCV Res)
57
diff --git a/packages/hmatrix/src/Numeric/GSL/Random.hs b/packages/hmatrix/src/Numeric/GSL/Random.hs
deleted file mode 100644
index 2872b17..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Random.hs
+++ /dev/null
@@ -1,84 +0,0 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Numeric.GSL.Random
4-- Copyright : (c) Alberto Ruiz 2009-14
5-- License : GPL
6--
7-- Maintainer : Alberto Ruiz
8-- Stability : provisional
9--
10-- Random vectors and matrices.
11--
12-----------------------------------------------------------------------------
13
14module Numeric.GSL.Random (
15 Seed,
16 RandDist(..),
17 randomVector,
18 gaussianSample,
19 uniformSample,
20 rand, randn
21) where
22
23import Numeric.GSL.Vector
24import Numeric.Container
25import Numeric.LinearAlgebra(Seed,RandDist(..),cholSH)
26import System.Random(randomIO)
27
28
29
30
31-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
32-- Gaussian distribution.
33gaussianSample :: Seed
34 -> Int -- ^ number of rows
35 -> Vector Double -- ^ mean vector
36 -> Matrix Double -- ^ covariance matrix
37 -> Matrix Double -- ^ result
38gaussianSample seed n med cov = m where
39 c = dim med
40 meds = konst' 1 n `outer` med
41 rs = reshape c $ randomVector seed Gaussian (c * n)
42 m = rs `mXm` cholSH cov `add` meds
43
44-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
45-- uniform distribution.
46uniformSample :: Seed
47 -> Int -- ^ number of rows
48 -> [(Double,Double)] -- ^ ranges for each column
49 -> Matrix Double -- ^ result
50uniformSample seed n rgs = m where
51 (as,bs) = unzip rgs
52 a = fromList as
53 cs = zipWith subtract as bs
54 d = dim a
55 dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
56 am = konst' 1 n `outer` a
57 m = fromColumns (zipWith scale cs dat) `add` am
58
59-- | pseudorandom matrix with uniform elements between 0 and 1
60randm :: RandDist
61 -> Int -- ^ rows
62 -> Int -- ^ columns
63 -> IO (Matrix Double)
64randm d r c = do
65 seed <- randomIO
66 return (reshape c $ randomVector seed d (r*c))
67
68-- | pseudorandom matrix with uniform elements between 0 and 1
69rand :: Int -> Int -> IO (Matrix Double)
70rand = randm Uniform
71
72{- | pseudorandom matrix with normal elements
73
74>>> x <- randn 3 5
75>>> disp 3 x
763x5
770.386 -1.141 0.491 -0.510 1.512
780.069 -0.919 1.022 -0.181 0.745
790.313 -0.670 -0.097 -1.575 -0.583
80
81-}
82randn :: Int -> Int -> IO (Matrix Double)
83randn = randm Gaussian
84
diff --git a/packages/hmatrix/src/Numeric/GSL/Root.hs b/packages/hmatrix/src/Numeric/GSL/Root.hs
deleted file mode 100644
index b9f3b94..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Root.hs
+++ /dev/null
@@ -1,199 +0,0 @@
1{- |
2Module : Numeric.GSL.Root
3Copyright : (c) Alberto Ruiz 2009
4License : GPL
5Maintainer : Alberto Ruiz
6Stability : provisional
7
8Multidimensional root finding.
9
10<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html>
11
12The example in the GSL manual:
13
14>>> let rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
15>>> let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
16>>> sol
17[1.0,1.0]
18>>> disp 3 path
1911x5
20 1.000 -10.000 -5.000 11.000 -1050.000
21 2.000 -3.976 24.827 4.976 90.203
22 3.000 -3.976 24.827 4.976 90.203
23 4.000 -3.976 24.827 4.976 90.203
24 5.000 -1.274 -5.680 2.274 -73.018
25 6.000 -1.274 -5.680 2.274 -73.018
26 7.000 0.249 0.298 0.751 2.359
27 8.000 0.249 0.298 0.751 2.359
28 9.000 1.000 0.878 -0.000 -1.218
2910.000 1.000 0.989 -0.000 -0.108
3011.000 1.000 1.000 0.000 0.000
31
32-}
33-----------------------------------------------------------------------------
34
35module Numeric.GSL.Root (
36 uniRoot, UniRootMethod(..),
37 uniRootJ, UniRootMethodJ(..),
38 root, RootMethod(..),
39 rootJ, RootMethodJ(..),
40) where
41
42import Data.Packed
43import Numeric.GSL.Internal
44import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
45import Foreign.C.Types
46import System.IO.Unsafe(unsafePerformIO)
47
48-------------------------------------------------------------------------
49type TVV = TV (TV Res)
50type TVM = TV (TM Res)
51
52
53data UniRootMethod = Bisection
54 | FalsePos
55 | Brent
56 deriving (Enum, Eq, Show, Bounded)
57
58uniRoot :: UniRootMethod
59 -> Double
60 -> Int
61 -> (Double -> Double)
62 -> Double
63 -> Double
64 -> (Double, Matrix Double)
65uniRoot method epsrel maxit fun xl xu = uniRootGen (fi (fromEnum method)) fun xl xu epsrel maxit
66
67uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
68 fp <- mkDoublefun f
69 rawpath <- createMIO maxit 4
70 (c_root m fp epsrel (fi maxit) xl xu)
71 "root"
72 let it = round (rawpath @@> (maxit-1,0))
73 path = takeRows it rawpath
74 [sol] = toLists $ dropRows (it-1) path
75 freeHaskellFunPtr fp
76 return (sol !! 1, path)
77
78foreign import ccall safe "root"
79 c_root:: CInt -> FunPtr (Double -> Double) -> Double -> CInt -> Double -> Double -> TM Res
80
81-------------------------------------------------------------------------
82data UniRootMethodJ = UNewton
83 | Secant
84 | Steffenson
85 deriving (Enum, Eq, Show, Bounded)
86
87uniRootJ :: UniRootMethodJ
88 -> Double
89 -> Int
90 -> (Double -> Double)
91 -> (Double -> Double)
92 -> Double
93 -> (Double, Matrix Double)
94uniRootJ method epsrel maxit fun dfun x = uniRootJGen (fi (fromEnum method)) fun
95 dfun x epsrel maxit
96
97uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do
98 fp <- mkDoublefun f
99 dfp <- mkDoublefun df
100 rawpath <- createMIO maxit 2
101 (c_rootj m fp dfp epsrel (fi maxit) x)
102 "rootj"
103 let it = round (rawpath @@> (maxit-1,0))
104 path = takeRows it rawpath
105 [sol] = toLists $ dropRows (it-1) path
106 freeHaskellFunPtr fp
107 return (sol !! 1, path)
108
109foreign import ccall safe "rootj"
110 c_rootj :: CInt -> FunPtr (Double -> Double) -> FunPtr (Double -> Double)
111 -> Double -> CInt -> Double -> TM Res
112
113-------------------------------------------------------------------------
114
115data RootMethod = Hybrids
116 | Hybrid
117 | DNewton
118 | Broyden
119 deriving (Enum,Eq,Show,Bounded)
120
121-- | Nonlinear multidimensional root finding using algorithms that do not require
122-- any derivative information to be supplied by the user.
123-- Any derivatives needed are approximated by finite differences.
124root :: RootMethod
125 -> Double -- ^ maximum residual
126 -> Int -- ^ maximum number of iterations allowed
127 -> ([Double] -> [Double]) -- ^ function to minimize
128 -> [Double] -- ^ starting point
129 -> ([Double], Matrix Double) -- ^ solution vector and optimization path
130
131root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit epsabs maxit
132
133rootGen m f xi epsabs maxit = unsafePerformIO $ do
134 let xiv = fromList xi
135 n = dim xiv
136 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
137 rawpath <- vec xiv $ \xiv' ->
138 createMIO maxit (2*n+1)
139 (c_multiroot m fp epsabs (fi maxit) // xiv')
140 "multiroot"
141 let it = round (rawpath @@> (maxit-1,0))
142 path = takeRows it rawpath
143 [sol] = toLists $ dropRows (it-1) path
144 freeHaskellFunPtr fp
145 return (take n $ drop 1 sol, path)
146
147
148foreign import ccall safe "multiroot"
149 c_multiroot:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
150
151-------------------------------------------------------------------------
152
153data RootMethodJ = HybridsJ
154 | HybridJ
155 | Newton
156 | GNewton
157 deriving (Enum,Eq,Show,Bounded)
158
159-- | Nonlinear multidimensional root finding using both the function and its derivatives.
160rootJ :: RootMethodJ
161 -> Double -- ^ maximum residual
162 -> Int -- ^ maximum number of iterations allowed
163 -> ([Double] -> [Double]) -- ^ function to minimize
164 -> ([Double] -> [[Double]]) -- ^ Jacobian
165 -> [Double] -- ^ starting point
166 -> ([Double], Matrix Double) -- ^ solution vector and optimization path
167
168rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun jac xinit epsabs maxit
169
170rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
171 let xiv = fromList xi
172 n = dim xiv
173 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
174 jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
175 rawpath <- vec xiv $ \xiv' ->
176 createMIO maxit (2*n+1)
177 (c_multirootj m fp jp epsabs (fi maxit) // xiv')
178 "multiroot"
179 let it = round (rawpath @@> (maxit-1,0))
180 path = takeRows it rawpath
181 [sol] = toLists $ dropRows (it-1) path
182 freeHaskellFunPtr fp
183 freeHaskellFunPtr jp
184 return (take n $ drop 1 sol, path)
185
186foreign import ccall safe "multirootj"
187 c_multirootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
188
189-------------------------------------------------------
190
191checkdim1 n v
192 | dim v == n = v
193 | otherwise = error $ "Error: "++ show n
194 ++ " components expected in the result of the function supplied to root"
195
196checkdim2 n m
197 | rows m == n && cols m == n = m
198 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
199 ++ " Jacobian expected in rootJ"
diff --git a/packages/hmatrix/src/Numeric/GSL/Vector.hs b/packages/hmatrix/src/Numeric/GSL/Vector.hs
deleted file mode 100644
index af79f32..0000000
--- a/packages/hmatrix/src/Numeric/GSL/Vector.hs
+++ /dev/null
@@ -1,107 +0,0 @@
1-----------------------------------------------------------------------------
2-- |
3-- Module : Numeric.GSL.Vector
4-- Copyright : (c) Alberto Ruiz 2007-14
5-- License : GPL
6-- Maintainer : Alberto Ruiz
7-- Stability : provisional
8--
9-----------------------------------------------------------------------------
10
11module Numeric.GSL.Vector (
12 randomVector,
13 saveMatrix,
14 fwriteVector, freadVector, fprintfVector, fscanfVector
15) where
16
17import Data.Packed
18import Numeric.LinearAlgebra(RandDist(..))
19import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
20
21import Foreign.Marshal.Alloc(free)
22import Foreign.Ptr(Ptr)
23import Foreign.C.Types
24import Foreign.C.String(newCString)
25import System.IO.Unsafe(unsafePerformIO)
26
27fromei x = fromIntegral (fromEnum x) :: CInt
28
29-----------------------------------------------------------------------
30
31-- | Obtains a vector of pseudorandom elements from the the mt19937 generator in GSL, with a given seed. Use randomIO to get a random seed.
32randomVector :: Int -- ^ seed
33 -> RandDist -- ^ distribution
34 -> Int -- ^ vector size
35 -> Vector Double
36randomVector seed dist n = unsafePerformIO $ do
37 r <- createVector n
38 app1 (c_random_vector_GSL (fi seed) ((fi.fromEnum) dist)) vec r "randomVectorGSL"
39 return r
40
41foreign import ccall unsafe "random_vector_GSL" c_random_vector_GSL :: CInt -> CInt -> TV
42
43--------------------------------------------------------------------------------
44
45-- | Saves a matrix as 2D ASCII table.
46saveMatrix :: FilePath
47 -> String -- ^ format (%f, %g, %e)
48 -> Matrix Double
49 -> IO ()
50saveMatrix filename fmt m = do
51 charname <- newCString filename
52 charfmt <- newCString fmt
53 let o = if orderOf m == RowMajor then 1 else 0
54 app1 (matrix_fprintf charname charfmt o) mat m "matrix_fprintf"
55 free charname
56 free charfmt
57
58foreign import ccall unsafe "matrix_fprintf" matrix_fprintf :: Ptr CChar -> Ptr CChar -> CInt -> TM
59
60--------------------------------------------------------------------------------
61
62-- | Loads a vector from an ASCII file (the number of elements must be known in advance).
63fscanfVector :: FilePath -> Int -> IO (Vector Double)
64fscanfVector filename n = do
65 charname <- newCString filename
66 res <- createVector n
67 app1 (gsl_vector_fscanf charname) vec res "gsl_vector_fscanf"
68 free charname
69 return res
70
71foreign import ccall unsafe "vector_fscanf" gsl_vector_fscanf:: Ptr CChar -> TV
72
73-- | Saves the elements of a vector, with a given format (%f, %e, %g), to an ASCII file.
74fprintfVector :: FilePath -> String -> Vector Double -> IO ()
75fprintfVector filename fmt v = do
76 charname <- newCString filename
77 charfmt <- newCString fmt
78 app1 (gsl_vector_fprintf charname charfmt) vec v "gsl_vector_fprintf"
79 free charname
80 free charfmt
81
82foreign import ccall unsafe "vector_fprintf" gsl_vector_fprintf :: Ptr CChar -> Ptr CChar -> TV
83
84-- | Loads a vector from a binary file (the number of elements must be known in advance).
85freadVector :: FilePath -> Int -> IO (Vector Double)
86freadVector filename n = do
87 charname <- newCString filename
88 res <- createVector n
89 app1 (gsl_vector_fread charname) vec res "gsl_vector_fread"
90 free charname
91 return res
92
93foreign import ccall unsafe "vector_fread" gsl_vector_fread:: Ptr CChar -> TV
94
95-- | Saves the elements of a vector to a binary file.
96fwriteVector :: FilePath -> Vector Double -> IO ()
97fwriteVector filename v = do
98 charname <- newCString filename
99 app1 (gsl_vector_fwrite charname) vec v "gsl_vector_fwrite"
100 free charname
101
102foreign import ccall unsafe "vector_fwrite" gsl_vector_fwrite :: Ptr CChar -> TV
103
104type PD = Ptr Double --
105type TV = CInt -> PD -> IO CInt --
106type TM = CInt -> CInt -> PD -> IO CInt --
107
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c b/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
deleted file mode 100644
index ffc5c20..0000000
--- a/packages/hmatrix/src/Numeric/GSL/gsl-aux.c
+++ /dev/null
@@ -1,945 +0,0 @@
1#include <gsl/gsl_complex.h>
2
3#define RVEC(A) int A##n, double*A##p
4#define RMAT(A) int A##r, int A##c, double* A##p
5#define KRVEC(A) int A##n, const double*A##p
6#define KRMAT(A) int A##r, int A##c, const double* A##p
7
8#define CVEC(A) int A##n, gsl_complex*A##p
9#define CMAT(A) int A##r, int A##c, gsl_complex* A##p
10#define KCVEC(A) int A##n, const gsl_complex*A##p
11#define KCMAT(A) int A##r, int A##c, const gsl_complex* A##p
12
13#define FVEC(A) int A##n, float*A##p
14#define FMAT(A) int A##r, int A##c, float* A##p
15#define KFVEC(A) int A##n, const float*A##p
16#define KFMAT(A) int A##r, int A##c, const float* A##p
17
18#define QVEC(A) int A##n, gsl_complex_float*A##p
19#define QMAT(A) int A##r, int A##c, gsl_complex_float* A##p
20#define KQVEC(A) int A##n, const gsl_complex_float*A##p
21#define KQMAT(A) int A##r, int A##c, const gsl_complex_float* A##p
22
23#include <gsl/gsl_blas.h>
24#include <gsl/gsl_math.h>
25#include <gsl/gsl_errno.h>
26#include <gsl/gsl_fft_complex.h>
27#include <gsl/gsl_integration.h>
28#include <gsl/gsl_deriv.h>
29#include <gsl/gsl_poly.h>
30#include <gsl/gsl_multimin.h>
31#include <gsl/gsl_multiroots.h>
32#include <gsl/gsl_min.h>
33#include <gsl/gsl_complex_math.h>
34#include <gsl/gsl_rng.h>
35#include <gsl/gsl_randist.h>
36#include <gsl/gsl_roots.h>
37#include <gsl/gsl_multifit_nlin.h>
38#include <string.h>
39#include <stdio.h>
40
41#define MACRO(B) do {B} while (0)
42#define ERROR(CODE) MACRO(return CODE;)
43#define REQUIRES(COND, CODE) MACRO(if(!(COND)) {ERROR(CODE);})
44#define OK return 0;
45
46#define MIN(A,B) ((A)<(B)?(A):(B))
47#define MAX(A,B) ((A)>(B)?(A):(B))
48
49#ifdef DBG
50#define DEBUGMSG(M) printf("*** calling aux C function: %s\n",M);
51#else
52#define DEBUGMSG(M)
53#endif
54
55#define CHECK(RES,CODE) MACRO(if(RES) return CODE;)
56
57#ifdef DBG
58#define DEBUGMAT(MSG,X) printf(MSG" = \n"); gsl_matrix_fprintf(stdout,X,"%f"); printf("\n");
59#else
60#define DEBUGMAT(MSG,X)
61#endif
62
63#ifdef DBG
64#define DEBUGVEC(MSG,X) printf(MSG" = \n"); gsl_vector_fprintf(stdout,X,"%f"); printf("\n");
65#else
66#define DEBUGVEC(MSG,X)
67#endif
68
69#define DVVIEW(A) gsl_vector_view A = gsl_vector_view_array(A##p,A##n)
70#define DMVIEW(A) gsl_matrix_view A = gsl_matrix_view_array(A##p,A##r,A##c)
71#define CVVIEW(A) gsl_vector_complex_view A = gsl_vector_complex_view_array((double*)A##p,A##n)
72#define CMVIEW(A) gsl_matrix_complex_view A = gsl_matrix_complex_view_array((double*)A##p,A##r,A##c)
73#define KDVVIEW(A) gsl_vector_const_view A = gsl_vector_const_view_array(A##p,A##n)
74#define KDMVIEW(A) gsl_matrix_const_view A = gsl_matrix_const_view_array(A##p,A##r,A##c)
75#define KCVVIEW(A) gsl_vector_complex_const_view A = gsl_vector_complex_const_view_array((double*)A##p,A##n)
76#define KCMVIEW(A) gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array((double*)A##p,A##r,A##c)
77
78#define FVVIEW(A) gsl_vector_float_view A = gsl_vector_float_view_array(A##p,A##n)
79#define FMVIEW(A) gsl_matrix_float_view A = gsl_matrix_float_view_array(A##p,A##r,A##c)
80#define QVVIEW(A) gsl_vector_complex_float_view A = gsl_vector_float_complex_view_array((float*)A##p,A##n)
81#define QMVIEW(A) gsl_matrix_complex_float_view A = gsl_matrix_float_complex_view_array((float*)A##p,A##r,A##c)
82#define KFVVIEW(A) gsl_vector_float_const_view A = gsl_vector_float_const_view_array(A##p,A##n)
83#define KFMVIEW(A) gsl_matrix_float_const_view A = gsl_matrix_float_const_view_array(A##p,A##r,A##c)
84#define KQVVIEW(A) gsl_vector_complex_float_const_view A = gsl_vector_complex_float_const_view_array((float*)A##p,A##n)
85#define KQMVIEW(A) gsl_matrix_complex_float_const_view A = gsl_matrix_complex_float_const_view_array((float*)A##p,A##r,A##c)
86
87#define V(a) (&a.vector)
88#define M(a) (&a.matrix)
89
90#define GCVEC(A) int A##n, gsl_complex*A##p
91#define KGCVEC(A) int A##n, const gsl_complex*A##p
92
93#define GQVEC(A) int A##n, gsl_complex_float*A##p
94#define KGQVEC(A) int A##n, const gsl_complex_float*A##p
95
96#define BAD_SIZE 2000
97#define BAD_CODE 2001
98#define MEM 2002
99#define BAD_FILE 2003
100
101
102void no_abort_on_error() {
103 gsl_set_error_handler_off();
104}
105
106
107
108int fft(int code, KCVEC(X), CVEC(R)) {
109 REQUIRES(Xn == Rn,BAD_SIZE);
110 DEBUGMSG("fft");
111 int s = Xn;
112 gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc (s);
113 gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc (s);
114 gsl_vector_const_view X = gsl_vector_const_view_array((double*)Xp, 2*Xn);
115 gsl_vector_view R = gsl_vector_view_array((double*)Rp, 2*Rn);
116 gsl_blas_dcopy(&X.vector,&R.vector);
117 if(code==0) {
118 gsl_fft_complex_forward ((double*)Rp, 1, s, wavetable, workspace);
119 } else {
120 gsl_fft_complex_inverse ((double*)Rp, 1, s, wavetable, workspace);
121 }
122 gsl_fft_complex_wavetable_free (wavetable);
123 gsl_fft_complex_workspace_free (workspace);
124 OK
125}
126
127
128int deriv(int code, double f(double, void*), double x, double h, double * result, double * abserr)
129{
130 gsl_function F;
131 F.function = f;
132 F.params = 0;
133
134 if(code==0) return gsl_deriv_central (&F, x, h, result, abserr);
135
136 if(code==1) return gsl_deriv_forward (&F, x, h, result, abserr);
137
138 if(code==2) return gsl_deriv_backward (&F, x, h, result, abserr);
139
140 return 0;
141}
142
143
144int integrate_qng(double f(double, void*), double a, double b, double aprec, double prec,
145 double *result, double*error) {
146 DEBUGMSG("integrate_qng");
147 gsl_function F;
148 F.function = f;
149 F.params = NULL;
150 size_t neval;
151 int res = gsl_integration_qng (&F, a,b, aprec, prec, result, error, &neval);
152 CHECK(res,res);
153 OK
154}
155
156int integrate_qags(double f(double,void*), double a, double b, double aprec, double prec, int w,
157 double *result, double* error) {
158 DEBUGMSG("integrate_qags");
159 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
160 gsl_function F;
161 F.function = f;
162 F.params = NULL;
163 int res = gsl_integration_qags (&F, a,b, aprec, prec, w,wk, result, error);
164 CHECK(res,res);
165 gsl_integration_workspace_free (wk);
166 OK
167}
168
169int integrate_qagi(double f(double,void*), double aprec, double prec, int w,
170 double *result, double* error) {
171 DEBUGMSG("integrate_qagi");
172 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
173 gsl_function F;
174 F.function = f;
175 F.params = NULL;
176 int res = gsl_integration_qagi (&F, aprec, prec, w,wk, result, error);
177 CHECK(res,res);
178 gsl_integration_workspace_free (wk);
179 OK
180}
181
182
183int integrate_qagiu(double f(double,void*), double a, double aprec, double prec, int w,
184 double *result, double* error) {
185 DEBUGMSG("integrate_qagiu");
186 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
187 gsl_function F;
188 F.function = f;
189 F.params = NULL;
190 int res = gsl_integration_qagiu (&F, a, aprec, prec, w,wk, result, error);
191 CHECK(res,res);
192 gsl_integration_workspace_free (wk);
193 OK
194}
195
196
197int integrate_qagil(double f(double,void*), double b, double aprec, double prec, int w,
198 double *result, double* error) {
199 DEBUGMSG("integrate_qagil");
200 gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
201 gsl_function F;
202 F.function = f;
203 F.params = NULL;
204 int res = gsl_integration_qagil (&F, b, aprec, prec, w,wk, result, error);
205 CHECK(res,res);
206 gsl_integration_workspace_free (wk);
207 OK
208}
209
210int integrate_cquad(double f(double,void*), double a, double b, double aprec, double prec,
211 int w, double *result, double* error, int *neval) {
212 DEBUGMSG("integrate_cquad");
213 gsl_integration_cquad_workspace * wk = gsl_integration_cquad_workspace_alloc (w);
214 gsl_function F;
215 F.function = f;
216 F.params = NULL;
217 size_t * sneval = NULL;
218 int res = gsl_integration_cquad (&F, a, b, aprec, prec, wk, result, error, sneval);
219 *neval = *sneval;
220 CHECK(res,res);
221 gsl_integration_cquad_workspace_free (wk);
222 OK
223}
224
225
226int polySolve(KRVEC(a), CVEC(z)) {
227 DEBUGMSG("polySolve");
228 REQUIRES(an>1,BAD_SIZE);
229 gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (an);
230 int res = gsl_poly_complex_solve ((double*)ap, an, w, (double*)zp);
231 CHECK(res,res);
232 gsl_poly_complex_workspace_free (w);
233 OK;
234}
235
236int vector_fscanf(char*filename, RVEC(a)) {
237 DEBUGMSG("gsl_vector_fscanf");
238 DVVIEW(a);
239 FILE * f = fopen(filename,"r");
240 CHECK(!f,BAD_FILE);
241 int res = gsl_vector_fscanf(f,V(a));
242 CHECK(res,res);
243 fclose (f);
244 OK
245}
246
247int vector_fprintf(char*filename, char*fmt, RVEC(a)) {
248 DEBUGMSG("gsl_vector_fprintf");
249 DVVIEW(a);
250 FILE * f = fopen(filename,"w");
251 CHECK(!f,BAD_FILE);
252 int res = gsl_vector_fprintf(f,V(a),fmt);
253 CHECK(res,res);
254 fclose (f);
255 OK
256}
257
258int vector_fread(char*filename, RVEC(a)) {
259 DEBUGMSG("gsl_vector_fread");
260 DVVIEW(a);
261 FILE * f = fopen(filename,"r");
262 CHECK(!f,BAD_FILE);
263 int res = gsl_vector_fread(f,V(a));
264 CHECK(res,res);
265 fclose (f);
266 OK
267}
268
269int vector_fwrite(char*filename, RVEC(a)) {
270 DEBUGMSG("gsl_vector_fwrite");
271 DVVIEW(a);
272 FILE * f = fopen(filename,"w");
273 CHECK(!f,BAD_FILE);
274 int res = gsl_vector_fwrite(f,V(a));
275 CHECK(res,res);
276 fclose (f);
277 OK
278}
279
280int matrix_fprintf(char*filename, char*fmt, int ro, RMAT(m)) {
281 DEBUGMSG("matrix_fprintf");
282 FILE * f = fopen(filename,"w");
283 CHECK(!f,BAD_FILE);
284 int i,j,sr,sc;
285 if (ro==1) { sr = mc; sc = 1;} else { sr = 1; sc = mr;}
286 #define AT(M,r,c) (M##p[(r)*sr+(c)*sc])
287 for (i=0; i<mr; i++) {
288 for (j=0; j<mc-1; j++) {
289 fprintf(f,fmt,AT(m,i,j));
290 fprintf(f," ");
291 }
292 fprintf(f,fmt,AT(m,i,j));
293 fprintf(f,"\n");
294 }
295 fclose (f);
296 OK
297}
298
299//---------------------------------------------------------------
300
301typedef double Trawfun(int, double*);
302
303double only_f_aux_min(const gsl_vector*x, void *pars) {
304 Trawfun * f = (Trawfun*) pars;
305 double* p = (double*)calloc(x->size,sizeof(double));
306 int k;
307 for(k=0;k<x->size;k++) {
308 p[k] = gsl_vector_get(x,k);
309 }
310 double res = f(x->size,p);
311 free(p);
312 return res;
313}
314
315double only_f_aux_root(double x, void *pars);
316int uniMinimize(int method, double f(double),
317 double epsrel, int maxit, double min,
318 double xl, double xu, RMAT(sol)) {
319 REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
320 DEBUGMSG("minimize_only_f");
321 gsl_function my_func;
322 my_func.function = only_f_aux_root;
323 my_func.params = f;
324 size_t iter = 0;
325 int status;
326 const gsl_min_fminimizer_type *T;
327 gsl_min_fminimizer *s;
328 // Starting point
329 switch(method) {
330 case 0 : {T = gsl_min_fminimizer_goldensection; break; }
331 case 1 : {T = gsl_min_fminimizer_brent; break; }
332 case 2 : {T = gsl_min_fminimizer_quad_golden; break; }
333 default: ERROR(BAD_CODE);
334 }
335 s = gsl_min_fminimizer_alloc (T);
336 gsl_min_fminimizer_set (s, &my_func, min, xl, xu);
337 do {
338 double current_min, current_lo, current_hi;
339 status = gsl_min_fminimizer_iterate (s);
340 current_min = gsl_min_fminimizer_x_minimum (s);
341 current_lo = gsl_min_fminimizer_x_lower (s);
342 current_hi = gsl_min_fminimizer_x_upper (s);
343 solp[iter*solc] = iter + 1;
344 solp[iter*solc+1] = current_min;
345 solp[iter*solc+2] = current_lo;
346 solp[iter*solc+3] = current_hi;
347 iter++;
348 if (status) /* check if solver is stuck */
349 break;
350
351 status =
352 gsl_min_test_interval (current_lo, current_hi, 0, epsrel);
353 }
354 while (status == GSL_CONTINUE && iter < maxit);
355 int i;
356 for (i=iter; i<solr; i++) {
357 solp[i*solc+0] = iter;
358 solp[i*solc+1]=0.;
359 solp[i*solc+2]=0.;
360 solp[i*solc+3]=0.;
361 }
362 gsl_min_fminimizer_free(s);
363 OK
364}
365
366
367
368// this version returns info about intermediate steps
369int minimize(int method, double f(int, double*), double tolsize, int maxit,
370 KRVEC(xi), KRVEC(sz), RMAT(sol)) {
371 REQUIRES(xin==szn && solr == maxit && solc == 3+xin,BAD_SIZE);
372 DEBUGMSG("minimizeList (nmsimplex)");
373 gsl_multimin_function my_func;
374 // extract function from pars
375 my_func.f = only_f_aux_min;
376 my_func.n = xin;
377 my_func.params = f;
378 size_t iter = 0;
379 int status;
380 double size;
381 const gsl_multimin_fminimizer_type *T;
382 gsl_multimin_fminimizer *s = NULL;
383 // Initial vertex size vector
384 KDVVIEW(sz);
385 // Starting point
386 KDVVIEW(xi);
387 // Minimizer nmsimplex, without derivatives
388 switch(method) {
389 case 0 : {T = gsl_multimin_fminimizer_nmsimplex; break; }
390#ifdef GSL110
391 case 1 : {T = gsl_multimin_fminimizer_nmsimplex; break; }
392#else
393 case 1 : {T = gsl_multimin_fminimizer_nmsimplex2; break; }
394#endif
395 default: ERROR(BAD_CODE);
396 }
397 s = gsl_multimin_fminimizer_alloc (T, my_func.n);
398 gsl_multimin_fminimizer_set (s, &my_func, V(xi), V(sz));
399 do {
400 status = gsl_multimin_fminimizer_iterate (s);
401 size = gsl_multimin_fminimizer_size (s);
402
403 solp[iter*solc+0] = iter+1;
404 solp[iter*solc+1] = s->fval;
405 solp[iter*solc+2] = size;
406
407 int k;
408 for(k=0;k<xin;k++) {
409 solp[iter*solc+k+3] = gsl_vector_get(s->x,k);
410 }
411 iter++;
412 if (status) break;
413 status = gsl_multimin_test_size (size, tolsize);
414 } while (status == GSL_CONTINUE && iter < maxit);
415 int i,j;
416 for (i=iter; i<solr; i++) {
417 solp[i*solc+0] = iter;
418 for(j=1;j<solc;j++) {
419 solp[i*solc+j]=0.;
420 }
421 }
422 gsl_multimin_fminimizer_free(s);
423 OK
424}
425
426// working with the gradient
427
428typedef struct {double (*f)(int, double*); int (*df)(int, double*, int, double*);} Tfdf;
429
430double f_aux_min(const gsl_vector*x, void *pars) {
431 Tfdf * fdf = ((Tfdf*) pars);
432 double* p = (double*)calloc(x->size,sizeof(double));
433 int k;
434 for(k=0;k<x->size;k++) {
435 p[k] = gsl_vector_get(x,k);
436 }
437 double res = fdf->f(x->size,p);
438 free(p);
439 return res;
440}
441
442
443void df_aux_min(const gsl_vector * x, void * pars, gsl_vector * g) {
444 Tfdf * fdf = ((Tfdf*) pars);
445 double* p = (double*)calloc(x->size,sizeof(double));
446 double* q = (double*)calloc(g->size,sizeof(double));
447 int k;
448 for(k=0;k<x->size;k++) {
449 p[k] = gsl_vector_get(x,k);
450 }
451
452 fdf->df(x->size,p,g->size,q);
453
454 for(k=0;k<x->size;k++) {
455 gsl_vector_set(g,k,q[k]);
456 }
457 free(p);
458 free(q);
459}
460
461void fdf_aux_min(const gsl_vector * x, void * pars, double * f, gsl_vector * g) {
462 *f = f_aux_min(x,pars);
463 df_aux_min(x,pars,g);
464}
465
466
467int minimizeD(int method, double f(int, double*), int df(int, double*, int, double*),
468 double initstep, double minimpar, double tolgrad, int maxit,
469 KRVEC(xi), RMAT(sol)) {
470 REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
471 DEBUGMSG("minimizeWithDeriv (conjugate_fr)");
472 gsl_multimin_function_fdf my_func;
473 // extract function from pars
474 my_func.f = f_aux_min;
475 my_func.df = df_aux_min;
476 my_func.fdf = fdf_aux_min;
477 my_func.n = xin;
478 Tfdf stfdf;
479 stfdf.f = f;
480 stfdf.df = df;
481 my_func.params = &stfdf;
482 size_t iter = 0;
483 int status;
484 const gsl_multimin_fdfminimizer_type *T;
485 gsl_multimin_fdfminimizer *s = NULL;
486 // Starting point
487 KDVVIEW(xi);
488 // conjugate gradient fr
489 switch(method) {
490 case 0 : {T = gsl_multimin_fdfminimizer_conjugate_fr; break; }
491 case 1 : {T = gsl_multimin_fdfminimizer_conjugate_pr; break; }
492 case 2 : {T = gsl_multimin_fdfminimizer_vector_bfgs; break; }
493 case 3 : {T = gsl_multimin_fdfminimizer_vector_bfgs2; break; }
494 case 4 : {T = gsl_multimin_fdfminimizer_steepest_descent; break; }
495 default: ERROR(BAD_CODE);
496 }
497 s = gsl_multimin_fdfminimizer_alloc (T, my_func.n);
498 gsl_multimin_fdfminimizer_set (s, &my_func, V(xi), initstep, minimpar);
499 do {
500 status = gsl_multimin_fdfminimizer_iterate (s);
501 solp[iter*solc+0] = iter+1;
502 solp[iter*solc+1] = s->f;
503 int k;
504 for(k=0;k<xin;k++) {
505 solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
506 }
507 iter++;
508 if (status) break;
509 status = gsl_multimin_test_gradient (s->gradient, tolgrad);
510 } while (status == GSL_CONTINUE && iter < maxit);
511 int i,j;
512 for (i=iter; i<solr; i++) {
513 solp[i*solc+0] = iter;
514 for(j=1;j<solc;j++) {
515 solp[i*solc+j]=0.;
516 }
517 }
518 gsl_multimin_fdfminimizer_free(s);
519 OK
520}
521
522//---------------------------------------------------------------
523
524double only_f_aux_root(double x, void *pars) {
525 double (*f)(double) = (double (*)(double)) pars;
526 return f(x);
527}
528
529int root(int method, double f(double),
530 double epsrel, int maxit,
531 double xl, double xu, RMAT(sol)) {
532 REQUIRES(solr == maxit && solc == 4,BAD_SIZE);
533 DEBUGMSG("root_only_f");
534 gsl_function my_func;
535 // extract function from pars
536 my_func.function = only_f_aux_root;
537 my_func.params = f;
538 size_t iter = 0;
539 int status;
540 const gsl_root_fsolver_type *T;
541 gsl_root_fsolver *s;
542 // Starting point
543 switch(method) {
544 case 0 : {T = gsl_root_fsolver_bisection; printf("7\n"); break; }
545 case 1 : {T = gsl_root_fsolver_falsepos; break; }
546 case 2 : {T = gsl_root_fsolver_brent; break; }
547 default: ERROR(BAD_CODE);
548 }
549 s = gsl_root_fsolver_alloc (T);
550 gsl_root_fsolver_set (s, &my_func, xl, xu);
551 do {
552 double best, current_lo, current_hi;
553 status = gsl_root_fsolver_iterate (s);
554 best = gsl_root_fsolver_root (s);
555 current_lo = gsl_root_fsolver_x_lower (s);
556 current_hi = gsl_root_fsolver_x_upper (s);
557 solp[iter*solc] = iter + 1;
558 solp[iter*solc+1] = best;
559 solp[iter*solc+2] = current_lo;
560 solp[iter*solc+3] = current_hi;
561 iter++;
562 if (status) /* check if solver is stuck */
563 break;
564
565 status =
566 gsl_root_test_interval (current_lo, current_hi, 0, epsrel);
567 }
568 while (status == GSL_CONTINUE && iter < maxit);
569 int i;
570 for (i=iter; i<solr; i++) {
571 solp[i*solc+0] = iter;
572 solp[i*solc+1]=0.;
573 solp[i*solc+2]=0.;
574 solp[i*solc+3]=0.;
575 }
576 gsl_root_fsolver_free(s);
577 OK
578}
579
580typedef struct {
581 double (*f)(double);
582 double (*jf)(double);
583} uniTfjf;
584
585double f_aux_uni(double x, void *pars) {
586 uniTfjf * fjf = ((uniTfjf*) pars);
587 return (fjf->f)(x);
588}
589
590double jf_aux_uni(double x, void * pars) {
591 uniTfjf * fjf = ((uniTfjf*) pars);
592 return (fjf->jf)(x);
593}
594
595void fjf_aux_uni(double x, void * pars, double * f, double * g) {
596 *f = f_aux_uni(x,pars);
597 *g = jf_aux_uni(x,pars);
598}
599
600int rootj(int method, double f(double),
601 double df(double),
602 double epsrel, int maxit,
603 double x, RMAT(sol)) {
604 REQUIRES(solr == maxit && solc == 2,BAD_SIZE);
605 DEBUGMSG("root_fjf");
606 gsl_function_fdf my_func;
607 // extract function from pars
608 my_func.f = f_aux_uni;
609 my_func.df = jf_aux_uni;
610 my_func.fdf = fjf_aux_uni;
611 uniTfjf stfjf;
612 stfjf.f = f;
613 stfjf.jf = df;
614 my_func.params = &stfjf;
615 size_t iter = 0;
616 int status;
617 const gsl_root_fdfsolver_type *T;
618 gsl_root_fdfsolver *s;
619 // Starting point
620 switch(method) {
621 case 0 : {T = gsl_root_fdfsolver_newton;; break; }
622 case 1 : {T = gsl_root_fdfsolver_secant; break; }
623 case 2 : {T = gsl_root_fdfsolver_steffenson; break; }
624 default: ERROR(BAD_CODE);
625 }
626 s = gsl_root_fdfsolver_alloc (T);
627
628 gsl_root_fdfsolver_set (s, &my_func, x);
629
630 do {
631 double x0;
632 status = gsl_root_fdfsolver_iterate (s);
633 x0 = x;
634 x = gsl_root_fdfsolver_root(s);
635 solp[iter*solc+0] = iter+1;
636 solp[iter*solc+1] = x;
637
638 iter++;
639 if (status) /* check if solver is stuck */
640 break;
641
642 status =
643 gsl_root_test_delta (x, x0, 0, epsrel);
644 }
645 while (status == GSL_CONTINUE && iter < maxit);
646
647 int i;
648 for (i=iter; i<solr; i++) {
649 solp[i*solc+0] = iter;
650 solp[i*solc+1]=0.;
651 }
652 gsl_root_fdfsolver_free(s);
653 OK
654}
655
656
657//---------------------------------------------------------------
658
659typedef void TrawfunV(int, double*, int, double*);
660
661int only_f_aux_multiroot(const gsl_vector*x, void *pars, gsl_vector*y) {
662 TrawfunV * f = (TrawfunV*) pars;
663 double* p = (double*)calloc(x->size,sizeof(double));
664 double* q = (double*)calloc(y->size,sizeof(double));
665 int k;
666 for(k=0;k<x->size;k++) {
667 p[k] = gsl_vector_get(x,k);
668 }
669 f(x->size,p,y->size,q);
670 for(k=0;k<y->size;k++) {
671 gsl_vector_set(y,k,q[k]);
672 }
673 free(p);
674 free(q);
675 return 0; //hmmm
676}
677
678int multiroot(int method, void f(int, double*, int, double*),
679 double epsabs, int maxit,
680 KRVEC(xi), RMAT(sol)) {
681 REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
682 DEBUGMSG("root_only_f");
683 gsl_multiroot_function my_func;
684 // extract function from pars
685 my_func.f = only_f_aux_multiroot;
686 my_func.n = xin;
687 my_func.params = f;
688 size_t iter = 0;
689 int status;
690 const gsl_multiroot_fsolver_type *T;
691 gsl_multiroot_fsolver *s;
692 // Starting point
693 KDVVIEW(xi);
694 switch(method) {
695 case 0 : {T = gsl_multiroot_fsolver_hybrids;; break; }
696 case 1 : {T = gsl_multiroot_fsolver_hybrid; break; }
697 case 2 : {T = gsl_multiroot_fsolver_dnewton; break; }
698 case 3 : {T = gsl_multiroot_fsolver_broyden; break; }
699 default: ERROR(BAD_CODE);
700 }
701 s = gsl_multiroot_fsolver_alloc (T, my_func.n);
702 gsl_multiroot_fsolver_set (s, &my_func, V(xi));
703
704 do {
705 status = gsl_multiroot_fsolver_iterate (s);
706
707 solp[iter*solc+0] = iter+1;
708
709 int k;
710 for(k=0;k<xin;k++) {
711 solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
712 }
713 for(k=xin;k<2*xin;k++) {
714 solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
715 }
716
717 iter++;
718 if (status) /* check if solver is stuck */
719 break;
720
721 status =
722 gsl_multiroot_test_residual (s->f, epsabs);
723 }
724 while (status == GSL_CONTINUE && iter < maxit);
725
726 int i,j;
727 for (i=iter; i<solr; i++) {
728 solp[i*solc+0] = iter;
729 for(j=1;j<solc;j++) {
730 solp[i*solc+j]=0.;
731 }
732 }
733 gsl_multiroot_fsolver_free(s);
734 OK
735}
736
737// working with the jacobian
738
739typedef struct {int (*f)(int, double*, int, double *);
740 int (*jf)(int, double*, int, int, double*);} Tfjf;
741
742int f_aux(const gsl_vector*x, void *pars, gsl_vector*y) {
743 Tfjf * fjf = ((Tfjf*) pars);
744 double* p = (double*)calloc(x->size,sizeof(double));
745 double* q = (double*)calloc(y->size,sizeof(double));
746 int k;
747 for(k=0;k<x->size;k++) {
748 p[k] = gsl_vector_get(x,k);
749 }
750 (fjf->f)(x->size,p,y->size,q);
751 for(k=0;k<y->size;k++) {
752 gsl_vector_set(y,k,q[k]);
753 }
754 free(p);
755 free(q);
756 return 0;
757}
758
759int jf_aux(const gsl_vector * x, void * pars, gsl_matrix * jac) {
760 Tfjf * fjf = ((Tfjf*) pars);
761 double* p = (double*)calloc(x->size,sizeof(double));
762 double* q = (double*)calloc((jac->size1)*(jac->size2),sizeof(double));
763 int i,j,k;
764 for(k=0;k<x->size;k++) {
765 p[k] = gsl_vector_get(x,k);
766 }
767
768 (fjf->jf)(x->size,p,jac->size1,jac->size2,q);
769
770 k=0;
771 for(i=0;i<jac->size1;i++) {
772 for(j=0;j<jac->size2;j++){
773 gsl_matrix_set(jac,i,j,q[k++]);
774 }
775 }
776 free(p);
777 free(q);
778 return 0;
779}
780
781int fjf_aux(const gsl_vector * x, void * pars, gsl_vector * f, gsl_matrix * g) {
782 f_aux(x,pars,f);
783 jf_aux(x,pars,g);
784 return 0;
785}
786
787int multirootj(int method, int f(int, double*, int, double*),
788 int jac(int, double*, int, int, double*),
789 double epsabs, int maxit,
790 KRVEC(xi), RMAT(sol)) {
791 REQUIRES(solr == maxit && solc == 1+2*xin,BAD_SIZE);
792 DEBUGMSG("root_fjf");
793 gsl_multiroot_function_fdf my_func;
794 // extract function from pars
795 my_func.f = f_aux;
796 my_func.df = jf_aux;
797 my_func.fdf = fjf_aux;
798 my_func.n = xin;
799 Tfjf stfjf;
800 stfjf.f = f;
801 stfjf.jf = jac;
802 my_func.params = &stfjf;
803 size_t iter = 0;
804 int status;
805 const gsl_multiroot_fdfsolver_type *T;
806 gsl_multiroot_fdfsolver *s;
807 // Starting point
808 KDVVIEW(xi);
809 switch(method) {
810 case 0 : {T = gsl_multiroot_fdfsolver_hybridsj;; break; }
811 case 1 : {T = gsl_multiroot_fdfsolver_hybridj; break; }
812 case 2 : {T = gsl_multiroot_fdfsolver_newton; break; }
813 case 3 : {T = gsl_multiroot_fdfsolver_gnewton; break; }
814 default: ERROR(BAD_CODE);
815 }
816 s = gsl_multiroot_fdfsolver_alloc (T, my_func.n);
817
818 gsl_multiroot_fdfsolver_set (s, &my_func, V(xi));
819
820 do {
821 status = gsl_multiroot_fdfsolver_iterate (s);
822
823 solp[iter*solc+0] = iter+1;
824
825 int k;
826 for(k=0;k<xin;k++) {
827 solp[iter*solc+k+1] = gsl_vector_get(s->x,k);
828 }
829 for(k=xin;k<2*xin;k++) {
830 solp[iter*solc+k+1] = gsl_vector_get(s->f,k-xin);
831 }
832
833 iter++;
834 if (status) /* check if solver is stuck */
835 break;
836
837 status =
838 gsl_multiroot_test_residual (s->f, epsabs);
839 }
840 while (status == GSL_CONTINUE && iter < maxit);
841
842 int i,j;
843 for (i=iter; i<solr; i++) {
844 solp[i*solc+0] = iter;
845 for(j=1;j<solc;j++) {
846 solp[i*solc+j]=0.;
847 }
848 }
849 gsl_multiroot_fdfsolver_free(s);
850 OK
851}
852
853//-------------- non linear least squares fitting -------------------
854
855int nlfit(int method, int f(int, double*, int, double*),
856 int jac(int, double*, int, int, double*),
857 double epsabs, double epsrel, int maxit, int p,
858 KRVEC(xi), RMAT(sol)) {
859 REQUIRES(solr == maxit && solc == 2+xin,BAD_SIZE);
860 DEBUGMSG("nlfit");
861 const gsl_multifit_fdfsolver_type *T;
862 gsl_multifit_fdfsolver *s;
863 gsl_multifit_function_fdf my_f;
864 // extract function from pars
865 my_f.f = f_aux;
866 my_f.df = jf_aux;
867 my_f.fdf = fjf_aux;
868 my_f.n = p;
869 my_f.p = xin; // !!!!
870 Tfjf stfjf;
871 stfjf.f = f;
872 stfjf.jf = jac;
873 my_f.params = &stfjf;
874 size_t iter = 0;
875 int status;
876
877 KDVVIEW(xi);
878 //DMVIEW(cov);
879
880 switch(method) {
881 case 0 : { T = gsl_multifit_fdfsolver_lmsder; break; }
882 case 1 : { T = gsl_multifit_fdfsolver_lmder; break; }
883 default: ERROR(BAD_CODE);
884 }
885
886 s = gsl_multifit_fdfsolver_alloc (T, my_f.n, my_f.p);
887 gsl_multifit_fdfsolver_set (s, &my_f, V(xi));
888
889 do { status = gsl_multifit_fdfsolver_iterate (s);
890
891 solp[iter*solc+0] = iter+1;
892 solp[iter*solc+1] = gsl_blas_dnrm2 (s->f);
893
894 int k;
895 for(k=0;k<xin;k++) {
896 solp[iter*solc+k+2] = gsl_vector_get(s->x,k);
897 }
898
899 iter++;
900 if (status) /* check if solver is stuck */
901 break;
902
903 status = gsl_multifit_test_delta (s->dx, s->x, epsabs, epsrel);
904 }
905 while (status == GSL_CONTINUE && iter < maxit);
906
907 int i,j;
908 for (i=iter; i<solr; i++) {
909 solp[i*solc+0] = iter;
910 for(j=1;j<solc;j++) {
911 solp[i*solc+j]=0.;
912 }
913 }
914
915 //gsl_multifit_covar (s->J, 0.0, M(cov));
916
917 gsl_multifit_fdfsolver_free (s);
918 OK
919}
920
921
922//////////////////////////////////////////////////////
923
924
925#define RAN(C,F) case C: { for(k=0;k<rn;k++) { rp[k]= F(gen); }; OK }
926
927int random_vector_GSL(int seed, int code, RVEC(r)) {
928 DEBUGMSG("random_vector_GSL")
929 static gsl_rng * gen = NULL;
930 if (!gen) { gen = gsl_rng_alloc (gsl_rng_mt19937);}
931 gsl_rng_set (gen, seed);
932 int k;
933 switch (code) {
934 RAN(0,gsl_rng_uniform)
935 RAN(1,gsl_ran_ugaussian)
936 default: ERROR(BAD_CODE);
937 }
938}
939#undef RAN
940
941//////////////////////////////////////////////////////
942
943#include "gsl-ode.c"
944
945//////////////////////////////////////////////////////
diff --git a/packages/hmatrix/src/Numeric/GSL/gsl-ode.c b/packages/hmatrix/src/Numeric/GSL/gsl-ode.c
deleted file mode 100644
index 3f2771b..0000000
--- a/packages/hmatrix/src/Numeric/GSL/gsl-ode.c
+++ /dev/null
@@ -1,182 +0,0 @@
1
2#ifdef GSLODE1
3
4////////////////////////////// ODE V1 //////////////////////////////////////////
5
6#include <gsl/gsl_odeiv.h>
7
8typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
9
10int odefunc (double t, const double y[], double f[], void *params) {
11 Tode * P = (Tode*) params;
12 (P->f)(t,P->n,y,P->n,f);
13 return GSL_SUCCESS;
14}
15
16int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
17 Tode * P = ((Tode*) params);
18 (P->j)(t,P->n,y,P->n,P->n,dfdy);
19 int j;
20 for (j=0; j< P->n; j++)
21 dfdt[j] = 0.0;
22 return GSL_SUCCESS;
23}
24
25
26int ode(int method, double h, double eps_abs, double eps_rel,
27 int f(double, int, const double*, int, double*),
28 int jac(double, int, const double*, int, int, double*),
29 KRVEC(xi), KRVEC(ts), RMAT(sol)) {
30
31 const gsl_odeiv_step_type * T;
32
33 switch(method) {
34 case 0 : {T = gsl_odeiv_step_rk2; break; }
35 case 1 : {T = gsl_odeiv_step_rk4; break; }
36 case 2 : {T = gsl_odeiv_step_rkf45; break; }
37 case 3 : {T = gsl_odeiv_step_rkck; break; }
38 case 4 : {T = gsl_odeiv_step_rk8pd; break; }
39 case 5 : {T = gsl_odeiv_step_rk2imp; break; }
40 case 6 : {T = gsl_odeiv_step_rk4imp; break; }
41 case 7 : {T = gsl_odeiv_step_bsimp; break; }
42 case 8 : { printf("Sorry: ODE rk1imp not available in this GSL version\n"); exit(0); }
43 case 9 : { printf("Sorry: ODE msadams not available in this GSL version\n"); exit(0); }
44 case 10: { printf("Sorry: ODE msbdf not available in this GSL version\n"); exit(0); }
45 default: ERROR(BAD_CODE);
46 }
47
48 gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, xin);
49 gsl_odeiv_control * c = gsl_odeiv_control_y_new (eps_abs, eps_rel);
50 gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (xin);
51
52 Tode P;
53 P.f = f;
54 P.j = jac;
55 P.n = xin;
56
57 gsl_odeiv_system sys = {odefunc, odejac, xin, &P};
58
59 double t = tsp[0];
60
61 double* y = (double*)calloc(xin,sizeof(double));
62 int i,j;
63 for(i=0; i< xin; i++) {
64 y[i] = xip[i];
65 solp[i] = xip[i];
66 }
67
68 for (i = 1; i < tsn ; i++)
69 {
70 double ti = tsp[i];
71 while (t < ti)
72 {
73 gsl_odeiv_evolve_apply (e, c, s,
74 &sys,
75 &t, ti, &h,
76 y);
77 // if (h < hmin) h = hmin;
78 }
79 for(j=0; j<xin; j++) {
80 solp[i*xin + j] = y[j];
81 }
82 }
83
84 free(y);
85 gsl_odeiv_evolve_free (e);
86 gsl_odeiv_control_free (c);
87 gsl_odeiv_step_free (s);
88 return 0;
89}
90
91#else
92
93///////////////////// ODE V2 ///////////////////////////////////////////////////
94
95#include <gsl/gsl_odeiv2.h>
96
97typedef struct {int n; int (*f)(double,int, const double*, int, double *); int (*j)(double,int, const double*, int, int, double*);} Tode;
98
99int odefunc (double t, const double y[], double f[], void *params) {
100 Tode * P = (Tode*) params;
101 (P->f)(t,P->n,y,P->n,f);
102 return GSL_SUCCESS;
103}
104
105int odejac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
106 Tode * P = ((Tode*) params);
107 (P->j)(t,P->n,y,P->n,P->n,dfdy);
108 int j;
109 for (j=0; j< P->n; j++)
110 dfdt[j] = 0.0;
111 return GSL_SUCCESS;
112}
113
114
115int ode(int method, double h, double eps_abs, double eps_rel,
116 int f(double, int, const double*, int, double*),
117 int jac(double, int, const double*, int, int, double*),
118 KRVEC(xi), KRVEC(ts), RMAT(sol)) {
119
120 const gsl_odeiv2_step_type * T;
121
122 switch(method) {
123 case 0 : {T = gsl_odeiv2_step_rk2; break; }
124 case 1 : {T = gsl_odeiv2_step_rk4; break; }
125 case 2 : {T = gsl_odeiv2_step_rkf45; break; }
126 case 3 : {T = gsl_odeiv2_step_rkck; break; }
127 case 4 : {T = gsl_odeiv2_step_rk8pd; break; }
128 case 5 : {T = gsl_odeiv2_step_rk2imp; break; }
129 case 6 : {T = gsl_odeiv2_step_rk4imp; break; }
130 case 7 : {T = gsl_odeiv2_step_bsimp; break; }
131 case 8 : {T = gsl_odeiv2_step_rk1imp; break; }
132 case 9 : {T = gsl_odeiv2_step_msadams; break; }
133 case 10: {T = gsl_odeiv2_step_msbdf; break; }
134 default: ERROR(BAD_CODE);
135 }
136
137 Tode P;
138 P.f = f;
139 P.j = jac;
140 P.n = xin;
141
142 gsl_odeiv2_system sys = {odefunc, odejac, xin, &P};
143
144 gsl_odeiv2_driver * d =
145 gsl_odeiv2_driver_alloc_y_new (&sys, T, h, eps_abs, eps_rel);
146
147 double t = tsp[0];
148
149 double* y = (double*)calloc(xin,sizeof(double));
150 int i,j;
151 int status=0;
152 for(i=0; i< xin; i++) {
153 y[i] = xip[i];
154 solp[i] = xip[i];
155 }
156
157 for (i = 1; i < tsn ; i++)
158 {
159 double ti = tsp[i];
160
161 status = gsl_odeiv2_driver_apply (d, &t, ti, y);
162
163 if (status != GSL_SUCCESS) {
164 printf ("error in ode, return value=%d\n", status);
165 break;
166 }
167
168// printf ("%.5e %.5e %.5e\n", t, y[0], y[1]);
169
170 for(j=0; j<xin; j++) {
171 solp[i*xin + j] = y[j];
172 }
173 }
174
175 free(y);
176 gsl_odeiv2_driver_free (d);
177
178 return status;
179}
180
181#endif
182
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-02 14:50:44 +0000