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