/* a collection of C routines for general purpose Simulated Annealing   intended to be the C equivalent of the C++ Simulated Annealing object   SimAnneal   Uses Cauchy training*/static char rcsid[] = "@(#)sa.c	1.2 15:54:31 3/30/93   EFC";#include <stdio.h>#include <math.h>#include <float.h>/* #include <stdlib.h>		/* for malloc *//* #include <rands.h> */#include <r250.h>#include "sa.h"/* #define DEBUG */#ifdef _R250_H_#define uniform(a,b)    ( a + (b - a) * dr250() )#endif#ifndef HUGE#define HUGE	HUGE_VAL#endif#ifndef PI#define PI		3.1415626536#endif#ifndef PI2#define PI2		(PI/2.0)#endiftypedef struct{	CostFunction func;	int dimension, maxit, dwell;	float *x, *xnew, *xbest;	float dt, c_jump, K, rho, t0, tscale, range;	double y, ybest;} SimAnneal;static SimAnneal s;/* iterate a few times at the present temperature *//* to get to thermal equilibrium */#ifdef NO_PROTOstatic int equilibrate(t, n)float t;int n;#elsestatic int equilibrate(float t,int n)#endif{	int i, j, equil = 0;	float xc, ynew, c, delta, p;	float *xtmp;	delta = 1.0;	for (i = 0; i < n; i++)	{		for (j = 0; j < s.dimension; j++)		{			xc = s.rho * t * tan ( uniform( -s.range, s.range ) ); 			s.xnew[j] = s.x[j] + xc;		}		/* "energy" */		ynew = s.func( s.xnew );		c = ynew - s.y;				if (c < 0.0)		/* keep xnew if energy is reduced */		{			xtmp = s.x;			s.x = s.xnew;			s.xnew = xtmp;			s.y = ynew;			if ( s.y < s.ybest )			{				for (j = 0; j < s.dimension; j++)					s.xbest[j] = s.x[j];				s.ybest = s.y;			}			delta = fabs( c );			delta = ( s.y != 0.0 ) ? delta / s.y : ( ynew != 0.0 ) ?					delta / ynew : delta;			/* equilibrium is defined as a 10% or smaller change			   in 10 iterations */			if ( delta < 0.10 )				equil++;			else				equil = 0;		}		else		{		/* keep xnew with probability, p, if ynew is increased *//*			p = exp( - (ynew - s.y) / (s.K * t) );			if ( p > number(0.0, 1.0) )			{				xtmp = s.x;				s.x = s.xnew;				s.xnew = xtmp;				s.y = ynew;				equil = 0;			}			else*/				equil++;		}				if (equil > 9)			break;	}	return i + 1;}/* initialize internal variables and define the cost function */#ifdef NO_PROTOint SAInit(f, d)CostFunction f;int d;#elseint SAInit(CostFunction f, int d)#endif{	int space;        	s.func = f;        s.dimension = d;        s.t0 = 0.0;        s.K  = 1.0;        s.rho = 0.5;        s.dt  = 0.1;        s.tscale = 0.1;        s.maxit = 400;        s.c_jump = 100.0;	s.range = PI2;	s.dwell = 20;	space = s.dimension * sizeof(float);                if ( (s.x = (float *)malloc( space ))  == NULL )        		return 0;        if ( (s.xnew = (float *)malloc( space )) == NULL )        		return 0;        if ( (s.xbest = (float *)malloc( space )) == NULL )        		return 0;        s.y = s.ybest = HUGE;#ifdef _R250_H_	r250_init( 12331 );#endif   	return 1;        }void SAfree(){	free( s.x );        free( s.xnew );        free( s.xbest );        s.dimension = 0;}#ifdef NO_PROTOint SAiterations(m)int m;#elseint SAiterations(int m)#endif{	if ( m > 0 )        	s.maxit = m;       return s.maxit;}#ifdef NO_PROTOint SAdwell(m)int m;#elseint SAdwell(int m)#endif{        if ( m > 0 )                s.dwell = m;       return s.dwell;}#ifdef NO_PROTOfloat SABoltzmann(k)float k;#elsefloat SABoltzmann(float k)#endif{	if ( k > 0.0 )        	s.K = k;       return s.K;}#ifdef NO_PROTOfloat SAtemperature(t)float t;#elsefloat SAtemperature(float t)#endif{	if ( t > 0.0 )        	s.t0 = t;       return s.t0;}#ifdef NO_PROTOfloat SAlearning_rate(r)float r;#elsefloat SAlearning_rate(float r)#endif{	if ( r > 0.0 )        	s.rho = r;       return s.rho;}#ifdef NO_PROTOfloat SAjump(j)float j;#elsefloat SAjump(float j)#endif{	if ( j > 0.0 )        	s.c_jump = j;       return s.c_jump;}#ifdef NO_PROTOfloat SArange(r)float r;#elsefloat SArange(float r)#endif{        if ( r > 0.0 )                s.range = r;       return s.range;}#ifdef NO_PROTOvoid SAinitial(xi)float* xi;#elsevoid SAinitial(float* xi)#endif{	int k;	for (k = 0; k < s.dimension; k++)		s.x[k] = xi[k];}#ifdef NO_PROTOvoid SAcurrent(xc)float* xc;#elsevoid SAcurrent(float* xc)#endif{	int k;        	for (k = 0; k < s.dimension; k++)		xc[k] = s.x[k];}#ifdef NO_PROTOvoid SAoptimum(xb)float* xb;#elsevoid SAoptimum(float* xb)#endif{	int k;        	for (k = 0; k < s.dimension; k++)		xb[k] = s.xbest[k];}/* increase the temperature until the system "melts" */#ifdef NO_PROTOfloat SAmelt(iters)int iters;#elsefloat SAmelt(int iters)	#endif{	int i, j, ok = 0;	float xc, ynew, t, cold, c = 0.0;	int n = iters;	if ( n < 1 )		n = s.maxit;	t = s.t0;	for (i = 0; i < n; i++)	{		if (i > 0 && c > 0.0)		{			cold = c;			ok = 1;		}				t += s.dt;		for (j = 0; j < s.dimension; j++)		{			xc = s.rho * t * tan ( uniform( -s.range, s.range ) ); 			s.x[j] += xc;		}		equilibrate( t, s.dwell);				/* "energy" */		ynew = s.func( s.x );		c = ynew - s.y;		if ( c < 0.0 && ynew < s.ybest)		{			for (j = 0; j < s.dimension; j++)				s.xbest[j] = s.x[j];			s.ybest = ynew;		}		s.y = ynew;		if ( ok && c > (s.c_jump * cold) )	/* phase transition */			break;	}	return s.t0 = t;	}/* cool the system with annealing */#ifdef NO_PROTOfloat SAanneal(iters)#elsefloat SAanneal(int iters)#endif{	int i, j;	float xc, p, ynew, t, c, dt, told;	float *xtmp;	int n = iters;	if ( n < 1 )		n = s.maxit;	equilibrate( s.t0, 10 * s.dwell );	told = s.t0;	for (i = 0; i < n; i++)	{		t = s.t0 /(1.0 + i * s.tscale);		dt = t - told;		told = t;		equilibrate(t, s.dwell);		for (j = 0; j < s.dimension; j++)	        {			xc = s.rho * t * tan ( uniform( -s.range, s.range ) ); 			s.xnew[j] = s.x[j] + xc;		}		/* "energy" */		ynew = s.func( s.xnew );		c = ynew - s.y;				if (ynew <= s.y)	/* keep xnew if energy is reduced */		{			xtmp = s.x;			s.x = s.xnew;			s.xnew = xtmp;			s.y = ynew;			if ( s.y < s.ybest )			{				for (j = 0; j < s.dimension; j++)					s.xbest[j] = s.x[j];				s.ybest = s.y;			}			continue;		}		else		{		/* keep xnew with probability, p, if ynew is increased */			p = exp( - (ynew - s.y) / (s.K * t) );			if ( p > uniform(0.0, 1.0) )			{				xtmp = s.x;				s.x = s.xnew;				s.xnew = xtmp;				s.y = ynew;			}		}	}	equilibrate( t, 10 * s.dwell );		return s.t0 = t;}