Page 144
Java: Simulated Annealing for the TSP
import java.io.*;
import java.util.*; // has class Random
public class sa
{
static String helptext[]={
"\n Simulated Annealing for a TSP using partial path reversal (ppr)",
"",
"Usage:\n\t prog [-v freq] [-i iterations] CityDataBase",
"",
"Note, the last argument is the city database file",
"Database consists of # cities as first line, followed by x,y",
"coordinates of the cities and that optionally followed by a",
"label of up to 3 characters.",
"freq = frequency of reports (e.g. 10000)",
"iterations= ",
};
private static final int NTRIALS=1000000; //for now
/*--- ggglobal var ---*/
static CITY[] city; // city array, start with 1
static int NC;
static int[] tour, modtour;
static int nTrials=NTRIALS;
static int vfreq = 0;
static Random rndm;
//mmm===
public static void main( String[] args)
{
getinputs(args);
long seed;
//seed ran. num. gen. against making random maps
seed = System.currentTimeMillis();
// seed = 999402363;
System.out.println("seed= "+seed);
rndm = new Random(seed);
/*
//this prints the labels in "identity" order
if( false )
{
for( int i=0; i < tour.length; ++i ) tour[i] = i;
print_tour(tour);
}
*/
int[] t = new int[NC];
double dE, E;
double r, bestE;
// double d=2200, a = 0.03, T=73330; //for avg |dE| = 737
double d=190, a=.001, T=99500; //for avg |dE|=1000, #its=10^6
// double d=3000, a=.03, T=99500; //for avg |dE|=1000, #its=10^6
mkperm(NC,tour); //makeup a random permutation
bestE = E = tourLength(tour);
int its=0;
/*
double mean = 0;
for( ; its < 40; ++its)
{
modtour = ppr(tour);
dE = tourLength(modtour) - E;
mean += Math.abs(dE);
System.out.println("dE= "+dE);
for( int i=0; i < NC; ++i) tour[i] = modtour[i];
}
System.out.println("avg absolute mean: "+mean/its);
System.exit(0);
*/
System.out.println("E0= "+bestE);
print_tour(tour);
//hhh
for( its=0; its < nTrials; ++its )
{
modtour = ppr(tour);
dE = tourLength(modtour) - E;
r = Math.exp(-dE/T);
if( rndm.nextDouble() < r )
{
for( int i=0; i < NC; ++i) tour[i] = modtour[i];
E = E + dE;
}
if( E < bestE)
{
bestE = Math.min(bestE,E);
for( int i=0; i < t.length; ++i) t[i] = tour[i];
System.out.print("iteration "+its+", ");
System.out.println(
"new best: "+bestE+" (cf "+tourLength(t)+"), tour:");
// print_tour(t);
}
T = d/(a + Math.log(its));
}
System.out.println("Done, #iterations= "+its+", best= "+bestE);
print_tour(t);
}
//iii===
static void
getinputs(String[] argv)
{
int i,j;
int oo;
int[] o = new int[1];
char s;
int pc;
String ipstring;
// StringBuffer ipstring = new StringBuffer();
boolean labelsset=false;
int k,p;
// see if&get command line arg
i = argv.length;
int v=0; //index into argv
while( i-- > 0 && argv[v].charAt(0) == '-') {
s = argv[v].charAt(1);
switch( s ) {
case 'v':
//System.out.println("v= "+v+" argv[v]= "+argv[v]+", argv[v+1]= "+argv[v+1]);
if( argv[v].length() > 2 ) { //argument packed with -v
vfreq = Integer.parseInt(argv[v].substring(2));
}
else { //-v might have a detached argument or none at all
if( i==0 ) { //no more arguments, an error
break;
}
//check that first char of next arg is digit
if( Character.isDigit(argv[v+1].charAt(0)) )
{ //read it
vfreq = Integer.parseInt(argv[++v]);
--i;
}
}
break;
case 'i':
if( argv[v].length() > 2 ) { //argument packed with -i
nTrials = Integer.parseInt(argv[v].substring(2));
}
else { //-i might have a detached argument or none at all
if( i==0 ) { //no more arguments, an error
break;
}
//check that first char of next arg is digit
if( Character.isDigit(argv[v+1].charAt(0)) )
{ //read it
nTrials = Integer.parseInt(argv[++v]);
--i;
}
}
break;
default:
if( ! argv[v].equals("-h") )
System.out.println("Illegal option "+argv[v]);
for( int l=0; l < helptext.length; ++l )
System.out.println(helptext[l]);
System.exit(0);
}
}
String inputFile = argv[argv.length-1];
try {
// Open a reader on the file.
FileInputStream stream = new FileInputStream(inputFile);
//InputStream stream = getClass().getResourceAsStream("/"+inputFile);
BufferedReader reader =
new BufferedReader(new InputStreamReader(stream));
ipstring = reader.readLine();
o[0]=0; oo = numpars(o,ipstring);
NC = Integer.parseInt(ipstring.substring(o[0],oo));
//----- memory for ------
city = new CITY[NC+1]; //+1 so first=last
char quoteChar;
//get city coordinates
for(j=0; j < NC; ){
city[j] = new CITY(); //default constructor
ipstring = reader.readLine();
for( o[0] = 0; ; ) {
if( (oo = numpars(o,ipstring)) < 0 ) break;
city[j].x = Double.parseDouble(ipstring.substring(o[0],oo));
o[0] = oo;
if( (oo = numpars(o,ipstring)) < 0 ) break;
city[j].y = Double.parseDouble(ipstring.substring(o[0],oo));
o[0] = oo;
//now scan for a label, otherwise make the label the index
if( (pc=ipstring.indexOf('"',oo)) >= 0 )
quoteChar = '"';
else if( (pc=ipstring.indexOf('\'',oo)) >= 0 )
quoteChar = '\'';
else quoteChar = 'z';
if( quoteChar != 'z' )
{ //find " or ' in ipstring
labelsset = true;
o[0] = pc+1; //start of label
if( (pc=ipstring.indexOf(quoteChar,o[0])) < 0 )
{ //error
System.out.println("Error, no end to label for: "+ipstring);
System.exit(0);
}
for(k=0; k < (pc-o[0]); ++k ) {
city[j].label[k] = ipstring.charAt(o[0]+k);
}
//city[j].label[k]=0;
o[0] = pc; //for next round
}
else { //make labels = index (or index +1) if -1 is set
if( labelsset ) {
System.out.println("no label for "+ipstring);
System.exit(0);
}
city[j].label = (Integer.toString(j+1)).toCharArray();
}
if (++j==NC) break;
}
}
// Close the file.
reader.close();
stream.close();
}
catch( IOException e ) {
System.out.println(e.getMessage());
System.exit(0);
}
//set last city equal to the first
city[NC] = new CITY();
city[NC].x = city[0].x; city[NC].y = city[0].y;
tour = new int[NC];
modtour = new int[NC];
tour[0] = modtour[0] = 0;
}//getinputs
//===
static void
print_tour(int[] tour)
{
int i,j;
int len = (city[0].label).length;
for( i=0; i < NC; ++i) {
System.out.print((new String(city[tour[i]].label))+" ");
}
System.out.println();
}
//===
static void
COPY_PERM(int[] modtour,int[] wbest)
{
for(int i=0;i < NC;++i) wbest[i]=modtour[i];
}
//===
static int
numpars(int[] pi, String s)
/*::::
scan string s for the next number starting from index *pi, advance
*pi to the start of the number. Return value is -1 if no number
encountered or the first location following the number otherwise.
mod: fix for scientific notation, treat embedded E as a "digit"
mod:(12/00) treat leading # as a comment, acts as if end of line
::::*/
{
int i;
//System.out.println("numpars: pi: "+pi[0]+", s: "+s);
//check for error input (negative index)
if( pi[0] < 0 || s.length() == 0 ) return -1;
//check for comment line, begins with #
if( pi[0]==0 && s.charAt(0) == '#' ) return -1;
/*--- find beginning of a number ---*/
for(i = pi[0];;++i) {
//System.out.print("(1)i now "+i);
//if( i < s.length() ) System.out.println(" reading "+s.charAt(i));
//else System.out.println(" last char read was "+s.charAt(i-1));
if( s.length() == i ||
s.charAt(i) == '\n' || s.charAt(i) == '\r' ) return -1;
if( s.charAt(i) >= '0' && s.charAt(i) <= '9' ) break;
if( s.charAt(i) == '+' || s.charAt(i) == '-' || s.charAt(i) == '.' ) break;
}
pi[0] = i;
/*--- find extent of the number, E followed by + or - is ok ---*/
for(++i;;++i) {
//System.out.print("(2)i now "+i);
//if( i < s.length() ) System.out.println(" reading "+s.charAt(i));
//else System.out.println(" last char read was "+s.charAt(i-1));
if( i >= s.length() ) return i;
if( s.charAt(i) >= '0' && s.charAt(i) <= '9' ) continue;
/*--- possible error here, should "continue" on - too ---*/
if( s.charAt(i) == '.' ) continue;
if( s.charAt(i) == 'E' || s.charAt(i) == 'e' )
if( s.charAt(i+1) == '+' || s.charAt(i+1) == '-' ) { ++i; continue;}
return i;
}
}
//sss===
static double
tourLength(int[] tour) // calculate tour length
{
double dx,dy,w=0;
int i, k=0;
//print_tour(tour);
for( i=1; i < tour.length; ++i ) {
dx = -city[tour[k]].x; dy = -city[tour[k]].y;
++k;
dx += city[tour[k]].x; dy += city[tour[k]].y;
w += Math.sqrt(dx*dx + dy*dy);
//System.out.println("link "+tour[k-1]+" to "+tour[k]+" d= "+
// Math.sqrt(dx*dx + dy*dy)+", sum= "+w);
} // pointer at N-1 when done
// now get back to first city
dx = -city[tour[k]].x; dy = -city[tour[k]].y;
dx += city[tour[0]].x; dy += city[tour[0]].y;
w += Math.sqrt(dx*dx + dy*dy);
//System.out.println("link "+tour[k]+" to "+tour[0]+" d= "+
// Math.sqrt(dx*dx + dy*dy)+", sum= "+w);
return w;
}
//===
static int[]
ppr(int[] tour)
{
int i, j;
int k, l;
i = 1 + (int)((NC-2)*rndm.nextDouble());
j = i+1 + (int)((NC-i-1)*rndm.nextDouble());
// make up modtour by reversing path from i to j
for( k=1; k < i; ++k ) modtour[k] = tour[k];
for( l=j; ; ) {
modtour[k++]=tour[l];
if( --l < i ) break;
}
for( ; k < NC; ++k ) modtour[k] = tour[k];
return modtour;
}
//===
static void
mkperm(int size, int[] perm)
{
int i,k,r;
for( i=1; i < size; ++i) perm[i] = i;
for( i=1; i < size-1; ++i)
{
k = 1+ (int)((size-i)*rndm.nextDouble());
r = perm[k];
for( ; k < size-i; ++k) perm[k] = perm[k+1];
perm[k] = r;
}
}
}
//kkk===
class CITY
{
double x,y;
char[] label = {' ',' ',' ',' '};
}
Page 147
Java: Genetic Algorithms, Complete Population Turnover Style
import java.io.*;
import java.util.*;
public class shekelGBS
{
static String[] helptext =
{
"Usage:\n\t <invoke program>",
""};
//command line variables
static int vfreq=5; //vfreq=0;
//hold overs we may yet implement
static final int NCEXIT=80000; // no change exit condition
//pertains to all instances
static int nIterations=50;
static int nTrials = 6;
static Random r;
static long seed;
static ArrayList minPerformance = new ArrayList();
static ArrayList maxPerformance = new ArrayList();
//=== problem specifics
static final int nBits = 10; //chromo length
static final int nStates = -1 +
(int)Math.rint(Math.exp(nBits*Math.log(2))); //0..2^nBits-1
static final double range = 10;
static double[][] S={
{7.382, 10.295, 0.1613},
{1.972, 4.722, 0.2327},
{8.111, 10.928, 0.2047},
{2.714, 11.052, 0.1579},
{4.022, 3.820, 0.1268},
{5.480, 12.911, 0.0519},
{8.316, 7.974, 0.1340},
{7.771, 8.085, 0.1040},
{1.618, 4.047, 0.1719},
{6.918, 5.310, 0.0665},
{8.571, 13.860, 0.0838},
{9.517, 6.163, 0.2204}};
//max value: 19.32926329141593 at 5.48
/*
{7.928, 11.380, 0.1515},
{7.480, 4.714, 0.1339},
{3.739, 7.692, 0.0900},
{1.535, 13.688, 0.0715},
{2.389, 11.687, 0.2139},
{8.236, 9.606, 0.1788},
{3.003, 12.746, 0.1783},
{5.249, 9.565, 0.1780},
{2.132, 9.786, 0.0551},
{3.773, 3.023, 0.1395},
{0.620, 13.985, 0.1181},
{3.929, 13.905, 0.2164}};
//max value: 18.349440998358663 at 2.132
*/
//=== genetic specifics
int popSize = 12;
double pm; //mutation probability, set via constructor
double geoRate = 1.0; //pm reduction rate, default to none
//ggg===
critterSG[] c0 = new critterSG[popSize]; //ptrs to old critters
critterSG[] c1 = new critterSG[popSize]; //ptrs to new critters
critterSG bestCrit = new critterSG(this);
long fevals=0L; /* keep track of fcn evaluations */
boolean logProg = false;
double logBase = 0; //only log fitnesses exceeding this
//sss===
double
runGA(ArrayList runRecord)
{
int t = 1; //algorithm time for cooling/stopping purposes
int noImprov=0; //no improvement count
bestCrit.fitness = 0;
//INITIALIZATION, create the population
for(int i=0; i < popSize; ++i )
{
c0[i] = new critterSG(this,r.nextInt(nStates+1)); //random init
c1[i] = new critterSG(this);
//System.out.println(
//"critter "+i+", x= "+c[i].decode()+" fitness: "+c[i].fitness);
}
bestCrit.state = c0[0].state; bestCrit.fitness = c0[0].fitness;
//EVOLUTION LOOP
//initial mututation rate is 1
double[] breakPts = new double[popSize];
double sum, w; //double helpers
int k, d, m, mc, x, y; //int helpers
critterSG[] temp; //needed for relay
int alIndex = 0; //record min/max progress periodically
for( t=1; t <= nIterations; ++t)
{
if( vfreq!=0 && t%vfreq == 0 )
{
double min =
((Double)shekelGBS.minPerformance.get(alIndex)).doubleValue();
double max =
((Double)shekelGBS.maxPerformance.get(alIndex)).doubleValue();
//System.out.format(
//"before: t= %2d, alIndex= %3d, min= %.3f, max= %.3f, best= %.3f\n",
//t,alIndex,min,max,bestCrit.fitness);
if( bestCrit.fitness < min )
shekelGBS.minPerformance.set(alIndex,(Double)bestCrit.fitness);
if( bestCrit.fitness > max )
shekelGBS.maxPerformance.set(alIndex,(Double)bestCrit.fitness);
/*
min = ((Double)shekelGBS.minPerformance.get(alIndex)).doubleValue();
max = ((Double)shekelGBS.maxPerformance.get(alIndex)).doubleValue();
System.out.format(
"after: t= %2d, alIndex= %3d, min= %.3f, max= %.3f, best= %.3f\n",
t,alIndex,min,max,bestCrit.fitness);
*/
++alIndex;
}
//do roulette wheel selection
sum=0;
for(int i=0; i < popSize; ++i) //make up wheel
{
sum += c0[i].fitness;
breakPts[i] = sum;
//System.out.println("breakPt["+i+"]= "+sum);
}
//select popSize times & add to new critters
for( int i=0; i < popSize; ++i)
{
w = sum*r.nextDouble(); // 0 <= w < sum, provided sum > 0
for( k=0; w>=breakPts[k]; ++k) {} //k selected
//System.out.println("w= "+w+" break selected "+(k));
c1[i].state = c0[k].state;
c1[i].fitness = c0[k].fitness;
c1[i].x = c0[k].x;
/*
System.out.println(" gen. "+t+", break points:");
for( int g=0; g < popSize; ++g)
System.out.println("breakPt["+g+"]= "+breakPts[g]);
System.out.println("w= "+w+" break selected "+(k));
*/
}
/*
if( vfreq!=0 && t%vfreq == 0 ){
System.out.println("post RWS pop. for generation "+(t));
for( int i=0; i < popSize; ++i) {
System.out.println(""+c1[i].state+" fit: "+c1[i].fitness);}
}
*/
//crossover critters 0 with 1, 2 with 3, etc
for( int i=0; i < popSize; i += 2 )
{
//select divide point, to right
d = 1+r.nextInt(nBits-1); //d=0 or d=nBits just name change
//bit mask and complement
m = (~0 << d) & nStates;
mc = ~m & nStates;
//System.out.println("Xover at "+d+" mask: "+m+" complement: "+mc);
x = c1[i].state; y = c1[i+1].state;
c1[i].state = (x & m) | (y & mc);
c1[i+1].state = (x & mc) | (y & m);
c1[i].setFitness();
c1[i+1].setFitness();
//System.out.println(
//"off1 state= "+c1[i].state+" fit: "+c1[i].fitness);
//System.out.println(
//"off2 state= "+c1[i+1].state+" fit: "+c1[i+1].fitness);
}
//check for improvement from mating
for( int i=0; i < popSize; ++i)
{
if( c1[i].fitness > bestCrit.fitness )
{
//copy c1[i] into best
bestCrit.state = c1[i].state;
bestCrit.x = c1[i].x;
bestCrit.fitness = c1[i].fitness;
noImprov=0;
//check whether to print solution
if( logProg && bestCrit.fitness > logBase )
{
//report(1,0,t); //1=only print value, 0=only for best
logBase = bestCrit.fitness;
}
}
}
//mutate each bit of each critter
//System.out.println("cur. pm "+pm);
for( int i=0; i < popSize; ++i )
{
for( int j=0; j < nBits; ++j)
{
if( r.nextDouble() < pm )
{
m = (1 << j);
//System.out.println("bit "+j+" selected, state before: "+c1[i].state+
//" state after: "+(m^c1[i].state));
c1[i].state ^= m;
}
}
c1[i].setFitness();
}
//check for improvement from mutation
for( int i=0; i < popSize; ++i)
{
if( c1[i].fitness > bestCrit.fitness )
{
//copy c1[i] into best
bestCrit.state = c1[i].state;
bestCrit.x = c1[i].x;
bestCrit.fitness = c1[i].fitness;
noImprov=0;
//check whether to print solution
if( logProg && bestCrit.fitness > logBase )
{
//report(1,0,t); //1=only print value, 0=only for best
logBase = bestCrit.fitness;
}
}
}
runRecord.add((Double)bestCrit.fitness);
//now c1 is the next population
temp = c0;
c0 = c1;
c1 = temp;
pm *= geoRate; //update the mutation rate
++noImprov; //update no improvement count
}
//report(2,0,t); //2= print value & matrix, 0= only for best
//System.out.println(
//"ending best: "+bestCrit.fitness+" @ x = "+bestCrit.x);
return bestCrit.fitness;
}
//===
double
objective( double x)
{
++fevals;
double y=0;
for( int i=0; i < S.length; ++i )
{
y += 1/(S[i][1]*S[i][1]*(x-S[i][0])*(x-S[i][0])+S[i][2]);
//System.out.println("i= "+i+", y= "+y);
}
return y;
}
//ccc===
public
shekelGBS(double pm)
{
this.pm = pm;
}
//ccc===
public
shekelGBS(double pm, double gr)
{
this.pm = pm;
geoRate = gr;
}
//ccc===
public
shekelGBS(double pm, double gr, double logBase)
{
this.pm = pm;
geoRate = gr;
logProg = true;
this.logBase = logBase;
}
//rrr===
void
report( int rLevel, int todo, long its)
{
/*::::::
rLevel= extent to report:
0 minimal (just dots),
1 just fitness values
2 fitness values and entire matrix
todo= # from top of gene pool to print
its iteration count
::::::*/
if( rLevel == 0 ) { System.out.print("."); return;}
if( rLevel == 1 )
{
System.out.format(
"best= %f, #its= %d\n",bestCrit.fitness, its);
for( int i=0; i < todo; ++i)
System.out.format("critter %2d fitness= %f\n",i,c1[i].fitness);
return;
}
if( rLevel == 2 )
{
System.out.format(
"best= %f, #its= %d\n", bestCrit.fitness,its);
for( int i=0; i < todo; ++i)
{
System.out.format("critter %2d fitness= %f\n",i,c1[i].fitness);
}
}
}
//mmm===
public static void
main(String[] args)
{
seed = System.currentTimeMillis();
// seed = 1190902532283L; // uncomment for testing
r = new Random(seed);
System.out.println("\nseed: "+seed);
double geoRate=1.00;
double pm = 0.1;
if( args.length >= 1 )
{
pm = Double.parseDouble(args[0]);
if( args.length >= 2 )
geoRate = Double.parseDouble(args[1]);
}
int successCount=0;
double bestOverAll = 0;
double runOut;
shekelGBS pprob;
ArrayList[] runLogs = new ArrayList[nTrials];
//initialize the min/max record
for(int i=1; i <= nIterations; ++i)
{
if( vfreq!=0 && i%vfreq == 0 )
{
minPerformance.add((Double)1000000.);
maxPerformance.add((Double)0.);
}
}
for( int i=0; i < nTrials; ++i)
{
if( i==0 )
pprob = new shekelGBS(pm,geoRate);
else
pprob = new shekelGBS(pm,geoRate,bestOverAll);
runLogs[i] = new ArrayList();
runOut = pprob.runGA(runLogs[i]);
if( runOut > bestOverAll ) bestOverAll = runOut;
System.out.println("trial "+i+", result= "+runOut);
//if( pprob.runGA() != 0 ) ++successCount;
}
System.out.println(
"in "+nTrials+" trials, best overall= "+bestOverAll);
for( int i=0; i < nTrials; ++i)
{
System.out.println("run "+i+" record");
for( int t=0; t < runLogs[0].size(); ++t)
{
System.out.format("t= %3d, best= %.3f\n",
(t+1),((Double)runLogs[i].get(t)).doubleValue());
}
}
System.out.println(
"min/max record every "+vfreq+" over "+nTrials+" runs");
for( int j=0; j < minPerformance.size(); ++j)
{
System.out.format(
"j= %4d, min %.3f, max %.3f\n",(j+1)*vfreq,
((Double)minPerformance.get(j)).doubleValue(),
((Double)maxPerformance.get(j)).doubleValue());
}
// System.out.println("geo. rate= "+geoRate);
//", success rate: "+successCount/(double)nTrials);
}
}
//kkk===
class critterSG //Shekel problem critter
{
//state holds the bits
int state=0; //0..2^nBits-1 = nStates, divide to get the range
double x; // x = range*(state/nStates)
double fitness=0;
shekelGBS sgb;
int[] chromo; //no longer needed
//ccc===
public
critterSG(shekelGBS s)
{
sgb = s; state = 0; setFitness();
}
//===
public
critterSG(shekelGBS s,int k)
{
sgb = s; state = k; setFitness();
}
//===
/**
* decode: from chromosome to state
* obsolete, no longer needed
*/
int
decode()
{
state = chromo[0];
//System.out.println("starting state = chromo[0] = "+state);
int base = 2;
for( int i=1; i < shekelGBS.nBits; ++i)
{
state += base*chromo[i];
//System.out.println("chromo["+i+"]= "+chromo[i]+", state now= "+state);
base *= 2;
}
return state;
}
//===
/**
* encode: from state to chromosome
* obsolete, no longer needed
*/
void
encode()
{
int k = state;
//from least significant to most significant
for( int i=0; i < shekelGBS.nBits; ++i)
{
chromo[i] = k%2;
k /= 2;
//System.out.println("encode: k= "+k+", chromo["+i+"]]= "+chromo[i]);
}
}
//===
void
setFitness()
{
x = sgb.range*(state/(double)shekelGBS.nStates);
fitness = sgb.objective(x);
}
}
Page 148
Java: Genetic Algorithms, Sorted Rank Style
import java.io.*;
import java.util.*;
public class shekelRSS
{
static String[] helptext =
{
"Implementation of the RS style GA for the shekel problem.",
"Usage:\n\t < invoke program>",
""};
//command line variables
static int vfreq=5; //vfreq=0;
//hold overs we may yet implement
static final int NCEXIT=80000; // no change exit condition
//pertains to all instances
static int nIterations=50;
static int nTrials = 6;
static Random r;
static long seed;
static ArrayList minPerformance = new ArrayList();
static ArrayList maxPerformance = new ArrayList();
static final int BUGINVALID=-1; //mark popspace free
static final double tempSTART=0.88;
static final int STDZ=12; //population size
static final int HALFSTDZ=STDZ/2; //half std pop size
static final int THREEQSTD=3*STDZ/4; // 3 quarters std pop size
static final int POPSPACE=3*STDZ; //space for critters
static final int MAXMATE=THREEQSTD; //max number of matings
static final int MAXMUTE=HALFSTDZ; // max number of mutations
//=== problem specifics
static final int nBits = 20; //chromo length
static final int nStates = -1 +
(int)Math.rint(Math.exp(nBits*Math.log(2))); //0..2^nBits-1
static final double range = 10;
static double[][] S={
{7.382, 10.295, 0.1613},
{1.972, 4.722, 0.2327},
{8.111, 10.928, 0.2047},
{2.714, 11.052, 0.1579},
{4.022, 3.820, 0.1268},
{5.480, 12.911, 0.0519},
{8.316, 7.974, 0.1340},
{7.771, 8.085, 0.1040},
{1.618, 4.047, 0.1719},
{6.918, 5.310, 0.0665},
{8.571, 13.860, 0.0838},
{9.517, 6.163, 0.2204}};
//max value: 19.32926329141593 at 5.48
//ggg===
critterSP[] c0 = new critterSP[POPSPACE]; //ptrs to critters
critterSP bestCrit = new critterSP(this);
int[] role = new int[POPSPACE];
long fevals=0L; /* keep track of fcn evaluations */
boolean logProg = false;
double logBase = 0; //only log fitnesses exceeding this
//sss===
double
runGA(ArrayList runRecord)
{
int t = 1; //algorithm time for cooling/stopping purposes
int z=0; //current population size
int noImprov=0; //no improvement count
bestCrit.fitness = 0;
double[] af = new double[1];
double temp; //temperature, using as RR living prob.
//mark population slots free
for( int i=0; i < POPSPACE; ++i )
{
c0[i] = new critterSP(this);
c0[i].fitness = BUGINVALID;
}
//INITIALIZATION, create the population
af[0]=0;
z=0;
for(int i=0; i < STDZ; ++i )
{
c0[i].state = r.nextInt(nStates+1); //random init
c0[i].setFitness();
z = merge(-1,i,af,z); /* sense = -1 for maximization */
}
//EVOLUTION LOOP
int nmate = MAXMATE; //currently 6
int nmute = MAXMUTE; //current 6
int dltz=0;
int ni, k; //int helpers
int alIndex = 0;
outerLoop: //block label
for( t=1; t < = nIterations; ++t)
{
if( vfreq!=0 && t%vfreq == 0 )
{
double min =
((Double)shekelRSS.minPerformance.get(alIndex)).doubleValue();
double max =
((Double)shekelRSS.maxPerformance.get(alIndex)).doubleValue();
//System.out.format(
//"before: t= %2d, alIndex= %3d, min= %.3f, max= %.3f, best= %.3f\n",
//t,alIndex,min,max,bestCrit.fitness);
if( bestCrit.fitness < min )
shekelRSS.minPerformance.set(alIndex,(Double)bestCrit.fitness);
if( bestCrit.fitness > max )
shekelRSS.maxPerformance.set(alIndex,(Double)bestCrit.fitness);
/*
min = ((Double)shekelRSS.minPerformance.get(alIndex)).doubleValue();
max = ((Double)shekelRSS.maxPerformance.get(alIndex)).doubleValue();
System.out.format(
"after: t= %2d, alIndex= %3d, min= %.3f, max= %.3f, best= %.3f\n",
t,alIndex,min,max,bestCrit.fitness);
*/
++alIndex;
}
//hhh
for( int i=0; i < nmute; ++i)
{
if( (ni = getslot()) < 0 ) break; //no free slots
//pick random one to mutate
k = r.nextInt(z);
z = mutate(c0[role[k]],ni,af,z);
}
//mating (cross-over) step
// mate top half, find mate at random (among current pop.)
for( int m=0; m < nmate; ++m )
{
if((ni = getslot()) < 0 ) break; //no free slots
//random population member
k = r.nextInt(z);
if( m < nmate-1 )
{
z = mate(c0[role[m]],c0[role[k]],ni,af,z);
}
else //do wrong-side-of-tracks for the last one
z = mate(c0[role[m]],new critterSP(this),ni,af,z);
}
//check for improvement
if( c0[role[0]].fitness > bestCrit.fitness )
{
//copy c0[role[0]] into best
bestCrit.state = c0[role[0]].state;
bestCrit.x = c0[role[0]].x;
bestCrit.fitness = c0[role[0]].fitness;
noImprov=0;
//check whether to print solution
if( logProg && bestCrit.fitness > logBase )
{
//report(1,0,t); //1=only print value, 0=only for best
logBase = bestCrit.fitness;
}
}
//check for no change exit
if( ++noImprov == NCEXIT ) break;
//select down the population size
for( dltz=STDZ-z; dltz < 0; ++dltz)
z = removal(gauntlet(tempSTART,z),af,z);
runRecord.add((Double)bestCrit.fitness);
}
//report(2,0,t); //2=print value & matrix, 0=only for best
//System.out.println(
//"ending best: "+bestCrit.fitness+" @ x = "+bestCrit.x);
return bestCrit.fitness;
}
//===
double
objective( double x)
{
++fevals;
double y=0;
for( int i=0; i < S.length; ++i )
{
y += 1/(S[i][1]*S[i][1]*(x-S[i][0])*(x-S[i][0])+S[i][2]);
//System.out.println("i= "+i+", y= "+y);
}
return y;
}
//ccc===
public
shekelRSS()
{
}
//ccc===
public
shekelRSS(double logBase)
{
logProg = true;
this.logBase = logBase;
}
//===
public int
merge(int sense,int i, double[] af, int cs)
/*::::::::
merge critter i into role according to rank via fitness
sense= +1, small to big, i.e. minimization
sense= -1, big to small, maximization
i= index of this critter
af= average fitness, cs= present population size
return new population size
:::::::*/
{
int m, k, j;
int z = cs;
//insert newbug according to its fitness
for( m=0; ; ++m) {
if( m == z ) break;
if( sense*c0[role[m]].fitness > sense*c0[i].fitness ) break;
}
//goes into slot m, see if room and update avg fitness
af[0] *= z; //sum of fitnesses
if( z < POPSPACE ) //have room, so add
{
++z;
af[0] = (af[0] + c0[i].fitness)/z;
}
else { //no room, knock off last
af[0] += c0[i].fitness - c0[POPSPACE-1].fitness;
af[0] = (af[0] + c0[i].fitness)/z;
}
for( j=z-2; j>=m; j--) role[j+1] = role[j];
role[m] = i;
return z;
}
//===
public int
getslot()
{
int i;
for(i=0; i < POPSPACE; ++i )
if( c0[i].fitness == BUGINVALID ) return i;
return -1;
}
//rrr===
void
report( int rLevel, int todo, long its)
{
/*::::::
rLevel= extent to report:
0 minimal (just dots),
1 just fitness values
2 fitness values and entire matrix
todo= # from top of gene pool to print
its iteration count
::::::*/
if( rLevel == 0 ) { System.out.print("."); return;}
if( rLevel == 1 )
{
System.out.format(
"best= %.3f, #its= %d\n",bestCrit.fitness,its);
for( int i=0; i < todo; ++i)
System.out.format("critter %4d fitness= %.3f\n",i,c0[i].fitness);
return;
}
if( rLevel == 2 )
{
System.out.format("best= %.3f, #its= %d\n",bestCrit.fitness,its);
for( int i=0; i < todo; ++i)
{
System.out.format(
"critter %4d fitness= %.3f, x= %.3f\n",i,c0[i].fitness,c0[i].x);
}
}
}
//mmm===
public static void
main(String[] args)
{
seed = System.currentTimeMillis();
// seed = 1174741248177L; // uncomment for testing
r = new Random(seed);
System.out.println("\nseed: "+seed);
int successCount=0;
double bestOverAll = 0;
double runOut;
shekelRSS rss;
ArrayList[] runLogs = new ArrayList[nTrials];
//initialize the min/max record
for(int i=1; i <= nIterations; ++i)
{
if( vfreq!=0 && i%vfreq == 0 )
{
minPerformance.add((Double)1000000.);
maxPerformance.add((Double)0.);
}
}
for( int i=0; i < nTrials; ++i)
{
if( i==0 )
rss = new shekelRSS();
else
rss = new shekelRSS(bestOverAll);
runLogs[i] = new ArrayList();
runOut = rss.runGA(runLogs[i]);
if( runOut > bestOverAll ) bestOverAll = runOut;
System.out.println("trial "+i+", result= "+runOut);
}
for( int i=0; i < nTrials; ++i)
{
System.out.println("run "+i+" record");
for( int t=0; t < runLogs[0].size(); ++t)
{
System.out.format("t= %3d, best= %.3f\n",
t,((Double)runLogs[i].get(t)).doubleValue());
}
}
System.out.println(
"min/max record every "+vfreq+" over "+nTrials+" runs");
for( int j=0; j < minPerformance.size(); ++j)
{
System.out.format(
"j= %4d, min %.3f, max %.3f\n",(j+1)*vfreq,
((Double)minPerformance.get(j)).doubleValue(),
((Double)maxPerformance.get(j)).doubleValue());
}
System.out.println(
"in "+nTrials+" trials, best overall= "+bestOverAll);
}
//===
/**
* crossover style mating, keep one offspring randomly
*/
int //return updated value of z
mate( critterSP c1, critterSP c2, int ni, double[] af, int z)
{
int k, d, m, mc, x, y; //int helpers
//select divide point, to right
d = 1+r.nextInt(nBits-1); //d=0 or d=nBits just name change
//bit mask and complement
m = (~0 << d) & nStates;
mc = ~m & nStates;
//System.out.println(
//"Xover at "+d+" mask: "+m+" complement: "+mc);
x = c1.state; y = c2.state;
//sometimes keep front half, sometimes the back half
if( r.nextDouble() < 0.5 )
{
c0[ni].state = (x & m) | (y & mc);
}
else
{
c0[ni].state = (x & mc) | (y & m);
}
c0[ni].setFitness();
//System.out.println(
//"off1 state= "+c1[i].state+" fit: "+c1[i].fitness);
//System.out.println(
//"off2 state= "+c1[i+1].state+" fit: "+c1[i+1].fitness);
return merge(-1,ni,af,z);
}
//===
/**
* flip one bit at random
*/
int
mutate(critterSP crit, int ni, double[] af, int z)
{
//copy crit into slot ni
c0[ni].state = crit.state;
for( int i=0; i < 1; ++i)
{
//bit to flip
int k = r.nextInt(nBits);
int m = (1 << k);
c0[ni].state ^= m;
}
c0[ni].setFitness();
//merge rank
return merge(-1,ni,af,z);
}
//===
int //return rank for reaping
gauntlet(double livprob,int z)
/*::::::
"Russian" roulette scheme, start at the bottom rank and work up.
Independent of max/min sense
::::::*/
{
int kr;
//Establish bug at risk, start at bottom of role
for( kr = z-1; ; )
{
if( r.nextDouble() > livprob) break; //livprob= prob. of surviving
if( kr == 0)
{
//cycled, tap random from bottom 2/3
kr = (int)(0.66*z*r.nextDouble()); //first top 2/3
break;
}
kr--;
}
return kr;
}
//===
int //return new population size
removal(int kr, double[] af, int cs)
/*:::::
sense independent, works for maximization or minimization, uses
rank according to the role.
kr; rank for reaping
cs; colony size
:::::*/
{
int ki,i;
int z=cs;
ki = role[kr]; // killed bug index, fitness rank kr
//remove bug ki = role[kr]
af[0] *= z; /* this is total fitness */
--z;
//following should work even if fitness is integer
af[0] = (af[0] - c0[ki].fitness)/z;
//must re-merge
c0[ki].fitness = BUGINVALID; //invalid mark for index ki
for( i=kr; i < z; ++i ) role[i] = role[i+1];
return z;
}
}
//kkk===
class critterSP //Shekel problem critter
{
//state holds the bits
int state=0; //0..2^nBits-1 = nStates, divide to get the range
double x; // x = range*(state/nStates)
double fitness=0;
shekelRSS srs;
int[] chromo;
//ccc===
public
critterSP(shekelRSS s)
{
srs = s;
state = 0; setFitness();
}
//===
public
critterSP(shekelRSS s,int k)
{
srs = s;
state = k; setFitness();
}
//===
/**
* decode: from chromosome to state
*/
int
decode()
{
state = chromo[0];
int base = 2;
for( int i=1; i < shekelRSS.nBits; ++i)
{
state += base*chromo[i];
base *= 2;
}
return state;
}
//===
/**
* encode: from state to chromosome
*/
void
encode()
{
int k = state;
//from least significant to most significant
for( int i=0; i < shekelRSS.nBits; ++i)
{
chromo[i] = k%2;
k /= 2;
}
//System.out.println("encode: k= "+k+", chromo["+i+"]]= "+chromo[i]);
}
//===
void
setFitness()
{
x = srs.range*(state/(double)shekelRSS.nStates);
fitness = srs.objective(x);
}
}
Page 154
Matlab: Computing the Permanent
% mat is an nxn 0/1 matrix
% used is an n-dim. array
function y=permrecurse(level,n,used,mat,v)
if level == n+1
y=v+1;
else
for i=1:n
if used(i) == 0 & mat(level,i) ~= 0
used(i)=1;
v = permrecurse(level+1,n,used,mat,v);
used(i)=0;
end
end
y=v;
end
Page 155
Fortran: Non-recursive 0/1 Permanent Algorithm
c23456789
PROGRAM PERMNR
C
INTEGER N,A(15,15)
INTEGER PER, PERM(15), OD(15)
COMMON A,PERM,OD,PER
C
N=3
DATA A(1,1),A(1,2),A(1,3)/1,0,1/
DATA A(2,1),A(2,2),A(2,3)/0,1,1/
DATA A(3,1),A(3,2),A(3,3)/1,1,1/
C
DO 10, I=1,N
10 PERM(I)=I
CALL ODOMETER(N)
C the permanent of A is now in PER
print *,'PER= ',PER
STOP
END
C
SUBROUTINE CP(N)
INTEGER N,A(15,15)
INTEGER PER, PERM(15), OD(15)
COMMON A,PERM,OD,PER
DO 10, I=1,N
IF (A(I,PERM(I)).EQ.0) RETURN
10 CONTINUE
PER = PER+1
RETURN
END
C
SUBROUTINE ODOMETER(N)
INTEGER N,A(15,15)
INTEGER PER, PERM(15), OD(15)
COMMON A,PERM,OD,PER
DO 10, I=1,N
10 OD(I) = 0
I = N-1
20 CONTINUE
IF (OD(I) .GT. N-I) THEN
OD(I) = 0
I=I-1
IF (I .LT. 1) GOTO 100
ELSE
K=PERM(I)
PERM(I)=PERM(I+1)
PERM(I+1)=K
IF (OD(I) .EQ. N-I) THEN
DO 30, J=I+1,N-I
K=PERM(J)
PERM(J) = PERM(J+1)
30 PERM(J+1) = K
END IF
CALL CP(N)
I=N-1
END IF
OD(I)=OD(I)+1
GOTO 20
100 RETURN
END
c234567