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