add generic algorithm parameters API

doc/docs/NLopt_Algorithms.md

@@ -311,6 +311,13 @@ At each point **x**, MMA forms a local approximation using the gradient of *f* a
 
 I also implemented another CCSA algorithm from the same paper, `NLOPT_LD_CCSAQ`: instead of constructing local MMA approximations, it constructs simple quadratic approximations (or rather, affine approximations plus a quadratic penalty term to stay conservative). This is the ccsa_quadratic code. It seems to have similar convergence rates to MMA for most problems, which is not surprising as they are both essentially similar. However, for the quadratic variant I implemented the possibility of [preconditioning](NLopt_Reference.md#preconditioning-with-approximate-hessians): including a user-supplied Hessian approximation in the local model. It is easy to incorporate this into the proof in Svanberg's paper, and to show that global convergence is still guaranteed as long as the user's "Hessian" is positive semidefinite, and it practice it can greatly improve convergence if the preconditioner is a good approximation for the real Hessian (at least for the eigenvectors of the largest eigenvalues).
 
+The `NLOPT_LD_MMA` and `NLOPT_LD_CCSAQ` algorithms support the following internal parameters, which can be
+specified using the [`nlopt_set_param` API](NLopt_Reference#Algorithm-specific_parameters.md):
+
+* `inner_maxeval`: If ≥ 0, gives maximum number of "inner" iterations of the algorithm where it tries to ensure that its approximatations are "conservative"; defaults to `0` (no limit).   It can be useful to specify a finite number (e.g. `5` or `10`) for this parameter if inaccuracies in your gradient or objective function are preventing the algorithm from making progress.
+* `dual_algorithm` (defaults to `NLOPT_LD_MMA`), `dual_ftol_rel` (defaults to `1e-14`), `dual_ftol_abs` (defaults to `0`), `dual_xtol_rel` (defaults to `0`), `dual_xtol_abs` (defaults to `0`), `dual_maxeval` (defaults to `100000`): These specify how the algorithm internally solves the "dual" optimization problem for its approximate objective.   Because this subsidiary solve requires no evaluations of the user's objective function, it is typically fast enough that we can solve it to high precision without worrying too much about the details.  Howeve,r in high-dimensional problems you may notice that MMA/CCSA is taking a long time between optimization steps, in which case you may want to increase `dual_ftol_rel` or make other changes.   If these parameters are not specified, NLopt takes them from the [subsidiary-optimizer algorithm](NLopt_Reference#localsubsidiary-optimization-algorithm) if that has been specified, and otherwise uses the defaults indicated here.
+* `verbosity`: If ≥ 0, causes the algorithm to print internal status information on each iteration.
+
 ### SLSQP
 
 Specified in NLopt as `NLOPT_LD_SLSQP`, this is a sequential quadratic programming (SQP) algorithm for nonlinearly constrained gradient-based optimization (supporting both inequality and equality constraints), based on the implementation by Dieter Kraft and described in:

doc/docs/NLopt_C-plus-plus_Reference.md

@@ -250,6 +250,21 @@ Request the number of evaluations.
 
 In certain cases, the caller may wish to *force* the optimization to halt, for some reason unknown to NLopt. For example, if the user presses Ctrl-C, or there is an error of some sort in the objective function. You can do this by throwing *any* exception inside your objective/constraint functions: the exception will be caught, the optimization will be halted gracefully, and another exception (possibly not the same one) will be rethrown. See [Exceptions](#Exceptions.md), below. The C++ equivalent of `nlopt_forced_stop` from the [C API](NLopt_Reference#Forced_termination.md) is to throw an `nlopt::forced_stop` exception.
 
+
+Algorithm-specific parameters
+-----------------------------
+
+Certain NLopt optimization algorithms allow you to specify additional parameters by calling
+```
+nlopt_result nlopt::opt::set_param(const char *name, double val);
+bool nlopt::opt::has_param(const char *name);
+double nlopt::opt::get_param(const char *name, double defaultval);
+unsigned nlopt::opt::num_params();
+const char *nlopt::opt::nth_param(unsigned n);
+```
+where the string `name` is the name of an algorithm-specific parameter and `val` is the value you are setting the parameter to.   These functions are equivalent to the [C API](NLopt_Reference#Algorithm-specific_parameters.md) functions of the corresponding names.
+
+
 Performing the optimization
 ---------------------------
 

doc/docs/NLopt_Python_Reference.md

@@ -236,6 +236,20 @@ Request the number of evaluations.
 
 In certain cases, the caller may wish to *force* the optimization to halt, for some reason unknown to NLopt. For example, if the user presses Ctrl-C, or there is an error of some sort in the objective function. You can do this by raising *any* exception inside your objective/constraint functions:the optimization will be halted gracefully, and the same exception will be raised to the caller. See [Exceptions](#exceptions), below. The Python equivalent of `nlopt_forced_stop` from the [C API](NLopt_Reference.md#forced-termination) is to throw an `nlopt.ForcedStop` exception.
 
+Algorithm-specific parameters
+-----------------------------
+
+Certain NLopt optimization algorithms allow you to specify additional parameters by calling
+```py
+opt.set_param("name", val);
+opt.has_param("name");
+opt.get_param("name", defaultval);
+opt.num_params();
+opt.nth_param(n);
+```
+where the string `"name"` is the name of an algorithm-specific parameter and `val` is the value you are setting the parameter to.   These functions are equivalent to the [C API](NLopt_Reference#Algorithm-specific_parameters.md) functions of the corresponding names.
+
+
 Performing the optimization
 ---------------------------
 

doc/docs/NLopt_Reference.md

@@ -300,6 +300,32 @@ int nlopt_get_force_stop(nlopt_opt opt)
 
 which returns the last force-stop value that was set since the last `nlopt_optimize`. The force-stop value is reset to zero at the beginning of `nlopt_optimize`. Passing `val=0` to `nlopt_set_force_stop` tells NLopt *not* to force a halt.
 
+
+Algorithm-specific parameters
+-----------------------------
+
+Certain NLopt optimization algorithms allow you to specify additional parameters by calling
+
+```c
+nlopt_result nlopt_set_param(nlopt_opt opt, const char *name, double val);
+```
+
+where the string `name` is the name of an algorithm-specific parameter and `val` is the value you are setting the parameter to.   For example, the MMA algorithm has a parameter `"inner_maxeval"`, an upper bound on the number of "inner" iterations of the algorithm, which you can set via `nlopt_set_param(opt, "inner_maxeval", 100)`.
+
+You can also check whether a parameter is set or get the current value of a parameter with
+```c
+double nlopt_has_param(const nlopt_opt opt, const char *name);
+double nlopt_get_param(const nlopt_opt opt, const char *name, double defaultval);
+```
+where `defaultval` is returned by `nlopt_get_param` if the parameter `name` has not been set.
+
+To inspect the list of currently set parameters, you can use:
+```c
+unsigned nlopt_num_params(const nlopt_opt opt);
+const char *nlopt_nth_param(const nlopt_opt opt, unsigned n);
+```
+which return the number of set parameters and the name of the `n`-th set parameters (from `0` to `num_params-1`), respectively.
+
 Performing the optimization
 ---------------------------
 

src/algs/mma/ccsa_quadratic.c

@@ -7,17 +7,17 @@
  * 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. 
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  */
 
 /* In this file we implement Svanberg's CCSA algorithm with the
@@ -26,7 +26,7 @@
    We also allow the user to specify an optional "preconditioner": an
    approximate Hessian (which must be symmetric & positive
    semidefinite) that can be added into the approximation.  [X. Liang
-   and I went through the convergence proof in Svanberg's paper 
+   and I went through the convergence proof in Svanberg's paper
    and it does not seem to be modified by such preconditioning, as
    long as the preconditioner eigenvalues are bounded above for all x.]
 
@@ -80,7 +80,7 @@ static double dual_func(unsigned m, const double *y, double *grad, void *d_)
 {
      dual_data *d = (dual_data *) d_;
      unsigned n = d->n;
-     const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma, 
+     const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma,
 	  *dfdx = d->dfdx;
      const double *dfcdx = d->dfcdx;
      double rho = d->rho, fval = d->fval;
@@ -94,7 +94,7 @@ static double dual_func(unsigned m, const double *y, double *grad, void *d_)
 
      val = d->gval = fval;
      d->wval = 0;
-     for (i = 0; i < m; ++i) 
+     for (i = 0; i < m; ++i)
 	  val += y[i] * (gcval[i] = fcval[i]);
 
      for (j = 0; j < n; ++j) {
@@ -128,7 +128,7 @@ static double dual_func(unsigned m, const double *y, double *grad, void *d_)
 	  if (xcur[j] > ub[j]) xcur[j] = ub[j];
 	  else if (xcur[j] < lb[j]) xcur[j] = lb[j];
 	  dx = xcur[j] - x[j];
-	  
+
 	  /* function value: */
 	  dx2 = dx * dx;
 	  val += v * dx + 0.5 * u * dx2 / sigma2;
@@ -152,7 +152,7 @@ static double dual_func(unsigned m, const double *y, double *grad, void *d_)
 /* compute g(x - x0) and its gradient */
 static double gfunc(unsigned n, double f, const double *dfdx,
 		    double rho, const double *sigma,
-		    const double *x0, 
+		    const double *x0,
 		    nlopt_precond pre, void *pre_data, double *scratch,
 		    const double *x, double *grad)
 {
@@ -199,7 +199,7 @@ static void gi(unsigned m, double *result,
 	  result[i] = gfunc(n, d->fcval[i], d->dfcdx + i*n, d->rhoc[i],
 			    d->sigma,
 			    d->x,
-			    d->prec ? d->prec[i] : NULL, 
+			    d->prec ? d->prec[i] : NULL,
 			    d->prec_data ? d->prec_data[i] : NULL,
 			    d->scratch,
 			    x, grad);
@@ -212,13 +212,13 @@ nlopt_result ccsa_quadratic_minimize(
      unsigned n, nlopt_func f, void *f_data,
      unsigned m, nlopt_constraint *fc,
 
-     nlopt_precond pre, 
+     nlopt_precond pre,
 
      const double *lb, const double *ub, /* bounds */
      double *x, /* in: initial guess, out: minimizer */
      double *minf,
      nlopt_stopping *stop,
-     nlopt_opt dual_opt)
+     nlopt_opt dual_opt, int inner_maxeval, unsigned verbose)
 {
      nlopt_result ret = NLOPT_SUCCESS;
      double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
@@ -233,6 +233,8 @@ nlopt_result ccsa_quadratic_minimize(
      unsigned no_precond;
      nlopt_opt pre_opt = NULL;
 
+	 verbose = MAX(ccsa_verbose, verbose);
+
      m = nlopt_count_constraints(mfc = m, fc);
      if (nlopt_get_dimension(dual_opt) != m) {
          nlopt_stop_msg(stop, "dual optimizer has wrong dimension %d != %d",
@@ -302,7 +304,7 @@ nlopt_result ccsa_quadratic_minimize(
 	  pre_ub = pre_lb + n;
 
 	  pre_opt = nlopt_create(nlopt_get_algorithm(dual_opt), n);
-	  if (!pre_opt) { 
+	  if (!pre_opt) {
               nlopt_stop_msg(stop, "failure creating precond. optimizer");
               ret = NLOPT_FAILURE;
               goto done;
@@ -318,7 +320,7 @@ nlopt_result ccsa_quadratic_minimize(
 	  ret = nlopt_set_maxeval(pre_opt, nlopt_get_maxeval(dual_opt));
 	  if (ret < 0) goto done;
      }
-     
+
      for (j = 0; j < n; ++j) {
 	  if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
 	       sigma[j] = 1.0; /* arbitrary default */
@@ -355,7 +357,7 @@ nlopt_result ccsa_quadratic_minimize(
 	constraint weighted by 1e40 -- basically, minimizing the
 	unfeasible constraint until it becomes feasible or until we at
 	least obtain a step towards a feasible point.
-	
+
 	Svanberg suggested a different approach in his 1987 paper, basically
 	introducing additional penalty variables for unfeasible constraints,
 	but this is easier to implement and at least as efficient. */
@@ -370,11 +372,12 @@ nlopt_result ccsa_quadratic_minimize(
      nlopt_remove_equality_constraints(dual_opt);
 
      while (1) { /* outer iterations */
+	  int inner_nevals = 0;
 	  double fprev = fcur;
 	  if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 	  else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 	  else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
-	  else if (feasible && *minf < stop->minf_max) 
+	  else if (feasible && *minf < stop->minf_max)
 	       ret = NLOPT_MINF_MAX_REACHED;
 	  if (ret != NLOPT_SUCCESS) goto done;
 	  if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
@@ -393,15 +396,15 @@ nlopt_result ccsa_quadratic_minimize(
 		    ccsa_verbose = 0; /* no recursive verbosity */
 		    reti = nlopt_optimize_limited(dual_opt, y, &min_dual,
 						  0,
-						  stop->maxtime 
-						  - (nlopt_seconds() 
+						  stop->maxtime
+						  - (nlopt_seconds()
 						     - stop->start));
 		    ccsa_verbose = save_verbose;
 		    if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
 			 ret = reti;
 			 goto done;
 		    }
-		    
+
 		    dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
 	       }
 	       else {
@@ -432,50 +435,53 @@ nlopt_result ccsa_quadratic_minimize(
 		    gi(m, dd.gcval, n, xcur, NULL, &dd);
 	       }
 
-	       if (ccsa_verbose) {
+	       if (verbose) {
 		    printf("CCSA dual converged in %d iters to g=%g:\n",
 			   dd.count, dd.gval);
-		    for (i = 0; i < MIN(ccsa_verbose, m); ++i)
+		    for (i = 0; i < MIN(verbose, m); ++i)
 			 printf("    CCSA y[%u]=%g, gc[%u]=%g\n",
 				i, y[i], i, dd.gcval[i]);
 	       }
 
 	       fcur = f(n, xcur, dfdx_cur, f_data);
 	       ++ *(stop->nevals_p);
-	       if (nlopt_stop_forced(stop)) { 
+		   ++inner_nevals;
+	       if (nlopt_stop_forced(stop)) {
 		    ret = NLOPT_FORCED_STOP; goto done; }
 	       feasible_cur = 1; infeasibility_cur = 0;
 	       inner_done = dd.gval >= fcur;
 	       for (i = ifc = 0; ifc < mfc; ++ifc) {
 		    nlopt_eval_constraint(fcval_cur + i, dfcdx_cur + i*n,
 					  fc + ifc, n, xcur);
 		    i += fc[ifc].m;
-		    if (nlopt_stop_forced(stop)) { 
+		    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;
 		    for (; i < inext; ++i) {
-			 feasible_cur = feasible_cur 
+			 feasible_cur = feasible_cur
 			      && fcval_cur[i] <= fc[ifc].tol[i-i0];
-			 inner_done = inner_done && 
+			 inner_done = inner_done &&
 			      (dd.gcval[i] >= fcval_cur[i]);
 			 if (fcval_cur[i] > infeasibility_cur)
 			      infeasibility_cur = fcval_cur[i];
 		    }
 	       }
 
+		   inner_done = inner_done || (inner_maxeval > 0 && inner_nevals == inner_maxeval);
+
 	       if ((fcur < *minf && (inner_done || feasible_cur || !feasible))
 		    || (!feasible && infeasibility_cur < infeasibility)) {
-		    if (ccsa_verbose && !feasible_cur)
+		    if (verbose && !feasible_cur)
 			 printf("CCSA - using infeasible point?\n");
 		    dd.fval = *minf = fcur;
 		    infeasibility = infeasibility_cur;
 		    memcpy(fcval, fcval_cur, sizeof(double)*m);
 		    memcpy(x, xcur, sizeof(double)*n);
 		    memcpy(dfdx, dfdx_cur, sizeof(double)*n);
 		    memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
-		    
+
 		    /* once we have reached a feasible solution, the
 		       algorithm should never make the solution infeasible
 		       again (if inner_done), although the constraints may
@@ -493,7 +499,7 @@ nlopt_result ccsa_quadratic_minimize(
 	       if (nlopt_stop_forced(stop)) ret = NLOPT_FORCED_STOP;
 	       else if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
 	       else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
-	       else if (feasible && *minf < stop->minf_max) 
+	       else if (feasible && *minf < stop->minf_max)
 		    ret = NLOPT_MINF_MAX_REACHED;
 	       if (ret != NLOPT_SUCCESS) goto done;
 
@@ -503,14 +509,14 @@ nlopt_result ccsa_quadratic_minimize(
 		    rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
 	       for (i = 0; i < m; ++i)
 		    if (fcval_cur[i] > dd.gcval[i])
-			 rhoc[i] = 
-			      MIN(10*rhoc[i], 
-				  1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i]) 
+			 rhoc[i] =
+			      MIN(10*rhoc[i],
+				  1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
 					 / dd.wval));
-	       
-	       if (ccsa_verbose)
+
+	       if (verbose)
 		    printf("CCSA inner iteration: rho -> %g\n", rho);
-	       for (i = 0; i < MIN(ccsa_verbose, m); ++i)
+	       for (i = 0; i < MIN(verbose, m); ++i)
 		    printf("                CCSA rhoc[%u] -> %g\n", i,rhoc[i]);
 	  }
 
@@ -519,14 +525,14 @@ nlopt_result ccsa_quadratic_minimize(
 	  if (nlopt_stop_x(stop, xcur, xprev))
 	       ret = NLOPT_XTOL_REACHED;
 	  if (ret != NLOPT_SUCCESS) goto done;
-	       
+
 	  /* update rho and sigma for iteration k+1 */
 	  rho = MAX(0.1 * rho, CCSA_RHOMIN);
-	  if (ccsa_verbose)
+	  if (verbose)
 	       printf("CCSA outer iteration: rho -> %g\n", rho);
 	  for (i = 0; i < m; ++i)
 	       rhoc[i] = MAX(0.1 * rhoc[i], CCSA_RHOMIN);
-	  for (i = 0; i < MIN(ccsa_verbose, m); ++i)
+	  for (i = 0; i < MIN(verbose, m); ++i)
 	       printf("                 CCSA rhoc[%u] -> %g\n", i, rhoc[i]);
 	  if (k > 1) {
 	       for (j = 0; j < n; ++j) {
@@ -540,8 +546,8 @@ nlopt_result ccsa_quadratic_minimize(
 			 sigma[j] = MAX(sigma[j], 1e-8*(ub[j]-lb[j]));
 		    }
 	       }
-	       for (j = 0; j < MIN(ccsa_verbose, n); ++j)
-		    printf("                 CCSA sigma[%u] -> %g\n", 
+	       for (j = 0; j < MIN(verbose, n); ++j)
+		    printf("                 CCSA sigma[%u] -> %g\n",
 			   j, sigma[j]);
 	  }
      }