1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
|
{- |
Module : Numeric.GSL.Root
Copyright : (c) Alberto Ruiz 2009
License : GPL
Maintainer : Alberto Ruiz (aruiz at um dot es)
Stability : provisional
Portability : uses ffi
Multidimensional root finding.
<http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Root_002dFinding.html>
The example in the GSL manual:
@import Numeric.GSL
import Numeric.LinearAlgebra(format)
import Text.Printf(printf)
rosenbrock a b [x,y] = [ a*(1-x), b*(y-x^2) ]
disp = putStrLn . format \" \" (printf \"%.3f\")
main = do
let (sol,path) = root Hybrids 1E-7 30 (rosenbrock 1 10) [-10,-5]
print sol
disp path
\> main
[1.0,1.0]
0.000 -10.000 -5.000 11.000 -1050.000
1.000 -3.976 24.827 4.976 90.203
2.000 -3.976 24.827 4.976 90.203
3.000 -3.976 24.827 4.976 90.203
4.000 -1.274 -5.680 2.274 -73.018
5.000 -1.274 -5.680 2.274 -73.018
6.000 0.249 0.298 0.751 2.359
7.000 0.249 0.298 0.751 2.359
8.000 1.000 0.878 -0.000 -1.218
9.000 1.000 0.989 -0.000 -0.108
10.000 1.000 1.000 0.000 0.000
@
-}
-----------------------------------------------------------------------------
module Numeric.GSL.Root (
root, RootMethod(..),
rootJ, RootMethodJ(..),
) where
import Data.Packed.Internal
import Data.Packed.Matrix
import Foreign
import Foreign.C.Types(CInt)
-------------------------------------------------------------------------
data RootMethod = Hybrids
| Hybrid
| DNewton
| Broyden
deriving (Enum,Eq,Show)
-- | Nonlinear multidimensional root finding using algorithms that do not require
-- any derivative information to be supplied by the user.
-- Any derivatives needed are approximated by finite differences.
root :: RootMethod
-> Double -- ^ maximum residual
-> Int -- ^ maximum number of iterations allowed
-> ([Double] -> [Double]) -- ^ function to minimize
-> [Double] -- ^ starting point
-> ([Double], Matrix Double) -- ^ solution vector and optimization path
root method epsabs maxit fun xinit = rootGen (fi (fromEnum method)) fun xinit epsabs maxit
rootGen m f xi epsabs maxit = unsafePerformIO $ do
let xiv = fromList xi
n = dim xiv
fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
rawpath <- withVector xiv $ \xiv' ->
createMIO maxit (2*n+1)
(c_root m fp epsabs (fi maxit) // xiv')
"root"
let it = round (rawpath @@> (maxit-1,0))
path = takeRows it rawpath
[sol] = toLists $ dropRows (it-1) path
freeHaskellFunPtr fp
return (take n $ drop 1 sol, path)
foreign import ccall "root"
c_root:: CInt -> FunPtr TVV -> Double -> CInt -> TVM
-------------------------------------------------------------------------
data RootMethodJ = HybridsJ
| HybridJ
| Newton
| GNewton
deriving (Enum,Eq,Show)
-- | Nonlinear multidimensional root finding using both the function and its derivatives.
rootJ :: RootMethodJ
-> Double -- ^ maximum residual
-> Int -- ^ maximum number of iterations allowed
-> ([Double] -> [Double]) -- ^ function to minimize
-> ([Double] -> [[Double]]) -- ^ Jacobian
-> [Double] -- ^ starting point
-> ([Double], Matrix Double) -- ^ solution vector and optimization path
rootJ method epsabs maxit fun jac xinit = rootJGen (fi (fromEnum method)) fun jac xinit epsabs maxit
rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
let xiv = fromList xi
n = dim xiv
fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
rawpath <- withVector xiv $ \xiv' ->
createMIO maxit (2*n+1)
(c_rootj m fp jp epsabs (fi maxit) // xiv')
"root"
let it = round (rawpath @@> (maxit-1,0))
path = takeRows it rawpath
[sol] = toLists $ dropRows (it-1) path
freeHaskellFunPtr fp
return (take n $ drop 1 sol, path)
foreign import ccall "rootj"
c_rootj:: CInt -> FunPtr TVV -> FunPtr TVM -> Double -> CInt -> TVM
---------------------------------------------------------------------
foreign import ccall "wrapper"
mkVecVecfun :: TVV -> IO (FunPtr TVV)
aux_vTov :: (Vector Double -> Vector Double) -> TVV
aux_vTov f n p nr r = g where
V {fptr = pr} = f x
x = createV (fromIntegral n) copy "aux_vTov"
copy n' q = do
copyArray q p (fromIntegral n')
return 0
g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral nr)
return 0
foreign import ccall "wrapper"
mkVecMatfun :: TVM -> IO (FunPtr TVM)
aux_vTom :: (Vector Double -> Matrix Double) -> TVM
aux_vTom f n p rr cr r = g where
V {fptr = pr} = flatten $ f x
x = createV (fromIntegral n) copy "aux_vTov"
copy n' q = do
copyArray q p (fromIntegral n')
return 0
g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
return 0
createV n fun msg = unsafePerformIO $ do
r <- createVector n
app1 fun vec r msg
return r
createMIO r c fun msg = do
res <- createMatrix RowMajor r c
app1 fun mat res msg
return res
checkdim1 n v
| dim v == n = v
| otherwise = error $ "Error: "++ show n
++ " results expected in the function supplied to root"
checkdim2 n m
| rows m == n && cols m == n = m
| otherwise = error $ "Error: "++ show n ++ "x" ++ show n
++ " Jacobian expected in root"
|
