summaryrefslogtreecommitdiff
path: root/packages/gsl
diff options
context:
space:
mode:
Diffstat (limited to 'packages/gsl')
-rw-r--r--packages/gsl/hmatrix-gsl.cabal1
-rw-r--r--packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs246
-rw-r--r--packages/gsl/src/Numeric/GSL/gsl-aux.c27
3 files changed, 273 insertions, 1 deletions
diff --git a/packages/gsl/hmatrix-gsl.cabal b/packages/gsl/hmatrix-gsl.cabal
index 6f983f2..f288c64 100644
--- a/packages/gsl/hmatrix-gsl.cabal
+++ b/packages/gsl/hmatrix-gsl.cabal
@@ -44,6 +44,7 @@ library
44 Numeric.GSL, 44 Numeric.GSL,
45 Numeric.GSL.LinearAlgebra, 45 Numeric.GSL.LinearAlgebra,
46 Numeric.GSL.Interpolation, 46 Numeric.GSL.Interpolation,
47 Numeric.GSL.SimulatedAnnealing,
47 Graphics.Plot 48 Graphics.Plot
48 other-modules: Numeric.GSL.Internal, 49 other-modules: Numeric.GSL.Internal,
49 Numeric.GSL.Vector, 50 Numeric.GSL.Vector,
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{- |
2Module : Numeric.GSL.Interpolation
3Copyright : (c) Matthew Peddie 2015
4License : GPL
5Maintainer : Alberto Ruiz
6Stability : provisional
7
8Simulated annealing routines.
9
10<https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing.html#Simulated-Annealing>
11
12Here 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
28The 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
36This global minimum is around 1.36.
37
38-}
39{-# OPTIONS_GHC -Wall #-}
40
41module Numeric.GSL.SimulatedAnnealing (
42 -- * Searching for minima
43 simanSolve
44 -- * Configuring the annealing process
45 , SimulatedAnnealingParams(..)
46 ) where
47
48import Numeric.GSL.Internal
49import Numeric.LinearAlgebra.HMatrix hiding(step)
50
51import Data.Vector.Storable(generateM)
52import Foreign.Storable(Storable(..))
53import Foreign.Marshal.Utils(with)
54import Foreign.Ptr(Ptr, FunPtr, nullFunPtr)
55import Foreign.StablePtr(StablePtr, newStablePtr, deRefStablePtr, freeStablePtr)
56import Foreign.C.Types
57import System.IO.Unsafe(unsafePerformIO)
58import Control.Applicative ((<*>), (<$>))
59
60import System.IO (hFlush, stdout)
61
62import 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>.
72data 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
82instance 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'.
119type P a = StablePtr (IORef a)
120
121copyConfig :: P a -> P a -> IO ()
122copyConfig src' dest' = do
123 dest <- deRefStablePtr dest'
124 src <- deRefStablePtr src'
125 readIORef src >>= writeIORef dest
126
127copyConstructConfig :: P a -> IO (P a)
128copyConstructConfig x = do
129 conf <- deRefRead x
130 newconf <- newIORef conf
131 newStablePtr newconf
132
133destroyConfig :: P a -> IO ()
134destroyConfig p = do
135 freeStablePtr p
136
137deRefRead :: P a -> IO a
138deRefRead p = deRefStablePtr p >>= readIORef
139
140wrapEnergy :: (a -> Double) -> P a -> Double
141wrapEnergy f p = unsafePerformIO $ f <$> deRefRead p
142
143wrapMetric :: (a -> a -> Double) -> P a -> P a -> Double
144wrapMetric f x y = unsafePerformIO $ f <$> deRefRead x <*> deRefRead y
145
146wrapStep :: Int
147 -> (Vector Double -> Double -> a -> a)
148 -> GSLRNG
149 -> P a
150 -> Double
151 -> IO ()
152wrapStep nrand f (GSLRNG rng) confptr stepSize = do
153 v <- generateM nrand (\_ -> gslRngUniform rng)
154 conf <- deRefStablePtr confptr
155 modifyIORef' conf $ f v stepSize
156
157wrapPrint :: (a -> String) -> P a -> IO ()
158wrapPrint pf ptr = deRefRead ptr >>= putStr . pf >> hFlush stdout
159
160foreign import ccall safe "wrapper"
161 mkEnergyFun :: (P a -> Double) -> IO (FunPtr (P a -> Double))
162
163foreign import ccall safe "wrapper"
164 mkMetricFun :: (P a -> P a -> Double) -> IO (FunPtr (P a -> P a -> Double))
165
166foreign import ccall safe "wrapper"
167 mkStepFun :: (GSLRNG -> P a -> Double -> IO ())
168 -> IO (FunPtr (GSLRNG -> P a -> Double -> IO ()))
169
170foreign import ccall safe "wrapper"
171 mkCopyFun :: (P a -> P a -> IO ()) -> IO (FunPtr (P a -> P a -> IO ()))
172
173foreign import ccall safe "wrapper"
174 mkCopyConstructorFun :: (P a -> IO (P a)) -> IO (FunPtr (P a -> IO (P a)))
175
176foreign import ccall safe "wrapper"
177 mkDestructFun :: (P a -> IO ()) -> IO (FunPtr (P a -> IO ()))
178
179newtype GSLRNG = GSLRNG (Ptr GSLRNG)
180
181foreign import ccall safe "gsl_rng.h gsl_rng_uniform"
182 gslRngUniform :: Ptr GSLRNG -> IO Double
183
184foreign 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.
219simanSolve :: 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
229simanSolve 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 480int 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
481int minimize(int method, double f(int, double*), double tolsize, int maxit, 506int minimize(int method, double f(int, double*), double tolsize, int maxit,