summaryrefslogtreecommitdiff
path: root/packages/gsl/src/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'packages/gsl/src/Numeric')
-rw-r--r--packages/gsl/src/Numeric/GSL/Fitting.hs16
-rw-r--r--packages/gsl/src/Numeric/GSL/Fourier.hs5
-rw-r--r--packages/gsl/src/Numeric/GSL/IO.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Internal.hs7
-rw-r--r--packages/gsl/src/Numeric/GSL/Interpolation.hs6
-rw-r--r--packages/gsl/src/Numeric/GSL/LinearAlgebra.hs2
-rw-r--r--packages/gsl/src/Numeric/GSL/Minimization.hs17
-rw-r--r--packages/gsl/src/Numeric/GSL/ODE.hs16
-rw-r--r--packages/gsl/src/Numeric/GSL/Polynomials.hs7
-rw-r--r--packages/gsl/src/Numeric/GSL/Random.hs16
-rw-r--r--packages/gsl/src/Numeric/GSL/Root.hs18
-rw-r--r--packages/gsl/src/Numeric/GSL/Vector.hs3
12 files changed, 63 insertions, 52 deletions
diff --git a/packages/gsl/src/Numeric/GSL/Fitting.hs b/packages/gsl/src/Numeric/GSL/Fitting.hs
index 0a92373..db9d82f 100644
--- a/packages/gsl/src/Numeric/GSL/Fitting.hs
+++ b/packages/gsl/src/Numeric/GSL/Fitting.hs
@@ -1,3 +1,5 @@
1{-# LANGUAGE FlexibleContexts #-}
2
1{- | 3{- |
2Module : Numeric.GSL.Fitting 4Module : Numeric.GSL.Fitting
3Copyright : (c) Alberto Ruiz 2010 5Copyright : (c) Alberto Ruiz 2010
@@ -50,7 +52,7 @@ module Numeric.GSL.Fitting (
50 fitModelScaled, fitModel 52 fitModelScaled, fitModel
51) where 53) where
52 54
53import Numeric.LinearAlgebra 55import Numeric.LinearAlgebra.HMatrix
54import Numeric.GSL.Internal 56import Numeric.GSL.Internal
55 57
56import Foreign.Ptr(FunPtr, freeHaskellFunPtr) 58import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
@@ -80,13 +82,13 @@ nlFitting :: FittingMethod
80nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit 82nlFitting method epsabs epsrel maxit fun jac xinit = nlFitGen (fi (fromEnum method)) fun jac xinit epsabs epsrel maxit
81 83
82nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do 84nlFitGen m f jac xiv epsabs epsrel maxit = unsafePerformIO $ do
83 let p = dim xiv 85 let p = size xiv
84 n = dim (f xiv) 86 n = size (f xiv)
85 fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f)) 87 fp <- mkVecVecfun (aux_vTov (checkdim1 n p . f))
86 jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac)) 88 jp <- mkVecMatfun (aux_vTom (checkdim2 n p . jac))
87 rawpath <- createMatrix RowMajor maxit (2+p) 89 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" 90 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)) 91 let it = round (rawpath `atIndex` (maxit-1,0))
90 path = takeRows it rawpath 92 path = takeRows it rawpath
91 [sol] = toRows $ dropRows (it-1) path 93 [sol] = toRows $ dropRows (it-1) path
92 freeHaskellFunPtr fp 94 freeHaskellFunPtr fp
@@ -99,7 +101,7 @@ foreign import ccall safe "nlfit"
99------------------------------------------------------- 101-------------------------------------------------------
100 102
101checkdim1 n _p v 103checkdim1 n _p v
102 | dim v == n = v 104 | size v == n = v
103 | otherwise = error $ "Error: "++ show n 105 | otherwise = error $ "Error: "++ show n
104 ++ " components expected in the result of the function supplied to nlFitting" 106 ++ " components expected in the result of the function supplied to nlFitting"
105 107
@@ -114,9 +116,9 @@ err (model,deriv) dat vsol = zip sol errs where
114 sol = toList vsol 116 sol = toList vsol
115 c = max 1 (chi/sqrt (fromIntegral dof)) 117 c = max 1 (chi/sqrt (fromIntegral dof))
116 dof = length dat - (rows cov) 118 dof = length dat - (rows cov)
117 chi = norm2 (fromList $ cost (resMs model) dat sol) 119 chi = norm_2 (fromList $ cost (resMs model) dat sol)
118 js = fromLists $ jacobian (resDs deriv) dat sol 120 js = fromLists $ jacobian (resDs deriv) dat sol
119 cov = inv $ trans js <> js 121 cov = inv $ tr js <> js
120 errs = toList $ scalar c * sqrt (takeDiag cov) 122 errs = toList $ scalar c * sqrt (takeDiag cov)
121 123
122 124
diff --git a/packages/gsl/src/Numeric/GSL/Fourier.hs b/packages/gsl/src/Numeric/GSL/Fourier.hs
index 734325b..d824b4f 100644
--- a/packages/gsl/src/Numeric/GSL/Fourier.hs
+++ b/packages/gsl/src/Numeric/GSL/Fourier.hs
@@ -16,14 +16,13 @@ module Numeric.GSL.Fourier (
16 ifft 16 ifft
17) where 17) where
18 18
19import Data.Packed 19import Numeric.LinearAlgebra.HMatrix
20import Numeric.GSL.Internal 20import Numeric.GSL.Internal
21import Data.Complex
22import Foreign.C.Types 21import Foreign.C.Types
23import System.IO.Unsafe (unsafePerformIO) 22import System.IO.Unsafe (unsafePerformIO)
24 23
25genfft code v = unsafePerformIO $ do 24genfft code v = unsafePerformIO $ do
26 r <- createVector (dim v) 25 r <- createVector (size v)
27 app2 (c_fft code) vec v vec r "fft" 26 app2 (c_fft code) vec v vec r "fft"
28 return r 27 return r
29 28
diff --git a/packages/gsl/src/Numeric/GSL/IO.hs b/packages/gsl/src/Numeric/GSL/IO.hs
index 0d6031a..936f6bf 100644
--- a/packages/gsl/src/Numeric/GSL/IO.hs
+++ b/packages/gsl/src/Numeric/GSL/IO.hs
@@ -14,7 +14,7 @@ module Numeric.GSL.IO (
14 fileDimensions, loadMatrix, fromFile 14 fileDimensions, loadMatrix, fromFile
15) where 15) where
16 16
17import Data.Packed 17import Numeric.LinearAlgebra.HMatrix hiding(saveMatrix, loadMatrix)
18import Numeric.GSL.Vector 18import Numeric.GSL.Vector
19import System.Process(readProcess) 19import System.Process(readProcess)
20 20
diff --git a/packages/gsl/src/Numeric/GSL/Internal.hs b/packages/gsl/src/Numeric/GSL/Internal.hs
index a1c4e0c..a269224 100644
--- a/packages/gsl/src/Numeric/GSL/Internal.hs
+++ b/packages/gsl/src/Numeric/GSL/Internal.hs
@@ -22,14 +22,13 @@ module Numeric.GSL.Internal(
22 aux_vTom, 22 aux_vTom,
23 createV, 23 createV,
24 createMIO, 24 createMIO,
25 module Data.Packed.Development, 25 module Numeric.LinearAlgebra.Devel,
26 check, 26 check,
27 Res,TV,TM,TCV,TCM 27 Res,TV,TM,TCV,TCM
28) where 28) where
29 29
30import Data.Packed 30import Numeric.LinearAlgebra.HMatrix
31import Data.Packed.Development hiding (check) 31import Numeric.LinearAlgebra.Devel hiding (check)
32import Data.Complex
33 32
34import Foreign.Marshal.Array(copyArray) 33import Foreign.Marshal.Array(copyArray)
35import Foreign.Ptr(Ptr, FunPtr) 34import Foreign.Ptr(Ptr, FunPtr)
diff --git a/packages/gsl/src/Numeric/GSL/Interpolation.hs b/packages/gsl/src/Numeric/GSL/Interpolation.hs
index 4d72ee2..fb49e40 100644
--- a/packages/gsl/src/Numeric/GSL/Interpolation.hs
+++ b/packages/gsl/src/Numeric/GSL/Interpolation.hs
@@ -32,8 +32,7 @@ module Numeric.GSL.Interpolation (
32 , evaluateIntegralV 32 , evaluateIntegralV
33) where 33) where
34 34
35import Data.Packed.Vector(Vector, fromList, dim) 35import Numeric.LinearAlgebra(Vector, fromList, size, Numeric)
36import Data.Packed.Foreign(appVector)
37import Foreign.C.Types 36import Foreign.C.Types
38import Foreign.Marshal.Alloc(alloca) 37import Foreign.Marshal.Alloc(alloca)
39import Foreign.Ptr(Ptr) 38import Foreign.Ptr(Ptr)
@@ -57,6 +56,9 @@ methodToInt CSplinePeriodic = 3
57methodToInt Akima = 4 56methodToInt Akima = 4
58methodToInt AkimaPeriodic = 5 57methodToInt AkimaPeriodic = 5
59 58
59dim :: Numeric t => Vector t -> Int
60dim = size
61
60applyCFun hsname cname fun mth xs ys x 62applyCFun hsname cname fun mth xs ys x
61 | dim xs /= dim ys = error $ 63 | dim xs /= dim ys = error $
62 "Error: Vectors of unequal sizes " ++ 64 "Error: Vectors of unequal sizes " ++
diff --git a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
index 17e2258..cb78bf4 100644
--- a/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
+++ b/packages/gsl/src/Numeric/GSL/LinearAlgebra.hs
@@ -15,7 +15,7 @@ module Numeric.GSL.LinearAlgebra (
15 fileDimensions, loadMatrix, fromFile 15 fileDimensions, loadMatrix, fromFile
16) where 16) where
17 17
18import Data.Packed 18import Numeric.LinearAlgebra.HMatrix hiding (RandDist,randomVector,saveMatrix,loadMatrix)
19import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) 19import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
20 20
21import Foreign.Marshal.Alloc(free) 21import Foreign.Marshal.Alloc(free)
diff --git a/packages/gsl/src/Numeric/GSL/Minimization.hs b/packages/gsl/src/Numeric/GSL/Minimization.hs
index 056d463..00e0619 100644
--- a/packages/gsl/src/Numeric/GSL/Minimization.hs
+++ b/packages/gsl/src/Numeric/GSL/Minimization.hs
@@ -1,3 +1,6 @@
1{-# LANGUAGE FlexibleContexts #-}
2
3
1{- | 4{- |
2Module : Numeric.GSL.Minimization 5Module : Numeric.GSL.Minimization
3Copyright : (c) Alberto Ruiz 2006-9 6Copyright : (c) Alberto Ruiz 2006-9
@@ -56,7 +59,7 @@ module Numeric.GSL.Minimization (
56) where 59) where
57 60
58 61
59import Data.Packed 62import Numeric.LinearAlgebra.HMatrix hiding(step)
60import Numeric.GSL.Internal 63import Numeric.GSL.Internal
61 64
62import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr) 65import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
@@ -99,7 +102,7 @@ uniMinimizeGen m f xmin xl xu epsrel maxit = unsafePerformIO $ do
99 rawpath <- createMIO maxit 4 102 rawpath <- createMIO maxit 4
100 (c_uniMinize m fp epsrel (fi maxit) xmin xl xu) 103 (c_uniMinize m fp epsrel (fi maxit) xmin xl xu)
101 "uniMinimize" 104 "uniMinimize"
102 let it = round (rawpath @@> (maxit-1,0)) 105 let it = round (rawpath `atIndex` (maxit-1,0))
103 path = takeRows it rawpath 106 path = takeRows it rawpath
104 [sol] = toLists $ dropRows (it-1) path 107 [sol] = toLists $ dropRows (it-1) path
105 freeHaskellFunPtr fp 108 freeHaskellFunPtr fp
@@ -137,13 +140,13 @@ minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList s
137ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 140ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
138 141
139minimizeV method eps maxit szv f xiv = unsafePerformIO $ do 142minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
140 let n = dim xiv 143 let n = size xiv
141 fp <- mkVecfun (iv f) 144 fp <- mkVecfun (iv f)
142 rawpath <- ww2 vec xiv vec szv $ \xiv' szv' -> 145 rawpath <- ww2 vec xiv vec szv $ \xiv' szv' ->
143 createMIO maxit (n+3) 146 createMIO maxit (n+3)
144 (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv') 147 (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
145 "minimize" 148 "minimize"
146 let it = round (rawpath @@> (maxit-1,0)) 149 let it = round (rawpath `atIndex` (maxit-1,0))
147 path = takeRows it rawpath 150 path = takeRows it rawpath
148 sol = flatten $ dropColumns 3 $ dropRows (it-1) path 151 sol = flatten $ dropColumns 3 $ dropRows (it-1) path
149 freeHaskellFunPtr fp 152 freeHaskellFunPtr fp
@@ -191,7 +194,7 @@ minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD
191 194
192 195
193minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do 196minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
194 let n = dim xiv 197 let n = size xiv
195 f' = f 198 f' = f
196 df' = (checkdim1 n . df) 199 df' = (checkdim1 n . df)
197 fp <- mkVecfun (iv f') 200 fp <- mkVecfun (iv f')
@@ -200,7 +203,7 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
200 createMIO maxit (n+2) 203 createMIO maxit (n+2)
201 (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv') 204 (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
202 "minimizeD" 205 "minimizeD"
203 let it = round (rawpath @@> (maxit-1,0)) 206 let it = round (rawpath `atIndex` (maxit-1,0))
204 path = takeRows it rawpath 207 path = takeRows it rawpath
205 sol = flatten $ dropColumns 2 $ dropRows (it-1) path 208 sol = flatten $ dropColumns 2 $ dropRows (it-1) path
206 freeHaskellFunPtr fp 209 freeHaskellFunPtr fp
@@ -217,6 +220,6 @@ foreign import ccall safe "gsl-aux.h minimizeD"
217--------------------------------------------------------------------- 220---------------------------------------------------------------------
218 221
219checkdim1 n v 222checkdim1 n v
220 | dim v == n = v 223 | size v == n = v
221 | otherwise = error $ "Error: "++ show n 224 | otherwise = error $ "Error: "++ show n
222 ++ " components expected in the result of the gradient supplied to minimizeD" 225 ++ " components expected in the result of the gradient supplied to minimizeD"
diff --git a/packages/gsl/src/Numeric/GSL/ODE.hs b/packages/gsl/src/Numeric/GSL/ODE.hs
index 7549a65..3258b83 100644
--- a/packages/gsl/src/Numeric/GSL/ODE.hs
+++ b/packages/gsl/src/Numeric/GSL/ODE.hs
@@ -1,3 +1,6 @@
1{-# LANGUAGE FlexibleContexts #-}
2
3
1{- | 4{- |
2Module : Numeric.GSL.ODE 5Module : Numeric.GSL.ODE
3Copyright : (c) Alberto Ruiz 2010 6Copyright : (c) Alberto Ruiz 2010
@@ -32,7 +35,7 @@ module Numeric.GSL.ODE (
32 odeSolve, odeSolveV, ODEMethod(..), Jacobian 35 odeSolve, odeSolveV, ODEMethod(..), Jacobian
33) where 36) where
34 37
35import Data.Packed 38import Numeric.LinearAlgebra.HMatrix
36import Numeric.GSL.Internal 39import Numeric.GSL.Internal
37 40
38import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr) 41import Foreign.Ptr(FunPtr, nullFunPtr, freeHaskellFunPtr)
@@ -68,7 +71,7 @@ odeSolve
68 -> Vector Double -- ^ desired solution times 71 -> Vector Double -- ^ desired solution times
69 -> Matrix Double -- ^ solution 72 -> Matrix Double -- ^ solution
70odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts 73odeSolve xdot xi ts = odeSolveV RKf45 hi epsAbs epsRel (l2v xdot) (fromList xi) ts
71 where hi = (ts@>1 - ts@>0)/100 74 where hi = (ts!1 - ts!0)/100
72 epsAbs = 1.49012e-08 75 epsAbs = 1.49012e-08
73 epsRel = 1.49012e-08 76 epsRel = 1.49012e-08
74 l2v f = \t -> fromList . f t . toList 77 l2v f = \t -> fromList . f t . toList
@@ -107,14 +110,14 @@ odeSolveV'
107 -> Vector Double -- ^ desired solution times 110 -> Vector Double -- ^ desired solution times
108 -> Matrix Double -- ^ solution 111 -> Matrix Double -- ^ solution
109odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do 112odeSolveV' method mbjac h epsAbs epsRel f xiv ts = unsafePerformIO $ do
110 let n = dim xiv 113 let n = size xiv
111 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t)) 114 fp <- mkDoubleVecVecfun (\t -> aux_vTov (checkdim1 n . f t))
112 jp <- case mbjac of 115 jp <- case mbjac of
113 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) 116 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
114 Nothing -> return nullFunPtr 117 Nothing -> return nullFunPtr
115 sol <- vec xiv $ \xiv' -> 118 sol <- vec xiv $ \xiv' ->
116 vec (checkTimes ts) $ \ts' -> 119 vec (checkTimes ts) $ \ts' ->
117 createMIO (dim ts) n 120 createMIO (size ts) n
118 (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' ) 121 (ode_c (method) h epsAbs epsRel fp jp // xiv' // ts' )
119 "ode" 122 "ode"
120 freeHaskellFunPtr fp 123 freeHaskellFunPtr fp
@@ -126,7 +129,7 @@ foreign import ccall safe "ode"
126------------------------------------------------------- 129-------------------------------------------------------
127 130
128checkdim1 n v 131checkdim1 n v
129 | dim v == n = v 132 | size v == n = v
130 | otherwise = error $ "Error: "++ show n 133 | otherwise = error $ "Error: "++ show n
131 ++ " components expected in the result of the function supplied to odeSolve" 134 ++ " components expected in the result of the function supplied to odeSolve"
132 135
@@ -135,6 +138,7 @@ checkdim2 n m
135 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n 138 | otherwise = error $ "Error: "++ show n ++ "x" ++ show n
136 ++ " Jacobian expected in odeSolve" 139 ++ " Jacobian expected in odeSolve"
137 140
138checkTimes ts | dim ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts 141checkTimes ts | size ts > 1 && all (>0) (zipWith subtract ts' (tail ts')) = ts
139 | otherwise = error "odeSolve requires increasing times" 142 | otherwise = error "odeSolve requires increasing times"
140 where ts' = toList ts 143 where ts' = toList ts
144
diff --git a/packages/gsl/src/Numeric/GSL/Polynomials.hs b/packages/gsl/src/Numeric/GSL/Polynomials.hs
index b1be85d..246e301 100644
--- a/packages/gsl/src/Numeric/GSL/Polynomials.hs
+++ b/packages/gsl/src/Numeric/GSL/Polynomials.hs
@@ -16,9 +16,8 @@ module Numeric.GSL.Polynomials (
16 polySolve 16 polySolve
17) where 17) where
18 18
19import Data.Packed 19import Numeric.LinearAlgebra.HMatrix
20import Numeric.GSL.Internal 20import Numeric.GSL.Internal
21import Data.Complex
22import System.IO.Unsafe (unsafePerformIO) 21import System.IO.Unsafe (unsafePerformIO)
23 22
24#if __GLASGOW_HASKELL__ >= 704 23#if __GLASGOW_HASKELL__ >= 704
@@ -47,8 +46,8 @@ polySolve :: [Double] -> [Complex Double]
47polySolve = toList . polySolve' . fromList 46polySolve = toList . polySolve' . fromList
48 47
49polySolve' :: Vector Double -> Vector (Complex Double) 48polySolve' :: Vector Double -> Vector (Complex Double)
50polySolve' v | dim v > 1 = unsafePerformIO $ do 49polySolve' v | size v > 1 = unsafePerformIO $ do
51 r <- createVector (dim v-1) 50 r <- createVector (size v-1)
52 app2 c_polySolve vec v vec r "polySolve" 51 app2 c_polySolve vec v vec r "polySolve"
53 return r 52 return r
54 | otherwise = error "polySolve on a polynomial of degree zero" 53 | otherwise = error "polySolve on a polynomial of degree zero"
diff --git a/packages/gsl/src/Numeric/GSL/Random.hs b/packages/gsl/src/Numeric/GSL/Random.hs
index f1f49e5..139c921 100644
--- a/packages/gsl/src/Numeric/GSL/Random.hs
+++ b/packages/gsl/src/Numeric/GSL/Random.hs
@@ -21,11 +21,13 @@ module Numeric.GSL.Random (
21) where 21) where
22 22
23import Numeric.GSL.Vector 23import Numeric.GSL.Vector
24import Numeric.LinearAlgebra(cholSH) 24import Numeric.LinearAlgebra.HMatrix hiding (
25import Numeric.Container hiding (
26 randomVector, 25 randomVector,
27 gaussianSample, 26 gaussianSample,
28 uniformSample 27 uniformSample,
28 Seed,
29 rand,
30 randn
29 ) 31 )
30import System.Random(randomIO) 32import System.Random(randomIO)
31 33
@@ -40,10 +42,10 @@ gaussianSample :: Seed
40 -> Matrix Double -- ^ covariance matrix 42 -> Matrix Double -- ^ covariance matrix
41 -> Matrix Double -- ^ result 43 -> Matrix Double -- ^ result
42gaussianSample seed n med cov = m where 44gaussianSample seed n med cov = m where
43 c = dim med 45 c = size med
44 meds = konst 1 n `outer` med 46 meds = konst 1 n `outer` med
45 rs = reshape c $ randomVector seed Gaussian (c * n) 47 rs = reshape c $ randomVector seed Gaussian (c * n)
46 m = rs `mXm` cholSH cov `add` meds 48 m = rs <> cholSH cov + meds
47 49
48-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate 50-- | Obtains a matrix whose rows are pseudorandom samples from a multivariate
49-- uniform distribution. 51-- uniform distribution.
@@ -55,10 +57,10 @@ uniformSample seed n rgs = m where
55 (as,bs) = unzip rgs 57 (as,bs) = unzip rgs
56 a = fromList as 58 a = fromList as
57 cs = zipWith subtract as bs 59 cs = zipWith subtract as bs
58 d = dim a 60 d = size a
59 dat = toRows $ reshape n $ randomVector seed Uniform (n*d) 61 dat = toRows $ reshape n $ randomVector seed Uniform (n*d)
60 am = konst 1 n `outer` a 62 am = konst 1 n `outer` a
61 m = fromColumns (zipWith scale cs dat) `add` am 63 m = fromColumns (zipWith scale cs dat) + am
62 64
63-- | pseudorandom matrix with uniform elements between 0 and 1 65-- | pseudorandom matrix with uniform elements between 0 and 1
64randm :: RandDist 66randm :: RandDist
diff --git a/packages/gsl/src/Numeric/GSL/Root.hs b/packages/gsl/src/Numeric/GSL/Root.hs
index b9f3b94..724f32f 100644
--- a/packages/gsl/src/Numeric/GSL/Root.hs
+++ b/packages/gsl/src/Numeric/GSL/Root.hs
@@ -1,3 +1,5 @@
1{-# LANGUAGE FlexibleContexts #-}
2
1{- | 3{- |
2Module : Numeric.GSL.Root 4Module : Numeric.GSL.Root
3Copyright : (c) Alberto Ruiz 2009 5Copyright : (c) Alberto Ruiz 2009
@@ -39,7 +41,7 @@ module Numeric.GSL.Root (
39 rootJ, RootMethodJ(..), 41 rootJ, RootMethodJ(..),
40) where 42) where
41 43
42import Data.Packed 44import Numeric.LinearAlgebra.HMatrix
43import Numeric.GSL.Internal 45import Numeric.GSL.Internal
44import Foreign.Ptr(FunPtr, freeHaskellFunPtr) 46import Foreign.Ptr(FunPtr, freeHaskellFunPtr)
45import Foreign.C.Types 47import Foreign.C.Types
@@ -69,7 +71,7 @@ uniRootGen m f xl xu epsrel maxit = unsafePerformIO $ do
69 rawpath <- createMIO maxit 4 71 rawpath <- createMIO maxit 4
70 (c_root m fp epsrel (fi maxit) xl xu) 72 (c_root m fp epsrel (fi maxit) xl xu)
71 "root" 73 "root"
72 let it = round (rawpath @@> (maxit-1,0)) 74 let it = round (rawpath `atIndex` (maxit-1,0))
73 path = takeRows it rawpath 75 path = takeRows it rawpath
74 [sol] = toLists $ dropRows (it-1) path 76 [sol] = toLists $ dropRows (it-1) path
75 freeHaskellFunPtr fp 77 freeHaskellFunPtr fp
@@ -100,7 +102,7 @@ uniRootJGen m f df x epsrel maxit = unsafePerformIO $ do
100 rawpath <- createMIO maxit 2 102 rawpath <- createMIO maxit 2
101 (c_rootj m fp dfp epsrel (fi maxit) x) 103 (c_rootj m fp dfp epsrel (fi maxit) x)
102 "rootj" 104 "rootj"
103 let it = round (rawpath @@> (maxit-1,0)) 105 let it = round (rawpath `atIndex` (maxit-1,0))
104 path = takeRows it rawpath 106 path = takeRows it rawpath
105 [sol] = toLists $ dropRows (it-1) path 107 [sol] = toLists $ dropRows (it-1) path
106 freeHaskellFunPtr fp 108 freeHaskellFunPtr fp
@@ -132,13 +134,13 @@ root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit ep
132 134
133rootGen m f xi epsabs maxit = unsafePerformIO $ do 135rootGen m f xi epsabs maxit = unsafePerformIO $ do
134 let xiv = fromList xi 136 let xiv = fromList xi
135 n = dim xiv 137 n = size xiv
136 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) 138 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
137 rawpath <- vec xiv $ \xiv' -> 139 rawpath <- vec xiv $ \xiv' ->
138 createMIO maxit (2*n+1) 140 createMIO maxit (2*n+1)
139 (c_multiroot m fp epsabs (fi maxit) // xiv') 141 (c_multiroot m fp epsabs (fi maxit) // xiv')
140 "multiroot" 142 "multiroot"
141 let it = round (rawpath @@> (maxit-1,0)) 143 let it = round (rawpath `atIndex` (maxit-1,0))
142 path = takeRows it rawpath 144 path = takeRows it rawpath
143 [sol] = toLists $ dropRows (it-1) path 145 [sol] = toLists $ dropRows (it-1) path
144 freeHaskellFunPtr fp 146 freeHaskellFunPtr fp
@@ -169,14 +171,14 @@ rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun ja
169 171
170rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do 172rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
171 let xiv = fromList xi 173 let xiv = fromList xi
172 n = dim xiv 174 n = size xiv
173 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) 175 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
174 jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList)) 176 jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
175 rawpath <- vec xiv $ \xiv' -> 177 rawpath <- vec xiv $ \xiv' ->
176 createMIO maxit (2*n+1) 178 createMIO maxit (2*n+1)
177 (c_multirootj m fp jp epsabs (fi maxit) // xiv') 179 (c_multirootj m fp jp epsabs (fi maxit) // xiv')
178 "multiroot" 180 "multiroot"
179 let it = round (rawpath @@> (maxit-1,0)) 181 let it = round (rawpath `atIndex` (maxit-1,0))
180 path = takeRows it rawpath 182 path = takeRows it rawpath
181 [sol] = toLists $ dropRows (it-1) path 183 [sol] = toLists $ dropRows (it-1) path
182 freeHaskellFunPtr fp 184 freeHaskellFunPtr fp
@@ -189,7 +191,7 @@ foreign import ccall safe "multirootj"
189------------------------------------------------------- 191-------------------------------------------------------
190 192
191checkdim1 n v 193checkdim1 n v
192 | dim v == n = v 194 | size v == n = v
193 | otherwise = error $ "Error: "++ show n 195 | otherwise = error $ "Error: "++ show n
194 ++ " components expected in the result of the function supplied to root" 196 ++ " components expected in the result of the function supplied to root"
195 197
diff --git a/packages/gsl/src/Numeric/GSL/Vector.hs b/packages/gsl/src/Numeric/GSL/Vector.hs
index af79f32..0cd99eb 100644
--- a/packages/gsl/src/Numeric/GSL/Vector.hs
+++ b/packages/gsl/src/Numeric/GSL/Vector.hs
@@ -14,8 +14,7 @@ module Numeric.GSL.Vector (
14 fwriteVector, freadVector, fprintfVector, fscanfVector 14 fwriteVector, freadVector, fprintfVector, fscanfVector
15) where 15) where
16 16
17import Data.Packed 17import Numeric.LinearAlgebra.HMatrix hiding(randomVector, saveMatrix)
18import Numeric.LinearAlgebra(RandDist(..))
19import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM) 18import Numeric.GSL.Internal hiding (TV,TM,TCV,TCM)
20 19
21import Foreign.Marshal.Alloc(free) 20import Foreign.Marshal.Alloc(free)
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-02 15:18:48 +0000