src/algs/stogo/tools.cc

   1
   2 #include <float.h>
   3 #include <iostream>
   4
   5 #include "stogo_config.h"
   6 #include "tools.h"
   7
   8 Trial::Trial():xvals(0) {
   9   objval=DBL_MAX;
  10 }
  11
  12 Trial::Trial(int n):xvals(n) {
  13   objval=DBL_MAX;
  14 }
  15
  16 Trial::Trial(RCTrial tr):xvals(tr.xvals) {
  17   objval=tr.objval ;
  18 }
  19
  20 RCTrial Trial::operator=(RCTrial tr) {
  21   xvals=tr.xvals ;
  22   objval=tr.objval ;
  23   return *this ;
  24 }
  25
  26 ostream & operator << (ostream & os, const Trial & T) {
  27   os << T.xvals << "  " << "(" << T.objval << ")" << endl ;
  28   return os;
  29 }
  30
  31 /********************* Vanilla Box **********************/
  32 VBox::VBox():lb(0),ub(0) {} // Constructor
  33
  34 VBox::VBox(int n):lb(n),ub(n) {} // Constructor
  35
  36 VBox::VBox(RCVBox box):lb(box.lb), ub(box.ub) {} // Constructor
  37
  38 RCVBox VBox::operator=(RCVBox box) {
  39   // Copy: Box1=Box2
  40   lb=box.lb; ub=box.ub;
  41   return *this;
  42 }
  43
  44 int VBox::GetDim() {
  45   return lb.GetLength();
  46 }
  47
  48 double VBox::Width(int i) {
  49   // Returns the width of the i-th interval, i between [0,...,dim-1]
  50   return ub(i)-lb(i);
  51 }
  52
  53 void VBox::Midpoint(RCRVector x) {
  54   // Returns the midpoint of Box in x
  55   int n=GetDim();
  56   for (int i=0 ; i<n ; i++)
  57     x(i)=fabs(ub(i)-lb(i))/2 + lb(i);
  58 }
  59
  60 ostream & operator << (ostream & os, const VBox & B) {
  61   int n=(B.lb).GetLength() ;
  62   for (int i=0 ; i<n ;i++)
  63     os << '[' << B.lb(i) << "," << B.ub(i) << "]";
  64   return os;
  65 }
  66
  67 /************************ Trial Box ***********************/
  68 TBox::TBox() : VBox() {
  69   // Constructor
  70   minf=DBL_MAX;
  71 }
  72
  73 TBox::TBox(int n) : VBox(n) {
  74   // Constructor
  75   minf=DBL_MAX;
  76 }
  77
  78 TBox::TBox(RCTBox box) : VBox(box) {
  79   // Constructor + Copy
  80   minf=box.minf;
  81   TList=box.TList ;
  82 }
  83
  84 RCTBox TBox::operator=(RCTBox box) {
  85   // Copy
  86   // NB We should 'somehow' use VBox to copy lb and ub
  87   // Note that assignment operators are _not_ inherited
  88   lb=box.lb ; ub=box.ub ;
  89   minf=box.minf ;
  90   TList=box.TList ;
  91   return *this ;
  92 }
  93
  94 double TBox::GetMin() {
  95   return minf;
  96 }
  97
  98 bool TBox::EmptyBox() {
  99   // Returns TRUE if the list of Trials is empty
 100   return TList.empty() ;
 101 }
 102
 103 void TBox::AddTrial(RCTrial T) {
 104   // Add a Trial to the (back of) TList
 105   TList.push_back(T) ;
 106   if (T.objval<minf)
 107     minf=T.objval ;
 108 }
 109
 110 void TBox::RemoveTrial(Trial &T) {
 111   // Remove a trial from the (back of) box
 112   T=TList.back() ;
 113   TList.pop_back() ;
 114 }
 115
 116 void TBox::GetLastTrial(Trial &T) {
 117 // Return a trial from the (back of) box
 118   T=TList.back() ;
 119 }
 120
 121 list<Trial>::const_iterator TBox::FirstTrial() {
 122   return TList.begin();
 123 }
 124
 125 list<Trial>::const_iterator TBox::LastTrial() {
 126   return TList.end();
 127 }
 128
 129 void TBox::GetTrial(list<Trial>::const_iterator itr, Trial &T) {
 130   T.xvals=(*itr).xvals;
 131   T.objval=(*itr).objval;
 132 }
 133
 134 void TBox::ClearBox() {
 135   TList.erase(TList.begin(), TList.end());
 136   minf=DBL_MAX;
 137 }
 138
 139 bool TBox::CloseToMin(RVector &vec, double *objval, double eps_cl) {
 140   // Returns TRUE if 'vec' is close to some of the trials in the box,
 141   // in this case, 'vec' and 'objval' are overwritten by the Trial data
 142   // otherwise 'vec' and 'objval' are not affected.
 143   //
 144   // Here "close" means norm2(T - some trial in box)<=eps_cl
 145   //
 146   // NB It might be better to accept Trial as argument instead of vector
 147
 148   int n=GetDim();
 149   RVector x(n), y(n);
 150   list<Trial>::const_iterator itr;
 151   for ( itr = TList.begin(); itr != TList.end(); ++itr ) {
 152     y=vec ; // NB Should be possible to avoid this copying
 153     x=(*itr).xvals;
 154     axpy(-1,x,y);
 155     if (norm2(y)<=eps_cl) {
 156       vec=x;
 157       *objval=(*itr).objval;
 158       return TRUE;
 159     }
 160   }
 161   return FALSE;
 162 }
 163
 164 unsigned int TBox::NStationary() {
 165   // Returns the number of trials in a box
 166   return TList.size() ;
 167 }
 168
 169 void TBox::split(RTBox B1, RTBox B2) {
 170   list<Trial>::const_iterator itr;
 171   double w,m,tmp;
 172   double fm1=DBL_MAX, fm2=DBL_MAX;
 173   int i, k, ns;
 174   int n=GetDim();
 175
 176   B1.lb=lb; B1.ub=ub;
 177   B2.lb=lb; B2.ub=ub;
 178   w=LongestSide(&i);
 179   ns=TList.size();
 180   switch (ns) {
 181   case 0: case 1:
 182     // Bisection
 183     w=ub(i)-lb(i); // Length of interval
 184     m=lb(i)+w/2;   // Midpoint
 185     B1.ub(i)=m; B2.lb(i)=m;
 186     break;
 187   default:
 188     // More elaborate split when there are more than one trials in B
 189     // See Serguies and Kaj's tech. report, page 11
 190     // Compute the center point of all stationary points
 191     RVector center(n), x(n), dispers(n);
 192     center=0; dispers=0;
 193     for ( itr = TList.begin(); itr != TList.end(); ++itr )
 194       axpy(1.0, (*itr).xvals, center);
 195     scal((double)(1.0/ns),center);
 196
 197     // Compute the relative deviations
 198     for ( itr = TList.begin(); itr != TList.end(); ++itr ) {
 199       for (k = 0; k<n; k++) {
 200         x=(*itr).xvals;
 201         dispers(k)=dispers(k)+pow(center(k)-x(k),2.0);
 202       }
 203     }
 204     scal((double)(1.0/ns),dispers);
 205
 206     // i=arg max(disp)
 207     tmp=dispers(0);i=0;
 208     for (k=1; k<n; k++) {
 209       if (dispers(k)>tmp) {
 210         tmp=dispers(k);i=k;
 211       }
 212     }
 213     B1.ub(i)=center(i) ; B2.lb(i)=center(i);
 214     break;
 215   }
 216
 217   // Split the set of trials accordingly
 218   for ( itr = TList.begin(); itr != TList.end(); ++itr ) {
 219     if ( B1.InsideBox((*itr).xvals) ) {
 220       fm1=min(fm1,(*itr).objval);
 221       B1.AddTrial(*itr);
 222     }
 223     else {
 224       B2.AddTrial(*itr);
 225       fm2=min(fm2,(*itr).objval);
 226     }
 227   }
 228   // Set minf of B1 and B2
 229   B1.minf=fm1 ; B2.minf=fm2;
 230 }
 231
 232 void TBox::dispTrials() {
 233   // Display information about box
 234 #ifdef KCC
 235   copy(TList.begin(), TList.end(), ostream_iterator<Trial,char>(cout));
 236 #else
 237   copy(TList.begin(), TList.end(), ostream_iterator<Trial>(cout));
 238 #endif
 239 }
 240
 241 ostream & operator << (ostream & os, const TBox & B) {
 242   int n=(B.lb).GetLength() ;
 243   for (int i=0 ; i<n ;i++)
 244     os << '[' << B.lb(i) << "," << B.ub(i) << "]";
 245   os << "   minf= " << B.minf << endl;
 246   return os;
 247 }
 248
 249 bool TBox::InsideBox(RCRVector x) {
 250   // Returns TRUE if the point X lies inside BOX, FALSE otherwise
 251   int n=GetDim();
 252   for (int i=0 ; i<n ; i++)
 253     if (x(i)<lb(i) || x(i)>ub(i)) return FALSE;
 254   return TRUE;
 255 }
 256
 257 int TBox::OutsideBox(RCRVector x, RCTBox domain) {
 258   // The function returns
 259   //    0 if X is inside both 'box' and 'domain'
 260   //    1 if X is inside 'domain' but outside 'box'
 261   //    2 if X is outside both 'domain' and 'box
 262
 263   int n=GetDim();
 264   int ins_box=1, ins_dom=1, outs=999;
 265   for (int i=0 ; i<n ; i++) {
 266     if (x(i)<lb(i) || x(i)>ub(i))
 267       ins_box=0;
 268     if (x(i)<domain.lb(i) || x(i)>domain.ub(i)) {
 269       ins_dom=0; break;
 270     }
 271   }
 272   if (ins_box==1 && ins_dom==1)
 273     outs=0;
 274   if (ins_box==0 && ins_dom==1)
 275     outs=1;
 276   if (ins_box==0 && ins_dom==0)
 277     outs=2;
 278   if (outs==999) {
 279     // Something has gone wrong!
 280     cout << "Error in OutsideBox, exiting\n";
 281     exit(1);
 282   }
 283   return outs;
 284 }
 285
 286 double TBox::ShortestSide(int *idx) {
 287   // Returns the shortest side of the box and it's index
 288   int n=GetDim(),j=0;
 289   double tmp=ub(0)-lb(0);
 290   for (int i=1 ; i<n ;i++)
 291     if ( (ub(i)-lb(i))<tmp ) {
 292       tmp=ub(i)-lb(i);
 293       j=i;
 294     }
 295   *idx=j;
 296   return tmp;
 297 }
 298
 299 double TBox::LongestSide(int *idx) {
 300   // Returns the longest side of the box and it's index
 301   int n=GetDim(),j=0;
 302   double tmp=ub(0)-lb(0);
 303   for (int i=1 ; i<n ;i++)
 304     if ( (ub(i)-lb(i))>tmp ) {
 305       tmp=ub(i)-lb(i);
 306       j=i;
 307     }
 308   *idx=j;
 309   return tmp;
 310 }
 311
 312 double TBox::ClosestSide(RCRVector x) {
 313   // Returns the shortest distance from point X to the box B.
 314
 315   //   Warning: The output of this functon is nonsense if the
 316   //   point X lies outside B. Should we try to detect this case?
 317
 318   double dist, tmp ;
 319   int n=GetDim();
 320   dist=DBL_MAX;
 321   for (int i=0 ; i<n ; i++) {
 322     tmp = min( x(i)-lb(i), ub(i)-x(i) );
 323     dist = min(dist, tmp);
 324   }
 325   return dist;
 326 }
 327
 328 double TBox::FarthestSide(RCRVector x) {
 329   // Returns the longest distance from point X to the box B.
 330   //   Same comment apply here as in ClosestSide(X)
 331
 332   double dist, tmp;
 333   int n=GetDim();
 334   dist=DBL_MIN;
 335   for (int i=0 ; i<n ; i++) {
 336     tmp = max( x(i)-lb(i), ub(i)-x(i) );
 337     dist = max(dist, tmp);
 338   }
 339   return dist;
 340 }
 341
 342 bool TBox::Intersection(RCRVector x, RCRVector h, RCRVector z) {
 343 // Intersection of a line with a box
 344 // The line is described by a vector 'h' and a point 'x' and
 345 // the point of intersection is returned in 'z'
 346 // Note that h+x must lie outside of box
 347 //
 348 // Known problem:
 349 //   Due to round of errors the algorithm can fail to find an intersection
 350 //   The caller is notified and should act accordingly
 351 //
 352 //  The routine returns FALSE if no intersection was found, TRUE otherwise
 353
 354   int n=GetDim();
 355   RVector tmpV(n);
 356   bool done;
 357   int i, j, k, isect;
 358   double alpha, gamma;
 359
 360   i=0; done=FALSE;
 361   while (i<n && done==FALSE) {
 362     if (h(i)==0) {
 363       z(i)=x(i);
 364       break;
 365     }
 366     for (k=1; k<=2; k++) {
 367       if (k==1)
 368         alpha=lb(i);
 369       else
 370         alpha=ub(i);
 371       gamma=(alpha-x(i))/h(i);
 372       z(i)=alpha;
 373       isect=1;
 374       for (j=0; j<n; j++) {
 375         if (j != i) {
 376           z(j)=x(j)+gamma*h(j);
 377           if (z(j)<lb(j) || z(j)>ub(j)) {
 378             isect=0;
 379             break;
 380           }
 381         }
 382       }
 383       copy(z,tmpV); axpy(-1.0,x,tmpV);  // tmpV=z-x
 384       if (isect==1 && dot(tmpV,h)>0) {
 385         done=TRUE; break;
 386       }
 387     }
 388     i++;
 389   }
 390   return done;
 391 }
 392
 393 double TBox::LowerBound(double maxgrad) {
 394 // Lower bound estimation
 395   double lb=minf ;
 396   double f1,f2,est ;
 397   list<Trial>::const_iterator itr1,itr2 ;
 398
 399   int n=GetDim();
 400   RVector x1(n), x2(n) ;
 401 #ifndef LB2
 402   for ( itr1 = TList.begin(); itr1 != TList.end(); ++itr1 ) {
 403     itr2=itr1 ;
 404     while (++itr2 != TList.end()) {
 405       x1=(*itr1).xvals ; f1=(*itr1).objval ;
 406       x2=(*itr2).xvals ; f2=(*itr2).objval ;
 407       axpy(-1.0,x2,x1) ; // x1=x1-x2
 408       est=0.5*(f1+f2-maxgrad*norm2(x1)) ;
 409       lb=min(lb,est) ;
 410       // cout << "est=" << est << " x1=" << x1 << " x2=" << x2 << endl ;
 411     }
 412   }
 413 #endif
 414 #ifdef LB2
 415   for ( itr1 = TList.begin(); itr1 != TList.end(); ++itr1 ) {
 416     // Use also max distance to border
 417     x1=(*itr1).xvals; f1=(*itr1).objval;
 418     est=f1-FarthestSide(x1)*maxgrad;
 419     lb=min(lb,est);
 420   }
 421 #endif
 422   return lb ;
 423 }