8
\$\begingroup\$

This simulated annealing program tries to look for the status that minimizes the energy value calculated by the energy function.

The status class, energy function and next function may be resource-intensive on future usage, so I would like to know if this is a suitable way to code it.

#pragma once

#include <cmath>
#include <random>
#include <utility>
#include <algorithm>

// To find a status with lower energy according to the given condition
template<typename status, typename count, typename energy_function, typename temperature_function, typename next_function, typename generator>
status simulated_annealing(status i_old, count c, const energy_function& ef, const temperature_function& tf, const next_function& nf, generator& g){

    auto   e_old  = ef(i_old);

    status i_best = i_old;
    auto   e_best = e_old;

    std::uniform_real_distribution<decltype(e_old)> rf(0, 1);

    for(; c > 0; --c){
        status i_new = nf(i_old, g);
        auto   e_new = ef(i_new);

        if(e_new < e_best){
            i_best =           i_new ;
            e_best =           e_new ;
        }

        if( e_new < e_old || std::exp( (e_old - e_new) / tf(c) ) > rf(g) ){
            i_old  = std::move(i_new);
            e_old  = std::move(e_new);
        }
    }
    return(i_best);
}

Example usage:

#include <cmath>
#include <random>

#include <chrono>
#include <cstdlib>
#include <iostream>

#include "simulated_annealing.hpp"

template<typename Clock = std::chrono::system_clock>
class timer{
private:
    typedef std::chrono::time_point<Clock> time_point;
    time_point _start;
    time_point _finish;
public:
    void start(){
        _start = Clock::now();
    }
    void finish(){
        _finish = Clock::now();
    }
    template<typename rtype>
    rtype count(){
        std::chrono::duration<rtype> elapsed = _finish - _start;
        return elapsed.count();
    }
};


auto ef = [](long double x){
    return std::abs( (x - 10)*(x + 20)*(x - 30) );
};

auto tf = [](long double k){
    return std::exp( -20 * k );
};

template<typename generator>
long double nf(long double x, generator& g){
    std::normal_distribution<decltype(x)> d(0, 1);
    return x + d(g);
}

int main(int argc, char const *argv[])
{
    std::random_device rd;
    std::mt19937_64 g(rd());
    timer<std::chrono::high_resolution_clock> t;

    long double root   = std::atof(argv[1]);
    std::size_t count  = std::atoi(argv[2]);

    t.start();
    root = simulated_annealing(root, count, ef, tf, nf<decltype(g)>, g);
    t.finish();

    std::cout << "Result: " << root << std::endl;
    std::cout << "Time elapaed: " << t.count<double>() << '\n';

    return 0;
}
\$\endgroup\$
  • \$\begingroup\$ One point (not worth an answer): Your counter c starts high and goes to 0. Your temperature function tf shrinks with k. Basically, you are raising the temperature with time instead of lowering it. \$\endgroup\$ – Avi Ginsburg Jan 22 '15 at 14:07
  • \$\begingroup\$ Was davidhigh's answer satisfying to you ? If yes please consider marking it as accepted. If anyone incorporated the suggestions into the code, it would be nice to share the resulting code (I would, but am considering Simulated Annealing Information, Taygeta Scientific Inc. first.) \$\endgroup\$ – Stéphane Gourichon May 2 '15 at 11:27
6
\$\begingroup\$

I think the algorithm is well abstracted and the code is nicely written, so there is not all too much you can do in order to make it more efficient. Still I have some minor suggestions for improvement:


First, regarding the algorithm itself:

  1. When the check if(e_new < e_best) is successful, you can directly treat the step as accepted. Although this might be negligible, you avoid one more check.

  2. In your check whether to accept the new step, you can think about choosing a cutoff (say ~ 10.0) parameter to bypass the relatively expensive calculation of std::exp as well as the random number generation.

Summarizing 1. and 2., one gets a slightly more verbose but hopefully more efficient code:

    if(e_new < e_best){
        i_best =           i_new ;
        e_best =           e_new ;
        i_old  = std::move(i_new);
        e_old  = std::move(e_new);
        continue;
    }

    auto t = tf(c);
    auto delta_e = e_new - e_old;

    if( delta_e > 10.0 * t) continue;  //as std::exp(-10.0) is a very small number

    if( delta_e <= 0.0 || std::exp( - delta_e / t ) > rf(g) ){
        i_old  = std::move(i_new);
        e_old  = std::move(e_new);
    }

Next, regarding your choice of functions: Here I would try to avoid expensive functions as well. I know, this might often not be possible as the functions are determined by your problem at hand. However, for example, in your temperature function you can think about replacing the std::exp(-20*k) by something more efficient such as 1/(k*k).


Regarding the language:

  1. You can get a bit more of a C++11 touch if you use universal references instead of lvalue references.

    template<typename status, typename count, typename energy_function
           , typename temperature_function, typename next_function, typename generator>
    status simulated_annealing(status&& i_old, count&& c, energy_function&& ef
                             , temperature_function&& tf, next_function&& nf, generator&& g)
    {
        //...
    }
    

    This is maybe also a bit more consistent with your std::move's of your status and count-type variables (which, as a side node, will resolve to simple copies when fundamental types such as double or size_t are used).

    Further, you can think about replacing your type count by size_t, since that is what this type is for. But this just because I personally try to avoid template parameters when possible.


Regarding the design:

  1. The nf function can be made more generic by defining a static generator inside

    long double nf(long double x){
        static std::mt19937_64 g(std::random_device());
        std::normal_distribution<decltype(x)> d(0, 1);
        return x + d(g);
    }
    

    I wouldn't pass a generator here, as it doesn't fit well into a generic desing. (For instance, the next state could also be found deterministically).

  2. The same holds for your main simulated_annealing function. I wouldn't bother the user the user to define a random number generator, but just define it inside -- moreover since it is used "only" to create random numbers in [0,1] (once you've followed suggestion 4.).

    template<typename status, typename count, typename energy_function
           , typename temperature_function, typename next_function>
    status simulated_annealing(status&& i_old, count&& c, energy_function&& ef
                         , temperature_function&& tf, next_function&& nf)
    {
        std::mt19937_64 g(std::random_device());            
        //...
    }
    

    This you can simply call the routine via

    root = simulated_annealing(root, count, ef, tf, nf);
    

    which is more concise and better reflects the main components of a simulated annealing algorithm.

\$\endgroup\$
  • \$\begingroup\$ I would create interfaces which allow the specification of the generator and allow use of a default. There are valid use cases for both approaches. \$\endgroup\$ – Brandon Kohn Dec 1 '16 at 1:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.