/* ga.c
This is search for minimum energy spin glass configurations,
using a genetic algorithm.
The output is solely text, since I wanted to concentrate on the
algorithm. It does, however, print out the lowest energy and the
mean and variance of the energy distribution among the population
for each step along the way.
usage: ga []
J. Watlington, 11/26/95
*/
#include
#include
#include
#define REPORT_PERIOD 0x01 /* num of time steps between prints */
#define POPULATION_SIZE 10
#define SPIN_GLASS_SIZE 100 /* number of spin glass lattice points */
#define EQUAL_PROB 0.5
#define NUM_BUFS 2 /* a double buffer is required */
/* The following variables represent the system being modeled. */
double glass_spin[ NUM_BUFS ][ POPULATION_SIZE ][ SPIN_GLASS_SIZE ];
double glass_J[ SPIN_GLASS_SIZE ]; /* The external stresses */
long rseed = -1; /* random generator state */
/* Quitting conditions */
double epsilon = 1.0;
#define MAX_STALE_GENERATIONS 12
#define MAX_ITERATIONS 4000
/* The following NR routines are linked in. */
extern float gasdev( long *idnum );
extern float ran1( long *idnum );
/* The following routines are declared locally */
int draw_parent( double *repro_prob, int exclude );
void reproduce( int buf, int adam, int eve, int child, double mutation );
double calc_total_energy( double *spin );
void display_glass( int time, int buf, int pop, int show_J );
void main( int argc, char **argv )
{
int i, member;
int adam, eve;
int best_member;
int time = 0;
int cur_gen = 0;
int next_gen = 1;
int no_progress_ctr = 0;
double total_prob;
double total_E, best_E, mean_E, var_E;
double minimum_E = 0.0;
double prev_best_E = 0.0;
double beta = 0.0;
double mutate = 0.0;
double stored_E[ POPULATION_SIZE ]; /* needed for variance calcs */
double repro_prob[ POPULATION_SIZE ]; /* probability of reproducing */
/* Parse the one or two arguments provided by the user. */
if( argc == 3 )
{
sscanf( argv[2], "%lf", &mutate );
}
else if( argc != 2 )
{
printf( "usage: ga \n" );
exit( -1 );
}
sscanf( argv[1], "%lf", &beta );
/* Initialize the spin glass */
for( i = 0; i < SPIN_GLASS_SIZE; i++ )
{
glass_J[ i ] = gasdev( &rseed );
minimum_E -= fabs( glass_J[ i ] );
}
for( member = 0; member < POPULATION_SIZE; member++ )
for( i = 0; i < SPIN_GLASS_SIZE; i++ )
if( ran1( &rseed ) > EQUAL_PROB )
glass_spin[ cur_gen ][ member ][ i ] = 1.0;
else
glass_spin[ cur_gen ][ member ][ i ] = -1.0;
/* Tell the user about the starting conditions */
printf( "Looking for minimum energy %lf (to within %lf)\n",
minimum_E, epsilon );
/* Now start evolving */
do {
/* Generate the probabilities of each population member
reproducing by calculating their energies. */
total_prob = 0.0;
best_E = 1000.0;
mean_E = 0;
best_member = 0;
for( member = 0; member < POPULATION_SIZE; member++ )
{
total_E = calc_total_energy( &glass_spin[cur_gen][member][0] );
repro_prob[ member ] = exp( -beta * ( total_E - minimum_E ) );
total_prob += repro_prob[ member ]; /* acc. the total prob. */
if( total_E < best_E ) /* and find member with least energy */
{
best_E = total_E;
best_member = member;
}
stored_E[ member ] = total_E; /* save for variance calcs */
mean_E += total_E; /* acc the total energy */
}
/* Calculate some statistics about the population energy levels,
since just looking at the raw numbers is useless ! */
mean_E = mean_E / (double)POPULATION_SIZE;
var_E = 0.0;
for( member = 0; member < POPULATION_SIZE; member++ )
{
total_E = stored_E[ member ] - mean_E;
var_E += total_E * total_E;
}
var_E = var_E / (double)POPULATION_SIZE;
/* Now build the probability array used by draw_parent() */
repro_prob[ 0 ] = repro_prob[ 0 ] / total_prob;
for( member = 1; member < POPULATION_SIZE; member++ )
{
repro_prob[ member ] = (repro_prob[ member ] / total_prob)
+ repro_prob[ member-1 ];
}
repro_prob[ POPULATION_SIZE - 1 ] = 1.0;
/* Tell the user about life... */
if( !(time & REPORT_PERIOD) )
printf( "t: %4d best: %3.1lf mean: %3.1lf var: %4.1lf\n", time,
best_E, mean_E, var_E );
/* For each member of the next generation, draw two parents from
the current generation and combine them into the offspring */
for( member = 0; member < POPULATION_SIZE; member++ )
{
adam = draw_parent( repro_prob, -1 );
eve = draw_parent( repro_prob, adam );
reproduce( cur_gen, adam, eve, member, mutate );
}
/* Check for forward progress ! */
if( prev_best_E == best_E )
no_progress_ctr++;
else
{
no_progress_ctr = 0;
prev_best_E = best_E;
}
/* Increment generations */
cur_gen = next_gen;
next_gen = 1 - next_gen;
} while( ((best_E - minimum_E) > epsilon) &&
(no_progress_ctr < MAX_STALE_GENERATIONS) &&
(time++ < MAX_ITERATIONS));
if( no_progress_ctr >= MAX_STALE_GENERATIONS )
printf( "Glass frozen in non-optimal state !\n" );
printf( "Best energy at end: %lf at time %d\n", best_E, time );
display_glass( time, next_gen, best_member, 1 );
}
/* draw_parent
Given a table which contains the cumulative probabilities
summed, (ie. a, a+b, a+b+c, a+b+c+d, ...), this function
randomly selects a parent. The exclude argument allows a
particular parent to be excluded from the drawing.
*/
int draw_parent( double *repro_prob, int exclude )
{
float uniform_dev;
int lower_limit, upper_limit;
int target;
do {
uniform_dev = ran1( &rseed );
lower_limit = 0;
upper_limit = POPULATION_SIZE - 1;
while( lower_limit != upper_limit )
{
target = lower_limit + ((upper_limit - lower_limit) >> 1);
if( uniform_dev > repro_prob[ target ] )
lower_limit = target + 1;
else
upper_limit = target;
}
} while( lower_limit == exclude );
return( lower_limit );
}
/* reproduce
This function takes two parents and creates a child from
their genetic information by randomly selecting a switchover
point, and copying information up to the switchover from
one parent and after the switchover from the other.
The mutate argument allows a mutation rate for the copying to
be specified.
*/
void reproduce( int buf, int adam, int eve, int child, double mutate )
{
int switchpoint;
int x;
double *parent;
double *baby;
parent = &glass_spin[ buf ][ adam ][ 0 ];
baby = &glass_spin[ 1 - buf ][ child ][ 0 ];
switchpoint = (int)( ran1( &rseed ) * (float)(SPIN_GLASS_SIZE - 1) );
for( x = 0; x < switchpoint; x++ )
*baby++ = *parent++;
parent = &glass_spin[ buf ][ eve ][ switchpoint ];
for( x = switchpoint; x < SPIN_GLASS_SIZE; x++ )
*baby++ = *parent++;
if( mutate > 0.0 )
{
baby = &glass_spin[ 1 - buf ][ child ][ 0 ];
for( x = 0; x < SPIN_GLASS_SIZE; x++ )
if( ran1( &rseed ) <= mutate )
baby[x] = -baby[x];
}
}
/* calc_total_energy
This function takes a pointer to a glass spin orientation
matrix, and returns the energy associated with that
configuration. The same external stress matrix (J) is used
for all glass spin energy evaluations.
*/
double calc_total_energy( double *spin )
{
int x;
double total = 0;
for( x = 0; x < SPIN_GLASS_SIZE-1; x++ )
{
total += (double)(glass_J[x] * spin[x] * spin[ x+1 ]);
}
/* Handle periodic boundary condition */
x = SPIN_GLASS_SIZE - 1;
total += (double)(glass_J[x] * spin[x] * spin[ 1 ]);
return( -total );
}
#define NUM_COLS 25
void display_glass( int time, int buf, int pop, int show_J )
{
int index = 0;
int row;
int col;
static int num_rows = SPIN_GLASS_SIZE / NUM_COLS;
printf( "At time %d:", time );
if( show_J )
{
printf( "\n**** J" );
for( row = 0; row < num_rows; row++ )
{
printf( "\n%2d>", row );
for( col = 0; col < NUM_COLS; col++ )
{
printf( " %2.1f", glass_J[ index++ ] );
}
}
printf( "\n**** S" );
index = 0;
}
for( row = 0; row < num_rows; row++ )
{
printf( "\n%2d>", row );
for( col = 0; col < NUM_COLS; col++ )
{
printf( " %2.0f", glass_spin[ buf ][ pop ][ index++ ] );
}
}
printf( "\n" );
}