stogo/global.cc

   1 /*
   2    Multi Dimensional Global Search.
   3
   4    Author: Steinn Gudmundsson
   5    Email: steinng@hotmail.com
   6
   7    This program is supplied without any warranty whatsoever.
   8
   9    NB The RNGs seed should be initialized using some timer
  10 */
  11
  12 #include <iostream>
  13 #include <time.h>
  14
  15 #include <iterator>
  16 #include <algorithm>
  17 #include <stack>
  18
  19 #include "stogo_config.h"
  20 #include "global.h"
  21 #include "local.h"
  22
  23 // Timer stuff
  24 time_t   StartTime;
  25
  26 double MacEpsilon ;
  27 int FC=0, GC=0 ;
  28
  29 int stogo_verbose = 0; /* set to nonzero for verbose output */
  30
  31 Global::Global(RTBox D, Pobj o, Pgrad g, GlobalParams P): Domain(D) {
  32
  33   dim=Domain.GetDim();
  34   Objective=o;
  35   Gradient=g;
  36
  37   // Initialize parameters
  38   maxtime=P.maxtime;
  39   maxeval = P.maxeval;
  40   numeval = 0;
  41   eps_cl=P.eps_cl; mu=P.mu; rshift=P.rshift;
  42   det_pnts=P.det_pnts; rnd_pnts=P.rnd_pnts;
  43   fbound=DBL_MAX;
  44 }
  45
  46 #if 0 // not necessary; default copy is sufficient
  47 Global& Global::operator=(const Global &G) {
  48   // Copy the problem info and parameter settings
  49   Domain=G.Domain; Objective=G.Objective;  Gradient=G.Gradient;
  50   maxtime=G.maxtime;
  51   maxeval = G.maxeval;
  52   numeval = G.numeval;
  53   eps_cl=G.eps_cl; mu=G.mu; rshift=G.rshift;
  54   det_pnts=G.det_pnts; rnd_pnts=G.rnd_pnts;
  55   return *this;
  56 }
  57 #endif
  58
  59 void Global::FillRegular(RTBox SampleBox, RTBox box) {
  60   // Generation of regular sampling points
  61   double w;
  62   int i, flag, dir;
  63   Trial tmpTrial(dim);
  64   RVector m(dim), x(dim);
  65
  66   if (det_pnts>0) {
  67     // Add midpoint
  68     box.Midpoint(m) ;
  69     tmpTrial.xvals=m ; tmpTrial.objval=DBL_MAX ;
  70     SampleBox.AddTrial(tmpTrial) ;
  71     // Add the rest
  72     i=1 ; flag=1 ; dir=0 ;
  73     x=m ; w=box.Width(dir) ;
  74     while (i<det_pnts) {
  75       x(dir)=m(dir)+flag*rshift*w ;
  76       tmpTrial.xvals=x ;
  77       SampleBox.AddTrial(tmpTrial) ;
  78       flag=-flag;
  79       if (flag==1 && dir<dim) {
  80         x(dir)=m(dir) ;
  81         dir++ ;
  82         w=box.Width(dir) ;
  83       }
  84       i++ ;
  85     }
  86   }
  87 }
  88
  89 void Global::FillRandom(RTBox SampleBox, RTBox box) {
  90   // Generation of stochastic sampling points
  91   Trial tmpTrial(dim);
  92
  93   tmpTrial.objval=DBL_MAX;
  94   for (int i=1 ; i<=rnd_pnts ; i++) {
  95     for (int dir=0 ; dir<dim ; dir++)
  96       tmpTrial.xvals(dir) =
  97         box.lb(dir)+(box.ub(dir)-box.lb(dir))*(double(rand())/RAND_MAX) ;
  98     SampleBox.AddTrial(tmpTrial) ;
  99   }
 100 }
 101
 102 double Global::NewtonTest(RTBox box, int axis, RCRVector x_av, int *noutside) {
 103   // Perform the Newton test
 104
 105   int info,nout=0;
 106   Trial tmpTrial(dim);
 107   TBox SampleBox(dim) ;
 108   double maxgrad=0 ;
 109
 110   // Create sampling points
 111   FillRegular(SampleBox, box);
 112   FillRandom(SampleBox, box);
 113
 114   // Perform the actual sampling
 115   while ( !SampleBox.EmptyBox() ) {
 116     SampleBox.RemoveTrial(tmpTrial) ;
 117     info = local(tmpTrial, box, Domain, eps_cl, &maxgrad, *this,
 118                  axis, x_av) ;
 119     // What should we do when info=LS_Unstable?
 120     if (info == LS_Out)
 121       nout++;
 122     if (info == LS_New ) {
 123       box.AddTrial(tmpTrial) ;
 124
 125       if (tmpTrial.objval<=fbound+mu && tmpTrial.objval<=box.fmin+mu) {
 126         if (stogo_verbose) {
 127           cout << "Found a candidate, x=" << tmpTrial.xvals;
 128           cout << " F=" <<tmpTrial.objval << " FC=" << FC << endl;
 129         }
 130         SolSet.push_back(tmpTrial);
 131       }
 132 #ifdef GS_DEBUG
 133       cout << "Found a stationary point, X= " << tmpTrial.xvals;
 134       cout <<" objval=" << tmpTrial.objval << endl;
 135 #endif
 136     }
 137
 138     if (!InTime())
 139       break;
 140   }
 141   *noutside=nout;
 142   return maxgrad;
 143 }
 144
 145 void Global::ReduceOrSubdivide(RTBox box, int axis, RCRVector x_av) {
 146   TBox B1(dim), B2(dim);
 147   Trial tmpTrial(dim);
 148   double maxgrad;
 149   int ns,nout;
 150
 151   // Monotonicity test has not been implemented yet
 152   maxgrad=NewtonTest(box, axis, x_av, &nout);
 153   ns=box.NStationary() ;
 154   if (ns==0) {
 155     // All iterates outside
 156     // NB result=Intersection(B,boundary(Domain))
 157     Garbage.push(box) ;
 158   }
 159   else
 160     if (ns==1 && nout==0) {
 161       // All iterates converge to same point
 162       Garbage.push(box) ;
 163     }
 164     else
 165       if ( (ns>1) && (box.LowerBound(maxgrad)>fbound) ) {
 166         // Several stationary points found and lower bound > fbound
 167         Garbage.push(box) ;
 168       }
 169       else {
 170         // Subdivision
 171         B1.ClearBox() ; B2.ClearBox() ;
 172         box.split(B1,B2) ;
 173         CandSet.push(B1) ; CandSet.push(B2) ;
 174       }
 175
 176   // Update fbound
 177   if (box.fmin < fbound) {
 178     fbound=box.fmin ;
 179 #ifdef GS_DEBUG
 180     cout <<"*** Improving fbound, fbound=" << fbound << endl;
 181 #endif
 182   }
 183 }
 184
 185 void Global::Search(int axis, RCRVector x_av){
 186   Trial tmpTrial(dim) ;
 187   TBox box(dim), B1(dim), B2(dim);
 188   RVector m(dim), x(dim);
 189   int inner_iter, outer_iter;
 190
 191   MacEpsilon=eps(); // Get machine precision
 192   if (det_pnts>2*dim+1) {
 193     det_pnts=2*dim+1;
 194     if (stogo_verbose)
 195       cout << "Warning: Reducing det_pnts to " << det_pnts << endl;
 196   }
 197
 198   // Initialize timer
 199   time(&StartTime);
 200
 201   // Clear priority_queues
 202   while (!Garbage.empty())
 203     Garbage.pop();
 204   while (!CandSet.empty())
 205     CandSet.pop();
 206
 207   box=Domain;
 208   CandSet.push(box);
 209   int done=0 ; outer_iter=0 ;
 210
 211   while (!done) {
 212     outer_iter++ ;
 213
 214     // Inner loop
 215     inner_iter=0 ;
 216     while (!CandSet.empty()) {
 217       inner_iter++ ;
 218       // Get best box from Candidate set
 219       box=CandSet.top() ; CandSet.pop() ;
 220
 221 #ifdef GS_DEBUG
 222       cout << "Iteration..." << inner_iter << " #CS=" << CandSet.size()+1 ;
 223       cout << " Processing " << box.NStationary() << " trials in the box " <<box;
 224 #endif
 225       ReduceOrSubdivide(box, axis, x_av);
 226
 227       if (!InTime()) {
 228         done=TRUE;
 229         if (stogo_verbose)
 230           cout << "The program has run out of time or function evaluations\n";
 231         break;
 232       }
 233
 234     } // inner while-loop
 235     if (stogo_verbose)
 236       cout << endl << "*** Inner loop completed ***" << endl ;
 237
 238     // Reduce SolSet if necessary
 239     SolSet.erase(remove_if(SolSet.begin(), SolSet.end(),
 240                            TrialGT(fbound+mu)),SolSet.end());
 241     if (InTime()) {
 242       if (stogo_verbose) {
 243         cout << "Current set of minimizers (" << SolSet.size() << ")" << endl ;
 244         DispMinimizers() ;
 245       }
 246
 247       while (!Garbage.empty()) {
 248         box=Garbage.top() ;
 249         Garbage.pop() ;
 250         // Split box
 251         B1.ClearBox() ; B2.ClearBox() ;
 252         box.split(B1,B2) ;
 253         // Add boxes to Candidate set
 254         CandSet.push(B1) ; CandSet.push(B2) ;
 255       }
 256     }
 257   } // Outer while-loop
 258
 259   if (stogo_verbose) {
 260     cout << "Number of outer iterations : " << outer_iter << endl;
 261     cout << "Number of unexplored boxes : " << CandSet.size() << endl;
 262     cout << "Number of boxes in garbage : " << Garbage.size() << endl;
 263     cout << "Number of elements in SolSet : " << SolSet.size() << endl;
 264     cout << "Number of function evaluations : " << FC << endl;
 265     cout << "Number of gradient evaluations : " << GC << endl;
 266   }
 267
 268   if (axis != -1) {
 269     // Return minimizer when doing the AV method
 270     tmpTrial=SolSet.back();
 271     x_av(axis)=tmpTrial.xvals(0);
 272   }
 273 }
 274
 275 /************* Various utility functions ****************/
 276 double Global::GetTime()
 277 {
 278  time_t ctime; time(&ctime);
 279  return difftime(ctime,StartTime);
 280 }
 281
 282 bool Global::InTime()
 283 {
 284  return (maxtime <= 0.0 || GetTime()<maxtime) && (!maxeval || numeval<maxeval);
 285 }
 286
 287 double Global::GetMinValue() {
 288   return fbound;
 289 }
 290
 291 void Global::SetMinValue(double new_fb) {
 292   fbound=new_fb;
 293 }
 294
 295 void Global::SetDomain(RTBox box) {
 296   Domain=box;
 297 }
 298
 299 void Global::GetDomain(RTBox box) {
 300   box=Domain;
 301 }
 302
 303 void Global::DispMinimizers() {
 304   copy(SolSet.begin(), SolSet.end(), ostream_iterator<Trial>(cout));
 305 }
 306
 307 double Global::OneMinimizer(RCRVector x) {
 308   if (NoMinimizers()) return 0.0;
 309   for (int i=0;i<x.GetLength();i++) x(i) = SolSet.front().xvals(i);
 310   return SolSet.front().objval;
 311 }
 312
 313 bool Global::NoMinimizers() {
 314   return SolSet.empty();
 315 }
 316
 317 void Global::ClearSolSet() {
 318   SolSet.erase(SolSet.begin(), SolSet.end()) ;
 319 }
 320
 321 void Global::AddPoint(RCRVector x, double f) {
 322   Trial T(dim);
 323   T.xvals=x; T.objval=f;
 324   Domain.AddTrial(T);
 325   SolSet.push_back(T);
 326 }