Added the bounds to all eoRealOp operators (and to the eoNormalMutation).

The eoESObjectiveBounds should soon dissapear!!

eo/src/es/eoNormalMutation.h

@@ -33,12 +33,16 @@
 #include <utils/eoUpdatable.h>
 #include <eoEvalFunc.h>
 #include <es/eoReal.h>
-
+#include <es/eoRealBounds.h>
 //-----------------------------------------------------------------------------
 
 /** 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 EOT> class eoNormalMutation: public eoMonOp<EOT>
@@ -46,11 +50,23 @@ template<class EOT> class eoNormalMutation: public eoMonOp<EOT>
  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 EOT> class eoNormalMutation: public eoMonOp<EOT>
 	  if (rng.flip(p_change)) 
 	    {
 	      _eo[lieu] += sigma*rng.normal();
+	      bounds.foldsInBounds(lieu, _eo[lieu]);
 	      hasChanged = true;
 	    }
 	}
@@ -77,6 +94,7 @@ template<class EOT> class eoNormalMutation: public eoMonOp<EOT>
 protected:
   double & sigma;
 private:
+  eoRealVectorBounds & bounds;
   double p_change;
 };
 

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<double, bool>
+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; }
+  // accessors  
-  double Maximum() { return repMaximum; }
+  virtual double minimum() { return repMinimum; }
-  double Range()   { return repRange; }
+  virtual double maximum() { return repMaximum; }
-  // for backward compatibility
+  virtual double range()   { return repRange; }
-  double minimum() { return repMinimum; }
+
-  double maximum() { return repMaximum; }
+  // description
-  double range()   { return repRange; }
+  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 
+class eoRealVectorBounds : public vector<eoRealBounds *>
 { 
 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<eoRealBounds *>(_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()), 
+    vector<eoRealBounds *>(_dim, &_bounds)
-    vecMaximum(_dim, _bounds.Maximum()), 
-    vecRange(_dim, _bounds.Range()) 
   {}
   
   // Ctor: different bonds for different variables, vectors of double
-  eoRealVectorBounds(vector<double> _min, vector<double> _max) : 
+  eoRealVectorBounds(vector<double> _min, vector<double> _max) 
-    vecMinimum(_min), vecMaximum(_max), vecRange(_min.size()) 
   {
     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];
+	push_back( new eoRealInterval(_min[i], _max[i]));
-	if (vecRange[i]<=0)
-	  throw std::logic_error("Void range in eoRealVectorBounds");
       }
   }
 
   // Ctor, particular case of dim-2
   eoRealVectorBounds(eoRealBounds & _xbounds, eoRealBounds & _ybounds) : 
-    vecMinimum(2), vecMaximum(2), vecRange(2) 
+    vector<eoRealBounds *>(0)
   {
-    vecMinimum[0] = _xbounds.Minimum();
+	push_back( &_xbounds);
-    vecMaximum[0] = _xbounds.Maximum();
+	push_back( &_ybounds);
-    vecRange[0] = _xbounds.Range();
-    vecMinimum[1] = _ybounds.Minimum();
-    vecMaximum[1] = _ybounds.Maximum();
-    vecRange[1] = _ybounds.Range();    
   }
   
-  // not a ctor, but usefull to initialize, too
+  virtual bool isBounded(unsigned _i) 
-  // is it safe to call it push_back? Maybe not, but it's meaningful!
-  void push_back(double _min=0, double _max=1)
   { 
-    vecMinimum.push_back(_min);
+    return (*this)[_i]->isBounded();
-    vecMaximum.push_back(_max);
-    if (_max-_min <= 0)
-      throw std::logic_error("Void range in eoRealVectorBounds::add");
-    vecRange.push_back(_max-_min);
   }
  
-  void push_back(eoRealBounds & _bounds)
+  // bounded iff all are bounded
+  virtual bool isBounded(void) 
   {
-    vecMinimum.push_back(_bounds.Minimum());
+    for (unsigned i=0; i<size(); i++)
-    vecMaximum.push_back(_bounds.Maximum());
+      if (! (*this)[i]->isBounded())
-    vecRange.push_back(_bounds.Range());
+	return false;
+    return true;
   }
 
-  // accessors - following rule that says that method start are capitalized
+  // these do not make any sense as vectors!
-  double Minimum(unsigned _i) { return vecMinimum[_i]; }
+  virtual bool isMinBounded(unsigned _i) 
-  double Maximum(unsigned _i) { return vecMaximum[_i]; }
+  { return (*this)[_i]->isMinBounded();} ;
-  double Range(unsigned _i)   { return vecRange[_i]; }
 
-  // accessors - for backward compatibility
+  virtual bool isMaxBounded(unsigned _i) 
-  double minimum(unsigned _i) { return vecMinimum[_i]; }
+  { return (*this)[_i]->isMaxBounded();} ;
-  double maximum(unsigned _i) { return vecMaximum[_i]; }
-  double range(unsigned _i)   { return vecRange[_i]; }
 
-  // handy: get the size
+  virtual void foldsInBounds(unsigned _i, double & _r)
-  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)
   {
-    return vecMinimum[_i] + _rng.uniform(vecRange[_i]);
+    (*this)[_i]->foldsInBounds(_r);
   }
 
-  // returns a vector of uniformly chosen variables in bounds
+  virtual void foldsInBounds(vector<double> & _v)
-  vector<double> uniform(eoRng & _rng = eo::rng)
   {
-    vector<double> v(vecMinimum.size());
+   for (unsigned i=0; i<size(); i++)
-    for (unsigned i=0; i<vecMinimum.size(); i++)
+     {
-      v[i] = vecMinimum[i] + _rng.uniform(vecRange[i]);
+       foldsInBounds(i, _v[i]);
+     }    
+  }
 
-    return v;
+  virtual bool isInBounds(unsigned _i, double _r)
+  { return (*this)[_i]->isInBounds(_r); }
+
+  // isInBounds iff all are in bouds
+  virtual bool isInBounds(vector<double> _v)
+  {
+    for (unsigned i=0; i<size(); i++)
+      if (! isInBounds(i, _v[i]))
+	return false;
+    return true;
+  }
+
+  // accessors  
+  virtual double minimum(unsigned _i) {return (*this)[_i]->minimum();}
+  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; i<size(); i++)
+      r += range(i);
+    return r/size();
+  }
+
+  // random generators
+  virtual double uniform(unsigned _i, eoRng & _rng = eo::rng)
+  { 
+    double r= (*this)[_i]->uniform();
+    return r;
   }
 
   // fills a vector with uniformly chosen variables in bounds
   void uniform(vector<double> & _v, eoRng & _rng = eo::rng)
   {
-    _v.resize(vecMinimum.size());
+    _v.resize(size());
-    for (unsigned i=0; i<vecMinimum.size(); i++)
+    for (unsigned i=0; i<size(); i++)
-      _v[i] = vecMinimum[i] + _rng.uniform(vecRange[i]);
-  }  
-
-  // says if a given double is within the bounds
-  bool operator()(unsigned _i, double _r)
       {
-    if (_r < vecMinimum[_i])
+      _v[i] = uniform(i, _rng);
-      return false;
-    if (_r > vecMaximum[_i])
-      return false;
-    return true;
       }
-
-  // check the bounds for a vector: true only if ALL ar ein bounds
-  bool operator()(vector<double> & _v)
-  {
-    for (unsigned i=0; i<_v.size(); i++)
-      if (! operator()(i, _v[i]) ) // out of bound
-	return false;
-    return true;
   }  
-private :
-  vector<double> vecMinimum;
-  vector<double> vecMaximum;
-  vector<double> 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<double> &) {return;}
+
+  virtual bool isInBounds(unsigned, double) {return true;}
+  virtual bool isInBounds(vector<double>) {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<double> &, eoRng & _rng = eo::rng)
+  {
+    throw logic_error("No uniform distribution on unbounded eoRealBounds");
+  }
+
+};
+
+// one object for all
+eoRealNoBounds eoDummyRealNoBounds;
+eoRealVectorNoBounds eoDummyVectorNoBounds;
 #endif

eo/src/es/eoRealOp.h

@@ -31,6 +31,7 @@
 #include <algorithm>    // swap_ranges
 #include <utils/eoRNG.h>
 #include <es/eoReal.h>
+#include <es/eoRealBounds.h>
 
 //-----------------------------------------------------------------------------
 
@@ -46,11 +47,23 @@ template<class EOT> class eoUniformMutation: public eoMonOp<EOT>
  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 EOT> class eoUniformMutation: public eoMonOp<EOT>
 	{
 	  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 EOT> class eoUniformMutation: public eoMonOp<EOT>
     }
   
 private:
+  eoRealVectorBounds & bounds;
   double epsilon;
   double p_change;
 };
@@ -94,7 +115,17 @@ template<class EOT> class eoDetUniformMutation: public eoMonOp<EOT>
    * @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 EOT> class eoDetUniformMutation: public eoMonOp<EOT>
 	{
 	  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 EOT> class eoSegmentCrossover: public eoQuadraticOp<EOT>
  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 EOT> class eoSegmentCrossover: public eoQuadraticOp<EOT>
     {
       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 EOT> class eoSegmentCrossover: public eoQuadraticOp<EOT>
     }
 
 protected:
+  eoRealVectorBounds & bounds;
   double alpha;
   double range;			   // == 1+2*alpha
 };
@@ -180,12 +264,35 @@ template<class EOT> class eoArithmeticCrossover: public eoQuadraticOp<EOT>
  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 EOT> class eoArithmeticCrossover: public eoQuadraticOp<EOT>
     {
       unsigned i;
       double r1, r2, fact;
+      if (alpha == 0.0)		   // no check to perform
 	  for (i=0; i<_eo1.size(); i++)
 	      {
 		r1=_eo1[i];
 		r2=_eo2[i];
-	  fact = rng.uniform(range);	   // in [0,range)
+	    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
 };