diff --git a/eo/src/es/eoNormalMutation.h b/eo/src/es/eoNormalMutation.h index cc3276e8..922ee6aa 100644 --- a/eo/src/es/eoNormalMutation.h +++ b/eo/src/es/eoNormalMutation.h @@ -33,12 +33,16 @@ #include #include #include - +#include //----------------------------------------------------------------------------- /** Simple normal mutation of a vector of real values. * The stDev is fixed - but it is passed ans stored as a reference, * to enable dynamic mutations (see eoOenFithMutation below). + * + * As for the bounds, the values are here folded back into the bounds. + * The other possiblity would be to iterate until we fall inside the bounds - + * but this sometimes takes a long time!!! */ template class eoNormalMutation: public eoMonOp @@ -46,11 +50,23 @@ template class eoNormalMutation: public eoMonOp public: /** * (Default) Constructor. + * The bounds are initialized with the global object that says: no bounds. + * * @param _sigma the range for uniform nutation * @param _p_change the probability to change a given coordinate */ eoNormalMutation(double & _sigma, const double& _p_change = 1.0): - sigma(_sigma), p_change(_p_change) {} + sigma(_sigma), bounds(eoDummyVectorNoBounds), p_change(_p_change) {} + + /** + * Constructor with bounds + * @param _bounds an eoRealVectorBounds that contains the bounds + * @param _sigma the range for uniform nutation + * @param _p_change the probability to change a given coordinate + */ + eoNormalMutation(eoRealVectorBounds & _bounds, + double & _sigma, const double& _p_change = 1.0): + sigma(_sigma), bounds(_bounds), p_change(_p_change) {} /// The class name. string className() const { return "eoNormalMutation"; } @@ -67,6 +83,7 @@ template class eoNormalMutation: public eoMonOp if (rng.flip(p_change)) { _eo[lieu] += sigma*rng.normal(); + bounds.foldsInBounds(lieu, _eo[lieu]); hasChanged = true; } } @@ -77,6 +94,7 @@ template class eoNormalMutation: public eoMonOp protected: double & sigma; private: + eoRealVectorBounds & bounds; double p_change; }; diff --git a/eo/src/es/eoRealBounds.h b/eo/src/es/eoRealBounds.h index 1ddf9983..97b5f551 100644 --- a/eo/src/es/eoRealBounds.h +++ b/eo/src/es/eoRealBounds.h @@ -30,11 +30,6 @@ /** \defgroup EvolutionStrategies - Various classes for the initialization and mutation of real valued vectors. - - Supports simple mutations and various more adaptable mutations, including - correlated mutations. - */ @@ -42,33 +37,113 @@ \class eoRealBounds eoRealBounds.h es/eoRealBounds.h \ingroup EvolutionStrategies - Defines the minima and maxima for real variables + Defines bound classes for real numbers. + +Scalar type: +------------ +Basic class is eoRealBounds, a pure virtual. + +The following pure virtual methods are to be used in mutations: +- void foldsInBounds(double &) that folds any value that falls out of + the bounds back into the bounds, by bouncing on the limit (if any) +- bool isInBounds(double &) that simply says whether or not the argument + is in the bounds + +So mutation can choose whetehr they want to iterate trying until +they fall in bounds, or only try once and "repair" by using +the foldsInBounds method + +There is also a uniform() method that generates a uniform value +(if possible, i.e. if bounded) in the interval. + +Derived class are + eoRealInterval, that holds a minimum and maximum value + eoRealNoBounds, that implements the "unbounded bounds" + +TODO: the eoRealMinBound and eoRealMaxBound that implement + the half-bounded intervals. + +Vector type: +------------ +Class eoRealVectorBounds implements the vectorized version: +it is basically a vector of eoRealBounds * and forwards all request +to the elements of the vector. */ -class eoBaseRealBounds : public eoUF -{ }; +class eoRealBounds +{ +public: + virtual bool isBounded(void) = 0; + virtual bool isMinBounded(void) = 0; + virtual bool isMaxBounded(void) = 0; + virtual void foldsInBounds(double &) = 0; + virtual bool isInBounds(double) = 0; -class eoRealBounds : public eoBaseRealBounds + // accessors + virtual double minimum() = 0; + virtual double maximum() = 0; + virtual double range() = 0; + + // random generators + virtual double uniform(eoRng & _rng = eo::rng) = 0; +}; + +class eoRealNoBounds : public eoRealBounds +{ +public: + virtual ~eoRealNoBounds(){} + + virtual bool isBounded(void) {return false;} + virtual bool isMinBounded(void) {return false;} + virtual bool isMaxBounded(void) {return false;} + virtual void foldsInBounds(double &) {return;} + virtual bool isInBounds(double) {return true;} + + // accessors + virtual double minimum() + { + throw logic_error("Trying to get minimum of unbounded eoRealBounds"); + } + virtual double maximum() + { + throw logic_error("Trying to get maximum of unbounded eoRealBounds"); + } + virtual double range() + { + throw logic_error("Trying to get range of unbounded eoRealBounds"); + } + + // random generators + virtual double uniform(eoRng & _rng = eo::rng) + { + throw logic_error("Trying to generate uniform values in unbounded eoRealBounds"); + } +}; + +/* fully bounded == interval */ +class eoRealInterval : public eoRealBounds { public : /** Simple bounds = minimum and maximum (allowed) */ - eoRealBounds(double _min=0, double _max=1) : + eoRealInterval(double _min=0, double _max=1) : repMinimum(_min), repMaximum(_max), repRange(_max-_min) { if (repRange<=0) throw std::logic_error("Void range in eoRealBounds"); } - - double Minimum() { return repMinimum; } - double Maximum() { return repMaximum; } - double Range() { return repRange; } - // for backward compatibility - double minimum() { return repMinimum; } - double maximum() { return repMaximum; } - double range() { return repRange; } + + // accessors + virtual double minimum() { return repMinimum; } + virtual double maximum() { return repMaximum; } + virtual double range() { return repRange; } + + // description + virtual bool isBounded(void) {return true;} + virtual bool isMinBounded(void) {return true;} + virtual bool isMaxBounded(void) {return true;} double uniform(eoRng & _rng = eo::rng) { @@ -76,7 +151,7 @@ public : } // says if a given double is within the bounds - bool operator()(double _r) + virtual bool isInBounds(double _r) { if (_r < repMinimum) return false; @@ -85,6 +160,35 @@ public : return true; } + // folds a value into bounds + void foldsInBounds(double & _r) + { + long iloc; + double dlargloc = 2 * range() ; + + if (fabs(_r) > 1.0E9) // iloc too large! + { + _r = uniform(); + return; + } + + if ( (_r > maximum()) ) + { + iloc = (long) ( (_r-minimum()) / dlargloc ) ; + _r -= dlargloc * iloc ; + if ( _r > maximum() ) + _r = 2*maximum() - _r ; + } + + if (_r < minimum()) + { + iloc = (long) ( (maximum()-_r) / dlargloc ) ; + _r += dlargloc * iloc ; + if (_r < minimum()) + _r = 2*minimum() - _r ; + } + } + private : double repMinimum; double repMaximum; @@ -94,16 +198,18 @@ private : // now the vectorized version -class eoRealVectorBounds -{ -public : - +class eoRealVectorBounds : public vector +{ +public: + // virtual desctructor (to avoid warining?) + virtual ~eoRealVectorBounds(){} + /** Simple bounds = minimum and maximum (allowed) */ // Ctor: same bonds for everybody, explicit eoRealVectorBounds(unsigned _dim, double _min=0, double _max=1) : - vecMinimum(_dim, _min), vecMaximum(_dim, _max), vecRange(_dim, _max-_min) + vector(_dim, new eoRealInterval(_min, _max)) { if (_max-_min<=0) throw std::logic_error("Void range in eoRealVectorBounds"); @@ -111,114 +217,165 @@ public : // Ctor: same bonds for everybody, given as a eoRealBounds eoRealVectorBounds(unsigned _dim, eoRealBounds & _bounds) : - vecMinimum(_dim, _bounds.Minimum()), - vecMaximum(_dim, _bounds.Maximum()), - vecRange(_dim, _bounds.Range()) + vector(_dim, &_bounds) {} // Ctor: different bonds for different variables, vectors of double - eoRealVectorBounds(vector _min, vector _max) : - vecMinimum(_min), vecMaximum(_max), vecRange(_min.size()) + eoRealVectorBounds(vector _min, vector _max) { if (_max.size() != _min.size()) throw std::logic_error("Dimensions don't match in eoRealVectorBounds"); for (unsigned i=0; i<_min.size(); i++) { - vecRange[i]=_max[i]-_min[i]; - if (vecRange[i]<=0) - throw std::logic_error("Void range in eoRealVectorBounds"); + push_back( new eoRealInterval(_min[i], _max[i])); } } // Ctor, particular case of dim-2 eoRealVectorBounds(eoRealBounds & _xbounds, eoRealBounds & _ybounds) : - vecMinimum(2), vecMaximum(2), vecRange(2) + vector(0) { - vecMinimum[0] = _xbounds.Minimum(); - vecMaximum[0] = _xbounds.Maximum(); - vecRange[0] = _xbounds.Range(); - vecMinimum[1] = _ybounds.Minimum(); - vecMaximum[1] = _ybounds.Maximum(); - vecRange[1] = _ybounds.Range(); + push_back( &_xbounds); + push_back( &_ybounds); } - // not a ctor, but usefull to initialize, too - // is it safe to call it push_back? Maybe not, but it's meaningful! - void push_back(double _min=0, double _max=1) + virtual bool isBounded(unsigned _i) + { + return (*this)[_i]->isBounded(); + } + + // bounded iff all are bounded + virtual bool isBounded(void) { - vecMinimum.push_back(_min); - vecMaximum.push_back(_max); - if (_max-_min <= 0) - throw std::logic_error("Void range in eoRealVectorBounds::add"); - vecRange.push_back(_max-_min); + for (unsigned i=0; iisBounded()) + return false; + return true; } - void push_back(eoRealBounds & _bounds) + // these do not make any sense as vectors! + virtual bool isMinBounded(unsigned _i) + { return (*this)[_i]->isMinBounded();} ; + + virtual bool isMaxBounded(unsigned _i) + { return (*this)[_i]->isMaxBounded();} ; + + virtual void foldsInBounds(unsigned _i, double & _r) { - vecMinimum.push_back(_bounds.Minimum()); - vecMaximum.push_back(_bounds.Maximum()); - vecRange.push_back(_bounds.Range()); + (*this)[_i]->foldsInBounds(_r); } - // accessors - following rule that says that method start are capitalized - double Minimum(unsigned _i) { return vecMinimum[_i]; } - double Maximum(unsigned _i) { return vecMaximum[_i]; } - double Range(unsigned _i) { return vecRange[_i]; } - - // accessors - for backward compatibility - double minimum(unsigned _i) { return vecMinimum[_i]; } - double maximum(unsigned _i) { return vecMaximum[_i]; } - double range(unsigned _i) { return vecRange[_i]; } - - // handy: get the size - unsigned int size() { return vecMinimum.size();} - - // returns a value uniformly chosen in bounds for a given variable - double uniform(unsigned _i, eoRng & _rng = eo::rng) + virtual void foldsInBounds(vector & _v) { - return vecMinimum[_i] + _rng.uniform(vecRange[_i]); - } + for (unsigned i=0; i uniform(eoRng & _rng = eo::rng) + virtual bool isInBounds(unsigned _i, double _r) + { return (*this)[_i]->isInBounds(_r); } + + // isInBounds iff all are in bouds + virtual bool isInBounds(vector _v) { - vector v(vecMinimum.size()); - for (unsigned i=0; iminimum();} + virtual double maximum(unsigned _i) {return (*this)[_i]->maximum();} + virtual double range(unsigned _i) {return (*this)[_i]->range();} + + virtual double averageRange() + { + double r=0.0; + for (unsigned i=0; iuniform(); + return r; + } // fills a vector with uniformly chosen variables in bounds void uniform(vector & _v, eoRng & _rng = eo::rng) { - _v.resize(vecMinimum.size()); - for (unsigned i=0; i vecMaximum[_i]) - return false; - return true; - } - - // check the bounds for a vector: true only if ALL ar ein bounds - bool operator()(vector & _v) - { - for (unsigned i=0; i<_v.size(); i++) - if (! operator()(i, _v[i]) ) // out of bound - return false; - return true; - } -private : - vector vecMinimum; - vector vecMaximum; - vector vecRange; // to minimize operations ??? }; +// the dummy unbounded eoRealVectorBounds: + +class eoRealVectorNoBounds: public eoRealVectorBounds +{ +public: + // virtual desctructor (to avoid warining?) + virtual ~eoRealVectorNoBounds(){} + + /** + Simple bounds = minimum and maximum (allowed) + */ + // Ctor: nothing to do! + eoRealVectorNoBounds(unsigned _dim=0) : eoRealVectorBounds(_dim) {} + + + virtual bool isBounded(unsigned) {return false;} + virtual bool isBounded(void) {return false;} + virtual bool isMinBounded(unsigned) {return false;} + virtual bool isMaxBounded(unsigned) {return false;} + + virtual void foldsInBounds(unsigned, double &) {return;} + virtual void foldsInBounds(vector &) {return;} + + virtual bool isInBounds(unsigned, double) {return true;} + virtual bool isInBounds(vector) {return true;} + + // accessors + virtual double minimum(unsigned) + { + throw logic_error("Trying to get minimum of unbounded eoRealBounds"); + } + virtual double maximum(unsigned) + { + throw logic_error("Trying to get maximum of unbounded eoRealBounds"); + } + virtual double range(unsigned) + { + throw logic_error("Trying to get range of unbounded eoRealBounds"); + } + + virtual double averageRange() + { + throw logic_error("Trying to get average range of unbounded eoRealBounds"); + } + + // random generators + virtual double uniform(unsigned, eoRng & _rng = eo::rng) + { + throw logic_error("No uniform distribution on unbounded eoRealBounds"); + } + + // fills a vector with uniformly chosen variables in bounds + void uniform(vector &, eoRng & _rng = eo::rng) + { + throw logic_error("No uniform distribution on unbounded eoRealBounds"); + } + +}; + +// one object for all +eoRealNoBounds eoDummyRealNoBounds; +eoRealVectorNoBounds eoDummyVectorNoBounds; #endif diff --git a/eo/src/es/eoRealOp.h b/eo/src/es/eoRealOp.h index 1cdef1cd..586453a5 100644 --- a/eo/src/es/eoRealOp.h +++ b/eo/src/es/eoRealOp.h @@ -31,6 +31,7 @@ #include // swap_ranges #include #include +#include //----------------------------------------------------------------------------- @@ -46,11 +47,23 @@ template class eoUniformMutation: public eoMonOp public: /** * (Default) Constructor. + * The bounds are initialized with the global object that says: no bounds. + * * @param _epsilon the range for uniform nutation * @param _p_change the probability to change a given coordinate */ eoUniformMutation(const double& _epsilon, const double& _p_change = 1.0): - epsilon(_epsilon), p_change(_p_change) {} + bounds(eoDummyVectorNoBounds), epsilon(_epsilon), p_change(_p_change) {} + + /** + * Constructor with bounds + * @param _bounds an eoRealVectorBounds that contains the bounds + * @param _epsilon the range for uniform nutation + * @param _p_change the probability to change a given coordinate + */ + eoUniformMutation(eoRealVectorBounds & _bounds, + const double& _epsilon, const double& _p_change = 1.0): + bounds(_bounds), epsilon(_epsilon), p_change(_p_change) {} /// The class name. string className() const { return "eoUniformMutation"; } @@ -66,7 +79,14 @@ template class eoUniformMutation: public eoMonOp { if (rng.flip(p_change)) { - _eo[lieu] += 2*epsilon*rng.uniform()-epsilon; + // check the bounds + double emin = _eo[lieu]-epsilon; + double emax = _eo[lieu]+epsilon; + if (bounds.isMinBounded(lieu)) + emin = max(bounds.minimum(lieu), emin); + if (bounds.isMaxBounded(lieu)) + emax = min(bounds.maximum(lieu), emax); + _eo[lieu] = emin + (emax-emin)*rng.uniform(); hasChanged = true; } } @@ -75,6 +95,7 @@ template class eoUniformMutation: public eoMonOp } private: + eoRealVectorBounds & bounds; double epsilon; double p_change; }; @@ -94,7 +115,17 @@ template class eoDetUniformMutation: public eoMonOp * @param number of coordinate to modify */ eoDetUniformMutation(const double& _epsilon, const unsigned& _no = 1): - epsilon(_epsilon), no(_no) {} + bounds(eoDummyVectorNoBounds), epsilon(_epsilon), no(_no) {} + + /** + * Constructor with bounds + * @param _bounds an eoRealVectorBounds that contains the bounds + * @param _epsilon the range for uniform nutation + * @param number of coordinate to modify + */ + eoDetUniformMutation(eoRealVectorBounds & _bounds, + const double& _epsilon, const unsigned& _no = 1): + bounds(_bounds), epsilon(_epsilon), no(_no) {} /// The class name. string className() const { return "eoDetUniformMutation"; } @@ -110,11 +141,20 @@ template class eoDetUniformMutation: public eoMonOp { unsigned lieu = rng.random(_eo.size()); // actually, we should test that we don't re-modify same variable! - _eo[lieu] += 2*epsilon*rng.uniform()-epsilon; + + // check the bounds + double emin = _eo[lieu]-epsilon; + double emax = _eo[lieu]+epsilon; + if (bounds.isMinBounded(lieu)) + emin = max(bounds.minimum(lieu), emin); + if (bounds.isMaxBounded(lieu)) + emax = min(bounds.maximum(lieu), emax); + _eo[lieu] = emin + (emax-emin)*rng.uniform(); } } private: + eoRealVectorBounds & bounds; double epsilon; unsigned no; }; @@ -133,12 +173,27 @@ template class eoSegmentCrossover: public eoQuadraticOp public: /** * (Default) Constructor. - * @param _alpha the amount of exploration OUTSIDE the parents + * The bounds are initialized with the global object that says: no bounds. + * + * @param _alphaMin the amount of exploration OUTSIDE the parents * as in BLX-alpha notation (Eshelman and Schaffer) * 0 == contractive application + * Must be positive */ eoSegmentCrossover(const double& _alpha = 0.0) : - alpha(_alpha), range(1+2*alpha) {} + bounds(eoDummyVectorNoBounds), alpha(_alpha), range(1+2*_alpha) {} + + /** + * Constructor with bounds + * @param _bounds an eoRealVectorBounds that contains the bounds + * @param _alphaMin the amount of exploration OUTSIDE the parents + * as in BLX-alpha notation (Eshelman and Schaffer) + * 0 == contractive application + * Must be positive + */ + eoSegmentCrossover(eoRealVectorBounds & _bounds, + const double& _alpha = 0.0) : + bounds(_bounds), alpha(_alpha), range(1+2*_alpha) {} /// The class name. string className() const { return "eoSegmentCrossover"; } @@ -152,7 +207,35 @@ template class eoSegmentCrossover: public eoQuadraticOp { unsigned i; double r1, r2, fact; - fact = rng.uniform(range); // in [0,range) + double alphaMin = -alpha; + double alphaMax = 1+alpha; + if (alpha == 0.0) // no check to perform + fact = -alpha + rng.uniform(range); // in [-alpha,1+alpha) + else // look for the bounds for fact + { + for (i=0; i<_eo1.size(); i++) + { + r1=_eo1[i]; + r2=_eo2[i]; + if (r1 != r2) { // otherwise you'll get NAN's + double rmin = min(r1, r2); + double rmax = max(r1, r2); + double length = rmax - rmin; + if (bounds.isMinBounded(i)) + { + alphaMin = max(alphaMin, (bounds.minimum(i)-rmin)/length); + alphaMin = max(alphaMin, (rmax-bounds.maximum(i))/length); + } + if (bounds.isMaxBounded(i)) + { + alphaMax = min(alphaMax, (bounds.maximum(i)-rmin)/length); + alphaMax = min(alphaMax, (rmax-bounds.minimum(i))/length); + } + } + } + fact = alphaMin + (alphaMax-alphaMin)*rng.uniform(); + } + for (i=0; i<_eo1.size(); i++) { r1=_eo1[i]; @@ -165,6 +248,7 @@ template class eoSegmentCrossover: public eoQuadraticOp } protected: + eoRealVectorBounds & bounds; double alpha; double range; // == 1+2*alpha }; @@ -180,12 +264,35 @@ template class eoArithmeticCrossover: public eoQuadraticOp public: /** * (Default) Constructor. - * @param _alpha the amount of exploration OUTSIDE the parents + * The bounds are initialized with the global object that says: no bounds. + * + * @param _alphaMin the amount of exploration OUTSIDE the parents * as in BLX-alpha notation (Eshelman and Schaffer) * 0 == contractive application + * Must be positive */ eoArithmeticCrossover(const double& _alpha = 0.0): - alpha(_alpha), range(1+2*alpha) {} + bounds(eoDummyVectorNoBounds), alpha(_alpha), range(1+2*_alpha) + { + if (_alpha < 0) + throw runtime_error("BLX coefficient should be positive"); + } + + /** + * Constructor with bounds + * @param _bounds an eoRealVectorBounds that contains the bounds + * @param _alphaMin the amount of exploration OUTSIDE the parents + * as in BLX-alpha notation (Eshelman and Schaffer) + * 0 == contractive application + * Must be positive + */ + eoArithmeticCrossover(eoRealVectorBounds & _bounds, + const double& _alpha = 0.0): + bounds(_bounds), alpha(_alpha), range(1+2*_alpha) + { + if (_alpha < 0) + throw runtime_error("BLX coefficient should be positive"); + } /// The class name. string className() const { return "eoArithmeticCrossover"; } @@ -199,19 +306,50 @@ template class eoArithmeticCrossover: public eoQuadraticOp { unsigned i; double r1, r2, fact; - for (i=0; i<_eo1.size(); i++) - { - r1=_eo1[i]; - r2=_eo2[i]; - fact = rng.uniform(range); // in [0,range) - _eo1[i] = fact * r1 + (1-fact) * r2; - _eo2[i] = (1-fact) * r1 + fact * r2; - } + if (alpha == 0.0) // no check to perform + for (i=0; i<_eo1.size(); i++) + { + r1=_eo1[i]; + r2=_eo2[i]; + fact = -alpha + rng.uniform(range); // in [-alpha,1+alpha) + _eo1[i] = fact * r1 + (1-fact) * r2; + _eo2[i] = (1-fact) * r1 + fact * r2; + } + else // check the bounds + for (i=0; i<_eo1.size(); i++) + { + r1=_eo1[i]; + r2=_eo2[i]; + if (r1 != r2) { // otherwise you'll get NAN's + double rmin = min(r1, r2); + double rmax = max(r1, r2); + double length = rmax - rmin; + double alphaMin = -alpha; + double alphaMax = 1+alpha; + // first find the limits on the alpha's + if (bounds.isMinBounded(i)) + { + alphaMin = max(alphaMin, (bounds.minimum(i)-rmin)/length); + alphaMin = max(alphaMin, (rmax-bounds.maximum(i))/length); + } + if (bounds.isMaxBounded(i)) + { + alphaMax = min(alphaMax, (bounds.maximum(i)-rmin)/length); + alphaMax = min(alphaMax, (rmax-bounds.minimum(i))/length); + } + fact = alphaMin + rng.uniform(alphaMax-alphaMin); + _eo1[i] = fact * rmin + (1-fact) * rmax; + _eo2[i] = (1-fact) * rmin + fact * rmax; + } + } + _eo1.invalidate(); + _eo2.invalidate(); } protected: + eoRealVectorBounds & bounds; double alpha; - double range; // == 1+2*alpha + double range; // == 1+2*alphaMin };