summaryrefslogtreecommitdiff
path: root/lib/Numeric
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric')
-rw-r--r--lib/Numeric/GSL/Minimization.hs93
1 files changed, 56 insertions, 37 deletions
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
index a595da9..3240ddf 100644
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -51,8 +51,9 @@ The nmsimplex2 version is a new O(N) implementation of the earlier O(N^2) nmsimp
51 51
52----------------------------------------------------------------------------- 52-----------------------------------------------------------------------------
53module Numeric.GSL.Minimization ( 53module Numeric.GSL.Minimization (
54 minimize, MinimizeMethod(..), 54 minimize, minimizeV, MinimizeMethod(..),
55 minimizeD, MinimizeMethodD(..), 55 minimizeD, minimizeVD, MinimizeMethodD(..),
56
56 minimizeNMSimplex, 57 minimizeNMSimplex,
57 minimizeConjugateGradient, 58 minimizeConjugateGradient,
58 minimizeVectorBFGS2 59 minimizeVectorBFGS2
@@ -82,7 +83,7 @@ data MinimizeMethod = NMSimplex
82 | NMSimplex2 83 | NMSimplex2
83 deriving (Enum,Eq,Show,Bounded) 84 deriving (Enum,Eq,Show,Bounded)
84 85
85-- | Minimization without derivatives. 86-- | Minimization without derivatives
86minimize :: MinimizeMethod 87minimize :: MinimizeMethod
87 -> Double -- ^ desired precision of the solution (size test) 88 -> Double -- ^ desired precision of the solution (size test)
88 -> Int -- ^ maximum number of iterations allowed 89 -> Int -- ^ maximum number of iterations allowed
@@ -91,7 +92,39 @@ minimize :: MinimizeMethod
91 -> [Double] -- ^ starting point 92 -> [Double] -- ^ starting point
92 -> ([Double], Matrix Double) -- ^ solution vector and optimization path 93 -> ([Double], Matrix Double) -- ^ solution vector and optimization path
93 94
94minimize method = minimizeGen (fi (fromEnum method)) 95-- | Minimization without derivatives (vector version)
96minimizeV :: MinimizeMethod
97 -> Double -- ^ desired precision of the solution (size test)
98 -> Int -- ^ maximum number of iterations allowed
99 -> Vector Double -- ^ sizes of the initial search box
100 -> (Vector Double -> Double) -- ^ function to minimize
101 -> Vector Double -- ^ starting point
102 -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
103
104minimize method eps maxit sz f xi = v2l $ minimizeV method eps maxit (fromList sz) (f.toList) (fromList xi)
105 where v2l (v,m) = (toList v, m)
106
107ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
108
109minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
110 let n = dim xiv
111 fp <- mkVecfun (iv f)
112 rawpath <- ww2 withVector xiv withVector szv $ \xiv' szv' ->
113 createMIO maxit (n+3)
114 (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
115 "minimize"
116 let it = round (rawpath @@> (maxit-1,0))
117 path = takeRows it rawpath
118 sol = cdat $ dropColumns 3 $ dropRows (it-1) path
119 freeHaskellFunPtr fp
120 return (sol, path)
121
122
123foreign import ccall "gsl-aux.h minimize"
124 c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM
125
126----------------------------------------------------------------------------------
127
95 128
96data MinimizeMethodD = ConjugateFR 129data MinimizeMethodD = ConjugateFR
97 | ConjugatePR 130 | ConjugatePR
@@ -111,49 +144,35 @@ minimizeD :: MinimizeMethodD
111 -> [Double] -- ^ starting point 144 -> [Double] -- ^ starting point
112 -> ([Double], Matrix Double) -- ^ solution vector and optimization path 145 -> ([Double], Matrix Double) -- ^ solution vector and optimization path
113 146
114minimizeD method = minimizeDGen (fi (fromEnum method)) 147-- | Minimization with derivatives (vector version)
115 148minimizeVD :: MinimizeMethodD
116------------------------------------------------------------------------- 149 -> Double -- ^ desired precision of the solution (gradient test)
117 150 -> Int -- ^ maximum number of iterations allowed
118ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2 151 -> Double -- ^ size of the first trial step
119 152 -> Double -- ^ tol (precise meaning depends on method)
120minimizeGen method eps maxit sz f xi = unsafePerformIO $ do 153 -> (Vector Double -> Double) -- ^ function to minimize
121 let xiv = fromList xi 154 -> (Vector Double -> Vector Double) -- ^ gradient
122 szv = fromList sz 155 -> Vector Double -- ^ starting point
123 n = dim xiv 156 -> (Vector Double, Matrix Double) -- ^ solution vector and optimization path
124 fp <- mkVecfun (iv (f.toList))
125 rawpath <- ww2 withVector xiv withVector szv $ \xiv' szv' ->
126 createMIO maxit (n+3)
127 (c_minimize method fp eps (fi maxit) // xiv' // szv')
128 "minimize"
129 let it = round (rawpath @@> (maxit-1,0))
130 path = takeRows it rawpath
131 [sol] = toLists $ dropRows (it-1) path
132 freeHaskellFunPtr fp
133 return (drop 3 sol, path)
134
135
136foreign import ccall "gsl-aux.h minimize"
137 c_minimize:: CInt -> FunPtr (CInt -> Ptr Double -> Double) -> Double -> CInt -> TVVM
138
139----------------------------------------------------------------------------------
140 157
158minimizeD method eps maxit istep tol f df xi = v2l $ minimizeVD
159 method eps maxit istep tol (f.toList) (fromList.df.toList) (fromList xi)
160 where v2l (v,m) = (toList v, m)
141 161
142 162
143minimizeDGen method eps maxit istep tol f df xi = unsafePerformIO $ do 163minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
144 let xiv = fromList xi 164 let n = dim xiv
145 n = dim xiv 165 f' = f
146 f' = f . toList 166 df' = (checkdim1 n . df)
147 df' = (checkdim1 n .fromList . df . toList)
148 fp <- mkVecfun (iv f') 167 fp <- mkVecfun (iv f')
149 dfp <- mkVecVecfun (aux_vTov df') 168 dfp <- mkVecVecfun (aux_vTov df')
150 rawpath <- withVector xiv $ \xiv' -> 169 rawpath <- withVector xiv $ \xiv' ->
151 createMIO maxit (n+2) 170 createMIO maxit (n+2)
152 (c_minimizeD method fp dfp istep tol eps (fi maxit) // xiv') 171 (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
153 "minimizeD" 172 "minimizeD"
154 let it = round (rawpath @@> (maxit-1,0)) 173 let it = round (rawpath @@> (maxit-1,0))
155 path = takeRows it rawpath 174 path = takeRows it rawpath
156 sol = toList $ cdat $ dropColumns 2 $ dropRows (it-1) path 175 sol = cdat $ dropColumns 2 $ dropRows (it-1) path
157 freeHaskellFunPtr fp 176 freeHaskellFunPtr fp
158 freeHaskellFunPtr dfp 177 freeHaskellFunPtr dfp
159 return (sol,path) 178 return (sol,path)