5 #include "stogo_config.h"
8 Trial::Trial():xvals(0) {
12 Trial::Trial(int n):xvals(n) {
16 Trial::Trial(RCTrial tr):xvals(tr.xvals) {
20 RCTrial Trial::operator=(RCTrial tr) {
26 ostream & operator << (ostream & os, const Trial & T) {
27 os << T.xvals << " " << "(" << T.objval << ")" << endl ;
31 /********************* Vanilla Box **********************/
32 VBox::VBox():lb(0),ub(0) {} // Constructor
34 VBox::VBox(int n):lb(n),ub(n) {} // Constructor
36 VBox::VBox(RCVBox box):lb(box.lb), ub(box.ub) {} // Constructor
38 RCVBox VBox::operator=(RCVBox box) {
45 return lb.GetLength();
48 double VBox::Width(int i) {
49 // Returns the width of the i-th interval, i between [0,...,dim-1]
53 void VBox::Midpoint(RCRVector x) {
54 // Returns the midpoint of Box in x
56 for (int i=0 ; i<n ; i++)
57 x(i)=fabs(ub(i)-lb(i))/2 + lb(i);
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) << "]";
67 /************************ Trial Box ***********************/
68 TBox::TBox() : VBox() {
73 TBox::TBox(int n) : VBox(n) {
78 TBox::TBox(RCTBox box) : VBox(box) {
84 RCTBox TBox::operator=(RCTBox box) {
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 ;
94 double TBox::GetMin() {
98 bool TBox::EmptyBox() {
99 // Returns TRUE if the list of Trials is empty
100 return TList.empty() ;
103 void TBox::AddTrial(RCTrial T) {
104 // Add a Trial to the (back of) TList
110 void TBox::RemoveTrial(Trial &T) {
111 // Remove a trial from the (back of) box
116 void TBox::GetLastTrial(Trial &T) {
117 // Return a trial from the (back of) box
121 list<Trial>::const_iterator TBox::FirstTrial() {
122 return TList.begin();
125 list<Trial>::const_iterator TBox::LastTrial() {
129 void TBox::GetTrial(list<Trial>::const_iterator itr, Trial &T) {
130 T.xvals=(*itr).xvals;
131 T.objval=(*itr).objval;
134 void TBox::ClearBox() {
135 TList.erase(TList.begin(), TList.end());
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.
144 // Here "close" means norm2(T - some trial in box)<=eps_cl
146 // NB It might be better to accept Trial as argument instead of vector
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
155 if (norm2(y)<=eps_cl) {
157 *objval=(*itr).objval;
164 unsigned int TBox::NStationary() {
165 // Returns the number of trials in a box
166 return TList.size() ;
169 void TBox::split(RTBox B1, RTBox B2) {
170 list<Trial>::const_iterator itr;
172 double fm1=DBL_MAX, fm2=DBL_MAX;
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;
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);
193 for ( itr = TList.begin(); itr != TList.end(); ++itr )
194 axpy(1.0, (*itr).xvals, center);
195 scal((double)(1.0/ns),center);
197 // Compute the relative deviations
198 for ( itr = TList.begin(); itr != TList.end(); ++itr ) {
199 for (k = 0; k<n; k++) {
201 dispers(k)=dispers(k)+pow(center(k)-x(k),2.0);
204 scal((double)(1.0/ns),dispers);
208 for (k=1; k<n; k++) {
209 if (dispers(k)>tmp) {
213 B1.ub(i)=center(i) ; B2.lb(i)=center(i);
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);
225 fm2=min(fm2,(*itr).objval);
228 // Set minf of B1 and B2
229 B1.minf=fm1 ; B2.minf=fm2;
232 void TBox::dispTrials() {
233 // Display information about box
235 copy(TList.begin(), TList.end(), ostream_iterator<Trial,char>(cout));
237 copy(TList.begin(), TList.end(), ostream_iterator<Trial>(cout));
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;
249 bool TBox::InsideBox(RCRVector x) {
250 // Returns TRUE if the point X lies inside BOX, FALSE otherwise
252 for (int i=0 ; i<n ; i++)
253 if (x(i)<lb(i) || x(i)>ub(i)) return FALSE;
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
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))
268 if (x(i)<domain.lb(i) || x(i)>domain.ub(i)) {
272 if (ins_box==1 && ins_dom==1)
274 if (ins_box==0 && ins_dom==1)
276 if (ins_box==0 && ins_dom==0)
279 // Something has gone wrong!
280 cout << "Error in OutsideBox, exiting\n";
286 double TBox::ShortestSide(int *idx) {
287 // Returns the shortest side of the box and it's index
289 double tmp=ub(0)-lb(0);
290 for (int i=1 ; i<n ;i++)
291 if ( (ub(i)-lb(i))<tmp ) {
299 double TBox::LongestSide(int *idx) {
300 // Returns the longest side of the box and it's index
302 double tmp=ub(0)-lb(0);
303 for (int i=1 ; i<n ;i++)
304 if ( (ub(i)-lb(i))>tmp ) {
312 double TBox::ClosestSide(RCRVector x) {
313 // Returns the shortest distance from point X to the box B.
315 // Warning: The output of this functon is nonsense if the
316 // point X lies outside B. Should we try to detect this case?
321 for (int i=0 ; i<n ; i++) {
322 tmp = min( x(i)-lb(i), ub(i)-x(i) );
323 dist = min(dist, tmp);
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)
335 for (int i=0 ; i<n ; i++) {
336 tmp = max( x(i)-lb(i), ub(i)-x(i) );
337 dist = max(dist, tmp);
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
349 // Due to round of errors the algorithm can fail to find an intersection
350 // The caller is notified and should act accordingly
352 // The routine returns FALSE if no intersection was found, TRUE otherwise
361 while (i<n && done==FALSE) {
366 for (k=1; k<=2; k++) {
371 gamma=(alpha-x(i))/h(i);
374 for (j=0; j<n; j++) {
376 z(j)=x(j)+gamma*h(j);
377 if (z(j)<lb(j) || z(j)>ub(j)) {
383 copy(z,tmpV); axpy(-1.0,x,tmpV); // tmpV=z-x
384 if (isect==1 && dot(tmpV,h)>0) {
393 double TBox::LowerBound(double maxgrad) {
394 // Lower bound estimation
397 list<Trial>::const_iterator itr1,itr2 ;
400 RVector x1(n), x2(n) ;
402 for ( itr1 = TList.begin(); itr1 != TList.end(); ++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)) ;
410 // cout << "est=" << est << " x1=" << x1 << " x2=" << x2 << endl ;
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;