From efa005713e44175c9b0ca20b5012b35dd6555ecf Mon Sep 17 00:00:00 2001 From: Joshua Nathaniel Pritikin Date: Fri, 8 May 2015 11:10:00 -0400 Subject: [PATCH] Separately track the best minor and major estimate This change ensures that we always move to the best point found during the line search even when the best point is not the last point. We also avoid 1 fit and gradient evaluation by moving the major iteration convergence check above the switch statement. --- slsqp/slsqp.c | 75 ++++++++++++++++++++++++++++++++++----------------- 1 file changed, 51 insertions(+), 24 deletions(-) diff --git a/slsqp/slsqp.c b/slsqp/slsqp.c index c58c3d4e..e8b29a22 100644 --- a/slsqp/slsqp.c +++ b/slsqp/slsqp.c @@ -2464,6 +2464,14 @@ static void estimate_return(struct estimate *est, double *minf, double *x) memcpy(x, est->par, sizeof(double)*est->n); } +static void estimate_copy(struct estimate *to, struct estimate *from) +{ + to->feasible = from->feasible; + to->infeasibility = from->infeasibility; + memcpy(to->par, from->par, sizeof(double)*to->n); + to->fval = from->fval; +} + nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data, unsigned m, nlopt_constraint *fc, unsigned p, nlopt_constraint *h, @@ -2484,7 +2492,8 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data, nlopt_result ret = NLOPT_SUCCESS; unsigned max_cdim; int want_grad = 1; - struct estimate cur, minor; + struct estimate cur, minor, major; + int foundBetterSpot = 0; max_cdim = MAX2(nlopt_max_constraint_dim(m, fc), nlopt_max_constraint_dim(p, h)); @@ -2493,7 +2502,7 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data, #define U(n) ((unsigned) (n)) work = (double *) malloc(sizeof(double) * (U(mpi1) * (n + 1) + U(mpi) - + n+1 + n + n + max_cdim*n + + n+1 + n + n + n + max_cdim*n + U(len_w)) + sizeof(int) * U(len_jw)); if (!work) return NLOPT_OUT_OF_MEMORY; @@ -2502,7 +2511,8 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data, grad = c + mpi; estimate_init(&cur, n, grad + n+1, theSpot); estimate_init(&minor, n, cur.par + n, theSpot); - cgradtmp = minor.par + n; + estimate_init(&major, n, minor.par + n, theSpot); + cgradtmp = major.par + n; w = cgradtmp + max_cdim*n; jw = (int *) (w + len_w); @@ -2516,6 +2526,29 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data, w, &len_w, jw, &len_jw, &state); + /* note: mode == -1 corresponds to the completion of a line search, + and is the only time we should check convergence (as in original slsqp code) */ + if (mode == -1 && !nlopt_isinf(minor.fval)) { + estimate_copy(&cur, &minor); + //printf("best minor %f %f\n", minor.fval, minor.infeasibility); + + if (minor.feasible || !foundBetterSpot) { + //printf("check major %f %f\n", minor.fval, minor.infeasibility); + if (!nlopt_isinf(major.fval)) { + if (nlopt_stop_ftol(stop, minor.fval, major.fval)) + ret = NLOPT_FTOL_REACHED; + else if (nlopt_stop_x(stop, minor.par, major.par)) + ret = NLOPT_XTOL_REACHED; + if (ret != NLOPT_SUCCESS) goto done; + } + estimate_copy(&major, &minor); + } + + /* prep for new minor iteration */ + minor.feasible = 0; + foundBetterSpot = 0; + } + switch (mode) { case -1: /* objective & gradient evaluation */ if (prev_mode == -2 && !want_grad) break; /* just evaluated this point */ @@ -2618,39 +2651,33 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data, } prev_mode = mode; + /* printf("check minor: fval %f < %f %d infe %f < %f %d\n", */ + /* cur.fval, minor.fval, cur.fval < minor.fval, */ + /* cur.infeasibility, minor.infeasibility, cur.infeasibility < minor.infeasibility); */ + /* update best point so far */ - if (nlopt_isfinite(cur.fval) && ((cur.fval < *minf && (cur.feasible || !minor.feasible)) - || (!minor.feasible && cur.infeasibility < minor.infeasibility))) { - *minf = cur.fval; - minor.feasible = cur.feasible; - minor.infeasibility = cur.infeasibility; - } + if (nlopt_isfinite(cur.fval) && + ((cur.fval < minor.fval && (cur.feasible || !minor.feasible)) + || (!minor.feasible && cur.infeasibility < minor.infeasibility))) { - /* note: mode == -1 corresponds to the completion of a line search, - and is the only time we should check convergence (as in original slsqp code) */ - if (mode == -1) { - if (!nlopt_isinf(minor.fval)) { - if (nlopt_stop_ftol(stop, cur.fval, minor.fval)) - ret = NLOPT_FTOL_REACHED; - else if (nlopt_stop_x(stop, cur.par, minor.par)) - ret = NLOPT_XTOL_REACHED; - } - minor.fval = cur.fval; - memcpy(minor.par, cur.par, sizeof(double)*n); + //printf("best eval so far %f %f\n", cur.fval, cur.infeasibility); + estimate_copy(&minor, &cur); + foundBetterSpot = 1; } /* do some additional termination tests */ if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED; else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED; - else if (minor.feasible && *minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED; + else if (major.feasible && major.fval < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED; } while (ret == NLOPT_SUCCESS); done: - if (!nlopt_isinf(minor.fval)) { + if (!nlopt_isinf(major.fval)) { + estimate_return(&major, minf, theSpot); + } else if (!nlopt_isinf(minor.fval)) { estimate_return(&minor, minf, theSpot); } else if (!nlopt_isinf(cur.fval)) { - /* didn't find any feasible points, just return last point evaluated */ - estimate_return(&cur, minf, theSpot); + estimate_return(&cur, minf, theSpot); } free(work);