From: stevenj Date: Fri, 24 Aug 2007 16:35:15 +0000 (-0400) Subject: added common random-number generation and timer utilities X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=4cdb51b10defaae0e766f84e4d784baf469f7e20;p=nlopt.git added common random-number generation and timer utilities darcs-hash:20070824163515-c8de0-a83acc3d0d58d0e818fa27585dcefbc0c6a578b4.gz --- diff --git a/Makefile.am b/Makefile.am index 0fa9377..a54a7fe 100644 --- a/Makefile.am +++ b/Makefile.am @@ -3,11 +3,11 @@ lib_LTLIBRARIES = libnlopt.la ACLOCAL_AMFLAGS=-I ./m4 -SUBDIRS=lbfgs subplex direct stogo api . test +SUBDIRS= util lbfgs subplex direct stogo api . test EXTRA_DIST=COPYRIGHT autogen.sh nlopt.pc.in m4 libnlopt_la_SOURCES = -libnlopt_la_LIBADD = lbfgs/liblbfgs.la subplex/libsubplex.la direct/libdirect.la stogo/libstogo.la api/libapi.la +libnlopt_la_LIBADD = lbfgs/liblbfgs.la subplex/libsubplex.la direct/libdirect.la stogo/libstogo.la api/libapi.la util/libutil.la libnlopt_la_LDFLAGS = -no-undefined -version-info @SHARED_VERSION_INFO@ diff --git a/api/Makefile.am b/api/Makefile.am index 130dc35..71a0fc2 100644 --- a/api/Makefile.am +++ b/api/Makefile.am @@ -1,4 +1,4 @@ -AM_CPPFLAGS = -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/lbfgs -I$(top_srcdir)/subplex +AM_CPPFLAGS = -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/lbfgs -I$(top_srcdir)/subplex -I$(top_srcdir)/util include_HEADERS = nlopt.h noinst_LTLIBRARIES = libapi.la diff --git a/api/nlopt.c b/api/nlopt.c index 8ee08a0..8935286 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -2,6 +2,7 @@ #include #include "nlopt.h" +#include "nlopt-util.h" #include "config.h" static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][128] = { @@ -18,6 +19,16 @@ const char *nlopt_algorithm_name(nlopt_algorithm a) return nlopt_algorithm_names[a]; } +static int nlopt_srand_called = 0; +void nlopt_srand(unsigned long seed) { + nlopt_srand_called = 1; + nlopt_init_genrand(seed); +} + +void nlopt_srand_time() { + nlopt_srand(nlopt_time_seed()); +} + static int my_isinf(double x) { return x == HUGE_VAL #ifdef HAVE_ISINF @@ -89,6 +100,10 @@ nlopt_result nlopt_minimize( d.lb = lb; d.ub = ub; + /* make sure rand generator is inited */ + if (!nlopt_srand_called) + nlopt_srand_time(); /* default is non-deterministic */ + /* check bound constraints */ for (i = 0; i < n; ++i) if (lb[i] > ub[i] || x[i] < lb[i] || x[i] > ub[i]) diff --git a/api/nlopt.h b/api/nlopt.h index 5676706..c480951 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -48,6 +48,9 @@ extern nlopt_result nlopt_minimize( double xtol_rel, const double *xtol_abs, int maxeval, double maxtime); +extern void nlopt_srand(unsigned long seed); +extern void nlopt_srand_time(void); + #ifdef __cplusplus } /* extern "C" */ #endif /* __cplusplus */ diff --git a/api/nlopt_minimize.3 b/api/nlopt_minimize.3 index 885a3c4..2a062bf 100644 --- a/api/nlopt_minimize.3 +++ b/api/nlopt_minimize.3 @@ -167,12 +167,14 @@ NaN is specified by the macro NAN in math.h (for ANSI C99). .SH ALGORITHMS The .I algorithm -parameter can take on any of the following values: +parameter specifies the optimization algorithm (for more detail on +these, see the README files in the source-code subdirectories), and +can take on any of the following values: .TP .B NLOPT_GLOBAL_DIRECT Perform a global derivative-free optimization using the DIRECT search -algorithm by Jones et al. See direct/README. Supports arbitrary -nonlinear constraints as described above. +algorithm by Jones et al. Supports arbitrary nonlinear constraints as +described above. .TP .B NLOPT_GLOBAL_DIRECT_L Perform a global derivative-free optimization using the DIRECT-L @@ -194,6 +196,12 @@ works (both for simple bound constraints via and .I ub as well as nonlinear constraints as described above). +.TP +.B NLOPT_GLOBAL_STOGO +Global optimization using the StoGO algorithm by Madsen et al. StoGO +exploits gradient information (which must be supplied by the +objective) for its local searches, and performs the global search by a +stochastic branch-and-bound technique. This stochastic .SH RETURN VALUE The value returned is one of the following enumerated constants. (Positive return values indicate successful termination, while negative diff --git a/configure.ac b/configure.ac index 92fcebe..c203bb4 100644 --- a/configure.ac +++ b/configure.ac @@ -22,13 +22,21 @@ AC_PROG_LIBTOOL dnl Checks for typedefs, structures, and compiler characteristics. AC_HEADER_STDC AC_HEADER_TIME -AC_CHECK_HEADERS([unistd.h getopt.h]) +AC_CHECK_HEADERS([unistd.h getopt.h stdint.h]) AC_C_CONST AC_C_INLINE +dnl find 32-bit unsigned integer type for random-number generator +AC_CHECK_SIZEOF(unsigned int) +AC_CHECK_SIZEOF(unsigned long) +AC_CHECK_TYPES(uint32_t, [], [], [$ac_includes_default +#ifdef HAVE_STDINT_H +# include +#endif]) + dnl Checks for libraries and functions AC_CHECK_LIB(m, sin) -AC_CHECK_FUNCS([BSDgettimeofday gettimeofday]) +AC_CHECK_FUNCS([BSDgettimeofday gettimeofday time]) AC_MSG_CHECKING([for isnan]) AC_TRY_LINK([#include @@ -79,6 +87,7 @@ AC_CONFIG_FILES([ Makefile nlopt.pc api/Makefile + util/Makefile direct/Makefile stogo/Makefile lbfgs/Makefile diff --git a/stogo/Makefile.am b/stogo/Makefile.am index f36a293..287f03a 100644 --- a/stogo/Makefile.am +++ b/stogo/Makefile.am @@ -1,3 +1,5 @@ +AM_CPPFLAGS = -I$(top_srcdir)/util + noinst_LTLIBRARIES = libstogo.la libstogo_la_SOURCES = global.cc linalg.cc local.cc stogo.cc tools.cc \ global.h linalg.h local.h stogo_config.h stogo.h tools.h diff --git a/stogo/global.cc b/stogo/global.cc index 80decbe..d133371 100644 --- a/stogo/global.cc +++ b/stogo/global.cc @@ -19,9 +19,10 @@ #include "stogo_config.h" #include "global.h" #include "local.h" +#include "nlopt-util.h" // Timer stuff -time_t StartTime; +double StartTime; double MacEpsilon ; int FC=0, GC=0 ; @@ -93,8 +94,7 @@ void Global::FillRandom(RTBox SampleBox, RTBox box) { tmpTrial.objval=DBL_MAX; for (int i=1 ; i<=rnd_pnts ; i++) { for (int dir=0 ; dir #endif -#if TIME_WITH_SYS_TIME -# include -# include -#else -# if HAVE_SYS_TIME_H -# include -# else -# include -# endif -#endif #include "nlopt.h" +#include "nlopt-util.h" #include "testfuncs.h" -static double urand(double a, double b) -{ - return a + (rand() * (b - a) / RAND_MAX); -} - static nlopt_algorithm algorithm = NLOPT_GLOBAL_DIRECT; static double ftol_rel = 0, ftol_abs = 0, xtol_rel = 0; static int maxeval = 1000; static double maxtime = 0.0; +static void listalgs(FILE *f) +{ + int i; + fprintf(f, "Available algorithms:\n"); + for (i = 0; i < NLOPT_NUM_ALGORITHMS; ++i) + fprintf(f, " %2d: %s\n", i, nlopt_algorithm_name((nlopt_algorithm) i)); +} + +static void listfuncs(FILE *f) +{ + int i; + fprintf(f, "Available objective functions:\n"); + for (i = 0; i < NTESTFUNCS; ++i) + fprintf(f, " %2d: %s (%d dims)\n", i, testfuncs[i].name, testfuncs[i].n); +} + static int test_function(int ifunc) { testfunc func; int i; double *x, fmin, f0; nlopt_result ret; + double start = nlopt_seconds(); if (ifunc < 0 || ifunc >= NTESTFUNCS) { fprintf(stderr, "testopt: invalid function %d\n", ifunc); + listfuncs(stderr); return 0; } func = testfuncs[ifunc]; @@ -61,7 +65,7 @@ static int test_function(int ifunc) func.name, func.n, nlopt_algorithm_name(algorithm)); printf("Starting guess x = ["); for (i = 0; i < func.n; ++i) - printf(" %g", x[i] = urand(func.lb[i], func.ub[i])); + printf(" %g", x[i] = nlopt_urand(func.lb[i], func.ub[i])); printf("]\n"); f0 = func.f(func.n, x, x + func.n, func.f_data); printf("Starting function value = %g\n", f0); @@ -85,6 +89,7 @@ static int test_function(int ifunc) x, &fmin, HUGE_VAL, ftol_rel, ftol_abs, xtol_rel, NULL, maxeval, maxtime); + printf("finished after %g seconds.\n", nlopt_seconds() - start); printf("return code %d from nlopt_minimize\n", ret); if (ret < 0) { fprintf(stderr, "testopt: error in nlopt_minimize\n"); @@ -110,6 +115,7 @@ static void usage(FILE *f) fprintf(f, "Usage: testopt [OPTIONS]\n" "Options:\n" " -h : print this help\n" + " -L : list available algorithms and objective functions\n" " -v : verbose mode\n" " -r : use random seed for starting guesses\n" " -a : use optimization algorithm \n" @@ -122,17 +128,21 @@ int main(int argc, char **argv) { int c; - srand((unsigned) time(NULL)); + nlopt_srand_time(); testfuncs_verbose = 0; if (argc <= 1) usage(stdout); - while ((c = getopt(argc, argv, "hvra:o:e:")) != -1) + while ((c = getopt(argc, argv, "hLvra:o:e:")) != -1) switch (c) { case 'h': usage(stdout); return EXIT_SUCCESS; + case 'L': + listalgs(stdout); + listfuncs(stdout); + return EXIT_SUCCESS; case 'v': testfuncs_verbose = 1; break; @@ -143,6 +153,7 @@ int main(int argc, char **argv) c = atoi(optarg); if (c < 0 || c >= NLOPT_NUM_ALGORITHMS) { fprintf(stderr, "testopt: invalid algorithm %d\n", c); + listalgs(stderr); return EXIT_FAILURE; } algorithm = (nlopt_algorithm) c; diff --git a/util/Makefile.am b/util/Makefile.am new file mode 100644 index 0000000..0bbe4b2 --- /dev/null +++ b/util/Makefile.am @@ -0,0 +1,2 @@ +noinst_LTLIBRARIES = libutil.la +libutil_la_SOURCES = mt19937ar.c timer.c nlopt-util.h diff --git a/util/mt19937ar.c b/util/mt19937ar.c new file mode 100644 index 0000000..c499573 --- /dev/null +++ b/util/mt19937ar.c @@ -0,0 +1,223 @@ +/* + A C-program for MT19937, with initialization improved 2002/1/26. + Coded by Takuji Nishimura and Makoto Matsumoto. + + Before using, initialize the state by using init_genrand(seed) + or init_by_array(init_key, key_length). + + Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, + All rights reserved. + + Modified 2007 by Steven G. Johnson for use with NLopt (to avoid + namespace pollution, use uint32_t instead of unsigned long, + and add the urand function). + + Redistribution and use in source and binary forms, with or without + modification, are permitted provided that the following conditions + are met: + + 1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + + 3. The names of its contributors may not be used to endorse or promote + products derived from this software without specific prior written + permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS + "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + + Any feedback is very welcome. + http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html + email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space) +*/ + +#include "nlopt-util.h" +#include "config.h" + +#if defined(HAVE_STDINT_H) +# include +#endif + +#ifndef HAVE_UINT32_T +# if SIZEOF_UNSIGNED_LONG == 4 + typedef unsigned long uint32_t; +# elif SIZEOF_UNSIGNED_INT == 4 + typedef unsigned int uint32_t; +# else +# error No 32-bit unsigned integer type +# endif +#endif + +/* Period parameters */ +#define N 624 +#define M 397 +#define MATRIX_A 0x9908b0dfUL /* constant vector a */ +#define UPPER_MASK 0x80000000UL /* most significant w-r bits */ +#define LOWER_MASK 0x7fffffffUL /* least significant r bits */ + +static uint32_t mt[N]; /* the array for the state vector */ +static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */ + +/* initializes mt[N] with a seed */ +void nlopt_init_genrand(unsigned long s) +{ + mt[0]= s & 0xffffffffUL; + for (mti=1; mti> 30)) + mti); + /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */ + /* In the previous versions, MSBs of the seed affect */ + /* only MSBs of the array mt[]. */ + /* 2002/01/09 modified by Makoto Matsumoto */ + mt[mti] &= 0xffffffffUL; + /* for >32 bit machines */ + } +} + +/* generates a random number on [0,0xffffffff]-interval */ +static uint32_t nlopt_genrand_int32(void) +{ + uint32_t y; + static uint32_t mag01[2]={0x0UL, MATRIX_A}; + /* mag01[x] = x * MATRIX_A for x=0,1 */ + + if (mti >= N) { /* generate N words at one time */ + int kk; + + if (mti == N+1) /* if init_genrand() has not been called, */ + nlopt_init_genrand(5489UL); /* a default initial seed is used */ + + for (kk=0;kk> 1) ^ mag01[y & 0x1UL]; + } + for (;kk> 1) ^ mag01[y & 0x1UL]; + } + y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK); + mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL]; + + mti = 0; + } + + y = mt[mti++]; + + /* Tempering */ + y ^= (y >> 11); + y ^= (y << 7) & 0x9d2c5680UL; + y ^= (y << 15) & 0xefc60000UL; + y ^= (y >> 18); + + return y; +} + +#if 0 /* not used in NLopt */ + +/* initialize by an array with array-length */ +/* init_key is the array for initializing keys */ +/* key_length is its length */ +/* slight change for C++, 2004/2/26 */ +static void nlopt_init_by_array(uint32_t init_key[], int key_length) +{ + int i, j, k; + nlopt_init_genrand(19650218UL); + i=1; j=0; + k = (N>key_length ? N : key_length); + for (; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL)) + + init_key[j] + j; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; j++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + if (j>=key_length) j=0; + } + for (k=N-1; k; k--) { + mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL)) + - i; /* non linear */ + mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */ + i++; + if (i>=N) { mt[0] = mt[N-1]; i=1; } + } + + mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */ +} + +/* generates a random number on [0,0x7fffffff]-interval */ +static long nlopt_genrand_int31(void) +{ + return (long)(nlopt_genrand_int32()>>1); +} + +/* generates a random number on [0,1]-real-interval */ +static double nlopt_genrand_real1(void) +{ + return nlopt_genrand_int32()*(1.0/4294967295.0); + /* divided by 2^32-1 */ +} + +/* generates a random number on [0,1)-real-interval */ +static double nlopt_genrand_real2(void) +{ + return nlopt_genrand_int32()*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +/* generates a random number on (0,1)-real-interval */ +static double nlopt_genrand_real3(void) +{ + return (((double)nlopt_genrand_int32()) + 0.5)*(1.0/4294967296.0); + /* divided by 2^32 */ +} + +#endif + +/* generates a random number on [0,1) with 53-bit resolution*/ +static double nlopt_genrand_res53(void) +{ + uint32_t a=nlopt_genrand_int32()>>5, b=nlopt_genrand_int32()>>6; + return(a*67108864.0+b)*(1.0/9007199254740992.0); +} +/* These real versions are due to Isaku Wada, 2002/01/09 added */ + +/* generate uniform random number in [a,b) with 53-bit resolution, + added by SGJ */ +double nlopt_urand(double a, double b) +{ + return(a + (b - a) * nlopt_genrand_res53()); +} + +#if 0 +#include +int main(void) +{ + int i; + uint32_t init[4]={0x123, 0x234, 0x345, 0x456}, length=4; + init_by_array(init, length); + printf("1000 outputs of nlopt_genrand_int32()\n"); + for (i=0; i<1000; i++) { + printf("%10lu ", nlopt_genrand_int32()); + if (i%5==4) printf("\n"); + } + printf("\n1000 outputs of genrand_real2()\n"); + for (i=0; i<1000; i++) { + printf("%10.8f ", genrand_real2()); + if (i%5==4) printf("\n"); + } + return 0; +} +#endif diff --git a/util/nlopt-util.h b/util/nlopt-util.h new file mode 100644 index 0000000..895db62 --- /dev/null +++ b/util/nlopt-util.h @@ -0,0 +1,21 @@ +#ifndef NLOPT_UTIL_H +#define NLOPT_UTIL_H + +#ifdef __cplusplus +extern "C" +{ +#endif /* __cplusplus */ + +/* seconds timer */ +extern double nlopt_seconds(void); +extern unsigned long nlopt_time_seed(void); + +/* pseudorandom number generation by Mersenne twister algorithm */ +extern void nlopt_init_genrand(unsigned long s); +extern double nlopt_urand(double a, double b); + +#ifdef __cplusplus +} /* extern "C" */ +#endif /* __cplusplus */ + +#endif diff --git a/util/timer.c b/util/timer.c new file mode 100644 index 0000000..843ba03 --- /dev/null +++ b/util/timer.c @@ -0,0 +1,72 @@ +#include "nlopt-util.h" +#include "config.h" + +#if TIME_WITH_SYS_TIME +# include +# include +#else +# if HAVE_SYS_TIME_H +# include +# else +# include +# endif +#endif + +#if defined(_WIN32) || defined(__WIN32__) +# include +#endif + +/* return time in seconds since some arbitrary point in the past */ +double nlopt_seconds(void) +{ + static int start_inited = 0; /* whether start time has been initialized */ +#if defined(HAVE_GETTIMEOFDAY) + static struct timeval start; + struct timeval tv; + if (!start_inited) { + start_inited = 1; + gettimeofday(&start, NULL); + } + gettimeofday(&tv, NULL); + return (tv.tv_sec - start.tv_sec) + 1.e-6 * (tv.tv_usec - start.tv_sec); +#elif defined(HAVE_TIME) + return time(NULL); +#elif defined(_WIN32) || defined(__WIN32__) + static ULONGLONG start; + FILETIME ft; + if (!start_inited) { + start_inited = 1; + GetSystemTimeAsFileTime(&ft); + start = (((ULONGLONG) ft.dwHighDateTime) << 32) + ft.dwLowDateTime; + } + GetSystemTimeAsFileTime(&ft); + return 100e-9 * (((((ULONGLONG) ft.dwHighDateTime) << 32) + ft.dwLowDateTime) - start); +#else + /* use clock() as a fallback... this is somewhat annoying + because clock() may wrap around with a fairly short period */ + static clock_t start; + if (!start_inited) { + start_inited = 1; + start = clock(); + } + return (clock() - start) * 1.0 / CLOCKS_PER_SEC; +#endif +} + +/* number based on time for use as random seed */ +unsigned long nlopt_time_seed(void) +{ +#if defined(HAVE_GETTIMEOFDAY) + struct timeval tv; + gettimeofday(&tv, NULL); + return (tv.tv_sec ^ tv.tv_usec); +#elif defined(HAVE_TIME) + return time(NULL); +#elif defined(_WIN32) || defined(__WIN32__) + FILETIME ft; + GetSystemTimeAsFileTime(&ft); + return ft.dwHighDateTime ^ ft.dwLowDateTime; +#else + return clock(); +#endif +}