Simulation annealing implemented in python.

sa.py

#!/usr/bin/env python
"""
Simulated annealing module.

Usage:
  sa.py [options]

Options:
  -h, --help  Show this message and exit.
"""

from __future__ import print_function
from docopt import docopt
import numpy as np
from random import random

def rosenbrock(x,y):
    return (1.0-x)*(1.0-x) +100.0*(y-x*x)*(y-x*x)

def run(func,varini,num_iteration,tempini,
        ptarget=0.5,update_interval=10,amin=0.1,amax=10.0,
        decay='liner',decay_coeff=5.0):
    """
    Run simulated annealing.

    func:
        Function to be minimized.
    varini:
        Variables to be optimized.
    num_iteration:
        Maximum number of iteration.
    tempini:
        Initial temperature. The temperature decreases to 0 in *num_iteration*.
        Degree of the temperature must be estimated heuristically by some
        preliminary calculations or your experience.
    ptarget:
        Target of probability to be achieved in each variable.
    update_interval:
        Update interval of max width of deviation to achieve probability target..
    amin, amax:
        Min and Max of coefficient update the max width of deviation.
    decay:
        Type of temperature decay: liner and exp.
    decay_coeff:
        Coefficient for temperature decay. Only used in exp decay.
    """

    ndim= len(varini)
    var= varini.copy()
    fini= func(varini[0],varini[1])
    fprev= fini
    temp= tempini

    #...set initial maximum deviation width
    width= np.zeros(ndim,dtype=float)
    history= []
    upcount= []
    for i in range(ndim):
        width[i]= abs(varini[i]) *0.1
        history.append([])
        upcount.append(update_interval)

    print('{0:6d} {1:12.5f} {2:8.4f} {3:8.4}'.format(0,fini,
                                                     var[0],var[1]),end='')
    print(' {0:6.3f} {1:6.3f}'.format(width[0],width[1]))
    for iteration in range(num_iteration):
        
        #...decide which var is changed
        idev= int(ndim *random())

        #...decide deviation width
        dev= width[idev] *(random()-0.5)

        #...add deviation
        var[idev]= var[idev] +dev

        #...evaluate function
        fnow= func(var[0],var[1])
        
        #...judge whether or not adopt the deviation
        prob= min(np.exp(-(fnow-fprev)/temp),1.0)
        if prob > random():
            fprev= fnow
            print('{0:6d} {1:12.5f} {2:8.4f} {3:8.4f}'.format(iteration+1,fnow,
                                                               var[0],var[1]),end='')
            print(' {0:6.3f} {1:6.3f}'.format(width[0],width[1]))
            history[idev].append(1)
        else:
            var[idev]= var[idev] -dev   # restore var
            history[idev].append(0)
        if len(history[idev]) > update_interval:
            history[idev].pop(0)

        #...adjust maximum deviation width
        upcount[idev] -= 1
        if upcount[idev] == 0:
            p= float(sum(history[idev]))/update_interval
            if p == 1.0:
                alpha= amax
            elif p == 0.0:
                alpha= amin
            else:
                alpha= np.sqrt(np.log(ptarget)/np.log(p))
            alpha= max(alpha,amin)
            width[idev]= alpha *width[idev]
            upcount[idev]= update_interval

        #...decrease temperature
        if decay == 'exp':
            temp= tempini*np.exp(-decay_coeff
                                 *iteration/num_iteration)
        else:
            temp= -float(iteration)/num_iteration*tempini +tempini


if __name__ == '__main__':

    args= docopt(__doc__)
    
    varini= np.array([-3.,-4.])
    run(rosenbrock,varini,10000,10.0,decay='exp')
    

Simulated annealing implemented with Python.