diff options
Diffstat (limited to 'packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs')
-rw-r--r-- | packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs | 245 |
1 files changed, 245 insertions, 0 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..11b22d3 --- /dev/null +++ b/packages/gsl/src/Numeric/GSL/SimulatedAnnealing.hs | |||
@@ -0,0 +1,245 @@ | |||
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 | |||
59 | import System.IO (hFlush, stdout) | ||
60 | |||
61 | import Data.IORef (IORef, newIORef, writeIORef, readIORef, modifyIORef') | ||
62 | |||
63 | -- | 'SimulatedAnnealingParams' is a translation of the | ||
64 | -- @gsl_siman_params_t@ structure documented in | ||
65 | -- <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing-functions.html#Simulated-Annealing-functions the GSL manual>, | ||
66 | -- which controls the simulated annealing algorithm. | ||
67 | -- | ||
68 | -- The annealing process is parameterized by the Boltzmann | ||
69 | -- distribution and the /cooling schedule/. For more details, see | ||
70 | -- <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing-algorithm.html#Simulated-Annealing-algorithm the relevant section of the manual>. | ||
71 | data SimulatedAnnealingParams = SimulatedAnnealingParams { | ||
72 | n_tries :: CInt -- ^ The number of points to try for each step. | ||
73 | , iters_fixed_T :: CInt -- ^ The number of iterations at each temperature | ||
74 | , step_size :: Double -- ^ The maximum step size in the random walk | ||
75 | , boltzmann_k :: Double -- ^ Boltzmann distribution parameter | ||
76 | , cooling_t_initial :: Double -- ^ Initial temperature | ||
77 | , cooling_mu_t :: Double -- ^ Cooling rate parameter | ||
78 | , cooling_t_min :: Double -- ^ Final temperature | ||
79 | } deriving (Eq, Show, Read) | ||
80 | |||
81 | instance Storable SimulatedAnnealingParams where | ||
82 | sizeOf p = sizeOf (n_tries p) + | ||
83 | sizeOf (iters_fixed_T p) + | ||
84 | sizeOf (step_size p) + | ||
85 | sizeOf (boltzmann_k p) + | ||
86 | sizeOf (cooling_t_initial p) + | ||
87 | sizeOf (cooling_mu_t p) + | ||
88 | sizeOf (cooling_t_min p) | ||
89 | -- TODO(MP): is this safe? | ||
90 | alignment p = alignment (step_size p) | ||
91 | -- TODO(MP): Is there a more automatic way to write these? | ||
92 | peek ptr = SimulatedAnnealingParams <$> | ||
93 | peekByteOff ptr 0 <*> | ||
94 | peekByteOff ptr i <*> | ||
95 | peekByteOff ptr (2*i) <*> | ||
96 | peekByteOff ptr (2*i + d) <*> | ||
97 | peekByteOff ptr (2*i + 2*d) <*> | ||
98 | peekByteOff ptr (2*i + 3*d) <*> | ||
99 | peekByteOff ptr (2*i + 4*d) | ||
100 | where | ||
101 | i = sizeOf (0 :: CInt) | ||
102 | d = sizeOf (0 :: Double) | ||
103 | poke ptr sap = do | ||
104 | pokeByteOff ptr 0 (n_tries sap) | ||
105 | pokeByteOff ptr i (iters_fixed_T sap) | ||
106 | pokeByteOff ptr (2*i) (step_size sap) | ||
107 | pokeByteOff ptr (2*i + d) (boltzmann_k sap) | ||
108 | pokeByteOff ptr (2*i + 2*d) (cooling_t_initial sap) | ||
109 | pokeByteOff ptr (2*i + 3*d) (cooling_mu_t sap) | ||
110 | pokeByteOff ptr (2*i + 4*d) (cooling_t_min sap) | ||
111 | where | ||
112 | i = sizeOf (0 :: CInt) | ||
113 | d = sizeOf (0 :: Double) | ||
114 | |||
115 | -- We use a StablePtr to an IORef so that we can keep hold of | ||
116 | -- StablePtr values but mutate their contents. A simple 'StablePtr a' | ||
117 | -- won't work, since we'd have no way to write 'copyConfig'. | ||
118 | type P a = StablePtr (IORef a) | ||
119 | |||
120 | copyConfig :: P a -> P a -> IO () | ||
121 | copyConfig src' dest' = do | ||
122 | dest <- deRefStablePtr dest' | ||
123 | src <- deRefStablePtr src' | ||
124 | readIORef src >>= writeIORef dest | ||
125 | |||
126 | copyConstructConfig :: P a -> IO (P a) | ||
127 | copyConstructConfig x = do | ||
128 | conf <- deRefRead x | ||
129 | newconf <- newIORef conf | ||
130 | newStablePtr newconf | ||
131 | |||
132 | destroyConfig :: P a -> IO () | ||
133 | destroyConfig p = do | ||
134 | freeStablePtr p | ||
135 | |||
136 | deRefRead :: P a -> IO a | ||
137 | deRefRead p = deRefStablePtr p >>= readIORef | ||
138 | |||
139 | wrapEnergy :: (a -> Double) -> P a -> Double | ||
140 | wrapEnergy f p = unsafePerformIO $ f <$> deRefRead p | ||
141 | |||
142 | wrapMetric :: (a -> a -> Double) -> P a -> P a -> Double | ||
143 | wrapMetric f x y = unsafePerformIO $ f <$> deRefRead x <*> deRefRead y | ||
144 | |||
145 | wrapStep :: Int | ||
146 | -> (Vector Double -> Double -> a -> a) | ||
147 | -> GSLRNG | ||
148 | -> P a | ||
149 | -> Double | ||
150 | -> IO () | ||
151 | wrapStep nrand f (GSLRNG rng) confptr stepSize = do | ||
152 | v <- generateM nrand (\_ -> gslRngUniform rng) | ||
153 | conf <- deRefStablePtr confptr | ||
154 | modifyIORef' conf $ f v stepSize | ||
155 | |||
156 | wrapPrint :: (a -> String) -> P a -> IO () | ||
157 | wrapPrint pf ptr = deRefRead ptr >>= putStr . pf >> hFlush stdout | ||
158 | |||
159 | foreign import ccall safe "wrapper" | ||
160 | mkEnergyFun :: (P a -> Double) -> IO (FunPtr (P a -> Double)) | ||
161 | |||
162 | foreign import ccall safe "wrapper" | ||
163 | mkMetricFun :: (P a -> P a -> Double) -> IO (FunPtr (P a -> P a -> Double)) | ||
164 | |||
165 | foreign import ccall safe "wrapper" | ||
166 | mkStepFun :: (GSLRNG -> P a -> Double -> IO ()) | ||
167 | -> IO (FunPtr (GSLRNG -> P a -> Double -> IO ())) | ||
168 | |||
169 | foreign import ccall safe "wrapper" | ||
170 | mkCopyFun :: (P a -> P a -> IO ()) -> IO (FunPtr (P a -> P a -> IO ())) | ||
171 | |||
172 | foreign import ccall safe "wrapper" | ||
173 | mkCopyConstructorFun :: (P a -> IO (P a)) -> IO (FunPtr (P a -> IO (P a))) | ||
174 | |||
175 | foreign import ccall safe "wrapper" | ||
176 | mkDestructFun :: (P a -> IO ()) -> IO (FunPtr (P a -> IO ())) | ||
177 | |||
178 | newtype GSLRNG = GSLRNG (Ptr GSLRNG) | ||
179 | |||
180 | foreign import ccall safe "gsl_rng.h gsl_rng_uniform" | ||
181 | gslRngUniform :: Ptr GSLRNG -> IO Double | ||
182 | |||
183 | foreign import ccall safe "gsl-aux.h siman" | ||
184 | siman :: CInt -- ^ RNG seed (for repeatability) | ||
185 | -> Ptr SimulatedAnnealingParams -- ^ params | ||
186 | -> P a -- ^ Configuration | ||
187 | -> FunPtr (P a -> Double) -- ^ Energy functional | ||
188 | -> FunPtr (P a -> P a -> Double) -- ^ Metric definition | ||
189 | -> FunPtr (GSLRNG -> P a -> Double -> IO ()) -- ^ Step evaluation | ||
190 | -> FunPtr (P a -> P a -> IO ()) -- ^ Copy config | ||
191 | -> FunPtr (P a -> IO (P a)) -- ^ Copy constructor for config | ||
192 | -> FunPtr (P a -> IO ()) -- ^ Destructor for config | ||
193 | -> FunPtr (P a -> IO ()) -- ^ Print function | ||
194 | -> IO CInt | ||
195 | |||
196 | -- | | ||
197 | -- Calling | ||
198 | -- | ||
199 | -- > simanSolve seed nrand params x0 e m step print | ||
200 | -- | ||
201 | -- performs a simulated annealing search through a given space. So | ||
202 | -- that any configuration type may be used, the space is specified by | ||
203 | -- providing the functions @e@ (the energy functional) and @m@ (the | ||
204 | -- metric definition). @x0@ is the initial configuration of the | ||
205 | -- system. The simulated annealing steps are generated using the | ||
206 | -- user-provided function @step@, which should randomly construct a | ||
207 | -- new system configuration. | ||
208 | -- | ||
209 | -- If 'Nothing' is passed instead of a printing function, no | ||
210 | -- incremental output will be generated. Otherwise, the GSL-formatted | ||
211 | -- output, including the configuration description the user function | ||
212 | -- generates, will be printed to stdout. | ||
213 | -- | ||
214 | -- Each time the step function is called, it is supplied with a random | ||
215 | -- vector containing @nrand@ 'Double' values, uniformly distributed in | ||
216 | -- @[0, 1)@. It should use these values to generate its new | ||
217 | -- configuration. | ||
218 | simanSolve :: Int -- ^ Seed for the random number generator | ||
219 | -> Int -- ^ @nrand@, the number of random 'Double's the | ||
220 | -- step function requires | ||
221 | -> SimulatedAnnealingParams -- ^ Parameters to configure the solver | ||
222 | -> a -- ^ Initial configuration @x0@ | ||
223 | -> (a -> Double) -- ^ Energy functional @e@ | ||
224 | -> (a -> a -> Double) -- ^ Metric definition @m@ | ||
225 | -> (Vector Double -> Double -> a -> a) -- ^ Stepping function @step@ | ||
226 | -> Maybe (a -> String) -- ^ Optional printing function | ||
227 | -> a -- ^ Best configuration the solver has found | ||
228 | simanSolve seed nrand params conf e m step printfun = | ||
229 | unsafePerformIO $ with params $ \paramptr -> do | ||
230 | ewrap <- mkEnergyFun $ wrapEnergy e | ||
231 | mwrap <- mkMetricFun $ wrapMetric m | ||
232 | stepwrap <- mkStepFun $ wrapStep nrand step | ||
233 | confptr <- newIORef conf >>= newStablePtr | ||
234 | cpwrap <- mkCopyFun copyConfig | ||
235 | ccwrap <- mkCopyConstructorFun copyConstructConfig | ||
236 | dwrap <- mkDestructFun destroyConfig | ||
237 | pwrap <- case printfun of | ||
238 | Nothing -> return nullFunPtr | ||
239 | Just pf -> mkDestructFun $ wrapPrint pf | ||
240 | siman (fromIntegral seed) | ||
241 | paramptr confptr | ||
242 | ewrap mwrap stepwrap cpwrap ccwrap dwrap pwrap // check "siman" | ||
243 | result <- deRefRead confptr | ||
244 | freeStablePtr confptr | ||
245 | return result | ||