From: stevenj Date: Wed, 22 Aug 2007 22:10:32 +0000 (-0400) Subject: initial checkin X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=0fd730a0ee0f4197821b5185351fd990ff4e3c40;p=nlopt.git initial checkin darcs-hash:20070822221032-c8de0-6fe4c3d6af0aa9a73f141cdbda9fabc08fe7805a.gz --- 0fd730a0ee0f4197821b5185351fd990ff4e3c40 diff --git a/direct/AUTHORS b/direct/AUTHORS new file mode 100644 index 0000000..b676a64 --- /dev/null +++ b/direct/AUTHORS @@ -0,0 +1,2 @@ +C conversion: Steven G. Johnson (stevenj@alum.mit.edu) +Original Fortran code: Joerg.M.Gablonsky (jmgablon@mailandnews.com) diff --git a/direct/COPYRIGHT b/direct/COPYRIGHT new file mode 100644 index 0000000..4b4a9c3 --- /dev/null +++ b/direct/COPYRIGHT @@ -0,0 +1,28 @@ +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. diff --git a/direct/DIRect.c b/direct/DIRect.c new file mode 100644 index 0000000..f208332 --- /dev/null +++ b/direct/DIRect.c @@ -0,0 +1,743 @@ +/* DIRect.f -- translated by f2c (version 20050501). + + f2c output hand-cleaned by SGJ (August 2007). +*/ + +#include +#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_ */ + diff --git a/direct/DIRparallel.c b/direct/DIRparallel.c new file mode 100644 index 0000000..aba092e --- /dev/null +++ b/direct/DIRparallel.c @@ -0,0 +1,381 @@ +/* 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_ */ diff --git a/direct/DIRserial.c b/direct/DIRserial.c new file mode 100644 index 0000000..bc45a51 --- /dev/null +++ b/direct/DIRserial.c @@ -0,0 +1,145 @@ +/* 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_ */ diff --git a/direct/DIRsubrout.c b/direct/DIRsubrout.c new file mode 100644 index 0000000..4d8a45e --- /dev/null +++ b/direct/DIRsubrout.c @@ -0,0 +1,1715 @@ +/* DIRsubrout.f -- translated by f2c (version 20050501). + + f2c output hand-cleaned by SGJ (August 2007). +*/ + +#include +#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__ */ + diff --git a/direct/direct-internal.h b/direct/direct-internal.h new file mode 100644 index 0000000..84fd63e --- /dev/null +++ b/direct/direct-internal.h @@ -0,0 +1,130 @@ +#ifndef DIRECT_INTERNAL_H +#define DIRECT_INTERNAL_H + +#include +#include +#include + +#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 */ diff --git a/direct/direct.h b/direct/direct.h new file mode 100644 index 0000000..41ed8c4 --- /dev/null +++ b/direct/direct.h @@ -0,0 +1,59 @@ +#ifndef DIRECT_H +#define DIRECT_H + +#include + +#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 */ diff --git a/direct/direct_wrap.c b/direct/direct_wrap.c new file mode 100644 index 0000000..406a44d --- /dev/null +++ b/direct/direct_wrap.c @@ -0,0 +1,97 @@ +/* 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; +} diff --git a/direct/tstc.c b/direct/tstc.c new file mode 100644 index 0000000..23a161a --- /dev/null +++ b/direct/tstc.c @@ -0,0 +1,41 @@ +#include +#include + +#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; +} diff --git a/direct/userguide.ps.gz b/direct/userguide.ps.gz new file mode 100644 index 0000000..c3a1655 Binary files /dev/null and b/direct/userguide.ps.gz differ diff --git a/stogo/README b/stogo/README new file mode 100644 index 0000000..71c6c62 --- /dev/null +++ b/stogo/README @@ -0,0 +1,21 @@ +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. diff --git a/stogo/global.cc b/stogo/global.cc new file mode 100644 index 0000000..9496358 --- /dev/null +++ b/stogo/global.cc @@ -0,0 +1,315 @@ +/* + 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 +#include + +#include +#include +#include + +#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(cout)); +} + +double Global::OneMinimizer(RCRVector x) { + if (NoMinimizers()) return 0.0; + for (int i=0;i +//#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 SolSet; + list::const_iterator titr; + priority_queue CandSet; + priority_queue 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 + diff --git a/stogo/linalg.cc b/stogo/linalg.cc new file mode 100644 index 0000000..04916cf --- /dev/null +++ b/stogo/linalg.cc @@ -0,0 +1,216 @@ +/* + Temporary implementation of vector and matrix classes + No attempt is made to check if the function arguments are valid +*/ + +#include +#include // 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;i0) 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 +using namespace std; +#include // for sqrt() +#include + +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 diff --git a/stogo/local.cc b/stogo/local.cc new file mode 100644 index 0000000..615a439 --- /dev/null +++ b/stogo/local.cc @@ -0,0 +1,322 @@ + +/* + Local search - A trust region algorithm with BFGS update. +*/ + +#include +#include + +#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 ; ik_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 ; +} diff --git a/stogo/local.h b/stogo/local.h new file mode 100644 index 0000000..21136dd --- /dev/null +++ b/stogo/local.h @@ -0,0 +1,25 @@ +/* + 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 diff --git a/stogo/prog.cc b/stogo/prog.cc new file mode 100644 index 0000000..1cb941a --- /dev/null +++ b/stogo/prog.cc @@ -0,0 +1,301 @@ +/* + 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=" <> 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=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 +#include + +#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::const_iterator TBox::FirstTrial() { + return TList.begin(); +} + +list::const_iterator TBox::LastTrial() { + return TList.end(); +} + +void TBox::GetTrial(list::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::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::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; ktmp) { + 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(cout)); +#else + copy(TList.begin(), TList.end(), ostream_iterator(cout)); +#endif +} + +ostream & operator << (ostream & os, const TBox & B) { + int n=(B.lb).GetLength() ; + for (int i=0 ; itmp ) { + 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 ; iub(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::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 ; +} diff --git a/stogo/tools.h b/stogo/tools.h new file mode 100644 index 0000000..81d8dc3 --- /dev/null +++ b/stogo/tools.h @@ -0,0 +1,119 @@ + +// Various datastructures and functions used by the global optimizer + +#ifndef TOOLS_H +#define TOOLS_H + +#include +#include + +#include +#include +#include + +#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 +// 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 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::const_iterator FirstTrial(); + list::const_iterator LastTrial(); + + void GetTrial(list::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 diff --git a/stogo/tst.cc b/stogo/tst.cc new file mode 100644 index 0000000..b2721ac --- /dev/null +++ b/stogo/tst.cc @@ -0,0 +1,62 @@ +#include + +#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"; +} diff --git a/stogo/tstc.c b/stogo/tstc.c new file mode 100644 index 0000000..c321b3e --- /dev/null +++ b/stogo/tstc.c @@ -0,0 +1,46 @@ +#include +#include + +#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; +}