diff options
author | Alberto Ruiz <aruiz@um.es> | 2015-06-18 13:43:46 +0200 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2015-06-18 13:43:46 +0200 |
commit | 54f7eb384bbfc1f478359dacd90f25759c44d0cf (patch) | |
tree | 1b5f0bf63e37a4f3f0d089d40a8e804e6e55bacf /packages/gsl/src/Numeric/GSL | |
parent | 52009006791ee2b71530a61f4bf9e1c065c04eae (diff) | |
parent | 992c8915ea4fb324a8d43203056d1f0224552236 (diff) |
Merge pull request #132 from peddie/siman
Add an interface to the GSL simulated annealing module.
Diffstat (limited to 'packages/gsl/src/Numeric/GSL')
-rw-r--r-- | packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs | 246 | ||||
-rw-r--r-- | packages/gsl/src/Numeric/GSL/gsl-aux.c | 27 |
2 files changed, 272 insertions, 1 deletions
diff --git a/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs b/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs new file mode 100644 index 0000000..9f9ed97 --- /dev/null +++ b/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs | |||
@@ -0,0 +1,246 @@ | |||
1 | {- | | ||
2 | Module : Numeric.GSL.Interpolation | ||
3 | Copyright : (c) Matthew Peddie 2015 | ||
4 | License : GPL | ||
5 | Maintainer : Alberto Ruiz | ||
6 | Stability : provisional | ||
7 | |||
8 | Simulated annealing routines. | ||
9 | |||
10 | <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing.html#Simulated-Annealing> | ||
11 | |||
12 | Here is a translation of the simple example given in | ||
13 | <https://www.gnu.org/software/gsl/manual/html_node/Trivial-example.html#Trivial-example the GSL manual>: | ||
14 | |||
15 | > import Numeric.GSL.SimulatedAnnealing | ||
16 | > import Numeric.LinearAlgebra.HMatrix | ||
17 | > | ||
18 | > main = print $ simanSolve 0 1 exampleParams 15.5 exampleE exampleM exampleS (Just show) | ||
19 | > | ||
20 | > exampleParams = SimulatedAnnealingParams 200 1000 1.0 1.0 0.008 1.003 2.0e-6 | ||
21 | > | ||
22 | > exampleE x = exp (-(x - 1)**2) * sin (8 * x) | ||
23 | > | ||
24 | > exampleM x y = abs $ x - y | ||
25 | > | ||
26 | > exampleS rands stepSize current = (rands ! 0) * 2 * stepSize - stepSize + current | ||
27 | |||
28 | The manual states: | ||
29 | |||
30 | > The first example, in one dimensional Cartesian space, sets up an | ||
31 | > energy function which is a damped sine wave; this has many local | ||
32 | > minima, but only one global minimum, somewhere between 1.0 and | ||
33 | > 1.5. The initial guess given is 15.5, which is several local minima | ||
34 | > away from the global minimum. | ||
35 | |||
36 | This global minimum is around 1.36. | ||
37 | |||
38 | -} | ||
39 | {-# OPTIONS_GHC -Wall #-} | ||
40 | |||
41 | module Numeric.GSL.SimulatedAnnealing ( | ||
42 | -- * Searching for minima | ||
43 | simanSolve | ||
44 | -- * Configuring the annealing process | ||
45 | , SimulatedAnnealingParams(..) | ||
46 | ) where | ||
47 | |||
48 | import Numeric.GSL.Internal | ||
49 | import Numeric.LinearAlgebra.HMatrix hiding(step) | ||
50 | |||
51 | import Data.Vector.Storable(generateM) | ||
52 | import Foreign.Storable(Storable(..)) | ||
53 | import Foreign.Marshal.Utils(with) | ||
54 | import Foreign.Ptr(Ptr, FunPtr, nullFunPtr) | ||
55 | import Foreign.StablePtr(StablePtr, newStablePtr, deRefStablePtr, freeStablePtr) | ||
56 | import Foreign.C.Types | ||
57 | import System.IO.Unsafe(unsafePerformIO) | ||
58 | import Control.Applicative ((<*>), (<$>)) | ||
59 | |||
60 | import System.IO (hFlush, stdout) | ||
61 | |||
62 | import Data.IORef (IORef, newIORef, writeIORef, readIORef, modifyIORef') | ||
63 | |||
64 | -- | 'SimulatedAnnealingParams' is a translation of the | ||
65 | -- @gsl_siman_params_t@ structure documented in | ||
66 | -- <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing-functions.html#Simulated-Annealing-functions the GSL manual>, | ||
67 | -- which controls the simulated annealing algorithm. | ||
68 | -- | ||
69 | -- The annealing process is parameterized by the Boltzmann | ||
70 | -- distribution and the /cooling schedule/. For more details, see | ||
71 | -- <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing-algorithm.html#Simulated-Annealing-algorithm the relevant section of the manual>. | ||
72 | data SimulatedAnnealingParams = SimulatedAnnealingParams { | ||
73 | n_tries :: CInt -- ^ The number of points to try for each step. | ||
74 | , iters_fixed_T :: CInt -- ^ The number of iterations at each temperature | ||
75 | , step_size :: Double -- ^ The maximum step size in the random walk | ||
76 | , boltzmann_k :: Double -- ^ Boltzmann distribution parameter | ||
77 | , cooling_t_initial :: Double -- ^ Initial temperature | ||
78 | , cooling_mu_t :: Double -- ^ Cooling rate parameter | ||
79 | , cooling_t_min :: Double -- ^ Final temperature | ||
80 | } deriving (Eq, Show, Read) | ||
81 | |||
82 | instance Storable SimulatedAnnealingParams where | ||
83 | sizeOf p = sizeOf (n_tries p) + | ||
84 | sizeOf (iters_fixed_T p) + | ||
85 | sizeOf (step_size p) + | ||
86 | sizeOf (boltzmann_k p) + | ||
87 | sizeOf (cooling_t_initial p) + | ||
88 | sizeOf (cooling_mu_t p) + | ||
89 | sizeOf (cooling_t_min p) | ||
90 | -- TODO(MP): is this safe? | ||
91 | alignment p = alignment (step_size p) | ||
92 | -- TODO(MP): Is there a more automatic way to write these? | ||
93 | peek ptr = SimulatedAnnealingParams <$> | ||
94 | peekByteOff ptr 0 <*> | ||
95 | peekByteOff ptr i <*> | ||
96 | peekByteOff ptr (2*i) <*> | ||
97 | peekByteOff ptr (2*i + d) <*> | ||
98 | peekByteOff ptr (2*i + 2*d) <*> | ||
99 | peekByteOff ptr (2*i + 3*d) <*> | ||
100 | peekByteOff ptr (2*i + 4*d) | ||
101 | where | ||
102 | i = sizeOf (0 :: CInt) | ||
103 | d = sizeOf (0 :: Double) | ||
104 | poke ptr sap = do | ||
105 | pokeByteOff ptr 0 (n_tries sap) | ||
106 | pokeByteOff ptr i (iters_fixed_T sap) | ||
107 | pokeByteOff ptr (2*i) (step_size sap) | ||
108 | pokeByteOff ptr (2*i + d) (boltzmann_k sap) | ||
109 | pokeByteOff ptr (2*i + 2*d) (cooling_t_initial sap) | ||
110 | pokeByteOff ptr (2*i + 3*d) (cooling_mu_t sap) | ||
111 | pokeByteOff ptr (2*i + 4*d) (cooling_t_min sap) | ||
112 | where | ||
113 | i = sizeOf (0 :: CInt) | ||
114 | d = sizeOf (0 :: Double) | ||
115 | |||
116 | -- We use a StablePtr to an IORef so that we can keep hold of | ||
117 | -- StablePtr values but mutate their contents. A simple 'StablePtr a' | ||
118 | -- won't work, since we'd have no way to write 'copyConfig'. | ||
119 | type P a = StablePtr (IORef a) | ||
120 | |||
121 | copyConfig :: P a -> P a -> IO () | ||
122 | copyConfig src' dest' = do | ||
123 | dest <- deRefStablePtr dest' | ||
124 | src <- deRefStablePtr src' | ||
125 | readIORef src >>= writeIORef dest | ||
126 | |||
127 | copyConstructConfig :: P a -> IO (P a) | ||
128 | copyConstructConfig x = do | ||
129 | conf <- deRefRead x | ||
130 | newconf <- newIORef conf | ||
131 | newStablePtr newconf | ||
132 | |||
133 | destroyConfig :: P a -> IO () | ||
134 | destroyConfig p = do | ||
135 | freeStablePtr p | ||
136 | |||
137 | deRefRead :: P a -> IO a | ||
138 | deRefRead p = deRefStablePtr p >>= readIORef | ||
139 | |||
140 | wrapEnergy :: (a -> Double) -> P a -> Double | ||
141 | wrapEnergy f p = unsafePerformIO $ f <$> deRefRead p | ||
142 | |||
143 | wrapMetric :: (a -> a -> Double) -> P a -> P a -> Double | ||
144 | wrapMetric f x y = unsafePerformIO $ f <$> deRefRead x <*> deRefRead y | ||
145 | |||
146 | wrapStep :: Int | ||
147 | -> (Vector Double -> Double -> a -> a) | ||
148 | -> GSLRNG | ||
149 | -> P a | ||
150 | -> Double | ||
151 | -> IO () | ||
152 | wrapStep nrand f (GSLRNG rng) confptr stepSize = do | ||
153 | v <- generateM nrand (\_ -> gslRngUniform rng) | ||
154 | conf <- deRefStablePtr confptr | ||
155 | modifyIORef' conf $ f v stepSize | ||
156 | |||
157 | wrapPrint :: (a -> String) -> P a -> IO () | ||
158 | wrapPrint pf ptr = deRefRead ptr >>= putStr . pf >> hFlush stdout | ||
159 | |||
160 | foreign import ccall safe "wrapper" | ||
161 | mkEnergyFun :: (P a -> Double) -> IO (FunPtr (P a -> Double)) | ||
162 | |||
163 | foreign import ccall safe "wrapper" | ||
164 | mkMetricFun :: (P a -> P a -> Double) -> IO (FunPtr (P a -> P a -> Double)) | ||
165 | |||
166 | foreign import ccall safe "wrapper" | ||
167 | mkStepFun :: (GSLRNG -> P a -> Double -> IO ()) | ||
168 | -> IO (FunPtr (GSLRNG -> P a -> Double -> IO ())) | ||
169 | |||
170 | foreign import ccall safe "wrapper" | ||
171 | mkCopyFun :: (P a -> P a -> IO ()) -> IO (FunPtr (P a -> P a -> IO ())) | ||
172 | |||
173 | foreign import ccall safe "wrapper" | ||
174 | mkCopyConstructorFun :: (P a -> IO (P a)) -> IO (FunPtr (P a -> IO (P a))) | ||
175 | |||
176 | foreign import ccall safe "wrapper" | ||
177 | mkDestructFun :: (P a -> IO ()) -> IO (FunPtr (P a -> IO ())) | ||
178 | |||
179 | newtype GSLRNG = GSLRNG (Ptr GSLRNG) | ||
180 | |||
181 | foreign import ccall safe "gsl_rng.h gsl_rng_uniform" | ||
182 | gslRngUniform :: Ptr GSLRNG -> IO Double | ||
183 | |||
184 | foreign import ccall safe "gsl-aux.h siman" | ||
185 | siman :: CInt -- ^ RNG seed (for repeatability) | ||
186 | -> Ptr SimulatedAnnealingParams -- ^ params | ||
187 | -> P a -- ^ Configuration | ||
188 | -> FunPtr (P a -> Double) -- ^ Energy functional | ||
189 | -> FunPtr (P a -> P a -> Double) -- ^ Metric definition | ||
190 | -> FunPtr (GSLRNG -> P a -> Double -> IO ()) -- ^ Step evaluation | ||
191 | -> FunPtr (P a -> P a -> IO ()) -- ^ Copy config | ||
192 | -> FunPtr (P a -> IO (P a)) -- ^ Copy constructor for config | ||
193 | -> FunPtr (P a -> IO ()) -- ^ Destructor for config | ||
194 | -> FunPtr (P a -> IO ()) -- ^ Print function | ||
195 | -> IO CInt | ||
196 | |||
197 | -- | | ||
198 | -- Calling | ||
199 | -- | ||
200 | -- > simanSolve seed nrand params x0 e m step print | ||
201 | -- | ||
202 | -- performs a simulated annealing search through a given space. So | ||
203 | -- that any configuration type may be used, the space is specified by | ||
204 | -- providing the functions @e@ (the energy functional) and @m@ (the | ||
205 | -- metric definition). @x0@ is the initial configuration of the | ||
206 | -- system. The simulated annealing steps are generated using the | ||
207 | -- user-provided function @step@, which should randomly construct a | ||
208 | -- new system configuration. | ||
209 | -- | ||
210 | -- If 'Nothing' is passed instead of a printing function, no | ||
211 | -- incremental output will be generated. Otherwise, the GSL-formatted | ||
212 | -- output, including the configuration description the user function | ||
213 | -- generates, will be printed to stdout. | ||
214 | -- | ||
215 | -- Each time the step function is called, it is supplied with a random | ||
216 | -- vector containing @nrand@ 'Double' values, uniformly distributed in | ||
217 | -- @[0, 1)@. It should use these values to generate its new | ||
218 | -- configuration. | ||
219 | simanSolve :: Int -- ^ Seed for the random number generator | ||
220 | -> Int -- ^ @nrand@, the number of random 'Double's the | ||
221 | -- step function requires | ||
222 | -> SimulatedAnnealingParams -- ^ Parameters to configure the solver | ||
223 | -> a -- ^ Initial configuration @x0@ | ||
224 | -> (a -> Double) -- ^ Energy functional @e@ | ||
225 | -> (a -> a -> Double) -- ^ Metric definition @m@ | ||
226 | -> (Vector Double -> Double -> a -> a) -- ^ Stepping function @step@ | ||
227 | -> Maybe (a -> String) -- ^ Optional printing function | ||
228 | -> a -- ^ Best configuration the solver has found | ||
229 | simanSolve seed nrand params conf e m step printfun = | ||
230 | unsafePerformIO $ with params $ \paramptr -> do | ||
231 | ewrap <- mkEnergyFun $ wrapEnergy e | ||
232 | mwrap <- mkMetricFun $ wrapMetric m | ||
233 | stepwrap <- mkStepFun $ wrapStep nrand step | ||
234 | confptr <- newIORef conf >>= newStablePtr | ||
235 | cpwrap <- mkCopyFun copyConfig | ||
236 | ccwrap <- mkCopyConstructorFun copyConstructConfig | ||
237 | dwrap <- mkDestructFun destroyConfig | ||
238 | pwrap <- case printfun of | ||
239 | Nothing -> return nullFunPtr | ||
240 | Just pf -> mkDestructFun $ wrapPrint pf | ||
241 | siman (fromIntegral seed) | ||
242 | paramptr confptr | ||
243 | ewrap mwrap stepwrap cpwrap ccwrap dwrap pwrap // check "siman" | ||
244 | result <- deRefRead confptr | ||
245 | freeStablePtr confptr | ||
246 | return result | ||
diff --git a/packages/gsl/src/Numeric/GSL/gsl-aux.c b/packages/gsl/src/Numeric/GSL/gsl-aux.c index e1b189c..1ca8199 100644 --- a/packages/gsl/src/Numeric/GSL/gsl-aux.c +++ b/packages/gsl/src/Numeric/GSL/gsl-aux.c | |||
@@ -36,6 +36,8 @@ | |||
36 | #include <gsl/gsl_roots.h> | 36 | #include <gsl/gsl_roots.h> |
37 | #include <gsl/gsl_spline.h> | 37 | #include <gsl/gsl_spline.h> |
38 | #include <gsl/gsl_multifit_nlin.h> | 38 | #include <gsl/gsl_multifit_nlin.h> |
39 | #include <gsl/gsl_siman.h> | ||
40 | |||
39 | #include <string.h> | 41 | #include <string.h> |
40 | #include <stdio.h> | 42 | #include <stdio.h> |
41 | 43 | ||
@@ -475,7 +477,30 @@ int uniMinimize(int method, double f(double), | |||
475 | OK | 477 | OK |
476 | } | 478 | } |
477 | 479 | ||
478 | 480 | int siman(int seed, | |
481 | gsl_siman_params_t *params, void *xp0, | ||
482 | double energy(void *), double metric(void *, void *), | ||
483 | void step(const gsl_rng *, void *, double), | ||
484 | void copy(void *, void *), void *copycons(void *), | ||
485 | void destroy(void *), void print(void *)) { | ||
486 | DEBUGMSG("siman"); | ||
487 | gsl_rng *gen = gsl_rng_alloc (gsl_rng_mt19937); | ||
488 | gsl_rng_set(gen, seed); | ||
489 | |||
490 | // The simulated annealing routine doesn't indicate with a return | ||
491 | // code how things went -- there's little notion of convergence for | ||
492 | // a randomized minimizer on a potentially non-convex problem, and I | ||
493 | // suppose it doesn't detect egregious failures like malloc errors | ||
494 | // in the copy-constructor. | ||
495 | gsl_siman_solve(gen, xp0, | ||
496 | energy, step, | ||
497 | metric, print, | ||
498 | copy, copycons, | ||
499 | destroy, 0, *params); | ||
500 | |||
501 | gsl_rng_free(gen); | ||
502 | OK | ||
503 | } | ||
479 | 504 | ||
480 | // this version returns info about intermediate steps | 505 | // this version returns info about intermediate steps |
481 | int minimize(int method, double f(int, double*), double tolsize, int maxit, | 506 | int minimize(int method, double f(int, double*), double tolsize, int maxit, |