C++11 class interface for GSL simulated annealing

siman.cpp

#include <functional>
#include <memory>
#include <gsl/gsl_siman.h>

namespace gsl {

template < typename State >
class Siman {

    static double Ef(void *xp) {
        auto self = static_cast<Siman*>(xp);
        return self->m_Ef(self->m_x);
    }

    static void take_step(gsl_rng const *, void *xp, double step_size) {
        auto self = static_cast<Siman*>(xp);

        self->m_x = self->m_take_step(self->m_x, step_size);
    }

    static double distance(void *xp, void *yp) {
        auto oldconf = static_cast<Siman*>(xp);
        auto newconf = static_cast<Siman*>(yp);

        return oldconf->m_distance(oldconf->m_x, newconf->m_x);
    }

    static void copyfunc(void *source, void *dest) {
        auto src = static_cast<Siman*>(source);
        auto dst = static_cast<Siman*>(dest);

        *dst = *src;
    }

    static void * copy_constructor(void *xp) {
        auto self = static_cast<Siman*>(xp);

        return static_cast<void*>(new Siman(*self));
    }

    static void destructor(void *xp) {
        auto self = static_cast<Siman*>(xp);

        delete self;

        self = nullptr;
    }
    
    std::function<double(State)> m_Ef;
    
    std::function<State(State, double)> m_take_step;
    
    std::function<double(State, State)> m_distance;
    
    State m_x;

    std::shared_ptr<gsl_rng> m_r;

    gsl_siman_params_t m_params;

public:
    
    template < typename EnergyFunctor, typename StepFunctor, typename DistanceFunctor, typename ... Args >
    Siman(State x, EnergyFunctor Ef, StepFunctor take_step, DistanceFunctor distance, Args ... params)
        : m_Ef{Ef}
        , m_take_step{take_step}
        , m_distance{distance}
        , m_x{x}
        , m_r{gsl_rng_alloc(gsl_rng_env_setup()),gsl_rng_free}
        , m_params{params...}
    {}

    State solve() {
        gsl_siman_solve(m_r.get(),
                        static_cast<void*>(this),
                        Ef,
                        take_step,
                        distance,
                        nullptr, // print_position
                        copyfunc,
                        copy_constructor,
                        destructor,
                        0, // element_size
                        m_params);
        return m_x;
    }
};

template < typename State, typename ... Args >
Siman<State> siman(State x, Args && ... args) {
    return Siman<State>{std::forward<State>(x), std::forward<Args>(args)...};
}

} // namespace gsl

// =====================================================
// =====================================================
// =====================================================
// =====================================================

#include <cmath>
#include <iostream>
#include <random>

int main()
{
    // set up parameters for this simulated annealing run
    int n_tries = 200;        // how many points do we try before stepping
    int iters_fixed_t = 1000; // how many iterations for each T?
    double step_size = 1.0;   // max step size in random walk
    double k = 1.0;           // Boltzmann constant
    double t_initial = 0.008; // initial temperature
    double mu_t = 1.003;      // damping factor for temperature
    double t_min = 2.0e-6;

    double x_initial = 15.5;

    auto energy = [] (double x) {
        return std::exp(-std::pow(x-1.0,2))*std::sin(8*x);
    };
    auto step = [] (double x, double step_size) {
        static thread_local std::default_random_engine eng{std::random_device{}()};
        double u = std::uniform_real_distribution<double>{0,1}(eng);
        return u * 2 * step_size - step_size + x;
    };
    auto distance = [] (double x, double y) {
        return std::abs(x - y);
    };

    auto annealer = gsl::siman(x_initial, energy, step, distance,
                               n_tries, iters_fixed_t, step_size,
                               k, t_initial, mu_t, t_min);
    auto min = annealer.solve();

    std::cout << "Minimum at " << min << '\n';
}