1 /* DIRect.f -- translated by f2c (version 20050501).
3 f2c output hand-cleaned by SGJ (August 2007).
7 #include "direct-internal.h"
9 /* Common Block Declarations */
11 /* Table of constant values */
13 /* +-----------------------------------------------------------------------+ */
14 /* | Program : Direct.f | */
15 /* | Last modified : 07-16-2001 | */
16 /* | Written by : Joerg Gablonsky (jmgablon@unity.ncsu.edu) | */
17 /* | North Carolina State University | */
18 /* | Dept. of Mathematics | */
19 /* | DIRECT is a method to solve problems of the form: | */
20 /* | min f: Q --> R, | */
21 /* | where f is the function to be minimized and Q is an n-dimensional | */
22 /* | hyperrectangle given by the the following equation: | */
24 /* | Q={ x : l(i) <= x(i) <= u(i), i = 1,...,n }. | */
25 /* | Note: This version of DIRECT can also handle hidden constraints. By | */
26 /* | this we mean that the function may not be defined over the whole| */
27 /* | hyperrectangle Q, but only over a subset, which is not given | */
28 /* | analytically. | */
30 /* | We now give a brief outline of the algorithm: | */
32 /* | The algorithm starts with mapping the hyperrectangle Q to the | */
33 /* | n-dimensional unit hypercube. DIRECT then samples the function at | */
34 /* | the center of this hypercube and at 2n more points, 2 in each | */
35 /* | coordinate direction. Uisng these function values, DIRECT then | */
36 /* | divides the domain into hyperrectangles, each having exactly one of | */
37 /* | the sampling points as its center. In each iteration, DIRECT chooses| */
38 /* | some of the existing hyperrectangles to be further divided. | */
39 /* | We provide two different strategies of how to decide which | */
40 /* | hyperrectangles DIRECT divides and several different convergence | */
43 /* | DIRECT was designed to solve problems where the function f is | */
44 /* | Lipschitz continues. However, DIRECT has proven to be effective on | */
45 /* | more complex problems than these. | */
46 /* +-----------------------------------------------------------------------+ */
47 /* Subroutine */ void direct_direct_(fp fcn, doublereal *x, integer *n, doublereal *eps, doublereal epsabs, integer *maxf, integer *maxt, doublereal *fmin, doublereal *l,
48 doublereal *u, integer *algmethod, integer *ierror, FILE *logfile,
49 doublereal *fglobal, doublereal *fglper, doublereal *volper,
50 doublereal *sigmaper, void *fcn_data)
52 /* System generated locals */
56 /* constants (FIXME: change to variable?) */
57 const integer MAXFUNC = 90000;
58 const integer MAXDEEP = 600;
59 const integer MAXDIV = 3000;
63 doublereal *c__ = 0 /* was [90000][64] */, *f = 0 /*
65 integer i__, j, *s = 0 /* was [3000][2] */, t;
68 integer ifeasiblef, iepschange, actmaxdeep;
69 integer actdeep_div__, iinfesiblef;
70 integer pos1, newtosample;
72 doublereal *oldl = 0, fmax;
74 doublereal kmax, *oldu = 0;
75 integer oops, *list2 = 0 /* was [64][2] */, cheat;
77 integer mdeep, *point = 0, start;
78 integer *anchor = 0, *length = 0 /* was [90000][64] */, *arrayi = 0;
79 doublereal *levels = 0, *thirds = 0;
82 integer oldpos, minpos, maxpos, tstart, actdeep, ifreeold, oldmaxf;
83 integer numfunc, version;
86 /* FIXME: change sizes dynamically? */
87 #define MY_ALLOC(p, t, n) p = (t *) malloc(sizeof(t) * (n)); \
88 if (!(p)) { *ierror = -100; goto cleanup; }
89 MY_ALLOC(c__, doublereal, MAXFUNC * (*n));
90 MY_ALLOC(f, doublereal, MAXFUNC * 2);
91 MY_ALLOC(s, integer, MAXDIV * 2);
92 MY_ALLOC(w, doublereal, (*n));
93 MY_ALLOC(oldl, doublereal, (*n));
94 MY_ALLOC(oldu, doublereal, (*n));
95 MY_ALLOC(list2, integer, (*n) * 2);
96 MY_ALLOC(point, integer, MAXFUNC);
97 MY_ALLOC(anchor, integer, MAXDEEP + 2);
98 MY_ALLOC(length, integer, MAXFUNC * (*n));
99 MY_ALLOC(arrayi, integer, (*n));
100 MY_ALLOC(levels, doublereal, MAXDEEP + 1);
101 MY_ALLOC(thirds, doublereal, MAXDEEP + 1);
103 /* +-----------------------------------------------------------------------+ */
104 /* | SUBROUTINE Direct | */
106 /* | fcn -- The argument containing the name of the user-supplied | */
107 /* | SUBROUTINE that returns values for the function to be | */
109 /* | n -- The dimension of the problem. | */
110 /* | eps -- Exceeding value. If eps > 0, we use the same epsilon for | */
111 /* | all iterations. If eps < 0, we use the update formula from | */
113 /* | eps = max(1.D-4*abs(fmin),epsfix), | */
114 /* | where epsfix = abs(eps), the absolute value of eps which is| */
115 /* | passed to the function. | */
116 /* | maxf -- The maximum number of function evaluations. | */
117 /* | maxT -- The maximum number of iterations. | */
118 /* | Direct stops when either the maximum number of iterations | */
119 /* | is reached or more than maxf function-evalutions were made.| */
120 /* | l -- The lower bounds of the hyperbox. | */
121 /* | u -- The upper bounds of the hyperbox. | */
122 /* |algmethod-- Choose the method, that is either use the original method | */
123 /* | as described by Jones et.al. (0) or use our modification(1)| */
124 /* | logfile -- File-Handle for the logfile. DIRECT expects this file to be| */
125 /* | opened and closed by the user outside of DIRECT. We moved | */
126 /* | this to the outside so the user can add extra informations | */
127 /* | to this file before and after the call to DIRECT. | */
128 /* | fglobal -- Function value of the global optimum. If this value is not | */
129 /* | known (that is, we solve a real problem, not a testproblem)| */
130 /* | set this value to -1.D100 and fglper (see below) to 0.D0. | */
131 /* | fglper -- Terminate the optimization when the percent error | */
132 /* | 100(f_min - fglobal)/max(1,abs(fglobal)) < fglper. | */
133 /* | volper -- Terminate the optimization when the volume of the | */
134 /* | hyperrectangle S with f(c(S)) = fmin is less then volper | */
135 /* | percent of the volume of the original hyperrectangle. | */
136 /* |sigmaper -- Terminate the optimization when the measure of the | */
137 /* | hyperrectangle S with f(c(S)) = fmin is less then sigmaper.| */
139 /* | User data that is passed through without being changed: | */
140 /* | fcn_data - opaque pointer to any user data | */
144 /* | x -- The final point obtained in the optimization process. | */
145 /* | X should be a good approximation to the global minimum | */
146 /* | for the function within the hyper-box. | */
148 /* | fmin -- The value of the function at x. | */
149 /* | Ierror -- Error flag. If Ierror is lower 0, an error has occured. The| */
150 /* | values of Ierror mean | */
151 /* | Fatal errors : | */
152 /* | -1 u(i) <= l(i) for some i. | */
153 /* | -2 maxf is too large. | */
154 /* | -3 Initialization in DIRpreprc failed. | */
155 /* | -4 Error in DIRSamplepoints, that is there was an error | */
156 /* | in the creation of the sample points. | */
157 /* | -5 Error in DIRSamplef, that is an error occured while | */
158 /* | the function was sampled. | */
159 /* | -6 Error in DIRDoubleInsert, that is an error occured | */
160 /* | DIRECT tried to add all hyperrectangles with the same| */
161 /* | size and function value at the center. Either | */
162 /* | increase maxdiv or use our modification (Jones = 1). | */
163 /* | Termination values : | */
164 /* | 1 Number of function evaluations done is larger then | */
166 /* | 2 Number of iterations is equal to maxT. | */
167 /* | 3 The best function value found is within fglper of | */
168 /* | the (known) global optimum, that is | */
169 /* | 100(fmin - fglobal/max(1,|fglobal|)) < fglper. | */
170 /* | Note that this termination signal only occurs when | */
171 /* | the global optimal value is known, that is, a test | */
172 /* | function is optimized. | */
173 /* | 4 The volume of the hyperrectangle with fmin at its | */
174 /* | center is less than volper percent of the volume of | */
175 /* | the original hyperrectangle. | */
176 /* | 5 The measure of the hyperrectangle with fmin at its | */
177 /* | center is less than sigmaper. | */
179 /* | SUBROUTINEs used : | */
181 /* | DIRheader, DIRInitSpecific, DIRInitList, DIRpreprc, DIRInit, DIRChoose| */
182 /* | DIRDoubleInsert, DIRGet_I, DIRSamplepoints, DIRSamplef, DIRDivide | */
183 /* | DIRInsertList, DIRreplaceInf, DIRWritehistbox, DIRsummary, Findareas | */
185 /* | Functions used : | */
187 /* | DIRgetMaxdeep, DIRgetlevel | */
188 /* +-----------------------------------------------------------------------+ */
189 /* +-----------------------------------------------------------------------+ */
191 /* +-----------------------------------------------------------------------+ */
192 /* +-----------------------------------------------------------------------+ */
193 /* | The maximum of function evaluations allowed. | */
194 /* | The maximum dept of the algorithm. | */
195 /* | The maximum number of divisions allowed. | */
196 /* | The maximal dimension of the problem. | */
197 /* +-----------------------------------------------------------------------+ */
198 /* +-----------------------------------------------------------------------+ */
199 /* | Global Variables. | */
200 /* +-----------------------------------------------------------------------+ */
201 /* +-----------------------------------------------------------------------+ */
202 /* | EXTERNAL Variables. | */
203 /* +-----------------------------------------------------------------------+ */
204 /* +-----------------------------------------------------------------------+ */
205 /* | User Variables. | */
206 /* | These can be used to pass user defined data to the function to be | */
208 /* +-----------------------------------------------------------------------+ */
209 /* +-----------------------------------------------------------------------+ */
210 /* | Place to define, if needed, some application-specific variables. | */
211 /* | Note: You should try to use the arrays defined above for this. | */
212 /* +-----------------------------------------------------------------------+ */
213 /* +-----------------------------------------------------------------------+ */
214 /* | End of application - specific variables ! | */
215 /* +-----------------------------------------------------------------------+ */
216 /* +-----------------------------------------------------------------------+ */
217 /* | Internal variables : | */
218 /* | f -- values of functions. | */
219 /* |divfactor-- Factor used for termination with known global minimum. | */
220 /* | anchor -- anchors of lists with deepness i, -1 is anchor for list of | */
221 /* | NaN - values. | */
222 /* | S -- List of potentially optimal points. | */
223 /* | point -- lists | */
224 /* | ifree -- first free position | */
225 /* | c -- midpoints of arrays | */
226 /* | thirds -- Precalculated values of 1/3^i. | */
227 /* | levels -- Length of intervals. | */
228 /* | length -- Length of intervall (index) | */
229 /* | t -- actual iteration | */
230 /* | j -- loop-variable | */
231 /* | actdeep -- the actual minimal interval-length index | */
232 /* | Minpos -- position of the actual minimum | */
233 /* | file -- The filehandle for a datafile. | */
234 /* | maxpos -- The number of intervalls, which are truncated. | */
235 /* | help -- A help variable. | */
236 /* | numfunc -- The actual number of function evaluations. | */
237 /* | file2 -- The filehandle for an other datafile. | */
238 /* | ArrayI -- Array with the indexes of the sides with maximum length. | */
239 /* | maxi -- Number of directions with maximal side length. | */
240 /* | oops -- Flag which shows if anything went wrong in the | */
241 /* | initialisation. | */
242 /* | cheat -- Obsolete. If equal 1, we don't allow Ktilde > kmax. | */
243 /* | writed -- If writed=1, store final division to plot with Matlab. | */
244 /* | List2 -- List of indicies of intervalls, which are to be truncated. | */
245 /* | i -- Another loop-variable. | */
246 /* |actmaxdeep-- The actual maximum (minimum) of possible Interval length. | */
247 /* | oldpos -- The old index of the minimum. Used to print only, if there | */
248 /* | is a new minimum found. | */
249 /* | tstart -- The start of the outer loop. | */
250 /* | start -- The postion of the starting point in the inner loop. | */
251 /* | Newtosample -- The total number of points to sample in the inner loop.| */
252 /* | w -- Array used to divide the intervalls | */
253 /* | kmax -- Obsolete. If cheat = 1, Ktilde was not allowed to be larger| */
254 /* | than kmax. If Ktilde > kmax, we set ktilde = kmax. | */
255 /* | delta -- The distance to new points from center of old hyperrec. | */
256 /* | pos1 -- Help variable used as an index. | */
257 /* | version -- Store the version number of DIRECT. | */
258 /* | oldmaxf -- Store the original function budget. | */
259 /* |increase -- Flag used to keep track if function budget was increased | */
260 /* | because no feasible point was found. | */
261 /* | ifreeold -- Keep track which index was free before. Used with | */
262 /* | SUBROUTINE DIRReplaceInf. | */
263 /* |actdeep_div-- Keep track of the current depths for divisions. | */
264 /* | oldl -- Array used to store the original bounds of the domain. | */
265 /* | oldu -- Array used to store the original bounds of the domain. | */
266 /* | epsfix -- If eps < 0, we use Jones update formula. epsfix stores the | */
267 /* | absolute value of epsilon. | */
268 /* |iepschange-- flag iepschange to store if epsilon stays fixed or is | */
270 /* |DIRgetMaxdeep-- Function to calculate the level of a hyperrectangle. | */
271 /* |DIRgetlevel-- Function to calculate the level and stage of a hyperrec. | */
272 /* | fmax -- Keep track of the maximum value of the function found. | */
273 /* |Ifeasiblef-- Keep track if a feasible point has been found so far. | */
274 /* | Ifeasiblef = 0 means a feasible point has been found, | */
275 /* | Ifeasiblef = 1 no feasible point has been found. | */
276 /* +-----------------------------------------------------------------------+ */
277 /* +-----------------------------------------------------------------------+ */
278 /* | JG 09/25/00 Version counter. | */
279 /* +-----------------------------------------------------------------------+ */
280 /* +-----------------------------------------------------------------------+ */
281 /* | JG 09/24/00 Add another actdeep to keep track of the current depths | */
282 /* | for divisions. | */
283 /* +-----------------------------------------------------------------------+ */
284 /* +-----------------------------------------------------------------------+ */
285 /* |JG 01/13/01 Added epsfix for epsilon update. If eps < 0, we use Jones | */
286 /* | update formula. epsfix stores the absolute value of epsilon| */
287 /* | then. Also added flag iepschange to store if epsilon stays | */
288 /* | fixed or is changed. | */
289 /* +-----------------------------------------------------------------------+ */
290 /* +-----------------------------------------------------------------------+ */
291 /* | JG 01/22/01 fmax is used to keep track of the maximum value found. | */
292 /* +-----------------------------------------------------------------------+ */
293 /* +-----------------------------------------------------------------------+ */
294 /* | JG 01/22/01 Ifeasiblef is used to keep track if a feasible point has | */
295 /* | been found so far. Ifeasiblef = 0 means a feasible point | */
296 /* | has been found, Ifeasiblef = 1 if not. | */
297 /* | JG 03/09/01 IInfeasible is used to keep track if an infeasible point | */
298 /* | has been found. IInfeasible > 0 means a infeasible point | */
299 /* | has been found, IInfeasible = 0 if not. | */
300 /* +-----------------------------------------------------------------------+ */
301 /* +-----------------------------------------------------------------------+ */
302 /* +-----------------------------------------------------------------------+ */
303 /* | Start of code. | */
304 /* +-----------------------------------------------------------------------+ */
305 /* +-----------------------------------------------------------------------+ */
306 /* Parameter adjustments */
314 /* +-----------------------------------------------------------------------+ */
315 /* | Save the upper and lower bounds. | */
316 /* +-----------------------------------------------------------------------+ */
318 for (i__ = 1; i__ <= i__1; ++i__) {
319 oldu[i__ - 1] = u[i__];
320 oldl[i__ - 1] = l[i__];
323 /* +-----------------------------------------------------------------------+ */
324 /* | Set version. | */
325 /* +-----------------------------------------------------------------------+ */
327 /* +-----------------------------------------------------------------------+ */
328 /* | Set parameters. | */
329 /* | If cheat > 0, we do not allow \tilde{K} to be larger than kmax, and| */
330 /* | set \tilde{K} to set value if necessary. Not used anymore. | */
331 /* +-----------------------------------------------------------------------+ */
335 /* +-----------------------------------------------------------------------+ */
336 /* | Write the header of the logfile. | */
337 /* +-----------------------------------------------------------------------+ */
338 direct_dirheader_(logfile, &version, &x[1], n, eps, maxf, maxt, &l[1], &u[1],
339 algmethod, &MAXFUNC, &MAXDEEP, fglobal, fglper, ierror, &epsfix, &
340 iepschange, volper, sigmaper);
341 /* +-----------------------------------------------------------------------+ */
342 /* | If an error has occured while writing the header (we do some checking | */
343 /* | of variables there), return to the main program. | */
344 /* +-----------------------------------------------------------------------+ */
348 /* +-----------------------------------------------------------------------+ */
349 /* | If the known global minimum is equal 0, we cannot divide by it. | */
350 /* | Therefore we set it to 1. If not, we set the divisionfactor to the | */
351 /* | absolute value of the global minimum. | */
352 /* +-----------------------------------------------------------------------+ */
353 if (*fglobal == 0.) {
356 divfactor = fabs(*fglobal);
358 /* +-----------------------------------------------------------------------+ */
359 /* | Save the budget given by the user. The variable maxf will be changed | */
360 /* | if in the beginning no feasible points are found. | */
361 /* +-----------------------------------------------------------------------+ */
364 /* +-----------------------------------------------------------------------+ */
365 /* | Initialiase the lists. | */
366 /* +-----------------------------------------------------------------------+ */
367 direct_dirinitlist_(anchor, &ifree, point, f, &MAXFUNC, &MAXDEEP);
368 /* +-----------------------------------------------------------------------+ */
369 /* | Call the routine to initialise the mapping of x from the n-dimensional| */
370 /* | unit cube to the hypercube given by u and l. If an error occured, | */
371 /* | give out a error message and return to the main program with the error| */
373 /* | JG 07/16/01 Changed call to remove unused data. | */
374 /* +-----------------------------------------------------------------------+ */
375 direct_dirpreprc_(&u[1], &l[1], n, &l[1], &u[1], &oops);
378 fprintf(logfile,"WARNING: Initialization in DIRpreprc failed.\n");
383 /* +-----------------------------------------------------------------------+ */
384 /* | Initialise the algorithm DIRECT. | */
385 /* +-----------------------------------------------------------------------+ */
386 /* +-----------------------------------------------------------------------+ */
387 /* | Added variable to keep track of the maximum value found. | */
388 /* +-----------------------------------------------------------------------+ */
389 direct_dirinit_(f, fcn, c__, length, &actdeep, point, anchor, &ifree,
390 logfile, arrayi, &maxi, list2, w, &x[1], &l[1], &u[1],
391 fmin, &minpos, thirds, levels, &MAXFUNC, &MAXDEEP, n, n, &
392 fmax, &ifeasiblef, &iinfesiblef, ierror, fcn_data, jones);
393 /* +-----------------------------------------------------------------------+ */
394 /* | Added error checking. | */
395 /* +-----------------------------------------------------------------------+ */
399 fprintf(logfile, "WARNING: Error occured in routine DIRsamplepoints.\n");
404 fprintf(logfile, "WARNING: Error occured in routine DIRsamplef..\n");
408 numfunc = maxi + 1 + maxi;
412 /* +-----------------------------------------------------------------------+ */
413 /* | If no feasible point has been found, give out the iteration, the | */
414 /* | number of function evaluations and a warning. Otherwise, give out | */
415 /* | the iteration, the number of function evaluations done and fmin. | */
416 /* +-----------------------------------------------------------------------+ */
417 if (ifeasiblef > 0) {
419 fprintf(logfile, "No feasible point found in %d iterations "
420 "and %d function evaluations.\n", tstart-1, numfunc);
423 fprintf(logfile, "%d, %d, %g, %g\n",
424 tstart-1, numfunc, *fmin, fmax);
426 /* +-----------------------------------------------------------------------+ */
427 /* +-----------------------------------------------------------------------+ */
429 /* +-----------------------------------------------------------------------+ */
430 /* +-----------------------------------------------------------------------+ */
432 for (t = tstart; t <= i__1; ++t) {
433 /* +-----------------------------------------------------------------------+ */
434 /* | Choose the sample points. The indices of the sample points are stored | */
435 /* | in the list S. | */
436 /* +-----------------------------------------------------------------------+ */
437 actdeep = actmaxdeep;
438 direct_dirchoose_(anchor, s, &MAXDEEP, f, fmin, *eps, epsabs, levels, &maxpos, length,
439 &MAXFUNC, &MAXDEEP, &MAXDIV, n, logfile, &cheat, &
440 kmax, &ifeasiblef, jones);
441 /* +-----------------------------------------------------------------------+ */
442 /* | Add other hyperrectangles to S, which have the same level and the same| */
443 /* | function value at the center as the ones found above (that are stored | */
444 /* | in S). This is only done if we use the original DIRECT algorithm. | */
445 /* | JG 07/16/01 Added Errorflag. | */
446 /* +-----------------------------------------------------------------------+ */
447 if (*algmethod == 0) {
448 direct_dirdoubleinsert_(anchor, s, &maxpos, point, f, &MAXDEEP, &MAXFUNC,
453 "WARNING: Capacity of array S in DIRDoubleInsert reached. Increase maxdiv.\n"
454 "This means that there are a lot of hyperrectangles with the same function\n"
455 "value at the center. We suggest to use our modification instead (Jones = 1)\n"
461 /* +-----------------------------------------------------------------------+ */
462 /* | Initialise the number of sample points in this outer loop. | */
463 /* +-----------------------------------------------------------------------+ */
466 for (j = 1; j <= i__2; ++j) {
467 actdeep = s[j + 2999];
468 /* +-----------------------------------------------------------------------+ */
469 /* | If the actual index is a point to sample, do it. | */
470 /* +-----------------------------------------------------------------------+ */
472 /* +-----------------------------------------------------------------------+ */
473 /* | JG 09/24/00 Calculate the value delta used for sampling points. | */
474 /* +-----------------------------------------------------------------------+ */
475 actdeep_div__ = direct_dirgetmaxdeep_(&s[j - 1], length, &MAXFUNC,
477 delta = thirds[actdeep_div__ + 1];
478 actdeep = s[j + 2999];
479 /* +-----------------------------------------------------------------------+ */
480 /* | If the current dept of division is only one under the maximal allowed | */
481 /* | dept, stop the computation. | */
482 /* +-----------------------------------------------------------------------+ */
483 if (actdeep + 1 >= mdeep) {
485 fprintf(logfile, "WARNING: Maximum number of levels reached. Increase maxdeep.\n");
489 actmaxdeep = MAX(actdeep,actmaxdeep);
491 if (! (anchor[actdeep + 1] == help)) {
492 pos1 = anchor[actdeep + 1];
493 while(! (point[pos1 - 1] == help)) {
494 pos1 = point[pos1 - 1];
496 point[pos1 - 1] = point[help - 1];
498 anchor[actdeep + 1] = point[help - 1];
501 actdeep = (integer) f[help - 1];
503 /* +-----------------------------------------------------------------------+ */
504 /* | Get the Directions in which to decrease the intervall-length. | */
505 /* +-----------------------------------------------------------------------+ */
506 direct_dirget_i__(length, &help, arrayi, &maxi, n, &MAXFUNC);
507 /* +-----------------------------------------------------------------------+ */
508 /* | Sample the function. To do this, we first calculate the points where | */
509 /* | we need to sample the function. After checking for errors, we then do | */
510 /* | the actual evaluation of the function, again followed by checking for | */
512 /* +-----------------------------------------------------------------------+ */
513 direct_dirsamplepoints_(c__, arrayi, &delta, &help, &start, length,
514 logfile, f, &ifree, &maxi, point, &x[
515 1], &l[1], fmin, &minpos, &u[1], n, &MAXFUNC, &
519 fprintf(logfile, "WARNING: Error occured in routine DIRsamplepoints.\n");
524 /* +-----------------------------------------------------------------------+ */
525 /* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
526 /* +-----------------------------------------------------------------------+ */
527 direct_dirsamplef_(c__, arrayi, &delta, &help, &start, length,
528 logfile, f, &ifree, &maxi, point, fcn, &x[
529 1], &l[1], fmin, &minpos, &u[1], n, &MAXFUNC, &
530 MAXDEEP, &oops, &fmax, &ifeasiblef, &iinfesiblef,
534 fprintf(logfile, "WARNING: Error occured in routine DIRsamplef.\n");
538 /* +-----------------------------------------------------------------------+ */
539 /* | Divide the intervalls. | */
540 /* +-----------------------------------------------------------------------+ */
541 direct_dirdivide_(&start, &actdeep_div__, length, point, arrayi, &
542 help, list2, w, &maxi, f, &MAXFUNC, &MAXDEEP, n);
543 /* +-----------------------------------------------------------------------+ */
544 /* | Insert the new intervalls into the list (sorted). | */
545 /* +-----------------------------------------------------------------------+ */
546 direct_dirinsertlist_(&start, anchor, point, f, &maxi, length, &
547 MAXFUNC, &MAXDEEP, n, &help, jones);
548 /* +-----------------------------------------------------------------------+ */
549 /* | Increase the number of function evaluations. | */
550 /* +-----------------------------------------------------------------------+ */
551 numfunc = numfunc + maxi + maxi;
553 /* +-----------------------------------------------------------------------+ */
554 /* | End of main loop. | */
555 /* +-----------------------------------------------------------------------+ */
558 /* +-----------------------------------------------------------------------+ */
559 /* | If there is a new minimum, show the actual iteration, the number of | */
560 /* | function evaluations, the minimum value of f (so far) and the position| */
561 /* | in the array. | */
562 /* +-----------------------------------------------------------------------+ */
563 if (oldpos < minpos) {
565 fprintf(logfile, "%d, %d, %g, %g\n",
566 t, numfunc, *fmin, fmax);
568 /* +-----------------------------------------------------------------------+ */
569 /* | If no feasible point has been found, give out the iteration, the | */
570 /* | number of function evaluations and a warning. | */
571 /* +-----------------------------------------------------------------------+ */
572 if (ifeasiblef > 0) {
574 fprintf(logfile, "No feasible point found in %d iterations "
575 "and %d function evaluations\n", t, numfunc);
577 /* +-----------------------------------------------------------------------+ */
578 /* +-----------------------------------------------------------------------+ */
579 /* | Termination Checks | */
580 /* +-----------------------------------------------------------------------+ */
581 /* +-----------------------------------------------------------------------+ */
582 /* | JG 01/22/01 Calculate the index for the hyperrectangle at which | */
583 /* | fmin is assumed. We then calculate the volume of this | */
584 /* | hyperrectangle and store it in delta. This delta can be | */
585 /* | used to stop DIRECT once the volume is below a certain | */
586 /* | percentage of the original volume. Since the original | */
587 /* | is 1 (scaled), we can stop once delta is below a certain | */
588 /* | percentage, given by volper. | */
589 /* +-----------------------------------------------------------------------+ */
592 actdeep_div__ = direct_dirgetlevel_(&minpos, length, &MAXFUNC, n, jones);
594 /* +-----------------------------------------------------------------------+ */
595 /* | JG 07/16/01 Use precalculated values to calculate volume. | */
596 /* +-----------------------------------------------------------------------+ */
597 delta = thirds[actdeep_div__] * 100;
598 if (delta <= *volper) {
601 fprintf(logfile, "DIRECT stopped: Volume of S_min is "
602 "%g%% < %g%% of the original volume.\n",
606 /* +-----------------------------------------------------------------------+ */
607 /* | JG 01/23/01 Calculate the measure for the hyperrectangle at which | */
608 /* | fmin is assumed. If this measure is smaller then sigmaper,| */
609 /* | we stop DIRECT. | */
610 /* +-----------------------------------------------------------------------+ */
611 actdeep_div__ = direct_dirgetlevel_(&minpos, length, &MAXFUNC, n, jones);
612 delta = levels[actdeep_div__];
613 if (delta <= *sigmaper) {
616 fprintf(logfile, "DIRECT stopped: Measure of S_min "
617 "= %g < %g.\n", delta, *sigmaper);
620 /* +-----------------------------------------------------------------------+ */
621 /* | If the best found function value is within fglper of the (known) | */
622 /* | global minimum value, terminate. This only makes sense if this optimal| */
623 /* | value is known, that is, in test problems. | */
624 /* +-----------------------------------------------------------------------+ */
625 if ((*fmin - *fglobal) * 100 / divfactor <= *fglper) {
628 fprintf(logfile, "DIRECT stopped: fmin within fglper of global minimum.\n");
631 /* +-----------------------------------------------------------------------+ */
632 /* | Find out if there are infeasible points which are near feasible ones. | */
633 /* | If this is the case, replace the function value at the center of the | */
634 /* | hyper rectangle by the lowest function value of a nearby function. | */
635 /* | If no infeasible points exist (IInfesiblef = 0), skip this. | */
636 /* +-----------------------------------------------------------------------+ */
637 if (iinfesiblef > 0) {
638 direct_dirreplaceinf_(&ifree, &ifreeold, f, c__, thirds, length, anchor,
639 point, &u[1], &l[1], &MAXFUNC, &MAXDEEP, n, n,
640 logfile, &fmax, jones);
643 /* +-----------------------------------------------------------------------+ */
644 /* | If iepschange = 1, we use the epsilon change formula from Jones. | */
645 /* +-----------------------------------------------------------------------+ */
646 if (iepschange == 1) {
648 d__1 = fabs(*fmin) * 1e-4;
649 *eps = MAX(d__1,epsfix);
651 /* +-----------------------------------------------------------------------+ */
652 /* | If no feasible point has been found yet, set the maximum number of | */
653 /* | function evaluations to the number of evaluations already done plus | */
654 /* | the budget given by the user. | */
655 /* | If the budget has already be increased, increase it again. If a | */
656 /* | feasible point has been found, remark that and reset flag. No further | */
657 /* | increase is needed. | */
658 /* +-----------------------------------------------------------------------+ */
660 *maxf = numfunc + oldmaxf;
661 if (ifeasiblef == 0) {
663 fprintf(logfile, "DIRECT found a feasible point. The "
664 "adjusted budget is now set to %d.\n", *maxf);
668 /* +-----------------------------------------------------------------------+ */
669 /* | Check if the number of function evaluations done is larger than the | */
670 /* | allocated budget. If this is the case, check if a feasible point was | */
671 /* | found. If this is a case, terminate. If no feasible point was found, | */
672 /* | increase the budget and set flag increase. | */
673 /* +-----------------------------------------------------------------------+ */
674 if (numfunc > *maxf) {
675 if (ifeasiblef == 0) {
678 fprintf(logfile, "DIRECT stopped: numfunc >= maxf.\n");
684 "DIRECT could not find a feasible point after %d function evaluations.\n"
685 "DIRECT continues until a feasible point is found.\n", numfunc);
686 *maxf = numfunc + oldmaxf;
691 /* +-----------------------------------------------------------------------+ */
692 /* +-----------------------------------------------------------------------+ */
693 /* | End of main loop. | */
694 /* +-----------------------------------------------------------------------+ */
695 /* +-----------------------------------------------------------------------+ */
696 /* +-----------------------------------------------------------------------+ */
697 /* | The algorithm stopped after maxT iterations. | */
698 /* +-----------------------------------------------------------------------+ */
701 fprintf(logfile, "DIRECT stopped: maxT iterations.\n");
704 /* +-----------------------------------------------------------------------+ */
705 /* | Store the position of the minimum in x. | */
706 /* +-----------------------------------------------------------------------+ */
708 for (i__ = 1; i__ <= i__1; ++i__) {
709 x[i__] = c__[minpos + i__ * 90000 - 90001] * l[i__] + l[i__] * u[i__];
710 u[i__] = oldu[i__ - 1];
711 l[i__] = oldl[i__ - 1];
714 /* +-----------------------------------------------------------------------+ */
715 /* | Store the number of function evaluations in maxf. | */
716 /* +-----------------------------------------------------------------------+ */
718 /* +-----------------------------------------------------------------------+ */
719 /* | Give out a summary of the run. | */
720 /* +-----------------------------------------------------------------------+ */
721 direct_dirsummary_(logfile, &x[1], &l[1], &u[1], n, fmin, fglobal, &numfunc,
723 /* +-----------------------------------------------------------------------+ */
724 /* | Format statements. | */
725 /* +-----------------------------------------------------------------------+ */
728 #define MY_FREE(p) if (p) free(p)