/*
Copyright 2006 by Sean Luke and George Mason University
Licensed under the Academic Free License version 3.0
See the file "LICENSE" for more information
*/
package ec.app.ecsuite;
import ec.util.*;
import ec.*;
import ec.simple.*;
import ec.vector.*;
/*
* ECSuite.java
*
* Created: Thu Mar 22 16:27:15 2001
* By: Liviu Panait and Sean Luke
*/
/*
* @author Liviu Panait and Sean Luke and Khaled Talukder
* @version 2.0
*/
/**
Several standard Evolutionary Computation functions are implemented.
As the SimpleFitness is used for maximization problems, the mapping f(x) --> -f(x) is used to transform
the problems into maximization ones.
Parameters
base.type
String, one of: rosenbrock rastrigin sphere step noisy-quartic kdj-f1 kdj-f2 kdj-f3 kdj-f4 booth griewank median sum product schwefel min rotated-rastrigin rotated-schwefel rotated-griewank langerman lennard-jones lunacek
(The vector problem to test against. Some of the types are synonyms: kdj-f1 = sphere, kdj-f2 = rosenbrock, kdj-f3 = step, kdj-f4 = noisy-quartic. "kdj" stands for "Ken DeJong", and the numbers are the problems in his test suite) base.seed
int > 0
(Random number seed for rotated problems)
*/
public class ECSuite extends Problem implements SimpleProblemForm
{
public static final String P_SEED = "seed";
public static final String P_WHICH_PROBLEM = "type";
public static final String V_ROSENBROCK = "rosenbrock";
public static final String V_RASTRIGIN = "rastrigin";
public static final String V_SPHERE = "sphere";
public static final String V_STEP = "step";
public static final String V_NOISY_QUARTIC = "noisy-quartic";
public static final String V_F1 = "kdj-f1";
public static final String V_F2 = "kdj-f2";
public static final String V_F3 = "kdj-f3";
public static final String V_F4 = "kdj-f4";
public static final String V_BOOTH = "booth";
public static final String V_GRIEWANGK = "griewangk";
public static final String V_GRIEWANK = "griewank";
public static final String V_MEDIAN = "median";
public static final String V_SUM = "sum";
public static final String V_PRODUCT = "product";
public static final String V_SCHWEFEL = "schwefel";
public static final String V_MIN = "min";
public static final String V_ROTATED_RASTRIGIN = "rotated-rastrigin";
public static final String V_ROTATED_SCHWEFEL = "rotated-schwefel";
public static final String V_ROTATED_GRIEWANK = "rotated-griewank";
public static final String V_LANGERMAN = "langerman" ;
public static final String V_LENNARDJONES = "lennard-jones" ;
public static final String V_LUNACEK = "lunacek" ;
public static final int PROB_ROSENBROCK = 0;
public static final int PROB_RASTRIGIN = 1;
public static final int PROB_SPHERE = 2;
public static final int PROB_STEP = 3;
public static final int PROB_NOISY_QUARTIC = 4;
public static final int PROB_BOOTH = 5;
public static final int PROB_GRIEWANK = 6;
public static final int PROB_MEDIAN = 7;
public static final int PROB_SUM = 8;
public static final int PROB_PRODUCT = 9;
public static final int PROB_SCHWEFEL = 10;
public static final int PROB_MIN = 11;
public static final int PROB_ROTATED_RASTRIGIN = 12;
public static final int PROB_ROTATED_SCHWEFEL = 13;
public static final int PROB_ROTATED_GRIEWANK = 14;
public static final int PROB_LANGERMAN = 15 ;
public static final int PROB_LENNARDJONES = 16 ;
public static final int PROB_LUNACEK = 17 ;
public int problemType = PROB_ROSENBROCK; // defaults on Rosenbrock
public static final String problemName[] = new String[]
{
V_ROSENBROCK,
V_RASTRIGIN,
V_SPHERE,
V_STEP,
V_NOISY_QUARTIC,
V_BOOTH,
V_GRIEWANK,
V_MEDIAN,
V_SUM,
V_PRODUCT,
V_SCHWEFEL,
V_MIN,
V_ROTATED_RASTRIGIN,
V_ROTATED_SCHWEFEL,
V_ROTATED_GRIEWANK,
V_LANGERMAN,
V_LENNARDJONES,
V_LUNACEK
};
public static final double minRange[] = new double[]
{
-2.048, // rosenbrock
-5.12, // rastrigin
-5.12, // sphere
-5.12, // step
-1.28, // noisy quartic
-5.12, // booth
-600.0, // griewank
0.0, // median
0.0, // sum
0.0, // product
-512.03, // schwefel
0.0, // min
-5.12, // rotated-rastrigin
-512.03, // rotated-schwefel
-600.0, // rotated-griewank
0, // langerman
-3.0, // lennard-jones
-5.0, // lunacek
};
public static final double maxRange[] = new double[]
{
2.048, // rosenbrock
5.12, // rastrigin
5.12, // sphere
5.12, // step
1.28, // noisy quartic
5.12, // booth
600.0, // griewank
1.0, // median
1.0, // sum
2.0, // product
511.97, // schwefel
1.0, // min
5.12, // rotated-rastrigin
511.97, // rotated-schwefel
600.0, // rotated-griewank
10, // langerman
3.0, // lennard-jones
5.0 // lunacek
};
public long seed; // rotation seed for rotation problems
boolean alreadyChecked = false;
public void checkRange(EvolutionState state, int problem, double[] genome)
{
if (alreadyChecked || state.generation > 0) return;
alreadyChecked = true;
for(int i = 0; i < state.population.subpops.length; i++)
{
if (!(state.population.subpops[i].species instanceof FloatVectorSpecies))
{
state.output.fatal("ECSuite requires species " + i + " to be a FloatVectorSpecies, but it is a: " + state.population.subpops[i].species);
}
FloatVectorSpecies species = (FloatVectorSpecies)(state.population.subpops[i].species);
for(int k = 0; k < genome.length; k++)
{
if (species.minGene(k) != minRange[problem] ||
species.maxGene(k) != maxRange[problem])
{
state.output.warning("Gene range is nonstandard for problem " + problemName[problem] + ".\nFirst occurrence: Subpopulation " + i + " Gene " + k +
" range was [" + species.minGene(k) + ", " + species.maxGene(k) +
"], expected [" + minRange[problem] + ", " + maxRange[problem] + "]");
return; // done here
}
}
}
if (problemType == PROB_LANGERMAN)
{
// Langerman has a maximum genome size of 10
if (genome.length > 10)
state.output.fatal("The Langerman function requires that the genome size be a value from 1 to 10 inclusive. It is presently " + genome.length);
}
else if (problemType == PROB_LENNARDJONES)
{
// Lennard-Jones requires that its genomes be multiples of 3
if (genome.length % 3 != 0)
state.output.fatal("The Lennard-Jones function requires that the genome size be a multiple of 3. It is presently " + genome.length);
}
}
// nothing....
public void setup(final EvolutionState state, final Parameter base)
{
super.setup(state, base);
String wp = state.parameters.getStringWithDefault( base.push( P_WHICH_PROBLEM ), null, "" );
if( wp.compareTo( V_ROSENBROCK ) == 0 || wp.compareTo (V_F2)==0 )
problemType = PROB_ROSENBROCK;
else if ( wp.compareTo( V_RASTRIGIN ) == 0 )
problemType = PROB_RASTRIGIN;
else if ( wp.compareTo( V_SPHERE ) == 0 || wp.compareTo (V_F1)==0)
problemType = PROB_SPHERE;
else if ( wp.compareTo( V_STEP ) == 0 || wp.compareTo (V_F3)==0)
problemType = PROB_STEP;
else if ( wp.compareTo( V_NOISY_QUARTIC ) == 0 || wp.compareTo (V_F4)==0)
problemType = PROB_NOISY_QUARTIC;
else if( wp.compareTo( V_BOOTH ) == 0 )
problemType = PROB_BOOTH;
else if( wp.compareTo( V_GRIEWANK ) == 0 )
problemType = PROB_GRIEWANK;
else if (wp.compareTo( V_GRIEWANGK ) == 0 )
{
state.output.warning("Incorrect parameter name (\"griewangk\") used, should be \"griewank\"", base.push( P_WHICH_PROBLEM ), null );
problemType = PROB_GRIEWANK;
}
else if( wp.compareTo( V_MEDIAN ) == 0 )
problemType = PROB_MEDIAN;
else if( wp.compareTo( V_SUM ) == 0 )
problemType = PROB_SUM;
else if( wp.compareTo( V_PRODUCT ) == 0 )
problemType = PROB_PRODUCT;
else if (wp.compareTo( V_SCHWEFEL ) == 0 )
problemType = PROB_SCHWEFEL;
else if (wp.compareTo( V_MIN ) == 0 )
problemType = PROB_MIN;
else if (wp.compareTo( V_ROTATED_RASTRIGIN) == 0)
problemType = PROB_ROTATED_RASTRIGIN;
else if (wp.compareTo( V_ROTATED_SCHWEFEL) == 0)
problemType = PROB_ROTATED_SCHWEFEL;
else if (wp.compareTo( V_ROTATED_GRIEWANK) == 0)
problemType = PROB_ROTATED_GRIEWANK;
else if (wp.compareTo(V_LANGERMAN) == 0)
problemType = PROB_LANGERMAN ;
else if (wp.compareTo(V_LENNARDJONES) == 0)
problemType = PROB_LENNARDJONES ;
else if (wp.compareTo(V_LUNACEK) == 0)
problemType = PROB_LUNACEK ;
else state.output.fatal(
"Invalid value for parameter, or parameter not found.\n" +
"Acceptable values are:\n" +
" " + V_ROSENBROCK + " (or " + V_F2 + ")\n" +
" " + V_RASTRIGIN + "\n" +
" " + V_SPHERE + " (or " + V_F1 + ")\n" +
" " + V_STEP + " (or " + V_F3 + ")\n" +
" " + V_NOISY_QUARTIC + " (or " + V_F4 + ")\n"+
" " + V_BOOTH + "\n" +
" " + V_GRIEWANK + "\n" +
" " + V_MEDIAN + "\n" +
" " + V_SUM + "\n" +
" " + V_PRODUCT + "\n" +
" " + V_SCHWEFEL + "\n"+
" " + V_MIN + "\n"+
" " + V_ROTATED_RASTRIGIN + "\n" +
" " + V_ROTATED_SCHWEFEL + "\n" +
" " + V_ROTATED_GRIEWANK + "\n" +
" " + V_LANGERMAN + "\n" +
" " + V_LENNARDJONES + "\n" +
" " + V_LUNACEK + "\n",
base.push( P_WHICH_PROBLEM ) );
seed = state.parameters.getLongWithDefault( base.push( P_SEED ), null, ROTATION_SEED );
if (seed <= 0)
state.output.fatal("If a rotation seed is provided, it must be > 0", base.push( P_SEED ), null);
}
public void evaluate(final EvolutionState state,
final Individual ind,
final int subpopulation,
final int threadnum)
{
if (ind.evaluated) // don't bother reevaluating
return;
if( !( ind instanceof DoubleVectorIndividual ) )
state.output.fatal( "The individuals for this problem should be DoubleVectorIndividuals." );
DoubleVectorIndividual temp = (DoubleVectorIndividual)ind;
double[] genome = temp.genome;
//int len = genome.length;
// this curious break-out makes it easy to use the isOptimal() and function() methods
// for other purposes, such as coevolutionary versions of this class.
// compute the fitness on a per-function basis
double fit = (function(state, problemType, temp.genome, threadnum));
// compute if we're optimal on a per-function basis
boolean isOptimal = isOptimal(problemType, fit);
// set the fitness appropriately
if (fit < (0.0 - Double.MAX_VALUE)) // uh oh -- can be caused by Product for example
{
((SimpleFitness)(ind.fitness)).setFitness( state, 0.0 - Double.MAX_VALUE, isOptimal );
state.output.warnOnce("'Product' type used: some fitnesses are negative infinity, setting to lowest legal negative number.");
}
else if (fit > Double.MAX_VALUE) // uh oh -- can be caused by Product for example
{
((SimpleFitness)(ind.fitness)).setFitness( state, Double.MAX_VALUE, isOptimal );
state.output.warnOnce("'Product' type used: some fitnesses are negative infinity, setting to lowest legal negative number.");
}
else
{
((SimpleFitness)(ind.fitness)).setFitness( state, fit, isOptimal );
}
ind.evaluated = true;
}
public boolean isOptimal(int function, double fitness)
{
switch(problemType)
{
case PROB_ROSENBROCK:
case PROB_RASTRIGIN:
case PROB_SPHERE:
case PROB_STEP:
return fitness == 0.0;
case PROB_NOISY_QUARTIC:
case PROB_BOOTH:
case PROB_GRIEWANK:
case PROB_MEDIAN:
case PROB_SUM:
case PROB_PRODUCT:
case PROB_SCHWEFEL:
case PROB_ROTATED_RASTRIGIN: // not sure
case PROB_ROTATED_SCHWEFEL:
case PROB_ROTATED_GRIEWANK:
case PROB_MIN:
case PROB_LANGERMAN: // may be around -1.4
case PROB_LENNARDJONES:
case PROB_LUNACEK:
default:
return false;
}
}
public double function(EvolutionState state, int function, double[] genome, int threadnum)
{
checkRange(state, function, genome);
double value = 0;
double len = genome.length;
switch(function)
{
case PROB_ROSENBROCK:
for( int i = 1 ; i < len ; i++ )
{
double gj = genome[i-1] ;
double gi = genome[i] ;
value += 100 * (gj*gj - gj) * (gj*gj - gj) + (1-gj) * (1-gj);
}
return -value;
case PROB_RASTRIGIN:
final double A = 10.0;
value = len * A;
for( int i = 0 ; i < len ; i++ )
{
double gi = genome[i] ;
value += ( gi*gi - A * Math.cos( 2 * Math.PI * gi ) );
}
return -value;
case PROB_SPHERE:
for( int i = 0 ; i < len ; i++ )
{
double gi = genome[i] ;
value += gi * gi;
}
return -value;
case PROB_STEP:
for( int i = 0 ; i < len ; i++ )
{
double gi = genome[i] ;
// The reason for the 6 is that this is the equation De Jong used in the De Jong Test Suite
value += 6 + Math.floor( gi );
}
return -value;
case PROB_NOISY_QUARTIC:
for( int i = 0 ; i < len ; i++ )
{
double gi = genome[i] ;
value += (i+1)*(gi*gi*gi*gi) + state.random[threadnum].nextDouble();
}
return -value;
case PROB_BOOTH:
if( len != 2 )
state.output.fatal( "The Booth problem is defined for only two terms, and as a consequence the genome of the DoubleVectorIndividual should have size 2." );
double g0 = genome[0] ;
double g1 = genome[1] ;
value = (g0 + 2*g1 - 7) * (g0 + 2*g1 - 7) +
(2*g0 + g1 - 5) * (2*g0 + g1 - 5);
return -value;
case PROB_GRIEWANK:
value = 1;
double prod = 1;
for( int i = 0 ; i < len ; i++ )
{
double gi = genome[i] ;
value += (gi*gi)/4000.0;
prod *= Math.cos( gi / Math.sqrt(i+1) );
}
value -= prod;
return -value;
case PROB_SCHWEFEL:
for( int i = 0 ; i < len ; i++ )
{
double gi = genome[i] ;
value += -gi * Math.sin(Math.sqrt(Math.abs(gi)));
}
return -value;
case PROB_MEDIAN: // FIXME, need to do a better median-finding algorithm, such as http://www.ics.uci.edu/~eppstein/161/960130.html
double[] sorted = new double[(int)len];
System.arraycopy(genome, 0, sorted, 0, sorted.length);
ec.util.QuickSort.qsort(sorted);
return sorted[sorted.length / 2] ; // note positive
case PROB_SUM:
value = 0.0;
for( int i = 0 ; i < len ; i++ )
{
double gi = genome[i] ;
value += gi;
}
return value; // note positive
case PROB_MIN:
value = genome[0] ;
for( int i = 1 ; i < len ; i++ )
{
double gi = genome[i] ;
if (value > gi) value = gi;
}
return value; // note positive
case PROB_PRODUCT:
value = 1.0;
for( int i = 0 ; i < len ; i++ )
{
double gi = genome[i] ;
value *= gi;
}
return value; // note positive
case PROB_ROTATED_RASTRIGIN:
{
synchronized(rotationMatrix) // synchronizations are rare in ECJ. :-(
{
if (rotationMatrix[0] == null)
rotationMatrix[0] = buildRotationMatrix(state, seed, (int)len);
}
// now we know the matrix exists rotate the matrix and return its value
double[] val = mul(rotationMatrix[0], genome);
return function(state, PROB_RASTRIGIN, val, threadnum);
}
case PROB_ROTATED_SCHWEFEL:
{
synchronized(rotationMatrix) // synchronizations are rare in ECJ. :-(
{
if (rotationMatrix[0] == null)
rotationMatrix[0] = buildRotationMatrix(state, seed, (int)len);
}
// now we know the matrix exists rotate the matrix and return its value
double[] val = mul(rotationMatrix[0], genome);
return function(state, PROB_SCHWEFEL, val, threadnum);
}
case PROB_ROTATED_GRIEWANK:
{
synchronized(rotationMatrix) // synchronizations are rare in ECJ. :-(
{
if (rotationMatrix[0] == null)
rotationMatrix[0] = buildRotationMatrix(state, seed, (int)len);
}
// now we know the matrix exists rotate the matrix and return its value
double[] val = mul(rotationMatrix[0], genome);
return function(state, PROB_GRIEWANK, val, threadnum);
}
case PROB_LANGERMAN:
{
return 0.0 - langerman(genome);
}
case PROB_LENNARDJONES:
{
int numAtoms = genome.length / 3;
double v = 0.0 ;
for(int i = 0 ; i < numAtoms - 1 ; i++ )
{
for(int j = i + 1 ; j < numAtoms ; j++ )
{
// double d = dist(genome, i, j);
double a = genome[i * 3] - genome[j * 3];
double b = genome[i * 3 + 1] - genome[j * 3 + 1];
double c = genome[i * 3 + 2] - genome[j * 3 + 2];
double d = Math.sqrt(a * a + b * b + c * c);
double r12 = Math.pow(d, -12.0);
double r6 = Math.pow(d, -6.0);
double e = r12 - r6 ;
v += e ;
}
}
v *= -4.0 ;
return v;
}
case PROB_LUNACEK:
{
// Lunacek function: for more information, please see --
// http://arxiv.org/pdf/1207.4318.pdf
// http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.154.1657
// // // //
double s = 0.7 ; // The shape of the boundary of the double sphere,
// could be like [0.2 - 1.4] but not 0.0.
// > 1.0 or < 1.0 means a parabolic shape, 1.0 means a linear boundary.
double d = 1.0 ; // depth of the sphere, could be 1, 2, 3, or 4. 1 is deeper than 4
// this could be also be a fraction I guess.
double mu1 = 0.0 ;
for(int i = 0 ; i < genome.length ; i++)
mu1 = genome[i] ;
double mu2 = -1.0 * Math.sqrt(((mu1 * mu1) - d)/s);
double sigma1 = 0.0;
double sigma2 = 0.0 ;
for(int i = 0 ; i < genome.length ; i++)
{
sigma1 = (genome[i] - mu1) * (genome[i] - mu1);
sigma2 = (genome[i] - mu2) * (genome[i] - mu2);
}
sigma2 = d * genome.length + s * sigma1 ;
double sphere = Math.min(sigma1, sigma2);
double rastrigin = function(state, PROB_RASTRIGIN, genome, threadnum);
// + or - ? not sure, I always get confused.
// Lunacek function is a combination of Rastrigin and a Double Sphere.
// As the Rastrigin is -, so this function should be.
return -1.0 * (sphere + rastrigin) ;
}
default:
state.output.fatal( "ec.app.ecsuite.ECSuite has an invalid problem -- how on earth did that happen?" );
return 0; // never happens
}
}
// magic arrays for the Langerman problem
private double[][] afox10 =
{
{9.681, 0.667, 4.783, 9.095, 3.517, 9.325, 6.544, 0.211, 5.122, 2.020},
{9.400, 2.041, 3.788, 7.931, 2.882, 2.672, 3.568, 1.284, 7.033, 7.374},
{8.025, 9.152, 5.114, 7.621, 4.564, 4.711, 2.996, 6.126, 0.734, 4.982},
{2.196, 0.415, 5.649, 6.979, 9.510, 9.166, 6.304, 6.054, 9.377, 1.426},
{8.074, 8.777, 3.467, 1.863, 6.708, 6.349, 4.534, 0.276, 7.633, 1.567},
{7.650, 5.658, 0.720, 2.764, 3.278, 5.283, 7.474, 6.274, 1.409, 8.208},
{1.256, 3.605, 8.623, 6.905, 0.584, 8.133, 6.071, 6.888, 4.187, 5.448},
{8.314, 2.261, 4.224, 1.781, 4.124, 0.932, 8.129, 8.658, 1.208, 5.762},
{0.226, 8.858, 1.420, 0.945, 1.622, 4.698, 6.228, 9.096, 0.972, 7.637},
{7.305, 2.228, 1.242, 5.928, 9.133, 1.826, 4.060, 5.204, 8.713, 8.247},
{0.652, 7.027, 0.508, 4.876, 8.807, 4.632, 5.808, 6.937, 3.291, 7.016},
{2.699, 3.516, 5.874, 4.119, 4.461, 7.496, 8.817, 0.690, 6.593, 9.789},
{8.327, 3.897, 2.017, 9.570, 9.825, 1.150, 1.395, 3.885, 6.354, 0.109},
{2.132, 7.006, 7.136, 2.641, 1.882, 5.943, 7.273, 7.691, 2.880, 0.564},
{4.707, 5.579, 4.080, 0.581, 9.698, 8.542, 8.077, 8.515, 9.231, 4.670},
{8.304, 7.559, 8.567, 0.322, 7.128, 8.392, 1.472, 8.524, 2.277, 7.826},
{8.632, 4.409, 4.832, 5.768, 7.050, 6.715, 1.711, 4.323, 4.405, 4.591},
{4.887, 9.112, 0.170, 8.967, 9.693, 9.867, 7.508, 7.770, 8.382, 6.740},
{2.440, 6.686, 4.299, 1.007, 7.008, 1.427, 9.398, 8.480, 9.950, 1.675},
{6.306, 8.583, 6.084, 1.138, 4.350, 3.134, 7.853, 6.061, 7.457, 2.258},
{0.652, 0.343, 1.370, 0.821, 1.310, 1.063, 0.689, 8.819, 8.833, 9.070},
{5.558, 1.272, 5.756, 9.857, 2.279, 2.764, 1.284, 1.677, 1.244, 1.234},
{3.352, 7.549, 9.817, 9.437, 8.687, 4.167, 2.570, 6.540, 0.228, 0.027},
{8.798, 0.880, 2.370, 0.168, 1.701, 3.680, 1.231, 2.390, 2.499, 0.064},
{1.460, 8.057, 1.336, 7.217, 7.914, 3.615, 9.981, 9.198, 5.292, 1.224},
{0.432, 8.645, 8.774, 0.249, 8.081, 7.461, 4.416, 0.652, 4.002, 4.644},
{0.679, 2.800, 5.523, 3.049, 2.968, 7.225, 6.730, 4.199, 9.614, 9.229},
{4.263, 1.074, 7.286, 5.599, 8.291, 5.200, 9.214, 8.272, 4.398, 4.506},
{9.496, 4.830, 3.150, 8.270, 5.079, 1.231, 5.731, 9.494, 1.883, 9.732},
{4.138, 2.562, 2.532, 9.661, 5.611, 5.500, 6.886, 2.341, 9.699, 6.500}
};
private double[] cfox10 =
{
0.806, 0.517, 1.5, 0.908, 0.965,
0.669, 0.524, 0.902, 0.531, 0.876,
0.462, 0.491, 0.463, 0.714, 0.352,
0.869, 0.813, 0.811, 0.828, 0.964,
0.789, 0.360, 0.369, 0.992, 0.332,
0.817, 0.632, 0.883, 0.608, 0.326
};
private double langerman(double genome[])
{
double sum = 0 ;
for ( int i = 0 ; i < 30 ; i++ )
{
// compute squared distance
double distsq = 0.0;
double t;
double[] afox10i = afox10[i];
for(int j = 0; j < genome.length; j++)
{
t = genome[j] - afox10i[j];
distsq += t * t;
}
sum += cfox10[i] * Math.exp(-distsq / Math.PI) * Math.cos(distsq * Math.PI);
}
return 0 - sum;
}
/*
-----------------
Rotation facility
-----------------
This code is just used by the Rotated Schwefel and Rotated Rastrigin functions to rotate their
functions by a certain amount. The code is largely based on the rotation scheme described in
"Completely Derandomized Self-Adaptation in Evolutionary Strategies",
Nikolaus Hansen and Andreas Ostermeier, Evolutionary Computation 9(2): 159--195.
We fix a hard-coded rotation matrix which is the same for all problems, in order to guarantee
correctness in gathering results over multiple jobs. But you can change that easily if you like.
*/
public static double[][][] rotationMatrix = new double[1][][]; // the actual matrix is stored in rotationMatrix[0] -- a hack
/** Dot product between two column vectors. Does not modify the original vectors. */
public static double dot(double[] x, double[] y)
{
double val = 0;
for(int i =0; i < x.length; i++)
val += x[i] * y[i];
return val;
}
/** Multiply a column vector against a matrix[row][column]. Does not modify the original vector or matrix. */
public static double[] mul(double [/* row */ ][ /* column */] matrix, double[] x)
{
double[] val = new double[matrix.length];
for(int i = 0; i < matrix.length; i++)
{
double sum = 0.0;
double[] m = matrix[i];
for(int j = 0; j < m.length; j++)
sum += m[j] * x[j];
val[i] = sum;
}
return val;
}
/** Scalar multiply against a column vector. Does not modify the original vector. */
public static double[] scalarMul(double scalar, double[] x)
{
double[] val = new double[x.length];
for(int i =0; i < x.length; i++)
val[i] = x[i] * scalar;
return val;
}
/** Subtract two column vectors. Does not modify the original vectors. */
public static double[] sub(double[] x, double[] y)
{
double[] val = new double[x.length];
for(int i =0; i < x.length; i++)
val[i] = x[i] - y[i];
return val;
}
/** Normalize a column vector. Does not modify the original vector. */
public static double[] normalize(double[] x)
{
double[] val = new double[x.length];
double sumsq = 0;
for(int i =0; i < x.length; i++)
sumsq += x[i] * x[i];
sumsq = Math.sqrt(sumsq);
for(int i =0; i < x.length; i++)
val[i] = x[i] / sumsq;
return val;
}
/** Fixed rotation seed so all jobs use exactly the same rotation space. */
public static final long ROTATION_SEED = 9731297;
/** Build an NxN rotation matrix[row][column] with a given seed. */
public static double[ /* row */ ][ /* column */] buildRotationMatrix(EvolutionState state, long rotationSeed, int N)
{
if (rotationSeed == ROTATION_SEED)
state.output.warnOnce("Default rotation seed being used (" + rotationSeed + ")");
MersenneTwisterFast rand = new MersenneTwisterFast(rotationSeed); // it's rare to need to do this, but we need to guarantee the same rotation space
for(int i = 0; i < 624 * 4; i++) // prime the MT for 4 full sample iterations to get it warmed up
rand.nextInt();
double o[ /* row */ ][ /* column */ ] = new double[N][N];
// make random values
for(int i = 0; i < N; i++)
for(int k = 0; k < N; k++)
o[i][k] = rand.nextGaussian();
// build random values
for(int i = 0; i < N; i++)
{
// extract o[i] -> no
double[] no = new double[N];
for(int k=0; k < N; k++)
no[k] = o[i][k];
// go through o[i] and o[j], modifying no
for(int j = 0; j < i; j++)
{
double d = dot(o[i], o[j]);
double[] val = scalarMul(d, o[j]);
no = sub(no, val);
}
o[i] = normalize(no);
}
return o;
}
}
Ads help maintain this website.