1 module minimized; 2 3 import std.algorithm : swapAt; 4 import std.math : isNaN, approxEqual; 5 import std.range : front, popFront, isRandomAccessRange; 6 import std.random : uniform, rndGen, isUniformRNG; 7 import std.traits : isCallable; 8 import std.functional : toDelegate; 9 import std.exception : enforce; 10 11 version( unittest ) 12 { 13 import std.stdio : writeln; 14 import std.algorithm : equal; 15 } 16 17 private void partialShuffle(Range, RandomGen)(Range r, in size_t n, ref RandomGen gen) 18 if(isRandomAccessRange!Range && isUniformRNG!RandomGen) 19 { 20 enforce(n <= r.length, "n must be <= r.length for partialShuffle."); 21 foreach (i; 0 .. n) 22 { 23 swapAt(r, i, i + uniform(0, r.length - i, gen)); 24 } 25 } 26 27 private void partialShuffle(Range)(Range r, in size_t n) 28 if(isRandomAccessRange!Range) 29 { 30 return partialShuffle(r, n, rndGen); 31 } 32 33 private struct Individual( RANGE ) 34 { 35 double temperature; 36 RANGE parameters; 37 38 /// Holds f and cr control parameters; 39 double[] evolutionControl; 40 } 41 42 class DifferentialEvolution( RANGE ) 43 { 44 double tau1 = 0.1; /// Tweaking parameter, see Brest et al for details 45 double tau2 = 0.1; /// Tweaking parameter, see Brest et al for details 46 47 @property void temperatureFunction(F)(auto ref F fp) 48 if (isCallable!F) 49 { 50 _temperatureFunction = toDelegate!F(fp); 51 } 52 53 double temperatureFunction( RANGE parameters ) 54 { 55 return _temperatureFunction( parameters ); 56 } 57 58 @property void randomIndividual(F)(auto ref F fp) 59 if (isCallable!F) 60 { 61 _randomIndividual = toDelegate!F(fp); 62 } 63 64 RANGE randomIndividual() 65 { 66 return _randomIndividual(); 67 } 68 69 void initialize() 70 { 71 foreach( i; 0..(10*randomIndividual().length) ) 72 { 73 auto pars = randomIndividual(); 74 auto individual = Individual!RANGE( temperatureFunction( pars ), 75 pars ); 76 77 individual.evolutionControl = [uniform( 0.1, 1.2 ), 78 uniform( 0.0, 1.0 )]; 79 80 population ~= individual; 81 82 if ( bestFit.temperature.isNaN || bestFit.temperature >= individual.temperature ) 83 bestFit = individual; 84 } 85 86 } 87 88 /// Perform one minimization step (generation) 89 /// 90 /// Return true if improvement was made 91 bool step() 92 { 93 if (population.length == 0) 94 initialize; 95 Individual!RANGE[] processed; 96 Individual!RANGE[] toProcess = population; 97 bool anyAccepted = false; 98 foreach( _; 0..toProcess.length ) 99 { 100 auto x = toProcess.front; 101 toProcess.popFront; 102 103 auto abc = (processed ~ toProcess); 104 105 abc.partialShuffle( 3 ); 106 107 size_t chosenR = uniform( 0, x.parameters.length-1 ); 108 109 Individual!RANGE y; 110 if (uniform(0.0,1.0)<tau1) 111 y.evolutionControl ~= uniform(0.1,1.2); 112 else 113 y.evolutionControl ~= x.evolutionControl[0]; 114 if (uniform(0.0,1.0)<tau2) 115 y.evolutionControl ~= uniform(0.0,1.0); 116 else 117 y.evolutionControl ~= x.evolutionControl[1]; 118 119 120 121 foreach( i; 0..(x.parameters.length ) ) 122 { 123 if ( i == chosenR 124 || uniform(0.0,1.0) < y.evolutionControl[1] ) 125 y.parameters ~= abc[0].parameters[i] 126 + y.evolutionControl[0] 127 *( abc[1].parameters[i] 128 - abc[2].parameters[i] ); 129 else 130 y.parameters ~= x.parameters[i]; 131 } 132 y.temperature = temperatureFunction( y.parameters ); 133 134 if ( y.temperature < x.temperature ) 135 { 136 anyAccepted = true; 137 processed ~= y; 138 if (y.temperature < bestFit.temperature) 139 bestFit = y; 140 } else { 141 processed ~= x; 142 } 143 } 144 145 if (anyAccepted) 146 { 147 population = processed; 148 } 149 return anyAccepted; 150 } 151 152 /// Return the current best fitting parameters 153 RANGE currentBestFit() 154 { 155 return bestFit.parameters; 156 } 157 158 RANGE minimize() 159 { 160 161 size_t sinceAccepted = 0; 162 163 while( sinceAccepted < 10 ) 164 { 165 if (step()) 166 { 167 sinceAccepted = 0; 168 } else { 169 sinceAccepted++; 170 } 171 } 172 173 return bestFit.parameters; 174 } 175 176 private: 177 Individual!RANGE[] population; 178 179 double delegate( RANGE parameters ) _temperatureFunction; 180 RANGE delegate() _randomIndividual; 181 Individual!RANGE bestFit; /// Current bestFit found 182 } 183 184 /// 185 unittest 186 { 187 // Function to minimize 188 auto fn = ( double[] xs ) { 189 auto p = [ 1.0, 2.0, 10, 20, 30 ]; 190 return p[2] * (xs[0] - p[0]) * (xs[0] - p[0]) + 191 p[3] * (xs[1] - p[1]) * (xs[1] - p[1]) + p[4]; 192 }; 193 194 // Function which will create random initial sets of parameters 195 auto initFunction = () 196 { 197 return [ uniform( 0.0, 10.0 ), uniform( 0.0, 10.0 )]; 198 }; 199 200 auto de = new DifferentialEvolution!(double[])(); 201 de.temperatureFunction = fn; 202 de.randomIndividual = initFunction; 203 204 auto min = de.minimize; 205 206 assert( equal!approxEqual( min, [ 1, 2 ] ) ); 207 } 208 209 unittest 210 { 211 // Function to minimize 212 auto p = [ 1.0, 2.0, 10, 20, 30 ]; 213 auto fn = ( double[] xs ) { 214 return p[2] * (xs[0] - p[0]) * (xs[0] - p[0]) + 215 p[3] * (xs[1] - p[1]) * (xs[1] - p[1]) + p[4]; 216 }; 217 218 // Function which will create random initial sets of parameters 219 auto initFunction = delegate double[]() 220 { 221 return [ uniform( 0.0, 10.0 ), uniform( 0.0, 10.0 )]; 222 }; 223 224 auto de = new DifferentialEvolution!(double[])(); 225 de.temperatureFunction = fn; 226 de.randomIndividual = initFunction; 227 228 auto min = de.minimize; 229 230 assert( equal!approxEqual( min, [ 1, 2 ] ) ); 231 assert( equal( min, de.currentBestFit ) ); 232 }