2 Multi Dimensional Global Search.
4 Author: Steinn Gudmundsson
5 Email: steinng@hotmail.com
7 This program is supplied without any warranty whatsoever.
9 NB The RNGs seed should be initialized using some timer
18 #include "stogo_config.h"
21 #include "nlopt-util.h"
29 int stogo_verbose = 0; /* set to nonzero for verbose output */
31 Global::Global(RTBox D, Pobj o, Pgrad g, GlobalParams P): Domain(D) {
37 // Initialize parameters
45 eps_cl=P.eps_cl; mu=P.mu; rshift=P.rshift;
46 det_pnts=P.det_pnts; rnd_pnts=P.rnd_pnts;
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;
61 eps_cl=G.eps_cl; mu=G.mu; rshift=G.rshift;
62 det_pnts=G.det_pnts; rnd_pnts=G.rnd_pnts;
67 void Global::FillRegular(RTBox SampleBox, RTBox box) {
68 // Generation of regular sampling points
72 RVector m(dim), x(dim);
76 tmpTrial.objval=DBL_MAX ;
78 i=1 ; flag=1 ; dir=0 ;
82 x(dir)=m(dir)+flag*rshift*w ;
84 SampleBox.AddTrial(tmpTrial) ;
86 if (flag==1 && dir<dim) {
94 SampleBox.AddTrial(tmpTrial) ;
98 void Global::FillRandom(RTBox SampleBox, RTBox box) {
99 // Generation of stochastic sampling points
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) ;
110 double Global::NewtonTest(RTBox box, int axis, RCRVector x_av, int *noutside) {
111 // Perform the Newton test
115 TBox SampleBox(dim) ;
118 // Create sampling points
119 FillRandom(SampleBox, box);
120 FillRegular(SampleBox, box);
122 // Perform the actual sampling
123 while ( !SampleBox.EmptyBox() ) {
124 SampleBox.RemoveTrial(tmpTrial) ;
125 info = local(tmpTrial, box, Domain, eps_cl, &maxgrad, *this,
131 // What should we do when info=LS_Unstable?
134 else if (info == LS_New ) {
135 box.AddTrial(tmpTrial) ;
137 if (tmpTrial.objval<=fbound+mu && tmpTrial.objval<=box.minf+mu) {
139 cout << "Found a candidate, x=" << tmpTrial.xvals;
140 cout << " F=" <<tmpTrial.objval << " FC=" << FC << endl;
142 SolSet.push_back(tmpTrial);
144 if (tmpTrial.objval < stop->minf_max)
149 cout << "Found a stationary point, X= " << tmpTrial.xvals;
150 cout <<" objval=" << tmpTrial.objval << endl;
154 if (!InTime() || info == LS_MaxEvalTime)
161 void Global::ReduceOrSubdivide(RTBox box, int axis, RCRVector x_av) {
162 TBox B1(dim), B2(dim);
167 // Monotonicity test has not been implemented yet
168 maxgrad=NewtonTest(box, axis, x_av, &nout);
169 ns=box.NStationary() ;
171 // All iterates outside
172 // NB result=Intersection(B,boundary(Domain))
176 if (ns==1 && nout==0) {
177 // All iterates converge to same point
181 if ( (ns>1) && (box.LowerBound(maxgrad)>fbound) ) {
182 // Several stationary points found and lower bound > fbound
187 B1.ClearBox() ; B2.ClearBox() ;
189 CandSet.push(B1) ; CandSet.push(B2) ;
193 if (box.minf < fbound) {
196 cout <<"*** Improving fbound, fbound=" << fbound << endl;
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;
207 MacEpsilon=eps(); // Get machine precision
208 if (det_pnts>2*dim+1) {
211 cout << "Warning: Reducing det_pnts to " << det_pnts << endl;
215 StartTime = nlopt_seconds();
217 // Clear priority_queues
218 while (!Garbage.empty())
220 while (!CandSet.empty())
225 int done=0 ; outer_iter=0 ;
232 while (!CandSet.empty()) {
234 // Get best box from Candidate set
235 box=CandSet.top() ; CandSet.pop() ;
238 cout << "Iteration..." << inner_iter << " #CS=" << CandSet.size()+1 ;
239 cout << " Processing " << box.NStationary() << " trials in the box " <<box;
241 ReduceOrSubdivide(box, axis, x_av);
244 if (!NoMinimizers() && OneMinimizer(x) < stop->minf_max) {
252 cout << "The program has run out of time or function evaluations\n";
256 } // inner while-loop
258 cout << endl << "*** Inner loop completed ***" << endl ;
260 // Reduce SolSet if necessary
261 SolSet.erase(remove_if(SolSet.begin(), SolSet.end(),
262 TrialGT(fbound+mu)),SolSet.end());
265 cout << "Current set of minimizers (" << SolSet.size() << ")" << endl ;
269 while (!Garbage.empty()) {
273 B1.ClearBox() ; B2.ClearBox() ;
275 // Add boxes to Candidate set
276 CandSet.push(B1) ; CandSet.push(B2) ;
279 } // Outer while-loop
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;
291 // Return minimizer when doing the AV method
292 tmpTrial=SolSet.back();
293 x_av(axis)=tmpTrial.xvals(0);
297 /************* Various utility functions ****************/
298 double Global::GetTime()
300 return nlopt_seconds() - StartTime;
303 bool Global::InTime()
306 return !nlopt_stop_evalstime(stop);
308 return (maxtime <= 0.0 || GetTime()<maxtime) && (!maxeval || numeval<maxeval);
312 double Global::GetMinValue() {
316 void Global::SetMinValue(double new_fb) {
320 void Global::SetDomain(RTBox box) {
324 void Global::GetDomain(RTBox box) {
328 void Global::DispMinimizers() {
329 copy(SolSet.begin(), SolSet.end(), ostream_iterator<Trial>(cout));
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;
338 bool Global::NoMinimizers() {
339 return SolSet.empty();
342 void Global::ClearSolSet() {
343 SolSet.erase(SolSet.begin(), SolSet.end()) ;
346 void Global::AddPoint(RCRVector x, double f) {
348 T.xvals=x; T.objval=f;