in functions with constraints, make sure forced_stop stops immediately, before evalua...

darcs-hash:20100708212353-c8de0-5060b668f7e91e93ca182f3c2a8ee7618b4d5032.gz

diff --git a/NEWS b/NEWS
index 5014cc1b7f3d38be332c18a808fc25c3531d9b5f..e8bd3b0abb8f7dc089861d239fdb5932e6f17540 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -4,6 +4,10 @@ NLopt 2.1.2 (8 July 2010)
   now pass a 2-dimensional array for the gradient argument, rather
   than a flattened 1d array.
 
+* Improved handling of exceptions and forced stops for constrained
+  optimization, making sure that no constraints are evaluated after
+  the stop.
+
 * Return an NLOPT_INVALID_ARGS error if more than n equality constraints
   are added in an n-dimensional problem.
 

diff --git a/auglag/auglag.c b/auglag/auglag.c
index 1b59aa83aaac9a6db883e9510c743f9a8b69312d..f43934e6f9037843f219e442a378f0fc0b509646 100644 (file)
--- a/auglag/auglag.c
+++ b/auglag/auglag.c
@@ -34,9 +34,12 @@ static double auglag(unsigned n, const double *x, double *grad, void *data)
      unsigned j, k;
 
      L = d->f(n, x, grad, d->f_data);
+     d->stop->nevals++;
+     if (nlopt_stop_forced(d->stop)) return L;
 
      for (ii = i = 0; i < d->p; ++i) {
          nlopt_eval_constraint(restmp, gradtmp, d->h + i, n, x);
+         if (nlopt_stop_forced(d->stop)) return L;
          for (k = 0; k < d->h[i].m; ++k) {
               double h = restmp[k] + lambda[ii++] / rho;
               L += 0.5 * rho * h*h;
@@ -47,6 +50,7 @@ static double auglag(unsigned n, const double *x, double *grad, void *data)
 
      for (ii = i = 0; i < d->m; ++i) {
          nlopt_eval_constraint(restmp, gradtmp, d->fc + i, n, x);
+         if (nlopt_stop_forced(d->stop)) return L;
          for (k = 0; k < d->fc[i].m; ++k) {
               double fc = restmp[k] + mu[ii++] / rho;
               if (fc > 0) {
@@ -56,8 +60,6 @@ static double auglag(unsigned n, const double *x, double *grad, void *data)
               }
          }
      }
-     
-     d->stop->nevals++;
 
      return L;
 }
@@ -142,10 +144,14 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data,
          double con2 = 0;
          d.stop->nevals++;
          fcur = f(n, xcur, NULL, f_data);
+         if (nlopt_stop_forced(stop)) {
+              ret = NLOPT_FORCED_STOP; goto done; }
          penalty = 0;
          feasible = 1;
          for (i = 0; i < d.p; ++i) {
               nlopt_eval_constraint(d.restmp, NULL, d.h + i, n, xcur);
+              if (nlopt_stop_forced(stop)) {
+                   ret = NLOPT_FORCED_STOP; goto done; }
               for (k = 0; k < d.h[i].m; ++k) {
                    double hi = d.restmp[k];
                    penalty += fabs(hi);
@@ -155,6 +161,8 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data,
          }
          for (i = 0; i < d.m; ++i) {
               nlopt_eval_constraint(d.restmp, NULL, d.fc + i, n, xcur);
+              if (nlopt_stop_forced(stop)) {
+                   ret = NLOPT_FORCED_STOP; goto done; }
               for (k = 0; k < d.fc[i].m; ++k) {
                    double fci = d.restmp[k];
                    penalty += fci > 0 ? fci : 0;
@@ -191,6 +199,8 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data,
          
          d.stop->nevals++;
          fcur = f(n, xcur, NULL, f_data);
+         if (nlopt_stop_forced(stop)) {
+              ret = NLOPT_FORCED_STOP; goto done; }
          if (auglag_verbose)
               printf("auglag: fcur = %g\n", fcur);
          
@@ -199,6 +209,8 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data,
          feasible = 1;
          for (i = ii = 0; i < d.p; ++i) {
               nlopt_eval_constraint(d.restmp, NULL, d.h + i, n, xcur);
+              if (nlopt_stop_forced(stop)) {
+                   ret = NLOPT_FORCED_STOP; goto done; }
               for (k = 0; k < d.h[i].m; ++k) {
                    double hi = d.restmp[k];
                    double newlam = d.lambda[ii] + d.rho * hi;
@@ -210,6 +222,8 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data,
          }
          for (i = ii = 0; i < d.m; ++i) {
               nlopt_eval_constraint(d.restmp, NULL, d.fc + i, n, xcur);
+              if (nlopt_stop_forced(stop)) {
+                   ret = NLOPT_FORCED_STOP; goto done; }
               for (k = 0; k < d.fc[i].m; ++k) {
                    double fci = d.restmp[k];
                    double newmu = d.mu[ii] + d.rho * fci;
@@ -269,6 +283,7 @@ nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data,
          if (ICM == 0) return NLOPT_FTOL_REACHED;
      } while (1);
 
+done:
      free(xcur);
      return ret;
 }

diff --git a/cobyla/cobyla.c b/cobyla/cobyla.c
index ea981cd6d6bb730aea36e1d43fb9a7e886dcd04c..990f38934038d11b29352ca788db95126497231b 100644 (file)
--- a/cobyla/cobyla.c
+++ b/cobyla/cobyla.c
@@ -71,6 +71,7 @@ typedef struct {
      double *xtmp;
      const double *lb, *ub;
      double *con_tol;
+     nlopt_stopping *stop;
 } func_wrap_state;
 
 static int func_wrap(int n, int m, double *x, double *f, double *con,
@@ -92,15 +93,18 @@ static int func_wrap(int n, int m, double *x, double *f, double *con,
      }
 
      *f = s->f(n, xtmp, NULL, s->f_data);
+     if (nlopt_stop_forced(s->stop)) return 1;
      i = 0;
      for (j = 0; j < s->m_orig; ++j) {
          nlopt_eval_constraint(con + i, NULL, s->fc+j, n, xtmp);
+         if (nlopt_stop_forced(s->stop)) return 1;
          for (k = 0; k < s->fc[j].m; ++k)
               con[i + k] = -con[i + k];
          i += s->fc[j].m;
      }
      for (j = 0; j < s->p; ++j) {
          nlopt_eval_constraint(con + i, NULL, s->h+j, n, xtmp);
+         if (nlopt_stop_forced(s->stop)) return 1;
          for (k = 0; k < s->h[j].m; ++k)
               con[(i + s->h[j].m) + k] = -con[i + k];
          i += 2 * s->h[j].m;
@@ -111,7 +115,6 @@ static int func_wrap(int n, int m, double *x, double *f, double *con,
          if (!nlopt_isinf(ub[j]))
               con[i++] = ub[j] - x[j];
      }
-     if (m != i) return 1; /* ... bug?? */
      return 0;
 }
 
@@ -185,6 +188,7 @@ nlopt_result cobyla_minimize(int n, nlopt_func f, void *f_data,
      s.p = p;
      s.h = h;
      s.lb = lb; s.ub = ub;
+     s.stop = stop;
      s.xtmp = (double *) malloc(sizeof(double) * n);
      if (!s.xtmp) return NLOPT_OUT_OF_MEMORY;
 
@@ -548,7 +552,7 @@ L40:
     if (*iprint >= 1) {
       fprintf(stderr, "cobyla: user requested end of minimization.\n");
     }
-    rc = NLOPT_FAILURE;
+    rc = NLOPT_FORCED_STOP;
     goto L600;
   }
 

diff --git a/isres/isres.c b/isres/isres.c
index 5b4fec3a50145550b24d5b5e26a8e21c1affb2d2..02b6f3c211ca76c4828c4fb83ae64059d4175b43 100644 (file)
--- a/isres/isres.c
+++ b/isres/isres.c
@@ -131,10 +131,14 @@ nlopt_result isres_minimize(int n, nlopt_func f, void *f_data,
               double gpenalty;
               stop->nevals++;
               fval[k] = f(n, xs + k*n, NULL, f_data);
+              if (nlopt_stop_forced(stop)) { 
+                   ret = NLOPT_FORCED_STOP; goto done; }
               penalty[k] = 0;
               for (c = 0; c < m; ++c) { /* inequality constraints */
                    nlopt_eval_constraint(results, NULL,
                                          fc + c, n, xs + k*n);
+                   if (nlopt_stop_forced(stop)) { 
+                        ret = NLOPT_FORCED_STOP; goto done; }
                    for (ires = 0; ires < fc[c].m; ++ires) {
                         double gval = results[ires];
                         if (gval > fc[c].tol[ires]) feasible = 0;
@@ -146,6 +150,8 @@ nlopt_result isres_minimize(int n, nlopt_func f, void *f_data,
               for (c = m; c < mp; ++c) { /* equality constraints */
                    nlopt_eval_constraint(results, NULL,
                                          h + (c-m), n, xs + k*n);
+                   if (nlopt_stop_forced(stop)) { 
+                        ret = NLOPT_FORCED_STOP; goto done; }
                    for (ires = 0; ires < h[c-m].m; ++ires) {
                         double hval = results[ires];
                         if (fabs(hval) > h[c-m].tol[ires]) feasible = 0;

diff --git a/mma/mma.c b/mma/mma.c
index f9fa5f2d4b47f19f5379506967d5550a204fe9dc..f3db3846f53da32b44ca446d9f654f003bdb9788 100644 (file)
--- a/mma/mma.c
+++ b/mma/mma.c
@@ -210,12 +210,14 @@ nlopt_result mma_minimize(unsigned n, nlopt_func f, void *f_data,
      dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
      stop->nevals++;
      memcpy(xcur, x, sizeof(double) * n);
+     if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
 
      feasible = 1; infeasibility = 0;
      for (i = ifc = 0; ifc < mfc; ++ifc) {
          nlopt_eval_constraint(fcval + i, dfcdx + i*n,
                                fc + ifc, n, x);
          i += fc[ifc].m;
+         if (nlopt_stop_forced(stop)) { ret = NLOPT_FORCED_STOP; goto done; }
      }
      for (i = 0; i < m; ++i) {
          feasible = feasible && (fcval[i] <= 0 || isnan(fcval[i]));
@@ -285,6 +287,8 @@ nlopt_result mma_minimize(unsigned n, nlopt_func f, void *f_data,
 
               fcur = f(n, xcur, dfdx_cur, f_data);
               stop->nevals++;
+              if (nlopt_stop_forced(stop)) { 
+                   ret = NLOPT_FORCED_STOP; goto done; }
               feasible_cur = 1; infeasibility_cur = 0;
               new_infeasible_constraint = 0;
               inner_done = dd.gval >= fcur;
@@ -292,6 +296,8 @@ nlopt_result mma_minimize(unsigned n, nlopt_func f, void *f_data,
                    nlopt_eval_constraint(fcval_cur + i, dfcdx_cur + i*n,
                                          fc + ifc, n, xcur);
                    i += fc[ifc].m;
+                   if (nlopt_stop_forced(stop)) { 
+                        ret = NLOPT_FORCED_STOP; goto done; }
               }
               for (i = ifc = 0; ifc < mfc; ++ifc) {
                    unsigned i0 = i, inext = i + fc[ifc].m;