{- |
Module      :  Numeric.GSL.Interpolation
Copyright   :  (c) Matthew Peddie 2015
License     :  GPL
Maintainer  :  Alberto Ruiz
Stability   :  provisional

Simulated annealing routines.

<https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing.html#Simulated-Annealing>

Here is a translation of the simple example given in
<https://www.gnu.org/software/gsl/manual/html_node/Trivial-example.html#Trivial-example the GSL manual>:

> import Numeric.GSL.SimulatedAnnealing
> import Numeric.LinearAlgebra.HMatrix
>
> main = print $ simanSolve 0 1 exampleParams 15.5 exampleE exampleM exampleS (Just show)
>
> exampleParams = SimulatedAnnealingParams 200 1000 1.0 1.0 0.008 1.003 2.0e-6
>
> exampleE x = exp (-(x - 1)**2) * sin (8 * x)
>
> exampleM x y = abs $ x - y
>
> exampleS rands stepSize current = (rands ! 0) * 2 * stepSize - stepSize + current

The manual states:

>     The first example, in one dimensional Cartesian space, sets up an
>     energy function which is a damped sine wave; this has many local
>     minima, but only one global minimum, somewhere between 1.0 and
>     1.5. The initial guess given is 15.5, which is several local minima
>     away from the global minimum.

This global minimum is around 1.36.

-}
{-# OPTIONS_GHC -Wall #-}

module Numeric.GSL.SimulatedAnnealing (
  -- * Searching for minima
  simanSolve
  -- * Configuring the annealing process
  , SimulatedAnnealingParams(..)
  ) where

import Numeric.GSL.Internal
import Numeric.LinearAlgebra.HMatrix hiding(step)

import Data.Vector.Storable(generateM)
import Foreign.Storable(Storable(..))
import Foreign.Marshal.Utils(with)
import Foreign.Ptr(Ptr, FunPtr, nullFunPtr)
import Foreign.StablePtr(StablePtr, newStablePtr, deRefStablePtr, freeStablePtr)
import Foreign.C.Types
import System.IO.Unsafe(unsafePerformIO)

import System.IO (hFlush, stdout)

import Data.IORef (IORef, newIORef, writeIORef, readIORef, modifyIORef')

-- | 'SimulatedAnnealingParams' is a translation of the
-- @gsl_siman_params_t@ structure documented in
-- <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing-functions.html#Simulated-Annealing-functions the GSL manual>,
-- which controls the simulated annealing algorithm.
--
-- The annealing process is parameterized by the Boltzmann
-- distribution and the /cooling schedule/.  For more details, see
-- <https://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing-algorithm.html#Simulated-Annealing-algorithm the relevant section of the manual>.
data SimulatedAnnealingParams = SimulatedAnnealingParams {
  SimulatedAnnealingParams -> CInt
n_tries :: CInt  -- ^ The number of points to try for each step.
  , SimulatedAnnealingParams -> CInt
iters_fixed_T :: CInt  -- ^ The number of iterations at each temperature
  , SimulatedAnnealingParams -> Double
step_size :: Double    -- ^ The maximum step size in the random walk
  , SimulatedAnnealingParams -> Double
boltzmann_k :: Double  -- ^ Boltzmann distribution parameter
  , SimulatedAnnealingParams -> Double
cooling_t_initial :: Double -- ^ Initial temperature
  , SimulatedAnnealingParams -> Double
cooling_mu_t :: Double      -- ^ Cooling rate parameter
  , SimulatedAnnealingParams -> Double
cooling_t_min :: Double     -- ^ Final temperature
  } deriving (SimulatedAnnealingParams -> SimulatedAnnealingParams -> Bool
(SimulatedAnnealingParams -> SimulatedAnnealingParams -> Bool)
-> (SimulatedAnnealingParams -> SimulatedAnnealingParams -> Bool)
-> Eq SimulatedAnnealingParams
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: SimulatedAnnealingParams -> SimulatedAnnealingParams -> Bool
$c/= :: SimulatedAnnealingParams -> SimulatedAnnealingParams -> Bool
== :: SimulatedAnnealingParams -> SimulatedAnnealingParams -> Bool
$c== :: SimulatedAnnealingParams -> SimulatedAnnealingParams -> Bool
Eq, Int -> SimulatedAnnealingParams -> ShowS
[SimulatedAnnealingParams] -> ShowS
SimulatedAnnealingParams -> String
(Int -> SimulatedAnnealingParams -> ShowS)
-> (SimulatedAnnealingParams -> String)
-> ([SimulatedAnnealingParams] -> ShowS)
-> Show SimulatedAnnealingParams
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [SimulatedAnnealingParams] -> ShowS
$cshowList :: [SimulatedAnnealingParams] -> ShowS
show :: SimulatedAnnealingParams -> String
$cshow :: SimulatedAnnealingParams -> String
showsPrec :: Int -> SimulatedAnnealingParams -> ShowS
$cshowsPrec :: Int -> SimulatedAnnealingParams -> ShowS
Show, ReadPrec [SimulatedAnnealingParams]
ReadPrec SimulatedAnnealingParams
Int -> ReadS SimulatedAnnealingParams
ReadS [SimulatedAnnealingParams]
(Int -> ReadS SimulatedAnnealingParams)
-> ReadS [SimulatedAnnealingParams]
-> ReadPrec SimulatedAnnealingParams
-> ReadPrec [SimulatedAnnealingParams]
-> Read SimulatedAnnealingParams
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [SimulatedAnnealingParams]
$creadListPrec :: ReadPrec [SimulatedAnnealingParams]
readPrec :: ReadPrec SimulatedAnnealingParams
$creadPrec :: ReadPrec SimulatedAnnealingParams
readList :: ReadS [SimulatedAnnealingParams]
$creadList :: ReadS [SimulatedAnnealingParams]
readsPrec :: Int -> ReadS SimulatedAnnealingParams
$creadsPrec :: Int -> ReadS SimulatedAnnealingParams
Read)

instance Storable SimulatedAnnealingParams where
  sizeOf :: SimulatedAnnealingParams -> Int
sizeOf p :: SimulatedAnnealingParams
p = CInt -> Int
forall a. Storable a => a -> Int
sizeOf (SimulatedAnnealingParams -> CInt
n_tries SimulatedAnnealingParams
p) Int -> Int -> Int
forall a. Num a => a -> a -> a
+
             CInt -> Int
forall a. Storable a => a -> Int
sizeOf (SimulatedAnnealingParams -> CInt
iters_fixed_T SimulatedAnnealingParams
p) Int -> Int -> Int
forall a. Num a => a -> a -> a
+
             Double -> Int
forall a. Storable a => a -> Int
sizeOf (SimulatedAnnealingParams -> Double
step_size SimulatedAnnealingParams
p) Int -> Int -> Int
forall a. Num a => a -> a -> a
+
             Double -> Int
forall a. Storable a => a -> Int
sizeOf (SimulatedAnnealingParams -> Double
boltzmann_k SimulatedAnnealingParams
p) Int -> Int -> Int
forall a. Num a => a -> a -> a
+
             Double -> Int
forall a. Storable a => a -> Int
sizeOf (SimulatedAnnealingParams -> Double
cooling_t_initial SimulatedAnnealingParams
p) Int -> Int -> Int
forall a. Num a => a -> a -> a
+
             Double -> Int
forall a. Storable a => a -> Int
sizeOf (SimulatedAnnealingParams -> Double
cooling_mu_t SimulatedAnnealingParams
p) Int -> Int -> Int
forall a. Num a => a -> a -> a
+
             Double -> Int
forall a. Storable a => a -> Int
sizeOf (SimulatedAnnealingParams -> Double
cooling_t_min SimulatedAnnealingParams
p)
  -- TODO(MP): is this safe?
  alignment :: SimulatedAnnealingParams -> Int
alignment p :: SimulatedAnnealingParams
p = Double -> Int
forall a. Storable a => a -> Int
alignment (SimulatedAnnealingParams -> Double
step_size SimulatedAnnealingParams
p)
  -- TODO(MP): Is there a more automatic way to write these?
  peek :: Ptr SimulatedAnnealingParams -> IO SimulatedAnnealingParams
peek ptr :: Ptr SimulatedAnnealingParams
ptr = CInt
-> CInt
-> Double
-> Double
-> Double
-> Double
-> Double
-> SimulatedAnnealingParams
SimulatedAnnealingParams (CInt
 -> CInt
 -> Double
 -> Double
 -> Double
 -> Double
 -> Double
 -> SimulatedAnnealingParams)
-> IO CInt
-> IO
     (CInt
      -> Double
      -> Double
      -> Double
      -> Double
      -> Double
      -> SimulatedAnnealingParams)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$>
             Ptr SimulatedAnnealingParams -> Int -> IO CInt
forall a b. Storable a => Ptr b -> Int -> IO a
peekByteOff Ptr SimulatedAnnealingParams
ptr 0 IO
  (CInt
   -> Double
   -> Double
   -> Double
   -> Double
   -> Double
   -> SimulatedAnnealingParams)
-> IO CInt
-> IO
     (Double
      -> Double
      -> Double
      -> Double
      -> Double
      -> SimulatedAnnealingParams)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*>
             Ptr SimulatedAnnealingParams -> Int -> IO CInt
forall a b. Storable a => Ptr b -> Int -> IO a
peekByteOff Ptr SimulatedAnnealingParams
ptr Int
i IO
  (Double
   -> Double
   -> Double
   -> Double
   -> Double
   -> SimulatedAnnealingParams)
-> IO Double
-> IO
     (Double -> Double -> Double -> Double -> SimulatedAnnealingParams)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*>
             Ptr SimulatedAnnealingParams -> Int -> IO Double
forall a b. Storable a => Ptr b -> Int -> IO a
peekByteOff Ptr SimulatedAnnealingParams
ptr (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i) IO
  (Double -> Double -> Double -> Double -> SimulatedAnnealingParams)
-> IO Double
-> IO (Double -> Double -> Double -> SimulatedAnnealingParams)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*>
             Ptr SimulatedAnnealingParams -> Int -> IO Double
forall a b. Storable a => Ptr b -> Int -> IO a
peekByteOff Ptr SimulatedAnnealingParams
ptr (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d) IO (Double -> Double -> Double -> SimulatedAnnealingParams)
-> IO Double -> IO (Double -> Double -> SimulatedAnnealingParams)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*>
             Ptr SimulatedAnnealingParams -> Int -> IO Double
forall a b. Storable a => Ptr b -> Int -> IO a
peekByteOff Ptr SimulatedAnnealingParams
ptr (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
d) IO (Double -> Double -> SimulatedAnnealingParams)
-> IO Double -> IO (Double -> SimulatedAnnealingParams)
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*>
             Ptr SimulatedAnnealingParams -> Int -> IO Double
forall a b. Storable a => Ptr b -> Int -> IO a
peekByteOff Ptr SimulatedAnnealingParams
ptr (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 3Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
d) IO (Double -> SimulatedAnnealingParams)
-> IO Double -> IO SimulatedAnnealingParams
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*>
             Ptr SimulatedAnnealingParams -> Int -> IO Double
forall a b. Storable a => Ptr b -> Int -> IO a
peekByteOff Ptr SimulatedAnnealingParams
ptr (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 4Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
d)
    where
      i :: Int
i = CInt -> Int
forall a. Storable a => a -> Int
sizeOf (0 :: CInt)
      d :: Int
d = Double -> Int
forall a. Storable a => a -> Int
sizeOf (0 :: Double)
  poke :: Ptr SimulatedAnnealingParams -> SimulatedAnnealingParams -> IO ()
poke ptr :: Ptr SimulatedAnnealingParams
ptr sap :: SimulatedAnnealingParams
sap = do
    Ptr SimulatedAnnealingParams -> Int -> CInt -> IO ()
forall a b. Storable a => Ptr b -> Int -> a -> IO ()
pokeByteOff Ptr SimulatedAnnealingParams
ptr 0 (SimulatedAnnealingParams -> CInt
n_tries SimulatedAnnealingParams
sap)
    Ptr SimulatedAnnealingParams -> Int -> CInt -> IO ()
forall a b. Storable a => Ptr b -> Int -> a -> IO ()
pokeByteOff Ptr SimulatedAnnealingParams
ptr Int
i (SimulatedAnnealingParams -> CInt
iters_fixed_T SimulatedAnnealingParams
sap)
    Ptr SimulatedAnnealingParams -> Int -> Double -> IO ()
forall a b. Storable a => Ptr b -> Int -> a -> IO ()
pokeByteOff Ptr SimulatedAnnealingParams
ptr (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i) (SimulatedAnnealingParams -> Double
step_size SimulatedAnnealingParams
sap)
    Ptr SimulatedAnnealingParams -> Int -> Double -> IO ()
forall a b. Storable a => Ptr b -> Int -> a -> IO ()
pokeByteOff Ptr SimulatedAnnealingParams
ptr (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d) (SimulatedAnnealingParams -> Double
boltzmann_k SimulatedAnnealingParams
sap)
    Ptr SimulatedAnnealingParams -> Int -> Double -> IO ()
forall a b. Storable a => Ptr b -> Int -> a -> IO ()
pokeByteOff Ptr SimulatedAnnealingParams
ptr (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
d) (SimulatedAnnealingParams -> Double
cooling_t_initial SimulatedAnnealingParams
sap)
    Ptr SimulatedAnnealingParams -> Int -> Double -> IO ()
forall a b. Storable a => Ptr b -> Int -> a -> IO ()
pokeByteOff Ptr SimulatedAnnealingParams
ptr (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 3Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
d) (SimulatedAnnealingParams -> Double
cooling_mu_t SimulatedAnnealingParams
sap)
    Ptr SimulatedAnnealingParams -> Int -> Double -> IO ()
forall a b. Storable a => Ptr b -> Int -> a -> IO ()
pokeByteOff Ptr SimulatedAnnealingParams
ptr (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 4Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
d) (SimulatedAnnealingParams -> Double
cooling_t_min SimulatedAnnealingParams
sap)
    where
      i :: Int
i = CInt -> Int
forall a. Storable a => a -> Int
sizeOf (0 :: CInt)
      d :: Int
d = Double -> Int
forall a. Storable a => a -> Int
sizeOf (0 :: Double)

-- We use a StablePtr to an IORef so that we can keep hold of
-- StablePtr values but mutate their contents.  A simple 'StablePtr a'
-- won't work, since we'd have no way to write 'copyConfig'.
type P a = StablePtr (IORef a)

copyConfig :: P a -> P a -> IO ()
copyConfig :: P a -> P a -> IO ()
copyConfig src' :: P a
src' dest' :: P a
dest' = do
  IORef a
dest <- P a -> IO (IORef a)
forall a. StablePtr a -> IO a
deRefStablePtr P a
dest'
  IORef a
src <- P a -> IO (IORef a)
forall a. StablePtr a -> IO a
deRefStablePtr P a
src'
  IORef a -> IO a
forall a. IORef a -> IO a
readIORef IORef a
src IO a -> (a -> IO ()) -> IO ()
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= IORef a -> a -> IO ()
forall a. IORef a -> a -> IO ()
writeIORef IORef a
dest

copyConstructConfig :: P a -> IO (P a)
copyConstructConfig :: P a -> IO (P a)
copyConstructConfig x :: P a
x = do
  a
conf <- P a -> IO a
forall a. P a -> IO a
deRefRead P a
x
  IORef a
newconf <- a -> IO (IORef a)
forall a. a -> IO (IORef a)
newIORef a
conf
  IORef a -> IO (P a)
forall a. a -> IO (StablePtr a)
newStablePtr IORef a
newconf

destroyConfig :: P a -> IO ()
destroyConfig :: P a -> IO ()
destroyConfig p :: P a
p = do
  P a -> IO ()
forall a. StablePtr a -> IO ()
freeStablePtr P a
p

deRefRead :: P a -> IO a
deRefRead :: P a -> IO a
deRefRead p :: P a
p = P a -> IO (IORef a)
forall a. StablePtr a -> IO a
deRefStablePtr P a
p IO (IORef a) -> (IORef a -> IO a) -> IO a
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= IORef a -> IO a
forall a. IORef a -> IO a
readIORef

wrapEnergy :: (a -> Double) -> P a -> Double
wrapEnergy :: (a -> Double) -> P a -> Double
wrapEnergy f :: a -> Double
f p :: P a
p = IO Double -> Double
forall a. IO a -> a
unsafePerformIO (IO Double -> Double) -> IO Double -> Double
forall a b. (a -> b) -> a -> b
$ a -> Double
f (a -> Double) -> IO a -> IO Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> P a -> IO a
forall a. P a -> IO a
deRefRead P a
p

wrapMetric :: (a -> a -> Double) -> P a -> P a -> Double
wrapMetric :: (a -> a -> Double) -> P a -> P a -> Double
wrapMetric f :: a -> a -> Double
f x :: P a
x y :: P a
y = IO Double -> Double
forall a. IO a -> a
unsafePerformIO (IO Double -> Double) -> IO Double -> Double
forall a b. (a -> b) -> a -> b
$ a -> a -> Double
f (a -> a -> Double) -> IO a -> IO (a -> Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> P a -> IO a
forall a. P a -> IO a
deRefRead P a
x IO (a -> Double) -> IO a -> IO Double
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> P a -> IO a
forall a. P a -> IO a
deRefRead P a
y

wrapStep :: Int
         -> (Vector Double -> Double -> a -> a)
         -> GSLRNG
         -> P a
         -> Double
         -> IO ()
wrapStep :: Int
-> (Vector Double -> Double -> a -> a)
-> GSLRNG
-> P a
-> Double
-> IO ()
wrapStep nrand :: Int
nrand f :: Vector Double -> Double -> a -> a
f (GSLRNG rng :: Ptr GSLRNG
rng) confptr :: P a
confptr stepSize :: Double
stepSize = do
  Vector Double
v <- Int -> (Int -> IO Double) -> IO (Vector Double)
forall (m :: * -> *) a.
(Monad m, Storable a) =>
Int -> (Int -> m a) -> m (Vector a)
generateM Int
nrand (\_ -> Ptr GSLRNG -> IO Double
gslRngUniform Ptr GSLRNG
rng)
  IORef a
conf <- P a -> IO (IORef a)
forall a. StablePtr a -> IO a
deRefStablePtr P a
confptr
  IORef a -> (a -> a) -> IO ()
forall a. IORef a -> (a -> a) -> IO ()
modifyIORef' IORef a
conf ((a -> a) -> IO ()) -> (a -> a) -> IO ()
forall a b. (a -> b) -> a -> b
$ Vector Double -> Double -> a -> a
f Vector Double
v Double
stepSize

wrapPrint :: (a -> String) -> P a -> IO ()
wrapPrint :: (a -> String) -> P a -> IO ()
wrapPrint pf :: a -> String
pf ptr :: P a
ptr = P a -> IO a
forall a. P a -> IO a
deRefRead P a
ptr IO a -> (a -> IO ()) -> IO ()
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= String -> IO ()
putStr (String -> IO ()) -> (a -> String) -> a -> IO ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. a -> String
pf IO () -> IO () -> IO ()
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Handle -> IO ()
hFlush Handle
stdout

foreign import ccall safe "wrapper"
  mkEnergyFun :: (P a -> Double) -> IO (FunPtr (P a -> Double))

foreign import ccall safe "wrapper"
  mkMetricFun :: (P a -> P a -> Double) -> IO (FunPtr (P a -> P a -> Double))

foreign import ccall safe "wrapper"
  mkStepFun :: (GSLRNG -> P a -> Double -> IO ())
            -> IO (FunPtr (GSLRNG -> P a -> Double -> IO ()))

foreign import ccall safe "wrapper"
  mkCopyFun :: (P a -> P a -> IO ()) -> IO (FunPtr (P a -> P a -> IO ()))

foreign import ccall safe "wrapper"
  mkCopyConstructorFun :: (P a -> IO (P a)) -> IO (FunPtr (P a -> IO (P a)))

foreign import ccall safe "wrapper"
  mkDestructFun :: (P a -> IO ()) -> IO (FunPtr (P a -> IO ()))

newtype GSLRNG = GSLRNG (Ptr GSLRNG)

foreign import ccall safe "gsl_rng.h gsl_rng_uniform"
  gslRngUniform :: Ptr GSLRNG -> IO Double

foreign import ccall safe "gsl-aux.h siman"
  siman :: CInt     -- ^ RNG seed (for repeatability)
        -> Ptr SimulatedAnnealingParams    -- ^ params
        -> P a                             -- ^ Configuration
        -> FunPtr (P a -> Double)          -- ^ Energy functional
        -> FunPtr (P a -> P a -> Double) -- ^ Metric definition
        -> FunPtr (GSLRNG -> P a -> Double -> IO ())  -- ^ Step evaluation
        -> FunPtr (P a -> P a -> IO ())  -- ^ Copy config
        -> FunPtr (P a -> IO (P a))      -- ^ Copy constructor for config
        -> FunPtr (P a -> IO ())           -- ^ Destructor for config
        -> FunPtr (P a -> IO ())           -- ^ Print function
        -> IO CInt

-- |
-- Calling
--
-- > simanSolve seed nrand params x0 e m step print
--
-- performs a simulated annealing search through a given space. So
-- that any configuration type may be used, the space is specified by
-- providing the functions @e@ (the energy functional) and @m@ (the
-- metric definition).  @x0@ is the initial configuration of the
-- system.  The simulated annealing steps are generated using the
-- user-provided function @step@, which should randomly construct a
-- new system configuration.
--
-- If 'Nothing' is passed instead of a printing function, no
-- incremental output will be generated.  Otherwise, the GSL-formatted
-- output, including the configuration description the user function
-- generates, will be printed to stdout.
--
-- Each time the step function is called, it is supplied with a random
-- vector containing @nrand@ 'Double' values, uniformly distributed in
-- @[0, 1)@.  It should use these values to generate its new
-- configuration.
simanSolve :: Int   -- ^ Seed for the random number generator
           -> Int   -- ^ @nrand@, the number of random 'Double's the
                    -- step function requires
           -> SimulatedAnnealingParams  -- ^ Parameters to configure the solver
           -> a                    -- ^ Initial configuration @x0@
           -> (a -> Double)        -- ^ Energy functional @e@
           -> (a -> a -> Double)   -- ^ Metric definition @m@
           -> (Vector Double -> Double -> a -> a)  -- ^ Stepping function @step@
           -> Maybe (a -> String)  -- ^ Optional printing function
           -> a          -- ^ Best configuration the solver has found
simanSolve :: Int
-> Int
-> SimulatedAnnealingParams
-> a
-> (a -> Double)
-> (a -> a -> Double)
-> (Vector Double -> Double -> a -> a)
-> Maybe (a -> String)
-> a
simanSolve seed :: Int
seed nrand :: Int
nrand params :: SimulatedAnnealingParams
params conf :: a
conf e :: a -> Double
e m :: a -> a -> Double
m step :: Vector Double -> Double -> a -> a
step printfun :: Maybe (a -> String)
printfun =
  IO a -> a
forall a. IO a -> a
unsafePerformIO (IO a -> a) -> IO a -> a
forall a b. (a -> b) -> a -> b
$ SimulatedAnnealingParams
-> (Ptr SimulatedAnnealingParams -> IO a) -> IO a
forall a b. Storable a => a -> (Ptr a -> IO b) -> IO b
with SimulatedAnnealingParams
params ((Ptr SimulatedAnnealingParams -> IO a) -> IO a)
-> (Ptr SimulatedAnnealingParams -> IO a) -> IO a
forall a b. (a -> b) -> a -> b
$ \paramptr :: Ptr SimulatedAnnealingParams
paramptr -> do
    FunPtr (P a -> Double)
ewrap <- (P a -> Double) -> IO (FunPtr (P a -> Double))
forall a. (P a -> Double) -> IO (FunPtr (P a -> Double))
mkEnergyFun ((P a -> Double) -> IO (FunPtr (P a -> Double)))
-> (P a -> Double) -> IO (FunPtr (P a -> Double))
forall a b. (a -> b) -> a -> b
$ (a -> Double) -> P a -> Double
forall a. (a -> Double) -> P a -> Double
wrapEnergy a -> Double
e
    FunPtr (P a -> P a -> Double)
mwrap <- (P a -> P a -> Double) -> IO (FunPtr (P a -> P a -> Double))
forall a.
(P a -> P a -> Double) -> IO (FunPtr (P a -> P a -> Double))
mkMetricFun ((P a -> P a -> Double) -> IO (FunPtr (P a -> P a -> Double)))
-> (P a -> P a -> Double) -> IO (FunPtr (P a -> P a -> Double))
forall a b. (a -> b) -> a -> b
$ (a -> a -> Double) -> P a -> P a -> Double
forall a. (a -> a -> Double) -> P a -> P a -> Double
wrapMetric a -> a -> Double
m
    FunPtr (GSLRNG -> P a -> Double -> IO ())
stepwrap <- (GSLRNG -> P a -> Double -> IO ())
-> IO (FunPtr (GSLRNG -> P a -> Double -> IO ()))
forall a.
(GSLRNG -> P a -> Double -> IO ())
-> IO (FunPtr (GSLRNG -> P a -> Double -> IO ()))
mkStepFun ((GSLRNG -> P a -> Double -> IO ())
 -> IO (FunPtr (GSLRNG -> P a -> Double -> IO ())))
-> (GSLRNG -> P a -> Double -> IO ())
-> IO (FunPtr (GSLRNG -> P a -> Double -> IO ()))
forall a b. (a -> b) -> a -> b
$ Int
-> (Vector Double -> Double -> a -> a)
-> GSLRNG
-> P a
-> Double
-> IO ()
forall a.
Int
-> (Vector Double -> Double -> a -> a)
-> GSLRNG
-> P a
-> Double
-> IO ()
wrapStep Int
nrand Vector Double -> Double -> a -> a
step
    P a
confptr <- a -> IO (IORef a)
forall a. a -> IO (IORef a)
newIORef a
conf IO (IORef a) -> (IORef a -> IO (P a)) -> IO (P a)
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= IORef a -> IO (P a)
forall a. a -> IO (StablePtr a)
newStablePtr
    FunPtr (P a -> P a -> IO ())
cpwrap <- (P a -> P a -> IO ()) -> IO (FunPtr (P a -> P a -> IO ()))
forall a.
(P a -> P a -> IO ()) -> IO (FunPtr (P a -> P a -> IO ()))
mkCopyFun P a -> P a -> IO ()
forall a. P a -> P a -> IO ()
copyConfig
    FunPtr (P a -> IO (P a))
ccwrap <- (P a -> IO (P a)) -> IO (FunPtr (P a -> IO (P a)))
forall a. (P a -> IO (P a)) -> IO (FunPtr (P a -> IO (P a)))
mkCopyConstructorFun P a -> IO (P a)
forall a. P a -> IO (P a)
copyConstructConfig
    FunPtr (P a -> IO ())
dwrap <- (P a -> IO ()) -> IO (FunPtr (P a -> IO ()))
forall a. (P a -> IO ()) -> IO (FunPtr (P a -> IO ()))
mkDestructFun P a -> IO ()
forall a. P a -> IO ()
destroyConfig
    FunPtr (P a -> IO ())
pwrap <- case Maybe (a -> String)
printfun of
      Nothing -> FunPtr (P a -> IO ()) -> IO (FunPtr (P a -> IO ()))
forall (m :: * -> *) a. Monad m => a -> m a
return FunPtr (P a -> IO ())
forall a. FunPtr a
nullFunPtr
      Just pf :: a -> String
pf -> (P a -> IO ()) -> IO (FunPtr (P a -> IO ()))
forall a. (P a -> IO ()) -> IO (FunPtr (P a -> IO ()))
mkDestructFun ((P a -> IO ()) -> IO (FunPtr (P a -> IO ())))
-> (P a -> IO ()) -> IO (FunPtr (P a -> IO ()))
forall a b. (a -> b) -> a -> b
$ (a -> String) -> P a -> IO ()
forall a. (a -> String) -> P a -> IO ()
wrapPrint a -> String
pf
    CInt
-> Ptr SimulatedAnnealingParams
-> P a
-> FunPtr (P a -> Double)
-> FunPtr (P a -> P a -> Double)
-> FunPtr (GSLRNG -> P a -> Double -> IO ())
-> FunPtr (P a -> P a -> IO ())
-> FunPtr (P a -> IO (P a))
-> FunPtr (P a -> IO ())
-> FunPtr (P a -> IO ())
-> IO CInt
forall a.
CInt
-> Ptr SimulatedAnnealingParams
-> P a
-> FunPtr (P a -> Double)
-> FunPtr (P a -> P a -> Double)
-> FunPtr (GSLRNG -> P a -> Double -> IO ())
-> FunPtr (P a -> P a -> IO ())
-> FunPtr (P a -> IO (P a))
-> FunPtr (P a -> IO ())
-> FunPtr (P a -> IO ())
-> IO CInt
siman (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
seed)
      Ptr SimulatedAnnealingParams
paramptr P a
confptr
      FunPtr (P a -> Double)
ewrap FunPtr (P a -> P a -> Double)
mwrap FunPtr (GSLRNG -> P a -> Double -> IO ())
stepwrap FunPtr (P a -> P a -> IO ())
cpwrap FunPtr (P a -> IO (P a))
ccwrap FunPtr (P a -> IO ())
dwrap FunPtr (P a -> IO ())
pwrap IO CInt -> (IO CInt -> IO ()) -> IO ()
forall x y. x -> (x -> y) -> y
// String -> IO CInt -> IO ()
check "siman"
    a
result <- P a -> IO a
forall a. P a -> IO a
deRefRead P a
confptr
    P a -> IO ()
forall a. StablePtr a -> IO ()
freeStablePtr P a
confptr
    a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return a
result