nlopt_set_maxeval(dual_opt, LO(maxeval, 100000));
#undef LO
- ret = mma_minimize(ni, f, f_data, (int) (opt->m), opt->fc,
+ ret = mma_minimize(n, f, f_data, opt->m, opt->fc,
lb, ub, x, minf, &stop, dual_opt);
nlopt_destroy(dual_opt);
return ret;
}
case NLOPT_LN_COBYLA: {
- double step;
- if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
- return NLOPT_OUT_OF_MEMORY;
- return cobyla_minimize(ni, f, f_data,
+ nlopt_result ret;
+ int freedx = 0;
+ if (!opt->dx) {
+ freedx = 1;
+ if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
+ return NLOPT_OUT_OF_MEMORY;
+ }
+ return cobyla_minimize(n, f, f_data,
opt->m, opt->fc,
opt->p, opt->h,
lb, ub, x, minf, &stop,
- step);
+ opt->dx);
+ if (freedx) { free(opt->dx); opt->dx = NULL; }
+ return ret;
}
case NLOPT_LN_NEWUOA: {
}
case NLOPT_LN_BOBYQA: {
- double step;
- if (initial_step(opt, x, &step) != NLOPT_SUCCESS)
- return NLOPT_OUT_OF_MEMORY;
- return bobyqa(ni, 2*n+1, x, lb, ub, step,
- &stop, minf, f_noderiv, opt);
+ nlopt_result ret;
+ int freedx = 0;
+ if (!opt->dx) {
+ freedx = 1;
+ if (nlopt_set_default_initial_step(opt, x) != NLOPT_SUCCESS)
+ return NLOPT_OUT_OF_MEMORY;
+ }
+ ret = bobyqa(ni, 2*n+1, x, lb, ub, opt->dx,
+ &stop, minf, opt->f, opt->f_data);
+ if (freedx) { free(opt->dx); opt->dx = NULL; }
+ return ret;
}
case NLOPT_LN_NELDERMEAD:
#include "bobyqa.h"
+typedef double (*bobyqa_func)(int n, const double *x, void *func_data);
+
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))
#define iabs(x) ((x) < 0 ? -(x) : (x))
/**************************************************************************/
+#define U(n) ((unsigned) (n))
+
+typedef struct {
+ double *s, *xs;
+ nlopt_func f; void *f_data;
+} rescale_fun_data;
+
+static double rescale_fun(int n, const double *x, void *d_)
+{
+ rescale_fun_data *d = (rescale_fun_data*) d_;
+ nlopt_unscale(U(n), d->s, x, d->xs);
+ return d->f(U(n), d->xs, NULL, d->f_data);
+}
+
nlopt_result bobyqa(int n, int npt, double *x,
- const double *xl, const double *xu, double rhobeg,
+ const double *xl, const double *xu,
+ const double *dx,
nlopt_stopping *stop, double *minf,
- bobyqa_func calfun, void *calfun_data)
+ nlopt_func f, void *f_data)
{
/* System generated locals */
int i__1;
double temp, zero;
int ibmat, izmat;
- double rhoend;
- double *w;
+ double rhobeg, rhoend;
+ double *w0 = NULL, *w;
nlopt_result ret;
-
-/* SGJ, 2009: compute rhoend from NLopt stop info */
+ double *s = NULL, *sxl = NULL, *sxu = NULL, *xs = NULL;
+ rescale_fun_data calfun_data;
+
+ /* SGJ 2010: rescale parameters to make the initial step sizes dx
+ equal in all directions */
+ s = nlopt_compute_rescaling(U(n), dx);
+ if (!s) return NLOPT_OUT_OF_MEMORY;
+
+ /* this statement must go before goto done, so that --x occurs */
+ nlopt_rescale(U(n), s, x, x); --x;
+
+ xs = (double *) malloc(sizeof(double) * (U(n)));
+
+ sxl = nlopt_new_rescaled(U(n), s, xl);
+ if (!sxl) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
+ xl = sxl;
+ sxu = nlopt_new_rescaled(U(n), s, xu);
+ if (!sxu) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
+ xu = sxu;
+
+ rhobeg = dx[0] / s[0]; /* equals all other dx[i] after rescaling */
+
+ calfun_data.s = s;
+ calfun_data.xs = xs;
+ calfun_data.f = f;
+ calfun_data.f_data = f_data;
+
+ /* SGJ, 2009: compute rhoend from NLopt stop info */
rhoend = stop->xtol_rel * (rhobeg);
for (j = 0; j < n; ++j)
- if (rhoend < stop->xtol_abs[j])
- rhoend = stop->xtol_abs[j];
+ if (rhoend < stop->xtol_abs[j] / s[j])
+ rhoend = stop->xtol_abs[j] / s[j];
/* This subroutine seeks the least value of a function of many variables, */
/* Parameter adjustments */
--xu;
--xl;
- --x;
/* Function Body */
np = n + 1;
if (npt < n + 2 || npt > (n + 2) * np / 2) {
/* Return from BOBYQA because NPT is not in the required interval */
- return NLOPT_INVALID_ARGS;
+ ret = NLOPT_INVALID_ARGS;
+ goto done;
}
/* Partition the working space array, so that different parts of it can */
ivl = id + n;
iw = ivl + ndim;
- w = (double *) malloc(sizeof(double) * ((npt+5)*(npt+n)+3*n*(n+5)/2));
- if (!w) return NLOPT_OUT_OF_MEMORY;
- --w;
+ w0 = (double *) malloc(sizeof(double) * U((npt+5)*(npt+n)+3*n*(n+5)/2));
+ if (!w0) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
+ w = w0 - 1;
/* Return if there is insufficient space between the bounds. Modify the */
/* initial X if necessary in order to avoid conflicts between the bounds */
if (temp < rhobeg + rhobeg) {
/* Return from BOBYQA because one of the differences
XU(I)-XL(I)s is less than 2*RHOBEG. */
- free(w+1);
- return NLOPT_INVALID_ARGS;
+ ret = NLOPT_INVALID_ARGS;
+ goto done;
}
jsl = isl + j - 1;
jsu = jsl + n;
/* Make the call of BOBYQB. */
ret = bobyqb_(&n, &npt, &x[1], &xl[1], &xu[1], &rhobeg, &rhoend,
- stop, calfun, calfun_data, minf,
+ stop, rescale_fun, &calfun_data, minf,
&w[ixb], &w[ixp], &w[ifv], &w[ixo], &w[igo], &w[ihq], &w[ipq],
&w[ibmat], &w[izmat], &ndim, &w[isl], &w[isu], &w[ixn], &w[ixa],
&w[id], &w[ivl], &w[iw]);
- free(w+1);
+
+done:
+ free(w0);
+ free(sxl);
+ free(sxu);
+ free(xs);
+ ++x; nlopt_unscale(U(n), s, x, x);
+ free(s);
return ret;
} /* bobyqa_ */
#include "nlopt-util.h"
#include "nlopt.h"
-typedef double (*bobyqa_func)(int n, const double *x, void *func_data);
-
extern nlopt_result bobyqa(int n, int npt, double *x,
const double *lb, const double *ub,
- double rhobeg, nlopt_stopping *stop, double *minf,
- bobyqa_func calfun, void *calfun_data);
+ const double *dx,
+ nlopt_stopping *stop, double *minf,
+ nlopt_func f, void *f_data);
#endif /* BOBYQA_H */
#define max(a,b) ((a) >= (b) ? (a) : (b))
#define abs(x) ((x) >= 0 ? (x) : -(x))
+#define U(n) ((unsigned) (n))
+
/**************************************************************************/
/* SGJ, 2008: NLopt-style interface function: */
typedef struct {
nlopt_func f;
void *f_data;
- int m_orig;
+ unsigned m_orig;
nlopt_constraint *fc;
- int p;
+ unsigned p;
nlopt_constraint *h;
double *xtmp;
- const double *lb, *ub;
- double *con_tol;
+ double *lb, *ub;
+ double *con_tol, *scale;
nlopt_stopping *stop;
} func_wrap_state;
-static int func_wrap(int n, int m, double *x, double *f, double *con,
+static int func_wrap(int ni, int mi, double *x, double *f, double *con,
func_wrap_state *s)
{
- int i, j, k;
+ unsigned n = U(ni);
+ unsigned i, j, k;
double *xtmp = s->xtmp;
const double *lb = s->lb, *ub = s->ub;
+ (void) mi; /* unused */
+
/* in nlopt, we guarante that the function is never evaluated outside
the lb and ub bounds, so we need force this with xtmp ... note
that this leads to discontinuity in the first derivative, which
else if (x[j] > ub[j]) xtmp[j] = ub[j];
else xtmp[j] = x[j];
}
+ nlopt_unscale(n, s->scale, xtmp, xtmp);
*f = s->f(n, xtmp, NULL, s->f_data);
if (nlopt_stop_forced(s->stop)) return 1;
* The cobyla function returns the usual nlopt_result codes.
*
*/
-extern nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, nlopt_stopping *stop, const double *lb, const double *ub,
+extern nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, double rhoend, nlopt_stopping *stop, const double *lb, const double *ub,
int message, cobyla_function *calcfc, func_wrap_state *state);
-nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data,
- int m, nlopt_constraint *fc,
- int p, nlopt_constraint *h,
+nlopt_result cobyla_minimize(unsigned n, nlopt_func f, void *f_data,
+ unsigned m, nlopt_constraint *fc,
+ unsigned p, nlopt_constraint *h,
const double *lb, const double *ub, /* bounds */
double *x, /* in: initial guess, out: minimizer */
double *minf,
nlopt_stopping *stop,
- double rhobegin)
+ const double *dx)
{
- int i, j;
+ unsigned i, j;
func_wrap_state s;
nlopt_result ret;
+ double rhobeg, rhoend;
s.f = f; s.f_data = f_data;
s.m_orig = m;
s.fc = fc;
s.p = p;
s.h = h;
- s.lb = lb; s.ub = ub;
s.stop = stop;
+ s.lb = s.ub = s.xtmp = s.con_tol = s.scale = NULL;
+
+ s.scale = nlopt_compute_rescaling(n, dx);
+ if (!s.scale) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
+
+ s.lb = nlopt_new_rescaled(n, s.scale, lb);
+ if (!s.lb) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
+ s.ub = nlopt_new_rescaled(n, s.scale, ub);
+ if (!s.ub) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
+
s.xtmp = (double *) malloc(sizeof(double) * n);
- if (!s.xtmp) return NLOPT_OUT_OF_MEMORY;
+ if (!s.xtmp) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
+
+ /* SGJ, 2008: compute rhoend from NLopt stop info */
+ rhobeg = dx[0] / s.scale[0];
+ rhoend = stop->xtol_rel * (rhobeg);
+ for (j = 0; j < n; ++j)
+ if (rhoend < stop->xtol_abs[j] / s.scale[j])
+ rhoend = stop->xtol_abs[j] / s.scale[j];
/* each equality constraint gives two inequality constraints */
m = nlopt_count_constraints(m, fc) + 2 * nlopt_count_constraints(p, h);
}
s.con_tol = (double *) malloc(sizeof(double) * m);
- if (m && !s.con_tol) {
- free(s.xtmp);
- return NLOPT_OUT_OF_MEMORY;
- }
+ if (m && !s.con_tol) { ret = NLOPT_OUT_OF_MEMORY; goto done; }
+
for (j = 0; j < m; ++j) s.con_tol[j] = 0;
for (j = i = 0; i < s.m_orig; ++i) {
- int j0 = j, jnext = j + fc[i].m;
- for (; j < jnext; ++j) s.con_tol[j] = fc[i].tol[j - j0];
+ unsigned ji = j, jnext = j + fc[i].m;
+ for (; j < jnext; ++j) s.con_tol[j] = fc[i].tol[j - ji];
}
for (i = 0; i < s.p; ++i) {
- int j0 = j, jnext = j + h[i].m;
- for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - j0];
- j0 = j; jnext = j + h[i].m;
- for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - j0];
+ unsigned ji = j, jnext = j + h[i].m;
+ for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - ji];
+ ji = j; jnext = j + h[i].m;
+ for (; j < jnext; ++j) s.con_tol[j] = h[i].tol[j - ji];
}
- ret = cobyla(n, m, x, minf, rhobegin, stop, lb, ub, COBYLA_MSG_NONE,
+ nlopt_rescale(n, s.scale, x, x);
+ ret = cobyla((int) n, (int) m, x, minf, rhobeg, rhoend,
+ stop, s.lb, s.ub, COBYLA_MSG_NONE,
func_wrap, &s);
+ nlopt_unscale(n, s.scale, x, x);
/* make sure e.g. rounding errors didn't push us slightly out of bounds */
for (j = 0; j < n; ++j) {
if (x[j] > ub[j]) x[j] = ub[j];
}
+done:
free(s.con_tol);
free(s.xtmp);
+ free(s.ub);
+ free(s.lb);
+ free(s.scale);
return ret;
}
/**************************************************************************/
-static nlopt_result cobylb(int *n, int *m, int *mpp, double *x, double *minf, double *rhobeg,
+static nlopt_result cobylb(int *n, int *m, int *mpp, double *x, double *minf, double *rhobeg, double rhoend,
nlopt_stopping *stop, const double *lb, const double *ub, int *iprint, double *con, double *sim,
double *simi, double *datmat, double *a, double *vsig, double *veta,
double *sigbar, double *dx, double *w, int *iact, cobyla_function *calcfc,
/* ------------------------------------------------------------------------ */
-nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, nlopt_stopping *stop, const double *lb, const double *ub, int iprint,
+nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, double rhoend, nlopt_stopping *stop, const double *lb, const double *ub, int iprint,
cobyla_function *calcfc, func_wrap_state *state)
{
int icon, isim, isigb, idatm, iveta, isimi, ivsig, iwork, ia, idx, mpp;
}
/* workspace allocation */
- w = (double*) malloc((n*(3*n+2*m+11)+4*m+6)*sizeof(*w));
+ w = (double*) malloc(U(n*(3*n+2*m+11)+4*m+6)*sizeof(*w));
if (w == NULL)
{
if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
return NLOPT_OUT_OF_MEMORY;
}
- iact = (int*)malloc((m+1)*sizeof(*iact));
+ iact = (int*)malloc(U(m+1)*sizeof(*iact));
if (iact == NULL)
{
if (iprint>=1) fprintf(stderr, "cobyla: memory allocation error.\n");
isigb = iveta + n;
idx = isigb + n;
iwork = idx + n;
- rc = cobylb(&n, &m, &mpp, &x[1], minf, &rhobeg, stop, &lb[1], &ub[1], &iprint,
+ rc = cobylb(&n, &m, &mpp, &x[1], minf, &rhobeg, rhoend, stop, &lb[1], &ub[1], &iprint,
&w[icon], &w[isim], &w[isimi], &w[idatm], &w[ia], &w[ivsig], &w[iveta],
&w[isigb], &w[idx], &w[iwork], &iact[1], calcfc, state);
/* ------------------------------------------------------------------------- */
static nlopt_result cobylb(int *n, int *m, int *mpp,
- double *x, double *minf, double *rhobeg,
+ double *x, double *minf, double *rhobeg, double rhoend,
nlopt_stopping *stop, const double *lb, const double *ub,
int *iprint, double *con, double *sim, double *simi,
double *datmat, double *a, double *vsig, double *veta,
int mp, np, iz, ibrnch;
int nbest, ifull, iptem, jdrop;
nlopt_result rc = NLOPT_SUCCESS;
- double rhoend;
uint32_t seed = (uint32_t) (*n + *m); /* arbitrary deterministic LCG seed */
int feasible;
- /* SGJ, 2008: compute rhoend from NLopt stop info */
- rhoend = stop->xtol_rel * (*rhobeg);
- for (j = 0; j < *n; ++j)
- if (rhoend < stop->xtol_abs[j])
- rhoend = stop->xtol_abs[j];
-
/* SGJ, 2008: added code to keep track of minimum feasible function val */
*minf = HUGE_VAL;
#endif /* __cplusplus */
/* NLopt-style interface function */
-nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data,
- int m, nlopt_constraint *fc,
- int p, nlopt_constraint *h,
+nlopt_result cobyla_minimize(unsigned n, nlopt_func f, void *f_data,
+ unsigned m, nlopt_constraint *fc,
+ unsigned p, nlopt_constraint *h,
const double *lb, const double *ub, /* bounds */
double *x, /* in: initial guess, out: minimizer */
double *minf,
nlopt_stopping *stop,
- double rhobegin);
+ const double *dx);
#ifdef __cplusplus
}
AM_CPPFLAGS = -I$(top_srcdir)/api
noinst_LTLIBRARIES = libutil.la
-libutil_la_SOURCES = mt19937ar.c sobolseq.c soboldata.h timer.c stop.c nlopt-util.h redblack.c redblack.h qsort_r.c
+libutil_la_SOURCES = mt19937ar.c sobolseq.c soboldata.h timer.c stop.c nlopt-util.h redblack.c redblack.h qsort_r.c rescale.c
noinst_PROGRAMS = redblack_test
redblack_test_SOURCES = redblack_test.c
const nlopt_constraint *c,
unsigned n, const double *x);
+/* rescale.c: */
+double *nlopt_compute_rescaling(unsigned n, const double *dx);
+double *nlopt_new_rescaled(unsigned n, const double *s, const double *x);
+void nlopt_rescale(unsigned n, const double *s, const double *x, double *xs);
+void nlopt_unscale(unsigned n, const double *s, const double *x, double *xs);
+
#ifdef __cplusplus
} /* extern "C" */
#endif /* __cplusplus */
--- /dev/null
+/* Copyright (c) 2007-2010 Massachusetts Institute of Technology
+ *
+ * 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.
+ */
+
+#include <stdlib.h>
+#include "nlopt-util.h"
+
+/* Return a new array of length n (> 0) that gives a rescaling factor
+ for each dimension, or NULL if out of memory, with dx being the
+ array of nonzero initial steps in each dimension. */
+double *nlopt_compute_rescaling(unsigned n, const double *dx)
+{
+ double *s = (double *) malloc(sizeof(double) * n);
+ unsigned i;
+
+ if (!s) return NULL;
+ for (i = 0; i < n; ++i) s[i] = 1.0; /* default: no rescaling */
+ if (n == 1) return s;
+
+ for (i = 1; i < n && dx[i] == dx[i-1]; ++i) ;
+ if (i < n) { /* unequal initial steps, rescale to make equal to dx[0] */
+ for (i = 1; i < n; ++i)
+ s[i] = dx[i] / dx[0];
+ }
+ return s;
+}
+
+void nlopt_rescale(unsigned n, const double *s, const double *x, double *xs)
+{
+ unsigned i;
+ if (!s) { for (i = 0; i < n;++i) xs[i] = x[i]; }
+ else { for (i = 0; i < n;++i) xs[i] = x[i] / s[i]; }
+}
+
+void nlopt_unscale(unsigned n, const double *s, const double *x, double *xs)
+{
+ unsigned i;
+ if (!s) { for (i = 0; i < n;++i) xs[i] = x[i]; }
+ else { for (i = 0; i < n;++i) xs[i] = x[i] * s[i]; }
+}
+
+/* return a new array of length n equal to the original array
+ x divided by the scale factors s, or NULL on a memory error */
+double *nlopt_new_rescaled(unsigned n, const double *s, const double *x)
+{
+ double *xs = (double *) malloc(sizeof(double) * n);
+ if (!xs) return NULL;
+ nlopt_rescale(n, s, x, xs);
+ return xs;
+}
~ianmdlvl / nlopt.git / commitdiff