--- /dev/null
+C conversion: Steven G. Johnson (stevenj@alum.mit.edu)
+Original Fortran code: Joerg.M.Gablonsky (jmgablon@mailandnews.com)
--- /dev/null
+This code is based on the DIRECT 2.0.4 Fortran code by Gablonsky et al. at
+ http://www4.ncsu.edu/~ctk/SOFTWARE/DIRECTv204.tar.gz
+The C version was initially converted via f2c and then cleaned up and
+reorganized by Steven G. Johnson (stevenj@alum.mit.edu), August 2007.
+
+******** Copyright and license for the original Fortran DIRECT code ********
+Copyright (c) 1999, 2000, 2001 North Carolina State University
+
+This program is distributed under the MIT License (see
+http://www.opensource.org/licenses/mit-license.php):
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
+of the Software, and to permit persons to whom the Software is furnished to do
+so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
--- /dev/null
+/* DIRect.f -- translated by f2c (version 20050501).
+
+ f2c output hand-cleaned by SGJ (August 2007).
+*/
+
+#include <math.h>
+#include "direct-internal.h"
+
+/* Common Block Declarations */
+
+/* Table of constant values */
+
+static integer c__90000 = 90000;
+static integer c__600 = 600;
+static integer c__64 = 64;
+static integer c__3000 = 3000;
+
+/* +-----------------------------------------------------------------------+ */
+/* | Program : Direct.f | */
+/* | Last modified : 07-16-2001 | */
+/* | Written by : Joerg Gablonsky (jmgablon@unity.ncsu.edu) | */
+/* | North Carolina State University | */
+/* | Dept. of Mathematics | */
+/* | DIRECT is a method to solve problems of the form: | */
+/* | min f: Q --> R, | */
+/* | where f is the function to be minimized and Q is an n-dimensional | */
+/* | hyperrectangle given by the the following equation: | */
+/* | | */
+/* | Q={ x : l(i) <= x(i) <= u(i), i = 1,...,n }. | */
+/* | Note: This version of DIRECT can also handle hidden constraints. By | */
+/* | this we mean that the function may not be defined over the whole| */
+/* | hyperrectangle Q, but only over a subset, which is not given | */
+/* | analytically. | */
+/* | | */
+/* | We now give a brief outline of the algorithm: | */
+/* | | */
+/* | The algorithm starts with mapping the hyperrectangle Q to the | */
+/* | n-dimensional unit hypercube. DIRECT then samples the function at | */
+/* | the center of this hypercube and at 2n more points, 2 in each | */
+/* | coordinate direction. Uisng these function values, DIRECT then | */
+/* | divides the domain into hyperrectangles, each having exactly one of | */
+/* | the sampling points as its center. In each iteration, DIRECT chooses| */
+/* | some of the existing hyperrectangles to be further divided. | */
+/* | We provide two different strategies of how to decide which | */
+/* | hyperrectangles DIRECT divides and several different convergence | */
+/* | criteria. | */
+/* | | */
+/* | DIRECT was designed to solve problems where the function f is | */
+/* | Lipschitz continues. However, DIRECT has proven to be effective on | */
+/* | more complex problems than these. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_direct_(fp fcn, doublereal *x, integer *n, doublereal *eps, doublereal epsabs, integer *maxf, integer *maxt, doublereal *fmin, doublereal *l,
+ doublereal *u, integer *algmethod, integer *ierror, FILE *logfile,
+ doublereal *fglobal, doublereal *fglper, doublereal *volper,
+ doublereal *sigmaper, void *fcn_data)
+{
+ /* System generated locals */
+ integer i__1, i__2;
+ doublereal d__1;
+
+ /* Local variables */
+ integer increase;
+ doublereal *c__ = 0 /* was [90000][64] */, *f = 0 /*
+ was [90000][2] */;
+ integer i__, j, *s = 0 /* was [3000][2] */, t;
+ doublereal *w = 0;
+ doublereal divfactor;
+ integer ifeasiblef, iepschange, actmaxdeep;
+ integer actdeep_div__, iinfesiblef;
+ integer pos1, newtosample;
+ integer ifree, help;
+ doublereal *oldl = 0, fmax;
+ integer maxi;
+ doublereal kmax, *oldu = 0;
+ integer oops, *list2 = 0 /* was [64][2] */, cheat;
+ doublereal delta;
+ integer mdeep, *point = 0, start;
+ integer *anchor = 0, *length = 0 /* was [90000][64] */, *arrayi = 0;
+ doublereal *levels = 0, *thirds = 0;
+ integer writed;
+ doublereal epsfix;
+ integer oldpos, minpos, maxpos, tstart, actdeep, ifreeold, oldmaxf;
+ integer numfunc, version;
+ integer jones;
+
+ /* FIXME: change sizes dynamically? */
+#define MY_ALLOC(p, t, n) p = (t *) malloc(sizeof(t) * (n)); \
+ if (!(p)) { *ierror = -100; goto cleanup; }
+ MY_ALLOC(c__, doublereal, 5760000);
+ MY_ALLOC(f, doublereal, 180000);
+ MY_ALLOC(s, integer, 6000);
+ MY_ALLOC(w, doublereal, 64);
+ MY_ALLOC(oldl, doublereal, 64);
+ MY_ALLOC(oldu, doublereal, 64);
+ MY_ALLOC(list2, integer, 128);
+ MY_ALLOC(point, integer, 90000);
+ MY_ALLOC(anchor, integer, 602);
+ MY_ALLOC(length, integer, 5760000);
+ MY_ALLOC(arrayi, integer, 64);
+ MY_ALLOC(levels, doublereal, 601);
+ MY_ALLOC(thirds, doublereal, 601);
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE Direct | */
+/* | On entry | */
+/* | fcn -- The argument containing the name of the user-supplied | */
+/* | SUBROUTINE that returns values for the function to be | */
+/* | minimized. | */
+/* | n -- The dimension of the problem. | */
+/* | eps -- Exceeding value. If eps > 0, we use the same epsilon for | */
+/* | all iterations. If eps < 0, we use the update formula from | */
+/* | Jones: | */
+/* | eps = max(1.D-4*abs(fmin),epsfix), | */
+/* | where epsfix = abs(eps), the absolute value of eps which is| */
+/* | passed to the function. | */
+/* | maxf -- The maximum number of function evaluations. | */
+/* | maxT -- The maximum number of iterations. | */
+/* | Direct stops when either the maximum number of iterations | */
+/* | is reached or more than maxf function-evalutions were made.| */
+/* | l -- The lower bounds of the hyperbox. | */
+/* | u -- The upper bounds of the hyperbox. | */
+/* |algmethod-- Choose the method, that is either use the original method | */
+/* | as described by Jones et.al. (0) or use our modification(1)| */
+/* | logfile -- File-Handle for the logfile. DIRECT expects this file to be| */
+/* | opened and closed by the user outside of DIRECT. We moved | */
+/* | this to the outside so the user can add extra informations | */
+/* | to this file before and after the call to DIRECT. | */
+/* | fglobal -- Function value of the global optimum. If this value is not | */
+/* | known (that is, we solve a real problem, not a testproblem)| */
+/* | set this value to -1.D100 and fglper (see below) to 0.D0. | */
+/* | fglper -- Terminate the optimization when the percent error | */
+/* | 100(f_min - fglobal)/max(1,abs(fglobal)) < fglper. | */
+/* | volper -- Terminate the optimization when the volume of the | */
+/* | hyperrectangle S with f(c(S)) = fmin is less then volper | */
+/* | percent of the volume of the original hyperrectangle. | */
+/* |sigmaper -- Terminate the optimization when the measure of the | */
+/* | hyperrectangle S with f(c(S)) = fmin is less then sigmaper.| */
+/* | | */
+/* | User data that is passed through without being changed: | */
+/* | fcn_data - opaque pointer to any user data | */
+/* | | */
+/* | On return | */
+/* | | */
+/* | x -- The final point obtained in the optimization process. | */
+/* | X should be a good approximation to the global minimum | */
+/* | for the function within the hyper-box. | */
+/* | | */
+/* | fmin -- The value of the function at x. | */
+/* | Ierror -- Error flag. If Ierror is lower 0, an error has occured. The| */
+/* | values of Ierror mean | */
+/* | Fatal errors : | */
+/* | -1 u(i) <= l(i) for some i. | */
+/* | -2 maxf is too large. | */
+/* | -3 Initialization in DIRpreprc failed. | */
+/* | -4 Error in DIRSamplepoints, that is there was an error | */
+/* | in the creation of the sample points. | */
+/* | -5 Error in DIRSamplef, that is an error occured while | */
+/* | the function was sampled. | */
+/* | -6 Error in DIRDoubleInsert, that is an error occured | */
+/* | DIRECT tried to add all hyperrectangles with the same| */
+/* | size and function value at the center. Either | */
+/* | increase maxdiv or use our modification (Jones = 1). | */
+/* | Termination values : | */
+/* | 1 Number of function evaluations done is larger then | */
+/* | maxf. | */
+/* | 2 Number of iterations is equal to maxT. | */
+/* | 3 The best function value found is within fglper of | */
+/* | the (known) global optimum, that is | */
+/* | 100(fmin - fglobal/max(1,|fglobal|)) < fglper. | */
+/* | Note that this termination signal only occurs when | */
+/* | the global optimal value is known, that is, a test | */
+/* | function is optimized. | */
+/* | 4 The volume of the hyperrectangle with fmin at its | */
+/* | center is less than volper percent of the volume of | */
+/* | the original hyperrectangle. | */
+/* | 5 The measure of the hyperrectangle with fmin at its | */
+/* | center is less than sigmaper. | */
+/* | | */
+/* | SUBROUTINEs used : | */
+/* | | */
+/* | DIRheader, DIRInitSpecific, DIRInitList, DIRpreprc, DIRInit, DIRChoose| */
+/* | DIRDoubleInsert, DIRGet_I, DIRSamplepoints, DIRSamplef, DIRDivide | */
+/* | DIRInsertList, DIRreplaceInf, DIRWritehistbox, DIRsummary, Findareas | */
+/* | | */
+/* | Functions used : | */
+/* | | */
+/* | DIRgetMaxdeep, DIRgetlevel | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Parameters | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | The maximum of function evaluations allowed. | */
+/* | The maximum dept of the algorithm. | */
+/* | The maximum number of divisions allowed. | */
+/* | The maximal dimension of the problem. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Global Variables. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | EXTERNAL Variables. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | User Variables. | */
+/* | These can be used to pass user defined data to the function to be | */
+/* | optimized. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Place to define, if needed, some application-specific variables. | */
+/* | Note: You should try to use the arrays defined above for this. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | End of application - specific variables ! | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Internal variables : | */
+/* | f -- values of functions. | */
+/* |divfactor-- Factor used for termination with known global minimum. | */
+/* | anchor -- anchors of lists with deepness i, -1 is anchor for list of | */
+/* | NaN - values. | */
+/* | S -- List of potentially optimal points. | */
+/* | point -- lists | */
+/* | ifree -- first free position | */
+/* | c -- midpoints of arrays | */
+/* | thirds -- Precalculated values of 1/3^i. | */
+/* | levels -- Length of intervals. | */
+/* | length -- Length of intervall (index) | */
+/* | t -- actual iteration | */
+/* | j -- loop-variable | */
+/* | actdeep -- the actual minimal interval-length index | */
+/* | Minpos -- position of the actual minimum | */
+/* | file -- The filehandle for a datafile. | */
+/* | maxpos -- The number of intervalls, which are truncated. | */
+/* | help -- A help variable. | */
+/* | numfunc -- The actual number of function evaluations. | */
+/* | file2 -- The filehandle for an other datafile. | */
+/* | ArrayI -- Array with the indexes of the sides with maximum length. | */
+/* | maxi -- Number of directions with maximal side length. | */
+/* | oops -- Flag which shows if anything went wrong in the | */
+/* | initialisation. | */
+/* | cheat -- Obsolete. If equal 1, we don't allow Ktilde > kmax. | */
+/* | writed -- If writed=1, store final division to plot with Matlab. | */
+/* | List2 -- List of indicies of intervalls, which are to be truncated. | */
+/* | i -- Another loop-variable. | */
+/* |actmaxdeep-- The actual maximum (minimum) of possible Interval length. | */
+/* | oldpos -- The old index of the minimum. Used to print only, if there | */
+/* | is a new minimum found. | */
+/* | tstart -- The start of the outer loop. | */
+/* | start -- The postion of the starting point in the inner loop. | */
+/* | Newtosample -- The total number of points to sample in the inner loop.| */
+/* | w -- Array used to divide the intervalls | */
+/* | kmax -- Obsolete. If cheat = 1, Ktilde was not allowed to be larger| */
+/* | than kmax. If Ktilde > kmax, we set ktilde = kmax. | */
+/* | delta -- The distance to new points from center of old hyperrec. | */
+/* | pos1 -- Help variable used as an index. | */
+/* | version -- Store the version number of DIRECT. | */
+/* | oldmaxf -- Store the original function budget. | */
+/* |increase -- Flag used to keep track if function budget was increased | */
+/* | because no feasible point was found. | */
+/* | ifreeold -- Keep track which index was free before. Used with | */
+/* | SUBROUTINE DIRReplaceInf. | */
+/* |actdeep_div-- Keep track of the current depths for divisions. | */
+/* | oldl -- Array used to store the original bounds of the domain. | */
+/* | oldu -- Array used to store the original bounds of the domain. | */
+/* | epsfix -- If eps < 0, we use Jones update formula. epsfix stores the | */
+/* | absolute value of epsilon. | */
+/* |iepschange-- flag iepschange to store if epsilon stays fixed or is | */
+/* | changed. | */
+/* |DIRgetMaxdeep-- Function to calculate the level of a hyperrectangle. | */
+/* |DIRgetlevel-- Function to calculate the level and stage of a hyperrec. | */
+/* | fmax -- Keep track of the maximum value of the function found. | */
+/* |Ifeasiblef-- Keep track if a feasible point has been found so far. | */
+/* | Ifeasiblef = 0 means a feasible point has been found, | */
+/* | Ifeasiblef = 1 no feasible point has been found. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 09/25/00 Version counter. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 09/24/00 Add another actdeep to keep track of the current depths | */
+/* | for divisions. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* |JG 01/13/01 Added epsfix for epsilon update. If eps < 0, we use Jones | */
+/* | update formula. epsfix stores the absolute value of epsilon| */
+/* | then. Also added flag iepschange to store if epsilon stays | */
+/* | fixed or is changed. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 fmax is used to keep track of the maximum value found. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Ifeasiblef is used to keep track if a feasible point has | */
+/* | been found so far. Ifeasiblef = 0 means a feasible point | */
+/* | has been found, Ifeasiblef = 1 if not. | */
+/* | JG 03/09/01 IInfeasible is used to keep track if an infeasible point | */
+/* | has been found. IInfeasible > 0 means a infeasible point | */
+/* | has been found, IInfeasible = 0 if not. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Start of code. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+ /* Parameter adjustments */
+ --u;
+ --l;
+ --x;
+
+ /* Function Body */
+ writed = 0;
+ jones = *algmethod;
+/* +-----------------------------------------------------------------------+ */
+/* | Save the upper and lower bounds. | */
+/* +-----------------------------------------------------------------------+ */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ oldu[i__ - 1] = u[i__];
+ oldl[i__ - 1] = l[i__];
+/* L150: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Set version. | */
+/* +-----------------------------------------------------------------------+ */
+ version = 204;
+/* +-----------------------------------------------------------------------+ */
+/* | Set parameters. | */
+/* | If cheat > 0, we do not allow \tilde{K} to be larger than kmax, and| */
+/* | set \tilde{K} to set value if necessary. Not used anymore. | */
+/* +-----------------------------------------------------------------------+ */
+ cheat = 0;
+ kmax = 1e10;
+ mdeep = 600;
+/* +-----------------------------------------------------------------------+ */
+/* | Write the header of the logfile. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirheader_(logfile, &version, &x[1], n, eps, maxf, maxt, &l[1], &u[1],
+ algmethod, &c__90000, &c__600, fglobal, fglper, ierror, &epsfix, &
+ iepschange, volper, sigmaper);
+/* +-----------------------------------------------------------------------+ */
+/* | If an error has occured while writing the header (we do some checking | */
+/* | of variables there), return to the main program. | */
+/* +-----------------------------------------------------------------------+ */
+ if (*ierror < 0) {
+ goto cleanup;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | If the known global minimum is equal 0, we cannot divide by it. | */
+/* | Therefore we set it to 1. If not, we set the divisionfactor to the | */
+/* | absolute value of the global minimum. | */
+/* +-----------------------------------------------------------------------+ */
+ if (*fglobal == 0.) {
+ divfactor = 1.;
+ } else {
+ divfactor = fabs(*fglobal);
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Save the budget given by the user. The variable maxf will be changed | */
+/* | if in the beginning no feasible points are found. | */
+/* +-----------------------------------------------------------------------+ */
+ oldmaxf = *maxf;
+ increase = 0;
+/* +-----------------------------------------------------------------------+ */
+/* | Initialiase the lists. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirinitlist_(anchor, &ifree, point, f, &c__90000, &c__600);
+/* +-----------------------------------------------------------------------+ */
+/* | Call the routine to initialise the mapping of x from the n-dimensional| */
+/* | unit cube to the hypercube given by u and l. If an error occured, | */
+/* | give out a error message and return to the main program with the error| */
+/* | flag set. | */
+/* | JG 07/16/01 Changed call to remove unused data. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirpreprc_(&u[1], &l[1], n, &l[1], &u[1], &oops);
+ if (oops > 0) {
+ if (logfile)
+ fprintf(logfile,"WARNING: Initialization in DIRpreprc failed.\n");
+ *ierror = -3;
+ goto cleanup;
+ }
+ tstart = 2;
+/* +-----------------------------------------------------------------------+ */
+/* | Initialise the algorithm DIRECT. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Added variable to keep track of the maximum value found. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirinit_(f, fcn, c__, length, &actdeep, point, anchor, &ifree,
+ logfile, arrayi, &maxi, list2, w, &x[1], &l[1], &u[1],
+ fmin, &minpos, thirds, levels, &c__90000, &c__600, n, &c__64, &
+ fmax, &ifeasiblef, &iinfesiblef, ierror, fcn_data, jones);
+/* +-----------------------------------------------------------------------+ */
+/* | Added error checking. | */
+/* +-----------------------------------------------------------------------+ */
+ if (*ierror < 0) {
+ if (*ierror == -4) {
+ if (logfile)
+ fprintf(logfile, "WARNING: Error occured in routine DIRsamplepoints.\n");
+ goto cleanup;
+ }
+ if (*ierror == -5) {
+ if (logfile)
+ fprintf(logfile, "WARNING: Error occured in routine DIRsamplef..\n");
+ goto cleanup;
+ }
+ }
+ numfunc = maxi + 1 + maxi;
+ actmaxdeep = 1;
+ oldpos = 0;
+ tstart = 2;
+/* +-----------------------------------------------------------------------+ */
+/* | If no feasible point has been found, give out the iteration, the | */
+/* | number of function evaluations and a warning. Otherwise, give out | */
+/* | the iteration, the number of function evaluations done and fmin. | */
+/* +-----------------------------------------------------------------------+ */
+ if (ifeasiblef > 0) {
+ if (logfile)
+ fprintf(logfile, "No feasible point found in %d iterations "
+ "and %d function evaluations.\n", tstart-1, numfunc);
+ } else {
+ if (logfile)
+ fprintf(logfile, "%d, %d, %g, %g\n",
+ tstart-1, numfunc, *fmin, fmax);
+ }
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Main loop! | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+ i__1 = *maxt;
+ for (t = tstart; t <= i__1; ++t) {
+/* +-----------------------------------------------------------------------+ */
+/* | Choose the sample points. The indices of the sample points are stored | */
+/* | in the list S. | */
+/* +-----------------------------------------------------------------------+ */
+ actdeep = actmaxdeep;
+ direct_dirchoose_(anchor, s, &c__600, f, fmin, *eps, epsabs, levels, &maxpos, length,
+ &c__90000, &c__600, &c__3000, n, logfile, &cheat, &
+ kmax, &ifeasiblef, jones);
+/* +-----------------------------------------------------------------------+ */
+/* | Add other hyperrectangles to S, which have the same level and the same| */
+/* | function value at the center as the ones found above (that are stored | */
+/* | in S). This is only done if we use the original DIRECT algorithm. | */
+/* | JG 07/16/01 Added Errorflag. | */
+/* +-----------------------------------------------------------------------+ */
+ if (*algmethod == 0) {
+ direct_dirdoubleinsert_(anchor, s, &maxpos, point, f, &c__600, &c__90000,
+ &c__3000, ierror);
+ if (*ierror == -6) {
+ if (logfile)
+ fprintf(logfile,
+"WARNING: Capacity of array S in DIRDoubleInsert reached. Increase maxdiv.\n"
+"This means that there are a lot of hyperrectangles with the same function\n"
+"value at the center. We suggest to use our modification instead (Jones = 1)\n"
+ );
+ goto cleanup;
+ }
+ }
+ oldpos = minpos;
+/* +-----------------------------------------------------------------------+ */
+/* | Initialise the number of sample points in this outer loop. | */
+/* +-----------------------------------------------------------------------+ */
+ newtosample = 0;
+ i__2 = maxpos;
+ for (j = 1; j <= i__2; ++j) {
+ actdeep = s[j + 2999];
+/* +-----------------------------------------------------------------------+ */
+/* | If the actual index is a point to sample, do it. | */
+/* +-----------------------------------------------------------------------+ */
+ if (s[j - 1] > 0) {
+/* +-----------------------------------------------------------------------+ */
+/* | JG 09/24/00 Calculate the value delta used for sampling points. | */
+/* +-----------------------------------------------------------------------+ */
+ actdeep_div__ = direct_dirgetmaxdeep_(&s[j - 1], length, &c__90000,
+ n);
+ delta = thirds[actdeep_div__ + 1];
+ actdeep = s[j + 2999];
+/* +-----------------------------------------------------------------------+ */
+/* | If the current dept of division is only one under the maximal allowed | */
+/* | dept, stop the computation. | */
+/* +-----------------------------------------------------------------------+ */
+ if (actdeep + 1 >= mdeep) {
+ if (logfile)
+ fprintf(logfile, "WARNING: Maximum number of levels reached. Increase maxdeep.\n");
+ *ierror = -6;
+ goto L100;
+ }
+ actmaxdeep = MAX(actdeep,actmaxdeep);
+ help = s[j - 1];
+ if (! (anchor[actdeep + 1] == help)) {
+ pos1 = anchor[actdeep + 1];
+ while(! (point[pos1 - 1] == help)) {
+ pos1 = point[pos1 - 1];
+ }
+ point[pos1 - 1] = point[help - 1];
+ } else {
+ anchor[actdeep + 1] = point[help - 1];
+ }
+ if (actdeep < 0) {
+ actdeep = (integer) f[help - 1];
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Get the Directions in which to decrease the intervall-length. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirget_i__(length, &help, arrayi, &maxi, n, &c__90000);
+/* +-----------------------------------------------------------------------+ */
+/* | Sample the function. To do this, we first calculate the points where | */
+/* | we need to sample the function. After checking for errors, we then do | */
+/* | the actual evaluation of the function, again followed by checking for | */
+/* | errors. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirsamplepoints_(c__, arrayi, &delta, &help, &start, length,
+ logfile, f, &ifree, &maxi, point, &x[
+ 1], &l[1], fmin, &minpos, &u[1], n, &c__90000, &
+ c__600, &oops);
+ if (oops > 0) {
+ if (logfile)
+ fprintf(logfile, "WARNING: Error occured in routine DIRsamplepoints.\n");
+ *ierror = -4;
+ goto cleanup;
+ }
+ newtosample += maxi;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirsamplef_(c__, arrayi, &delta, &help, &start, length,
+ logfile, f, &ifree, &maxi, point, fcn, &x[
+ 1], &l[1], fmin, &minpos, &u[1], n, &c__90000, &
+ c__600, &oops, &fmax, &ifeasiblef, &iinfesiblef,
+ fcn_data);
+ if (oops > 0) {
+ if (logfile)
+ fprintf(logfile, "WARNING: Error occured in routine DIRsamplef.\n");
+ *ierror = -5;
+ goto cleanup;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Divide the intervalls. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirdivide_(&start, &actdeep_div__, length, point, arrayi, &
+ help, list2, w, &maxi, f, &c__90000, &c__600, n);
+/* +-----------------------------------------------------------------------+ */
+/* | Insert the new intervalls into the list (sorted). | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirinsertlist_(&start, anchor, point, f, &maxi, length, &
+ c__90000, &c__600, n, &help, jones);
+/* +-----------------------------------------------------------------------+ */
+/* | Increase the number of function evaluations. | */
+/* +-----------------------------------------------------------------------+ */
+ numfunc = numfunc + maxi + maxi;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | End of main loop. | */
+/* +-----------------------------------------------------------------------+ */
+/* L20: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | If there is a new minimum, show the actual iteration, the number of | */
+/* | function evaluations, the minimum value of f (so far) and the position| */
+/* | in the array. | */
+/* +-----------------------------------------------------------------------+ */
+ if (oldpos < minpos) {
+ if (logfile)
+ fprintf(logfile, "%d, %d, %g, %g\n",
+ t, numfunc, *fmin, fmax);
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | If no feasible point has been found, give out the iteration, the | */
+/* | number of function evaluations and a warning. | */
+/* +-----------------------------------------------------------------------+ */
+ if (ifeasiblef > 0) {
+ if (logfile)
+ fprintf(logfile, "No feasible point found in %d iterations "
+ "and %d function evaluations\n", t, numfunc);
+ }
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Termination Checks | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Calculate the index for the hyperrectangle at which | */
+/* | fmin is assumed. We then calculate the volume of this | */
+/* | hyperrectangle and store it in delta. This delta can be | */
+/* | used to stop DIRECT once the volume is below a certain | */
+/* | percentage of the original volume. Since the original | */
+/* | is 1 (scaled), we can stop once delta is below a certain | */
+/* | percentage, given by volper. | */
+/* +-----------------------------------------------------------------------+ */
+ *ierror = jones;
+ jones = 0;
+ actdeep_div__ = direct_dirgetlevel_(&minpos, length, &c__90000, n, jones);
+ jones = *ierror;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Use precalculated values to calculate volume. | */
+/* +-----------------------------------------------------------------------+ */
+ delta = thirds[actdeep_div__] * 100;
+ if (delta <= *volper) {
+ *ierror = 4;
+ if (logfile)
+ fprintf(logfile, "DIRECT stopped: Volume of S_min is "
+ "%g%% < %g%% of the original volume.\n",
+ delta, *volper);
+ goto L100;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/23/01 Calculate the measure for the hyperrectangle at which | */
+/* | fmin is assumed. If this measure is smaller then sigmaper,| */
+/* | we stop DIRECT. | */
+/* +-----------------------------------------------------------------------+ */
+ actdeep_div__ = direct_dirgetlevel_(&minpos, length, &c__90000, n, jones);
+ delta = levels[actdeep_div__];
+ if (delta <= *sigmaper) {
+ *ierror = 5;
+ if (logfile)
+ fprintf(logfile, "DIRECT stopped: Measure of S_min "
+ "= %g < %g.\n", delta, *sigmaper);
+ goto L100;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | If the best found function value is within fglper of the (known) | */
+/* | global minimum value, terminate. This only makes sense if this optimal| */
+/* | value is known, that is, in test problems. | */
+/* +-----------------------------------------------------------------------+ */
+ if ((*fmin - *fglobal) * 100 / divfactor <= *fglper) {
+ *ierror = 3;
+ if (logfile)
+ fprintf(logfile, "DIRECT stopped: fmin within fglper of global minimum.\n");
+ goto L100;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Find out if there are infeasible points which are near feasible ones. | */
+/* | If this is the case, replace the function value at the center of the | */
+/* | hyper rectangle by the lowest function value of a nearby function. | */
+/* | If no infeasible points exist (IInfesiblef = 0), skip this. | */
+/* +-----------------------------------------------------------------------+ */
+ if (iinfesiblef > 0) {
+ direct_dirreplaceinf_(&ifree, &ifreeold, f, c__, thirds, length, anchor,
+ point, &u[1], &l[1], &c__90000, &c__600, &c__64, n,
+ logfile, &fmax, jones);
+ }
+ ifreeold = ifree;
+/* +-----------------------------------------------------------------------+ */
+/* | If iepschange = 1, we use the epsilon change formula from Jones. | */
+/* +-----------------------------------------------------------------------+ */
+ if (iepschange == 1) {
+/* Computing MAX */
+ d__1 = fabs(*fmin) * 1e-4;
+ *eps = MAX(d__1,epsfix);
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | If no feasible point has been found yet, set the maximum number of | */
+/* | function evaluations to the number of evaluations already done plus | */
+/* | the budget given by the user. | */
+/* | If the budget has already be increased, increase it again. If a | */
+/* | feasible point has been found, remark that and reset flag. No further | */
+/* | increase is needed. | */
+/* +-----------------------------------------------------------------------+ */
+ if (increase == 1) {
+ *maxf = numfunc + oldmaxf;
+ if (ifeasiblef == 0) {
+ if (logfile)
+ fprintf(logfile, "DIRECT found a feasible point. The "
+ "adjusted budget is now set to %d.\n", *maxf);
+ increase = 0;
+ }
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Check if the number of function evaluations done is larger than the | */
+/* | allocated budget. If this is the case, check if a feasible point was | */
+/* | found. If this is a case, terminate. If no feasible point was found, | */
+/* | increase the budget and set flag increase. | */
+/* +-----------------------------------------------------------------------+ */
+ if (numfunc > *maxf) {
+ if (ifeasiblef == 0) {
+ *ierror = 1;
+ if (logfile)
+ fprintf(logfile, "DIRECT stopped: numfunc >= maxf.\n");
+ goto L100;
+ } else {
+ increase = 1;
+ if (logfile)
+ fprintf(logfile,
+"DIRECT could not find a feasible point after %d function evaluations.\n"
+"DIRECT continues until a feasible point is found.\n", numfunc);
+ *maxf = numfunc + oldmaxf;
+ }
+ }
+/* L10: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | End of main loop. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | The algorithm stopped after maxT iterations. | */
+/* +-----------------------------------------------------------------------+ */
+ *ierror = 2;
+ if (logfile)
+ fprintf(logfile, "DIRECT stopped: maxT iterations.\n");
+
+L100:
+/* +-----------------------------------------------------------------------+ */
+/* | Store the position of the minimum in x. | */
+/* +-----------------------------------------------------------------------+ */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ x[i__] = c__[minpos + i__ * 90000 - 90001] * l[i__] + l[i__] * u[i__];
+ u[i__] = oldu[i__ - 1];
+ l[i__] = oldl[i__ - 1];
+/* L50: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Store the number of function evaluations in maxf. | */
+/* +-----------------------------------------------------------------------+ */
+ *maxf = numfunc;
+/* +-----------------------------------------------------------------------+ */
+/* | Give out a summary of the run. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirsummary_(logfile, &x[1], &l[1], &u[1], n, fmin, fglobal, &numfunc,
+ ierror);
+/* +-----------------------------------------------------------------------+ */
+/* | Format statements. | */
+/* +-----------------------------------------------------------------------+ */
+
+ cleanup:
+#define MY_FREE(p) if (p) free(p)
+ MY_FREE(c__);
+ MY_FREE(f);
+ MY_FREE(s);
+ MY_FREE(w);
+ MY_FREE(oldl);
+ MY_FREE(oldu);
+ MY_FREE(list2);
+ MY_FREE(point);
+ MY_FREE(anchor);
+ MY_FREE(length);
+ MY_FREE(arrayi);
+ MY_FREE(levels);
+ MY_FREE(thirds);
+} /* direct_ */
+
--- /dev/null
+/* DIRparallel.f -- translated by f2c (version 20050501).
+
+ f2c output hand-cleaned by SGJ (August 2007).
+*/
+
+#include "direct-internal.h"
+
+/* Table of constant values */
+
+static integer c__0 = 0;
+static integer c_n1 = -1;
+
+/* +-----------------------------------------------------------------------+ */
+/* | Program : Direct.f (subfile DIRseriell.f) | */
+/* | Last modified : 02-22-01 | */
+/* | Written by : Joerg Gablonsky | */
+/* | Subroutines, which differ depending on the serial or parallel version.| */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Parallel Direct. This routine replaces the normal main routine DIRect.| */
+/* | In it, we find out if this pe is the master or slave. If it is the | */
+/* | master, it calls the serial DIRect main routine. The only routine that| */
+/* | has to change for parallel Direct is DIRSamplef, where the actual | */
+/* | sampling of the function is done. If we are on the slave, wait for | */
+/* | either the coordinates of a point to sample the function or the | */
+/* | termination signal. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ int direct_pardirect_(fp fcn, doublereal *x, integer *n,
+ doublereal *eps, integer *maxf, integer *maxt, doublereal *fmin,
+ doublereal *l, doublereal *u, integer *algmethod, integer *ierror,
+ FILE *logfile, doublereal *fglobal, doublereal *fglper, doublereal
+ *volper, doublereal *sigmaper, void *fcn_data)
+{
+ /* System generated locals */
+ integer i__1;
+
+ /* Local variables */
+ integer i__, k;
+ integer tid;
+ integer flag__;
+ doublereal fval;
+ integer tids[360], kret;
+ integer mytid;
+ doublereal fscale;
+ integer nprocs;
+
+/* +-----------------------------------------------------------------------+ */
+/* | Parameters | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | The maximum of function evaluations allowed. | */
+/* | The maximum dept of the algorithm. | */
+/* | The maximum number of divisions allowed. | */
+/* | The maximal dimension of the problem. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Global Variables. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | External Variables. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | User Variables. | */
+/* | These can be used to pass user defined data to the function to be | */
+/* | optimized. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Parallel programming variables | */
+/* +-----------------------------------------------------------------------+ */
+/* maxprocs should be >= the number of processes used for DIRECT */
+/* +-----------------------------------------------------------------------+ */
+/* | End of parallel programming variables | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Internal variables | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 02/28/01 Begin of parallel additions | */
+/* | DETERMINE MASTER PROCESSOR. GET TIDS OF ALL PROCESSORS. | */
+/* +-----------------------------------------------------------------------+ */
+ /* Parameter adjustments */
+ --u;
+ --l;
+ --x;
+
+ /* Function Body */
+ getmytidif_(&mytid);
+ getnprocsif_(&nprocs);
+ gettidif_(&c__0, tids);
+/* +-----------------------------------------------------------------------+ */
+/* | If I am the master get the other tids and start running DIRECT. | */
+/* | Otherwise, branch off to do function evaluations. | */
+/* +-----------------------------------------------------------------------+ */
+ if (mytid == tids[0]) {
+ i__1 = nprocs - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ gettidif_(&i__, &tids[i__]);
+/* L46: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Call Direct main routine. This routine calls DIRSamplef for the | */
+/* | function evaluations, which are then done in parallel. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_direct_(fcn, &x[1], n, eps, maxf, maxt, fmin, &l[1], &u[1],
+ algmethod, ierror, logfile, fglobal, fglper, volper, sigmaper,
+ fcn_data);
+/* +-----------------------------------------------------------------------+ */
+/* | Send exit message to rest of pe's. | */
+/* +-----------------------------------------------------------------------+ */
+ flag__ = 0;
+ i__1 = nprocs;
+ for (tid = 2; tid <= i__1; ++tid) {
+ mastersendif_(&tids[tid - 1], &tids[tid - 1], n, &flag__, &flag__,
+ &x[1], &u[1], &l[1], &x[1]);
+/* L200: */
+ }
+ } else {
+/* +-----------------------------------------------------------------------+ */
+/* | This is what the slaves do!! | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Receive the first point from the master processor. | */
+/* +-----------------------------------------------------------------------+ */
+ slaverecvif_(tids, &c_n1, n, &flag__, &k, &fscale, &u[1], &l[1], &x[1]
+ );
+/* +-----------------------------------------------------------------------+ */
+/* | Repeat until master signals to stop. | */
+/* +-----------------------------------------------------------------------+ */
+ while(flag__ > 0) {
+/* +-----------------------------------------------------------------------+ */
+/* | Evaluate f(x). | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &fval, &kret, &
+ fcn_data);
+/* +-----------------------------------------------------------------------+ */
+/* | Send result and wait for next point / message with signal to stop. | */
+/* +-----------------------------------------------------------------------+ */
+ slavesendif_(tids, &mytid, &k, &mytid, &fval, &kret);
+ slaverecvif_(tids, &c_n1, n, &flag__, &k, &fscale, &u[1], &l[1], &
+ x[1]);
+ }
+ }
+ return 0;
+} /* pardirect_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | Subroutine for sampling. This sampling is done in parallel, the master| */
+/* | prozessor is also evaluating the function sometimes. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirsamplef_(doublereal *c__, integer *arrayi, doublereal
+ *delta, integer *sample, integer *new__, integer *length,
+ FILE *logfile, doublereal *f, integer *free, integer *maxi,
+ integer *point, fp fcn, doublereal *x, doublereal *l, doublereal *
+ fmin, integer *minpos, doublereal *u, integer *n, integer *maxfunc,
+ integer *maxdeep, integer *oops, doublereal *fmax, integer *
+ ifeasiblef, integer *iinfesiblef, void *fcn_data)
+{
+ /* System generated locals */
+ integer length_dim1, length_offset, c_dim1, c_offset, f_dim1, f_offset,
+ i__1;
+ doublereal d__1;
+
+ /* Local variables */
+ integer i__, j, k, helppoint, tid, pos;
+ integer flag__, tids[360], kret, npts;
+ doublereal fhelp;
+ integer oldpos, nprocs, datarec;
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 fcn must be declared external. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Removed fcn. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
+/* | Added variable to keep track if feasible point was found. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Variables to pass user defined data to the function to be optimized. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Parallel programming variables. | */
+/* +-----------------------------------------------------------------------+ */
+/* JG 09/05/00 Increase this if more processors are used. */
+/* +-----------------------------------------------------------------------+ */
+/* | Find out the id's of all processors. | */
+/* +-----------------------------------------------------------------------+ */
+ /* Parameter adjustments */
+ --u;
+ --l;
+ --x;
+ --arrayi;
+ --point;
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+ c_dim1 = *maxfunc;
+ c_offset = 1 + c_dim1;
+ c__ -= c_offset;
+
+ /* Function Body */
+ getnprocsif_(&nprocs);
+ i__1 = nprocs - 1;
+ for (i__ = 0; i__ <= i__1; ++i__) {
+ gettidif_(&i__, &tids[i__]);
+/* L46: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Set the pointer to the first function to be evaluated, | */
+/* | store this position also in helppoint. | */
+/* +-----------------------------------------------------------------------+ */
+ pos = *new__;
+ helppoint = pos;
+/* +-----------------------------------------------------------------------+ */
+/* | Iterate over all points, where the function should be | */
+/* | evaluated. | */
+/* +-----------------------------------------------------------------------+ */
+ flag__ = 1;
+ npts = *maxi + *maxi;
+ k = 1;
+ while(k <= npts && k < nprocs) {
+/* +-----------------------------------------------------------------------+ */
+/* | tid is the id of the prozessor the next points is send to. | */
+/* +-----------------------------------------------------------------------+ */
+ tid = k + 1;
+/* +-----------------------------------------------------------------------+ */
+/* | Copy the position into the helparray x. | */
+/* +-----------------------------------------------------------------------+ */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ x[i__] = c__[pos + i__ * c_dim1];
+/* L60: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Send the point. | */
+/* +-----------------------------------------------------------------------+ */
+ mastersendif_(&tids[tid - 1], &tids[tid - 1], n, &flag__, &pos, &x[1],
+ &u[1], &l[1], &x[1]);
+ ++k;
+ pos = point[pos];
+/* +-----------------------------------------------------------------------+ */
+/* | Get the next point. | */
+/* +-----------------------------------------------------------------------+ */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Get data until it is all received. | */
+/* +-----------------------------------------------------------------------+ */
+ datarec = 0;
+ while(datarec < npts) {
+ if ((doublereal) datarec / (doublereal) nprocs - datarec / nprocs <
+ 1e-5 && k <= npts) {
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ x[i__] = c__[pos + i__ * c_dim1];
+/* L165: */
+ }
+ direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &fhelp, &kret,
+ fcn_data);
+ oldpos = pos;
+ f[oldpos + f_dim1] = fhelp;
+ ++datarec;
+/* +-----------------------------------------------------------------------+ */
+/* | Remember if an infeasible point has been found. | */
+/* +-----------------------------------------------------------------------+ */
+ *iinfesiblef = MAX(*iinfesiblef,kret);
+ if (kret == 0) {
+/* +-----------------------------------------------------------------------+ */
+/* | if the function evaluation was O.K., set the flag in | */
+/* | f(pos,2). | */
+/* +-----------------------------------------------------------------------+ */
+ f[oldpos + (f_dim1 << 1)] = 0.;
+ *ifeasiblef = 0;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
+/* +-----------------------------------------------------------------------+ */
+/* Computing MAX */
+ d__1 = f[pos + f_dim1];
+ *fmax = MAX(d__1,*fmax);
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Remember if an infeasible point has been found. | */
+/* +-----------------------------------------------------------------------+ */
+ *iinfesiblef = MAX(*iinfesiblef,kret);
+ if (kret == 1) {
+/* +-----------------------------------------------------------------------+ */
+/* | If the function could not be evaluated at the given point, | */
+/* | set flag to mark this (f(pos,2) and store the maximum | */
+/* | box-sidelength in f(pos,1). | */
+/* +-----------------------------------------------------------------------+ */
+ f[oldpos + (f_dim1 << 1)] = 2.;
+ f[oldpos + f_dim1] = *fmax;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | If the function could not be evaluated due to a failure in | */
+/* | the setup, mark this. | */
+/* +-----------------------------------------------------------------------+ */
+ if (kret == -1) {
+ f[oldpos + (f_dim1 << 1)] = -1.;
+ }
+ ++k;
+ pos = point[pos];
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Recover where to store the value. | */
+/* +-----------------------------------------------------------------------+ */
+ masterrecvif_(&c_n1, &c_n1, &oldpos, &tid, &fhelp, &kret);
+ f[oldpos + f_dim1] = fhelp;
+ ++datarec;
+/* +-----------------------------------------------------------------------+ */
+/* | Remember if an infeasible point has been found. | */
+/* +-----------------------------------------------------------------------+ */
+ *iinfesiblef = MAX(*iinfesiblef,kret);
+ if (kret == 0) {
+/* +-----------------------------------------------------------------------+ */
+/* | if the function evaluation was O.K., set the flag in | */
+/* | f(pos,2). | */
+/* +-----------------------------------------------------------------------+ */
+ f[oldpos + (f_dim1 << 1)] = 0.;
+ *ifeasiblef = 0;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
+/* +-----------------------------------------------------------------------+ */
+/* Computing MAX */
+ d__1 = f[oldpos + f_dim1];
+ *fmax = MAX(d__1,*fmax);
+ }
+ if (kret == 1) {
+/* +-----------------------------------------------------------------------+ */
+/* | If the function could not be evaluated at the given point, | */
+/* | set flag to mark this (f(pos,2) and store the maximum | */
+/* | box-sidelength in f(pos,1). | */
+/* +-----------------------------------------------------------------------+ */
+ f[oldpos + (f_dim1 << 1)] = 2.;
+ f[oldpos + f_dim1] = *fmax;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | If the function could not be evaluated due to a failure in | */
+/* | the setup, mark this. | */
+/* +-----------------------------------------------------------------------+ */
+ if (kret == -1) {
+ f[oldpos + (f_dim1 << 1)] = -1.;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Send data until it is all sent. | */
+/* +-----------------------------------------------------------------------+ */
+ if (k <= npts) {
+/* +-----------------------------------------------------------------------+ */
+/* | Copy the position into the helparray x. | */
+/* +-----------------------------------------------------------------------+ */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ x[i__] = c__[pos + i__ * c_dim1];
+/* L160: */
+ }
+ mastersendif_(&tid, &tid, n, &flag__, &pos, &x[1], &u[1], &l[1], &
+ x[1]);
+ ++k;
+ pos = point[pos];
+ }
+ }
+ pos = helppoint;
+/* +-----------------------------------------------------------------------+ */
+/* | Iterate over all evaluated points and see, if the minimal | */
+/* | value of the function has changed. If this has happend, | */
+/* | store the minimal value and its position in the array. | */
+/* | Attention: Only valied values are checked!! | */
+/* +-----------------------------------------------------------------------+ */
+ i__1 = *maxi + *maxi;
+ for (j = 1; j <= i__1; ++j) {
+ if (f[pos + f_dim1] < *fmin && f[pos + (f_dim1 << 1)] == 0.) {
+ *fmin = f[pos + f_dim1];
+ *minpos = pos;
+ }
+ pos = point[pos];
+/* L50: */
+ }
+} /* dirsamplef_ */
--- /dev/null
+/* DIRserial.f -- translated by f2c (version 20050501).
+
+ f2c output hand-cleaned by SGJ (August 2007).
+*/
+
+#include "direct-internal.h"
+
+/* +-----------------------------------------------------------------------+ */
+/* | Program : Direct.f (subfile DIRserial.f) | */
+/* | Last modified : 04-12-2001 | */
+/* | Written by : Joerg Gablonsky | */
+/* | SUBROUTINEs, which differ depENDing on the serial or parallel version.| */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE for sampling. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirsamplef_(doublereal *c__, integer *arrayi, doublereal
+ *delta, integer *sample, integer *new__, integer *length,
+ FILE *logfile, doublereal *f, integer *free, integer *maxi,
+ integer *point, fp fcn, doublereal *x, doublereal *l, doublereal *
+ fmin, integer *minpos, doublereal *u, integer *n, integer *maxfunc,
+ integer *maxdeep, integer *oops, doublereal *fmax, integer *
+ ifeasiblef, integer *iinfesiblef, void *fcn_data)
+{
+ /* System generated locals */
+ integer length_dim1, length_offset, c_dim1, c_offset, f_dim1, f_offset,
+ i__1, i__2;
+ doublereal d__1;
+
+ /* Local variables */
+ integer i__, j, helppoint, pos, kret;
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 fcn must be declared external. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Removed fcn. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
+/* | Added variable to keep track IF feasible point was found. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Variables to pass user defined data to the function to be optimized. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Set the pointer to the first function to be evaluated, | */
+/* | store this position also in helppoint. | */
+/* +-----------------------------------------------------------------------+ */
+ /* Parameter adjustments */
+ --u;
+ --l;
+ --x;
+ --arrayi;
+ --point;
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+ c_dim1 = *maxfunc;
+ c_offset = 1 + c_dim1;
+ c__ -= c_offset;
+
+ /* Function Body */
+ pos = *new__;
+ helppoint = pos;
+/* +-----------------------------------------------------------------------+ */
+/* | Iterate over all points, where the function should be | */
+/* | evaluated. | */
+/* +-----------------------------------------------------------------------+ */
+ i__1 = *maxi + *maxi;
+ for (j = 1; j <= i__1; ++j) {
+/* +-----------------------------------------------------------------------+ */
+/* | Copy the position into the helparrayy x. | */
+/* +-----------------------------------------------------------------------+ */
+ i__2 = *n;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ x[i__] = c__[pos + i__ * c_dim1];
+/* L60: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Call the function. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &f[pos + f_dim1], &kret,
+ fcn_data);
+/* +-----------------------------------------------------------------------+ */
+/* | Remember IF an infeasible point has been found. | */
+/* +-----------------------------------------------------------------------+ */
+ *iinfesiblef = MAX(*iinfesiblef,kret);
+ if (kret == 0) {
+/* +-----------------------------------------------------------------------+ */
+/* | IF the function evaluation was O.K., set the flag in | */
+/* | f(pos,2). Also mark that a feasible point has been found. | */
+/* +-----------------------------------------------------------------------+ */
+ f[pos + (f_dim1 << 1)] = 0.;
+ *ifeasiblef = 0;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
+/* +-----------------------------------------------------------------------+ */
+/* Computing MAX */
+ d__1 = f[pos + f_dim1];
+ *fmax = MAX(d__1,*fmax);
+ }
+ if (kret >= 1) {
+/* +-----------------------------------------------------------------------+ */
+/* | IF the function could not be evaluated at the given point, | */
+/* | set flag to mark this (f(pos,2) and store the maximum | */
+/* | box-sidelength in f(pos,1). | */
+/* +-----------------------------------------------------------------------+ */
+ f[pos + (f_dim1 << 1)] = 2.;
+ f[pos + f_dim1] = *fmax;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | IF the function could not be evaluated due to a failure in | */
+/* | the setup, mark this. | */
+/* +-----------------------------------------------------------------------+ */
+ if (kret == -1) {
+ f[pos + (f_dim1 << 1)] = -1.;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Set the position to the next point, at which the function | */
+/* | should e evaluated. | */
+/* +-----------------------------------------------------------------------+ */
+ pos = point[pos];
+/* L40: */
+ }
+ pos = helppoint;
+/* +-----------------------------------------------------------------------+ */
+/* | Iterate over all evaluated points and see, IF the minimal | */
+/* | value of the function has changed. IF this has happEND, | */
+/* | store the minimal value and its position in the array. | */
+/* | Attention: Only valied values are checked!! | */
+/* +-----------------------------------------------------------------------+ */
+ i__1 = *maxi + *maxi;
+ for (j = 1; j <= i__1; ++j) {
+ if (f[pos + f_dim1] < *fmin && f[pos + (f_dim1 << 1)] == 0.) {
+ *fmin = f[pos + f_dim1];
+ *minpos = pos;
+ }
+ pos = point[pos];
+/* L50: */
+ }
+} /* dirsamplef_ */
--- /dev/null
+/* DIRsubrout.f -- translated by f2c (version 20050501).
+
+ f2c output hand-cleaned by SGJ (August 2007).
+*/
+
+#include <math.h>
+#include "direct-internal.h"
+
+/* Table of constant values */
+
+static integer c__1 = 1;
+static integer c__32 = 32;
+static integer c__0 = 0;
+
+/* +-----------------------------------------------------------------------+ */
+/* | INTEGER Function DIRGetlevel | */
+/* | Returns the level of the hyperrectangle. Depending on the value of the| */
+/* | global variable JONES. IF JONES equals 0, the level is given by | */
+/* | kN + p, where the rectangle has p sides with a length of| */
+/* | 1/3^(k+1), and N-p sides with a length of 1/3^k. | */
+/* | If JONES equals 1, the level is the power of 1/3 of the length of the | */
+/* | longest side hyperrectangle. | */
+/* | | */
+/* | On Return : | */
+/* | the maximal length | */
+/* | | */
+/* | pos -- the position of the midpoint in the array length | */
+/* | length -- the array with the dimensions | */
+/* | maxfunc -- the leading dimension of length | */
+/* | n -- the dimension of the problem | */
+/* | | */
+/* +-----------------------------------------------------------------------+ */
+integer direct_dirgetlevel_(integer *pos, integer *length, integer *maxfunc, integer
+ *n, integer jones)
+{
+ /* System generated locals */
+ integer length_dim1, length_offset, ret_val, i__1;
+
+ /* Local variables */
+ integer i__, k, p, help;
+
+/* JG 09/15/00 Added variable JONES (see above) */
+ /* Parameter adjustments */
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+
+ /* Function Body */
+ if (jones == 0) {
+ help = length[*pos + length_dim1];
+ k = help;
+ p = 1;
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+ if (length[*pos + i__ * length_dim1] < k) {
+ k = length[*pos + i__ * length_dim1];
+ }
+ if (length[*pos + i__ * length_dim1] == help) {
+ ++p;
+ }
+/* L100: */
+ }
+ if (k == help) {
+ ret_val = k * *n + *n - p;
+ } else {
+ ret_val = k * *n + p;
+ }
+ } else {
+ help = length[*pos + length_dim1];
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+ if (length[*pos + i__ * length_dim1] < help) {
+ help = length[*pos + i__ * length_dim1];
+ }
+/* L10: */
+ }
+ ret_val = help;
+ }
+ return ret_val;
+} /* dirgetlevel_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | Program : Direct.f (subfile DIRsubrout.f) | */
+/* | Last modified : 07-16-2001 | */
+/* | Written by : Joerg Gablonsky | */
+/* | Subroutines used by the algorithm DIRECT. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRChoose | */
+/* | Decide, which is the next sampling point. | */
+/* | Changed 09/25/00 JG | */
+/* | Added maxdiv to call and changed S to size maxdiv. | */
+/* | Changed 01/22/01 JG | */
+/* | Added Ifeasiblef to call to keep track if a feasible point has| */
+/* | been found. | */
+/* | Changed 07/16/01 JG | */
+/* | Changed if statement to prevent run-time errors. |
+ | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirchoose_(integer *anchor, integer *s, integer *actdeep,
+ doublereal *f, doublereal *fmin, doublereal epsrel, doublereal epsabs, doublereal *thirds,
+ integer *maxpos, integer *length, integer *maxfunc, integer *maxdeep,
+ integer *maxdiv, integer *n, FILE *logfile,
+ integer *cheat, doublereal *kmax, integer *ifeasiblef, integer jones)
+{
+ /* System generated locals */
+ integer s_dim1, s_offset, length_dim1, length_offset, f_dim1, f_offset,
+ i__1;
+
+ /* Local variables */
+ integer i__, j, k;
+ doublereal helplower;
+ integer i___, j___;
+ doublereal helpgreater;
+ integer novaluedeep;
+ doublereal help2;
+ integer novalue;
+
+ /* Parameter adjustments */
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ ++anchor;
+ s_dim1 = *maxdiv;
+ s_offset = 1 + s_dim1;
+ s -= s_offset;
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+
+ /* Function Body */
+ helplower = 1e20;
+ helpgreater = 0.;
+ k = 1;
+ if (*ifeasiblef >= 1) {
+ i__1 = *actdeep;
+ for (j = 0; j <= i__1; ++j) {
+ if (anchor[j] > 0) {
+ s[k + s_dim1] = anchor[j];
+ s[k + (s_dim1 << 1)] = direct_dirgetlevel_(&s[k + s_dim1], &length[
+ length_offset], maxfunc, n, jones);
+ goto L12;
+ }
+/* L1001: */
+ }
+L12:
+ ++k;
+ *maxpos = 1;
+ return;
+ } else {
+ i__1 = *actdeep;
+ for (j = 0; j <= i__1; ++j) {
+ if (anchor[j] > 0) {
+ s[k + s_dim1] = anchor[j];
+ s[k + (s_dim1 << 1)] = direct_dirgetlevel_(&s[k + s_dim1], &length[
+ length_offset], maxfunc, n, jones);
+ ++k;
+ }
+/* L10: */
+ }
+ }
+ novalue = 0;
+ if (anchor[-1] > 0) {
+ novalue = anchor[-1];
+ novaluedeep = direct_dirgetlevel_(&novalue, &length[length_offset], maxfunc,
+ n, jones);
+ }
+ *maxpos = k - 1;
+ i__1 = *maxdeep;
+ for (j = k - 1; j <= i__1; ++j) {
+ s[k + s_dim1] = 0;
+/* L11: */
+ }
+ for (j = *maxpos; j >= 1; --j) {
+ helplower = 1e20;
+ helpgreater = 0.;
+ j___ = s[j + s_dim1];
+ i__1 = j - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ i___ = s[i__ + s_dim1];
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Changed IF statement into two to prevent run-time errors | */
+/* | which could occur if the compiler checks the second | */
+/* | expression in an .AND. statement although the first | */
+/* | statement is already not true. | */
+/* +-----------------------------------------------------------------------+ */
+ if (i___ > 0 && ! (i__ == j)) {
+ if (f[i___ + (f_dim1 << 1)] <= 1.) {
+ help2 = thirds[s[i__ + (s_dim1 << 1)]] - thirds[s[j + (
+ s_dim1 << 1)]];
+ help2 = (f[i___ + f_dim1] - f[j___ + f_dim1]) / help2;
+ if (help2 <= 0.) {
+ if (logfile)
+ fprintf(logfile, "thirds > 0, help2 <= 0\n");
+ goto L60;
+ }
+ if (help2 < helplower) {
+ if (logfile)
+ fprintf(logfile, "helplower = %g\n", help2);
+ helplower = help2;
+ }
+ }
+ }
+/* L30: */
+ }
+ i__1 = *maxpos;
+ for (i__ = j + 1; i__ <= i__1; ++i__) {
+ i___ = s[i__ + s_dim1];
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Changed IF statement into two to prevent run-time errors | */
+/* | which could occur if the compiler checks the second | */
+/* | expression in an .AND. statement although the first | */
+/* | statement is already not true. | */
+/* +-----------------------------------------------------------------------+ */
+ if (i___ > 0 && ! (i__ == j)) {
+ if (f[i___ + (f_dim1 << 1)] <= 1.) {
+ help2 = thirds[s[i__ + (s_dim1 << 1)]] - thirds[s[j + (
+ s_dim1 << 1)]];
+ help2 = (f[i___ + f_dim1] - f[j___ + f_dim1]) / help2;
+ if (help2 <= 0.) {
+ if (logfile)
+ fprintf(logfile, "thirds < 0, help2 <= 0\n");
+ goto L60;
+ }
+ if (help2 > helpgreater) {
+ if (logfile)
+ fprintf(logfile, "helpgreater = %g\n", help2);
+ helpgreater = help2;
+ }
+ }
+ }
+/* L31: */
+ }
+ if (helplower > 1e20 && helpgreater > 0.) {
+ helplower = helpgreater;
+ helpgreater += -1.;
+ }
+ if (helpgreater <= helplower) {
+ if (*cheat == 1 && helplower > *kmax) {
+ helplower = *kmax;
+ }
+ if (f[j___ + f_dim1] - helplower * thirds[s[j + (s_dim1 << 1)]] >
+ MIN(*fmin - epsrel * fabs(*fmin),
+ *fmin - epsabs)) {
+ if (logfile)
+ fprintf(logfile, "> fmin - epslfminl\n");
+ goto L60;
+ }
+ } else {
+ if (logfile)
+ fprintf(logfile, "helpgreater > helplower: %g %g %g\n",
+ helpgreater, helplower, helpgreater - helplower);
+ goto L60;
+ }
+ goto L40;
+L60:
+ s[j + s_dim1] = 0;
+L40:
+ ;
+ }
+ if (novalue > 0) {
+ ++(*maxpos);
+ s[*maxpos + s_dim1] = novalue;
+ s[*maxpos + (s_dim1 << 1)] = novaluedeep;
+ }
+} /* dirchoose_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRDoubleInsert | */
+/* | Routine to make sure that if there are several potential optimal | */
+/* | hyperrectangles of the same level (i.e. hyperrectangles that have| */
+/* | the same level and the same function value at the center), all of| */
+/* | them are divided. This is the way as originally described in | */
+/* | Jones et.al. | */
+/* | JG 07/16/01 Added errorflag to calling sequence. We check if more | */
+/* | we reach the capacity of the array S. If this happens, we | */
+/* | return to the main program with an error. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirdoubleinsert_(integer *anchor, integer *s, integer *
+ maxpos, integer *point, doublereal *f, integer *maxdeep, integer *
+ maxfunc, integer *maxdiv, integer *ierror)
+{
+ /* System generated locals */
+ integer s_dim1, s_offset, f_dim1, f_offset, i__1;
+
+ /* Local variables */
+ integer i__, oldmaxpos, pos, help, iflag, actdeep;
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Added flag to prevent run time-errors on some systems. | */
+/* +-----------------------------------------------------------------------+ */
+ /* Parameter adjustments */
+ ++anchor;
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ --point;
+ s_dim1 = *maxdiv;
+ s_offset = 1 + s_dim1;
+ s -= s_offset;
+
+ /* Function Body */
+ oldmaxpos = *maxpos;
+ i__1 = oldmaxpos;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (s[i__ + s_dim1] > 0) {
+ actdeep = s[i__ + (s_dim1 << 1)];
+ help = anchor[actdeep];
+ pos = point[help];
+ iflag = 0;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Added flag to prevent run time-errors on some systems. On | */
+/* | some systems the second conditions in an AND statement is | */
+/* | evaluated even if the first one is already not true. | */
+/* +-----------------------------------------------------------------------+ */
+ while(pos > 0 && iflag == 0) {
+ if (f[pos + f_dim1] - f[help + f_dim1] <= 1e-13) {
+ if (*maxpos < *maxdiv) {
+ ++(*maxpos);
+ s[*maxpos + s_dim1] = pos;
+ s[*maxpos + (s_dim1 << 1)] = actdeep;
+ pos = point[pos];
+ } else {
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Maximum number of elements possible in S has been reached!| */
+/* +-----------------------------------------------------------------------+ */
+ *ierror = -6;
+ return;
+ }
+ } else {
+ iflag = 1;
+ }
+ }
+ }
+/* L10: */
+ }
+} /* dirdoubleinsert_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | INTEGER Function GetmaxDeep | */
+/* | function to get the maximal length (1/length) of the n-dimensional | */
+/* | rectangle with midpoint pos. | */
+/* | | */
+/* | On Return : | */
+/* | the maximal length | */
+/* | | */
+/* | pos -- the position of the midpoint in the array length | */
+/* | length -- the array with the dimensions | */
+/* | maxfunc -- the leading dimension of length | */
+/* | n -- the dimension of the problem | */
+/* | | */
+/* +-----------------------------------------------------------------------+ */
+integer direct_dirgetmaxdeep_(integer *pos, integer *length, integer *maxfunc,
+ integer *n)
+{
+ /* System generated locals */
+ integer length_dim1, length_offset, i__1, i__2, i__3;
+
+ /* Local variables */
+ integer i__, help;
+
+ /* Parameter adjustments */
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+
+ /* Function Body */
+ help = length[*pos + length_dim1];
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+/* Computing MIN */
+ i__2 = help, i__3 = length[*pos + i__ * length_dim1];
+ help = MIN(i__2,i__3);
+/* L10: */
+ }
+ return help;
+} /* dirgetmaxdeep_ */
+
+static integer isinbox_(doublereal *x, doublereal *a, doublereal *b, integer *n,
+ integer *lmaxdim)
+{
+ /* System generated locals */
+ integer ret_val, i__1;
+
+ /* Local variables */
+ integer outofbox, i__;
+
+ /* Parameter adjustments */
+ --b;
+ --a;
+ --x;
+
+ /* Function Body */
+ outofbox = 1;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (a[i__] > x[i__] || b[i__] < x[i__]) {
+ outofbox = 0;
+ goto L1010;
+ }
+/* L1000: */
+ }
+L1010:
+ ret_val = outofbox;
+ return ret_val;
+} /* isinbox_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG Added 09/25/00 | */
+/* | | */
+/* | SUBROUTINE DIRResortlist | */
+/* | | */
+/* | Resort the list so that the infeasible point is in the list with the | */
+/* | replaced value. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ static void dirresortlist_(integer *replace, integer *anchor,
+ doublereal *f, integer *point, integer *length, integer *n, integer *
+ maxfunc, integer *maxdim, integer *maxdeep, FILE *logfile,
+ integer jones)
+{
+ /* System generated locals */
+ integer f_dim1, f_offset, length_dim1, length_offset, i__1;
+
+ /* Local variables */
+ integer i__, l, pos;
+ integer start;
+
+/* +-----------------------------------------------------------------------+ */
+/* | Get the length of the hyper rectangle with infeasible mid point and | */
+/* | Index of the corresponding list. | */
+/* +-----------------------------------------------------------------------+ */
+/* JG 09/25/00 Replaced with DIRgetlevel */
+/* l = DIRgetmaxDeep(replace,length,maxfunc,n) */
+ /* Parameter adjustments */
+ --point;
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+ ++anchor;
+
+ /* Function Body */
+ l = direct_dirgetlevel_(replace, &length[length_offset], maxfunc, n, jones);
+ start = anchor[l];
+/* +-----------------------------------------------------------------------+ */
+/* | If the hyper rectangle with infeasibel midpoint is already the start | */
+/* | of the list, give out message, nothing to do. | */
+/* +-----------------------------------------------------------------------+ */
+ if (*replace == start) {
+/* write(logfile,*) 'No resorting of list necessarry, since new ', */
+/* + 'point is already anchor of list .',l */
+ } else {
+/* +-----------------------------------------------------------------------+ */
+/* | Take the hyper rectangle with infeasible midpoint out of the list. | */
+/* +-----------------------------------------------------------------------+ */
+ pos = start;
+ i__1 = *maxfunc;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (point[pos] == *replace) {
+ point[pos] = point[*replace];
+ goto L20;
+ } else {
+ pos = point[pos];
+ }
+ if (pos == 0) {
+ if (logfile)
+ fprintf(logfile, "Error in DIRREsortlist: "
+ "We went through the whole list\n"
+ "and could not find the point to replace!!\n");
+ goto L20;
+ }
+/* L10: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | If the anchor of the list has a higher value than the value of a | */
+/* | nearby point, put the infeasible point at the beginning of the list. | */
+/* +-----------------------------------------------------------------------+ */
+L20:
+ if (f[start + f_dim1] > f[*replace + f_dim1]) {
+ anchor[l] = *replace;
+ point[*replace] = start;
+/* write(logfile,*) 'Point is replacing current anchor for ' */
+/* + , 'this list ',l,replace,start */
+ } else {
+/* +-----------------------------------------------------------------------+ */
+/* | Insert the point into the list according to its (replaced) function | */
+/* | value. | */
+/* +-----------------------------------------------------------------------+ */
+ pos = start;
+ i__1 = *maxfunc;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+/* +-----------------------------------------------------------------------+ */
+/* | The point has to be added at the end of the list. | */
+/* +-----------------------------------------------------------------------+ */
+ if (point[pos] == 0) {
+ point[*replace] = point[pos];
+ point[pos] = *replace;
+/* write(logfile,*) 'Point is added at the end of the ' */
+/* + , 'list ',l, replace */
+ goto L40;
+ } else {
+ if (f[point[pos] + f_dim1] > f[*replace + f_dim1]) {
+ point[*replace] = point[pos];
+ point[pos] = *replace;
+/* write(logfile,*) 'There are points with a higher ' */
+/* + ,'f-value in the list ',l,replace, pos */
+ goto L40;
+ }
+ pos = point[pos];
+ }
+/* L30: */
+ }
+L40:
+ pos = pos;
+ }
+ }
+} /* dirresortlist_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG Added 09/25/00 | */
+/* | SUBROUTINE DIRreplaceInf | */
+/* | | */
+/* | Find out if there are infeasible points which are near feasible ones. | */
+/* | If this is the case, replace the function value at the center of the | */
+/* | hyper rectangle by the lowest function value of a nearby function. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirreplaceinf_(integer *free, integer *freeold,
+ doublereal *f, doublereal *c__, doublereal *thirds, integer *length,
+ integer *anchor, integer *point, doublereal *c1, doublereal *c2,
+ integer *maxfunc, integer *maxdeep, integer *maxdim, integer *n,
+ FILE *logfile, doublereal *fmax, integer jones)
+{
+ /* System generated locals */
+ integer f_dim1, f_offset, c_dim1, c_offset, length_dim1, length_offset,
+ i__1, i__2, i__3;
+ doublereal d__1, d__2;
+
+ /* Local variables */
+ doublereal a[32], b[32];
+ integer i__, j, k, l;
+ doublereal x[32], sidelength;
+ integer help;
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
+/* +-----------------------------------------------------------------------+ */
+ /* Parameter adjustments */
+ --point;
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ ++anchor;
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+ c_dim1 = *maxfunc;
+ c_offset = 1 + c_dim1;
+ c__ -= c_offset;
+ --c2;
+ --c1;
+
+ /* Function Body */
+ i__1 = *free - 1;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (f[i__ + (f_dim1 << 1)] > 0.) {
+/* +-----------------------------------------------------------------------+ */
+/* | Get the maximum side length of the hyper rectangle and then set the | */
+/* | new side length to this lengths times the growth factor. | */
+/* +-----------------------------------------------------------------------+ */
+ help = direct_dirgetmaxdeep_(&i__, &length[length_offset], maxfunc, n);
+ sidelength = thirds[help] * 2.;
+/* +-----------------------------------------------------------------------+ */
+/* | Set the Center and the upper and lower bounds of the rectangles. | */
+/* +-----------------------------------------------------------------------+ */
+ i__2 = *n;
+ for (j = 1; j <= i__2; ++j) {
+ sidelength = thirds[length[i__ + j * length_dim1]];
+ a[j - 1] = c__[i__ + j * c_dim1] - sidelength;
+ b[j - 1] = c__[i__ + j * c_dim1] + sidelength;
+/* L20: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | The function value is reset to 'Inf', since it may have been changed | */
+/* | in an earlier iteration and now the feasible point which was close | */
+/* | is not close anymore (since the hyper rectangle surrounding the | */
+/* | current point may have shrunk). | */
+/* +-----------------------------------------------------------------------+ */
+ f[i__ + f_dim1] = 1e6f;
+ f[i__ + (f_dim1 << 1)] = 2.;
+/* +-----------------------------------------------------------------------+ */
+/* | Check if any feasible point is near this infeasible point. | */
+/* +-----------------------------------------------------------------------+ */
+ i__2 = *free - 1;
+ for (k = 1; k <= i__2; ++k) {
+/* +-----------------------------------------------------------------------+ */
+/* | If the point k is feasible, check if it is near. | */
+/* +-----------------------------------------------------------------------+ */
+ if (f[k + (f_dim1 << 1)] == 0.) {
+/* +-----------------------------------------------------------------------+ */
+/* | Copy the coordinates of the point k into x. | */
+/* +-----------------------------------------------------------------------+ */
+ i__3 = *n;
+ for (l = 1; l <= i__3; ++l) {
+ x[l - 1] = c__[k + l * c_dim1];
+/* L40: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Check if the point k is near the infeasible point, if so, replace the | */
+/* | value at */
+/* +-----------------------------------------------------------------------+ */
+ if (isinbox_(x, a, b, n, &c__32) == 1) {
+/* Computing MIN */
+ d__1 = f[i__ + f_dim1], d__2 = f[k + f_dim1];
+ f[i__ + f_dim1] = MIN(d__1,d__2);
+ f[i__ + (f_dim1 << 1)] = 1.;
+ }
+ }
+/* L30: */
+ }
+ if (f[i__ + (f_dim1 << 1)] == 1.) {
+ f[i__ + f_dim1] += (d__1 = f[i__ + f_dim1], fabs(d__1)) *
+ 1e-6f;
+ i__2 = *n;
+ for (l = 1; l <= i__2; ++l) {
+ x[l - 1] = c__[i__ + l * c_dim1] * c1[l] + c__[i__ + l *
+ c_dim1] * c2[l];
+/* L200: */
+ }
+ dirresortlist_(&i__, &anchor[-1], &f[f_offset], &point[1], &
+ length[length_offset], n, maxfunc, maxdim, maxdeep,
+ logfile, jones);
+ } else {
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 | */
+/* | Replaced fixed value for infeasible points with maximum value found, | */
+/* | increased by 1. | */
+/* +-----------------------------------------------------------------------+ */
+ if (! (*fmax == f[i__ + f_dim1])) {
+/* Computing MAX */
+ d__1 = *fmax + 1., d__2 = f[i__ + f_dim1];
+ f[i__ + f_dim1] = MAX(d__1,d__2);
+ }
+ }
+ }
+/* L10: */
+ }
+/* L1000: */
+} /* dirreplaceinf_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRInsert | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ static void dirinsert_(integer *start, integer *ins, integer *point,
+ doublereal *f, integer *maxfunc)
+{
+ /* System generated locals */
+ integer f_dim1, f_offset, i__1;
+
+ /* Local variables */
+ integer i__, help;
+
+/* JG 09/17/00 Rewrote this routine. */
+/* DO 10,i = 1,maxfunc */
+/* IF (f(ins,1) .LT. f(point(start),1)) THEN */
+/* help = point(start) */
+/* point(start) = ins */
+/* point(ins) = help */
+/* GOTO 20 */
+/* END IF */
+/* IF (point(start) .EQ. 0) THEN */
+/* point(start) = ins */
+/* point(ins) = 0 */
+/* GOTO 20 */
+/* END IF */
+/* start = point(start) */
+/* 10 CONTINUE */
+/* 20 END */
+ /* Parameter adjustments */
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ --point;
+
+ /* Function Body */
+ i__1 = *maxfunc;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (point[*start] == 0) {
+ point[*start] = *ins;
+ point[*ins] = 0;
+ return;
+ } else if (f[*ins + f_dim1] < f[point[*start] + f_dim1]) {
+ help = point[*start];
+ point[*start] = *ins;
+ point[*ins] = help;
+ return;
+ }
+ *start = point[*start];
+/* L10: */
+ }
+} /* dirinsert_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRInsertList | */
+/* | Changed 02-24-2000 | */
+/* | Got rid of the distinction between feasible and infeasible points| */
+/* | I could do this since infeasible points get set to a high | */
+/* | function value, which may be replaced by a function value of a | */
+/* | nearby function at the end of the main loop. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirinsertlist_(integer *new__, integer *anchor, integer *
+ point, doublereal *f, integer *maxi, integer *length, integer *
+ maxfunc, integer *maxdeep, integer *n, integer *samp,
+ integer jones)
+{
+ /* System generated locals */
+ integer length_dim1, length_offset, f_dim1, f_offset, i__1;
+
+ /* Local variables */
+ integer j;
+ integer pos;
+ integer pos1, pos2, deep;
+
+/* JG 09/24/00 Changed this to Getlevel */
+ /* Parameter adjustments */
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ --point;
+ ++anchor;
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+
+ /* Function Body */
+ i__1 = *maxi;
+ for (j = 1; j <= i__1; ++j) {
+ pos1 = *new__;
+ pos2 = point[pos1];
+ *new__ = point[pos2];
+/* JG 09/24/00 Changed this to Getlevel */
+/* deep = DIRGetMaxdeep(pos1,length,maxfunc,n) */
+ deep = direct_dirgetlevel_(&pos1, &length[length_offset], maxfunc, n, jones);
+ if (anchor[deep] == 0) {
+ if (f[pos2 + f_dim1] < f[pos1 + f_dim1]) {
+ anchor[deep] = pos2;
+ point[pos2] = pos1;
+ point[pos1] = 0;
+ } else {
+ anchor[deep] = pos1;
+ point[pos2] = 0;
+ }
+ } else {
+ pos = anchor[deep];
+ if (f[pos2 + f_dim1] < f[pos1 + f_dim1]) {
+ if (f[pos2 + f_dim1] < f[pos + f_dim1]) {
+ anchor[deep] = pos2;
+/* JG 08/30/00 Fixed bug. Sorting was not correct when */
+/* f(pos2,1) < f(pos1,1) < f(pos,1) */
+ if (f[pos1 + f_dim1] < f[pos + f_dim1]) {
+ point[pos2] = pos1;
+ point[pos1] = pos;
+ } else {
+ point[pos2] = pos;
+ dirinsert_(&pos, &pos1, &point[1], &f[f_offset],
+ maxfunc);
+ }
+ } else {
+ dirinsert_(&pos, &pos2, &point[1], &f[f_offset], maxfunc);
+ dirinsert_(&pos, &pos1, &point[1], &f[f_offset], maxfunc);
+ }
+ } else {
+ if (f[pos1 + f_dim1] < f[pos + f_dim1]) {
+/* JG 08/30/00 Fixed bug. Sorting was not correct when */
+/* f(pos1,1) < f(pos2,1) < f(pos,1) */
+ anchor[deep] = pos1;
+ if (f[pos + f_dim1] < f[pos2 + f_dim1]) {
+ point[pos1] = pos;
+ dirinsert_(&pos, &pos2, &point[1], &f[f_offset],
+ maxfunc);
+ } else {
+ point[pos1] = pos2;
+ point[pos2] = pos;
+ }
+ } else {
+ dirinsert_(&pos, &pos1, &point[1], &f[f_offset], maxfunc);
+ dirinsert_(&pos, &pos2, &point[1], &f[f_offset], maxfunc);
+ }
+ }
+ }
+/* L10: */
+ }
+/* JG 09/24/00 Changed this to Getlevel */
+/* deep = DIRGetMaxdeep(samp,length,maxfunc,n) */
+ deep = direct_dirgetlevel_(samp, &length[length_offset], maxfunc, n, jones);
+ pos = anchor[deep];
+ if (f[*samp + f_dim1] < f[pos + f_dim1]) {
+ anchor[deep] = *samp;
+ point[*samp] = pos;
+ } else {
+ dirinsert_(&pos, samp, &point[1], &f[f_offset], maxfunc);
+ }
+} /* dirinsertlist_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRInsertList2 (Old way to do it.) | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ static void dirinsertlist_2__(integer *start, integer *j, integer *k,
+ integer *list2, doublereal *w, integer *maxi, integer *n)
+{
+ /* System generated locals */
+ integer list2_dim1, list2_offset, i__1;
+
+ /* Local variables */
+ integer i__, pos;
+
+ /* Parameter adjustments */
+ --w;
+ list2_dim1 = *n;
+ list2_offset = 1 + list2_dim1;
+ list2 -= list2_offset;
+
+ /* Function Body */
+ pos = *start;
+ if (*start == 0) {
+ list2[*j + list2_dim1] = 0;
+ *start = *j;
+ goto L50;
+ }
+ if (w[*start] > w[*j]) {
+ list2[*j + list2_dim1] = *start;
+ *start = *j;
+ } else {
+ i__1 = *maxi;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (list2[pos + list2_dim1] == 0) {
+ list2[*j + list2_dim1] = 0;
+ list2[pos + list2_dim1] = *j;
+ goto L50;
+ } else {
+ if (w[*j] < w[list2[pos + list2_dim1]]) {
+ list2[*j + list2_dim1] = list2[pos + list2_dim1];
+ list2[pos + list2_dim1] = *j;
+ goto L50;
+ }
+ }
+ pos = list2[pos + list2_dim1];
+/* L10: */
+ }
+ }
+L50:
+ list2[*j + (list2_dim1 << 1)] = *k;
+} /* dirinsertlist_2__ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRSearchmin | */
+/* | Search for the minimum in the list. ! */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ static void dirsearchmin_(integer *start, integer *list2, integer *
+ pos, integer *k, integer *n)
+{
+ /* System generated locals */
+ integer list2_dim1, list2_offset;
+
+ /* Parameter adjustments */
+ list2_dim1 = *n;
+ list2_offset = 1 + list2_dim1;
+ list2 -= list2_offset;
+
+ /* Function Body */
+ *k = *start;
+ *pos = list2[*start + (list2_dim1 << 1)];
+ *start = list2[*start + list2_dim1];
+} /* dirsearchmin_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRSamplepoints | */
+/* | Subroutine to sample the new points. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirsamplepoints_(doublereal *c__, integer *arrayi,
+ doublereal *delta, integer *sample, integer *start, integer *length,
+ FILE *logfile, doublereal *f, integer *free,
+ integer *maxi, integer *point, doublereal *x, doublereal *l,
+ doublereal *fmin, integer *minpos, doublereal *u, integer *n,
+ integer *maxfunc, integer *maxdeep, integer *oops)
+{
+ /* System generated locals */
+ integer length_dim1, length_offset, c_dim1, c_offset, f_dim1, f_offset,
+ i__1, i__2;
+
+ /* Local variables */
+ integer j, k, pos;
+
+
+ /* Parameter adjustments */
+ --u;
+ --l;
+ --x;
+ --arrayi;
+ --point;
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+ c_dim1 = *maxfunc;
+ c_offset = 1 + c_dim1;
+ c__ -= c_offset;
+
+ /* Function Body */
+ *oops = 0;
+ pos = *free;
+ *start = *free;
+ i__1 = *maxi + *maxi;
+ for (k = 1; k <= i__1; ++k) {
+ i__2 = *n;
+ for (j = 1; j <= i__2; ++j) {
+ length[*free + j * length_dim1] = length[*sample + j *
+ length_dim1];
+ c__[*free + j * c_dim1] = c__[*sample + j * c_dim1];
+/* L20: */
+ }
+ pos = *free;
+ *free = point[*free];
+ if (*free == 0) {
+ if (logfile)
+ fprintf(logfile, "Error, no more free positions! "
+ "Increase maxfunc!\n");
+ *oops = 1;
+ return;
+ }
+/* L10: */
+ }
+ point[pos] = 0;
+ pos = *start;
+ i__1 = *maxi;
+ for (j = 1; j <= i__1; ++j) {
+ c__[pos + arrayi[j] * c_dim1] = c__[*sample + arrayi[j] * c_dim1] + *
+ delta;
+ pos = point[pos];
+ c__[pos + arrayi[j] * c_dim1] = c__[*sample + arrayi[j] * c_dim1] - *
+ delta;
+ pos = point[pos];
+/* L30: */
+ }
+ ASRT(pos <= 0);
+} /* dirsamplepoints_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRDivide | */
+/* | Subroutine to divide the hyper rectangles according to the rules. | */
+/* | Changed 02-24-2000 | */
+/* | Replaced if statement by min (line 367) | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirdivide_(integer *new__, integer *currentlength,
+ integer *length, integer *point, integer *arrayi, integer *sample,
+ integer *list2, doublereal *w, integer *maxi, doublereal *f, integer *
+ maxfunc, integer *maxdeep, integer *n)
+{
+ /* System generated locals */
+ integer length_dim1, length_offset, list2_dim1, list2_offset, f_dim1,
+ f_offset, i__1, i__2;
+ doublereal d__1, d__2;
+
+ /* Local variables */
+ integer i__, j, k, pos, pos2;
+ integer start;
+
+
+ /* Parameter adjustments */
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ --point;
+ --w;
+ list2_dim1 = *n;
+ list2_offset = 1 + list2_dim1;
+ list2 -= list2_offset;
+ --arrayi;
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+
+ /* Function Body */
+ start = 0;
+ pos = *new__;
+ i__1 = *maxi;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ j = arrayi[i__];
+ w[j] = f[pos + f_dim1];
+ k = pos;
+ pos = point[pos];
+/* Computing MIN */
+ d__1 = f[pos + f_dim1], d__2 = w[j];
+ w[j] = MIN(d__1,d__2);
+ pos = point[pos];
+ dirinsertlist_2__(&start, &j, &k, &list2[list2_offset], &w[1], maxi,
+ n);
+/* L10: */
+ }
+ ASRT(pos <= 0);
+ i__1 = *maxi;
+ for (j = 1; j <= i__1; ++j) {
+ dirsearchmin_(&start, &list2[list2_offset], &pos, &k, n);
+ pos2 = start;
+ length[*sample + k * length_dim1] = *currentlength + 1;
+ i__2 = *maxi - j + 1;
+ for (i__ = 1; i__ <= i__2; ++i__) {
+ length[pos + k * length_dim1] = *currentlength + 1;
+ pos = point[pos];
+ length[pos + k * length_dim1] = *currentlength + 1;
+/* JG 07/10/01 pos2 = 0 at the end of the 30-loop. Since we end */
+/* the loop now, we do not need to reassign pos and pos2. */
+ if (pos2 > 0) {
+ pos = list2[pos2 + (list2_dim1 << 1)];
+ pos2 = list2[pos2 + list2_dim1];
+ }
+/* L30: */
+ }
+/* L20: */
+ }
+} /* dirdivide_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | | */
+/* | SUBROUTINE DIRINFCN | */
+/* | | */
+/* | Subroutine DIRinfcn unscales the variable x for use in the | */
+/* | user-supplied function evaluation subroutine fcn. After fcn returns | */
+/* | to DIRinfcn, DIRinfcn then rescales x for use by DIRECT. | */
+/* | | */
+/* | On entry | */
+/* | | */
+/* | fcn -- The argument containing the name of the user-supplied | */
+/* | subroutine that returns values for the function to be | */
+/* | minimized. | */
+/* | | */
+/* | x -- A double-precision vector of length n. The point at | */
+/* | which the derivative is to be evaluated. | */
+/* | | */
+/* | xs1 -- A double-precision vector of length n. Used for | */
+/* | scaling and unscaling the vector x by DIRinfcn. | */
+/* | | */
+/* | xs2 -- A double-precision vector of length n. Used for | */
+/* | scaling and unscaling the vector x by DIRinfcn. | */
+/* | | */
+/* | n -- An integer. The dimension of the problem. | */
+/* | kret -- An Integer. If kret = 1, the point is infeasible, | */
+/* | kret = -1, bad problem set up, | */
+/* | kret = 0, feasible. | */
+/* | | */
+/* | On return | */
+/* | | */
+/* | f -- A double-precision scalar. | */
+/* | | */
+/* | Subroutines and Functions | */
+/* | | */
+/* | The subroutine whose name is passed through the argument fcn. | */
+/* | | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirinfcn_(fp fcn, doublereal *x, doublereal *c1,
+ doublereal *c2, integer *n, doublereal *f, integer *flag__,
+ void *fcn_data)
+{
+ /* System generated locals */
+ integer i__1;
+
+ /* Local variables */
+ integer i__;
+
+/* +-----------------------------------------------------------------------+ */
+/* | Variables to pass user defined data to the function to be optimized. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | Unscale the variable x. | */
+/* +-----------------------------------------------------------------------+ */
+ /* Parameter adjustments */
+ --c2;
+ --c1;
+ --x;
+
+ /* Function Body */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ x[i__] = (x[i__] + c2[i__]) * c1[i__];
+/* L20: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Call the function-evaluation subroutine fcn. | */
+/* +-----------------------------------------------------------------------+ */
+ *flag__ = 0;
+ *f = fcn(*n, &x[1], flag__, fcn_data);
+/* +-----------------------------------------------------------------------+ */
+/* | Rescale the variable x. | */
+/* +-----------------------------------------------------------------------+ */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ x[i__] = x[i__] / c1[i__] - c2[i__];
+/* L30: */
+ }
+} /* dirinfcn_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRGet_I | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirget_i__(integer *length, integer *pos, integer *
+ arrayi, integer *maxi, integer *n, integer *maxfunc)
+{
+ /* System generated locals */
+ integer length_dim1, length_offset, i__1;
+
+ /* Local variables */
+ integer i__, j, help;
+
+ /* Parameter adjustments */
+ --arrayi;
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+
+ /* Function Body */
+ j = 1;
+ help = length[*pos + length_dim1];
+ i__1 = *n;
+ for (i__ = 2; i__ <= i__1; ++i__) {
+ if (length[*pos + i__ * length_dim1] < help) {
+ help = length[*pos + i__ * length_dim1];
+ }
+/* L10: */
+ }
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (length[*pos + i__ * length_dim1] == help) {
+ arrayi[j] = i__;
+ ++j;
+ }
+/* L20: */
+ }
+ *maxi = j - 1;
+} /* dirget_i__ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRInit | */
+/* | Initialise all needed variables and do the first run of the | */
+/* | algorithm. | */
+/* | Changed 02/24/2000 | */
+/* | Changed fcn Double precision to fcn external! | */
+/* | Changed 09/15/2000 | */
+/* | Added distinction between Jones way to characterize rectangles | */
+/* | and our way. Common variable JONES controls which way we use. | */
+/* | JONES = 0 Jones way (Distance from midpoint to corner) | */
+/* | JONES = 1 Our way (Length of longest side) | */
+/* | Changed 09/24/00 | */
+/* | Added array levels. Levels contain the values to characterize | */
+/* | the hyperrectangles. | */
+/* | Changed 01/22/01 | */
+/* | Added variable fmax to keep track of maximum value found. | */
+/* | Added variable Ifeasiblef to keep track if feasibel point has | */
+/* | been found. | */
+/* | Changed 01/23/01 | */
+/* | Added variable Ierror to keep track of errors. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirinit_(doublereal *f, fp fcn, doublereal *c__,
+ integer *length, integer *actdeep, integer *point, integer *anchor,
+ integer *free, FILE *logfile, integer *arrayi,
+ integer *maxi, integer *list2, doublereal *w, doublereal *x,
+ doublereal *l, doublereal *u, doublereal *fmin, integer *minpos,
+ doublereal *thirds, doublereal *levels, integer *maxfunc, integer *
+ maxdeep, integer *n, integer *maxor, doublereal *fmax, integer *
+ ifeasiblef, integer *iinfeasible, integer *ierror, void *fcndata,
+ integer jones)
+{
+ /* System generated locals */
+ integer f_dim1, f_offset, c_dim1, c_offset, length_dim1, length_offset,
+ list2_dim1, list2_offset, i__1, i__2;
+
+ /* Local variables */
+ integer i__, j;
+ integer new__, help, oops;
+ doublereal help2, delta;
+ doublereal costmin;
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
+/* +-----------------------------------------------------------------------+ */
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable Ifeasiblef to keep track if feasibel point | */
+/* | has been found. | */
+/* | JG 01/23/01 Added variable Ierror to keep track of errors. | */
+/* | JG 03/09/01 Added IInfeasible to keep track if an infeasible point has| */
+/* | been found. | */
+/* +-----------------------------------------------------------------------+ */
+/* JG 09/15/00 Added variable JONES (see above) */
+/* +-----------------------------------------------------------------------+ */
+/* | Variables to pass user defined data to the function to be optimized. | */
+/* +-----------------------------------------------------------------------+ */
+ /* Parameter adjustments */
+ --point;
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ ++anchor;
+ --u;
+ --l;
+ --x;
+ --w;
+ list2_dim1 = *maxor;
+ list2_offset = 1 + list2_dim1;
+ list2 -= list2_offset;
+ --arrayi;
+ length_dim1 = *maxfunc;
+ length_offset = 1 + length_dim1;
+ length -= length_offset;
+ c_dim1 = *maxfunc;
+ c_offset = 1 + c_dim1;
+ c__ -= c_offset;
+
+ /* Function Body */
+ *fmin = 1e20;
+ costmin = *fmin;
+/* JG 09/15/00 If Jones way of characterising rectangles is used, */
+/* initialise thirds to reflect this. */
+ if (jones == 0) {
+ i__1 = *n - 1;
+ for (j = 0; j <= i__1; ++j) {
+ w[j + 1] = sqrt(*n - j + j / 9.) * .5;
+/* L5: */
+ }
+ help2 = 1.;
+ i__1 = *maxdeep / *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ i__2 = *n - 1;
+ for (j = 0; j <= i__2; ++j) {
+ levels[(i__ - 1) * *n + j] = w[j + 1] / help2;
+/* L8: */
+ }
+ help2 *= 3.;
+/* L10: */
+ }
+ } else {
+/* JG 09/15/00 Initialiase levels to contain 1/j */
+ help2 = 3.;
+ i__1 = *maxdeep;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ levels[i__] = 1. / help2;
+ help2 *= 3.;
+/* L11: */
+ }
+ levels[0] = 1.;
+ }
+ help2 = 3.;
+ i__1 = *maxdeep;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ thirds[i__] = 1. / help2;
+ help2 *= 3.;
+/* L21: */
+ }
+ thirds[0] = 1.;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ c__[i__ * c_dim1 + 1] = .5;
+ x[i__] = .5;
+ length[i__ * length_dim1 + 1] = 0;
+/* L20: */
+ }
+ direct_dirinfcn_(fcn, &x[1], &l[1], &u[1], n, &f[f_dim1 + 1], &help, fcndata);
+ f[(f_dim1 << 1) + 1] = (doublereal) help;
+ *iinfeasible = help;
+ *fmax = f[f_dim1 + 1];
+/* 09/25/00 Added this */
+/* if (f(1,1) .ge. 1.E+6) then */
+ if (f[(f_dim1 << 1) + 1] > 0.) {
+ f[f_dim1 + 1] = 1e6;
+ *fmax = f[f_dim1 + 1];
+ *ifeasiblef = 1;
+ } else {
+ *ifeasiblef = 0;
+ }
+/* JG 09/25/00 Remove IF */
+ *fmin = f[f_dim1 + 1];
+ costmin = f[f_dim1 + 1];
+ *minpos = 1;
+ *actdeep = 2;
+ point[1] = 0;
+ *free = 2;
+ delta = thirds[1];
+ direct_dirget_i__(&length[length_offset], &c__1, &arrayi[1], maxi, n, maxfunc);
+ new__ = *free;
+ direct_dirsamplepoints_(&c__[c_offset], &arrayi[1], &delta, &c__1, &new__, &
+ length[length_offset], logfile, &f[f_offset], free, maxi, &
+ point[1], &x[1], &l[1], fmin, minpos, &u[1], n,
+ maxfunc, maxdeep, &oops);
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/23/01 Added error checking. | */
+/* +-----------------------------------------------------------------------+ */
+ if (oops > 0) {
+ *ierror = -4;
+ return;
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/22/01 Added variable to keep track of the maximum value found. | */
+/* | Added variable to keep track if feasible point was found. | */
+/* +-----------------------------------------------------------------------+ */
+ direct_dirsamplef_(&c__[c_offset], &arrayi[1], &delta, &c__1, &new__, &length[
+ length_offset], logfile, &f[f_offset], free, maxi, &point[
+ 1], fcn, &x[1], &l[1], fmin, minpos, &u[1], n, maxfunc,
+ maxdeep, &oops, fmax, ifeasiblef, iinfeasible, fcndata);
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/23/01 Added error checking. | */
+/* +-----------------------------------------------------------------------+ */
+ if (oops > 0) {
+ *ierror = -5;
+ return;
+ }
+ direct_dirdivide_(&new__, &c__0, &length[length_offset], &point[1], &arrayi[1], &
+ c__1, &list2[list2_offset], &w[1], maxi, &f[f_offset], maxfunc,
+ maxdeep, n);
+ direct_dirinsertlist_(&new__, &anchor[-1], &point[1], &f[f_offset], maxi, &
+ length[length_offset], maxfunc, maxdeep, n, &c__1, jones);
+} /* dirinit_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRInitList | */
+/* | Initialise the list. | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirinitlist_(integer *anchor, integer *free, integer *
+ point, doublereal *f, integer *maxfunc, integer *maxdeep)
+{
+ /* System generated locals */
+ integer f_dim1, f_offset, i__1;
+
+ /* Local variables */
+ integer i__;
+
+/* f -- values of functions. */
+/* anchor -- anchors of lists with deep i */
+/* point -- lists */
+/* free -- first free position */
+ /* Parameter adjustments */
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ --point;
+ ++anchor;
+
+ /* Function Body */
+ i__1 = *maxdeep;
+ for (i__ = -1; i__ <= i__1; ++i__) {
+ anchor[i__] = 0;
+/* L10: */
+ }
+ i__1 = *maxfunc;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ f[i__ + f_dim1] = 0.;
+ f[i__ + (f_dim1 << 1)] = 0.;
+ point[i__] = i__ + 1;
+/* point(i) = 0 */
+/* L20: */
+ }
+ point[*maxfunc] = 0;
+ *free = 1;
+} /* dirinitlist_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRSort3 | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ static void dirsort3_(integer *pos1, integer *pos2, integer *pos3,
+ doublereal *f, integer *maxfunc)
+{
+ /* System generated locals */
+ integer f_dim1, f_offset;
+
+ /* Local variables */
+ integer help;
+
+ /* Parameter adjustments */
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+
+ /* Function Body */
+ if (f[*pos1 + f_dim1] < f[*pos2 + f_dim1]) {
+ if (f[*pos1 + f_dim1] < f[*pos3 + f_dim1]) {
+ if (f[*pos3 + f_dim1] < f[*pos2 + f_dim1]) {
+ help = *pos2;
+ *pos2 = *pos3;
+ *pos3 = help;
+ }
+ } else {
+ help = *pos1;
+ *pos1 = *pos3;
+ *pos3 = *pos2;
+ *pos2 = help;
+ }
+ } else {
+ if (f[*pos2 + f_dim1] < f[*pos3 + f_dim1]) {
+ if (f[*pos3 + f_dim1] < f[*pos1 + f_dim1]) {
+ help = *pos1;
+ *pos1 = *pos2;
+ *pos2 = *pos3;
+ *pos3 = help;
+ } else {
+ help = *pos1;
+ *pos1 = *pos2;
+ *pos2 = help;
+ }
+ } else {
+ help = *pos1;
+ *pos1 = *pos3;
+ *pos3 = help;
+ }
+ }
+} /* dirsort3_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | SUBROUTINE DIRInsert3 | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirinsert3_(integer *pos1, integer *pos2, integer *pos3,
+ integer *deep, integer *anchor, integer *point, integer *free,
+ doublereal *f, doublereal *fmin, integer *minpos, integer *maxfunc,
+ integer *maxdeep)
+{
+ /* System generated locals */
+ integer f_dim1, f_offset;
+
+ /* Local variables */
+ integer pos;
+
+ /* Parameter adjustments */
+ f_dim1 = *maxfunc;
+ f_offset = 1 + f_dim1;
+ f -= f_offset;
+ --point;
+ ++anchor;
+
+ /* Function Body */
+ dirsort3_(pos1, pos2, pos3, &f[f_offset], maxfunc);
+ if (anchor[*deep] == 0) {
+ anchor[*deep] = *pos1;
+ point[*pos1] = *pos2;
+ point[*pos2] = *pos3;
+ point[*pos3] = 0;
+ } else {
+ pos = anchor[*deep];
+ if (f[*pos1 + f_dim1] < f[pos + f_dim1]) {
+ anchor[*deep] = *pos1;
+ point[*pos1] = pos;
+ } else {
+ dirinsert_(&pos, pos1, &point[1], &f[f_offset], maxfunc);
+ }
+ dirinsert_(&pos, pos2, &point[1], &f[f_offset], maxfunc);
+ dirinsert_(&pos, pos3, &point[1], &f[f_offset], maxfunc);
+ }
+ if (f[*pos1 + f_dim1] < *fmin && f[*pos1 + (f_dim1 << 1)] == 0.) {
+ *fmin = f[*pos1 + f_dim1];
+ *minpos = *pos1;
+ }
+} /* dirinsert3_ */
+
+/* +-----------------------------------------------------------------------+ */
+/* | | */
+/* | SUBROUTINE DIRPREPRC | */
+/* | | */
+/* | Subroutine DIRpreprc uses an afine mapping to map the hyper-box given | */
+/* | by the constraints on the variable x onto the n-dimensional unit cube.| */
+/* | This mapping is done using the following equation: | */
+/* | | */
+/* | x(i)=x(i)/(u(i)-l(i))-l(i)/(u(i)-l(i)). | */
+/* | | */
+/* | DIRpreprc checks if the bounds l and u are well-defined. That is, if | */
+/* | | */
+/* | l(i) < u(i) forevery i. | */
+/* | | */
+/* | On entry | */
+/* | | */
+/* | u -- A double-precision vector of length n. The vector | */
+/* | containing the upper bounds for the n independent | */
+/* | variables. | */
+/* | | */
+/* | l -- A double-precision vector of length n. The vector | */
+/* | containing the lower bounds for the n independent | */
+/* | variables. | */
+/* | | */
+/* | n -- An integer. The dimension of the problem. | */
+/* | | */
+/* | On return | */
+/* | | */
+/* | xs1 -- A double-precision vector of length n, used for scaling | */
+/* | and unscaling the vector x. | */
+/* | | */
+/* | xs2 -- A double-precision vector of length n, used for scaling | */
+/* | and unscaling the vector x. | */
+/* | | */
+/* | | */
+/* | oops -- An integer. If an upper bound is less than a lower | */
+/* | bound or if the initial point is not in the | */
+/* | hyper-box oops is set to 1 and iffco terminates. | */
+/* | | */
+/* +-----------------------------------------------------------------------+ */
+/* Subroutine */ void direct_dirpreprc_(doublereal *u, doublereal *l, integer *n,
+ doublereal *xs1, doublereal *xs2, integer *oops)
+{
+ /* System generated locals */
+ integer i__1;
+
+ /* Local variables */
+ integer i__;
+ doublereal help;
+
+ /* Parameter adjustments */
+ --xs2;
+ --xs1;
+ --l;
+ --u;
+
+ /* Function Body */
+ *oops = 0;
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+/* +-----------------------------------------------------------------------+ */
+/* | Check if the hyper-box is well-defined. | */
+/* +-----------------------------------------------------------------------+ */
+ if (u[i__] <= l[i__]) {
+ *oops = 1;
+ return;
+ }
+/* L20: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | Scale the initial iterate so that it is in the unit cube. | */
+/* +-----------------------------------------------------------------------+ */
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ help = u[i__] - l[i__];
+ xs2[i__] = l[i__] / help;
+ xs1[i__] = help;
+/* L50: */
+ }
+} /* dirpreprc_ */
+
+/* Subroutine */ void direct_dirheader_(FILE *logfile, integer *version,
+ doublereal *x, integer *n, doublereal *eps, integer *maxf, integer *
+ maxt, doublereal *l, doublereal *u, integer *algmethod, integer *
+ maxfunc, integer *maxdeep, doublereal *fglobal, doublereal *fglper,
+ integer *ierror, doublereal *epsfix, integer *iepschange, doublereal *
+ volper, doublereal *sigmaper)
+{
+ /* System generated locals */
+ integer i__1;
+
+ /* Local variables */
+ integer imainver, i__, numerrors, isubsubver, ihelp, isubver;
+
+
+/* +-----------------------------------------------------------------------+ */
+/* | Variables to pass user defined data to the function to be optimized. | */
+/* +-----------------------------------------------------------------------+ */
+ /* Parameter adjustments */
+ --u;
+ --l;
+ --x;
+
+ /* Function Body */
+ if (logfile)
+ fprintf(logfile, "------------------- Log file ------------------\n");
+
+ numerrors = 0;
+ *ierror = 0;
+ imainver = *version / 100;
+ ihelp = *version - imainver * 100;
+ isubver = ihelp / 10;
+ ihelp -= isubver * 10;
+ isubsubver = ihelp;
+/* +-----------------------------------------------------------------------+ */
+/* | JG 01/13/01 Added check for epsilon. If epsilon is smaller 0, we use | */
+/* | the update formula from Jones. We then set the flag | */
+/* | iepschange to 1, and store the absolute value of eps in | */
+/* | epsfix. epsilon is then changed after each iteration. | */
+/* +-----------------------------------------------------------------------+ */
+ if (*eps < 0.) {
+ *iepschange = 1;
+ *epsfix = -(*eps);
+ *eps = -(*eps);
+ } else {
+ *iepschange = 0;
+ *epsfix = 1e100;
+ }
+
+/* +-----------------------------------------------------------------------+ */
+/* | JG 07/16/01 Removed printout of contents in cdata(1). | */
+/* +-----------------------------------------------------------------------+ */
+/* write(logfile,*) cdata(1) */
+
+ if (logfile) {
+ fprintf(logfile, "DIRECT Version %d.%d.%d\n"
+ " Problem dimension n: %d\n"
+ " Eps value: %e\n"
+ " Maximum number of f-evaluations (maxf): %d\n"
+ " Maximum number of iterations (MaxT): %d\n"
+ " Value of f_global: %e\n"
+ " Global percentage wanted: %e\n"
+ " Volume percentage wanted: %e\n"
+ " Measure percentage wanted: %e\n",
+ imainver, isubver, isubsubver, *n, *eps, *maxf, *maxt,
+ *fglobal, *fglper, *volper, *sigmaper);
+ fprintf(logfile, *iepschange == 1
+ ? "Epsilon is changed using the Jones formula.\n"
+ : "Epsilon is constant.\n");
+ fprintf(logfile, *algmethod == 0
+ ? "Jones original DIRECT algorithm is used.\n"
+ : "Our modification of the DIRECT algorithm is used.\n");
+ }
+
+ i__1 = *n;
+ for (i__ = 1; i__ <= i__1; ++i__) {
+ if (u[i__] <= l[i__]) {
+ *ierror = -1;
+ if (logfile)
+ fprintf(logfile, "WARNING: bounds on variable x%d: "
+ "%g <= xi <= %g\n", i__, l[i__], u[i__]);
+ ++numerrors;
+ } else {
+ if (logfile)
+ fprintf(logfile, "Bounds on variable x%d: "
+ "%g <= xi <= %g\n", i__, l[i__], u[i__]);
+ }
+/* L1010: */
+ }
+/* +-----------------------------------------------------------------------+ */
+/* | If there are to many function evaluations or to many iteration, note | */
+/* | this and set the error flag accordingly. Note: If more than one error | */
+/* | occurred, we give out an extra message. | */
+/* +-----------------------------------------------------------------------+ */
+ if (*maxf + 20 > *maxfunc) {
+ if (logfile)
+ fprintf(logfile,
+"WARNING: The maximum number of function evaluations (%d) is higher than\n"
+" the constant maxfunc (%d). Increase maxfunc in subroutine DIRECT\n"
+" or decrease the maximum number of function evaluations.\n",
+ *maxf, *maxfunc);
+ ++numerrors;
+ *ierror = -2;
+ }
+ if (*ierror < 0) {
+ if (logfile) fprintf(logfile, "----------------------------------\n");
+ if (numerrors == 1) {
+ if (logfile)
+ fprintf(logfile, "WARNING: One error in the input!\n");
+ } else {
+ if (logfile)
+ fprintf(logfile, "WARNING: %d errors in the input!\n",
+ numerrors);
+ }
+ }
+ if (logfile) fprintf(logfile, "----------------------------------\n");
+ if (*ierror >= 0) {
+ if (logfile)
+ fprintf(logfile, "Iteration # of f-eval. fmin\n");
+ }
+/* L10005: */
+} /* dirheader_ */
+
+/* Subroutine */ void direct_dirsummary_(FILE *logfile, doublereal *x, doublereal *
+ l, doublereal *u, integer *n, doublereal *fmin, doublereal *fglobal,
+ integer *numfunc, integer *ierror)
+{
+ /* Local variables */
+ integer i__;
+
+ /* Parameter adjustments */
+ --u;
+ --l;
+ --x;
+
+ /* Function Body */
+ if (logfile) {
+ fprintf(logfile, "-----------------------Summary------------------\n"
+ "Final function value: %g\n"
+ "Number of function evaluations: %d\n", *fmin, *numfunc);
+ if (*fglobal > -1e99)
+ fprintf(logfile, "Final function value is within %g%% of global optimum\n", 100*(*fmin - *fglobal) / MAX(1.0, fabs(*fglobal)));
+ fprintf(logfile, "Index, final solution, x(i)-l(i), u(i)-x(i)\n");
+ for (i__ = 1; i__ <= *n; ++i__)
+ fprintf(logfile, "%d, %g, %g, %g\n", i__, x[i__],
+ x[i__]-l[i__], u[i__] - x[i__]);
+ fprintf(logfile, "-----------------------------------------------\n");
+
+ }
+} /* dirsummary_ */
+
+/* Subroutine */ void direct_dirmaxf_to_high1__(integer *maxf, integer *maxfunc,
+ FILE *logfile)
+{
+ if (logfile) {
+ fprintf(logfile,
+"The maximum number of function evaluations (%d)\n"
+"is higher than the constant maxfunc (%d). Increase maxfunc\n"
+"in the SUBROUTINE DIRECT or decrease the maximum number\n"
+ "of function evaluations.\n", *maxf, *maxfunc);
+ }
+} /* dirmaxf_to_high1__ */
+
+/* Subroutine */ void direct_dirmaxt_to_high1__(integer *maxt, integer *maxdeep,
+ FILE *logfile)
+{
+ if (logfile) {
+ fprintf(logfile,
+"The maximum number of iterations (%d) is higher than the constant\n"
+"maxdeep (%d). Increase maxdeep or decrease the number of iterations.\n",
+ *maxt, *maxdeep);
+ }
+} /* dirmaxt_to_high1__ */
+
--- /dev/null
+#ifndef DIRECT_INTERNAL_H
+#define DIRECT_INTERNAL_H
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "direct.h"
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif /* __cplusplus */
+
+typedef int integer;
+typedef double doublereal;
+typedef direct_objective_func fp;
+
+#define ASRT(c) if (!(c)) { fprintf(stderr, "DIRECT assertion failure at " __FILE__ ":%d -- " #c "\n", __LINE__); exit(EXIT_FAILURE); }
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
+/* DIRect.c */
+extern void direct_dirsamplef_(
+ doublereal *c__, integer *arrayi, doublereal
+ *delta, integer *sample, integer *new__, integer *length,
+ FILE *logfile, doublereal *f, integer *free, integer *maxi,
+ integer *point, fp fcn, doublereal *x, doublereal *l, doublereal *
+ fmin, integer *minpos, doublereal *u, integer *n, integer *maxfunc,
+ integer *maxdeep, integer *oops, doublereal *fmax, integer *
+ ifeasiblef, integer *iinfesiblef, void *fcndata);
+
+/* DIRsubrout.c */
+
+extern void direct_dirheader_(
+ FILE *logfile, integer *version,
+ doublereal *x, integer *n, doublereal *eps, integer *maxf, integer *
+ maxt, doublereal *l, doublereal *u, integer *algmethod, integer *
+ maxfunc, integer *maxdeep, doublereal *fglobal, doublereal *fglper,
+ integer *ierror, doublereal *epsfix, integer *iepschange, doublereal *
+ volper, doublereal *sigmaper);
+extern void direct_dirinit_(
+ doublereal *f, fp fcn, doublereal *c__,
+ integer *length, integer *actdeep, integer *point, integer *anchor,
+ integer *free, FILE *logfile, integer *arrayi,
+ integer *maxi, integer *list2, doublereal *w, doublereal *x,
+ doublereal *l, doublereal *u, doublereal *fmin, integer *minpos,
+ doublereal *thirds, doublereal *levels, integer *maxfunc, integer *
+ maxdeep, integer *n, integer *maxor, doublereal *fmax, integer *
+ ifeasiblef, integer *iinfeasible, integer *ierror, void *fcndata,
+ integer jones);
+extern void direct_dirinitlist_(
+ integer *anchor, integer *free, integer *
+ point, doublereal *f, integer *maxfunc, integer *maxdeep);
+extern void direct_dirpreprc_(doublereal *u, doublereal *l, integer *n,
+ doublereal *xs1, doublereal *xs2, integer *oops);
+extern void direct_dirchoose_(
+ integer *anchor, integer *s, integer *actdeep,
+ doublereal *f, doublereal *fmin, doublereal epsrel, doublereal epsabs, doublereal *thirds,
+ integer *maxpos, integer *length, integer *maxfunc, integer *maxdeep,
+ integer *maxdiv, integer *n, FILE *logfile,
+ integer *cheat, doublereal *kmax, integer *ifeasiblef, integer jones);
+extern void direct_dirdoubleinsert_(
+ integer *anchor, integer *s, integer *maxpos, integer *point,
+ doublereal *f, integer *maxdeep, integer *maxfunc,
+ integer *maxdiv, integer *ierror);
+extern integer direct_dirgetmaxdeep_(integer *pos, integer *length, integer *maxfunc,
+ integer *n);
+extern void direct_dirget_i__(
+ integer *length, integer *pos, integer *arrayi, integer *maxi,
+ integer *n, integer *maxfunc);
+extern void direct_dirsamplepoints_(
+ doublereal *c__, integer *arrayi,
+ doublereal *delta, integer *sample, integer *start, integer *length,
+ FILE *logfile, doublereal *f, integer *free,
+ integer *maxi, integer *point, doublereal *x, doublereal *l,
+ doublereal *fmin, integer *minpos, doublereal *u, integer *n,
+ integer *maxfunc, integer *maxdeep, integer *oops);
+extern void direct_dirdivide_(
+ integer *new__, integer *currentlength,
+ integer *length, integer *point, integer *arrayi, integer *sample,
+ integer *list2, doublereal *w, integer *maxi, doublereal *f,
+ integer *maxfunc, integer *maxdeep, integer *n);
+extern void direct_dirinsertlist_(
+ integer *new__, integer *anchor, integer *point, doublereal *f,
+ integer *maxi, integer *length, integer *maxfunc,
+ integer *maxdeep, integer *n, integer *samp, integer jones);
+extern void direct_dirreplaceinf_(
+ integer *free, integer *freeold,
+ doublereal *f, doublereal *c__, doublereal *thirds, integer *length,
+ integer *anchor, integer *point, doublereal *c1, doublereal *c2,
+ integer *maxfunc, integer *maxdeep, integer *maxdim, integer *n,
+ FILE *logfile, doublereal *fmax, integer jones);
+extern void direct_dirsummary_(
+ FILE *logfile, doublereal *x, doublereal *l, doublereal *u,
+ integer *n, doublereal *fmin, doublereal *fglobal,
+ integer *numfunc, integer *ierror);
+extern integer direct_dirgetlevel_(
+ integer *pos, integer *length,
+ integer *maxfunc, integer *n, integer jones);
+extern void direct_dirinfcn_(
+ fp fcn, doublereal *x, doublereal *c1,
+ doublereal *c2, integer *n, doublereal *f, integer *flag__,
+ void *fcn_data);
+
+/* DIRserial.c / DIRparallel.c */
+extern void direct_dirsamplef_(
+ doublereal *c__, integer *arrayi, doublereal
+ *delta, integer *sample, integer *new__, integer *length,
+ FILE *logfile, doublereal *f, integer *free, integer *maxi,
+ integer *point, fp fcn, doublereal *x, doublereal *l, doublereal *
+ fmin, integer *minpos, doublereal *u, integer *n, integer *maxfunc,
+ integer *maxdeep, integer *oops, doublereal *fmax, integer *
+ ifeasiblef, integer *iinfesiblef, void *fcn_data);
+
+/* DIRect.c */
+
+extern void direct_direct_(
+ fp fcn, doublereal *x, integer *n, doublereal *eps, doublereal epsabs,
+ integer *maxf, integer *maxt, doublereal *fmin, doublereal *l,
+ doublereal *u, integer *algmethod, integer *ierror, FILE *logfile,
+ doublereal *fglobal, doublereal *fglper, doublereal *volper,
+ doublereal *sigmaper, void *fcn_data);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif /* __cplusplus */
+
+#endif /* DIRECT_INTERNAL_H */
--- /dev/null
+#ifndef DIRECT_H
+#define DIRECT_H
+
+#include <math.h>
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif /* __cplusplus */
+
+typedef double (*direct_objective_func)(int n, const double *x,
+ int *undefined_flag,
+ void *data);
+
+typedef enum {
+ DIRECT_ORIGINAL, DIRECT_GABLONSKY
+} direct_algorithm;
+
+typedef enum {
+ DIRECT_INVALID_BOUNDS = -1,
+ DIRECT_MAXFEVAL_TOOBIG = -2,
+ DIRECT_INIT_FAILED = -3,
+ DIRECT_SAMPLEPOINTS_FAILED = -4,
+ DIRECT_SAMPLE_FAILED = -5,
+ DIRECT_MAXFEVAL_EXCEEDED = 1,
+ DIRECT_MAXITER_EXCEEDED = 2,
+ DIRECT_GLOBAL_FOUND = 3,
+ DIRECT_VOLTOL = 4,
+ DIRECT_SIGMATOL = 5,
+
+ DIRECT_OUT_OF_MEMORY = -100,
+ DIRECT_INVALID_ARGS = -101
+} direct_return_code;
+
+#define DIRECT_UNKNOWN_FGLOBAL (-HUGE_VAL)
+#define DIRECT_UNKNOWN_FGLOBAL_RELTOL (0.0)
+
+extern direct_return_code direct_optimize(
+ direct_objective_func f, void *f_data,
+ int dimension,
+ const double *lower_bounds, const double *upper_bounds,
+
+ double *x, double *fmin,
+
+ int max_feval, int max_iter,
+ double reltol, double abstol,
+ double volume_reltol, double sigma_reltol,
+
+ double fglobal,
+ double fglobal_reltol,
+
+ FILE *logfile,
+ direct_algorithm algorithm);
+
+#ifdef __cplusplus
+} /* extern "C" */
+#endif /* __cplusplus */
+
+#endif /* DIRECT_H */
--- /dev/null
+/* C-style API for DIRECT functions. SGJ (August 2007). */
+
+#include "direct-internal.h"
+
+/* Perform global minimization using (Gablonsky implementation of) DIRECT
+ algorithm. Arguments:
+
+ f, f_data: the objective function and any user data
+ -- the objective function f(n, x, undefined_flag, data) takes 4 args:
+ int n: the dimension, same as dimension arg. to direct_optimize
+ const double *x: array x[n] of point to evaluate
+ int *undefined_flag: set to 1 on return if x violates constraints
+ or don't touch otherwise
+ void *data: same as f_data passed to direct_optimize
+ return value = value of f(x)
+
+ dimension: the number of minimization variable dimensions
+ lower_bounds, upper_bounds: arrays of length dimension of variable bounds
+
+ x: an array of length dimension, set to optimum variables upon return
+ fmin: on return, set to minimum f value
+
+ max_feval, max_iter: maximum number of function evaluations & DIRECT iters
+ reltol, abstol: relative and absolute tolerances (0 if none)
+ volume_reltol: relative tolerance on hypercube volume (0 if none)
+ sigma_reltol: relative tolerance on hypercube "measure" (??) (0 if none)
+
+ fglobal: the global minimum of f, if known ahead of time
+ -- this is mainly for benchmarking, in most cases it
+ is not known and you should pass DIRECT_UNKNOWN_FGLOBAL
+ fglobal_reltol: relative tolerance on how close we should find fglobal
+ -- ignored if fglobal is DIRECT_UNKNOWN_FGLOBAL
+
+ logfile: an output file to write diagnostic info to (NULL for no I/O)
+
+ algorithm: whether to use the original DIRECT algorithm (DIRECT_ORIGINAL)
+ or Gablonsky's "improved" version (DIRECT_GABLONSKY)
+*/
+direct_return_code direct_optimize(
+ direct_objective_func f, void *f_data,
+ int dimension,
+ const double *lower_bounds, const double *upper_bounds,
+
+ double *x, double *fmin,
+
+ int max_feval, int max_iter,
+ double reltol, double abstol,
+ double volume_reltol, double sigma_reltol,
+
+ double fglobal,
+ double fglobal_reltol,
+
+ FILE *logfile,
+ direct_algorithm algorithm)
+{
+ integer algmethod = algorithm == DIRECT_GABLONSKY;
+ integer ierror;
+ doublereal *l, *u;
+ int i;
+
+ /* convert to percentages: */
+ volume_reltol *= 100;
+ sigma_reltol *= 100;
+ fglobal_reltol *= 100;
+
+ /* make sure these are ignored if <= 0 */
+ if (volume_reltol <= 0) volume_reltol = -1;
+ if (sigma_reltol <= 0) sigma_reltol = -1;
+
+ if (fglobal == DIRECT_UNKNOWN_FGLOBAL)
+ fglobal_reltol = DIRECT_UNKNOWN_FGLOBAL_RELTOL;
+
+ if (dimension < 1) return DIRECT_INVALID_ARGS;
+
+ l = (doublereal *) malloc(sizeof(doublereal) * dimension * 2);
+ if (!l) return DIRECT_OUT_OF_MEMORY;
+ u = l + dimension;
+ for (i = 0; i < dimension; ++i) {
+ l[i] = lower_bounds[i];
+ u[i] = upper_bounds[i];
+ }
+
+ direct_direct_(f, x, &dimension, &reltol, abstol,
+ &max_feval, &max_iter,
+ fmin,
+ l, u,
+ &algmethod,
+ &ierror,
+ logfile,
+ &fglobal, &fglobal_reltol,
+ &volume_reltol, &sigma_reltol,
+ f_data);
+
+ free(l);
+
+ return (direct_return_code) ierror;
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "direct.h"
+
+/* has two global minima at (0.09,-0.71) and (-0.09,0.71), plus
+ 4 additional local minima */
+static int cnt=0;
+double tst_obj(int n, const double *xy, int *undefined_flag, void *unused)
+{
+ double x, y, f;
+ x = xy[0];
+ y = xy[1];
+ f = ((x*x)*(4-2.1*(x*x)+((x*x)*(x*x))/3) + x*y + (y*y)*(-4+4*(y*y)));
+ printf("feval:, %d, %g, %g, %g\n", ++cnt, x,y, f);
+ return f;
+}
+
+int main(int argc, char **argv)
+{
+ int n = 2;
+ double x[2], l[2], u[2];
+ long int maxits = 0;
+ int info;
+ double fmin;
+
+ maxits = argc < 2 ? 100 : atoi(argv[1]);
+
+ l[0] = -3; l[1] = -3;
+ u[0] = 3; u[1] = 3;
+
+ info = direct_optimize(tst_obj, NULL, n, l, u, x, &fmin,
+ maxits, 500,
+ 0, 0, 0, 0, DIRECT_UNKNOWN_FGLOBAL, 0,
+ stdout, DIRECT_GABLONSKY);
+
+ printf("min f = %g at (%g,%g) after %d evals, return value %d\n",
+ fmin, x[0], x[1], cnt, info);
+
+ return EXIT_SUCCESS;
+}
--- /dev/null
+This code is modified from the original StoGO Global Optimization library
+by Madsen et al., downloaded from:
+
+ http://www2.imm.dtu.dk/~km/GlobOpt/opt.html
+
+It was modified to allow C-callable wrappers and some other niceties by
+Steven G. Johnson (stevenj@alum.mit.edu) in 2007.
+
+StoGO uses a gradient-based direct-search branch-and-bound algorithm,
+described in:
+
+S. Gudmundsson, "Parallel Global Optimization," M.Sc. Thesis, IMM,
+ Technical University of Denmark, 1998.
+
+K. Madsen, S. Zertchaninov, and A. Zilinskas, "Global Optimization
+ using Branch-and-Bound," Submitted to the Journal of Global
+ Optimization, 1998. [ does not seem to have been published? ]
+
+S. Zertchaninov and K. Madsen, "A C++ Programme for Global Optimization,"
+ IMM-REP-1998-04, Department of Mathematical Modelling,
+ Technical University of Denmark, DK-2800 Lyngby, Denmark, 1998.
--- /dev/null
+/*
+ Multi Dimensional Global Search.
+
+ Author: Steinn Gudmundsson
+ Email: steinng@hotmail.com
+
+ This program is supplied without any warranty whatsoever.
+
+ NB The RNGs seed should be initialized using some timer
+*/
+
+#include <iostream>
+#include <time.h>
+
+#include <iterator>
+#include <algorithm>
+#include <stack>
+
+#include "stogo_config.h"
+#include "global.h"
+#include "local.h"
+
+// Timer stuff
+time_t StartTime;
+
+double MacEpsilon ;
+int FC=0, GC=0 ;
+
+Global::Global(RTBox D, Pobj o, Pgrad g, GlobalParams P): Domain(D) {
+
+ dim=Domain.GetDim();
+ Objective=o;
+ Gradient=g;
+
+ // Initialize parameters
+ maxtime=P.maxtime;
+ maxeval = P.maxeval;
+ numeval = 0;
+ eps_cl=P.eps_cl; mu=P.mu; rshift=P.rshift;
+ det_pnts=P.det_pnts; rnd_pnts=P.rnd_pnts;
+ fbound=DBL_MAX;
+}
+
+#if 0 // not necessary; default copy is sufficient
+Global& Global::operator=(const Global &G) {
+ // Copy the problem info and parameter settings
+ Domain=G.Domain; Objective=G.Objective; Gradient=G.Gradient;
+ maxtime=G.maxtime;
+ maxeval = G.maxeval;
+ numeval = G.numeval;
+ eps_cl=G.eps_cl; mu=G.mu; rshift=G.rshift;
+ det_pnts=G.det_pnts; rnd_pnts=G.rnd_pnts;
+ return *this;
+}
+#endif
+
+void Global::FillRegular(RTBox SampleBox, RTBox box) {
+ // Generation of regular sampling points
+ double w;
+ int i, flag, dir;
+ Trial tmpTrial(dim);
+ RVector m(dim), x(dim);
+
+ if (det_pnts>0) {
+ // Add midpoint
+ box.Midpoint(m) ;
+ tmpTrial.xvals=m ; tmpTrial.objval=DBL_MAX ;
+ SampleBox.AddTrial(tmpTrial) ;
+ // Add the rest
+ i=1 ; flag=1 ; dir=0 ;
+ x=m ; w=box.Width(dir) ;
+ while (i<det_pnts) {
+ x(dir)=m(dir)+flag*rshift*w ;
+ tmpTrial.xvals=x ;
+ SampleBox.AddTrial(tmpTrial) ;
+ flag=-flag;
+ if (flag==1 && dir<dim) {
+ x(dir)=m(dir) ;
+ dir++ ;
+ w=box.Width(dir) ;
+ }
+ i++ ;
+ }
+ }
+}
+
+void Global::FillRandom(RTBox SampleBox, RTBox box) {
+ // Generation of stochastic sampling points
+ Trial tmpTrial(dim);
+
+ tmpTrial.objval=DBL_MAX;
+ for (int i=1 ; i<=rnd_pnts ; i++) {
+ for (int dir=0 ; dir<dim ; dir++)
+ tmpTrial.xvals(dir) =
+ box.lb(dir)+(box.ub(dir)-box.lb(dir))*(double(rand())/RAND_MAX) ;
+ SampleBox.AddTrial(tmpTrial) ;
+ }
+}
+
+double Global::NewtonTest(RTBox box, int axis, RCRVector x_av, int *noutside) {
+ // Perform the Newton test
+
+ int info,nout=0;
+ Trial tmpTrial(dim);
+ TBox SampleBox(dim) ;
+ double maxgrad=0 ;
+
+ // Create sampling points
+ FillRegular(SampleBox, box);
+ FillRandom(SampleBox, box);
+
+ // Perform the actual sampling
+ while ( !SampleBox.EmptyBox() ) {
+ SampleBox.RemoveTrial(tmpTrial) ;
+ info = local(tmpTrial, box, Domain, eps_cl, &maxgrad, *this,
+ axis, x_av) ;
+ // What should we do when info=LS_Unstable?
+ if (info == LS_Out)
+ nout++;
+ if (info == LS_New ) {
+ box.AddTrial(tmpTrial) ;
+
+ if (tmpTrial.objval<=fbound+mu && tmpTrial.objval<=box.fmin+mu) {
+ cout << "Found a candidate, x=" << tmpTrial.xvals;
+ cout << " F=" <<tmpTrial.objval << " FC=" << FC << endl;
+ SolSet.push_back(tmpTrial);
+ }
+#ifdef GS_DEBUG
+ cout << "Found a stationary point, X= " << tmpTrial.xvals;
+ cout <<" objval=" << tmpTrial.objval << endl;
+#endif
+ }
+
+ if (!InTime())
+ break;
+ }
+ *noutside=nout;
+ return maxgrad;
+}
+
+void Global::ReduceOrSubdivide(RTBox box, int axis, RCRVector x_av) {
+ TBox B1(dim), B2(dim);
+ Trial tmpTrial(dim);
+ double maxgrad;
+ int ns,nout;
+
+ // Monotonicity test has not been implemented yet
+ maxgrad=NewtonTest(box, axis, x_av, &nout);
+ ns=box.NStationary() ;
+ if (ns==0) {
+ // All iterates outside
+ // NB result=Intersection(B,boundary(Domain))
+ Garbage.push(box) ;
+ }
+ else
+ if (ns==1 && nout==0) {
+ // All iterates converge to same point
+ Garbage.push(box) ;
+ }
+ else
+ if ( (ns>1) && (box.LowerBound(maxgrad)>fbound) ) {
+ // Several stationary points found and lower bound > fbound
+ Garbage.push(box) ;
+ }
+ else {
+ // Subdivision
+ B1.ClearBox() ; B2.ClearBox() ;
+ box.split(B1,B2) ;
+ CandSet.push(B1) ; CandSet.push(B2) ;
+ }
+
+ // Update fbound
+ if (box.fmin < fbound) {
+ fbound=box.fmin ;
+#ifdef GS_DEBUG
+ cout <<"*** Improving fbound, fbound=" << fbound << endl;
+#endif
+ }
+}
+
+void Global::Search(int axis, RCRVector x_av){
+ Trial tmpTrial(dim) ;
+ TBox box(dim), B1(dim), B2(dim);
+ RVector m(dim), x(dim);
+ int inner_iter, outer_iter;
+
+ MacEpsilon=eps(); // Get machine precision
+ if (det_pnts>2*dim+1) {
+ det_pnts=2*dim+1;
+ cout << "Warning: Reducing det_pnts to " << det_pnts << endl;
+ }
+
+ // Initialize timer
+ time(&StartTime);
+
+ // Clear priority_queues
+ while (!Garbage.empty())
+ Garbage.pop();
+ while (!CandSet.empty())
+ CandSet.pop();
+
+ box=Domain;
+ CandSet.push(box);
+ int done=0 ; outer_iter=0 ;
+
+ while (!done) {
+ outer_iter++ ;
+
+ // Inner loop
+ inner_iter=0 ;
+ while (!CandSet.empty()) {
+ inner_iter++ ;
+ // Get best box from Candidate set
+ box=CandSet.top() ; CandSet.pop() ;
+
+#ifdef GS_DEBUG
+ cout << "Iteration..." << inner_iter << " #CS=" << CandSet.size()+1 ;
+ cout << " Processing " << box.NStationary() << " trials in the box " <<box;
+#endif
+ ReduceOrSubdivide(box, axis, x_av);
+
+ if (!InTime()) {
+ done=TRUE;
+ cout << "The program has run out of time or function evaluations\n";
+ break;
+ }
+
+ } // inner while-loop
+ cout << endl << "*** Inner loop completed ***" << endl ;
+
+ // Reduce SolSet if necessary
+ SolSet.erase(remove_if(SolSet.begin(), SolSet.end(),
+ TrialGT(fbound+mu)),SolSet.end());
+ if (InTime()) {
+ cout << "Current set of minimizers (" << SolSet.size() << ")" << endl ;
+ DispMinimizers() ;
+
+ while (!Garbage.empty()) {
+ box=Garbage.top() ;
+ Garbage.pop() ;
+ // Split box
+ B1.ClearBox() ; B2.ClearBox() ;
+ box.split(B1,B2) ;
+ // Add boxes to Candidate set
+ CandSet.push(B1) ; CandSet.push(B2) ;
+ }
+ }
+ } // Outer while-loop
+
+ cout << "Number of outer iterations : " << outer_iter << endl;
+ cout << "Number of unexplored boxes : " << CandSet.size() << endl;
+ cout << "Number of boxes in garbage : " << Garbage.size() << endl;
+ cout << "Number of elements in SolSet : " << SolSet.size() << endl;
+ cout << "Number of function evaluations : " << FC << endl;
+ cout << "Number of gradient evaluations : " << GC << endl;
+
+ if (axis != -1) {
+ // Return minimizer when doing the AV method
+ tmpTrial=SolSet.back();
+ x_av(axis)=tmpTrial.xvals(0);
+ }
+}
+
+/************* Various utility functions ****************/
+long int Global::GetTime()
+{
+ time_t ctime; time(&ctime);
+ return (long int)difftime(ctime,StartTime);
+}
+
+bool Global::InTime()
+{
+ return (!maxtime || GetTime()<maxtime) && (!maxeval || numeval<maxeval);
+}
+
+double Global::GetMinValue() {
+ return fbound;
+}
+
+void Global::SetMinValue(double new_fb) {
+ fbound=new_fb;
+}
+
+void Global::SetDomain(RTBox box) {
+ Domain=box;
+}
+
+void Global::GetDomain(RTBox box) {
+ box=Domain;
+}
+
+void Global::DispMinimizers() {
+ copy(SolSet.begin(), SolSet.end(), ostream_iterator<Trial>(cout));
+}
+
+double Global::OneMinimizer(RCRVector x) {
+ if (NoMinimizers()) return 0.0;
+ for (int i=0;i<x.GetLength();i++) x(i) = SolSet.front().xvals(i);
+ return SolSet.front().objval;
+}
+
+bool Global::NoMinimizers() {
+ return SolSet.empty();
+}
+
+void Global::ClearSolSet() {
+ SolSet.erase(SolSet.begin(), SolSet.end()) ;
+}
+
+void Global::AddPoint(RCRVector x, double f) {
+ Trial T(dim);
+ T.xvals=x; T.objval=f;
+ Domain.AddTrial(T);
+ SolSet.push_back(T);
+}
--- /dev/null
+#ifndef GLOBAL_H
+#define GLOBAL_H
+
+#include <queue>
+//#include "function.h"
+#include "tools.h"
+using namespace std;
+
+typedef void dom(RTBox) ;
+typedef dom* Pdom ;
+
+typedef double obj(RCRVector) ;
+typedef obj* Pobj ;
+
+typedef void grad(RCRVector,RVector&) ;
+typedef grad* Pgrad ;
+
+typedef enum { OBJECTIVE_ONLY, GRADIENT_ONLY, OBJECTIVE_AND_GRADIENT } whichO;
+
+typedef double objgrad(RCRVector,RCRVector,whichO) ;
+typedef objgrad* Pobjgrad ;
+
+class GlobalParams {
+public:
+ long int maxtime, maxeval;
+ double eps_cl, mu, rshift;
+ int det_pnts, rnd_pnts;
+};
+
+class Global: public GlobalParams {
+public:
+ // Problem specification
+ int dim ;
+ Pobj Objective ;
+ Pgrad Gradient ;
+ long int numeval;
+
+ virtual double ObjectiveGradient(RCRVector xy, RVector&grad, whichO which){
+ ++numeval;
+ switch (which) {
+ case OBJECTIVE_AND_GRADIENT:
+ Gradient(xy, grad);
+ case OBJECTIVE_ONLY:
+ return Objective(xy);
+ case GRADIENT_ONLY:
+ Gradient(xy, grad);
+ return 0.0;
+ }
+ }
+
+ Global(RTBox, Pobj, Pgrad, GlobalParams);
+// Global& operator=(const Global &);
+
+ void Search(int, RCRVector);
+ void DispMinimizers();
+ double OneMinimizer(RCRVector);
+ bool NoMinimizers();
+ void SetDomain(RTBox);
+ void GetDomain(RTBox);
+ double GetMinValue();
+ void SetMinValue(double);
+ void ClearSolSet();
+ void AddPoint(RCRVector, double);
+
+ long int GetTime();
+ bool InTime();
+
+private:
+ list<Trial> SolSet;
+ list<Trial>::const_iterator titr;
+ priority_queue<TBox> CandSet;
+ priority_queue<TBox> Garbage;
+
+ double fbound;
+ TBox Domain;
+
+ void FillRegular(RTBox, RTBox);
+ void FillRandom(RTBox, RTBox);
+ double NewtonTest(RTBox, int, RCRVector, int*);
+ void ReduceOrSubdivide(RTBox, int, RCRVector);
+};
+#endif
+
--- /dev/null
+/*
+ Temporary implementation of vector and matrix classes
+ No attempt is made to check if the function arguments are valid
+*/
+
+#include <iostream>
+#include <math.h> // for sqrt()
+
+#include "linalg.h"
+
+double eps() {
+ /* Returns the machine precision : (min { x >= 0 : 1 + x > 1 })
+ NB This routine should be replaced by LAPACK_LAMCH */
+ double Current, Last, OnePlusCurrent ;
+ Current = 1.0 ;
+ do
+ {
+ Last = Current ;
+ Current /= 2.0 ;
+ OnePlusCurrent = 1.0 + Current ;
+ } while (OnePlusCurrent > 1.0) ;
+ return Last ;
+}
+
+
+RVector::RVector() {
+ // Constructor
+ len=0;
+ elements=0; (*this)=0.;
+}
+
+RVector::RVector(int n) {
+ // Constructor
+ len=n;
+ elements=new double[len]; (*this)=0.;
+}
+
+RVector::RVector(RCRVector vect)
+{
+ // Constructor + Copy
+ len=vect.len;
+ elements=new double[len]; (*this)=vect;
+}
+
+RCRVector RVector::operator=(RCRVector vect)
+{
+ // Copy
+ for (int i=0;i<len;i++) elements[i]=vect.elements[i];
+ return *this;
+}
+
+RCRVector RVector::operator=(double num) {
+ // Assignment
+ for (int i=0;i<len;i++) elements[i]=num;
+ return *this;
+}
+
+double RVector::nrm2() {
+ double sum=0 ;
+ for (int i = 0; i < len; i++)
+ sum+=elements[i]*elements[i] ;
+ return sqrt(sum) ;
+}
+
+void scal(double alpha, RCRVector x) {
+ int n=x.len ;
+ for (int i = 0; i < n; i++)
+ x.elements[i] = alpha*x.elements[i] ;
+}
+
+double norm2(RCRVector x) {
+ // Euclidian norm
+ double sum=0 ;
+ double* pa=x.elements ;
+ int n=x.len ;
+ for (int i = 0; i < n; i++) {
+ sum+=(*pa)*(*pa) ;
+ pa++ ;
+ }
+ return sqrt(sum) ;
+}
+
+double normInf(RCRVector x) {
+ // Infinity norm
+ int n=x.len ;
+ double tmp=DBL_MIN ;
+ for (int i=0 ; i<n ; i++)
+ tmp=max(tmp,fabs(x.elements[i])) ;
+ return tmp ;
+}
+
+double dot(RCRVector x, RCRVector y) {
+ // dot <- x'y
+ int n=x.len ;
+ double sum=0 ;
+ for (int i=0 ; i<n ; i++)
+ sum+= x.elements[i]*y.elements[i] ;
+ return sum ;
+}
+
+void copy(RCRVector x, RCRVector y) {
+ // y <- x
+ double *px=x.elements ;
+ double *py=y.elements ;
+ int n=x.len ;
+ for (int i=0 ; i<n ; i++)
+ (*py++)=(*px++) ;
+}
+
+void axpy(double alpha, RCRVector x, RCRVector y) {
+ // y <- alpha*x + y
+ double *px=x.elements ;
+ double *py=y.elements ;
+ int n=x.len ;
+ for (int i=0 ; i<n ; i++) {
+ *py=alpha*(*px) + *py ;
+ px++ ; py++ ;
+ }
+}
+
+void gemv(char trans, double alpha, RCRMatrix A,RCRVector x,
+ double beta, RCRVector y) {
+ // Matrix-vector multiplication
+ int i, j, dim=A.Dim;
+ double sum ;
+
+ if (trans=='N') {
+ // y=alpha*A*x + beta*y
+ for (i=0;i<dim;i++) {
+ sum=0.;
+ for (j=0;j<dim;j++)
+ sum+=A.Vals[j+i*dim]*x.elements[j]*alpha;
+ y.elements[i]=y.elements[i]*beta + sum ;
+ }
+ }
+ else {
+ // y=alpha*transpose(A)*x + beta*y
+ for (i=0;i<dim;i++) {
+ sum=0.;
+ for (j=0;j<dim;j++) {
+ sum+=A.Vals[i+j*dim]*x.elements[j]*alpha ;
+ }
+ y.elements[i]=y.elements[i]*beta + sum ;
+ }
+ }
+}
+
+ostream & operator << (ostream & os, const RVector & v) {
+ os << '[';
+ for (int i = 0; i < v.len; i++) {
+ if (i>0) os << "," ;
+ os << v.elements[i] ;
+ }
+ return os << ']';
+}
+
+/*************************** Matrix Class ***********************/
+
+RMatrix::RMatrix() {
+ // Constructor
+ Dim=0 ; Vals=0 ; (*this)=0 ;
+}
+
+RMatrix::RMatrix(int dim) {
+ // Constructor
+ Dim=dim;
+ Vals=new double[long(Dim)*long(Dim)]; (*this)=0.;
+}
+
+RMatrix::RMatrix(RCRMatrix matr) {
+ // Constructor + Copy
+ Dim=matr.Dim;
+ Vals=new double[long(Dim)*long(Dim)]; (*this)=matr;
+}
+
+RCRMatrix RMatrix::operator=(RCRMatrix mat)
+{ // Assignment, A=B
+ long int Len=long(Dim)*long(Dim);
+ for (int i=0;i<Len;i++) Vals[i]=mat.Vals[i] ;
+ return *this;
+}
+
+RCRMatrix RMatrix::operator=(double num) {
+ long int Len=long(Dim)*long(Dim);
+ for (long int i=0;i<Len;i++) Vals[i]=num;
+ return *this;
+}
+
+double& RMatrix::operator()(int vidx,int hidx) {
+ return Vals[vidx*Dim+hidx];
+}
+
+void ger(double alpha, RCRVector x,RCRVector y, RCRMatrix A) {
+ // Rank one update : A=alpha*xy'+A
+ int dim=x.len;
+ double* pa=A.Vals ;
+
+ for (int i=0;i<dim;i++)
+ for (int j=0;j<dim;j++) {
+ *pa=alpha*x.elements[i]*y.elements[j] + *pa;
+ pa++ ;
+ }
+}
+
+ostream & operator << (ostream & os, const RMatrix & A) {
+ int n=A.Dim ;
+ double* pa=A.Vals ;
+ os << endl ;
+ for (int i = 0; i < n; i++) {
+ for (int j=0 ; j< n ; j++) {
+ os << (*(pa++)) << " " ;
+ }
+ os << endl ;
+ }
+ return os ;
+}
--- /dev/null
+/*
+ Temporary implementation of vector and matrix classes
+ This is more or less borrowed from Serguei's program
+*/
+
+#ifndef LINALG_H
+#define LINALG_H
+
+#include <iostream>
+using namespace std;
+#include <math.h> // for sqrt()
+#include <float.h>
+
+typedef const class RVector CRVector;
+typedef CRVector& RCRVector;
+typedef const class RMatrix CRMatrix ;
+typedef CRMatrix& RCRMatrix;
+
+double eps() ;
+
+#define max(A,B) ((A) > (B) ? (A):(B))
+#define min(A,B) ((A) < (B) ? (A):(B))
+
+/********************* Class RVector *********************/
+
+class RVector{
+protected:
+ double* elements; // array of values
+ int len; // size of array
+
+ public:
+ RVector() ;
+ RVector(int); // Constructor
+ RVector(RCRVector); // copy constructor
+ ~RVector() { delete[] elements; elements=0 ; len=0; }
+
+ RCRVector operator=(double) ;
+ RCRVector operator=(RCRVector);
+
+ double & operator () (int i) const {return elements[i] ; }
+ double nrm2() ; // Euclidian norm
+
+ double *raw_data() { return elements; }
+ const double *raw_data_const() const { return elements; }
+
+ friend ostream & operator << (ostream &, const RVector &);
+
+ friend double norm2(RCRVector) ;
+ friend double normInf(RCRVector) ;
+ friend double dot(RCRVector, RCRVector) ;
+ friend void scal(double, RCRVector) ;
+ friend void copy(RCRVector, RCRVector) ;
+ friend void axpy(double, RCRVector, RCRVector) ;
+ friend void gemv(char,double, RCRMatrix, RCRVector, double, RCRVector);
+ friend void ger(double alpha, RCRVector, RCRVector, RCRMatrix);
+
+ int GetLength() const { return len; }; // get vector size
+};
+
+/******************* Class RMatrix *************************/
+
+class RMatrix
+{
+ protected:
+ double* Vals; // array of values
+ int Dim; // dimension
+
+ public:
+ RMatrix() ;
+ RMatrix(int); // dimension
+ ~RMatrix() { delete[] Vals; Vals=0 ; Dim=0; }
+
+ RMatrix(RCRMatrix); // copy constructor
+ RCRMatrix operator=(double num) ;
+ RCRMatrix operator=(RCRMatrix) ; // (needed for template stuff)
+
+ double& operator()(int vidx,int hidx) ;
+ friend ostream & operator << (ostream &, const RMatrix &);
+
+ friend void gemv(char,double, RCRMatrix, RCRVector, double, RCRVector);
+ friend void ger(double alpha,RCRVector,RCRVector,RCRMatrix);
+
+ int GetDim() { return Dim; }; // get dimension
+};
+
+#endif
--- /dev/null
+
+/*
+ Local search - A trust region algorithm with BFGS update.
+*/
+
+#include <iostream>
+#include <stdlib.h>
+
+#include "stogo_config.h"
+#include "global.h"
+#include "local.h"
+#include "tools.h"
+
+int local(Trial &T, TBox &box, TBox &domain, double eps_cl, double *mgr,
+ Global &glob, int axis, RCRVector x_av) {
+
+ int k_max, info, outside ;
+ int k, i, good_enough, iTmp ;
+ int n=box.GetDim();
+
+ double maxgrad, delta, f, f_new, tmp ;
+ double alpha, gamma, beta, d2, s2, nom, den, ro ;
+ double nrm_sd, nrm_hn, snrm_hn, nrm_dl ;
+ RVector x(n), g(n), h_sd(n), h_dl(n), h_n(n), x_new(n), g_new(n) ;
+ RVector s(n),y(n),z(n),w(n) ; // Temporary vectors
+ RMatrix B(n), H(n) ; // Hessian and it's inverse
+
+ k_max = max_iter*n ;
+ x=T.xvals ;
+
+#ifdef LS_DEBUG
+ cout << "Local Search, x=" << x << endl;
+#endif
+
+ if (box.OutsideBox(x, domain) != 0) {
+ cout << "Starting point is not inside the boundary. Exiting...\n" ;
+ exit(1) ;
+ return LS_Out ;
+ }
+
+ // Check if we are close to a stationary point located previously
+ if (box.CloseToMin(x, &tmp, eps_cl)) {
+#ifdef LS_DEBUG
+ cout << "Close to a previously located stationary point, exiting" << endl;
+#endif
+ T.objval=tmp;
+ return LS_Old ;
+ }
+
+ // Initially B and H are equal to the identity matrix
+ B=0 ; H=0 ;
+ for (i=0 ; i<n ; i++) {
+ B(i,i)=1 ;
+ H(i,i)=1 ;
+ }
+
+ RVector g_av(x_av.GetLength());
+ if (axis==-1) {
+ f=glob.ObjectiveGradient(x,g,OBJECTIVE_AND_GRADIENT);
+ }
+ else {
+ x_av(axis)=x(0);
+ f=glob.ObjectiveGradient(x_av,g_av,OBJECTIVE_AND_GRADIENT);
+ g(0)=g_av(axis);
+ }
+ FC++;GC++;
+
+ if (axis == -1) {
+ // Skipping AV
+#ifdef INI3
+ // Elaborate scheme to initalize delta
+ delta=delta_coef*norm2(g) ;
+ copy(g,z) ;
+ axpy(1.0,x,z) ;
+ if (!box.InsideBox(z)) {
+ if (box.Intersection(x,g,z)==TRUE) {
+ axpy(-1.0,x,z) ;
+ delta=min(delta,delta_coef*norm2(z)) ;
+ }
+ else {
+ // Algorithm broke down, use INI1
+ delta = (1.0/7)*box.ShortestSide(&iTmp) ;
+ }
+ }
+#endif
+#ifdef INI2
+ // Use INI2 scheme
+ delta = box.ClosestSide(x)*delta_coef ;
+ if (delta<MacEpsilon)
+ // Patch to avoid trust region with radius close to zero
+ delta = (1.0/7)*box.ShortestSide(&iTmp) ;
+#endif
+#ifdef INI1
+ delta = delta_coef*box.ShortestSide(&iTmp) ;
+#endif
+ }
+ else {
+ // Use a simple scheme for the 1D minimization (INI1)
+ delta = (1.0/7.0)*box.ShortestSide(&iTmp) ;
+ }
+
+ k=0 ; good_enough = 0 ; info=LS_New ; outside=0 ;
+ maxgrad=*mgr ;
+ while (good_enough == 0) {
+ k++ ;
+ if (k>k_max) {
+#ifdef LS_DEBUG
+ cout << "Maximum number of iterations reached\n" ;
+#endif
+ info=LS_MaxIter ;
+ break ;
+ }
+
+ // Update maximal gradient value
+ maxgrad=max(maxgrad,normInf(g)) ;
+
+ // Steepest descent, h_sd = -g
+ copy(g,h_sd) ;
+ scal(-1.0,h_sd) ;
+ nrm_sd=norm2(h_sd) ;
+
+ if (nrm_sd < epsilon) {
+ // Stop criterion (gradient) fullfilled
+#ifdef LS_DEBUG
+ cout << "Gradient small enough" << endl ;
+#endif
+ good_enough = 1 ;
+ break ;
+ }
+
+ // Compute Newton step, h_n = -H*g
+ gemv('N',-1.0, H, g, 0.0, h_n) ;
+ nrm_hn = norm2(h_n) ;
+
+ if (nrm_hn < delta) {
+ // Pure Newton step
+ copy(h_n, h_dl) ;
+#ifdef LS_DEBUG
+ cout << "[Newton step] " ;
+#endif
+ }
+ else {
+ gemv('N',1.0,B,g,0.0,z) ;
+ tmp=dot(g,z) ;
+ if (tmp==0) {
+ info = LS_Unstable ;
+ break ;
+ }
+ alpha=(nrm_sd*nrm_sd)/tmp ; // Normalization (N38,eq. 3.30)
+ scal(alpha,h_sd) ;
+ nrm_sd=fabs(alpha)*nrm_sd ;
+
+ if (nrm_sd >= delta) {
+ gamma = delta/nrm_sd ; // Normalization (N38, eq. 3.33)
+ copy(h_sd,h_dl) ;
+ scal(gamma,h_dl) ;
+#ifdef LS_DEBUG
+ cout << "[Steepest descent] " ;
+#endif
+ }
+ else {
+ // Combination of Newton and SD steps
+ d2 = delta*delta ;
+ copy(h_sd,s) ;
+ s2=nrm_sd*nrm_sd ;
+ nom = d2 - s2 ;
+ snrm_hn=nrm_hn*nrm_hn ;
+ tmp = dot(h_n,s) ;
+ den = tmp-s2 + sqrt((tmp-d2)*(tmp-d2)+(snrm_hn-d2)*(d2-s2)) ;
+ if (den==0) {
+ info = LS_Unstable ;
+ break ;
+ }
+ // Normalization (N38, eq. 3.31)
+ beta = nom/den ;
+ copy(h_n,h_dl) ;
+ scal(beta,h_dl) ;
+ axpy((1-beta),h_sd,h_dl) ;
+#ifdef LS_DEBUG
+ cout << "[Mixed step] " ;
+#endif
+ }
+ }
+ nrm_dl=norm2(h_dl) ;
+
+ //x_new = x+h_dl ;
+ copy(x,x_new) ;
+ axpy(1.0,h_dl,x_new) ;
+
+ // Check if x_new is inside the box
+ iTmp=box.OutsideBox(x_new, domain) ;
+ if (iTmp == 1) {
+#ifdef LS_DEBUG
+ cout << "x_new is outside the box " << endl ;
+#endif
+ outside++ ;
+ if (outside>max_outside_steps) {
+ // Previous point was also outside, exit
+ break ;
+ }
+ }
+ else if (iTmp == 2) {
+#ifdef LS_DEBUG
+ cout << " x_new is outside the domain" << endl ;
+#endif
+ info=LS_Out ;
+ break ;
+ }
+ else {
+ outside=0 ;
+ }
+
+ // Compute the gain
+ if (axis==-1)
+ f_new=glob.ObjectiveGradient(x_new,x_new,OBJECTIVE_ONLY);
+ else {
+ x_av(axis)=x_new(0);
+ f_new=glob.ObjectiveGradient(x_av,x_new,OBJECTIVE_ONLY);
+ }
+ FC++;
+ gemv('N',0.5,B,h_dl,0.0,z);
+ ro = (f_new-f) / (dot(g,h_dl) + dot(h_dl,z)); // Quadratic model
+ if (ro > 0.75) {
+ delta = delta*2;
+ }
+ if (ro < 0.25) {
+ delta = delta/3;
+ }
+ if (ro > 0) {
+ // Update the Hessian and it's inverse using the BFGS formula
+ if (axis==-1)
+ glob.ObjectiveGradient(x_new,g_new,GRADIENT_ONLY);
+ else {
+ x_av(axis)=x_new(0);
+ glob.ObjectiveGradient(x_av,g_av,GRADIENT_ONLY);
+ g_new(0)=g_av(axis);
+ }
+ GC++;
+
+ // y=g_new-g
+ copy(g_new,y);
+ axpy(-1.0,g,y);
+
+ // Check curvature condition
+ alpha=dot(y,h_dl);
+ if (alpha <= sqrt(MacEpsilon)*nrm_dl*norm2(y)) {
+#ifdef LS_DEBUG
+ cout << "Curvature condition violated " ;
+#endif
+ }
+ else {
+ // Update Hessian
+ gemv('N',1.0,B,h_dl,0.0,z) ; // z=Bh_dl
+ beta=-1/dot(h_dl,z) ;
+ ger(1/alpha,y,y,B) ;
+ ger(beta,z,z,B) ;
+
+ // Update Hessian inverse
+ gemv('N',1.0,H,y,0.0,z) ; // z=H*y
+ gemv('T',1.0,H,y,0.0,w) ; // w=y'*H
+ beta=dot(y,z) ;
+ beta=(1+beta/alpha)/alpha ;
+
+ // It should be possible to do this updating more efficiently, by
+ // exploiting the fact that (h_dl*y'*H) = transpose(H*y*h_dl')
+ ger(beta,h_dl,h_dl,H) ;
+ ger(-1/alpha,z,h_dl,H) ;
+ ger(-1/alpha,h_dl,w,H) ;
+ }
+
+ if (nrm_dl < norm2(x)*epsilon) {
+ // Stop criterion (iteration progress) fullfilled
+#ifdef LS_DEBUG
+ cout << "Progress is marginal" ;
+#endif
+ good_enough = 1 ;
+ }
+
+ // Check if we are close to a stationary point located previously
+ if (box.CloseToMin(x_new, &f_new, eps_cl)) {
+ // Note that x_new and f_new may be overwritten on exit from CloseToMin
+#ifdef LS_DEBUG
+ cout << "Close to a previously located stationary point, exiting" << endl;
+#endif
+ info = LS_Old ;
+ good_enough = 1 ;
+ }
+
+ // Update x, g and f
+ copy(x_new,x) ; copy(g_new,g) ; f=f_new ;
+
+#ifdef LS_DEBUG
+ cout << " x=" << x << endl ;
+#endif
+
+ }
+ else {
+#ifdef LS_DEBUG
+ cout << "Step is no good, ro=" << ro << " delta=" << delta << endl ;
+#endif
+ }
+
+ } // wend
+
+ // Make sure the routine returns correctly...
+ // Check if last iterate is outside the boundary
+ if (box.OutsideBox(x, domain) != 0) {
+ info=LS_Out; f=DBL_MAX;
+ }
+
+ if (info == LS_Unstable) {
+ cout << "Local search became unstable. No big deal but exiting anyway\n" ;
+ exit(1);
+ }
+
+ T.xvals=x ; T.objval=f ;
+ *mgr=maxgrad ;
+ if (outside>0)
+ return LS_Out ;
+ else
+ return info ;
+}
--- /dev/null
+/*
+ Definitions of various variables used in the local search
+*/
+
+#ifndef LOCAL_H
+#define LOCAL_H
+
+#include "tools.h"
+#include "global.h"
+
+extern int FC, GC ;
+
+// Results of local search
+enum {LS_Unstable, LS_MaxIter, LS_Old, LS_New,LS_Out} ;
+
+const double delta_coef = 1.0/2.0; // Initialization of trust region
+const double epsilon = 1.0E-4 ; // Stopping criterion, var 1e-4
+const int max_outside_steps=1 ; // Maximum number of steps outside the box
+const int max_iter=50 ; // Max iterations = max_iter*dim. of problem
+
+extern double MacEpsilon ; // min {x >= 0 : 1 + x > 1}
+
+int local(Trial &, TBox &, TBox &, double, double*, Global&, int, RCRVector) ;
+
+#endif
--- /dev/null
+/*
+ A simple program to test the global optimizer.
+*/
+
+#include "global.h"
+#include "tools.h"
+#include "testfun.h"
+
+#define STRLEN_MAX 80
+
+int main() {
+ bool AVfail, AVflag;
+ int testfnc, dim, axis, i;
+ double AVbest;
+ Pdom Dom;
+ Pgrad Grad;
+ Pobj Obj;
+
+ cout << "Global Optimizer v0.1" << endl;
+ cout << "1. Rosenbrock (min=0)" << endl;
+ cout << "2. McCormick (min=-1.91)" << endl;
+ cout << "3. Box & Betts (min=0)" << endl;
+ cout << "4. Paviani (min=-45.7)" << endl;
+ cout << "5. Extended Beale (min=0)" << endl;
+ cout << "6. Three-hump Camel (min=0)" << endl;
+ cout << "7. Jacobsen & Pedersen (min=3)" << endl;
+ cout << "8. Poly1 (min=0)" << endl;
+ cout << "9. Six-hump Camel (min=-1.03) " << endl;
+ cout << "10. One-dimensional (min=-23.58)" << endl;
+ cout << "11. Kowalik (min=0.0003075)" << endl;
+ cout << "12. Hansen (min=-176.54)" << endl;
+ cout << "13. Shubert (min=-24.06)" << endl;
+ cout << "14. Griewank (min=0)" << endl;
+ cout << "15. Levy4 (min=-21.502)" << endl;
+ cout << "16. Levy5 (min=-11.504)" << endl;
+ cout << "17. Levy6 (min=-11.504)" << endl;
+ cout << "18. Levy7 (min=-11.504)" << endl;
+ cout << "19. Rastrigin (min=0)" << endl;
+ cout << "20. Trid" << endl;
+ cout << "21. Perm(4,50)" << endl;
+ cout << "22. Perm(4,0.5)" << endl;
+ cout << "23. Powersum(8,18,44,114)" << endl;
+ cout << "24. Schwefel (min=-12569.5)" << endl;
+ cout << "25. Shekel (-10.0)" << endl;
+ cout << endl << "Select test function :";
+ cin >> testfnc;
+
+ switch (testfnc) {
+ case 10:
+ dim=1;
+ Dom=Domain_OneDim;
+ Obj=Objective_OneDim;
+ Grad=Gradient_OneDim;
+ break;
+ case 1:
+ dim=2;
+ Dom=Domain_Rosenbrock;
+ Obj=Objective_Rosenbrock;
+ Grad=Gradient_Rosenbrock;
+ break;
+ case 8:
+ dim=2;
+ Dom=Domain_Poly1;
+ Obj=Objective_Poly1;
+ Grad=Gradient_Poly1;
+ break;
+ case 7:
+ dim=2;
+ Dom=Domain_Jacobsen;
+ Obj=Objective_Jacobsen;
+ Grad=Gradient_Jacobsen;
+ break;
+ case 6:
+ dim=2;
+ Dom=Domain_Camel3;
+ Obj=Objective_Camel3;
+ Grad=Gradient_Camel3;
+ break;
+ case 5:
+ dim=2;
+ Dom=Domain_Beale3;
+ Obj=Objective_Beale3;
+ Grad=Gradient_Beale3;
+ break;
+ case 4:
+ dim=10;
+ Dom=Domain_Paviani;
+ Obj=Objective_Paviani;
+ Grad=Gradient_Paviani;
+ break;
+ case 12:
+ dim=2;
+ Dom=Domain_Hansen;
+ Obj=Objective_Hansen;
+ Grad=Gradient_Hansen;
+ break;
+ case 13:
+ dim=2;
+ Dom=Domain_Shubert;
+ Obj=Objective_Shubert;
+ Grad=Gradient_Shubert;
+ break;
+ case 11:
+ dim=4;
+ Dom=Domain_Kowalik;
+ Obj=Objective_Kowalik;
+ Grad=Gradient_Kowalik;
+ break;
+ case 9:
+ dim=2;
+ Dom=Domain_Camel6;
+ Obj=Objective_Camel6;
+ Grad=Gradient_Camel6;
+ break;
+ case 2:
+ dim=2;
+ Dom=Domain_McCormick;
+ Obj=Objective_McCormick;
+ Grad=Gradient_McCormick;
+ break;
+ case 3:
+ dim=3;
+ Dom=Domain_BoxBetts;
+ Obj=Objective_BoxBetts;
+ Grad=Gradient_BoxBetts;
+ break;
+ case 14:
+ dim=10;
+ Dom=Domain_Griewank;
+ Obj=Objective_Griewank;
+ Grad=Gradient_Griewank;
+ break;
+ case 15:
+ dim=4;
+ Dom=Domain_Levy;
+ Obj=Objective_Levy;
+ Grad=Gradient_Levy;
+ break;
+ case 16:
+ dim=5;
+ Dom=Domain_Levy;
+ Obj=Objective_Levy;
+ Grad=Gradient_Levy;
+ break;
+ case 17:
+ dim=6;
+ Dom=Domain_Levy;
+ Obj=Objective_Levy;
+ Grad=Gradient_Levy;
+ break;
+ case 18:
+ dim=7;
+ Dom=Domain_Levy;
+ Obj=Objective_Levy;
+ Grad=Gradient_Levy;
+ break;
+ case 19:
+ cout << "Enter problem dimension ";
+ int rast_dim;
+ cin >> rast_dim;
+ dim=rast_dim;
+ Dom=Domain_Rastrigin;
+ Obj=Objective_Rastrigin;
+ Grad=Gradient_Rastrigin;
+ break;
+ case 20:
+ cout << "Enter problem dimension (two or larger) ";
+ int trid_dim;
+ cin >> trid_dim;
+ dim=trid_dim;
+ Dom=Domain_Trid;
+ Obj=Objective_Trid;
+ Grad=Gradient_Trid;
+ break;
+ case 21:
+ dim=4;
+ Dom=Domain_Perm;
+ Obj=Objective_Perm_4_50;
+ Grad=Gradient_Perm_4_50;
+ break;
+ case 22:
+ dim=4;
+ Dom=Domain_Perm;
+ Obj=Objective_Perm_4_05;
+ Grad=Gradient_Perm_4_05;
+ break;
+ case 23:
+ dim=4;
+ Dom=Domain_Powersum;
+ Obj=Objective_Powersum;
+ Grad=Gradient_Powersum;
+ break;
+ case 24:
+ cout << "Enter problem dimension ";
+ int schwef_dim;
+ cin >> schwef_dim;
+ dim=schwef_dim;
+ Dom=Domain_Schwefel;
+ Obj=Objective_Schwefel;
+ Grad=Gradient_Schwefel;
+ break;
+ case 25:
+ dim=4;
+ Dom=Domain_Shekel;
+ Obj=Objective_Shekel;
+ Grad=Gradient_Shekel;
+ break;
+ default:
+ cout << "Error : Function not defined" << endl;
+ exit(1);
+ }
+
+ cout << "Dimension=" <<dim << endl;
+ TBox D(dim);
+ Dom(D);
+ cout << "Domain=";
+ for (i=0; i<dim; i++)
+ cout << "[" << D.lb(i) << "," << D.ub(i) << "]";
+ cout << endl << endl;
+
+ GlobalParams params;
+ cout << "Enter time limit (seconds) ";
+ cin >> params.maxtime;
+ if (params.maxtime<1) {
+ cout << "Warning: time limit set to 1 second\n";
+ }
+ params.maxeval = 0;
+ cout << "Use factory settings (y/n) ";
+ char str[STRLEN_MAX]; cin >> str;
+ if (str[0]=='y') {
+ params.det_pnts=2*dim+1; params.rnd_pnts=0;
+ params.eps_cl=0.1; params.rshift=0.3;
+ params.mu=1.0E-4; AVflag=FALSE;
+ }
+ else {
+ cout << "Number of deterministic points ";
+ cin >> params.det_pnts;
+ cout << "Numer of stochastic points ";
+ cin >> params.rnd_pnts;
+ cout << "Radius of attraction ";
+ cin >> params.eps_cl;
+ cout << "Parameter rshift ";
+ cin >> params.rshift;
+ cout << "Parameter mu ";
+ cin >> params.mu;
+ cout << "Use the AV initialization (y/n) ";
+ cin >> str;
+ if (str[0]=='y') AVflag=TRUE; else AVflag=FALSE;
+ }
+
+ Global Problem(D,Obj, Grad, params);
+ RVector x_av(dim);
+ if (AVflag==TRUE) {
+ cout << "Enter time limit for each coordinate direction (seconds) ";
+ cin >> params.maxtime;
+ if (params.maxtime<1) {
+ cout << "Warning: time limit set to 1 second\n";
+ }
+ params.maxeval = 0;
+ params.det_pnts=3;
+ TBox I(1);
+ Global AV(I, Obj, Grad, params);
+
+ x_av=0.0; AVfail=FALSE;
+ for (axis=0; axis<Problem.dim; axis++) {
+ cout << "### axis=" << axis << " ###" << endl;
+ I.lb=(D.lb)(axis); I.ub=(D.ub)(axis);
+ AV.SetDomain(I);
+ AV.ClearSolSet();
+ AV.Search(axis, x_av);
+
+ if (AV.NoMinimizers()) {
+ cout << "AV failed with axis=" << axis << endl;
+ AVfail=TRUE; break;
+ }
+ }
+
+ if (AVfail==FALSE) {
+ AVbest=AV.GetMinValue();
+ cout << "### AV Located x=" << x_av << " fbound=" << AVbest << endl;
+ RVector AVx(Problem.dim);
+ AVx=x_av;
+
+ // Use result from AV for new fbound
+ Problem.SetMinValue(AVbest);
+
+ // Add the best point found to the initial box (domain)
+ Problem.AddPoint(x_av, AVbest);
+ }
+ }
+
+ // Perform the main search
+ cout << "### Starting main search ###" << endl;
+ Problem.Search(-1, x_av);
+
+ cout << "Optimization terminated. Current set of minimizers is" << endl;
+ if (Problem.NoMinimizers() && AVflag==FALSE)
+ cout << "### No improvement found ###" << endl;
+ else
+ Problem.DispMinimizers();
+}
--- /dev/null
+#include "linalg.h"
+#include "tools.h"
+#include "stogo_config.h"
+
+void Domain_Rosenbrock(RTBox box) {
+ box.lb=-10.0 ; box.ub=10.0;
+}
+
+double Objective_Rosenbrock(RCRVector x) {
+ double a=x(1)-x(0)*x(0);
+ double b=1-x(0);
+ return 100*a*a + b*b;
+}
+
+void Gradient_Rosenbrock(RCRVector x, RCRVector grad) {
+ grad(0)=200*(x(1)-x(0)*x(0))*(-2*x(0))-2*(1-x(0));
+ grad(1)=200*(x(1)-x(0)*x(0));
+}
--- /dev/null
+// A C-callable front-end to the StoGO global-optimization library.
+// -- Steven G. Johnson
+
+#include "stogo.h"
+#include "global.h"
+
+class MyGlobal : public Global {
+protected:
+ objective_func my_func;
+ void *my_data;
+
+public:
+
+ MyGlobal(RTBox D, GlobalParams P, objective_func func, void *data) : Global(D, 0, 0, P), my_func(func), my_data(data) {}
+
+ virtual double ObjectiveGradient(RCRVector xy, RVector &grad, whichO which){
+ ++numeval;
+ switch (which) {
+ case GRADIENT_ONLY:
+ case OBJECTIVE_AND_GRADIENT:
+ return my_func(xy.GetLength(), xy.raw_data_const(), grad.raw_data(), my_data);
+ case OBJECTIVE_ONLY:
+ return my_func(xy.GetLength(), xy.raw_data_const(), NULL, my_data);
+ }
+ }
+};
+
+int stogo_minimize(int n,
+ objective_func fgrad, void *data,
+ double *x, double *fmin,
+ const double *l, const double *u,
+ long int maxeval, long int maxtime)
+{
+ GlobalParams params;
+
+ // FIXME: WTF do these parameters mean?
+ params.det_pnts=2*n+1; params.rnd_pnts=0;
+ params.eps_cl=0.1; params.rshift=0.3;
+ params.mu=1.0E-4;
+
+ params.maxtime = maxtime;
+ params.maxeval = maxeval;
+
+ TBox D(n);
+ for (int i = 0; i < n; ++i) {
+ D.lb(i) = l[i];
+ D.ub(i) = u[i];
+ }
+
+ MyGlobal Problem(D, params, fgrad, data);
+ RVector dummyvec(n);
+ Problem.Search(-1, dummyvec);
+
+ if (Problem.NoMinimizers())
+ return 0;
+
+ *fmin = Problem.OneMinimizer(dummyvec);
+ for (int i = 0; i < n; ++i) x[i] = dummyvec(i);
+ return 1;
+}
--- /dev/null
+/* A C-callable front-end to the StoGO global-optimization library.
+ -- Steven G. Johnson */
+
+#ifndef STOGO_H
+#define STOGO_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef double (*objective_func)(int n, const double *x, double *grad,
+ void *data);
+
+/* search for the global minimum of the function fgrad(n, x, grad, data)
+ inside a simple n-dimensional hyperrectangle.
+
+ Input:
+
+ n: dimension of search space (number of decision variables)
+
+ fgrad: the objective function of the form fgrad(n, x, grad, data),
+ returning the objective function at x, where
+ n: dimension of search space
+ x: pointer to array of length n, point to evaluate
+ grad: if non-NULL, an array of length n which
+ should on return be the gradient d(fgrad)/dx
+ [ if NULL, grad should be ignored ]
+ data: arbitrary data pointer, whose value is the
+ data argument of stogo_minimize
+
+ data: arbitrary pointer to any auxiliary data needed by fgrad
+
+ l, u: arrays of length n giving the lower and upper bounds of the
+ search space
+
+ maxeval: if nonzero, a maximum number of fgrad evaluations
+
+ maxtime: if nonzero, a maximum time (in seconds)
+
+ Output:
+
+ fmin: the minimum value of the objective function found
+
+ x: pointer to array of length n, giving the location of the minimum
+
+ Return value: 0 if no minimum found, 1 otherwise.
+
+ */
+
+int stogo_minimize(int n,
+ objective_func fgrad, void *data,
+ double *x, double *fmin,
+ const double *l, const double *u,
+ long int maxeval, long int maxtime);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
--- /dev/null
+#ifndef STOGO_CONFIG_H
+#define STOGO_CONFIG_H
+
+// Compiler flags to enable/disable various options in the code
+
+// To obtain extra information on the global search process:
+//#define GS_DEBUG
+
+// To obtain info on the local search process:
+//#define LS_DEBUG
+
+// Initialization strategy in the local search (INI1,INI2 or INI3)
+#define INI2
+
+// Use the more pessimistic strategy for lower bound estimation
+//#define LB2
+
+
+#endif
--- /dev/null
+#ifndef TESTFUN_H
+#define TESTFUN_H
+
+#include "linalg.h"
+#include "tools.h"
+#include "stogo_config.h"
+
+const double pi=fabs(acos(-1.));
+
+/* The Matrix a and vector c are needed in the Shekel function */
+static double a[10][4]={ { 4,4,4,4 } ,
+ { 1,1,1,1 } ,
+ { 8,8,8,8 } ,
+ { 6,6,6,6 } ,
+ { 3,7,3,7 } ,
+ { 2,9,2,9 } ,
+ { 5,5,3,3 } ,
+ { 8,1,8,1 } ,
+ { 6,2,6,2 } ,
+ {7,3.6,7,3.6} };
+static double c[10]= { .1 , .2 , .2 , .4 , .4 , .6 , .3, .7 , .5 , .5 };
+
+void Domain_Shekel(RTBox box) {
+ box.lb=0.0 ; box.ub=10.0 ;
+}
+
+double Objective_Shekel(RCRVector x) {
+ int n=x.GetLength() ;
+ double R=0.0, S;
+ for(int i=0;i<10;i++) {
+ S=0;
+ for(int j=0;j<n;j++) S+=pow(x(j)-a[i][j],2);
+ R-=1/(S+c[i]);
+ }
+ return R;
+}
+
+void Gradient_Shekel(RCRVector x, RVector &grad) {
+ int n=x.GetLength() ;
+ double R;
+ for(int k=0;k<n;k++) {
+ R=0.0;
+ for(int i=0;i<10;i++) {
+ R+=(2.0*x(k)-2.0*a[i][k])/(pow(pow(x(0)-a[i][0],2.0)+pow(x(1)-a[i][1],2.0)
+ +pow(x(2)-a[i][2],2.0)+pow(x(3)-a[i][3],2.0)+c[i],2.0));
+ }
+ grad(k)=R;
+ }
+}
+
+
+/******************** Unimodal functions ******************/
+
+void Domain_Rosenbrock(RTBox box) {
+ box.lb=-10.0 ; box.ub=10.0 ;
+}
+
+double Objective_Rosenbrock(RCRVector x) {
+ double a=x(1)-x(0)*x(0) ;
+ double b=1-x(0) ;
+ return 100*a*a + b*b ;
+}
+
+void Gradient_Rosenbrock(RCRVector x, RVector &grad) {
+ grad(0)=200*(x(1)-x(0)*x(0))*(-2*x(0))-2*(1-x(0)) ;
+ grad(1)=200*(x(1)-x(0)*x(0)) ;
+}
+
+
+void Domain_McCormick(RTBox box) {
+ box.lb(0)=-1.5 ; box.ub(0)=4.0 ;
+ box.lb(1)=-3.0 ; box.ub(1)=4.0 ;
+}
+
+double Objective_McCormick(RCRVector x) {
+ return sin(x(0)+x(1)) + pow(x(0)-x(1),2.0) - 1.5*x(0) + 2.5*x(1) + 1.0 ;
+}
+
+void Gradient_McCormick(RCRVector x, RVector &grad) {
+ grad(0)=cos(x(0)+x(1)) + 2*(x(0)-x(1)) - 1.5 ;
+ grad(1)=cos(x(0)+x(1)) - 2*(x(0)-x(1)) + 2.5 ;
+}
+
+
+void Domain_BoxBetts(RTBox box) {
+ box.lb(0)=0.9 ; box.ub(0)=1.2 ;
+ box.lb(1)=9.0 ; box.ub(1)=11.2 ;
+ box.lb(2)=0.9 ; box.ub(2)=1.2 ;
+}
+
+double Objective_BoxBetts(RCRVector x) {
+ double x0=x(0),x1=x(1),x2=x(2) ;
+ double sum=0.0 ;
+ for (int i=1 ; i<=10 ; i++)
+ sum+=pow(exp(-0.1*i*x0)-exp(-0.1*i*x1)-(exp(-0.1*i)-exp(-1.0*i))*x2,2.0);
+ return sum ;
+}
+
+void Gradient_BoxBetts(RCRVector x, RVector &grad) {
+ double x0=x(0),x1=x(1),x2=x(2) ;
+ double g0=0.0, g1=0.0, g2=0.0 ;
+ for (int i=1 ; i<=10 ; i++) {
+ g0 += -0.2*(exp(-0.1*i*x0)-exp(-0.1*i*x1)
+ -(exp(-0.1*i)-exp(-1.0*i))*x2)*i*exp(-0.1*i*x0);
+ g1 += 0.2*(exp(-0.1*i*x0)-exp(-0.1*i*x1)-(exp(-0.1*i)
+ -exp(-1.0*i))*x2)*i*exp(-0.1*i*x1);
+ g2 += 2.0*(exp(-0.1*i*x0)-exp(-0.1*i*x1)
+ -(exp(-0.1*i)-exp(-1.0*i))*x2)*(-exp(-0.1*i)+exp(-1.0*i));
+ }
+ grad(0)=g0 ; grad(1)=g1 ; grad(2)=g2 ;
+}
+
+
+void Domain_Paviani(RTBox box) {
+ box.lb=2.001 ; box.ub=9.999 ;
+}
+
+double Objective_Paviani(RCRVector x) {
+ double a,b,sum=0.0, mul=1.0 ;
+ int n=x.GetLength() ;
+ for (int i=0 ; i<n ; i++) {
+ a=log(x(i)-2.0) ; b=log(10.0-x(i)) ;
+ sum+= a*a + b*b ;
+ mul*= x(i) ;
+ }
+ return sum - pow(mul,0.2) ;
+}
+
+void Gradient_Paviani(RCRVector x, RVector &grad) {
+ double sum, mul=1.0 ;
+ int n=10 ;
+ for (int i=0 ; i<n ; i++) {
+ mul*= x(i) ;
+ }
+
+ for (int j=0 ; j<n ; j++) {
+ sum=2*log(x(j)-2.0)/(x(j)-2.0) - 2*log(10.0-x(j))/(10.0-x(j)) ;
+ grad(j) = sum - 0.2*(mul/x(j))/pow(mul,0.8) ;
+ }
+}
+
+
+void Domain_Beale3(RTBox box) { // NB Make sure that there is only one minima
+ box.lb=-0.5 ; box.ub=4.0 ;
+}
+
+double Objective_Beale3(RCRVector x) {
+ double x1=x(0), x2=x(1) ;
+ double b1=1.500 + (x2 - 1.)*x1 ;
+ double b2=2.250 + (x2*x2 - 1.)*x1 ;
+ double b3=2.625 + (x2*x2*x2 - 1.)*x1 ;
+ return b1*b1 + b2*b2 + b3*b3 ;
+}
+
+void Gradient_Beale3(RCRVector x, RVector &grad) {
+ double x1=x(0), x2=x(1) ;
+ grad(0) = 2*(1.500 + (x2 - 1.)*x1)*(x2 - 1.)
+ + 2*(2.250 + (x2*x2 - 1.)*x1)*(x2*x2 - 1.)
+ + 2*(2.625 + (x2*x2*x2 - 1.)*x1)*(x2*x2*x2 - 1.) ;
+
+ grad(1) = 2*(1.500 + (x2 - 1.)*x1)*x1
+ + 4*(2.250 + (x2*x2 - 1.)*x1)*x2*x1
+ + 6*(2.625 + (x2*x2*x2 - 1.)*x1)*x2*x2*x1 ;
+}
+
+/************** Difficult unimodal problems ****************/
+void Domain_Trid(RTBox box) {
+ int n=(box.lb).GetLength() ;
+ double tmp=pow(n,2.0);
+ box.lb=-tmp ; box.ub=tmp ; // [-n^2,n^2]
+}
+
+double Objective_Trid(RCRVector x) {
+ int n=x.GetLength();
+ double sum1=0.0, sum2=0.0;
+ for (int i=1 ; i<=n ; i++)
+ sum1+=pow(x(i-1)-1,2.0);
+ for (int i=2 ; i<=n ; i++)
+ sum2+=x(i-1)*x(i-2);
+ return sum1 - sum2;
+}
+
+void Gradient_Trid(RCRVector x, RVector &grad) {
+ int n=x.GetLength();
+ grad(0)=2*(x(0)-1)-x(1) ;
+ for (int i=1 ; i<=n-2 ; i++)
+ grad(i)=2*(x(i)-1) - (x(i-1) + x(i+1)) ;
+ grad(n-1)=2*(x(n-1)-1) - x(n-2) ;
+}
+
+/******************** Multimodal functions ****************/
+
+void Domain_OneDim(RTBox box) {
+ box.lb=-25.0 ; box.ub=25.0 ;
+}
+
+double Objective_OneDim(RCRVector x) {
+ return x(0)*sin(x(0)) ;
+}
+
+void Gradient_OneDim(RCRVector x, RVector &grad) {
+ grad(0) = x(0)*cos(x(0)) + sin(x(0)) ;
+}
+
+void Domain_Poly1(RTBox box) {
+// box.lb=-5.0 ; box.ub=5.0 ;
+ box.lb=-4.0;box.ub=4.01;
+}
+
+double Objective_Poly1(RCRVector x) {
+ double a=(x(0)*x(0) + x(1)*x(1) - 11) ;
+ double b=(x(0)+x(1)*x(1) - 7) ;
+ return a*a + b*b ;
+}
+
+void Gradient_Poly1(RCRVector x, RVector &grad) {
+ double a=(x(0)*x(0) + x(1)*x(1) - 11) ;
+ double b=(x(0)+x(1)*x(1) - 7) ;
+ grad(0) = 4*x(0)* a + 2*b ;
+ grad(1) = 4*x(1)* a + 4*x(1)*b ;
+}
+
+void Domain_Jacobsen(RTBox box) {
+ box.lb=-100.0 ; box.ub=100.0 ;
+}
+
+double Objective_Jacobsen(RCRVector x) {
+ double a=(x(0)*x(0) + x(1) - 11) ;
+ double b=(x(0)+x(1)*x(1) - 7) ;
+ return a*a + b*b + 3;
+}
+
+void Gradient_Jacobsen(RCRVector x, RVector &grad) {
+ double a=(x(0)*x(0) + x(1) - 11) ;
+ double b=(x(0)+x(1)*x(1) - 7) ;
+ grad(0) = 4*x(0)* a + 2*b ;
+ grad(1) = 2*x(1)* a + 4*x(1)*b ;
+}
+
+void Domain_Camel3(RTBox box) {
+ box.lb=-1.0 ; box.ub=2.0 ;
+}
+
+double Objective_Camel3(RCRVector x) {
+ double a=x(0)*x(0) ;
+ return 12*a - 6.3*a*a + a*a*a + 6*x(1)*(x(1)-x(0)) ;
+}
+
+void Gradient_Camel3(RCRVector x, RVector &grad) {
+ double a=x(0)*x(0) ;
+ grad(0) = 24*x(0) - 25.2*a*x(0) + 6*a*a*x(0) - 6*x(1) ;
+ grad(1) = 12*x(1) - 6*x(0) ;
+}
+
+
+void Domain_Hansen(RTBox box) {
+ box.lb=-10.0 ; box.ub=10.0 ;
+}
+
+double Objective_Hansen(RCRVector x) {
+ return (cos(1.0)+2.0*cos(x(0)+2.0)+3.0*cos(2.0*x(0)+3.0)+4.0*cos(3.0*x(0)
+ +4.0)+5.0*cos(4.0*x(0)+5.0))*(cos(2.0*x(1)+1.0)+2.0*cos(3.0*x(1)+2.0)
+ +3.0*cos(4.0*x(1)+3.0)+4.0*cos(5.0*x(1)+4.0)+5.0*cos(6.0*x(1)+5.0));
+}
+
+void Gradient_Hansen(RCRVector x, RVector &grad) {
+ grad(0) = (-2.0*sin(x(0)+2.0)-6.0*sin(2.0*x(0)+3.0)-12.0*sin(3.0*x(0)+4.0)
+ -20.0*sin(4.0*x(0)+5.0))*(cos(2.0*x(1)+1.0)+2.0*cos(3.0*x(1)+2.0)
+ +3.0*cos(4.0*x(1)+3.0)+4.0*cos(5.0*x(1)+4.0)+5.0*cos(6.0*x(1)+5.0));
+
+ grad(1) = (cos(1.0)+2.0*cos(x(0)+2.0)+3.0*cos(2.0*x(0)+3.0)+4.0*cos(3.0*x(0)
+ +4.0)+5.0*cos(4.0*x(0)+5.0))*(-2.0*sin(2.0*x(1)+1.0)
+ -6.0*sin(3.0*x(1)+2.0)-12.0*sin(4.0*x(1)+3.0)
+ -20.0*sin(5.0*x(1)+4.0)-30.0*sin(6.0*x(1)+5.0));
+}
+
+
+void Domain_Shubert(RTBox box) {
+ box.lb=-10.0 ; box.ub=10.0 ;
+}
+
+double Objective_Shubert(RCRVector x) {
+ return -sin(2.0*x(0)+1.0)-2.0*sin(3.0*x(0)+2.0)-3.0*sin(4.0*x(0)+3.0)
+ -4.0*sin(5.0*x(0)+4.0)-5.0*sin(6.0*x(0)+5.0)-sin(2.0*x(1)+1.0)
+ -2.0*sin(3.0*x(1)+2.0)-3.0*sin(4.0*x(1)+3.0)-4.0*sin(5.0*x(1)+4.0)
+ -5.0*sin(6.0*x(1)+5.0);
+}
+
+void Gradient_Shubert(RCRVector x, RVector &grad) {
+ grad(0) = -2.0*cos(2.0*x(0)+1.0)-6.0*cos(3.0*x(0)+2.0)-12.0*cos(4.0*x(0)+3.0)
+ -20.0*cos(5.0*x(0)+4.0)-30.0*cos(6.0*x(0)+5.0);
+ grad(1) = -2.0*cos(2.0*x(1)+1.0)-6.0*cos(3.0*x(1)+2.0)-12.0*cos(4.0*x(1)+3.0)
+ -20.0*cos(5.0*x(1)+4.0)-30.0*cos(6.0*x(1)+5.0);
+}
+
+
+void Domain_Griewank(RTBox box) {
+ box.lb=-500 ; box.ub=700 ;
+}
+
+double Objective_Griewank(RCRVector x) {
+ double sum=0 ;
+ double prod=1 ;
+ for (int i=0 ; i<10 ; i++) {
+ sum+=x(i)*x(i) ;
+ prod*=cos(x(i)/sqrt(double(i+1))) ;
+ }
+ return sum/4000.0-prod + 1 ;
+}
+
+void Gradient_Griewank(RCRVector x, RVector &grad) {
+ double prod=1 ;
+ for (int i=0 ; i<10 ; i++) {
+ prod*=cos(x(i)/sqrt(double(i+1))) ;
+ }
+ for (int i=0 ; i<10 ; i++) {
+ grad(i)=x(i)/2000.0 + 1/sqrt(double(i+1))
+ *prod/cos(x(i)/sqrt(double(i+1)))*sin(x(i)/sqrt(double(i+1))) ;
+ }
+}
+
+
+void Domain_Levy(RTBox box) {
+ int n=(box.lb).GetLength() ;
+ switch (n) {
+ case 4 :
+ box.lb=-10.0 ; box.ub=10.0 ;
+ break ;
+ default:
+ box.lb=-5.0 ; box.ub=5.0 ;
+ }
+}
+
+double Objective_Levy(RCRVector x) {
+ int n=x.GetLength();
+ double sum=0.0;
+
+ for (int i=0 ; i<=n-2 ; i++)
+ sum+=pow(x(i)-1,2.0)*(1+pow(sin(3*pi*x(i+1)),2.0));
+ return pow(sin(3*pi*x(0)),2.0) + sum + (x(n-1)-1)*(1+pow(sin(2*pi*x(n-1)),2.0));
+}
+
+
+void Gradient_Levy(RCRVector x, RVector &grad) {
+ int n=x.GetLength();
+
+ grad(0)=6*sin(3*pi*x(0))*cos(3*pi*x(0))*pi + 2*(x(0)-1)*(1+pow(sin(3*pi*x(1)),2.0));
+
+ for (int i=1 ; i<=n-2 ; i++)
+ grad(i)=6*pow(x(i-1)-1,2.0)*sin(3*pi*x(i))*cos(3*pi*x(i))*pi
+ + 2*(x(i)-1)*(1+pow(sin(3*pi*x(i+1)),2.0)) ;
+
+ grad(n-1)=6*pow(x(n-2)-1,2.0)*sin(3*pi*x(n-1))*cos(3*pi*x(n-1))*pi
+ + 1 + pow(sin(2*pi*x(n-1)),2.0) + 4*(x(n-1)-1)*sin(2*pi*x(n-1))*cos(2*pi*x(n-1))*pi;
+}
+
+
+void Domain_Kowalik(RTBox box) {
+ box.lb=0 ; box.ub=5.0 ;
+}
+
+double Objective_Kowalik(RCRVector x) {
+ // First element in a and b is not used
+ double a[]={999,0.1957,0.1947,0.1735,0.1600,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246};
+ double b[]={999,1./0.25,1./0.5,1.,1./2,1./4,1./6,1./8,1./10,1./12,1./14,1./16};
+ double x1=x(0),x2=x(1),x3=x(2),x4=x(3);
+
+ double s=0;
+ for (int i=1 ; i<=11 ; i++)
+ s+=pow(a[i]- (x1*(b[i]*b[i] + b[i]*x2))/(b[i]*b[i] + b[i]*x3 + x4),2.0) ;
+ return s;
+}
+
+void Gradient_Kowalik(RCRVector x, RVector &grad) {
+ double a[]={999,0.1957,0.1947,0.1735,0.1600,0.0844,0.0627,0.0456,0.0342,0.0323,0.0235,0.0246};
+ double b[]={999,1./0.25,1./0.5,1.,1./2,1./4,1./6,1./8,1./10,1./12,1./14,1./16};
+ double x1=x(0),x2=x(1),x3=x(2),x4=x(3);
+ double g1=0,g2=0,g3=0,g4=0,tmp;
+ for (int i=1 ; i<=11 ; i++) {
+ tmp=a[i]- (x1*(b[i]*b[i] + b[i]*x2))/(b[i]*b[i] + b[i]*x3 + x4) ;
+ g1=g1-2*tmp*(b[i]*b[i]+b[i]*x2)/(b[i]*b[i]+b[i]*x3+x4) ;
+ g2=g2-2*tmp*(x1*b[i])/(b[i]*b[i]+b[i]*x3+x4) ;
+ g3=g3+2*tmp*x1*(b[i]*b[i]+b[i]*x2)*b[i]/pow(b[i]*b[i]+b[i]*x3+x4,2.0) ;
+ g4=g4+2*tmp*x1*(b[i]*b[i]+b[i]*x2)/pow(b[i]*b[i]+b[i]*x3+x4,2.0) ;
+ }
+ grad(0)=g1; grad(1)=g2; grad(2)=g3; grad(3)=g4;
+}
+
+
+void Domain_Camel6(RTBox box) {
+ box.lb=-5.0 ; box.ub=10.0 ;
+}
+
+double Objective_Camel6(RCRVector x) {
+ double x1=x(0),x2=x(1) ;
+ return 4.0*x1*x1-0.21E1*pow(x1,4.0)+pow(x1,6.0)/3+x1*x2-4.0*x2*x2
+ + 4.0*pow(x2,4.0);
+}
+
+void Gradient_Camel6(RCRVector x, RVector &grad) {
+ double x1=x(0),x2=x(1) ;
+ grad(0) = 8.0*x1-0.84E1*x1*x1*x1+2.0*pow(x1,5.0)+x2;
+ grad(1) = x1-8.0*x2+16.0*x2*x2*x2;
+}
+
+/** Multimodal function with a huge number of local optima **/
+
+void Domain_Rastrigin(RTBox box) {
+ box.lb=-4.12 ; box.ub=6.12;
+}
+
+double Objective_Rastrigin(RCRVector x) {
+ double sum=0.0;
+ int n=x.GetLength();
+
+ for (int i=1 ; i<=n ; i++)
+ sum+=pow(x(i-1),2.0)-10*cos(2*pi*x(i-1))+10 ;
+ return sum ;
+}
+
+void Gradient_Rastrigin(RCRVector x, RVector &grad) {
+ int n=x.GetLength();
+
+ for (int i=1 ; i<=n ; i++)
+ grad(i-1)= 2*x(i-1) + 20*sin(2*pi*x(i-1))*pi ;
+}
+
+void Domain_Schwefel(RTBox box) {
+ box.lb=-500.0;box.ub=500.0;
+}
+
+double Objective_Schwefel(RCRVector x) {
+ double sum=0.0;
+ int n=x.GetLength();
+
+ for (int i=0 ; i<n ; i++)
+ sum+=x(i)*sin(sqrt(fabs(x(i))));
+ return -sum;
+}
+
+void Gradient_Schwefel(RCRVector x, RVector &grad) {
+ int n=x.GetLength();
+
+ for (int i=0; i<n ; i++) {
+ if (x(i)>=0)
+ grad(i)=-( sin(sqrt(x(i)))+sqrt(x(i))/2.0*cos(sqrt(x(i))) );
+ else
+ grad(i)=-( sin(sqrt(-x(i)))+sqrt(-x(i))/2.0*cos(sqrt(-x(i))) );
+ }
+}
+
+/******* Difficult multimodal problem, PERM(n) *********/
+
+void Domain_Perm(RTBox box) {
+ int n=(box.lb).GetLength();
+ box.lb=double(-n) ; box.ub=double(n) ;
+}
+
+double ObjPerm(RCRVector x, double beta) {
+ int n=x.GetLength();
+ double s1=0.0;
+ for (int k=1; k<=n ; k++) {
+ double s2=0.0;
+ for (int i=1 ; i<=n; i++)
+ // s2+=(i^k+beta)*((x(i-1)/i)^k-1);
+ s2+=(pow(1.0*i,1.0*k)+beta)*(pow(x(i-1)/i,1.0*k)-1) ;
+ s1+=s2*s2;
+ }
+ return s1;
+}
+
+void GradPerm(RCRVector x, RVector &grad, double beta) {
+ int n=x.GetLength();
+ for (int j=1 ; j<=n ; j++) {
+ double s1=0.0;
+ for (int k=1 ; k<=n ; k++) {
+ double s2=0.0;
+ for (int i=1 ; i<=n; i++)
+ // s2+=(i^k+beta)*((x(i-1)/i)^k-1);
+ s2+=(pow(1.0*i,1.0*k)+beta)*(pow(x(i-1)/i,1.0*k)-1);
+ // s1+=2*s2*(j^k+beta)/(x(j-1)*j^k)*k*x(j-1)^k;
+ s1+=2*s2*(pow(1.0*j,1.0*k)+beta)/pow(1.0*j,1.0*k)*k*pow(x(j-1),k-1.0);
+ }
+ grad(j-1)=s1;
+ }
+}
+
+double Objective_Perm_4_50(RCRVector x) {
+ return ObjPerm(x,50);
+}
+
+void Gradient_Perm_4_50(RCRVector x, RVector &grad) {
+ GradPerm(x,grad,50);
+}
+
+double Objective_Perm_4_05(RCRVector x) {
+ return ObjPerm(x,0.5);
+}
+
+void Gradient_Perm_4_05(RCRVector x, RVector &grad) {
+ GradPerm(x,grad,0.5);
+}
+
+/******************* Powersum **************/
+void Domain_Powersum(RTBox box) {
+ box.lb=0.0 ; box.ub=4.0 ;
+}
+
+double Objective_Powersum(RCRVector x) {
+ int n=x.GetLength();
+ double b[]={8,18,44,114};
+
+ double s1=0.0;
+ for (int k=1 ; k<=n ; k++) {
+ double s2=0.0;
+ for (int i=1 ; i<=n ; i++)
+ s2+=pow(x(i-1),1.0*k);
+ s1+=pow(s2-b[k-1],2.0);
+ }
+ return s1;
+}
+
+void Gradient_Powersum(RCRVector x, RVector &grad) {
+ int n=x.GetLength();
+ double b[]={8,18,44,114};
+
+ for (int j=0 ; j<n ; j++) {
+ double s1=0.0;
+ for (int k=1 ; k<=n ; k++) {
+ double s2=0.0;
+ for (int i=1 ; i<=n ; i++)
+ s2+=pow(x(i-1),1.0*k);
+ s1+=2*(s2-b[k-1])*k*pow(x(j),k-1.0);
+ }
+ grad(j)=s1;
+ }
+}
+
+#endif
--- /dev/null
+#include "global.h"
+#include "tools.h"
+#include "rosen.h"
+
+int main() {
+ int dim=2;
+
+ Pdom Dom=Domain_Rosenbrock;
+ Pobj Obj=Objective_Rosenbrock;
+ Pgrad Grad=Gradient_Rosenbrock;
+
+ GlobalParams params;
+ params.det_pnts=2*dim+1; params.rnd_pnts=0;
+ params.eps_cl=0.1; params.rshift=0.3;
+ params.mu=1.0E-4; params.maxtime=5;
+
+ TBox D(dim);
+ Dom(D);
+ Global Problem(D,Obj, Grad, params);
+ RVector dummyvec(dim);
+ Problem.Search(-1, dummyvec);
+
+ cout << "Minimizers found\n";
+ Problem.DispMinimizers();
+}
--- /dev/null
+
+#include <float.h>
+#include <iostream>
+
+#include "stogo_config.h"
+#include "tools.h"
+
+Trial::Trial():xvals(0) {
+ objval=DBL_MAX;
+}
+
+Trial::Trial(int n):xvals(n) {
+ objval=DBL_MAX;
+}
+
+Trial::Trial(RCTrial tr):xvals(tr.xvals) {
+ objval=tr.objval ;
+}
+
+RCTrial Trial::operator=(RCTrial tr) {
+ xvals=tr.xvals ;
+ objval=tr.objval ;
+ return *this ;
+}
+
+ostream & operator << (ostream & os, const Trial & T) {
+ os << T.xvals << " " << "(" << T.objval << ")" << endl ;
+ return os;
+}
+
+/********************* Vanilla Box **********************/
+VBox::VBox():lb(0),ub(0) {} // Constructor
+
+VBox::VBox(int n):lb(n),ub(n) {} // Constructor
+
+VBox::VBox(RCVBox box):lb(box.lb), ub(box.ub) {} // Constructor
+
+RCVBox VBox::operator=(RCVBox box) {
+ // Copy: Box1=Box2
+ lb=box.lb; ub=box.ub;
+ return *this;
+}
+
+int VBox::GetDim() {
+ return lb.GetLength();
+}
+
+double VBox::Width(int i) {
+ // Returns the width of the i-th interval, i between [0,...,dim-1]
+ return ub(i)-lb(i);
+}
+
+void VBox::Midpoint(RCRVector x) {
+ // Returns the midpoint of Box in x
+ int n=GetDim();
+ for (int i=0 ; i<n ; i++)
+ x(i)=fabs(ub(i)-lb(i))/2 + lb(i);
+}
+
+ostream & operator << (ostream & os, const VBox & B) {
+ int n=(B.lb).GetLength() ;
+ for (int i=0 ; i<n ;i++)
+ os << '[' << B.lb(i) << "," << B.ub(i) << "]";
+ return os;
+}
+
+/************************ Trial Box ***********************/
+TBox::TBox() : VBox() {
+ // Constructor
+ fmin=DBL_MAX;
+}
+
+TBox::TBox(int n) : VBox(n) {
+ // Constructor
+ fmin=DBL_MAX;
+}
+
+TBox::TBox(RCTBox box) : VBox(box) {
+ // Constructor + Copy
+ fmin=box.fmin;
+ TList=box.TList ;
+}
+
+RCTBox TBox::operator=(RCTBox box) {
+ // Copy
+ // NB We should 'somehow' use VBox to copy lb and ub
+ // Note that assignment operators are _not_ inherited
+ lb=box.lb ; ub=box.ub ;
+ fmin=box.fmin ;
+ TList=box.TList ;
+ return *this ;
+}
+
+double TBox::GetMin() {
+ return fmin;
+}
+
+bool TBox::EmptyBox() {
+ // Returns TRUE if the list of Trials is empty
+ return TList.empty() ;
+}
+
+void TBox::AddTrial(RCTrial T) {
+ // Add a Trial to the (back of) TList
+ TList.push_back(T) ;
+ if (T.objval<fmin)
+ fmin=T.objval ;
+}
+
+void TBox::RemoveTrial(Trial &T) {
+ // Remove a trial from the (back of) box
+ T=TList.back() ;
+ TList.pop_back() ;
+}
+
+void TBox::GetLastTrial(Trial &T) {
+// Return a trial from the (back of) box
+ T=TList.back() ;
+}
+
+list<Trial>::const_iterator TBox::FirstTrial() {
+ return TList.begin();
+}
+
+list<Trial>::const_iterator TBox::LastTrial() {
+ return TList.end();
+}
+
+void TBox::GetTrial(list<Trial>::const_iterator itr, Trial &T) {
+ T.xvals=(*itr).xvals;
+ T.objval=(*itr).objval;
+}
+
+void TBox::ClearBox() {
+ TList.erase(TList.begin(), TList.end());
+ fmin=DBL_MAX;
+}
+
+bool TBox::CloseToMin(RVector &vec, double *objval, double eps_cl) {
+ // Returns TRUE if 'vec' is close to some of the trials in the box,
+ // in this case, 'vec' and 'objval' are overwritten by the Trial data
+ // otherwise 'vec' and 'objval' are not affected.
+ //
+ // Here "close" means norm2(T - some trial in box)<=eps_cl
+ //
+ // NB It might be better to accept Trial as argument instead of vector
+
+ int n=GetDim();
+ RVector x(n), y(n);
+ list<Trial>::const_iterator itr;
+ for ( itr = TList.begin(); itr != TList.end(); ++itr ) {
+ y=vec ; // NB Should be possible to avoid this copying
+ x=(*itr).xvals;
+ axpy(-1,x,y);
+ if (norm2(y)<=eps_cl) {
+ vec=x;
+ *objval=(*itr).objval;
+ return TRUE;
+ }
+ }
+ return FALSE;
+}
+
+unsigned int TBox::NStationary() {
+ // Returns the number of trials in a box
+ return TList.size() ;
+}
+
+void TBox::split(RTBox B1, RTBox B2) {
+ list<Trial>::const_iterator itr;
+ double w,m,tmp;
+ double fm1=DBL_MAX, fm2=DBL_MAX;
+ int i, k, ns;
+ int n=GetDim();
+
+ B1.lb=lb; B1.ub=ub;
+ B2.lb=lb; B2.ub=ub;
+ w=LongestSide(&i);
+ ns=TList.size();
+ switch (ns) {
+ case 0: case 1:
+ // Bisection
+ w=ub(i)-lb(i); // Length of interval
+ m=lb(i)+w/2; // Midpoint
+ B1.ub(i)=m; B2.lb(i)=m;
+ break;
+ default:
+ // More elaborate split when there are more than one trials in B
+ // See Serguies and Kaj's tech. report, page 11
+ // Compute the center point of all stationary points
+ RVector center(n), x(n), dispers(n);
+ center=0; dispers=0;
+ for ( itr = TList.begin(); itr != TList.end(); ++itr )
+ axpy(1.0, (*itr).xvals, center);
+ scal((double)(1.0/ns),center);
+
+ // Compute the relative deviations
+ for ( itr = TList.begin(); itr != TList.end(); ++itr ) {
+ for (k = 0; k<n; k++) {
+ x=(*itr).xvals;
+ dispers(k)=dispers(k)+pow(center(k)-x(k),2.0);
+ }
+ }
+ scal((double)(1.0/ns),dispers);
+
+ // i=arg max(disp)
+ tmp=dispers(0);i=0;
+ for (k=1; k<n; k++) {
+ if (dispers(k)>tmp) {
+ tmp=dispers(k);i=k;
+ }
+ }
+ B1.ub(i)=center(i) ; B2.lb(i)=center(i);
+ break;
+ }
+
+ // Split the set of trials accordingly
+ for ( itr = TList.begin(); itr != TList.end(); ++itr ) {
+ if ( B1.InsideBox((*itr).xvals) ) {
+ fm1=min(fm1,(*itr).objval);
+ B1.AddTrial(*itr);
+ }
+ else {
+ B2.AddTrial(*itr);
+ fm2=min(fm2,(*itr).objval);
+ }
+ }
+ // Set fmin of B1 and B2
+ B1.fmin=fm1 ; B2.fmin=fm2;
+}
+
+void TBox::dispTrials() {
+ // Display information about box
+#ifdef KCC
+ copy(TList.begin(), TList.end(), ostream_iterator<Trial,char>(cout));
+#else
+ copy(TList.begin(), TList.end(), ostream_iterator<Trial>(cout));
+#endif
+}
+
+ostream & operator << (ostream & os, const TBox & B) {
+ int n=(B.lb).GetLength() ;
+ for (int i=0 ; i<n ;i++)
+ os << '[' << B.lb(i) << "," << B.ub(i) << "]";
+ os << " fmin= " << B.fmin << endl;
+ return os;
+}
+
+bool TBox::InsideBox(RCRVector x) {
+ // Returns TRUE if the point X lies inside BOX, FALSE otherwise
+ int n=GetDim();
+ for (int i=0 ; i<n ; i++)
+ if (x(i)<lb(i) || x(i)>ub(i)) return FALSE;
+ return TRUE;
+}
+
+int TBox::OutsideBox(RCRVector x, RCTBox domain) {
+ // The function returns
+ // 0 if X is inside both 'box' and 'domain'
+ // 1 if X is inside 'domain' but outside 'box'
+ // 2 if X is outside both 'domain' and 'box
+
+ int n=GetDim();
+ int ins_box=1, ins_dom=1, outs=999;
+ for (int i=0 ; i<n ; i++) {
+ if (x(i)<lb(i) || x(i)>ub(i))
+ ins_box=0;
+ if (x(i)<domain.lb(i) || x(i)>domain.ub(i)) {
+ ins_dom=0; break;
+ }
+ }
+ if (ins_box==1 && ins_dom==1)
+ outs=0;
+ if (ins_box==0 && ins_dom==1)
+ outs=1;
+ if (ins_box==0 && ins_dom==0)
+ outs=2;
+ if (outs==999) {
+ // Something has gone wrong!
+ cout << "Error in OutsideBox, exiting\n";
+ exit(1);
+ }
+ return outs;
+}
+
+double TBox::ShortestSide(int *idx) {
+ // Returns the shortest side of the box and it's index
+ int n=GetDim(),j=0;
+ double tmp=ub(0)-lb(0);
+ for (int i=1 ; i<n ;i++)
+ if ( (ub(i)-lb(i))<tmp ) {
+ tmp=ub(i)-lb(i);
+ j=i;
+ }
+ *idx=j;
+ return tmp;
+}
+
+double TBox::LongestSide(int *idx) {
+ // Returns the longest side of the box and it's index
+ int n=GetDim(),j=0;
+ double tmp=ub(0)-lb(0);
+ for (int i=1 ; i<n ;i++)
+ if ( (ub(i)-lb(i))>tmp ) {
+ tmp=ub(i)-lb(i);
+ j=i;
+ }
+ *idx=j;
+ return tmp;
+}
+
+double TBox::ClosestSide(RCRVector x) {
+ // Returns the shortest distance from point X to the box B.
+
+ // Warning: The output of this functon is nonsense if the
+ // point X lies outside B. Should we try to detect this case?
+
+ double dist, tmp ;
+ int n=GetDim();
+ dist=DBL_MAX;
+ for (int i=0 ; i<n ; i++) {
+ tmp = min( x(i)-lb(i), ub(i)-x(i) );
+ dist = min(dist, tmp);
+ }
+ return dist;
+}
+
+double TBox::FarthestSide(RCRVector x) {
+ // Returns the longest distance from point X to the box B.
+ // Same comment apply here as in ClosestSide(X)
+
+ double dist, tmp;
+ int n=GetDim();
+ dist=DBL_MIN;
+ for (int i=0 ; i<n ; i++) {
+ tmp = max( x(i)-lb(i), ub(i)-x(i) );
+ dist = max(dist, tmp);
+ }
+ return dist;
+}
+
+bool TBox::Intersection(RCRVector x, RCRVector h, RCRVector z) {
+// Intersection of a line with a box
+// The line is described by a vector 'h' and a point 'x' and
+// the point of intersection is returned in 'z'
+// Note that h+x must lie outside of box
+//
+// Known problem:
+// Due to round of errors the algorithm can fail to find an intersection
+// The caller is notified and should act accordingly
+//
+// The routine returns FALSE if no intersection was found, TRUE otherwise
+
+ int n=GetDim();
+ RVector tmpV(n);
+ bool done;
+ int i, j, k, isect;
+ double alpha, gamma;
+
+ i=0; done=FALSE;
+ while (i<n && done==FALSE) {
+ if (h(i)==0) {
+ z(i)=x(i);
+ break;
+ }
+ for (k=1; k<=2; k++) {
+ if (k==1)
+ alpha=lb(i);
+ else
+ alpha=ub(i);
+ gamma=(alpha-x(i))/h(i);
+ z(i)=alpha;
+ isect=1;
+ for (j=0; j<n; j++) {
+ if (j != i) {
+ z(j)=x(j)+gamma*h(j);
+ if (z(j)<lb(j) || z(j)>ub(j)) {
+ isect=0;
+ break;
+ }
+ }
+ }
+ copy(z,tmpV); axpy(-1.0,x,tmpV); // tmpV=z-x
+ if (isect==1 && dot(tmpV,h)>0) {
+ done=TRUE; break;
+ }
+ }
+ i++;
+ }
+ return done;
+}
+
+double TBox::LowerBound(double maxgrad) {
+// Lower bound estimation
+ double lb=fmin ;
+ double f1,f2,est ;
+ list<Trial>::const_iterator itr1,itr2 ;
+
+ int n=GetDim();
+ RVector x1(n), x2(n) ;
+#ifndef LB2
+ for ( itr1 = TList.begin(); itr1 != TList.end(); ++itr1 ) {
+ itr2=itr1 ;
+ while (++itr2 != TList.end()) {
+ x1=(*itr1).xvals ; f1=(*itr1).objval ;
+ x2=(*itr2).xvals ; f2=(*itr2).objval ;
+ axpy(-1.0,x2,x1) ; // x1=x1-x2
+ est=0.5*(f1+f2-maxgrad*norm2(x1)) ;
+ lb=min(lb,est) ;
+ // cout << "est=" << est << " x1=" << x1 << " x2=" << x2 << endl ;
+ }
+ }
+#endif
+#ifdef LB2
+ for ( itr1 = TList.begin(); itr1 != TList.end(); ++itr1 ) {
+ // Use also max distance to border
+ x1=(*itr1).xvals; f1=(*itr1).objval;
+ est=f1-FarthestSide(x1)*maxgrad;
+ lb=min(lb,est);
+ }
+#endif
+ return lb ;
+}
--- /dev/null
+
+// Various datastructures and functions used by the global optimizer
+
+#ifndef TOOLS_H
+#define TOOLS_H
+
+#include <float.h>
+#include <iostream>
+
+#include <algorithm>
+#include <iterator>
+#include <list>
+
+#include "linalg.h"
+
+#ifndef FALSE
+const int FALSE=(1==0); // boolean FALSE
+#endif
+#ifndef FALSE
+const int TRUE=(1==1); // boolean TRUE
+#endif
+
+typedef const class Trial CTrial;
+typedef CTrial& RCTrial;
+typedef CTrial* PCTrial;
+
+class Trial{
+public:
+ RVector xvals;
+ double objval;
+
+ Trial();
+ Trial(int);
+ Trial(RCTrial); // Copy constructor
+ RCTrial operator=(RCTrial) ; // assignment operator
+ friend ostream & operator << (ostream &, RCTrial) ;
+};
+
+class TrialGT : public unary_function<Trial, bool>
+// Predicate for Trial (needed for remove_if)
+{
+public:
+ explicit TrialGT(double val) : _val(val) {}
+ bool operator()(Trial& foo) {
+ return foo.objval > _val;
+ }
+private:
+ double _val;
+};
+
+typedef class VBox& RVBox;
+typedef const class VBox CVBox;
+typedef CVBox* PCVBox;
+typedef CVBox& RCVBox;
+
+class VBox{
+public:
+ RVector lb, ub;
+
+ VBox(); // Construct a box
+ VBox(int);
+ VBox(RCVBox); // Copy constructor
+ RCVBox operator=(RCVBox); // assignment operator
+
+ int GetDim(); // Returns the dimension of the box
+ double Width(int) ; // Returns the width of the i-th interval
+ void Midpoint(RCRVector); // Returns the midpoint
+
+ friend ostream & operator << (ostream &, const VBox &);
+};
+
+typedef class TBox& RTBox;
+typedef const class TBox CTBox;
+typedef CTBox* PCTBox;
+typedef CTBox& RCTBox;
+
+class TBox: public VBox {
+public:
+ double fmin; // Smallest function value found so far
+ list<Trial> TList; // List of trials
+
+ TBox(); // Construct a box
+ TBox(int);
+ TBox(RCTBox); // Copy constructor
+
+ RCTBox operator=(RCTBox); // assignment operator
+
+ double GetMin(); // Returns 'fmin'
+ bool EmptyBox(); // Returns TRUE if Box contains no trials
+ void AddTrial(RCTrial); // Add a trial to the (back of) box
+ void RemoveTrial(Trial &); // Remove a trial from the (back of) box
+ void GetLastTrial(Trial &); // Return a trial from the back of the box
+
+ list<Trial>::const_iterator FirstTrial();
+ list<Trial>::const_iterator LastTrial();
+
+ void GetTrial(list<Trial>::const_iterator, Trial&);
+ void ClearBox();
+ bool CloseToMin(RVector&, double*, double);
+
+ unsigned int NStationary(); // Returns the number of stationary points
+
+ void split(RTBox, RTBox); // Split a box
+ void dispTrials();
+
+ bool InsideBox(RCRVector);
+ int OutsideBox(RCRVector, RCTBox);
+ double ShortestSide(int*); // Returns the shortest side
+ double LongestSide(int*); // Returns the longest side
+ double ClosestSide(RCRVector x);
+ double FarthestSide(RCRVector x);
+ bool Intersection(RCRVector, RCRVector, RCRVector);
+ double LowerBound(double);
+
+ bool operator<(const TBox & x) const {return (fmin>x.fmin);}
+ friend ostream & operator << (ostream &, const TBox &);
+};
+
+#endif
--- /dev/null
+#include <unistd.h>
+
+#include "global.h"
+#include "tools.h"
+
+#include "linalg.h"
+#include "tools.h"
+#include "stogo_config.h"
+
+void Domain_Mine(RTBox box) {
+ box.lb=-3.0 ; box.ub=3.0;
+}
+
+static int cnt = 0;
+double Objective_Mine(RCRVector xy) {
+ double x = xy(0);
+ double y = xy(1);
+ double f = ((x*x)*(4-2.1*(x*x)+((x*x)*(x*x))/3) + x*y + (y*y)*(-4+4*(y*y)));
+ printf("feval:, %d, %g, %g, %g\n", ++cnt, x,y, f);
+ return f;
+}
+
+void Gradient_Mine(RCRVector xy, RVector &grad) {
+ double x = xy(0);
+ double y = xy(1);
+ grad(0) = /* df/dx */
+ ((2*x)*(4-2.1*(x*x)+((x*x)*(x*x))/3)
+ + (x*x)*(-4.2*x+4*(x*(x*x))/3)
+ + y);
+ grad(1) = /* df/dy */
+ (x + (2*y)*(-4+4*(y*y)) + (y*y)*(8*(y)));
+}
+
+
+int main(int argc, char **argv) {
+ int dim=2;
+
+ Pdom Dom=Domain_Mine;
+ Pobj Obj=Objective_Mine;
+ Pgrad Grad=Gradient_Mine;
+
+ GlobalParams params;
+ params.det_pnts=2*dim+1; params.rnd_pnts=0;
+ params.eps_cl=0.1; params.rshift=0.3;
+ params.mu=1.0E-4; params.maxtime=0;
+
+ params.maxeval = argc < 2 ? 100 : atoi(argv[1]);
+
+ TBox D(dim);
+ Dom(D);
+ Global Problem(D,Obj, Grad, params);
+ RVector dummyvec(dim);
+ Problem.Search(-1, dummyvec);
+
+ cout << "Minimizers found\n";
+ Problem.DispMinimizers();
+
+ cout << "dummyvec: " << dummyvec << "\n";
+
+ double val = Problem.OneMinimizer(dummyvec);
+ cout << "one minimizer: " << val << ": " << dummyvec << "\n";
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "stogo.h"
+
+/* has two global minima at (0.09,-0.71) and (-0.09,0.71), plus
+ 4 additional local minima */
+static int cnt=0;
+double tst_obj(int n, const double *xy, double *g, void *unused)
+{
+ double x, y, f;
+ x = xy[0];
+ y = xy[1];
+ f = ((x*x)*(4-2.1*(x*x)+((x*x)*(x*x))/3) + x*y + (y*y)*(-4+4*(y*y)));
+ printf("feval:, %d, %g, %g, %g\n", ++cnt, x,y, f);
+ if (g) {
+ g[0] = /* df/dx */
+ ((2*x)*(4-2.1*(x*x)+((x*x)*(x*x))/3)
+ + (x*x)*(-4.2*x+4*(x*(x*x))/3)
+ + y);
+ g[1] = /* df/dy */
+ (x + (2*y)*(-4+4*(y*y)) + (y*y)*(8*(y)));
+ }
+ return f;
+}
+
+int main(int argc, char **argv)
+{
+ int n = 2;
+ double x[2], l[2], u[2];
+ long int maxits = 0, maxtim = 0;
+ int info;
+ double fmin;
+
+ maxits = argc < 2 ? 100 : atoi(argv[1]);
+
+ l[0] = -3; l[1] = -3;
+ u[0] = 3; u[1] = 3;
+
+ info = stogo_minimize(n, tst_obj, NULL, x, &fmin, l, u, maxits, maxtim);
+
+ printf("min f = %g at (%g,%g) after %d evals, return value %d\n",
+ fmin, x[0], x[1], cnt, info);
+
+ return EXIT_SUCCESS;
+}
~ianmdlvl / nlopt.git / commitdiff