From ed93948c785f42c084229d5b6501f026056fe667 Mon Sep 17 00:00:00 2001
From: stevenj <stevenj@alum.mit.edu>
Date: Mon, 18 Aug 2008 17:13:32 -0400
Subject: [PATCH] more verbose output from MMA; more paranoia required in
 feasibility check because feasibility constraints may be violated by rounding
 errors

darcs-hash:20080818211332-c8de0-e4e18ad0909671df4f63813f52c7f586ad93da86.gz
---
 api/nlopt.c |  2 +-
 mma/mma.c   | 26 +++++++++++++++++++++-----
 2 files changed, 22 insertions(+), 6 deletions(-)

diff --git a/api/nlopt.c b/api/nlopt.c
index e10448f..b370894 100644
--- a/api/nlopt.c
+++ b/api/nlopt.c
@@ -415,7 +415,7 @@ static nlopt_result nlopt_minimize_(
 	 case NLOPT_LD_MMA:
 	      return mma_minimize(n, f, f_data, m, fc, fc_data, fc_datum_size,
 				  lb, ub, x, minf, &stop,
-				  local_search_alg_deriv, 1e-8, 100000);
+				  local_search_alg_deriv, 1e-12, 100000);
 
 	 default:
 	      return NLOPT_INVALID_ARGS;
diff --git a/mma/mma.c b/mma/mma.c
index cf3a214..f1b72d1 100644
--- a/mma/mma.c
+++ b/mma/mma.c
@@ -21,6 +21,7 @@ int mma_verbose = 0; /* > 0 for verbose output */
 /* function for MMA's dual solution of the approximate problem */
 
 typedef struct {
+     int count; /* evaluation count, incremented each call */
      int n; /* must be set on input to dimension of x */
      const double *x, *lb, *ub, *sigma, *dfdx; /* arrays of length n */
      const double *dfcdx; /* m-by-n array of fc gradients */
@@ -46,6 +47,8 @@ static double dual_func(int m, const double *y, double *grad, void *d_)
      int i, j;
      double val;
 
+     d->count++;
+
      val = d->gval = fval;
      d->wval = 0;
      for (i = 0; i < m; ++i) val += y[i] * (gcval[i] = fcval[i]);
@@ -191,7 +194,7 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
 	       nlopt_result reti;
 
 	       /* solve dual problem */
-	       dd.rho = rho;
+	       dd.rho = rho; dd.count = 0;
 	       reti = nlopt_minimize(
 		    dual_alg, m, dual_func, &dd,
 		    dual_lb, dual_ub, y, &min_dual,
@@ -203,20 +206,33 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
 	       }
 
 	       dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
+	       if (mma_verbose) {
+		    printf("MMA dual converged in %d iterations to g=%g:\n",
+			   dd.count, dd.gval);
+		    for (i = 0; i < mma_verbose; ++i)
+			 printf("    MMA y[%d]=%g, gc[%d]=%g\n",
+				i, y[i], i, dd.gcval[i]);
+	       }
 
 	       fcur = f(n, xcur, dfdx_cur, f_data);
 	       stop->nevals++;
-	       inner_done = dd.gval >= fcur;
 	       feasible_cur = 1;
+	       inner_done = dd.gval >= fcur;
 	       for (i = 0; i < m; ++i) {
 		    fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n, 
 				      fc_data + fc_datum_size * i);
-		    feasible_cur = feasible_cur && (fcval_cur[i] <= 0);
+		    feasible_cur = feasible_cur && (fcval[i] <= 0);
 		    inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]);
 	       }
 
-	       if (fcur < *minf && (feasible_cur || !feasible)) {
-		    feasible = feasible_cur;
+	       /* once we have reached a feasible solution, the
+		  algorithm should never make the solution infeasible
+		  again, although the constraints may be violated
+		  slightly by rounding errors etc. so we must be a
+		  little careful about checking feasibility */
+	       if (feasible_cur) feasible = 1;
+
+	       if (fcur < *minf) {
 		    dd.fval = *minf = fcur;
 		    memcpy(fcval, fcval_cur, sizeof(double)*m);
 		    memcpy(x, xcur, sizeof(double)*n);
-- 
2.30.2