From: ajayshah@almaak.usc.edu (Ajay Shah)
Newsgroups: comp.sources.misc
Subject: v22i062:  gademo - Genetic Algorithm for unconstrained nonlinear maximisation, Part01/01
Date: 25 Aug 91 22:44:15 GMT

Submitted-by: Ajay Shah <ajayshah@almaak.usc.edu>
Posting-number: Volume 22, Issue 62
Archive-name: gademo/part01

I wrote this program to help me understand a simple GA.  It is 
not optimised for speed.

Out of the box, it will find the _global_ maximum of a univariate function.
That function need not be continuous.  Both these aspects:

   - global maximum as opposed to local, and 
   - ability to deal with non-differentiable and discontinuous functions

are plus points when compared with gradient-oriented algorithms.

The code currently contains two functions as demos, one of which
has a unique local (global) maximum, the other has one local 
maximum and a different global maximum.

Ajay Shah
-------
#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of shell archive."
# Contents:  README gademo.c makefile mylib.c mylib.h
# Wrapped by ajayshah@max on Fri Aug 23 15:37:34 1991
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'README' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'README'\"
else
echo shar: Extracting \"'README'\" \(1961 characters\)
sed "s/^X//" >'README' <<'END_OF_FILE'
X
XI wrote this program to help me understand a simple GA.  It is 
Xnot optimised for speed.
X
XOut of the box, it will find the _global_ maximum of a univariate function.
XThat function need not be continuous.  Both these aspects:
X
X	- global maximum as opposed to local, and 
X	- ability to deal with non-differentiable and 
X		discontinuous functions
X
Xare plus points when compared with gradient-oriented algorithms.
X
XThe code currently contains two functions as demos, one of which
Xhas a unique local (global) maximum, the other has one local 
Xmaximum and one global maximum.
X
XWhat you do:
X
X	1. in makefile, the CFLAGS contains cpp defines which guide
X	the code.  
X
X	You say -DVERBOSE to get a detailed trace of the programs
X	activities, and skip it (do not put -DVERBOSE) to get a more
X	terse summary.
X
X	You say -DTRIVIAL to maximise the trivial function and
X	-DHAIRY to maximise the multiple-local-maxima function.
X
X	2. Once makefile is setup to suit your tastes, say "make".
X	You will get an executable called "main".  Run it.  It's
X	interesting working through different combinations of these two
X	cpp variables.
X
X	3. Now you can get to playing with tuning parameters
X	embedded inside the code.  Early in the file "gademo.c",
X	you will find these hashdefines: 
X
X	SIZE_POPULATION
X	COSMIC_RAY_DENSITY
X	TIME0_MIN
X	TIME0_MAX
X	CONVERGED
X
X	I won't document them here; there are comments alongside
X	which explain each of them.  You can modify each of these
X	and see the effects.  Perhaps the most interesting is the
X	COSMIC_RAY_DENSITY -- it influences the speed to
X	convergence in dramatic ways.
X
X	4. Finally, right after these hashdefines, is the code for
X	the function being maximised.  It's just 
X
X		double f(double x)
X
X	While you must stay with this prototype, you can replace
X	the innards with anything you want (i.e., any f:R^1-->R^1).
X	Be sure to choose TIME0_MIN and TIME0_MAX appropriately.
X
X	5. Deeper in gademo.c, after the code for f(x), is the GA
X	code.
X
END_OF_FILE
if test 1961 -ne `wc -c <'README'`; then
    echo shar: \"'README'\" unpacked with wrong size!
fi
# end of 'README'
fi
if test -f 'gademo.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'gademo.c'\"
else
echo shar: Extracting \"'gademo.c'\" \(5869 characters\)
sed "s/^X//" >'gademo.c' <<'END_OF_FILE'
X
X#include <stdio.h>
X#include <math.h>
X#include <malloc.h>
X#include "mylib.h"
X
X/* This is a simple demo program where we use a genetic algorithm
Xto solve a nonlinear unconstrained maximisation problem. 
X
XThe program takes no arguments.  All tuning parameters are in
Xhash-defines at the top of file.  */
X
X/*------------------------- Start of tuning parameters */
X#define SIZE_POPULATION 10
X	/* each generation normally (not always) contains so many animals */
X#define COSMIC_RAY_DENSITY 0.4
X	/* It's about the amount of mutation in going from 
X	generation n to (n+1). */
X#define TIME0_MIN -20.0
X#define TIME0_MAX 20.0
X	/* first generation animals will be uniformly distributed
X	between these limits.  Be sure to make 'em generous! */
X#define CONVERGED 0.001
X	/* convergence criterion. */
X/* ------------------------- End of tuning parameters */
X
X/* Objective function */
X#ifdef TRIVIAL
Xdouble f(double x)
X{
X	return (-x*x - 2.0*x - 1.0);
X}	/* a trivial demo; maximised at -1 */
X#endif TRIVIAL
X
X#ifdef HAIRY
Xdouble f(double x)
X{
X	double x2, x3, x4;
X
X	x2 = x*x;
X	x3 = x2*x;
X	x4 = x3*x;
X	return ((47*x/3) - (239*x2/12) + (25*x3/3) - (13*x4/12));
X}
X	/* This function has two local maxima and one local minimum.
X	f(0.6) = 3.8896
X	f(1.85) = 0.8924
X	f(3.35) = 5.8234	global maximum.
X	*/  
X#endif HAIRY
X
X/* Algorithm from here on works for any scalar function f(x). */
X
X/* Data structures: at time t, the population will be
Xcharacterised by three things:
X	N, the number of animals
X	x, a vector with N elements where x_i is the attribute of animal i
X	y, a vector where y_i = f(x_i)
X*/
X
X#define uniform(min, max, seed) (min + ((max - min)*ran1(seed)))
X
Xdouble *first_gen()
X/* This creates a vector of animals in the first generation. */
X{
X	double *t;
X	int i;
X	int seed = 0;
X
X	t = (double *) malloc(SIZE_POPULATION*sizeof(double));
X	for (i=0; i<SIZE_POPULATION; i++) 
X		t[i] = uniform(TIME0_MIN, TIME0_MAX, &seed);
X	return t;
X}
X
Xvoid ofn2probs(int N, double **y, int *locbest, double *bestsofar)
X/* Given a vector of values of the objective function, we
Xcreate probabilities of sex, in-place */
X{
X	double avg, min, max, scale, sum, tmp;
X	int i;
X
X	min = 1.0e30; max = -1.0e30;
X	for (i=0; i<N; i++) {
X		tmp = (*y)[i];
X		if (tmp < min) min = tmp;
X		if (tmp > max) {
X			max = tmp;
X			*locbest = i;
X		}
X	}
X	*bestsofar = max;	
X	/* in this y vector, this is the best value of the objective
X	function found */
X
X	/* now replace by fitness as follows: let avg=(min+max)/2.
X		<= avg --> 0
X		max-->1.0
X	*/
X	avg = (max + min)/2.0; scale = max - avg; sum = 0.0;
X	for (i=0; i<N; i++) {
X		tmp = (*y)[i];
X		if (tmp <= avg) tmp = 0.0;
X		else tmp = (tmp - min)/scale;
X		sum += tmp;
X		(*y)[i] = tmp;
X	}
X
X	/* these aren't marginal probabilities yet.  They don't add
X	up to one (they addup to sum). */
X	for (i=0; i<N; i++) (*y)[i] = (*y)[i]/sum;
X
X	/* Now we have the vector of marginal probabilities
X	mapping every animal into probability of having sex */
X}
X/* The ofn2probs operation is three passes over the table of 
Xsize N.  It does 2N double comparisons takes 4N (double) flops. */
X
Xvoid print_generation(int time, int N, double *x, double *y, double best)
X{
X	int i;
X
X	printf("\nGeneration %d, best f(x) = %f\n", time, best);
X#ifdef VERBOSE
X	for (i=0; i<N; i++) 
X		printf("Pr(x = %7.3f) = %7.3f\n", x[i], y[i]);
X#else
X	for (i=0; i<N; i++) printf("%5.1f ", x[i]);
X#endif
X}
X
Xvoid sex(double *oldx, double *oldy, double **newx, int *seed)
X/* Given one generation characterised by vectors oldx and oldy,
Xcreate a new generation in vector newx. */
X{
X	double *cdf;
X	int *outcomes, N, i, j, numdarts;
X	double avg, p1, p2, child;
X
X	N = SIZE_POPULATION;
X	cdf = dvector(1, N);
X	numdarts = 2*N;
X	outcomes = ivector(1, numdarts);
X	pdf2cdf(oldy-1, cdf, N);
X	throw_darts(cdf, N, seed, outcomes, numdarts);
X#ifdef VERBOSE
X	printf("\nSex:\t");
X	for (i=1; i<=numdarts; i+=2) 
X		printf("%d+%d ", outcomes[i], outcomes[i+1]);
X	printf("\n\n");
X#endif VERBOSE
X	for (i=1, j=0; j<N; j++) {
X		p1 = oldx[outcomes[i]-1]; i++;
X		p2 = oldx[outcomes[i]-1]; i++;
X		avg = (0.5*(p1 + p2));
X		child = avg + (COSMIC_RAY_DENSITY*avg*gasdev(seed));
X#ifdef VERBOSE
X		printf("parents %f %f, child %f\n", p1, p2, child);
X#endif
X		(*newx)[j] = child;
X	}
X
X	free_dvector(cdf, 1, N);
X	free_ivector(outcomes, 1, numdarts);
X}
X
Xint main()
X{
X	int N, i, done, time;
X	double *x, *y, best;
X	double *oldx, *oldy, *newx, *newy;
X	double oldbest, newbest;
X	double *_dptr;
X	int seed, which_is_best;
X
X	/* setup forest at time 0 */
X	N = SIZE_POPULATION; 
X	x = first_gen();
X	y = (double *) malloc(SIZE_POPULATION*sizeof(double));
X	for (i=0; i<N; i++) y[i] = f(x[i]);
X	ofn2probs(N, &y, &which_is_best, &best); 
X		/* values of objective fn --> probs */
X
X	/* We now have one generation: N, x and y (in probabilities).
X	In the main loop of the GA, we're going to constantly
X	manipulate two generations: an old and a new.  We will
X	share the same storage by shuffling pointers.  No memory
X	allocation will happen while GA is running.  */
X
X	done = FALSE; time = 0; oldx = x; oldy = y; 
X	oldbest = best; seed = 0;
X
X	newx = (double *) malloc(SIZE_POPULATION*sizeof(double));
X	newy = (double *) malloc(SIZE_POPULATION*sizeof(double));
X
X	while (done == FALSE) {
X		print_generation(time, N, oldx, oldy, oldbest);
X
X		do {
X			printf(".");
X			sex(oldx, oldy, &newx, &seed);
X			for (i=0; i<N; i++) newy[i] = f(newx[i]);
X			ofn2probs(N, &newy, &which_is_best, &newbest);
X			time++;
X		} while (newbest <= oldbest);
X
X		if (fabs(oldbest - newbest) < CONVERGED) done = TRUE;
X
X		/* now get ready for next iteration. */
X		_dptr = oldx; /* points to storage used by oldx */
X		oldx = newx; newx = _dptr;
X
X		_dptr = oldy; oldy = newy; newy = _dptr;
X		oldbest = newbest;
X	}
X
X	printf("\n\nConverged!\n");
X	printf("Best x = %10.6f, f(x) here is %10.6f\n\n", 
X		oldx[which_is_best], f(oldx[which_is_best]));
X	print_generation(time, N, oldx, oldy, oldbest);
X	printf("\n");
X
X	return 0;
X}
X
END_OF_FILE
if test 5869 -ne `wc -c <'gademo.c'`; then
    echo shar: \"'gademo.c'\" unpacked with wrong size!
fi
# end of 'gademo.c'
fi
if test -f 'makefile' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'makefile'\"
else
echo shar: Extracting \"'makefile'\" \(199 characters\)
sed "s/^X//" >'makefile' <<'END_OF_FILE'
X
XL = mylib.o
XCC = gcc
XCFLAGS = -DHAIRY
X
Xmain : gademo.c makefile mylib.o
X	$(CC) $(CFLAGS) -o main gademo.c $L -lm
X
Xmylib.o : mylib.c
X	$(CC) $(CFLAGS) -c mylib.c
X
Xclean :
X	rm -f main mylib.o gademo.o
END_OF_FILE
if test 199 -ne `wc -c <'makefile'`; then
    echo shar: \"'makefile'\" unpacked with wrong size!
fi
# end of 'makefile'
fi
if test -f 'mylib.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'mylib.c'\"
else
echo shar: Extracting \"'mylib.c'\" \(3194 characters\)
sed "s/^X//" >'mylib.c' <<'END_OF_FILE'
X
X#include <stdio.h>
X#include <math.h>
X#include <malloc.h>
X
Xvoid FatalError(char *error_text)
X{
X	fprintf(stderr,"Run-time error...\n");
X	fprintf(stderr,"%s\n",error_text);
X	fprintf(stderr,"...now exiting to system...\n");
X	exit(1);
X}
X
Xdouble gasdev(int *idum)
X{
X	static int iset=0;
X	static double gset;
X	double fac,r,v1,v2;
X	double ran1();
X
X	if  (iset == 0) {
X		do {
X			v1=2.0*ran1(idum)-1.0;
X			v2=2.0*ran1(idum)-1.0;
X			r=v1*v1+v2*v2;
X		} while (r >= 1.0 || r == 0.0);
X		fac=sqrt(-2.0*log(r)/r);
X		gset=v1*fac;
X		iset=1;
X		return v2*fac;
X	} else {
X		iset=0;
X		return gset;
X	}
X}
X
X#define M1 259200
X#define IA1 7141
X#define IC1 54773
X#define RM1 (1.0/M1)
X#define M2 134456
X#define IA2 8121
X#define IC2 28411
X#define RM2 (1.0/M2)
X#define M3 243000
X#define IA3 4561
X#define IC3 51349
X
Xdouble ran1(idum)
Xint *idum;
X{
X	static long ix1,ix2,ix3;
X	static double r[98];
X	double temp;
X	static int iff=0;
X	int j;
X	void nrerror();
X
X	if (*idum < 0 || iff == 0) {
X		iff=1;
X		ix1=(IC1-(*idum)) % M1;
X		ix1=(IA1*ix1+IC1) % M1;
X		ix2=ix1 % M2;
X		ix1=(IA1*ix1+IC1) % M1;
X		ix3=ix1 % M3;
X		for (j=1;j<=97;j++) {
X			ix1=(IA1*ix1+IC1) % M1;
X			ix2=(IA2*ix2+IC2) % M2;
X			r[j]=(ix1+ix2*RM2)*RM1;
X		}
X		*idum=1;
X	}
X	ix1=(IA1*ix1+IC1) % M1;
X	ix2=(IA2*ix2+IC2) % M2;
X	ix3=(IA3*ix3+IC3) % M3;
X	j=1 + ((97*ix3)/M3);
X	if (j > 97 || j < 1) FatalError("RAN1: This cannot happen.\n");
X	temp=r[j];
X	r[j]=(ix1+ix2*RM2)*RM1;
X	return temp;
X}
X
X#undef M1
X#undef IA1
X#undef IC1
X#undef RM1
X#undef M2
X#undef IA2
X#undef IC2
X#undef RM2
X#undef M3
X#undef IA3
X#undef IC3
X
X/* Allocate and deallocate routines which help you think of
Xarrays from index 1: */
X
Xfloat *vector(int nl, int nh)
X{
X	float *v;
X
X	v=(float *)malloc((unsigned) (nh-nl+1)*sizeof(float));
X	if (!v) FatalError("allocation failure in vector()");
X	return v-nl;
X}
X
Xint *ivector(nl,nh)
Xint nl,nh;
X{
X	int *v;
X
X	v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
X	if (!v) FatalError("allocation failure in ivector()");
X	return v-nl;
X}
X
Xdouble *dvector(nl,nh)
Xint nl,nh;
X{
X	double *v;
X
X	v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
X	if (!v) FatalError("allocation failure in dvector()");
X	return v-nl;
X}
X
Xvoid free_vector(v,nl,nh)
Xfloat *v;
Xint nl,nh;
X{
X	free((char*) (v+nl));
X}
X
Xvoid free_ivector(v,nl,nh)
Xint *v,nl,nh;
X{
X	free((char*) (v+nl));
X}
X
Xvoid free_dvector(v,nl,nh)
Xdouble *v;
Xint nl,nh;
X{
X	free((char*) (v+nl));
X}
X
Xvoid mk_freq(int *outcomes, int N, int *frequencies, int CATEGORIES)
X{
X	int i;
X
X	for (i=1; i<=CATEGORIES; i++) frequencies[i] = 0;
X	for (i=1; i<=N; i++) 
X		frequencies[outcomes[i]]++;
X}
X
Xvoid throw_darts(double *cdf, int n, 
X	int *seed, int *outcomes, int numdarts)
X{
X	int i, j;
X	double U;
X
X	for (i=1; i<=numdarts; i++) {
X		U = ran1(seed);
X		for (j=1; j<=n; j++) {
X			if (cdf[j] > U) {
X				outcomes[i] = j;
X				break;
X			}
X		}
X	}
X}
X
Xint pdf2cdf(double *pdf, double *cdf, int n)
X{
X	double sum;
X	int i;
X
X	for (sum = 0.0, i=1; i<=n; i++) {
X		sum += pdf[i];
X		cdf[i] = sum;
X	}
X	if (fabs(sum - 1.0) > 1.0e-2) {
X		printf("Sums to %20.18f, not a pdf!\n", sum);
X		return 1;
X	}
X	return 0;
X}
X
Xvoid cdf2pdf(double *cdf, double *pdf, int n)
X{
X	int i; 
X	double last;
X
X	last = 0.0;
X	for (i=1; i<=n; i++) {
X		pdf[i] = cdf[i] - last;
X		last = cdf[i];
X	}
X}
X
END_OF_FILE
if test 3194 -ne `wc -c <'mylib.c'`; then
    echo shar: \"'mylib.c'\" unpacked with wrong size!
fi
# end of 'mylib.c'
fi
if test -f 'mylib.h' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'mylib.h'\"
else
echo shar: Extracting \"'mylib.h'\" \(590 characters\)
sed "s/^X//" >'mylib.h' <<'END_OF_FILE'
X
X#define FALSE 0
X#define TRUE 1
X
Xvoid FatalError(char *error_text);
X
Xdouble gasdev(int *idum);
Xdouble ran1(int *idum);
X
Xfloat *vector(int nl, int nh);
Xint *ivector(int nl, int nh);
Xdouble *dvector(int nl, int nh);
Xvoid free_vector(float *v, int nl, int nh);
Xvoid free_ivector(int *v, int nl, int nh);
Xvoid free_dvector(double *v, int nl, int nh);
X
Xvoid mk_freq(int *outcomes, int N, int *frequencies, int CATEGORIES);
Xvoid throw_darts(double *cdf, int n, int *seed, int *outcomes, int numdarts);
Xint pdf2cdf(double *pdf, double *cdf, int n);
Xvoid cdf2pdf(double *cdf, double *pdf, int n);
X
END_OF_FILE
if test 590 -ne `wc -c <'mylib.h'`; then
    echo shar: \"'mylib.h'\" unpacked with wrong size!
fi
# end of 'mylib.h'
fi
echo shar: End of shell archive.
exit 0

-- 
_______________________________________________________________________________
Ajay Shah, (213)734-3930, ajayshah@usc.edu
                             The more things change, the more they stay insane.
_______________________________________________________________________________

exit 0 # Just in case...