From 2dd33a7cbf438f3545794ff27415f5824541022a Mon Sep 17 00:00:00 2001 From: stevenj Date: Thu, 15 Jul 2010 12:20:51 -0400 Subject: [PATCH] allow BOBYQA and COBYLA to use unequal initial step sizes in different dimensions (internally rescaling the coordinates as needed); thanks to Tom Fiddaman for pointing out the problem when different coordinates have very different units darcs-hash:20100715162051-c8de0-5685b685d14f5042f824bbe5bff7f3aeed0fb808.gz --- api/optimize.c | 34 +++++++++++------ bobyqa/bobyqa.c | 84 +++++++++++++++++++++++++++++++--------- bobyqa/bobyqa.h | 7 ++-- cobyla/cobyla.c | 97 ++++++++++++++++++++++++++++------------------- cobyla/cobyla.h | 8 ++-- util/Makefile.am | 2 +- util/nlopt-util.h | 6 +++ util/rescale.c | 68 +++++++++++++++++++++++++++++++++ 8 files changed, 231 insertions(+), 75 deletions(-) create mode 100644 util/rescale.c diff --git a/api/optimize.c b/api/optimize.c index 6221b24..ef9402a 100644 --- a/api/optimize.c +++ b/api/optimize.c @@ -426,21 +426,27 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) 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: { @@ -460,11 +466,17 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) } 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: diff --git a/bobyqa/bobyqa.c b/bobyqa/bobyqa.c index 55f4f7c..9c3a8b7 100644 --- a/bobyqa/bobyqa.c +++ b/bobyqa/bobyqa.c @@ -9,6 +9,8 @@ #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)) @@ -3054,10 +3056,25 @@ L720: /**************************************************************************/ +#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; @@ -3069,15 +3086,41 @@ nlopt_result bobyqa(int n, int npt, double *x, 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, */ @@ -3120,13 +3163,13 @@ nlopt_result bobyqa(int n, int npt, double *x, /* 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 */ @@ -3152,9 +3195,9 @@ nlopt_result bobyqa(int n, int npt, double *x, 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 */ @@ -3170,8 +3213,8 @@ nlopt_result bobyqa(int n, int npt, double *x, 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; @@ -3208,11 +3251,18 @@ nlopt_result bobyqa(int n, int npt, double *x, /* 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_ */ diff --git a/bobyqa/bobyqa.h b/bobyqa/bobyqa.h index a5917cb..598753a 100644 --- a/bobyqa/bobyqa.h +++ b/bobyqa/bobyqa.h @@ -4,11 +4,10 @@ #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 */ diff --git a/cobyla/cobyla.c b/cobyla/cobyla.c index 990f389..71916f0 100644 --- a/cobyla/cobyla.c +++ b/cobyla/cobyla.c @@ -58,29 +58,34 @@ static char const rcsid[] = #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 @@ -91,6 +96,7 @@ static int func_wrap(int n, int m, double *x, double *f, double *con, 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; @@ -166,31 +172,48 @@ typedef int cobyla_function(int n, int m, double *x, double *f, double *con, * 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); @@ -204,24 +227,25 @@ nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data, } 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) { @@ -229,8 +253,12 @@ nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data, 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; } @@ -274,7 +302,7 @@ static double lcg_urand(uint32_t *seed, double a, double b) /**************************************************************************/ -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, @@ -285,7 +313,7 @@ static nlopt_result trstlp(int *n, int *m, double *a, double *b, double *rho, /* ------------------------------------------------------------------------ */ -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; @@ -365,13 +393,13 @@ nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, nlopt_ } /* 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"); @@ -397,7 +425,7 @@ nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, nlopt_ 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); @@ -413,7 +441,7 @@ nlopt_result cobyla(int n, int m, double *x, double *minf, double rhobeg, nlopt_ /* ------------------------------------------------------------------------- */ 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, @@ -447,16 +475,9 @@ static nlopt_result cobylb(int *n, int *m, int *mpp, 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; diff --git a/cobyla/cobyla.h b/cobyla/cobyla.h index 824c3bd..fd71ee1 100644 --- a/cobyla/cobyla.h +++ b/cobyla/cobyla.h @@ -46,14 +46,14 @@ extern "C" #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 } diff --git a/util/Makefile.am b/util/Makefile.am index 04ec0c9..ca4dff9 100644 --- a/util/Makefile.am +++ b/util/Makefile.am @@ -1,7 +1,7 @@ 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 diff --git a/util/nlopt-util.h b/util/nlopt-util.h index 484c3b7..9b4f51f 100644 --- a/util/nlopt-util.h +++ b/util/nlopt-util.h @@ -115,6 +115,12 @@ extern void nlopt_eval_constraint(double *result, double *grad, 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 */ diff --git a/util/rescale.c b/util/rescale.c new file mode 100644 index 0000000..48c8c3a --- /dev/null +++ b/util/rescale.c @@ -0,0 +1,68 @@ +/* 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 +#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; +} -- 2.30.2